Actual source code: meshMovement.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: meshMovement.c,v 1.19 2000/10/17 13:48:57 knepley Exp $";
3: #endif
5: #include "src/gvec/meshMovementImpl.h" /*I "meshMovement.h" I*/
6: #include "src/grid/gridimpl.h" /*I "grid.h" I*//*I "gvec.h" I*/
8: extern int MeshPreCopy_Private(Mesh);
9: extern int MeshPostCopy_Private(Mesh, Mesh);
11: /*-------------------------------------------------- Mesh Functions --------------------------------------------------*/
12: #undef __FUNCT__
14: /*@
15: MeshSetMover - This function provides the MeshMover object for this mesh.
17: Not collective
19: Input Parameters:
20: + mesh - The mesh
21: - mover - The mover, or PETSC_NULL
23: Level: intermediate
25: .keywords: Mesh, movement, mover
26: .seealso: MeshGetMover(), MeshSetMovement()
27: @*/
28: int MeshSetMover(Mesh mesh, MeshMover mover) {
33: PetscObjectCompose((PetscObject) mesh, "MeshMover", (PetscObject) mover);
34: if (mover != PETSC_NULL) {
36: PetscObjectDereference((PetscObject) mover);
37: MeshMoverSetMesh(mover, mesh);
38: }
39: return(0);
40: }
42: #undef __FUNCT__
44: /*@
45: MeshGetMover - This function returns the MeshMover object for this mesh.
47: Not collective
49: Input Parameter:
50: . mesh - The mesh
52: Output Parameters:
53: . mover - The mover, or PETSC_NULL if it does not exist
55: Level: intermediate
57: .keywords: Mesh, movement, mover
58: .seealso: MeshSetMover(), MeshGetMovement()
59: @*/
60: int MeshGetMover(Mesh mesh, MeshMover *mover) {
65: PetscObjectQuery((PetscObject) mesh, "MeshMover", (PetscObject *) mover);
66: return(0);
67: }
69: #undef __FUNCT__
71: /*@
72: MeshMoveMesh - Recalculates the node coordinates based on the
73: current velocities and accelerations.
75: Collective on Mesh
77: Input Parameters:
78: . mesh - The mesh
79: . dt - The timestep
81: Level: advanced
83: .keywords: mesh, movement
84: .seealso: MeshCalcNodeVelocitiess(), MeshCalcNodeAccelerations()
85: @*/
86: int MeshMoveMesh(Mesh mesh, double dt)
87: {
88: MeshMover mover;
89: PetscTruth isMoving;
90: int ierr;
95: MeshGetMovement(mesh, &isMoving);
96: if (isMoving == PETSC_FALSE) return(0);
97: PetscLogEventBegin(MESH_MoveMesh, mesh, 0, 0, 0);
98: MeshGetMover(mesh, &mover);
99: (*mover->ops->movemesh)(mover, dt);
100: PetscLogEventEnd(MESH_MoveMesh, mesh, 0, 0, 0);
101: return(0);
102: }
104: #undef __FUNCT__
106: /*@
107: MeshUpdateNodeValues - Updates the mesh velocity, acceleration, or any other node-based quantity.
109: Collective on Mesh
111: Input Parameters:
112: + mesh - The mesh
113: - sol - The solution to an auxilliary problem
115: Output Parameters:
116: + vec - The vector which will hold he solution
117: - ghostVec - The corresponding ghost vector for vec
119: Level: developer
121: .keywords: mesh, movement
122: .seealso: MeshMoverSetNodeVelocities()
123: @*/
124: int MeshUpdateNodeValues(Mesh mesh, GVec sol, Vec vec, Vec ghostVec)
125: {
126: MeshMover mover;
127: int ierr;
134: MeshGetMover(mesh, &mover);
135: (*mover->ops->updatenodevalues)(mover, sol, vec, ghostVec);
136: return(0);
137: }
139: EXTERN_C_BEGIN
140: #undef __FUNCT__
142: /*
143: MeshPostReform_Private - Create a suitable MeshMover after reforming a Mesh.
145: Collective on Mesh
147: Input Parameters:
148: + mesh - The mesh
149: - newMesh - The reformed mesh
151: Level: developer
153: .keywords: mesh, refine, reform
154: .seealso: MeshReform()
155: */
156: int MeshPostReform_Private(Mesh mesh, Mesh newMesh)
157: {
158: MeshMover mover, newMover;
159: int ierr;
164: MeshGetMover(mesh, &mover);
165: if (mover != PETSC_NULL) {
166: MeshMoverConvert(mover, newMesh, &newMover);
167: }
168: return(0);
169: }
170: EXTERN_C_END
172: #undef __FUNCT__
174: /*@
175: MeshReform - Reform a mesh using the original boundary.
177: Collective on Mesh
179: Input Parameters:
180: + mesh - The previous mesh
181: . refine - Flag indicating whether to refine or recreate the mesh
182: . area - A function which gives an area constraint when evaluated inside an element
183: - newBd - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage
185: Output Parameter:
186: . newMesh - The new mesh
188: Level: beginner
190: .keywords: mesh, refinement, reform
191: .seealso: MeshReform(), MeshSetBoundary()
192: @*/
193: int MeshReform(Mesh mesh, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *newMesh)
194: {
195: PetscObject obj = (PetscObject) mesh;
196: Mesh refinedMesh;
197: MeshMover mover;
198: PetscReal maxArea, startX, endX, startY, endY;
199: void *tempCtx;
200: int ierr;
205: PetscObjectComposeFunction(obj, "PostCopy", "MeshPostReform_Private", (void (*)(void)) MeshPostReform_Private);
206: if (refine == PETSC_TRUE) {
207: /* Refine grid with tolerance equal to the entire size
208: This just "flips" the current mesh and does the necessary point insertion
209: */
210: MeshGetBoundingBox(mesh, &startX, &startY, PETSC_NULL, &endX, &endY, PETSC_NULL);
211: maxArea = (endX - startX)*(endY - startY);
212: MeshGetUserContext(mesh, &tempCtx);
213: MeshSetUserContext(mesh, (void *) &maxArea);
214: MeshRefine(mesh, PointFunctionConstant, newMesh);
215: MeshSetUserContext(mesh, tempCtx);
216: MeshSetUserContext(*newMesh, tempCtx);
217: } else {
218: MeshGetMover(mesh, &mover);
219: if (mover->ops->reform == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " ");
220: PetscLogEventBegin(MESH_Reform, mesh, 0, 0, 0);
221: MeshPreCopy_Private(mesh);
222: (*mover->ops->reform)(mesh, refine, area, newBd, newMesh);
223: MeshPostCopy_Private(mesh, *newMesh);
224: PetscLogEventEnd(MESH_Reform, mesh, 0, 0, 0);
225: }
227: /* Refine new mesh if necessary */
228: if (area != PETSC_NULL) {
230: PetscObjectComposeFunction((PetscObject) *newMesh, "PostCopy", "MeshPostReform_Private", (void (*)(void)) MeshPostReform_Private);
231: MeshRefine(*newMesh, area, &refinedMesh);
232: MeshDestroy(*newMesh);
233: *newMesh = refinedMesh;
234: PetscObjectComposeFunction((PetscObject) *newMesh, "PostCopy", "", PETSC_NULL);
235: }
236: PetscObjectComposeFunction(obj, "PostCopy", "", PETSC_NULL);
238: return(0);
239: }
241: /*----------------------------------------------- MeshMover Functions ------------------------------------------------*/
242: #undef __FUNCT__
244: /*@
245: MeshMoverCreate - This function creates an empty mesh mover.
247: Collective on Mesh
249: Input Parameter:
250: . mesh - The Mesh to be moved
252: Output Parameter:
253: . mover - The mover
255: Level: beginner
257: .keywords: Mesh, create, mover
258: .seealso: MeshSetUp(), MeshDestroy(), GridCreate()
259: @*/
260: int MeshMoverCreate(Mesh mesh, MeshMover *mover) {
261: MeshMover m;
262: MPI_Comm comm;
263: int ierr;
267: *mover = PETSC_NULL;
269: PetscObjectGetComm((PetscObject) mesh, &comm);
270: PetscHeaderCreate(m, _MeshMover, struct _MeshMoverOps, MESH_COOKIE, -1, "MeshMover", comm, MeshMoverDestroy, MeshMoverView);
271: PetscLogObjectCreate(m);
272: PetscLogObjectMemory(m, sizeof(struct _MeshMover));
273: PetscMemzero(m->ops, sizeof(struct _MeshMoverOps));
274: m->bops->publish = PETSC_NULL /* MeshMoverPublish_Petsc */;
275: m->type_name = PETSC_NULL;
276: m->serialize_name = PETSC_NULL;
278: m->meshVelocityMethod = MESH_LAPLACIAN;
279: m->meshVelocityGrid = PETSC_NULL;
280: m->nodeVelocities = PETSC_NULL;
281: m->nodeVelocitiesGhost = PETSC_NULL;
282: m->meshVelocityMat = PETSC_NULL;
283: m->meshVelocityRhs = PETSC_NULL;
284: m->meshVelocitySol = PETSC_NULL;
285: m->meshVelocitySles = PETSC_NULL;
286: m->meshAccelerationMethod = MESH_LAPLACIAN;
287: m->meshAccelerationGrid = PETSC_NULL;
288: m->nodeAccelerations = PETSC_NULL;
289: m->nodeAccelerationsGhost = PETSC_NULL;
290: m->meshAccelerationMat = PETSC_NULL;
291: m->meshAccelerationRhs = PETSC_NULL;
292: m->meshAccelerationSol = PETSC_NULL;
293: m->meshAccelerationSles = PETSC_NULL;
294: m->meshVelocityScatter = PETSC_NULL;
295: m->movingMeshViewer = PETSC_NULL;
296: m->movingMeshViewerCaption = "Moving Mesh";
297: m->movingMeshCtx = PETSC_NULL;
299: MeshSetMover(mesh, m);
300: *mover = m;
301: return(0);
302: }
304: #undef __FUNCT__
306: /*@
307: MeshMoverSetup - Creates all the necessary structures for mesh movement.
309: Collective on MeshMover
311: Input Parameter:
312: . mover - The mover
314: Note:
315: This function should be called directly before beginning the
316: calculation, and after all options have been set. It is
317: normally called automatically by the solver.
319: Level: advanced
321: .keywords: mesh, movement
322: .seealso: MeshMoverGetMovement()
323: @*/
324: int MeshMoverSetup(MeshMover mover)
325: {
326: Mesh mesh;
327: Partition part;
328: IS globalIS, localIS;
329: int dim;
330: int *indices;
331: int numNodes, numOverlapNodes;
332: int node, gNode, j;
333: PetscTruth reduceVel, reduceAcc;
334: PetscTruth opt;
335: int ierr;
339: /* Destroy old structures */
340: if (mover->nodeVelocities != PETSC_NULL) {
341: VecDestroy(mover->nodeVelocities);
342: VecDestroy(mover->nodeVelocitiesGhost);
343: }
344: if (mover->meshVelocityGrid != PETSC_NULL) {
345: GridDestroy(mover->meshVelocityGrid);
346: mover->meshVelocityGrid = PETSC_NULL;
347: GMatDestroy(mover->meshVelocityMat);
348: mover->meshVelocityMat = PETSC_NULL;
349: GVecDestroy(mover->meshVelocityRhs);
350: mover->meshVelocityRhs = PETSC_NULL;
351: GVecDestroy(mover->meshVelocitySol);
352: mover->meshVelocitySol = PETSC_NULL;
353: SLESDestroy(mover->meshVelocitySles);
354: mover->meshVelocitySles = PETSC_NULL;
355: }
356: if (mover->nodeAccelerations != PETSC_NULL) {
357: VecDestroy(mover->nodeAccelerations);
358: VecDestroy(mover->nodeAccelerationsGhost);
359: }
360: if (mover->meshAccelerationGrid != PETSC_NULL) {
361: GridDestroy(mover->meshAccelerationGrid);
362: mover->meshAccelerationGrid = PETSC_NULL;
363: GMatDestroy(mover->meshAccelerationMat);
364: mover->meshAccelerationMat = PETSC_NULL;
365: GVecDestroy(mover->meshAccelerationRhs);
366: mover->meshAccelerationRhs = PETSC_NULL;
367: GVecDestroy(mover->meshAccelerationSol);
368: mover->meshAccelerationSol = PETSC_NULL;
369: SLESDestroy(mover->meshAccelerationSles);
370: mover->meshAccelerationSles = PETSC_NULL;
371: }
372: if (mover->meshVelocityScatter != PETSC_NULL) {
373: VecScatterDestroy(mover->meshVelocityScatter);
374: }
375: /* Determine solution method */
376: ierr = PetscOptionsHasName(mover->prefix, "-mesh_move_lap", &opt);
377: if (opt == PETSC_TRUE) {
378: mover->meshVelocityMethod = MESH_LAPLACIAN;
379: mover->meshAccelerationMethod = MESH_LAPLACIAN;
380: }
381: ierr = PetscOptionsHasName(mover->prefix, "-mesh_move_weighted_lap", &opt);
382: if (opt == PETSC_TRUE) {
383: mover->meshVelocityMethod = MESH_WEIGHTED_LAPLACIAN;
384: mover->meshAccelerationMethod = MESH_WEIGHTED_LAPLACIAN;
385: }
386: /* Setup mesh node velocities vector */
387: MeshMoverGetMesh(mover, &mesh);
388: MeshGetPartition(mesh, &part);
389: MeshGetDimension(mesh, &dim);
390: MeshGetInfo(mesh, PETSC_NULL, &numNodes, PETSC_NULL, PETSC_NULL);
391: VecCreateMPI(mover->comm, numNodes*dim, PETSC_DECIDE, &mover->nodeVelocities);
392: PartitionGetNumOverlapNodes(part, &numOverlapNodes);
393: VecCreateSeq(PETSC_COMM_SELF, numOverlapNodes*dim, &mover->nodeVelocitiesGhost);
394: /* Setup mesh node accelerations vector */
395: VecCreateMPI(mover->comm, numNodes*dim, PETSC_DECIDE, &mover->nodeAccelerations);
396: VecCreateSeq(PETSC_COMM_SELF, numOverlapNodes*dim, &mover->nodeAccelerationsGhost);
397: /* Setup grids and other implementation specific stuff */
398: (*mover->ops->setup)(mover);
399: /* Setup scatter from node values vector to ghost vector */
400: ISCreateStride(mover->comm, numOverlapNodes*dim, 0, 1, &localIS);
401: PetscMalloc(numOverlapNodes*dim * sizeof(int), &indices);
402: for(node = 0; node < numOverlapNodes; node++) {
403: PartitionLocalToGlobalNodeIndex(part, node, &gNode);
404: for(j = 0; j < dim; j++)
405: indices[node*dim+j] = gNode*dim+j;
406: }
407: ISCreateGeneral(mover->comm, numOverlapNodes*dim, indices, &globalIS);
408: VecScatterCreate(mover->nodeVelocities, globalIS, mover->nodeVelocitiesGhost, localIS, &mover->meshVelocityScatter);
409:
410: ISDestroy(localIS);
411: ISDestroy(globalIS);
412: PetscFree(indices);
413: /* Setup viewer */
414: PetscOptionsHasName(mover->prefix, "-mesh_moving_view", &opt);
415: if ((opt == PETSC_TRUE) && (mover->movingMeshViewer == PETSC_NULL)) {
416: mover->movingMeshViewer = PETSC_VIEWER_STDOUT_WORLD;
417: }
418: PetscOptionsHasName(mover->prefix, "-mesh_moving_view_draw", &opt);
419: if ((opt == PETSC_TRUE) && (mover->movingMeshViewer == PETSC_NULL)) {
420: PetscReal startX, startY, endX, endY, sizeX, sizeY;
421: int xsize, ysize;
423: ierr = MeshGetBoundingBox(mesh, &startX, &startY, 0, &endX, &endY, 0);
424: sizeX = endX - startX;
425: sizeY = endY - startY;
426: if (sizeX > sizeY) {
427: ysize = 300;
428: xsize = ysize * (int) (sizeX/sizeY);
429: } else {
430: xsize = 300;
431: ysize = xsize * (int) (sizeY/sizeX);
432: }
433: PetscViewerDrawOpen(mover->comm, 0, 0, 0, 315, xsize, ysize, &mover->movingMeshViewer);
434: }
435: /* Check grids */
436: GridGetReduceSystem(mover->meshVelocityGrid, &reduceVel);
437: GridGetReduceSystem(mover->meshAccelerationGrid, &reduceAcc);
438: if (reduceVel != reduceAcc) {
439: SETERRQ(PETSC_ERR_ARG_WRONG, "Velocity and acceleration grids must have same BC reduction");
440: }
441: return(0);
442: }
444: #undef __FUNCT__
446: /*
447: MeshMoverSetTypeFromOptions - Sets the type of mesh movement from user options.
449: Collective on MeshMover
451: Input Parameter:
452: . mover - The mover
454: Level: intermediate
456: .keywords: MeshMover, set, options, database, type
457: .seealso: MeshMoverSetFromOptions(), MeshMoverSetType()
458: */
459: static int MeshMoverSetTypeFromOptions(MeshMover mover)
460: {
461: Mesh mesh;
462: PetscTruth opt;
463: char *defaultType;
464: char typeName[256];
465: int dim;
466: int ierr;
469: MeshMoverGetMesh(mover, &mesh);
470: MeshGetDimension(mesh, &dim);
471: if (mover->type_name != PETSC_NULL) {
472: defaultType = mover->type_name;
473: } else {
474: switch(dim)
475: {
476: case 2:
477: defaultType = MESH_TRIANGULAR_2D;
478: break;
479: default:
480: SETERRQ1(PETSC_ERR_SUP, "MeshMover dimension %d is not supported", dim);
481: }
482: }
484: if (!MeshMoverRegisterAllCalled) {
485: MeshMoverRegisterAll(PETSC_NULL);
486: }
487: PetscOptionsList("-mesh_mover_type", "Mesh movement method", "MeshMoverSetType", MeshMoverList, defaultType, typeName, 256, &opt);
488:
489: if (opt == PETSC_TRUE) {
490: MeshMoverSetType(mover, typeName);
491: } else {
492: MeshMoverSetType(mover, defaultType);
493: }
494: return(0);
495: }
497: #undef __FUNCT__
499: /*@
500: MeshMoverSetFromOptions - Sets various MeshMover parameters from user options.
502: Collective on MeshMover
504: Input Parameter:
505: . mover - The mover
507: Notes: To see all options, run your program with the -help option, or consult the users manual.
508: Must be called after MeshMoverCreate() but before the vector is used.
510: Level: intermediate
512: .keywords: MeshMover, set, options, database
513: .seealso: MeshMoverCreate(), MeshMoverPrintHelp(), MeshMoverSetOptionsPrefix()
514: @*/
515: int MeshMoverSetFromOptions(MeshMover mover)
516: {
517: PetscTruth opt;
518: int ierr;
522: PetscOptionsBegin(mover->comm, mover->prefix, "MeshMover options", "MeshMover");
524: /* Handle generic mesh options */
525: PetscOptionsHasName(PETSC_NULL, "-help", &opt);
526: if (opt == PETSC_TRUE) {
527: MeshMoverPrintHelp(mover);
528: }
530: /* Handle mesh type options */
531: MeshMoverSetTypeFromOptions(mover);
533: /* Handle specific mesh options */
534: if (mover->ops->setfromoptions != PETSC_NULL) {
535: (*mover->ops->setfromoptions)(mover);
536: }
538: PetscOptionsEnd();
540: /* Handle subobject options */
541: return(0);
542: }
544: #undef __FUNCT__
546: /*@
547: MeshMoverPrintHelp - Prints all options for the MeshMover.
549: Input Parameter:
550: . mover - The mover
552: Options Database Keys:
553: $ -help, -h
555: Level: intermediate
557: .keywords: MeshMover, help
558: .seealso: MeshMoverSetFromOptions()
559: @*/
560: int MeshMoverPrintHelp(MeshMover mover)
561: {
562: char p[64];
563: int ierr;
568: PetscStrcpy(p, "-");
569: if (mover->prefix != PETSC_NULL) {
570: PetscStrcat(p, mover->prefix);
571: }
573: (*PetscHelpPrintf)(mover->comm, "MeshMover options ------------------------------------------------n");
574: (*PetscHelpPrintf)(mover->comm," %smesh_mover_type <typename> : Sets the mesh movement typen", p);
575: return(0);
576: }
578: #undef __FUNCT__
580: /*@
581: MeshMoverDuplicate - This function returns an identical MeshMover.
583: Not collective
585: Input Parameter:
586: . mover - The mover
588: Output Parameter:
589: . newMover - The new mover
591: Level: intermediate
593: .keywords: Mesh, movement, mover
594: .seealso: MeshMoverCopy()
595: @*/
596: int MeshMoverDuplicate(MeshMover mover, MeshMover *newMover) {
597: Mesh mesh;
598: int ierr;
601: MeshMoverGetMesh(mover, &mesh);
602: MeshMoverConvert(mover, mesh, newMover);
603: return(0);
604: }
606: #undef __FUNCT__
608: /*@
609: MeshMoverConvert - This function returns an identical MeshMover on the new Mesh.
611: Not collective
613: Input Parameters:
614: + mover - The mover
615: - newMesh - The new Mesh
617: Output Parameter:
618: . newMover - The new mover
620: Level: intermediate
622: .keywords: Mesh, movement, mover
623: .seealso: MeshMoverCopy()
624: @*/
625: int MeshMoverConvert(MeshMover mover, Mesh newMesh, MeshMover *newMover) {
626: MeshMover m;
627: Mesh mesh;
628: int ierr;
631: MeshMoverGetMesh(mover, &mesh);
632: if (mesh == newMesh) {
633: /* Need to play reference games since MeshMoverCreate automatically sets it in the Mesh */
634: PetscObjectReference((PetscObject) mover);
635: MeshMoverCreate(newMesh, &m);
636: PetscObjectReference((PetscObject) m);
637: MeshSetMover(mesh, mover);
638: } else {
639: MeshMoverCreate(newMesh, &m);
640: }
642: MeshMoverSetFromOptions(m);
643: m->meshVelocityMethod = mover->meshVelocityMethod;
644: m->meshAccelerationMethod = mover->meshAccelerationMethod;
645: m->movingMeshViewer = mover->movingMeshViewer;
646: if (mover->movingMeshViewer != PETSC_NULL) PetscObjectReference((PetscObject) mover->movingMeshViewer);
647: m->movingMeshViewerCaption = mover->movingMeshViewerCaption;
648: m->movingMeshCtx = mover->movingMeshCtx;
649: MeshMoverSetup(m);
650: GridDuplicateBC(mover->meshVelocityGrid, m->meshVelocityGrid);
651: GridDuplicateBC(mover->meshAccelerationGrid, m->meshAccelerationGrid);
652: *newMover = m;
653: return(0);
654: }
656: #undef __FUNCT__
658: /*@
659: MeshMoverView - Views a mesh mover object.
661: Collective on MeshMover
663: Input Parameters:
664: + mover - The mover
665: - viewer - The viewer with which to view the mesh mover
667: Level: beginner
669: .keywords: mesh, view, mover
670: .seealso: MeshMoverDestroy(), PetscViewerDrawOpen()
671: @*/
672: int MeshMoverView(MeshMover mover, PetscViewer viewer)
673: {
676: if (viewer == PETSC_NULL) {
677: viewer = PETSC_VIEWER_STDOUT_SELF;
678: } else {
680: }
681: return(0);
682: }
684: #undef __FUNCT__
686: /*@
687: MeshMoverDestroy - Destroys a mesh mover object.
689: Collective on MeshMover
691: Input Parameter:
692: . mover - The mover
694: Level: beginner
696: .keywords: mesh, destroy, mover
697: .seealso MeshMoverCreate()
698: @*/
699: int MeshMoverDestroy(MeshMover mover) {
704: if (--mover->refct > 0) return(0);
706: if (mover->nodeVelocities != PETSC_NULL) {
707: GVecDestroy(mover->nodeVelocities);
708: GVecDestroy(mover->nodeVelocitiesGhost);
709: }
710: if (mover->meshVelocitySles != PETSC_NULL) {
711: GMatDestroy(mover->meshVelocityMat);
712: GVecDestroy(mover->meshVelocityRhs);
713: GVecDestroy(mover->meshVelocitySol);
714: SLESDestroy(mover->meshVelocitySles);
715: }
716: if (mover->meshVelocityGrid != PETSC_NULL) {
717: GridFinalizeBC(mover->meshVelocityGrid);
718: GridDestroy(mover->meshVelocityGrid);
719: }
720: if (mover->nodeAccelerations != PETSC_NULL) {
721: GVecDestroy(mover->nodeAccelerations);
722: GVecDestroy(mover->nodeAccelerationsGhost);
723: }
724: if (mover->meshAccelerationSles != PETSC_NULL) {
725: GMatDestroy(mover->meshAccelerationMat);
726: GVecDestroy(mover->meshAccelerationRhs);
727: GVecDestroy(mover->meshAccelerationSol);
728: SLESDestroy(mover->meshAccelerationSles);
729: }
730: if (mover->meshAccelerationGrid != PETSC_NULL) {
731: GridFinalizeBC(mover->meshAccelerationGrid);
732: GridDestroy(mover->meshAccelerationGrid);
733: }
734: if (mover->meshVelocityScatter != PETSC_NULL) {
735: VecScatterDestroy(mover->meshVelocityScatter);
736: }
737: if (mover->movingMeshViewer != PETSC_NULL) {
738: PetscViewerDestroy(mover->movingMeshViewer);
739: }
741: (*mover->ops->destroy)(mover);
743: PetscLogObjectDestroy(mover);
744: PetscHeaderDestroy(mover);
745: return(0);
746: }
748: #undef __FUNCT__
750: /*@
751: MeshMoverSetMesh - This function returns the associated Mesh.
753: Not collective
755: Input Parameter:
756: . mover - The mover
758: Output Parameter:
759: . mesh - The mesh, or PETSC_NULL if it does not exist
761: Level: intermediate
763: .keywords: Mesh, movement, mover
764: .seealso: MeshSetMover()
765: @*/
766: int MeshMoverSetMesh(MeshMover mover, Mesh mesh) {
770: mover->mesh = mesh;
771: return(0);
772: }
774: #undef __FUNCT__
776: /*@
777: MeshMoverGetMesh - This function returns the associated Mesh.
779: Not collective
781: Input Parameter:
782: . mover - The mover
784: Output Parameter:
785: . mesh - The mesh, or PETSC_NULL if it does not exist
787: Level: intermediate
789: .keywords: Mesh, movement, mover
790: .seealso: MeshGetMover()
791: @*/
792: int MeshMoverGetMesh(MeshMover mover, Mesh *mesh) {
796: *mesh = mover->mesh;
797: return(0);
798: }
800: #undef __FUNCT__
802: /*@
803: MeshMoverGetMovement - Returns the type of mesh movement calculation.
805: Not collective
807: Input Parameter:
808: . mover - The mover
810: Output Parameters:
811: + vtype - The type of mesh velocity calculation, like MESH_LAPLACIAN
812: . atype - The type of mesh acceleration calculation, like MESH_LAPLACIAN
813: - ctx - A user context
815: Level: intermediate
817: .keywords: mesh, movement
818: .seealso: MeshMoverSetMovement()
819: @*/
820: int MeshMoverGetMovement(MeshMover mover, MeshSolveMethod *vtype, MeshSolveMethod *atype, PetscObject *ctx)
821: {
824: if (vtype != PETSC_NULL) {
826: *vtype = mover->meshVelocityMethod;
827: }
828: if (atype != PETSC_NULL) {
830: *atype = mover->meshAccelerationMethod;
831: }
832: if (ctx != PETSC_NULL) {
834: *ctx = mover->movingMeshCtx;
835: }
836: return(0);
837: }
839: #undef __FUNCT__
841: /*@
842: MeshMoverSetMovement - Determines the type of mesh mvoement calculation
844: Collective on MeshMover
846: Input Parameters:
847: + mover - The mover
848: . vtype - The type of mesh velocity calculation, like MESH_LAPLACIAN
849: . atype - The type of mesh acceleration calculation, like MESH_LAPLACIAN
850: - ctx - A user context
852: Level: intermediate
854: .keywords: mesh, movement
855: .seealso: MeshMoverGetMovement()
856: @*/
857: int MeshMoverSetMovement(MeshMover mover, MeshSolveMethod vtype, MeshSolveMethod atype, PetscObject ctx)
858: {
863: if (vtype != PETSC_DECIDE) mover->meshVelocityMethod = vtype;
864: if (atype != PETSC_DECIDE) mover->meshAccelerationMethod = atype;
865: if (ctx != PETSC_NULL) mover->movingMeshCtx = ctx;
866: if (mover->meshVelocityGrid == PETSC_NULL) {
867: MeshMoverSetup(mover);
868: }
869: return(0);
870: }
872: #undef __FUNCT__
874: /*@
875: MeshMoverGetVelocityGrid - Returns the grid which defines the mesh velocity.
877: Not collective
879: Input Parameter:
880: . mover - The mover
882: Output Parameter:
883: . velGrid - The mesh velocity grid
885: Level: intermediate
887: .keywords: grid, mesh, mesh velocity, movement
888: .seealso: MeshMoverGetAccelerationGrid(), GridGetMeshMovement()
889: @*/
890: int MeshMoverGetVelocityGrid(MeshMover mover, Grid *velGrid)
891: {
894: *velGrid = mover->meshVelocityGrid;
895: return(0);
896: }
898: #undef __FUNCT__
900: /*@
901: MeshMoverGetAccelerationGrid - Returns the grid which defines the mesh acceleration.
903: Not collective
905: Input Parameter:
906: . mover - The mover
908: Output Parameter:
909: . accGrid - The mesh acceleration grid
911: Level: intermediate
913: .keywords: grid, mesh, mesh acceleration, movement
914: .seealso: MeshMoverGetVelocityGrid(), GridGetMeshMovement()
915: @*/
916: int MeshMoverGetAccelerationGrid(MeshMover mover, Grid *accGrid)
917: {
920: *accGrid = mover->meshAccelerationGrid;
921: return(0);
922: }
924: #undef __FUNCT__
926: /*@
927: MeshMoverCalcVelocityBCValues - This function calculates the boundary values for the
928: mesh velocity. It is normally called once a timestep when using time dependent boundary
929: conditions.
931: Collective on MeshMover
933: Input Parameters:
934: + mover - The MeshMover
935: . save - A flag used to store old values, usually for timestepping
936: - ctx - The context
938: Level: advanced
940: .keywords: MeshMover, reduction, boundary conditions
941: .seealso: MeshMoverSetVelocityBCContext(), MeshMoverSetVelocityBC()
942: @*/
943: int MeshMoverCalcVelocityBCValues(MeshMover mover)
944: {
945: Grid g;
946: PetscTruth reduce;
947: void *reduceCtx;
948: int ierr;
952: MeshMoverGetMovement(mover, PETSC_NULL, PETSC_NULL, PETSC_NULL);
953: MeshMoverGetVelocityGrid(mover, &g);
955: GridGetReduceSystem(g, &reduce);
956: if (reduce == PETSC_FALSE) return(0);
957: GridSetUp(g);
958: GridGetBCContext(g, &reduceCtx);
959: GridCalcBCValues(g, PETSC_FALSE, reduceCtx);
960: return(0);
961: }
963: #undef __FUNCT__
965: /*@
966: MeshMoverCalcAccelerationBCValues - This function calculates the boundary values for the
967: mesh acceleration. It is normally called once a timestep when using time dependent boundary
968: conditions.
970: Collective on MeshMover
972: Input Parameters:
973: + mover - The MeshMover
974: . save - A flag used to store old values, usually for timestepping
975: - ctx - The context
977: Level: advanced
979: .keywords: MeshMover, reduction, boundary conditions
980: .seealso: MeshMoverSetAccelerationBCContext(), MeshMoverSetAccelerationBC()
981: @*/
982: int MeshMoverCalcAccelerationBCValues(MeshMover mover)
983: {
984: Grid g;
985: PetscTruth reduce;
986: void *reduceCtx;
987: int ierr;
991: MeshMoverGetMovement(mover, PETSC_NULL, PETSC_NULL, PETSC_NULL);
992: MeshMoverGetAccelerationGrid(mover, &g);
994: GridGetReduceSystem(g, &reduce);
995: if (reduce == PETSC_FALSE) return(0);
996: GridSetUp(g);
997: GridGetBCContext(g, &reduceCtx);
998: GridCalcBCValues(g, PETSC_FALSE, reduceCtx);
999: return(0);
1000: }
1002: #undef __FUNCT__
1004: /*@
1005: MeshMoverCalcNodeVelocities - Recalculates the node velocities.
1007: Collective on MeshMover
1009: Input Parameters:
1010: + mover - The mover
1011: - flag - Flag for reuse of old operator in calculation
1013: Level: advanced
1015: .keywords: mesh, movement
1016: .seealso: MeshMoverCalcNodeAccelerations()
1017: @*/
1018: int MeshMoverCalcNodeVelocities(MeshMover mover, PetscTruth flag)
1019: {
1020: Grid grid = mover->meshVelocityGrid;
1021: MatStructure matStructFlag = DIFFERENT_NONZERO_PATTERN;
1022: Mesh mesh;
1023: SLES sles;
1024: GMat mat;
1025: GVec rhs, sol;
1026: PetscScalar zero = 0.0;
1027: int its;
1028: int ierr;
1033: PetscLogEventBegin(MESH_CalcNodeVelocities, mover, 0, 0, 0);
1035: /* Allocation */
1036: if (mover->meshVelocitySles == PETSC_NULL) {
1037: GVecCreate(grid, &mover->meshVelocityRhs);
1038: PetscLogObjectParent(mover, mover->meshVelocityRhs);
1039: GVecCreate(grid, &mover->meshVelocitySol);
1040: PetscLogObjectParent(mover, mover->meshVelocitySol);
1041: GMatCreate(grid, &mover->meshVelocityMat);
1042: PetscLogObjectParent(mover, mover->meshVelocityMat);
1043: SLESCreate(mover->comm, &mover->meshVelocitySles);
1044: SLESAppendOptionsPrefix(mover->meshVelocitySles, "mesh_vel_");
1045: SLESSetFromOptions(mover->meshVelocitySles);
1046: PetscLogObjectParent(mover, mover->meshVelocitySles);
1047: }
1048: sles = mover->meshVelocitySles;
1049: mat = mover->meshVelocityMat;
1050: rhs = mover->meshVelocityRhs;
1051: sol = mover->meshVelocitySol;
1053: /* Form System Matrix */
1054: if (flag == PETSC_TRUE) {
1055: MatZeroEntries(mat);
1056: GridEvaluateSystemMatrix(grid, PETSC_NULL, &mat, &mat, &matStructFlag, mover->movingMeshCtx);
1057: GMatSetBoundary(mat, 1.0, mover->movingMeshCtx);
1058: MatSetOption(mat, MAT_NO_NEW_NONZERO_LOCATIONS);
1059: }
1060: /* Form Rhs */
1061: VecSet(&zero, rhs);
1062: GridEvaluateRhs(grid, PETSC_NULL, rhs, mover->movingMeshCtx);
1063: GVecSetBoundary(rhs, mover->movingMeshCtx);
1064: /* Solve system */
1065: SLESSetOperators(sles, mat, mat, SAME_NONZERO_PATTERN);
1066: SLESSolve(sles, rhs, sol, &its);
1067: /* Copy values into the mesh vector */
1068: MeshMoverGetMesh(mover, &mesh);
1069: MeshUpdateNodeValues(mesh, sol, mover->nodeVelocities, mover->nodeVelocitiesGhost);
1070: PetscLogEventEnd(MESH_CalcNodeVelocities, mover, 0, 0, 0);
1071: return(0);
1072: }
1074: #undef __FUNCT__
1076: /*@
1077: MeshMoverCalcNodeAccelerations - Recalculates the node accelerations.
1079: Collective on MeshMover
1081: Input Parameters:
1082: + mover - The mover
1083: - flag - Flag for reuse of old operator in calculation
1085: Level: advanced
1087: .keywords: mesh, movement
1088: .seealso: MeshMoverCalcNodeVelocities()
1089: @*/
1090: int MeshMoverCalcNodeAccelerations(MeshMover mover, PetscTruth flag)
1091: {
1092: Grid grid = mover->meshAccelerationGrid;
1093: MatStructure matStructFlag = DIFFERENT_NONZERO_PATTERN;
1094: Mesh mesh;
1095: SLES sles;
1096: GMat mat;
1097: GVec rhs, sol;
1098: PetscScalar zero = 0.0;
1099: int its;
1100: int ierr;
1105: PetscLogEventBegin(MESH_CalcNodeAccelerations, mover, 0, 0, 0);
1107: /* Allocation */
1108: if (mover->meshAccelerationSles == PETSC_NULL) {
1109: GVecCreate(grid, &mover->meshAccelerationRhs);
1110: PetscLogObjectParent(mover, mover->meshAccelerationRhs);
1111: GVecCreate(grid, &mover->meshAccelerationSol);
1112: PetscLogObjectParent(mover, mover->meshAccelerationSol);
1113: GMatCreate(grid, &mover->meshAccelerationMat);
1114: PetscLogObjectParent(mover, mover->meshAccelerationMat);
1115: SLESCreate(mover->comm, &mover->meshAccelerationSles);
1116: SLESAppendOptionsPrefix(mover->meshAccelerationSles, "mesh_accel_");
1117: SLESSetFromOptions(mover->meshAccelerationSles);
1118: PetscLogObjectParent(mover, mover->meshAccelerationSles);
1119: }
1120: sles = mover->meshAccelerationSles;
1121: mat = mover->meshAccelerationMat;
1122: rhs = mover->meshAccelerationRhs;
1123: sol = mover->meshAccelerationSol;
1125: /* Form System Matrix */
1126: if (flag == PETSC_TRUE) {
1127: MatZeroEntries(mat);
1128: GridEvaluateSystemMatrix(grid, PETSC_NULL, &mat, &mat, &matStructFlag, mover->movingMeshCtx);
1129: GMatSetBoundary(mat, 1.0, mover->movingMeshCtx);
1130: MatSetOption(mat, MAT_NO_NEW_NONZERO_LOCATIONS);
1131: }
1132: /* Form Rhs */
1133: VecSet(&zero, rhs);
1134: GridEvaluateRhs(grid, PETSC_NULL, rhs, mover->movingMeshCtx);
1135: GVecSetBoundary(rhs, mover->movingMeshCtx);
1136: /* Solve system */
1137: SLESSetOperators(sles, mat, mat, SAME_NONZERO_PATTERN);
1138: SLESSolve(sles, rhs, sol, &its);
1139: /* Copy values into the mesh vector */
1140: MeshMoverGetMesh(mover, &mesh);
1141: MeshUpdateNodeValues(mesh, sol, mover->nodeAccelerations, mover->nodeAccelerationsGhost);
1142: PetscLogEventEnd(MESH_CalcNodeAccelerations, mover, 0, 0, 0);
1143: return(0);
1144: }
1146: #undef __FUNCT__
1148: /*@
1149: MeshMoverSetNodeVelocities - Allows the user to specify node velocities.
1151: Collective on MeshMover
1153: Input Parameters:
1154: + mover - The mover
1155: . func - A PointFunction defining the node velocities
1156: . alpha - A scalar multiplier
1157: - ctx - A user context
1159: Level: advanced
1161: .keywords: mesh, movement
1162: .seealso: MeshMoverSetNodeAccelerations()
1163: @*/
1164: int MeshMoverSetNodeVelocities(MeshMover mover, PointFunction func, PetscScalar alpha, void *ctx)
1165: {
1166: #if 0
1167: int field = 0;
1171: GVecEvaluateFunction(mover->nodeVelocities, 1, &field, func, alpha, ctx);
1172: /* Now need a way to get mesh velocities back to the system */
1173: return(0);
1174: #else
1175: SETERRQ(PETSC_ERR_SUP, "Not yet a grid vector");
1176: #endif
1177: }
1179: #undef __FUNCT__
1181: /*@
1182: MeshMoverSetNodeAccelerations - Allows the user to specify node accelerations.
1184: Collective on MeshMover
1186: Input Parameters:
1187: + mover - The mover
1188: . func - A function defining the node accelerations
1189: . alpha - A scalar multiplier
1190: - ctx - A user context
1192: Level: intermediate
1194: .keywords: mesh, movement
1195: .seealso: MeshMoverSetNodeVelocities()
1196: @*/
1197: int MeshMoverSetNodeAccelerations(MeshMover mover, PointFunction func, PetscScalar alpha, void *ctx)
1198: {
1199: #if 0
1200: int field = 0;
1204: GVecEvaluateFunction(mover->nodeAccelerations, 1, &field, func, alpha, ctx);
1205: /* Now need a way to get mesh accelerations back to the system */
1206: return(0);
1207: #else
1208: SETERRQ(PETSC_ERR_SUP, "Not yet a grid vector");
1209: #endif
1210: }
1212: #undef __FUNCT__
1214: /*@C
1215: MeshMoverSetVelocityBC - This function sets the boundary conditions to use
1216: for the mesh velocity calculation.
1218: Collective on MeshMover
1220: Input Parameters:
1221: + mover - The mover
1222: . bd - The marker for the boundary along which conditions are applied
1223: . f - The function which defines the boundary condition
1224: - reduce - The flag for explicit reduction of the system
1226: Note:
1227: This function clears any previous conditions, whereas MeshMoverAddVelocityBC()
1228: appends the condition.
1230: Level: intermediate
1232: .keywords mesh, velocity, boundary condition
1233: .seealso MeshMoverAddVelocityBC(), MeshMoverSetAccelerationBC(), GridSetBC(), GridAddBC()
1234: @*/
1235: int MeshMoverSetVelocityBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1236: {
1242: if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created");
1243: GridSetBC(mover->meshVelocityGrid, bd, 0, f, reduce);
1244: return(0);
1245: }
1247: #undef __FUNCT__
1249: /*@C
1250: MeshMoverAddVelocityBC - This function adds a boundary condition to use
1251: for the mesh velocity calculation.
1253: Collective on MeshMover
1255: Input Parameters:
1256: + mover - The mover
1257: . bd - The marker for the boundary along which conditions are applied
1258: . f - The function which defines the boundary condition
1259: - reduce - The flag for explicit reduction of the system
1261: Note:
1262: This function appends the condition, whereas MeshMoverSetVelocityBC()
1263: clears any previous conditions.
1265: Level: intermediate
1267: .keywords mesh, velocity, boundary condition
1268: .seealso MeshMoverSetVelocityBC(), MeshMoverAddAccelerationBC(), GridSetBC(), GridAddBC()
1269: @*/
1270: int MeshMoverAddVelocityBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1271: {
1277: if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created");
1278: GridAddBC(mover->meshVelocityGrid, bd, 0, f, reduce);
1279: return(0);
1280: }
1282: #undef __FUNCT__
1284: /*@C
1285: MeshMoverSetVelocityBCContext - This function sets the boundary condition context
1286: to use for the mesh velocity calculation.
1288: Collective on MeshMover
1290: Input Parameters:
1291: + mover - The mover
1292: - ctx - The user context
1294: Level: intermediate
1296: .keywords mesh, velocity, boundary condition
1297: .seealso MeshMoverAddVelocityBC(), MeshMoverSetVelocityBC()
1298: @*/
1299: int MeshMoverSetVelocityBCContext(MeshMover mover, void *ctx)
1300: {
1306: if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created");
1307: GridSetBCContext(mover->meshVelocityGrid, ctx);
1308: return(0);
1309: }
1311: #undef __FUNCT__
1313: /*@C
1314: MeshMoverSetAccelerationBC - This function sets the boundary condition to use
1315: for the mesh acceleration calculation.
1317: Collective on MeshMover
1319: Input Parameters:
1320: + mover - The mover
1321: . bd - The marker for the boundary along which conditions are applied
1322: . f - The function which defines the boundary condition
1323: - reduce - The flag for explicit reduction of the system
1325: Note:
1326: This function clears any previous conditions, whereas MeshMoverAddAccelerationBC()
1327: appends the condition.
1329: Level: intermediate
1331: .keywords mesh, velocity, boundary condition
1332: .seealso MeshMoverAddAccelerationBC(), MeshMoverSetVelocityBC(), GridSetBC(), GridAddBC()
1333: @*/
1334: int MeshMoverSetAccelerationBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1335: {
1341: if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created");
1342: GridSetBC(mover->meshAccelerationGrid, bd, 0, f, reduce);
1343: return(0);
1344: }
1346: #undef __FUNCT__
1348: /*@C
1349: MeshMoverAddAccelerationBC - This function adds a boundary condition to use
1350: for the mesh acceleration calculation.
1352: Collective on MeshMover
1354: Input Parameters:
1355: + mover - The mover
1356: . bd - The marker for the boundary along which conditions are applied
1357: . f - The function which defines the boundary condition
1358: - reduce - The flag for explicit reduction of the system
1360: Note:
1361: This function appends the condition, whereas MeshMoverSetVelocityBC()
1362: clears any previous conditions.
1364: Level: intermediate
1366: .keywords mesh, velocity, boundary condition
1367: .seealso MeshMoverSetAccelerationBC(), MeshMoverAddVelocityBC(), GridSetBC(), GridAddBC()
1368: @*/
1369: int MeshMoverAddAccelerationBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1370: {
1376: if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created");
1377: GridAddBC(mover->meshAccelerationGrid, bd, 0, f, reduce);
1378: return(0);
1379: }
1381: #undef __FUNCT__
1383: /*@C
1384: MeshMoverSetAccelerationBCContext - This function sets the boundary condition context
1385: to use for the mesh acceleration calculation.
1387: Collective on MeshMover
1389: Input Parameter:
1390: + mover - The mover
1391: - ctx - The user context
1393: Level: intermediate
1395: .keywords mesh, velocity, boundary condition
1396: .seealso MeshMoverAddAccelerationBC(), MeshMoverSetAccelerationBC()
1397: @*/
1398: int MeshMoverSetAccelerationBCContext(MeshMover mover, void *ctx)
1399: {
1405: if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created");
1406: GridSetBCContext(mover->meshAccelerationGrid, ctx);
1407: return(0);
1408: }
1410: #undef __FUNCT__
1412: /*@
1413: MeshMoverSetMovementCaption - Sets the caption of the moving mesh
1414: viewer. Can be used to provide auxilliary information, e.g. time.
1416: Collective on MeshMover
1418: Input Parameters:
1419: + mover - The mover
1420: - caption - The caption
1422: Level: intermediate
1424: .keywords: mesh, movement
1425: .seealso: MeshMoverMoveMesh()
1426: @*/
1427: int MeshMoverSetMovementCaption(MeshMover mover, char *caption)
1428: {
1432: mover->movingMeshViewerCaption = caption;
1433: return(0);
1434: }
1436: /*----------------------------------------- MeshMover Registration Functions -----------------------------------------*/
1437: PetscFList MeshMoverList = 0;
1438: int MeshMoverRegisterAllCalled = 0;
1440: #undef __FUNCT__
1442: /*@C
1443: MeshMoverSetType - Builds a MeshMover, for a particular Mesh implementation.
1445: Collective on MeshMover
1447: Input Parameters:
1448: + mover - The MeshMover
1449: - type - The type name
1451: Options Database Key:
1452: . -mesh_mover_type <type> - Sets the mesh mvoer type
1454: Level: intermediate
1456: .keywords: vector, set, type
1457: .seealso: MeshMoverCreate()
1458: @*/
1459: int MeshMoverSetType(MeshMover mover, MeshMoverType method) {
1460: int (*r)(MeshMover);
1461: PetscTruth match;
1462: int ierr;
1466: PetscTypeCompare((PetscObject) mover, method, &match);
1467: if (match == PETSC_TRUE) return(0);
1469: /* Get the function pointers for the vector requested */
1470: if (!MeshMoverRegisterAllCalled) {
1471: MeshMoverRegisterAll(PETSC_NULL);
1472: }
1473: PetscFListFind(mover->comm, MeshMoverList, method,(void (**)(void)) &r);
1474: if (!r) SETERRQ1(PETSC_ERR_ARG_WRONG, "Unknown mesh mover type: %s", method);
1476: if (mover->ops->destroy != PETSC_NULL) {
1477: (*mover->ops->destroy)(mover);
1478: }
1479: (*r)(mover);
1481: PetscObjectChangeTypeName((PetscObject) mover, method);
1482: return(0);
1483: }
1485: #undef __FUNCT__
1487: /*@C
1488: MeshMoverGetType - Gets the MeshMover method type name (as a string).
1490: Not collective
1492: Input Parameter:
1493: . mover - The mover
1495: Output Parameter:
1496: . type - The name of MeshMover method
1498: Level: intermediate
1500: .keywords: mover, get, type, name
1501: .seealso MeshMoverSetType()
1502: @*/
1503: int MeshMoverGetType(MeshMover mover, MeshMoverType *type)
1504: {
1510: if (!MeshMoverRegisterAllCalled) {
1511: MeshMoverRegisterAll(PETSC_NULL);
1512: }
1513: *type = mover->type_name;
1514: return(0);
1515: }
1517: /*@C
1518: MeshMoverRegister - Adds a creation method to the mesh movement package.
1520: Synopsis:
1522: MeshMoverRegister(char *name, char *path, char *func_name, int (*create_func)(MeshMover))
1524: Not Collective
1526: Input Parameters:
1527: + name - The name of a new user-defined creation routine
1528: . path - The path (either absolute or relative) of the library containing this routine
1529: . func_name - The name of the creation routine
1530: - create_func - The creation routine itself
1532: Notes:
1533: MeshMoverRegister() may be called multiple times to add several user-defined meshes.
1535: If dynamic libraries are used, then the fourth input argument (create_func) is ignored.
1537: Sample usage:
1538: .vb
1539: MeshMoverRegisterDynamic("my_mesh", "/home/username/my_lib/lib/libO/solaris/libmy.a", "MyMeshMoverCreate", MyMeshMoverCreate);
1540: .ve
1542: Then, your mesh type can be chosen with the procedural interface via
1543: .vb
1544: MeshMoverCreate(MPI_Comm, MeshMover *);
1545: MeshMoverSetType(vec, "my_mesh_mover")
1546: .ve
1547: or at runtime via the option
1548: .vb
1549: -mesh_mover_type my_mesh_mover
1550: .ve
1552: Note: $PETSC_ARCH and $BOPT occuring in pathname will be replaced with appropriate values.
1554: Level: advanced
1556: .keywords: MeshMover, register
1557: .seealso: MeshMoverRegisterAll(), MeshMoverRegisterDestroy()
1558: @*/
1559: #undef __FUNCT__
1561: int MeshMoverRegister(const char sname[], const char path[], const char name[], int (*function)(MeshMover))
1562: {
1563: char fullname[256];
1564: int ierr;
1567: PetscStrcpy(fullname, path);
1568: PetscStrcat(fullname, ":");
1569: PetscStrcat(fullname, name);
1570: PetscFListAdd(&MeshMoverList, sname, fullname, (void (*)(void)) function);
1571: return(0);
1572: }
1574: #undef __FUNCT__
1576: /*@C
1577: MeshMoverRegisterDestroy - Frees the list of creation routines for
1578: meshes that were registered by MeshMoverRegister().
1580: Not collective
1582: Level: advanced
1584: .keywords: mesh, register, destroy
1585: .seealso: MeshMoverRegister(), MeshMoverRegisterAll()
1586: @*/
1587: int MeshMoverRegisterDestroy() {
1591: if (MeshMoverList != PETSC_NULL) {
1592: PetscFListDestroy(&MeshMoverList);
1593: MeshMoverList = PETSC_NULL;
1594: }
1595: MeshMoverRegisterAllCalled = 0;
1596: return(0);
1597: }
1599: EXTERN_C_BEGIN
1600: extern int MeshMoverCreate_Triangular_2D(MeshMover);
1601: EXTERN_C_END
1603: #undef __FUNCT__
1605: /*@C
1606: MeshMoverRegisterAll - Registers all of the generation routines in the MeshMover package.
1608: Not Collective
1610: Input parameter:
1611: . path - The dynamic library path
1613: Level: advanced
1615: .keywords: MeshMover, register, all
1616: .seealso: MeshMoverRegister(), MeshMoverRegisterDestroy()
1617: @*/
1618: int MeshMoverRegisterAll(const char path[])
1619: {
1623: MeshMoverRegisterAllCalled = 1;
1625: MeshMoverRegister(MESH_MOVER_TRIANGULAR_2D, path, "MeshMoverCreate_Triangular_2D", MeshMoverCreate_Triangular_2D);
1626:
1627: return(0);
1628: }
1630: /*-------------------------------------------------- Grid Functions --------------------------------------------------*/
1631: #undef __FUNCT__
1633: /*@
1634: GridReformMesh - Reconstruct the grid from the current description
1635: and underlying mesh boundary. The old mesh is returned if requested.
1636: This allows postprocessing of the mesh before GridReform() is called.
1638: Collective on Grid
1640: Input Parameters:
1641: + grid - The grid
1642: . refine - Flag indicating whether to refine or recreate the mesh
1643: . area - [Optional] A function which gives area constraints for refinement after reform, PETSC_NULL for no refinement
1644: - newBd - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage
1646: Output Parameter:
1647: . oldMesh - The old mesh
1649: Level: advanced
1651: .keywords: Grid, mesh, reform, refine
1652: .seealso: MeshReform(), GridReform(), GridRefineMesh()
1653: @*/
1654: int GridReformMesh(Grid grid, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *oldMesh)
1655: {
1656: Mesh mesh;
1657: Mesh newMesh;
1658: int ierr;
1663: /* Reform mesh */
1664: GridGetMesh(grid, &mesh);
1665: MeshReform(mesh, refine, area, newBd, &newMesh);
1667: /* Save old information */
1668: if (oldMesh != PETSC_NULL){
1670: *oldMesh = mesh;
1671: }
1672: grid->mesh = newMesh;
1674: /* Recalculate grid structures for mesh */
1675: GridCalcPointBCNodes(grid);
1677: /* Recreate class structure */
1678: FieldClassMapDestroy(grid->cm);
1679: FieldClassMapCreateTriangular2D(grid, grid->numActiveFields, grid->defaultFields, &grid->cm);
1681: /* Specific setup */
1682: if (grid->ops->gridreformmesh != PETSC_NULL) {
1683: (*grid->ops->gridreformmesh)(grid);
1684: }
1685: return(0);
1686: }
1688: #undef __FUNCT__
1690: /*@
1691: GridReform - Reconstruct the variable structures, usually after a change
1692: in mesh structure. The optional arguments allow projection of the
1693: previous solution onto the new mesh.
1695: Collective on Grid
1697: Input Parameters:
1698: + grid - The grid
1699: . oldMesh - [Optional] The old mesh
1700: - x - [Optional] A grid vector to project onto the new mesh
1702: Output Parameter:
1703: . y - [Optional] The projected grid vector
1705: Level: advanced
1707: .keywords: grid, mesh, reform
1708: .seealso: GridReformMesh(), MeshReform()
1709: @*/
1710: int GridReform(Grid grid, Mesh oldMesh, GVec x, GVec *y)
1711: {
1712: MPI_Comm comm;
1713: VarOrdering order;
1714: FieldClassMap map;
1715: int numFields;
1716: int *fields;
1717: PetscTruth isConstrained, explicitConstraints;
1718: InterpolationContext iCtx; /* Used to form the rhs */
1719: /* L_2 projection variables */
1720: #if 0
1721: GMat M; /* The (local) mass matrix M */
1722: GMat m; /* The restriction of M to the boundary m */
1723: GMat S; /* The Schur complement S = m^T M^{-1} m */
1724: GVec sol, rhs; /* The solution x and rhs f */
1725: GVec projRhs; /* The projection rhs g */
1726: GVec uzawaRhs; /* The Uzawa rhs m^T M^{-1} f - g */
1727: GVec uzawaSol; /* The Uzawa solution lambda = (m^T M^{-1} m)^{-1} (m^T M^{-1} f - g) */
1728: KSP ksp; /* The Krylov solver for M */
1729: PC pc; /* The preconditioner for M */
1730: SLES sles; /* The solver for M */
1731: SLES uzawaSles; /* The solver for S */
1732: PetscScalar minusOne = -1.0;
1733: PetscScalar zero = 0.0;
1734: int its;
1735: #endif
1736: /* Local Interpolation variables */
1737: VarOrdering reduceOrder, reduceOrderOld; /* TEMPORARY, until I fix the interpolation */
1738: GVec bdReduceVec, prevBdReduceVec; /* TEMPORARY, until I fix the interpolation */
1739: GVec bdReduceVecOld, prevBdReduceVecOld; /* TEMPORARY, until I fix the interpolation */
1740: GVec bdReduceVecDiff, prevBdReduceVecDiff; /* TEMPORARY, until I fix the interpolation */
1741: PetscReal norm;
1742: int maxComp;
1743: int field, f;
1744: int ierr;
1748: PetscLogEventBegin(GRID_Reform, grid, 0, 0, 0);
1749: /* Make sure the vector is saved in local storage */
1750: GridGlobalToLocal(grid, INSERT_VALUES, x);
1751: /* Need to setup grid again */
1752: grid->setupcalled = PETSC_FALSE;
1753: grid->bdSetupCalled = PETSC_FALSE;
1754: /* Check arguments */
1755: if ((x == PETSC_NULL) || (y == PETSC_NULL)) {
1756: GridSetupGhostScatter(grid);
1757: return(0);
1758: }
1762: PetscObjectGetComm((PetscObject) oldMesh, &comm);
1763: GridIsConstrained(grid, &isConstrained);
1764: GridGetExplicitConstraints(grid, &explicitConstraints);
1765: if (explicitConstraints == PETSC_FALSE) {
1766: order = grid->order;
1767: } else {
1768: order = grid->constraintOrder;
1769: }
1770: VarOrderingGetClassMap(order, &map);
1771: #ifdef INTERFACE_UPDATE
1772: FieldClassMapGetNumFields(map, &numFields);
1773: FieldClassMapGetFields(map, &fields);
1774: #endif
1775: numFields = map->numFields;
1776: fields = map->fields;
1778: /* Project fields onto new mesh */
1779: switch(grid->interpolationType)
1780: {
1781: #if 0
1782: case INTERPOLATION_L2:
1783: /* Interpolation Procedure:
1784: We must locally interpolate the value of the old field at each Gaussian intergration
1785: point in an element of the new field. The function PointFunctionInterpolateField()
1786: calls GVecInterpolateField(), which calls DiscretizationInterpolateField(). The vector, field,
1787: and old mesh are stored in a context so that PointFunctionInterpolateField() can pass
1788: them to GVecInterpolateField(). The indices and values for DiscretizationInterpolateField() are
1789: calculated by GridCalcLocalElementVecIndices() and GridLocalToElement(). These
1790: functions need the old mesh, classes, variable ordering, and local vector ghostVec in
1791: order to correctly calculate the old vector indices and values. However, the
1792: GVecEvaluateFunctionGalerkin() uses the new values to get the weak form. The old mesh
1793: we have, and the classes are stored in classesOld[]. We will keep the variable
1794: ordering and local vector in temporary veriables. The old copies must be set to NULL
1795: or they will be deleted.
1796: */
1798: /* Save old structures */
1799: iCtx.vec = x;
1800: iCtx.mesh = oldMesh;
1801: iCtx.order = order;
1802: iCtx.ghostVec = grid->ghostVec;
1804: /* Setup context */
1805: for(f = 0, maxComp = 1; f < numFields; f++) {
1806: maxComp = PetscMax(maxComp, grid->comp[fields[f]]);
1807: }
1808: MPI_Comm_size(comm, &iCtx.numProcs);
1809: MPI_Comm_rank(comm, &iCtx.rank);
1810: iCtx.batchMode = PETSC_FALSE;
1811: iCtx.curPoint = 0;
1812: iCtx.maxPoints = iCtx.numProcs;
1813: PetscMalloc((iCtx.numProcs+1) * sizeof(int), &iCtx.numPoints);
1814: PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.activeProcs);
1815: PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.calcProcs);
1816: PetscMalloc(iCtx.maxPoints*3 * sizeof(PetscScalar), &iCtx.coords);
1817: PetscMalloc(iCtx.numProcs*maxComp * sizeof(PetscScalar), &iCtx.values);
1819: /* Create new variable numbering and ghost structures */
1820: if (explicitConstraints == PETSC_TRUE) {
1821: grid->constraintOrder = PETSC_NULL;
1822: } else {
1823: grid->order = PETSC_NULL;
1824: }
1825: grid->ghostVec = PETSC_NULL;
1826: reduceOrderOld = grid->reduceOrder;
1827: prevBdReduceVec = grid->bdReduceVec;
1828: prevBdReduceVecOld = grid->bdReduceVecOld;
1829: prevBdReduceVecDiff = grid->bdReduceVecDiff;
1830: grid->reduceOrder = PETSC_NULL;
1831: grid->bdReduceVec = PETSC_NULL;
1832: grid->bdReduceVecOld = PETSC_NULL;
1833: grid->bdReduceVecDiff = PETSC_NULL;
1834: GridSetupGhostScatter(grid);
1835: reduceOrder = grid->reduceOrder;
1836: bdReduceVec = grid->bdReduceVec;
1837: bdReduceVecOld = grid->bdReduceVecOld;
1838: bdReduceVecDiff = grid->bdReduceVecDiff;
1839: grid->reduceOrder = reduceOrderOld;
1840: grid->bdReduceVec = prevBdReduceVec;
1841: grid->bdReduceVecOld = prevBdReduceVecOld;
1842: grid->bdReduceVecDiff = prevBdReduceVecDiff;
1843: if (grid->bdReduceVecCur == bdReduceVec) {
1844: grid->bdReduceVecCur = grid->bdReduceVec;
1845: } else if (grid->bdReduceVecCur == bdReduceVecOld) {
1846: grid->bdReduceVecCur = grid->bdReduceVecOld;
1847: } else if (grid->bdReduceVecCur == bdReduceVecDiff) {
1848: grid->bdReduceVecCur = grid->bdReduceVecDiff;
1849: } else {
1850: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
1851: }
1853: /* Allocate structures */
1854: if (explicitConstraints) {
1855: GVecCreateConstrained(grid, y);
1856: } else {
1857: GVecCreate(grid, y);
1858: }
1859: GVecDuplicate(*y, &rhs);
1860: GVecDuplicate(*y, &sol);
1861: GVecCreateBoundaryRestriction(grid, &uzawaSol);
1862: GVecDuplicate(uzawaSol, &projRhs);
1863: GVecDuplicate(uzawaSol, &uzawaRhs);
1865: /* Create the rhs */
1866: for(f = 0; f < numFields; f++) {
1867: field = fields[f];
1868: iCtx.field = field;
1869: if (iCtx.batchMode == PETSC_FALSE) {
1870: GVecEvaluateFunctionGalerkinCollective(rhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1871:
1872: #if 1
1873: GVecEvaluateBoundaryFunctionGalerkinCollective(projRhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1874:
1875: #else
1876: VecSet(&zero, projRhs);
1877: #endif
1878: } else {
1879: GVecEvaluateFunctionGalerkin(rhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1880: GVecInterpolateFieldBatch(rhs, field, &iCtx);
1881: iCtx.curPoint = iCtx.numPoints[iCtx.rank];
1882: GVecEvaluateFunctionGalerkin(rhs, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx);
1883:
1884: if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing");
1885: GVecEvaluateBoundaryFunctionGalerkin(projRhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1886:
1887: GVecInterpolateFieldBatch(projRhs, field, &iCtx);
1888: iCtx.curPoint = iCtx.numPoints[iCtx.rank];
1889: GVecEvaluateFunctionGalerkin(projRhs, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx);
1890:
1891: if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing");
1892: }
1893: }
1894: if (grid->bdReduceVecCur == grid->bdReduceVec) {
1895: grid->bdReduceVecCur = bdReduceVec;
1896: } else if (grid->bdReduceVecCur == grid->bdReduceVecOld) {
1897: grid->bdReduceVecCur = bdReduceVecOld;
1898: } else if (grid->bdReduceVecCur == grid->bdReduceVecDiff) {
1899: grid->bdReduceVecCur = bdReduceVecDiff;
1900: } else {
1901: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
1902: }
1903: grid->reduceOrder = reduceOrder;
1904: grid->bdReduceVec = bdReduceVec;
1905: grid->bdReduceVecOld = bdReduceVecOld;
1906: grid->bdReduceVecDiff = bdReduceVecDiff;
1907: if (grid->reduceSystem == PETSC_TRUE) {
1908: VarOrderingDestroy(reduceOrderOld);
1909: VecDestroy(prevBdReduceVec);
1910: VecDestroy(prevBdReduceVecOld);
1911: VecDestroy(prevBdReduceVecDiff);
1912: }
1914: /* Destroy old structures */
1915: VarOrderingDestroy(iCtx.order);
1916: VecDestroy(iCtx.ghostVec);
1918: /* Cleanup */
1919: PetscFree(iCtx.numPoints);
1920: PetscFree(iCtx.activeProcs);
1921: PetscFree(iCtx.calcProcs);
1922: PetscFree(iCtx.values);
1923: PetscFree(iCtx.coords);
1925: /* Create the inner solver */
1926: SLESCreate(PETSC_COMM_SELF, &sles);
1927: SLESAppendOptionsPrefix(sles, "grid_proj_");
1928: SLESGetKSP(sles, &ksp);
1929: SLESGetPC(sles, &pc);
1930: KSPSetType(ksp, KSPCG);
1931: PCSetType(pc, PCJACOBI);
1932: SLESSetFromOptions(sles);
1933: /* Create the outer solver */
1934: SLESCreate(PETSC_COMM_SELF, &uzawaSles);
1935: SLESAppendOptionsPrefix(uzawaSles, "grid_proj_uzawa_");
1936: SLESGetKSP(uzawaSles, &ksp);
1937: SLESGetPC(uzawaSles, &pc);
1938: KSPSetType(ksp, KSPGMRES);
1939: PCSetType(pc, PCNONE);
1940: SLESSetFromOptions(uzawaSles);
1942: /* Create M, the mass matrix */
1943: GMatCreate(grid, &M);
1944: GMatEvaluateOperatorGalerkin(M, PETSC_NULL, numFields, fields, fields, IDENTITY, 1.0, MAT_FINAL_ASSEMBLY, PETSC_NULL);
1945:
1946: MatNorm(M, NORM_FROBENIUS, &norm);
1947: PetscPrintf(PETSC_COMM_SELF, "Mass matrix norm %lfn", norm);
1949: /* Create m, the restriction of M to the boundary */
1950: GMatCreateBoundaryRestriction(grid, &m);
1951: #if 1
1952: GMatEvaluateBoundaryOperatorGalerkin(m, numFields, grid->cm->fields, grid->bdOrder, grid->bdLocOrder,
1953: fields, grid->order, grid->locOrder, IDENTITY, 1.0, MAT_FINAL_ASSEMBLY,
1954: PETSC_NULL);
1955:
1956: #else
1957: MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY);
1958: MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY);
1959: MatZeroEntries(m);
1960: #endif
1961: MatNorm(m, NORM_FROBENIUS, &norm);
1962: PetscPrintf(PETSC_COMM_SELF, "Restricted mass matrix norm %lfn", norm);
1964: VecNorm(rhs, NORM_2, &norm);
1965: PetscPrintf(PETSC_COMM_SELF, "Rhs norm %lfn", norm);
1966: VecNorm(projRhs, NORM_2, &norm);
1967: PetscPrintf(PETSC_COMM_SELF, "Restricted rhs norm %lfn", norm);
1969: /* Create the Schur complement */
1970: SLESSetOperators(sles, M, M, SAME_NONZERO_PATTERN);
1971: GMatCreateUzawa(M, m, sles, &S);
1973: /* Create the restriction of the rhs */
1974: SLESSolve(sles, rhs, sol, &its);
1975: MatMultTranspose(m, sol, uzawaRhs);
1976: VecAXPY(&minusOne, projRhs, uzawaRhs);
1977: VecNorm(uzawaRhs, NORM_2, &norm);
1978: PetscPrintf(PETSC_COMM_SELF, "Uzawa rhs norm %lfn", norm);
1980: /* Solve system */
1981: SLESSetOperators(uzawaSles, S, S, SAME_NONZERO_PATTERN);
1982: SLESSolve(uzawaSles, uzawaRhs, uzawaSol, &its);
1983: PetscPrintf(PETSC_COMM_SELF, "Uzawa solution took %d iterationsn", its);
1984: VecNorm(uzawaSol, NORM_2, &norm);
1985: PetscPrintf(PETSC_COMM_SELF, "Uzawa solution norm %lfn", norm);
1987: /* Get solution */
1988: MatMult(m, uzawaSol, sol);
1989: VecNorm(sol, NORM_2, &norm);
1990: PetscPrintf(PETSC_COMM_SELF, "Boundary solution update norm %lfn", norm);
1991: VecAYPX(&minusOne, rhs, sol);
1992: VecNorm(sol, NORM_2, &norm);
1993: PetscPrintf(PETSC_COMM_SELF, "Solution update norm %lfn", norm);
1994: SLESSolve(sles, sol, *y, &its);
1995: PetscPrintf(PETSC_COMM_SELF, "Projection solution took %d iterationsn", its);
1997: VecNorm(*y, NORM_2, &norm);
1998: PetscPrintf(PETSC_COMM_SELF, "Projected solution norm %lfn", norm);
2000: /* Cleanup */
2001: SLESDestroy(sles);
2002: SLESDestroy(uzawaSles);
2003: GMatDestroyUzawa(S);
2004: GMatDestroy(M);
2005: GMatDestroy(m);
2006: GVecDestroy(rhs);
2007: GVecDestroy(sol);
2008: GVecDestroy(projRhs);
2009: GVecDestroy(uzawaRhs);
2010: GVecDestroy(uzawaSol);
2012: break;
2013: #endif
2014: case INTERPOLATION_LOCAL:
2015: /* Interpolation Procedure:
2016: We must locally interpolate the value of the old field at each node in an element of
2017: the new field. The function PointFunctionInterpolateField()calls GVecInterpolateField(),
2018: which calls DiscretizationInterpolateField(). The vector, field, and old mesh are stored in a
2019: context so that PointFunctionInterpolateField() can pass them to GVecInterpolateField().
2020: The indices and values for DiscretizationInterpolateField() are calculated by
2021: GridCalcLocalElementVecIndices() and GridLocalToElement(). These functions need the old
2022: mesh, classes, variable ordering, and local vector ghostVec in order to correctly calculate
2023: the old vector indices and values. The old mesh we have, and the classes are stored in
2024: classesOld[]. We will keep the variable ordering and local vector in temporary veriables.
2025: The old copies must be set to NULL or they will be deleted.
2026: */
2028: /* Save old structures */
2029: iCtx.vec = x;
2030: iCtx.mesh = oldMesh;
2031: iCtx.order = order;
2032: iCtx.ghostVec = grid->ghostVec;
2034: /* Setup context */
2035: for(f = 0, maxComp = 1; f < numFields; f++) {
2036: maxComp = PetscMax(maxComp, grid->fields[fields[f]].numComp);
2037: }
2038: MPI_Comm_size(comm, &iCtx.numProcs);
2039: MPI_Comm_rank(comm, &iCtx.rank);
2040: iCtx.batchMode = PETSC_FALSE;
2041: iCtx.curPoint = 0;
2042: iCtx.maxPoints = iCtx.numProcs;
2043: PetscMalloc((iCtx.numProcs+1) * sizeof(int), &iCtx.numPoints);
2044: PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.activeProcs);
2045: PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.calcProcs);
2046: PetscMalloc(iCtx.maxPoints*3 * sizeof(PetscScalar), &iCtx.coords);
2047: PetscMalloc(iCtx.numProcs*maxComp * sizeof(PetscScalar), &iCtx.values);
2048: PetscMemzero(iCtx.numPoints, (iCtx.numProcs+1) * sizeof(int));
2049: PetscOptionsHasName(PETSC_NULL, "-grid_int_mode", &iCtx.batchMode);
2051: /* Create new variable numbering and ghost structures */
2052: if (explicitConstraints == PETSC_TRUE) {
2053: grid->constraintOrder = PETSC_NULL;
2054: } else {
2055: grid->order = PETSC_NULL;
2056: }
2057: grid->ghostVec = PETSC_NULL;
2058: reduceOrderOld = grid->reduceOrder;
2059: prevBdReduceVec = grid->bdReduceVec;
2060: prevBdReduceVecOld = grid->bdReduceVecOld;
2061: prevBdReduceVecDiff = grid->bdReduceVecDiff;
2062: grid->reduceOrder = PETSC_NULL;
2063: grid->bdReduceVec = PETSC_NULL;
2064: grid->bdReduceVecOld = PETSC_NULL;
2065: grid->bdReduceVecDiff = PETSC_NULL;
2066: GridSetupGhostScatter(grid);
2067: reduceOrder = grid->reduceOrder;
2068: bdReduceVec = grid->bdReduceVec;
2069: bdReduceVecOld = grid->bdReduceVecOld;
2070: bdReduceVecDiff = grid->bdReduceVecDiff;
2071: grid->reduceOrder = reduceOrderOld;
2072: grid->bdReduceVec = prevBdReduceVec;
2073: grid->bdReduceVecOld = prevBdReduceVecOld;
2074: grid->bdReduceVecDiff = prevBdReduceVecDiff;
2075: if (grid->bdReduceVecCur == bdReduceVec) {
2076: grid->bdReduceVecCur = grid->bdReduceVec;
2077: } else if (grid->bdReduceVecCur == bdReduceVecOld) {
2078: grid->bdReduceVecCur = grid->bdReduceVecOld;
2079: } else if (grid->bdReduceVecCur == bdReduceVecDiff) {
2080: grid->bdReduceVecCur = grid->bdReduceVecDiff;
2081: } else {
2082: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
2083: }
2085: /* Allocate structures */
2086: if (explicitConstraints) {
2087: GVecCreateConstrained(grid, y);
2088: } else {
2089: GVecCreate(grid, y);
2090: }
2092: /* Locally interpolate the solution values */
2093: for(f = 0; f < numFields; f++) {
2094: field = fields[f];
2095: iCtx.field = field;
2096: if (iCtx.batchMode == PETSC_FALSE) {
2097: GVecEvaluateFunctionCollective(*y, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
2098: } else {
2099: GVecEvaluateFunction(*y, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
2100: GVecInterpolateFieldBatch(*y, field, &iCtx);
2101: iCtx.curPoint = iCtx.numPoints[iCtx.rank];
2102: GVecEvaluateFunction(*y, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx);
2103: if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1])
2104: SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing");
2105: }
2106: }
2107: if (grid->bdReduceVecCur == grid->bdReduceVec) {
2108: grid->bdReduceVecCur = bdReduceVec;
2109: } else if (grid->bdReduceVecCur == grid->bdReduceVecOld) {
2110: grid->bdReduceVecCur = bdReduceVecOld;
2111: } else if (grid->bdReduceVecCur == grid->bdReduceVecDiff) {
2112: grid->bdReduceVecCur = bdReduceVecDiff;
2113: } else {
2114: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
2115: }
2116: grid->reduceOrder = reduceOrder;
2117: grid->bdReduceVec = bdReduceVec;
2118: grid->bdReduceVecOld = bdReduceVecOld;
2119: grid->bdReduceVecDiff = bdReduceVecDiff;
2120: if (grid->reduceSystem == PETSC_TRUE) {
2121: VarOrderingDestroy(reduceOrderOld);
2122: VecDestroy(prevBdReduceVec);
2123: VecDestroy(prevBdReduceVecOld);
2124: VecDestroy(prevBdReduceVecDiff);
2125: }
2127: VecNorm(*y, NORM_2, &norm);
2128: PetscPrintf(PETSC_COMM_SELF, "Projected solution norm %gn", norm);
2130: /* Destroy old structures */
2131: VarOrderingDestroy(iCtx.order);
2132: VecDestroy(iCtx.ghostVec);
2134: /* Cleanup */
2135: PetscFree(iCtx.numPoints);
2136: PetscFree(iCtx.activeProcs);
2137: PetscFree(iCtx.calcProcs);
2138: PetscFree(iCtx.values);
2139: PetscFree(iCtx.coords);
2140: break;
2141: case INTERPOLATION_HALO:
2142: case INTERPOLATION_L2_LOCAL:
2143: SETERRQ1(PETSC_ERR_SUP, "Interpolation type %d not yet supported", grid->interpolationType);
2144: default:
2145: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid interpolation type %d", grid->interpolationType);
2146: }
2147: /* Anything specific to the implementation */
2148: if (grid->ops->gridreform != PETSC_NULL) {
2149: (*grid->ops->gridreform)(grid, oldMesh, x, y);
2150: }
2151: PetscLogEventEnd(GRID_Reform, grid, 0, 0, 0);
2152: return(0);
2153: }