Actual source code: grid.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: grid.c,v 1.24 2000/07/16 05:40:26 knepley Exp $";
3: #endif
5: /*
6: This file provides routines for grid vectors (vectors that are associated
7: with grids, possibly with multiple degrees of freedom per node).
9: This component is new; changes in the user interface will be occuring in the
10: near future!
11: */
13: #include "petscsles.h" /* For ALE Mesh Support *//*I "petscsles.h" I*/
14: #include "src/grid/gridimpl.h" /*I "grid.h" I*//*I "gvec.h" I*/
16: /* Logging support */
17: int GRID_COOKIE;
18: int GRID_Reform, GRID_SetUp;
20: /*---------------------------------------------- Standard Grid Functions --------------------------------------------*/
21: #undef __FUNCT__
23: /*@
24: GridDestroy - Destroys a grid.
26: Collective on Grid
28: Input Parameter:
29: . grid - The grid
31: Level: beginner
33: .keywrods: grid, destroy
34: .seealso: GridView(), GridCreateGVec(), GridSetUp(), GridRefine()
35: @*/
36: int GridDestroy(Grid grid)
37: {
42: if (--grid->refct > 0) return(0);
43: (*grid->ops->destroy)(grid);
44: PetscLogObjectDestroy(grid);
45: PetscHeaderDestroy(grid);
46: return(0);
47: }
49: #undef __FUNCT__
51: /*@
52: GridView - Views a grid.
54: Collective on Grid
56: Input Parameters:
57: + grid - The grid
58: - viewer - The viewer, PETSC_VIEWER_STDOUT_WORLD by default
60: Level: beginner
62: .keywords: grid, view
63: .seealso: GVecView()
64: @*/
65: int GridView(Grid grid, PetscViewer viewer)
66: {
71: /* Do not check return from setup, since we may not have a problem yet */
72: if (grid->setupcalled == PETSC_FALSE) GridSetUp(grid);
73: if (!viewer) {
74: viewer = PETSC_VIEWER_STDOUT_SELF;
75: } else {
77: }
78: (*grid->ops->view)(grid, viewer);
79: return(0);
80: }
82: EXTERN_C_BEGIN
83: #undef __FUNCT__
85: int GridSerialize_Generic(MPI_Comm comm, Grid *grid, PetscViewer viewer, PetscTruth store) {
86: Grid g;
87: Mesh mesh;
88: int fd, hasParent;
89: int one = 1;
90: int zero = 0;
91: int field, len;
92: char *type_name = PETSC_NULL;
93: PetscTruth match;
94: int ierr;
100: PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match);
101: if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer");
102: PetscViewerBinaryGetDescriptor(viewer, &fd);
104: if (store) {
105: g = *grid;
107: MeshSerialize(comm, &g->mesh, viewer, store);
108: PetscStrlen(g->type_name, &len);
109: PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0);
110: PetscBinaryWrite(fd, g->type_name, len, PETSC_CHAR, 0);
112: PetscBinaryWrite(fd, &g->numFields, 1, PETSC_INT, 0);
113: for(field = 0; field < g->numFields; field++) {
114: if (g->fields[field].name != PETSC_NULL) {
115: PetscStrlen(g->fields[field].name, &len);
116: } else {
117: len = 0;
118: }
119: PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0);
120: PetscBinaryWrite(fd, g->fields[field].name, len, PETSC_CHAR, 0);
121: PetscBinaryWrite(fd, &g->fields[field].numComp, 1, PETSC_INT, 0);
122: PetscStrlen(g->fields[field].discType, &len);
123: PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0);
124: PetscBinaryWrite(fd, g->fields[field].discType, len, PETSC_CHAR, 0);
125: DiscretizationSerialize(g->comm, &g->fields[field].disc, viewer, store);
126: PetscBinaryWrite(fd, &g->fields[field].isActive, 1, PETSC_INT, 0);
127: PetscBinaryWrite(fd, &g->fields[field].isConstrained, 1, PETSC_INT, 0);
128: PetscBinaryWrite(fd, &g->fields[field].constraintCompDiff, 1, PETSC_INT, 0);
129: }
131: if (g->gridparent != PETSC_NULL) {
132: PetscBinaryWrite(fd, &one, 1, PETSC_INT, 0);
133: GridSerialize(comm, &g->gridparent, viewer, store);
134: } else {
135: PetscBinaryWrite(fd, &zero, 1, PETSC_INT, 0);
136: }
137: PetscBinaryWrite(fd, &g->isConstrained, 1, PETSC_INT, 0);
138: #ifdef NEW_GRID_SERIALIZE
139: /* This stuff is all made in GridSetUp */
140: FieldClassSerialize(g->comm, &g->cm, viewer, store);
141: VarOrderingSerialize(g->cm, &g->order, viewer, store);
142: FieldClassSerialize(g->comm, &g->constraintCM, viewer, store);
143: VarOrderingSerialize(g->constraintCM, &g->constraintOrder, viewer, store);
144: ISSerialize(g->comm, &g->constraintOrdering, viewer, store);
145: MatSerialize(g->comm, &g->constraintMatrix, viewer, store);
146: #endif
148: PetscBinaryWrite(fd, &g->reduceAlpha, 1, PETSC_SCALAR, 0);
149: PetscBinaryWrite(fd, &g->interpolationType, 1, PETSC_INT, 0);
150: PetscBinaryWrite(fd, &g->viewField, 1, PETSC_INT, 0);
151: PetscBinaryWrite(fd, &g->viewComp, 1, PETSC_INT, 0);
152: } else {
153: MeshSerialize(comm, &mesh, viewer, store);
154: GridCreate(mesh, &g);
155: MeshDestroy(mesh);
156: PetscBinaryRead(fd, &len, 1, PETSC_INT);
157: if (len) {
158: PetscMalloc((len+1) * sizeof(char), &type_name);
159: type_name[len] = 0;
160: }
161: PetscBinaryRead(fd, type_name, len, PETSC_CHAR);
162: PetscBinaryRead(fd, &g->numFields, 1, PETSC_INT);
163: g->maxFields = g->numFields;
164: PetscFree(g->fields);
165: PetscMalloc(g->numFields * sizeof(Field), &g->fields);
166: for(field = 0; field < g->numFields; field++) {
167: PetscBinaryRead(fd, &len, 1, PETSC_INT);
168: if (len) {
169: PetscMalloc((len+1) * sizeof(char), &g->fields[field].name);
170: g->fields[field].name[len] = 0;
171: }
172: PetscBinaryRead(fd, g->fields[field].name, len, PETSC_CHAR);
173: PetscBinaryRead(fd, &g->fields[field].numComp, 1, PETSC_INT);
174: PetscBinaryRead(fd, &len, 1, PETSC_INT);
175: PetscMalloc((len+1) * sizeof(DiscretizationType), &g->fields[field].discType);
176: g->fields[field].discType[len] = 0;
177: PetscBinaryRead(fd, g->fields[field].discType, len, PETSC_CHAR);
178: DiscretizationSerialize(g->comm, &g->fields[field].disc, viewer, store);
179: PetscTypeCompare((PetscObject) g->fields[field].disc, g->fields[field].discType, &match);
180: if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Incorrect serialization of Discretization");
181: PetscBinaryRead(fd, &g->fields[field].isActive, 1, PETSC_INT);
182: PetscBinaryRead(fd, &g->fields[field].isConstrained, 1, PETSC_INT);
183: PetscBinaryRead(fd, &g->fields[field].constraintCompDiff, 1, PETSC_INT);
184: }
185: GridSetType(g, type_name);
186: if (len) {
187: PetscFree(type_name);
188: }
190: PetscBinaryRead(fd, &hasParent, 1, PETSC_INT);
191: if (hasParent) {
192: GridSerialize(comm, &g->gridparent, viewer, store);
193: }
194: PetscBinaryRead(fd, &g->isConstrained, 1, PETSC_INT);
195: #ifdef NEW_GRID_SERIALIZE
196: /* This stuff is all made in GridSetUp */
197: FieldClassSerialize(g->comm, &g->cm, viewer, store);
198: VarOrderingSerialize(g->cm, &g->order, viewer, store);
199: FieldClassSerialize(g->comm, &g->constraintCM, viewer, store);
200: VarOrderingSerialize(g->constraintCM, &g->constraintOrder, viewer, store);
201: ISSerialize(g->comm, &g->constraintOrdering, viewer, store);
202: MatSerialize(g->comm, &g->constraintMatrix, viewer, store);
203: #endif
205: PetscBinaryRead(fd, &g->reduceAlpha, 1, PETSC_SCALAR);
206: PetscBinaryRead(fd, &g->interpolationType, 1, PETSC_INT);
207: PetscBinaryRead(fd, &g->viewField, 1, PETSC_INT);
208: PetscBinaryRead(fd, &g->viewComp, 1, PETSC_INT);
210: *grid = g;
211: }
213: return(0);
214: }
215: EXTERN_C_END
217: #undef __FUNCT__
219: /*@
220: GridSetUp - Set up any required internal data structures for the grid.
222: Collective on Grid
224: Input Parameter:
225: . grid - The grid
227: Level: intermediate
229: .keywords: grid, setup
230: .seealso: GridDestroy(), GridCreateGVec()
231: @*/
232: int GridSetUp(Grid grid)
233: {
234: LocalVarOrdering reduceLocOrder;
235: LocalVarOrdering locOrder;
236: int sField, tField;
237: int bd, bc, op;
238: int ierr;
242: if (grid->setupcalled == PETSC_TRUE) return(0);
243: grid->setupcalled = PETSC_TRUE;
245: PetscLogEventBegin(GRID_SetUp, grid, 0, 0, 0);
246: if (grid->numFields <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "At least one field must be defined on the Grid");
247: /* Setup boundary sizes */
248: for(bd = 0; bd < grid->numBd; bd++) {
249: if (grid->bdSize[bd] != PETSC_NULL) {
250: PetscFree(grid->bdSize[bd]);
251: }
252: PetscMalloc(grid->numFields * sizeof(int), &grid->bdSize[bd]);
253: PetscLogObjectMemory(grid, grid->numFields * sizeof(int));
254: }
255: /* Setup type-specific stuff */
256: if (grid->ops->gridsetup) {
257: (*grid->ops->gridsetup)(grid);
258: if (ierr != 0) PetscFunctionReturn(ierr);
259: }
260: /* Setup constraint structures */
261: (*grid->ops->gridsetupconstraints)(grid, grid->constraintCtx);
262: /* Setup ghost variable structures */
263: GridSetupGhostScatter(grid);
264: /* Setup boundary conditions reductions */
265: if (grid->reduceSystem == PETSC_TRUE) {
266: /* Create storage for reduction of Rhs */
267: if (grid->bdReduceVec != PETSC_NULL) {
268: GVecDestroy(grid->bdReduceVec);
269: }
270: GVecCreateRectangularGhost(grid, grid->reduceOrder, &grid->bdReduceVec);
271: grid->bdReduceVecCur = grid->bdReduceVec;
273: /* Hopefully, the user specified the right context by now */
274: GridCalcBCValues(grid, PETSC_FALSE, grid->reduceContext);
276: if (grid->reduceElement == PETSC_FALSE) {
277: /* Create storage for explicit reduction of Rhs */
278: GVecCreate(grid, &grid->reduceVec);
279: GMatCreateRectangular(grid, grid->reduceOrder, grid->order, &grid->bdReduceMat);
280: MatSetOption(grid->bdReduceMat, MAT_NEW_NONZERO_ALLOCATION_ERR);
282: /* Initialize storage */
283: MatZeroEntries(grid->bdReduceMat);
285: /* Evaluate the elimination block */
286: for(bc = 0; bc < grid->numBC; bc++) {
287: for(op = 0; op < grid->numMatOps; op++) {
288: sField = grid->matOps[op].field;
289: tField = grid->fields[sField].disc->operators[grid->matOps[op].op]->test->field;
290: if (sField == grid->bc[bc].field) {
291: LocalVarOrderingCreate(grid, 1, &tField, &locOrder);
292: LocalVarOrderingCreate(grid, 1, &sField, &reduceLocOrder);
293: GMatEvaluateOperatorGalerkin(grid->bdReduceMat, PETSC_NULL, 1, &sField, &tField, grid->matOps[op].op,
294: grid->matOps[op].alpha, MAT_FLUSH_ASSEMBLY, grid->reduceContext);
295:
296: LocalVarOrderingDestroy(locOrder);
297: LocalVarOrderingDestroy(reduceLocOrder);
298: }
299: #ifdef PETSC_USE_BOPT_g
300: PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
301: #endif
302: }
303: }
304: for(bc = 0; bc < grid->numPointBC; bc++) {
305: for(op = 0; op < grid->numMatOps; op++) {
306: sField = grid->matOps[op].field;
307: tField = grid->fields[sField].disc->operators[grid->matOps[op].op]->test->field;
308: if (sField == grid->pointBC[bc].field) {
309: LocalVarOrderingCreate(grid, 1, &tField, &locOrder);
310: LocalVarOrderingCreate(grid, 1, &sField, &reduceLocOrder);
311: GMatEvaluateOperatorGalerkin(grid->bdReduceMat, PETSC_NULL, 1, &sField, &tField, grid->matOps[op].op,
312: grid->matOps[op].alpha, MAT_FLUSH_ASSEMBLY, grid->reduceContext);
313:
314: LocalVarOrderingDestroy(locOrder);
315: LocalVarOrderingDestroy(reduceLocOrder);
316: }
317: #ifdef PETSC_USE_BOPT_g
318: PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
319: #endif
320: }
321: }
322: MatAssemblyBegin(grid->bdReduceMat, MAT_FINAL_ASSEMBLY);
323: MatAssemblyEnd(grid->bdReduceMat, MAT_FINAL_ASSEMBLY);
324: }
325: }
326: PetscLogEventEnd(GRID_SetUp, grid, 0, 0, 0);
327: return(0);
328: }
330: #undef __FUNCT__
332: /*@
333: GridSetupGhostScatter - Set up the scatter from a global vector to local values, including ghosts.
335: Collective on Grid
337: Input Parameter:
338: . grid - The grid
340: Level: advanced
342: .keywords: grid, ghost, scatter
343: .seealso: GridDestroy(), GridCreateGVec()
344: @*/
345: int GridSetupGhostScatter(Grid grid)
346: {
347: PetscTruth isConstrained, explicitConstraints;
348: int ierr;
352: GridSetUp(grid);
353: /* Destroy old structures */
354: if (grid->ghostVec != PETSC_NULL) {
355: VecDestroy(grid->ghostVec);
356: }
357: if (grid->ghostScatter != PETSC_NULL) {
358: VecScatterDestroy(grid->ghostScatter);
359: }
360: /* Setup scatter from global vector to local ghosted vector */
361: GridIsConstrained(grid, &isConstrained);
362: GridGetExplicitConstraints(grid, &explicitConstraints);
363: if ((isConstrained == PETSC_TRUE) && (explicitConstraints == PETSC_TRUE)) {
364: (*grid->ops->gridsetupghostscatter)(grid, grid->constraintOrder, &grid->ghostVec, &grid->ghostScatter);
365:
366: } else {
367: (*grid->ops->gridsetupghostscatter)(grid, grid->order, &grid->ghostVec, &grid->ghostScatter);
368: }
369: return(0);
370: }
372: #undef __FUNCT__
374: /*@
375: GridSetupBoundary - Set up any required internal data structures for the boundary.
377: Collective on Grid
379: Input Parameter:
380: . grid - The grid
382: Level: advanced
384: .keywords: grid, setup, boundary
385: .seealso: GridDestroy(), GridCreateGVec()
386: @*/
387: int GridSetupBoundary(Grid grid)
388: {
393: if (grid->bdSetupCalled == PETSC_TRUE) return(0);
394: GridSetUp(grid);
395: if (grid->ops->gridsetupboundary) {
396: (*grid->ops->gridsetupboundary)(grid);
397: }
398: return(0);
399: }
401: #undef __FUNCT__
403: /*@
404: GridGetMesh - Returns the context for the underlying mesh.
406: Not collective
408: Input Parameter:
409: . grid - The initial grid
411: Output Parameter:
412: . mesh - The mesh
414: Level: intermediate
416: .keywords: grid, mesh
417: .seealso: GridDestroy()
418: @*/
419: int GridGetMesh(Grid grid, Mesh *mesh)
420: {
424: *mesh = grid->mesh;
425: return(0);
426: }
428: #undef __FUNCT__
430: /*@
431: GridGetClassMap - This function returns the default field class map for the grid.
433: Not collective
435: Input Parameter:
436: . grid - The grid
438: Output Parameter:
439: . order - The default field class map
441: Level: intermediate
443: .keywords: Grid, get, class map
444: .seealso: GridGetMesh(), GridGetOrder(), FieldClassMapCreate()
445: @*/
446: int GridGetClassMap(Grid grid, FieldClassMap *map)
447: {
451: if (grid->isConstrained == PETSC_TRUE) {
452: *map = grid->constraintCM;
453: } else {
454: *map = grid->cm;
455: }
456: return(0);
457: }
459: #undef __FUNCT__
461: /*@
462: GridGetOrder - This function returns the default ordering for the grid.
464: Not collective
466: Input Parameter:
467: . grid - The grid
469: Output Parameter:
470: . order - The default variable ordering
472: Level: intermediate
474: .keywords: Grid, get, order
475: .seealso: GridGetLocalOrder(), VarOrderingCreate()
476: @*/
477: int GridGetOrder(Grid grid, VarOrdering *order)
478: {
482: if (grid->isConstrained == PETSC_TRUE) {
483: *order = grid->constraintOrder;
484: } else {
485: *order = grid->order;
486: }
487: return(0);
488: }
490: #undef __FUNCT__
492: /*@
493: GridGetLocalOrder - This function returns the local ordering for the grid.
495: Not collective
497: Input Parameter:
498: . grid - The grid
500: Output Parameter:
501: . order - The local variable ordering
503: Level: intermediate
505: .keywords: Grid, get, local, order
506: .seealso: GridGetOrder(), LocalVarOrderingCreate()
507: @*/
508: int GridGetLocalOrder(Grid grid, LocalVarOrdering *order)
509: {
513: *order = grid->locOrder;
514: return(0);
515: }
517: #undef __FUNCT__
519: /*@
520: GridGetDiscretization - This function returns the discretization for the given field.
522: Not collective
524: Input Parameters:
525: + grid - The grid
526: - field - The field
528: Output Parameter:
529: . disc - The Discretization
531: Level: intermediate
533: .keywords: Grid, get, Discretization, field
534: .seealso: GridGetOrder(), DiscretizationCreate()
535: @*/
536: int GridGetDiscretization(Grid grid, int field, Discretization *disc) {
540: GridValidField(grid, field);
541: *disc = grid->fields[field].disc;
542: return(0);
543: }
545: #undef __FUNCT__
547: /*
548: GridSetTypeFromOptions - Sets the type of grid from user options.
550: Collective on Grid
552: Input Parameter:
553: . grid - The grid
555: Level: intermediate
557: .keywords: Grid, set, options, database, type
558: .seealso: GridSetFromOptions(), GridSetType()
559: */
560: static int GridSetTypeFromOptions(Grid grid)
561: {
562: PetscTruth opt;
563: char *defaultType;
564: char typeName[256];
565: int dim = grid->dim;
566: int ierr;
569: if (grid->type_name != PETSC_NULL) {
570: defaultType = grid->type_name;
571: } else {
572: switch(dim) {
573: case 1:
574: defaultType = GRID_TRIANGULAR_1D;
575: break;
576: case 2:
577: defaultType = GRID_TRIANGULAR_2D;
578: break;
579: default:
580: SETERRQ1(PETSC_ERR_SUP, "Grid dimension %d is not supported", dim);
581: }
582: }
584: if (GridRegisterAllCalled == PETSC_FALSE) {
585: GridRegisterAll(PETSC_NULL);
586: }
587: PetscOptionsList("-grid_type", "Grid method"," GridSetType", GridList, defaultType, typeName, 256, &opt);
588: if (opt == PETSC_TRUE) {
589: GridSetType(grid, typeName);
590: } else {
591: GridSetType(grid, defaultType);
592: }
593: return(0);
594: }
596: #undef __FUNCT__
598: /*
599: GridSetSerializeTypeFromOptions - Sets the type of grid serialization from user options.
601: Collective on Grid
603: Input Parameter:
604: . grid - The grid
606: Level: intermediate
608: .keywords: Grid, set, options, database, type
609: .seealso: GridSetFromOptions(), GridSetSerializeType()
610: */
611: static int GridSetSerializeTypeFromOptions(Grid grid)
612: {
613: PetscTruth opt;
614: char *defaultType;
615: char typeName[256];
616: int ierr;
619: if (grid->serialize_name != PETSC_NULL) {
620: defaultType = grid->serialize_name;
621: } else {
622: defaultType = GRID_SER_GENERIC;
623: }
625: if (GridSerializeRegisterAllCalled == PETSC_FALSE) {
626: GridSerializeRegisterAll(PETSC_NULL);
627: }
628: PetscOptionsList("-grid_serialize_type", "Grid serialization method"," GridSetSerializeType", GridSerializeList,
629: defaultType, typeName, 256, &opt);
630: if (opt == PETSC_TRUE) {
631: GridSetSerializeType(grid, typeName);
632: } else {
633: GridSetSerializeType(grid, defaultType);
634: }
635: return(0);
636: }
638: #undef __FUNCT__
640: /*@
641: GridSetFromOptions - Sets various Grid parameters from user options.
643: Collective on Grid
645: Input Parameter:
646: . grid - The grid
648: Notes: To see all options, run your program with the -help option, or consult the users manual.
649: Must be called after GridCreate() but before the Grid is used.
651: Level: intermediate
653: .keywords: Grid, set, options, database
654: .seealso: GridCreate(), GridPrintHelp(), GridSetOptionsPrefix()
655: @*/
656: int GridSetFromOptions(Grid grid)
657: {
658: Mesh mesh;
659: char *intNames[NUM_INTERPOLATIONS] = {"Local", "Halo", "L2", "L2Local"};
660: char intType[256];
661: PetscTruth opt;
662: int ierr;
666: PetscOptionsBegin(grid->comm, grid->prefix, "Grid options", "Grid");
668: /* Handle generic grid options */
669: PetscOptionsHasName(PETSC_NULL, "-help", &opt);
670: if (opt == PETSC_TRUE) {
671: GridPrintHelp(grid);
672: }
674: /* Constraint options */
675: PetscOptionsName("-grid_explicit_constraints", "Explicitly form constrained system", "GridGetExplicitConstraints", &grid->explicitConstraints);
677: /* Interpolation options */
678: PetscOptionsEList("-grid_int_type", "Interpolation Method", "None", intNames, NUM_INTERPOLATIONS,
679: intNames[grid->interpolationType], intType, 256, &opt);
680: if (opt == PETSC_TRUE) {
681: PetscTruth islocal, ishalo, isl2, isl2local;
683: PetscStrcasecmp(intType, intNames[INTERPOLATION_LOCAL], &islocal);
684: PetscStrcasecmp(intType, intNames[INTERPOLATION_HALO], &ishalo);
685: PetscStrcasecmp(intType, intNames[INTERPOLATION_L2], &isl2);
686: PetscStrcasecmp(intType, intNames[INTERPOLATION_L2_LOCAL], &isl2local);
687: if (islocal == PETSC_TRUE) {
688: grid->interpolationType = INTERPOLATION_LOCAL;
689: } else if (ishalo == PETSC_TRUE) {
690: grid->interpolationType = INTERPOLATION_HALO;
691: } else if (isl2 == PETSC_TRUE) {
692: grid->interpolationType = INTERPOLATION_L2;
693: } else if (isl2local == PETSC_TRUE) {
694: grid->interpolationType = INTERPOLATION_L2_LOCAL;
695: }
696: }
698: /* Graphics options */
699: PetscOptionsInt("-grid_view_field", "View Field", "None", grid->viewField, &grid->viewField, &opt);
700: PetscOptionsInt("-grid_view_comp", "View Component", "None", grid->viewComp, &grid->viewComp, &opt);
702: /* Handle grid type options */
703: GridSetTypeFromOptions(grid);
705: /* Handle grid serialization options */
706: GridSetSerializeTypeFromOptions(grid);
708: /* Handle specific grid options */
709: if (grid->ops->setfromoptions != PETSC_NULL) {
710: (*grid->ops->setfromoptions)(grid);
711: }
713: PetscOptionsEnd();
715: /* Handle subobject options */
716: GridGetMesh(grid, &mesh);
717: MeshSetFromOptions(mesh);
719: GridViewFromOptions(grid, grid->name);
720: return(0);
721: }
723: #undef __FUNCT__
725: /*@
726: GridViewFromOptions - This function visualizes the grid based upon user options.
728: Collective on Grid
730: Input Parameter:
731: . grid - The grid
733: Level: intermediate
735: .keywords: Grid, view, options, database
736: .seealso: GridSetFromOptions(), GridView()
737: @*/
738: int GridViewFromOptions(Grid grid, char *title)
739: {
740: PetscViewer viewer;
741: PetscDraw draw;
742: PetscTruth opt;
743: char *titleStr;
744: char typeName[1024];
745: char fileName[PETSC_MAX_PATH_LEN];
746: int len;
747: int ierr;
750: PetscOptionsHasName(grid->prefix, "-grid_view", &opt);
751: if (opt == PETSC_TRUE) {
752: PetscOptionsGetString(grid->prefix, "-grid_view", typeName, 1024, &opt);
753: PetscStrlen(typeName, &len);
754: if (len > 0) {
755: PetscViewerCreate(grid->comm, &viewer);
756: PetscViewerSetType(viewer, typeName);
757: PetscOptionsGetString(grid->prefix, "-grid_view_file", fileName, 1024, &opt);
758: if (opt == PETSC_TRUE) {
759: PetscViewerSetFilename(viewer, fileName);
760: } else {
761: PetscViewerSetFilename(viewer, grid->name);
762: }
763: GridView(grid, viewer);
764: PetscViewerFlush(viewer);
765: PetscViewerDestroy(viewer);
766: } else {
767: GridView(grid, PETSC_NULL);
768: }
769: }
770: PetscOptionsHasName(grid->prefix, "-grid_view_draw", &opt);
771: if (opt == PETSC_TRUE) {
772: PetscViewerDrawOpen(grid->comm, 0, 0, 0, 0, 300, 300, &viewer);
773: PetscViewerDrawGetDraw(viewer, 0, &draw);
774: if (title != PETSC_NULL) {
775: titleStr = title;
776: } else {
777: PetscObjectName((PetscObject) grid); CHKERRQ(ierr) ;
778: titleStr = grid->name;
779: }
780: PetscDrawSetTitle(draw, titleStr);
781: GridView(grid, viewer);
782: PetscViewerFlush(viewer);
783: PetscDrawPause(draw);
784: PetscViewerDestroy(viewer);
785: }
786: return(0);
787: }
789: #undef __FUNCT__
791: /*@
792: GridPrintHelp - Prints all options for the grid.
794: Input Parameter:
795: . grid - The grid
797: Options Database Keys:
798: $ -help, -h
800: Level: intermediate
802: .keywords: Grid, help
803: .seealso: GridSetFromOptions()
804: @*/
805: int GridPrintHelp(Grid grid)
806: {
807: char p[64];
808: int ierr;
813: PetscStrcpy(p, "-");
814: if (grid->prefix != PETSC_NULL) {
815: PetscStrcat(p, grid->prefix);
816: }
818: (*PetscHelpPrintf)(grid->comm, "Grid options ------------------------------------------------n");
819: (*PetscHelpPrintf)(grid->comm," %sgrid_type <typename> : Sets the grid typen", p);
820: (*PetscHelpPrintf)(grid->comm," %sgrid_view_field <field>: Sets the field to visualizen", p);
821: (*PetscHelpPrintf)(grid->comm," %sgrid_view_comp <field>: Sets the component to visualizen", p);
822: return(0);
823: }
825: #undef __FUNCT__
827: /*@C
828: GridSetOptionsPrefix - Sets the prefix used for searching for all grid options in the database.
830: Input Parameters:
831: + grid - The grid
832: - prefix - The prefix to prepend to all option names
834: Notes:
835: A hyphen (-) must NOT be given at the beginning of the prefix name.
836: The first character of all runtime options is AUTOMATICALLY the hyphen.
838: Level: intermediate
840: .keywords: grid, set, options, prefix, database
841: .seealso: GridSetFromOptions()
842: @*/
843: int GridSetOptionsPrefix(Grid grid, char *prefix)
844: {
849: PetscObjectSetOptionsPrefix((PetscObject) grid, prefix);
850: return(0);
851: }
853: #undef __FUNCT__
855: /*@C
856: GridAppendOptionsPrefix - Appends to the prefix used for searching for all grid options in the database.
858: Collective on Grid
860: Input Parameters:
861: + grid - The grid
862: - prefix - The prefix to prepend to all option names
864: Notes:
865: A hyphen (-) must NOT be given at the beginning of the prefix name.
866: The first character of all runtime options is AUTOMATICALLY the hyphen.
868: Level: intermediate
870: .keywords: grid, append, options, prefix, database
871: .seealso: GridGetOptionsPrefix()
872: @*/
873: int GridAppendOptionsPrefix(Grid grid, char *prefix)
874: {
876:
879: PetscObjectAppendOptionsPrefix((PetscObject) grid, prefix);
880: return(0);
881: }
883: #undef __FUNCT__
885: /*@
886: GridGetOptionsPrefix - Sets the prefix used for searching for all grid options in the database.
888: Input Parameter:
889: . grid - The grid
891: Output Parameter:
892: . prefix - A pointer to the prefix string used
894: Level: intermediate
896: .keywords: grid, get, options, prefix, database
897: .seealso: GridAppendOptionsPrefix()
898: @*/
899: int GridGetOptionsPrefix(Grid grid, char **prefix)
900: {
905: PetscObjectGetOptionsPrefix((PetscObject) grid, prefix);
906: return(0);
907: }
909: /*----------------------------------------------- Graphics Functions ------------------------------------------------*/
910: #undef __FUNCT__
912: /*@C GridSetViewField
913: This function sets the field component to be used in visualization.
915: Collective on Grid
917: Input Parameters:
918: + grid - The grid
919: . field - The field
920: - comp - The component
922: Level: intermediate
924: .keywords grid, view, field
925: .seealso GridGetMeshInfo
926: @*/
927: int GridSetViewField(Grid grid, int field, int comp)
928: {
931: GridValidField(grid, field);
932: if ((comp < 0) || (comp >= grid->fields[field].numComp)) {
933: SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid field component %d should be in [0,%d)", comp, grid->fields[field].numComp);
934: }
935: grid->viewField = field;
936: grid->viewComp = comp;
937: return(0);
938: }
940: /*------------------------------------------ Interface to Mesh Functions --------------------------------------------*/
941: #undef __FUNCT__
943: /*@C
944: GridGetNodeClass - This function returns the class of a given node
946: Not collective
948: Input Parameters:
949: + grid - The grid
950: - node - The node
952: Output Parameter:
953: . nclass - The node class
955: Level: intermediate
957: .keywords grid, mesh, node, class
958: .seealso GridGetMeshInfo
959: @*/
960: int GridGetNodeClass(Grid grid, int node, int *nclass)
961: {
962: FieldClassMap map;
963: int ierr;
968: GridGetClassMap(grid, &map);
969: FieldClassMapGetNodeClass(map, node, nclass);
970: return(0);
971: }
973: #undef __FUNCT__
975: /*@
976: GridGetNearestNode - This function returns the node nearest to the given point on which
977: the field is defined, or -1 if the node is not contained in the local mesh.
979: Not collective
981: Input Parameters:
982: + mesh - The mesh
983: . field - The field
984: - x,y,z - The node coordinates
986: Output Parameter:
987: . node - The nearest node
989: Level: beginner
991: .keywords: grid, node, point location
992: .seealso MeshNearestNode(), MeshLocatePoint()
993: @*/
994: int GridGetNearestNode(Grid grid, int field, double x, double y, double z, int *node)
995: {
996: FieldClassMap map;
997: Mesh mesh;
998: Partition p;
999: PetscTruth defined;
1000: double minDist, dist, nx, ny, nz;
1001: int minNode, globalMinNode, allMinNode;
1002: int numCorners, startNode, endNode;
1003: int elem, corner, nearNode, nclass;
1004: int ierr;
1009: FieldClassMapCreateTriangular2D(grid, 1, &field, &map);
1010: GridGetMesh(grid, &mesh);
1011: MeshGetNumCorners(mesh, &numCorners);
1012: MeshGetPartition(mesh, &p);
1013: PartitionGetStartNode(p, &startNode);
1014: PartitionGetEndNode(p, &endNode);
1015: MeshLocatePoint(mesh, x, y, z, &elem);
1016: if (elem >= 0) {
1017: /* Find first node with field defined */
1018: for(corner = 0, minDist = PETSC_MAX; corner < numCorners; corner++) {
1019: MeshGetNodeFromElement(mesh, elem, corner, &minNode);
1020: FieldClassMapGetNodeClass(map, minNode, &nclass);
1021: FieldClassMapIsDefined(map, field, nclass, &defined);
1022: if (defined == PETSC_TRUE) {
1023: MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz);
1024: minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y));
1025: break;
1026: }
1027: }
1028: if (corner == numCorners) SETERRQ2(PETSC_ERR_ARG_WRONG, "Element %d has no node with field %d defined", elem, field);
1029: /* Find closest node with field defined */
1030: for(corner = 1; corner < numCorners; corner++) {
1031: MeshGetNodeFromElement(mesh, elem, corner, &nearNode);
1032: MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz);
1033: FieldClassMapGetNodeClass(map, nearNode, &nclass);
1034: FieldClassMapIsDefined(map, field, nclass, &defined);
1035: dist = PetscSqr(nx - x) + PetscSqr(ny - y);
1036: if ((defined == PETSC_TRUE) && (dist < minDist)) {
1037: minDist = dist;
1038: minNode = nearNode;
1039: }
1040: }
1041: PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
1042: } else {
1043: minNode = -1;
1044: globalMinNode = -1;
1045: }
1046: FieldClassMapDestroy(map);
1048: /* The minimum node might still be a ghost node */
1049: MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, grid->comm);
1050: if ((allMinNode >= startNode) && (allMinNode < endNode)) {
1051: *node = allMinNode - startNode;
1052: } else {
1053: *node = -1;
1054: }
1055: if (allMinNode < 0) PetscFunctionReturn(1);
1056: return(0);
1057: }
1059: #undef __FUNCT__
1061: /*@
1062: GridGetNearestBdNode - This function returns the boundary node nearest to the given point
1063: on which the field is defined, or -1 if the node is not contained in the local mesh.
1065: Not collective
1067: Input Parameters:
1068: + mesh - The mesh
1069: . field - The field
1070: - x,y,z - The node coordinates
1072: Output Parameter:
1073: . node - The nearest boundary node
1075: Level: beginner
1077: .keywords: grid, node, point location
1078: .seealso MeshNearestBdNode(), MeshLocatePoint()
1079: @*/
1080: int GridGetNearestBdNode(Grid grid, int field, double x, double y, double z, int *node)
1081: {
1082: FieldClassMap map;
1083: Mesh mesh;
1084: Partition p;
1085: PetscTruth defined;
1086: double minDist, dist, nx, ny, nz;
1087: int minNode, globalMinNode, allMinNode;
1088: int numCorners, startNode, endNode;
1089: int elem, corner, nearNode, nclass, bd;
1090: int ierr;
1095: FieldClassMapCreateTriangular2D(grid, 1, &field, &map);
1096: GridGetMesh(grid, &mesh);
1097: MeshGetNumCorners(mesh, &numCorners);
1098: MeshGetPartition(mesh, &p);
1099: PartitionGetStartNode(p, &startNode);
1100: PartitionGetEndNode(p, &endNode);
1101: MeshLocatePoint(mesh, x, y, z, &elem);
1102: if (elem >= 0) {
1103: /* Find first node with field defined */
1104: for(corner = 0, minDist = PETSC_MAX; corner < numCorners; corner++) {
1105: MeshGetNodeFromElement(mesh, elem, corner, &minNode);
1106: MeshGetNodeBoundary(mesh, minNode, &bd);
1107: FieldClassMapGetNodeClass(map, minNode, &nclass);
1108: FieldClassMapIsDefined(map, field, nclass, &defined);
1109: if ((bd != 0) && (defined == PETSC_TRUE)) {
1110: MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz);
1111: minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y));
1112: break;
1113: }
1114: }
1115: if (corner == numCorners) SETERRQ2(PETSC_ERR_ARG_WRONG, "Element %d has no node with field %d defined", elem, field);
1116: /* Find closest node with field defined */
1117: for(corner = 1; corner < numCorners; corner++) {
1118: MeshGetNodeFromElement(mesh, elem, corner, &nearNode);
1119: MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz);
1120: MeshGetNodeBoundary(mesh, nearNode, &bd);
1121: FieldClassMapGetNodeClass(map, nearNode, &nclass);
1122: FieldClassMapIsDefined(map, field, nclass, &defined);
1123: dist = PetscSqr(nx - x) + PetscSqr(ny - y);
1124: if ((bd != 0) && (defined == PETSC_TRUE) && (dist < minDist)) {
1125: minDist = dist;
1126: minNode = nearNode;
1127: }
1128: }
1129: PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
1130: } else {
1131: minNode = -1;
1132: globalMinNode = -1;
1133: }
1134: FieldClassMapDestroy(map);
1136: /* The minimum node might still be a ghost node */
1137: MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, grid->comm);
1138: if ((allMinNode >= startNode) && (allMinNode < endNode)) {
1139: *node = allMinNode - startNode;
1140: } else {
1141: *node = -1;
1142: }
1143: if (allMinNode < 0)
1144: PetscFunctionReturn(1);
1145: return(0);
1146: }
1148: #undef __FUNCT__
1150: /*@
1151: GridGetBoundarySize - Retrieves the size of the specified boundary for
1152: a given field.
1154: Not collective
1156: Input Parameters:
1157: + grid - The grid
1158: . bd - The boundary marker
1159: - field - The field
1161: Output Parameter:
1162: . size - The size of the boundary
1164: Level: intermediate
1166: .keywords: grid, boundary, size
1167: .seealso: GridGetBoundaryNext(), GridGetBoundaryStart()
1168: @*/
1169: int GridGetBoundarySize(Grid grid, int bd, int field, int *size)
1170: {
1171: int bdIndex;
1177: GridValidField(grid, field);
1178: if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first")
1179: ierr = MeshGetBoundaryIndex(grid->mesh, bd, &bdIndex);
1180: *size = grid->bdSize[bdIndex][field];
1181: return(0);
1182: }
1184: #undef __FUNCT__
1186: /*@
1187: GridGetBoundaryStart - Retrieves the canonical node number of the first node
1188: on the specified boundary in a class contained in a given field.
1190: Not collective
1192: Input Parameters:
1193: + grid - The grid
1194: . bd - The boundary marker
1195: . field - The field
1196: - ghost - Flag for including ghost nodes
1198: Output Parameters:
1199: + node - The canonical node number
1200: - nclass - The node class
1202: Level: intermediate
1204: .keywords: grid, boundary, start
1205: .seealso: GridGetBoundaryNext(), MeshGetBoundaryStart
1206: @*/
1207: int GridGetBoundaryStart(Grid grid, int bd, int field, PetscTruth ghost, int *node, int *nclass)
1208: {
1209: int f;
1216: if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first");
1217: /* Look for field in default class map */
1218: for(f = 0; f < grid->cm->numFields; f++) {
1219: if (grid->cm->fields[f] == field) break;
1220: }
1221: if (f == grid->cm->numFields) {
1222: SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Field %d not present in the default class map", field);
1223: }
1224: (*grid->ops->getboundarystart)(grid, bd, f, ghost, grid->cm, node, nclass);
1225: return(0);
1226: }
1228: #undef __FUNCT__
1230: /*@
1231: GridGetBoundaryNext - Retrieves the canonical node number of the next node
1232: on the specified boundary in a class contained in a given field, with -1
1233: indicating boundary termination.
1235: Not collective
1237: Input Parameters:
1238: + grid - The grid
1239: . bd - The boundary marker
1240: . field - The field
1241: - ghost - Flag for including ghost nodes
1243: Output Parameters:
1244: + node - The canonical node number
1245: - nclass - The node class
1247: Level: intermediate
1249: .keywords: grid, boundary, next
1250: .seealso: GridGetBoundaryStart(), MeshGetBoundaryNext
1251: @*/
1252: int GridGetBoundaryNext(Grid grid, int bd, int field, PetscTruth ghost, int *node, int *nclass)
1253: {
1254: int f;
1261: if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first");
1262: /* Look for field in default class map */
1263: for(f = 0; f < grid->cm->numFields; f++) {
1264: if (grid->cm->fields[f] == field) break;
1265: }
1266: if (f == grid->cm->numFields) {
1267: SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Field %d not present in the default class map", field);
1268: }
1269: (*grid->ops->getboundarynext)(grid, bd, field, ghost, grid->cm, node, nclass);
1270: return(0);
1271: }
1273: #undef __FUNCT__
1275: /*@
1276: GridRefineMesh - Refine a grid based on area constraints.
1278: Collective on Grid
1280: Input Parameters:
1281: + grid - The initial grid
1282: - area - A function which gives an area constraint when evaluated inside an element
1284: Output Parameter:
1285: . newgrid - The refined grid
1287: Level: advanced
1289: .keywords: grid, refinement
1290: .seealso: MeshRefine(), GridReformMesh()
1291: @*/
1292: int GridRefineMesh(Grid grid, PointFunction area, Grid *newgrid)
1293: {
1294: Mesh newMesh;
1295: int field;
1296: int ierr;
1301: MeshRefine(grid->mesh, area, &newMesh);
1302: GridCreate(newMesh, newgrid);
1303: MeshDestroy(newMesh);
1304: for(field = 0; field < grid->numFields; field++) {
1305: GridAddField(*newgrid, grid->fields[field].name, grid->fields[field].discType, grid->fields[field].numComp, 0);
1306: }
1307: GridSetType(*newgrid, grid->type_name);
1308: /* Copy registered operators in discretizations */
1309: /* Set parent: (*newgrid)->gridparent = grid; */
1311: return(0);
1312: }
1314: /*------------------------------------------ Distributed Data Functions ---------------------------------------------*/
1315: #undef __FUNCT__
1317: /*@
1318: GridGlobalToLocalGeneral - Scatters values including ghost variables
1319: from the given global vector into the given ghost vector.
1321: Collective on Grid
1323: Input Parameters:
1324: + grid - The grid
1325: . vec - The global vector
1326: . mode - The insertion mode, either INSERT_VALUES or ADD_VALUES
1327: - scatter - The scatter
1329: Output Parameters:
1330: . ghostVec - The local vector
1332: Level: intermediate
1334: .keywords: grid, scatter, global, local
1335: .seealso: GridGlobalToLocal(), GridLocalToGlobal(), GVecGlobalToLocal()
1336: @*/
1337: int GridGlobalToLocalGeneral(Grid grid, GVec vec, Vec ghostVec, InsertMode mode, VecScatter scatter)
1338: {
1342: if (mode == ADD_VALUES) SETERRQ(PETSC_ERR_SUP, "No support for adding ghost values");
1343: VecScatterBegin(vec, ghostVec, mode, SCATTER_FORWARD, scatter);
1344: VecScatterEnd(vec, ghostVec, mode, SCATTER_FORWARD, scatter);
1345: return(0);
1346: }
1348: #undef __FUNCT__
1350: /*@
1351: GridGlobalToLocal - Scatters values including ghost variables
1352: from the given global vector into the local grid ghost vector.
1354: Collective on Grid
1356: Input Parameters:
1357: + grid - The grid
1358: . mode - The insertion mode, either INSERT_VALUES or ADD_VALUES
1359: - vec - The grid vector
1361: Level: intermediate
1363: .keywords: grid, scatter, global, local
1364: .seealso: GridGlobalToLocalGeneral(), GridLocalToGlobal(), GVecGlobalToLocal()
1365: @*/
1366: int GridGlobalToLocal(Grid grid, InsertMode mode, GVec vec)
1367: {
1368: Grid g;
1369: int ierr;
1375: GVecGetGrid(vec, &g);
1376: if (grid != g) SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector does not arise from this grid");
1377: GridGlobalToLocalGeneral(grid, vec, grid->ghostVec, mode, grid->ghostScatter);
1378: return(0);
1379: }
1381: #undef __FUNCT__
1383: /*@
1384: GridLocalToGlobal - Scatters local processor values including ghost variables
1385: from the local grid ghost vector into the given global vector.
1387: Collective on Grid
1389: Input Parameters:
1390: + grid - The grid
1391: . mode - The insertion mode, either SET_VALUES or ADD_VALUES
1392: - vec - The grid vector
1394: Level: intermediate
1396: .seealso: GridGlobalToLocal(), GVecLocalToGlobal()
1397: @*/
1398: int GridLocalToGlobal(Grid grid, InsertMode mode, GVec vec)
1399: {
1400: Grid g;
1401: int ierr;
1407: GVecGetGrid(vec, &g);
1408: if (grid != g) SETERRQ(PETSC_ERR_ARG_WRONG, "Vector does not arise from this grid");
1409: VecScatterBegin(grid->ghostVec, vec, mode, SCATTER_REVERSE, grid->ghostScatter);
1410: VecScatterEnd(grid->ghostVec, vec, mode, SCATTER_REVERSE, grid->ghostScatter);
1411: return(0);
1412: }
1414: /*---------------------------------------------- InterGrid Functions ------------------------------------------------*/
1415: #undef __FUNCT__
1417: /*@
1418: GridInterpolateElementVec - Interpolates a given element vector from one discretization to another.
1420: Not collective
1422: Input Parameters:
1423: + grid - The original grid
1424: . field - The original field
1425: . vec - The original element vector
1426: . newGrid - The grid defining the new element vector
1427: - newField - The new field
1429: Output Parameter:
1430: . newVec - The interpolated element vector
1432: Level: advanced
1434: .keywords: grid, interpolation
1435: .seealso: ElementVecCreate()
1436: @*/
1437: int GridInterpolateElementVec(Grid grid, int field, ElementVec vec, Grid newGrid, int newField, ElementVec newVec)
1438: {
1446: GridValidField(grid, field);
1447: GridValidField(newGrid, newField);
1448: if (grid->fields[field].numComp != newGrid->fields[newField].numComp) {
1449: SETERRQ2(PETSC_ERR_ARG_INCOMP, "Fields have a different number of components %d != %d",
1450: grid->fields[field].numComp, newGrid->fields[newField].numComp);
1451: }
1452: DiscretizationInterpolateElementVec(grid->fields[field].disc, vec, newGrid->fields[newField].disc, newVec);
1453: return(0);
1454: }
1456: #undef __FUNCT__
1458: /*@
1459: GridCreateRestriction - Generates restriction matrix between two grids.
1461: Input Parameters:
1462: + vf - The fine grid
1463: - vc - The coarse grid
1465: Output Parameter:
1466: . gmat - A sparse matrix containing restriction
1468: Level: advanced
1470: @*/
1471: int GridCreateRestriction(Grid vf, Grid vc, GMat *gmat)
1472: {
1480: if (vf->setupcalled == PETSC_FALSE) {
1481: GridSetUp(vf);
1482: }
1483: if (vc->setupcalled == PETSC_FALSE) {
1484: GridSetUp(vc);
1485: }
1486: if (!vf->ops->gridcreaterestriction) SETERRQ(PETSC_ERR_SUP, " ");
1487: (*vf->ops->gridcreaterestriction)(vf, vc, gmat);
1488: return(0);
1489: }
1491: /*---------------------------------------------- Assembly Functions -------------------------------------------------*/
1492: #undef __FUNCT__
1494: /*@
1495: GridLocalToElementGeneral - Scatters values including ghost variables from the given ghost vector
1496: into the given element vector.
1498: Not collective
1500: Input Parameters:
1501: + grid - The grid
1502: . ghostVec - The vector of values (including ghsot points)
1503: . reduceVec - The vector of boundary values
1504: . reduceSystem - The flag for reducing boundary conditions
1505: . reduceElement - The flag for putting boundary values in the element vector
1506: - vec - The element vector
1508: WARNING:
1509: Make sure that the indices in the element vector are local indices.
1511: Note:
1512: If reduceSystem and reduceElement are PTESC_TRUE, then boundary values are
1513: placed in vec. If reduceElement is PETSC_FALSE, then zero is used instead.
1515: Level: advanced
1517: .keywords: grid, element, scatter
1518: .seealso: GridLocalToElement(), GridGlobalToLocal()
1519: @*/
1520: int GridLocalToElementGeneral(Grid grid, Vec ghostVec, Vec reduceVec, PetscTruth reduceSystem, PetscTruth reduceElement, ElementVec vec)
1521: {
1522: PetscScalar *array, *reduceArray;
1523: PetscScalar alpha = -grid->reduceAlpha;
1524: int elemSize = vec->reduceSize;
1525: int numOverlapVars, numReduceOverlapVars;
1526: int var;
1527: int ierr;
1530: VecGetSize(ghostVec, &numOverlapVars);
1531: VecGetArray(ghostVec, &array);
1532: if (reduceSystem == PETSC_TRUE) {
1533: VecGetSize(reduceVec, &numReduceOverlapVars);
1534: VecGetArray(reduceVec, &reduceArray);
1535: for(var = 0; var < elemSize; var++) {
1536: #ifdef PETSC_USE_BOPT_g
1537: if (vec->indices[var] >= numOverlapVars) {
1538: SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)",
1539: vec->indices[var], numOverlapVars);
1540: }
1541: #endif
1542: if (vec->indices[var] < 0) {
1543: if (reduceElement == PETSC_TRUE) {
1544: #ifdef PETSC_USE_BOPT_g
1545: if (-(vec->indices[var]+1) >= numReduceOverlapVars) {
1546: SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)",
1547: -(vec->indices[var]+1), numReduceOverlapVars);
1548: }
1549: #endif
1550: vec->array[var] = alpha*reduceArray[-(vec->indices[var]+1)];
1551: } else {
1552: vec->array[var] = 0.0;
1553: }
1554: } else {
1555: vec->array[var] = array[vec->indices[var]];
1556: }
1557: }
1558: VecRestoreArray(reduceVec, &reduceArray);
1559: } else {
1560: for(var = 0; var < elemSize; var++) {
1561: #ifdef PETSC_USE_BOPT_g
1562: if ((vec->indices[var] < 0) || (vec->indices[var] >= numOverlapVars)) {
1563: SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)",
1564: vec->indices[var], numOverlapVars);
1565: }
1566: #endif
1567: vec->array[var] = array[vec->indices[var]];
1568: }
1569: }
1570: VecRestoreArray(ghostVec, &array);
1571: return(0);
1572: }
1574: #undef __FUNCT__
1576: /*@
1577: GridLocalToElement - Scatters values including ghost variables from the local grid ghost vector
1578: into the given element vector.
1580: Not collective
1582: Input Parameters:
1583: + grid - The grid
1584: - vec - The element vector
1586: WARNING:
1587: Make sure that the indices in the element vector are local indices.
1589: Level: advanced
1591: .keywords: grid, element, scatter
1592: .seealso: GridLocalToElementGeneral(), GridGlobalToLocal()
1593: @*/
1594: int GridLocalToElement(Grid grid, ElementVec vec)
1595: {
1601: GridLocalToElementGeneral(grid, grid->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec);
1602: return(0);
1603: }
1605: /*-------------------------------------------- Evaluation Functions -------------------------------------------------*/
1606: #undef __FUNCT__
1608: /*@C
1609: GridEvaluateRhs - This function constructs the weak form of the functions and
1610: operators specified with GridSetRhsFunction(), GridSetRhsOperator(),
1611: and GridSetRhsNonlinearOperator().
1613: Collective on Grid
1615: Input Parameters:
1616: + grid - The grid
1617: . x - The current iterate
1618: - ctx - The optional user context
1620: Output Parameters
1621: . f - The function value
1623: Level: beginner
1625: .keywords: grid, rhs
1626: .seealso: GridEvaluateSystemMatrix(), GridSetRhsFunction(), GridSetRhsOperator(), GridSetRhsNonlinearOperator()
1627: @*/
1628: int GridEvaluateRhs(Grid grid, GVec x, GVec f, PetscObject ctx)
1629: {
1630: Mesh mesh;
1631: MeshMover mover;
1632: Grid meshVelGrid;
1633: PetscScalar alpha = -grid->reduceAlpha;
1634: int ierr;
1638: /* x can be PETSC_NULL */
1640: /* Setup mesh velocity calculation for ALE variables */
1641: if (grid->ALEActive == PETSC_TRUE) {
1642: GridGetMesh(grid, &mesh);
1643: MeshGetMover(mesh, &mover);
1644: MeshMoverGetVelocityGrid(mover, &meshVelGrid);
1645: GridSetUp(meshVelGrid);
1646: }
1647: /* Calculate Rhs */
1648: (*grid->ops->gridevaluaterhs)(grid, x, f, ctx);
1649: /* Add boundary values */
1650: if ((grid->reduceSystem == PETSC_TRUE) && (grid->reduceElement == PETSC_FALSE) && (grid->bdReduceMat != PETSC_NULL)) {
1651: /* Calculate contribution of constrained variables to the Rhs */
1652: MatMult(grid->bdReduceMat, grid->bdReduceVecCur, grid->reduceVec);
1653: /* Update Rhs */
1654: VecAXPY(&alpha, grid->reduceVec, f);
1655: }
1656: return(0);
1657: }
1659: #undef __FUNCT__
1661: /*@C GridEvaluateRhsFunction
1662: This function constructs the weak form of the functions specified with GridSetRhsFunction().
1664: Collective on Grid
1666: Input Parameters:
1667: + grid - The grid
1668: . x - The current iterate
1669: - ctx - The optional user context
1671: Output Parameter:
1672: . f - The function value
1674: Level: beginner
1676: .keywords grid, rhs, function
1677: .seealso GridEvaluateRhsOperator(), GridEvaluateRhs(), GridEvaluateSystemMatrix()
1678: @*/
1679: int GridEvaluateRhsFunction(Grid grid, GVec x, GVec f, void *ctx)
1680: {
1681: PetscTruth oldLinear = grid->activeOpTypes[1];
1682: PetscTruth oldNonlinear = grid->activeOpTypes[2];
1683: int ierr;
1687: /* x can be PETSC_NULL */
1689: /* Turn off computation of operators */
1690: grid->activeOpTypes[1] = PETSC_FALSE;
1691: grid->activeOpTypes[2] = PETSC_FALSE;
1692: GridEvaluateRhs(grid, x, f, (PetscObject) ctx);
1693: grid->activeOpTypes[1] = oldLinear;
1694: grid->activeOpTypes[2] = oldNonlinear;
1695: return(0);
1696: }
1698: #undef __FUNCT__
1700: /*@C GridEvaluateRhsOperator
1701: This function constructs the weak form of the operators specified with GridSetRhsOperator().
1703: Collective on Grid
1705: Input Parameters:
1706: + grid - The grid
1707: . x - The current iterate
1708: . linear - The flag for including linear operators
1709: . nonlinear - The flag for including nonlinear operators
1710: - ctx - The optional user context
1712: Output Parameter:
1713: . f - The operator value
1715: Level: beginner
1717: .keywords grid, rhs, operator
1718: .seealso GridEvaluateRhsFunction(), GridEvaluateRhs(), GridEvaluateSystemMatrix()
1719: @*/
1720: int GridEvaluateRhsOperator(Grid grid, GVec x, GVec f, PetscTruth linear, PetscTruth nonlinear, void *ctx)
1721: {
1722: PetscTruth oldFunction = grid->activeOpTypes[0];
1723: PetscTruth oldLinear = grid->activeOpTypes[1];
1724: PetscTruth oldNonlinear = grid->activeOpTypes[2];
1725: int ierr;
1729: /* x can be PETSC_NULL */
1731: /* Turn off computation of operators */
1732: grid->activeOpTypes[0] = PETSC_FALSE;
1733: grid->activeOpTypes[1] = linear;
1734: grid->activeOpTypes[2] = nonlinear;
1735: GridEvaluateRhs(grid, x, f, (PetscObject) ctx);
1736: grid->activeOpTypes[0] = oldFunction;
1737: grid->activeOpTypes[1] = oldLinear;
1738: grid->activeOpTypes[2] = oldNonlinear;
1739: return(0);
1740: }
1742: #undef __FUNCT__
1744: /*@C GridEvaluateSystemMatrix
1745: This function constructs the weak form of the linear operators
1746: specified with GridSetMatOperator().
1748: Collective on Grid
1750: Input Parameters:
1751: + grid - The grid
1752: . x - The current iterate
1753: - ctx - The optional user context
1755: Output Parameter:
1756: . f - The function value
1758: Level: beginner
1760: .keywords: grid, matrix
1761: .seealso: GridEvaluateRhs()
1762: @*/
1763: int GridEvaluateSystemMatrix(Grid grid, GVec x, GMat *J, GMat *M, MatStructure *flag, PetscObject ctx)
1764: {
1765: Mesh mesh;
1766: MeshMover mover;
1767: Grid meshVelGrid;
1768: #ifdef PETSC_HAVE_MATHEMATICA
1769: PetscTruth opt;
1770: #endif
1771: int ierr;
1775: if (grid->ALEActive == PETSC_TRUE) {
1776: GridGetMesh(grid, &mesh);
1777: MeshGetMover(mesh, &mover);
1778: MeshMoverGetVelocityGrid(mover, &meshVelGrid);
1779: GridSetUp(meshVelGrid);
1780: }
1781: if (grid->isMatrixFree == PETSC_TRUE) {
1782: /* Should probably set this as the matrix context, but there is no MatShellSetContext() */
1783: grid->matrixFreeArg = x;
1784: grid->matrixFreeContext = ctx;
1785: } else {
1786: (*grid->ops->gridevaluatesystemmatrix)(grid, x, J, M, flag, ctx);
1787: }
1788: #ifdef PETSC_HAVE_MATHEMATICA
1789: /* Debugging */
1790: PetscOptionsHasName(PETSC_NULL, "-trace_grid_math", &opt);
1791: if (opt == PETSC_TRUE) {
1792: ViewerMathematicaSetName(PETSC_VIEWER_MATHEMATICA_WORLD, "jac");
1793: MatView(*J, PETSC_VIEWER_MATHEMATICA_WORLD);
1794: ViewerMathematicaClearName(PETSC_VIEWER_MATHEMATICA_WORLD);
1795: }
1796: #endif
1797: return(0);
1798: }
1800: #if 0
1801: /*----------------------------------------- Parallel Communication Functions ----------------------------------------*/
1802: #undef __FUNCT__
1804: /*@C
1805: GridGhostExchange - This functions transfers data between local and ghost storage without a predefined mapping.
1807: Collective on MPI_Comm
1809: Input Parameters:
1810: + comm - The communicator
1811: . numGhosts - The number of ghosts in this domain
1812: . ghostProcs - The processor from which to obtain each ghost
1813: . ghostIndices - The global index for each ghost
1814: . dataType - The type of the variables
1815: . firstVar - The first variable on each processor
1816: . addv - The insert mode, INSERT_VALUES or ADD_VALUES
1817: - mode - The direction of the transfer, SCATTER_FORWARD or SCATTER_REVERSE
1819: Output Paramters:
1820: + locVars - The local variable array
1821: - ghostVars - The ghost variables
1823: Note:
1824: The data in ghostVars is assumed contiguous and implicitly indexed by the order of
1825: ghostProcs and ghostIndices. The SCATTER_FORWARD mode will take the requested data
1826: from locVars and copy it to ghostVars in the order specified by ghostIndices. The
1827: SCATTER_REVERSE mode will take data from ghostVars and copy it to locVars.
1829: Level: advanced
1831: .keywords: ghost, exchange, grid
1832: .seealso: GridGlobalToLocal(), GridLocalToGlobal()
1833: @*/
1834: int GridGhostExchange(MPI_Comm comm, int numGhosts, int *ghostProcs, int *ghostIndices, PetscDataType dataType,
1835: int *firstVar, InsertMode addv, ScatterMode mode, void *locVars, void *ghostVars)
1836: {
1837: int *numSendGhosts; /* The number of ghosts from each domain */
1838: int *numRecvGhosts; /* The number of local variables which are ghosts in each domain */
1839: int *sumSendGhosts; /* The prefix sums of numSendGhosts */
1840: int *sumRecvGhosts; /* The prefix sums of numRecvGhosts */
1841: int *offsets; /* The offset into the send array for each domain */
1842: int totSendGhosts; /* The number of ghosts to request variables for */
1843: int totRecvGhosts; /* The number of nodes to provide class info about */
1844: int *sendIndices; /* The canonical indices of ghosts in this domain */
1845: int *recvIndices; /* The canonical indices of ghosts to return variables for */
1846: char *tempVars; /* The variables of the requested or submitted ghosts */
1847: char *locBytes = (char *) locVars;
1848: MPI_Datatype MPIType;
1849: int typeSize;
1850: #ifdef PETSC_USE_BOPT_g
1851: int numLocVars;
1852: #endif
1853: int numProcs, rank;
1854: int proc, ghost, locIndex, byte;
1855: int ierr;
1858: /* Initialize communication */
1859: MPI_Comm_size(comm, &numProcs);
1860: MPI_Comm_rank(comm, &rank);
1861: PetscMalloc(numProcs * sizeof(int), &numSendGhosts);
1862: PetscMalloc(numProcs * sizeof(int), &numRecvGhosts);
1863: PetscMalloc(numProcs * sizeof(int), &sumSendGhosts);
1864: PetscMalloc(numProcs * sizeof(int), &sumRecvGhosts);
1865: PetscMalloc(numProcs * sizeof(int), &offsets);
1866: PetscMemzero(numSendGhosts, numProcs * sizeof(int));
1867: PetscMemzero(numRecvGhosts, numProcs * sizeof(int));
1868: PetscMemzero(sumSendGhosts, numProcs * sizeof(int));
1869: PetscMemzero(sumRecvGhosts, numProcs * sizeof(int));
1870: PetscMemzero(offsets, numProcs * sizeof(int));
1871: #ifdef PETSC_USE_BOPT_g
1872: numLocVars = firstVar[rank+1] - firstVar[rank];
1873: #endif
1875: /* Get number of ghosts needed from each processor */
1876: for(ghost = 0; ghost < numGhosts; ghost++)
1877: numSendGhosts[ghostProcs[ghost]]++;
1879: /* Get number of ghosts to provide variables for */
1880: MPI_Alltoall(numSendGhosts, 1, MPI_INT, numRecvGhosts, 1, MPI_INT, comm);
1881: for(proc = 1; proc < numProcs; proc++) {
1882: sumSendGhosts[proc] = sumSendGhosts[proc-1] + numSendGhosts[proc-1];
1883: sumRecvGhosts[proc] = sumRecvGhosts[proc-1] + numRecvGhosts[proc-1];
1884: offsets[proc] = sumSendGhosts[proc];
1885: }
1886: totSendGhosts = sumSendGhosts[numProcs-1] + numSendGhosts[numProcs-1];
1887: totRecvGhosts = sumRecvGhosts[numProcs-1] + numRecvGhosts[numProcs-1];
1888: if (numGhosts != totSendGhosts) {
1889: SETERRQ2(PETSC_ERR_PLIB, "Invalid number of ghosts %d in send, should be %d", totSendGhosts, numGhosts);
1890: }
1892: PetscDataTypeGetSize(dataType, &typeSize);
1893: if (totSendGhosts) {
1894: sendIndices = (int *) PetscMalloc(totSendGhosts * sizeof(int)); CHKERRQ(sendIndices);
1895: }
1896: if (totRecvGhosts) {
1897: recvIndices = (int *) PetscMalloc(totRecvGhosts * sizeof(int)); CHKERRQ(recvIndices);
1898: tempVars = (char *) PetscMalloc(totRecvGhosts * typeSize); CHKERRQ(tempVars);
1899: }
1901: /* Must order ghosts by processor */
1902: for(ghost = 0; ghost < numGhosts; ghost++)
1903: sendIndices[offsets[ghostProcs[ghost]]++] = ghostIndices[ghost];
1905: /* Get canonical indices of ghosts to provide variables for */
1906: MPI_Alltoallv(sendIndices, numSendGhosts, sumSendGhosts, MPI_INT,
1907: recvIndices, numRecvGhosts, sumRecvGhosts, MPI_INT, comm);
1908:
1910: switch(mode)
1911: {
1912: case SCATTER_FORWARD:
1913: /* Get ghost variables */
1914: if (addv == INSERT_VALUES) {
1915: for(ghost = 0; ghost < totRecvGhosts; ghost++)
1916: {
1917: locIndex = recvIndices[ghost] - firstVar[rank];
1918: #ifdef PETSC_USE_BOPT_g
1919: if ((locIndex < 0) || (locIndex >= numLocVars)) {
1920: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1921: }
1922: #endif
1923: for(byte = 0; byte < typeSize; byte++)
1924: tempVars[ghost*typeSize+byte] = locBytes[locIndex*typeSize+byte];
1925: }
1926: } else {
1927: for(ghost = 0; ghost < totRecvGhosts; ghost++)
1928: {
1929: locIndex = recvIndices[ghost] - firstVar[rank];
1930: #ifdef PETSC_USE_BOPT_g
1931: if ((locIndex < 0) || (locIndex >= numLocVars)) {
1932: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1933: }
1934: #endif
1935: for(byte = 0; byte < typeSize; byte++)
1936: tempVars[ghost*typeSize+byte] += locBytes[locIndex*typeSize+byte];
1937: }
1938: }
1940: /* Communicate local variables to ghost storage */
1941: PetscDataTypeToMPIDataType(dataType, &MPIType);
1942: MPI_Alltoallv(tempVars, numRecvGhosts, sumRecvGhosts, MPIType,
1943: ghostVars, numSendGhosts, sumSendGhosts, MPIType, comm);
1944:
1945: break;
1946: case SCATTER_REVERSE:
1947: /* Communicate ghost variables to local storage */
1948: PetscDataTypeToMPIDataType(dataType, &MPIType);
1949: MPI_Alltoallv(ghostVars, numSendGhosts, sumSendGhosts, MPIType,
1950: tempVars, numRecvGhosts, sumRecvGhosts, MPIType, comm);
1951:
1953: /* Get ghost variables */
1954: if (addv == INSERT_VALUES) {
1955: for(ghost = 0; ghost < totRecvGhosts; ghost++)
1956: {
1957: locIndex = recvIndices[ghost] - firstVar[rank];
1958: #ifdef PETSC_USE_BOPT_g
1959: if ((locIndex < 0) || (locIndex >= numLocVars)) {
1960: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1961: }
1962: #endif
1963: for(byte = 0; byte < typeSize; byte++)
1964: locBytes[locIndex*typeSize+byte] = tempVars[ghost*typeSize+byte];
1965: }
1966: } else {
1967: /* There must be a better way to do this -- Ask Bill */
1968: if (typeSize == sizeof(int)) {
1969: int *tempInt = (int *) tempVars;
1970: int *locInt = (int *) locVars;
1972: for(ghost = 0; ghost < totRecvGhosts; ghost++) {
1973: locIndex = recvIndices[ghost] - firstVar[rank];
1974: #ifdef PETSC_USE_BOPT_g
1975: if ((locIndex < 0) || (locIndex >= numLocVars)) {
1976: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1977: }
1978: #endif
1979: locInt[locIndex] += tempInt[ghost];
1980: }
1981: } else if (typeSize == sizeof(long int)) {
1982: long int *tempLongInt = (long int *) tempVars;
1983: long int *locLongInt = (long int *) locVars;
1985: for(ghost = 0; ghost < totRecvGhosts; ghost++) {
1986: locIndex = recvIndices[ghost] - firstVar[rank];
1987: #ifdef PETSC_USE_BOPT_g
1988: if ((locIndex < 0) || (locIndex >= numLocVars)) {
1989: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1990: }
1991: #endif
1992: locLongInt[locIndex] += tempLongInt[ghost];
1993: }
1994: } else {
1995: for(ghost = 0; ghost < totRecvGhosts; ghost++) {
1996: locIndex = recvIndices[ghost] - firstVar[rank];
1997: #ifdef PETSC_USE_BOPT_g
1998: if ((locIndex < 0) || (locIndex >= numLocVars)) {
1999: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
2000: }
2001: #endif
2002: for(byte = 0; byte < typeSize; byte++)
2003: locBytes[locIndex*typeSize+byte] += tempVars[ghost*typeSize+byte];
2004: }
2005: }
2006: }
2007: break;
2008: default:
2009: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid scatter mode %d", mode);
2010: }
2012: /* Cleanup */
2013: PetscFree(numSendGhosts);
2014: PetscFree(numRecvGhosts);
2015: PetscFree(sumSendGhosts);
2016: PetscFree(sumRecvGhosts);
2017: PetscFree(offsets);
2018: if (totSendGhosts) {
2019: PetscFree(sendIndices);
2020: }
2021: if (totRecvGhosts) {
2022: PetscFree(recvIndices);
2023: PetscFree(tempVars);
2024: }
2025: return(0);
2026: }
2027: #endif