Actual source code: fieldClassMap2d.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: fieldClassMap2d.c,v 1.6 2000/01/30 18:30:46 huangp Exp $";
3: #endif
5: #include "src/grid/gridimpl.h" /*I "grid.h" I*/
6: #include "fieldClassMap2d.h"
8: /*---------------------------------------- Creation and Destruction Functions ---------------------------------------*/
9: #undef __FUNCT__
11: /*@C
12: FieldClassMapCreateTriangular2D - This function creates a class map for the given grid.
14: Collective on Grid
16: Input Parameters:
17: + grid - The underlying grid for the ordering
18: . numFields - The number of fields in the ordering
19: - fields - The fields in the ordering
21: Output Parameter:
22: . map - The map
24: Level: beginner
26: .keywords: class, field class, class map, create
27: .seealso: FieldClassMapCreate(), FieldClassMapDestroy()
28: @*/
29: int FieldClassMapCreateTriangular2D(Grid grid, int numFields, int *fields, FieldClassMap *map)
30: {
31: MPI_Comm comm;
32: ParameterDict dict;
33: int numTotalFields, f;
34: int ierr;
37: GridGetNumFields(grid, &numTotalFields);
38: if (numFields > numTotalFields) {
39: SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid number %d of fields, grid only has %d", numFields, numTotalFields);
40: }
41: for(f = 0; f < numFields; f++) {
42: if ((fields[f] < 0) || (fields[f] >= numTotalFields)) {
43: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", fields[f]);
44: }
45: }
46: PetscObjectGetComm((PetscObject) grid, &comm);
47: FieldClassMapCreate(comm, map);
48: PetscObjectGetComm((PetscObject) grid, &comm);
49: ParameterDictCreate(comm, &dict);
50: ParameterDictSetInteger(dict, "numFields", numFields);
51: ParameterDictSetObject(dict, "fields", fields);
52: ParameterDictSetObject(dict, "grid", grid);
53: PetscObjectSetParameterDict((PetscObject) *map, dict);
54: ParameterDictDestroy(dict);
55: FieldClassMapSetType(*map, CLASS_MAP_TRIANGULAR_2D);
56: return(0);
57: }
59: #undef __FUNCT__
61: int FieldClassMapDestroy_Triangular_2D(FieldClassMap map)
62: {
64: return(0);
65: }
67: #undef __FUNCT__
69: int FieldClassMapConstrain_Triangular_2D(FieldClassMap map, Grid grid, PetscTruth useBC, PetscTruth useConst,
70: FieldClassMap *newMap)
71: {
72: FieldClassMap cm;
73: Mesh mesh;
74: Partition part;
75: int numBd = grid->numBd;
76: int numBC = grid->numBC;
77: GridBC *gBC = grid->bc;
78: int numPointBC = grid->numPointBC;
79: GridBC *pBC = grid->pointBC;
80: int numFields = map->numFields;
81: int *fields = map->fields;
82: int numNodes = map->numNodes;
83: int numOverlapNodes = map->numOverlapNodes;
84: int numClasses = map->numClasses;
85: int *classes = map->classes;
86: int ***localClassMap, ***classMap;
87: PetscTruth isConst;
88: int f, g, field, comp, node, bd, oldBd, bc, marker, nclass, newNClass, oldNClass;
89: int ierr;
92: FieldClassMapCreate(map->comm, &cm);
93: PetscMemcpy(cm->ops, map->ops, sizeof(FieldClassMapOps));
94: PetscStrallocpy(map->serialize_name, &cm->serialize_name);
96: GridGetMesh(grid, &mesh);
97: cm->numFields = numFields;
98: cm->numNodes = numNodes;
99: cm->numOverlapNodes = numOverlapNodes;
100: cm->numGhostNodes = cm->numOverlapNodes - cm->numNodes;
101: PetscMalloc(numFields * sizeof(int), &cm->fields);
102: PetscMemcpy(cm->fields, fields, numFields * sizeof(int));
103: /*
104: classMap[numBd+numPointBC+1][numClasses][numFields]:
106: The classMap[][][] structure is organized as follows. For every boundary, we determine the classes which are
107: present on that boundary. For each class, we create a new class which does not contain the constrained fields.
108: The classMap[][][] structure has for each boundary
110: flag_0, flag_1, ldots, flag_numFields-1
111: flag_0, flag_1, ldots, flag_numFields-1
112: ...
113: flag_0, flag_1, ldots, flag_numFields-1
115: where a positive flag indicates that the field is constrained in that class, and a zero indicates no change.
116: Point boundary conditions are treated like additional boundaries, and since we only allow one constraint,
117: the last position is used for the constraint.
118: */
119: cm->mapSize = numBd + numPointBC + 1;
120: cm->numOldClasses = numClasses;
121: PetscMalloc(cm->mapSize * sizeof(int **), &localClassMap);
122: PetscMalloc(cm->mapSize * sizeof(int **), &classMap);
123: for(bd = 0; bd < cm->mapSize; bd++) {
124: PetscMalloc(numClasses * sizeof(int *), &localClassMap[bd]);
125: PetscMalloc(numClasses * sizeof(int *), &classMap[bd]);
126: for(nclass = 0; nclass < numClasses; nclass++) {
127: PetscMalloc(numFields * sizeof(int), &localClassMap[bd][nclass]);
128: PetscMalloc(numFields * sizeof(int), &classMap[bd][nclass]);
129: PetscMemzero(localClassMap[bd][nclass], numFields * sizeof(int));
130: PetscMemzero(classMap[bd][nclass], numFields * sizeof(int));
131: }
132: }
133: /* Mark constrained classes for reassignment */
134: for(node = 0; node < numNodes; node++) {
135: MeshGetNodeBoundary(mesh, node, &marker);
136: if (marker == 0)
137: continue;
138: nclass = classes[node];
139: if ((useConst == PETSC_TRUE) && (marker < 0)) {
140: /* All constrained nodes have negative markers */
141: for(f = 0; f < numFields; f++) {
142: field = fields[f];
143: if (grid->fields[field].isConstrained == PETSC_TRUE) {
144: localClassMap[cm->mapSize-1][nclass][f] = 1;
145: }
146: }
147: } else if ((useBC == PETSC_TRUE) && (marker > 0)) {
148: for(bc = 0; bc < numBC; bc++) {
149: if (gBC[bc].reduce == PETSC_FALSE) continue;
150: if (gBC[bc].boundary == marker) {
151: MeshGetBoundaryIndex(mesh, marker, &bd);
152: for(f = 0; f < numFields; f++) {
153: if (fields[f] == gBC[bc].field) {
154: localClassMap[bd][nclass][f] = 1;
155: break;
156: }
157: }
158: }
159: }
160: for(bc = 0; bc < numPointBC; bc++) {
161: if (pBC[bc].reduce == PETSC_FALSE) continue;
162: if (pBC[bc].node == node) {
163: for(f = 0; f < numFields; f++) {
164: if (fields[f] == pBC[bc].field) {
165: localClassMap[numBd+bc][nclass][f] = 1;
166: break;
167: }
168: }
169: }
170: }
171: }
172: }
173: /* Synchronize localClassMap */
174: for(bd = 0; bd < cm->mapSize; bd++) {
175: for(nclass = 0; nclass < numClasses; nclass++) {
176: MPI_Allreduce(localClassMap[bd][nclass], classMap[bd][nclass], numFields, MPI_INT, MPI_LOR, cm->comm);
177:
178: }
179: }
180: /* Assign new classes */
181: PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);
182: for(bd = 0; bd < cm->mapSize; bd++) {
183: PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);
184: for(nclass = 0; nclass < numClasses; nclass++) {
185: cm->classMap[bd][nclass] = nclass;
186: }
187: }
188: for(bd = 0, cm->numClasses = numClasses; bd < cm->mapSize; bd++) {
189: for(nclass = 0; nclass < numClasses; nclass++) {
190: for(f = 0, isConst = PETSC_FALSE; f < numFields; f++) {
191: if (classMap[bd][nclass][f] != 0) {
192: isConst = PETSC_TRUE;
193: field = fields[f];
194: /* Map old classes to new with the condition which created it for a given field */
195: classMap[bd][nclass][f] = cm->numClasses;
196: /* Check that field should be constrained */
197: if ((bd == cm->mapSize-1) && (grid->fields[field].isConstrained == PETSC_FALSE)) {
198: SETERRQ1(PETSC_ERR_PLIB, "Invalid constrained field %d", field);
199: }
200: /* Must include previous constraints on the node */
201: for(g = 0; g < numFields; g++) {
202: for(oldBd = 0; oldBd < bd; oldBd++) {
203: for(oldNClass = 0; oldNClass < numClasses; oldNClass++) {
204: if (classMap[oldBd][oldNClass][g] != 0) {
205: classMap[bd][nclass][g] = cm->numClasses;
206: break;
207: }
208: }
209: }
210: }
211: break;
212: }
213: }
214: /* Create a new class */
215: if (isConst == PETSC_TRUE) {
216: cm->classMap[bd][nclass] = cm->numClasses++;
217: }
218: }
219: }
220: /* Set new classes on nodes */
221: PetscMalloc(numOverlapNodes * sizeof(int), &cm->classes);
222: PetscMemcpy(cm->classes, classes, numOverlapNodes * sizeof(int));
223: for(node = 0; node < numOverlapNodes; node++) {
224: MeshGetNodeBoundary(mesh, node, &marker);
225: if (marker == 0)
226: continue;
227: nclass = classes[node];
228: if (marker < 0) {
229: cm->classes[node] = cm->classMap[cm->mapSize-1][nclass];
230: } else {
231: for(bc = 0; bc < numBC; bc++) {
232: if (gBC[bc].reduce == PETSC_FALSE) continue;
233: if (gBC[bc].boundary == marker) {
234: MeshGetBoundaryIndex(mesh, marker, &bd);
235: if (cm->classMap[bd][nclass] != nclass) {
236: cm->classes[node] = cm->classMap[bd][nclass];
237: }
238: }
239: }
240: /* Notice that point BC must come last since these nodes may lie on a boundary which was also constrained */
241: for(bc = 0; bc < numPointBC; bc++) {
242: if (pBC[bc].reduce == PETSC_FALSE) continue;
243: if (pBC[bc].node == node) {
244: if (cm->classMap[numBd+bc][nclass] != nclass) {
245: cm->classes[node] = cm->classMap[numBd+bc][nclass];
246: }
247: }
248: }
249: }
250: }
251: /* Calculate new class structure */
252: cm->isReduced = map->isReduced;
253: cm->isConstrained = map->isConstrained;
254: PetscMalloc(numFields * sizeof(int *), &cm->fieldClasses);
255: PetscMalloc(numFields * sizeof(int *), &cm->reduceFieldClasses);
256: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizes);
257: PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);
258: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);
259: PetscMemzero(cm->isClassConstrained, cm->numClasses * sizeof(int));
260: PetscMemzero(cm->classSizeDiffs, cm->numClasses * sizeof(int));
261: /* Copy old structure */
262: for(nclass = 0; nclass < numClasses; nclass++) {
263: cm->classSizes[nclass] = map->classSizes[nclass];
264: cm->isClassConstrained[nclass] = map->isClassConstrained[nclass];
265: cm->classSizeDiffs[nclass] = map->classSizeDiffs[nclass];
266: }
267: for(newNClass = numClasses; newNClass < cm->numClasses; newNClass++) {
268: /* Invert the class map */
269: for(bd = 0; bd < cm->mapSize; bd++) {
270: for(nclass = 0; nclass < numClasses; nclass++)
271: if (cm->classMap[bd][nclass] == newNClass)
272: break;
273: if (nclass < numClasses)
274: break;
275: }
276: if ((nclass >= numClasses) || (bd >= cm->mapSize)) {
277: SETERRQ1(PETSC_ERR_PLIB, "Unable to locate class %d in class map", newNClass);
278: }
280: cm->classSizes[newNClass] = map->classSizes[nclass];
281: cm->classSizeDiffs[newNClass] = map->classSizeDiffs[nclass];
282: }
283: /* Apply constraints */
284: for(f = 0; f < numFields; f++) {
285: PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);
286: PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);
287: PetscMemzero(cm->reduceFieldClasses[f], cm->numClasses * sizeof(int));
288: /* Copy old structure */
289: for(nclass = 0; nclass < numClasses; nclass++) {
290: cm->fieldClasses[f][nclass] = map->fieldClasses[f][nclass];
291: if (map->reduceFieldClasses) {
292: cm->reduceFieldClasses[f][nclass] = map->reduceFieldClasses[f][nclass];
293: }
294: }
296: /* Remove constrained fields from new classes */
297: for(newNClass = numClasses; newNClass < cm->numClasses; newNClass++) {
298: /* Invert the class map */
299: for(bd = 0; bd < cm->mapSize; bd++) {
300: for(nclass = 0; nclass < numClasses; nclass++)
301: if (cm->classMap[bd][nclass] == newNClass)
302: break;
303: if (nclass < numClasses)
304: break;
305: }
306: if ((nclass >= numClasses) || (bd >= cm->mapSize)) {
307: SETERRQ1(PETSC_ERR_PLIB, "Unable to locate class %d in class map", newNClass);
308: }
309: /* Check every BC for constraints on this field/class */
310: cm->fieldClasses[f][newNClass] = map->fieldClasses[f][nclass];
311: for(bd = 0; bd < cm->mapSize; bd++) {
312: if (classMap[bd][nclass][f] == newNClass) {
313: /* New classes can be redefined by another BC to be a higher number, but do not
314: operate on the existing new class with that BC (hence the second check) */
315: if (cm->classMap[bd][nclass] != newNClass) {
316: SETERRQ3(PETSC_ERR_PLIB, "Invalid bc %d constraining class %d to %d", bd, nclass, newNClass);
317: }
318: /* If the field is constrained, it has no values on these nodes */
319: cm->fieldClasses[f][newNClass] = 0;
320: if (bd < cm->mapSize-1) {
321: cm->isReduced = PETSC_TRUE;
322: cm->reduceFieldClasses[f][newNClass] = 1;
323: }
324: /* Subtract the degrees of freedom of the constrained field for each BC */
325: GridGetFieldComponents(grid, fields[f], &comp);
326: cm->classSizes[newNClass] -= comp;
327: if (cm->classSizes[newNClass] < 0) {
328: SETERRQ2(PETSC_ERR_PLIB, "Invalid size for class %d (derived from class %d)", newNClass, nclass);
329: }
330: /* Handle user constraints separately from BC */
331: if (bd == cm->mapSize-1) {
332: cm->isConstrained = PETSC_TRUE;
333: /* This class is constrained */
334: cm->isClassConstrained[newNClass] = 1;
335: /* There are this many components in the new field */
336: cm->classSizeDiffs[newNClass] += grid->fields[fields[f]].constraintCompDiff + comp;
337: }
338: }
339: }
340: }
341: }
343: /* Cleanup */
344: if (cm->isReduced == PETSC_FALSE) {
345: for(f = 0; f < map->numFields; f++) {
346: PetscFree(cm->reduceFieldClasses[f]);
347: }
348: PetscFree(cm->reduceFieldClasses);
349: }
350: for(bd = 0; bd < cm->mapSize; bd++) {
351: for(nclass = 0; nclass < numClasses; nclass++) {
352: PetscFree(localClassMap[bd][nclass]);
353: PetscFree(classMap[bd][nclass]);
354: }
355: PetscFree(localClassMap[bd]);
356: PetscFree(classMap[bd]);
357: }
358: PetscFree(localClassMap);
359: PetscFree(classMap);
361: /* Communicate ghost classes */
362: MeshGetPartition(mesh, &part);
363: PartitionGhostNodeExchange(part, INSERT_VALUES, SCATTER_FORWARD, cm->classes, &cm->classes[cm->numNodes]);
364:
366: PetscLogObjectMemory(cm, (cm->mapSize + numFields*2) * sizeof(int *) + (numFields + cm->numOldClasses*cm->mapSize +
367: cm->numClasses*(numFields*2 + 3) + numOverlapNodes) * sizeof(int));
368: *newMap = cm;
369: return(0);
370: }
372: #undef __FUNCT__
374: int FieldClassMapReduce_Triangular_2D(FieldClassMap map, Grid grid, FieldClassMap *newMap) {
375: FieldClassMap cm;
376: Mesh mesh;
377: Partition part;
378: int numBd = grid->numBd;
379: int numBC = grid->numBC;
380: GridBC *gBC = grid->bc;
381: int numPointBC = grid->numPointBC;
382: GridBC *pBC = grid->pointBC;
383: int numFields = map->numFields;
384: int *fields = map->fields;
385: int numNodes = map->numNodes;
386: int numOverlapNodes = map->numOverlapNodes;
387: int numClasses = map->numClasses;
388: int *classes = map->classes;
389: int ***localClassMap, ***classMap;
390: PetscTruth isConst;
391: int f, g, field, comp, node, bd, bc, marker, nclass, newNClass;
392: int ierr;
395: FieldClassMapCreate(map->comm, &cm);
396: PetscMemcpy(cm->ops, map->ops, sizeof(FieldClassMapOps));
397: PetscStrallocpy(map->serialize_name, &cm->serialize_name);
399: GridGetMesh(grid, &mesh);
400: cm->numNodes = numNodes;
401: cm->numOverlapNodes = numOverlapNodes;
402: cm->numGhostNodes = cm->numOverlapNodes - cm->numNodes;
403: PetscMalloc(numFields * sizeof(int), &cm->fields);
404: /*
405: classMap[numBd+numPointBC][numClasses][numFields]:
407: The classMap[][][] structure is organized as follows. For every boundary, we determine the classes which are
408: present on that boundary. For each class, we create a new class which does not contain the constrained fields.
409: The classMap[][][] structure has for each boundary
411: flag_0, flag_1, ldots, flag_numFields-1
412: flag_0, flag_1, ldots, flag_numFields-1
413: ...
414: flag_0, flag_1, ldots, flag_numFields-1
416: where a positive flag indicates that the field is constrained in that class, and a zero indicates no change.
417: Point boundary conditions are treated like additional boundaries.
418: */
419: cm->mapSize = numBd + numPointBC;
420: cm->numOldClasses = numClasses;
421: PetscMalloc(cm->mapSize * sizeof(int **), &localClassMap);
422: PetscMalloc(cm->mapSize * sizeof(int **), &classMap);
423: for(bd = 0; bd < cm->mapSize; bd++) {
424: PetscMalloc(numClasses * sizeof(int *), &localClassMap[bd]);
425: PetscMalloc(numClasses * sizeof(int *), &classMap[bd]);
426: for(nclass = 0; nclass < numClasses; nclass++) {
427: PetscMalloc(numFields * sizeof(int), &localClassMap[bd][nclass]);
428: PetscMalloc(numFields * sizeof(int), &classMap[bd][nclass]);
429: PetscMemzero(localClassMap[bd][nclass], numFields * sizeof(int));
430: PetscMemzero(classMap[bd][nclass], numFields * sizeof(int));
431: }
432: }
433: /* Mark constrained classes for reassignment */
434: for(node = 0; node < numNodes; node++) {
435: MeshGetNodeBoundary(mesh, node, &marker);
436: if (marker == 0)
437: continue;
438: nclass = classes[node];
439: if (marker > 0) {
440: for(bc = 0; bc < numBC; bc++) {
441: if (gBC[bc].reduce == PETSC_FALSE)
442: continue;
443: if (gBC[bc].boundary == marker) {
444: MeshGetBoundaryIndex(mesh, marker, &bd);
445: for(f = 0; f < numFields; f++) {
446: if (fields[f] == gBC[bc].field) {
447: localClassMap[bd][nclass][f] = 1;
448: break;
449: }
450: }
451: }
452: }
453: for(bc = 0; bc < numPointBC; bc++) {
454: if (pBC[bc].reduce == PETSC_FALSE) continue;
455: if (pBC[bc].node == node) {
456: for(f = 0; f < numFields; f++) {
457: if (fields[f] == pBC[bc].field) {
458: localClassMap[numBd+bc][nclass][f] = 1;
459: break;
460: }
461: }
462: }
463: }
464: }
465: }
466: /* Synchronize localClassMap */
467: for(bd = 0; bd < cm->mapSize; bd++) {
468: for(nclass = 0; nclass < numClasses; nclass++) {
469: MPI_Allreduce(localClassMap[bd][nclass], classMap[bd][nclass], numFields, MPI_INT, MPI_LOR, cm->comm);
470:
471: }
472: }
474: /* Assign new classes and build field list */
475: PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);
476: for(bd = 0; bd < cm->mapSize; bd++) {
477: PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);
478: for(nclass = 0; nclass < numClasses; nclass++) {
479: cm->classMap[bd][nclass] = nclass;
480: }
481: }
482: cm->numFields = 0;
483: for(bd = 0, cm->numClasses = numClasses; bd < cm->mapSize; bd++) {
484: for(nclass = 0; nclass < numClasses; nclass++) {
485: for(f = 0, isConst = PETSC_FALSE; f < numFields; f++) {
486: if (classMap[bd][nclass][f] != 0) {
487: isConst = PETSC_TRUE;
488: field = fields[f];
489: /* Add to list of reduced fields */
490: for(g = 0; g < cm->numFields; g++)
491: if (field == cm->fields[g])
492: break;
493: if (g == cm->numFields)
494: cm->fields[cm->numFields++] = field;
495: }
496: }
497: /* Create a new class */
498: if (isConst == PETSC_TRUE) {
499: cm->classMap[bd][nclass] = cm->numClasses++;
500: }
501: }
502: }
504: /* Set new classes on nodes */
505: PetscMalloc(numOverlapNodes * sizeof(int), &cm->classes);
506: PetscMemcpy(cm->classes, classes, numOverlapNodes * sizeof(int));
507: for(node = 0; node < numOverlapNodes; node++) {
508: MeshGetNodeBoundary(mesh, node, &marker);
509: if (marker == 0)
510: continue;
511: nclass = classes[node];
512: for(bc = 0; bc < numBC; bc++) {
513: if (gBC[bc].reduce == PETSC_FALSE) continue;
514: if (gBC[bc].boundary == marker) {
515: MeshGetBoundaryIndex(mesh, marker, &bd);
516: cm->classes[node] = cm->classMap[bd][nclass];
517: }
518: }
519: /* Notice that point BC must come last since these nodes may lie on a boundary which was also constrained */
520: for(bc = 0; bc < numPointBC; bc++) {
521: if (pBC[bc].reduce == PETSC_FALSE) continue;
522: if (pBC[bc].node == node) {
523: cm->classes[node] = cm->classMap[numBd+bc][nclass];
524: }
525: }
526: }
527: /* Calculate new class structure */
528: cm->isReduced = PETSC_TRUE;
529: cm->isConstrained = map->isConstrained;
530: PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);
531: PetscMalloc(cm->numFields * sizeof(int *), &cm->reduceFieldClasses);
532: PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);
533: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);
534: PetscMemzero(cm->isClassConstrained, cm->numClasses * sizeof(int));
535: PetscMemzero(cm->classSizeDiffs, cm->numClasses * sizeof(int));
536: /* Copy old structure */
537: for(nclass = 0; nclass < numClasses; nclass++) {
538: cm->isClassConstrained[nclass] = map->isClassConstrained[nclass];
539: cm->classSizeDiffs[nclass] = map->classSizeDiffs[nclass];
540: }
541: /* Apply constraints */
542: for(f = 0; f < cm->numFields; f++) {
543: PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);
544: PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);
545: PetscMemzero(cm->fieldClasses[f], cm->numClasses * sizeof(int));
546: /* Copy old structure */
547: for(nclass = 0; nclass < numClasses; nclass++) {
548: cm->reduceFieldClasses[f][nclass] = map->fieldClasses[f][nclass];
549: if (map->reduceFieldClasses) {
550: cm->fieldClasses[f][nclass] = map->reduceFieldClasses[f][nclass];
551: }
552: }
554: /* Remove constrained fields from new classes */
555: for(newNClass = numClasses; newNClass < cm->numClasses; newNClass++) {
556: /* Invert the class map */
557: for(bd = 0; bd < cm->mapSize; bd++) {
558: for(nclass = 0; nclass < numClasses; nclass++)
559: if (cm->classMap[bd][nclass] == newNClass)
560: break;
561: if (nclass < numClasses)
562: break;
563: }
564: if ((nclass >= numClasses) || (bd >= cm->mapSize)) {
565: SETERRQ1(PETSC_ERR_PLIB, "Unable to locate class %d in class map", newNClass);
566: }
567: /* Check every BC for constraints on this field/class */
568: cm->reduceFieldClasses[f][newNClass] = map->fieldClasses[f][nclass];
569: for(bd = 0; bd < cm->mapSize; bd++) {
570: /* New classes can be redefined by another BC to be a higher number, but do not
571: operate on the existing new class with that BC (hence the second check)
572: */
573: if ((classMap[bd][nclass][f] == 1) && (cm->classMap[bd][nclass] <= newNClass)) {
574: /* If the field is constrained, it has no values on these nodes */
575: cm->fieldClasses[f][newNClass] = 1;
576: cm->reduceFieldClasses[f][newNClass] = 0;
577: }
578: }
579: }
580: }
582: /* Setup class sizes */
583: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizes);
584: PetscMemzero(cm->classSizes, cm->numClasses * sizeof(int));
585: for(nclass = 0; nclass < cm->numClasses; nclass++) {
586: for(f = 0; f < cm->numFields; f++) {
587: field = cm->fields[f];
588: if (cm->fieldClasses[f][nclass]) {
589: GridGetFieldComponents(grid, field, &comp);
590: cm->classSizes[nclass] += comp;
591: }
592: }
593: }
595: /* Cleanup */
596: for(bd = 0; bd < cm->mapSize; bd++) {
597: for(nclass = 0; nclass < numClasses; nclass++) {
598: PetscFree(localClassMap[bd][nclass]);
599: PetscFree(classMap[bd][nclass]);
600: }
601: PetscFree(localClassMap[bd]);
602: PetscFree(classMap[bd]);
603: }
604: PetscFree(localClassMap);
605: PetscFree(classMap);
607: /* Communicate ghost classes */
608: MeshGetPartition(mesh, &part);
609: PartitionGhostNodeExchange(part, INSERT_VALUES, SCATTER_FORWARD, cm->classes, &cm->classes[cm->numNodes]);
610:
612: PetscLogObjectMemory(cm, (cm->mapSize + cm->numFields*2) * sizeof(int *) + (numFields + cm->numOldClasses*cm->mapSize +
613: cm->numClasses*(cm->numFields*2 + 3) + numOverlapNodes) * sizeof(int));
614: *newMap = cm;
615: return(0);
616: }
618: static FieldClassMapOps cmOps = {0,
619: FieldClassMapConstrain_Triangular_2D,
620: FieldClassMapReduce_Triangular_2D,
621: FieldClassMapDestroy_Triangular_2D,
622: 0};
624: EXTERN_C_BEGIN
625: #undef __FUNCT__
627: int FieldClassMapSerialize_Triangular_2D(MPI_Comm comm, FieldClassMap *map, PetscViewer viewer, PetscTruth store)
628: {
629: FieldClassMap cm;
630: int fd;
631: int zero = 0;
632: int one = 0;
633: int hasClassMap;
634: int bd, f;
635: int ierr;
638: PetscViewerBinaryGetDescriptor(viewer, &fd);
639: if (store) {
640: cm = *map;
641: PetscBinaryWrite(fd, &cm->numFields, 1, PETSC_INT, 0);
642: PetscBinaryWrite(fd, cm->fields, cm->numFields, PETSC_INT, 0);
643: PetscBinaryWrite(fd, &cm->numNodes, 1, PETSC_INT, 0);
644: PetscBinaryWrite(fd, &cm->numGhostNodes, 1, PETSC_INT, 0);
645: PetscBinaryWrite(fd, &cm->numOverlapNodes, 1, PETSC_INT, 0);
646: PetscBinaryWrite(fd, &cm->numClasses, 1, PETSC_INT, 0);
647: for(f = 0; f < cm->numFields; f++) {
648: PetscBinaryWrite(fd, cm->fieldClasses[f], cm->numClasses, PETSC_INT, 0);
649: }
650: PetscBinaryWrite(fd, cm->classes, cm->numOverlapNodes, PETSC_INT, 0);
651: PetscBinaryWrite(fd, cm->classSizes, cm->numClasses, PETSC_INT, 0);
652: PetscBinaryWrite(fd, &cm->isReduced, 1, PETSC_INT, 0);
653: if (cm->isReduced == PETSC_TRUE) {
654: for(f = 0; f < cm->numFields; f++) {
655: PetscBinaryWrite(fd, cm->reduceFieldClasses[f], cm->numClasses, PETSC_INT, 0);
656: }
657: }
658: PetscBinaryWrite(fd, &cm->isConstrained, 1, PETSC_INT, 0);
659: PetscBinaryWrite(fd, cm->isClassConstrained, cm->numClasses, PETSC_INT, 0);
660: PetscBinaryWrite(fd, cm->classSizeDiffs, cm->numClasses, PETSC_INT, 0);
661: PetscBinaryWrite(fd, &cm->mapSize, 1, PETSC_INT, 0);
662: PetscBinaryWrite(fd, &cm->numOldClasses, 1, PETSC_INT, 0);
663: if (cm->classMap != PETSC_NULL) {
664: PetscBinaryWrite(fd, &one, 1, PETSC_INT, 0);
665: for(bd = 0; bd < cm->mapSize; bd++) {
666: PetscBinaryWrite(fd, cm->classMap[bd], cm->numOldClasses, PETSC_INT, 0);
667: }
668: } else {
669: PetscBinaryWrite(fd, &zero, 1, PETSC_INT, 0);
670: }
671: } else {
672: /* Create the mesh context */
673: FieldClassMapCreate(comm, &cm);
674: PetscMemcpy(cm->ops, &cmOps, sizeof(FieldClassMapOps));
675: PetscStrallocpy(CLASS_MAP_TRIANGULAR_2D, &cm->type_name);
676: cm->data = PETSC_NULL;
678: PetscBinaryRead(fd, &cm->numFields, 1, PETSC_INT);
679: PetscMalloc(cm->numFields * sizeof(int), &cm->fields);
680: PetscBinaryRead(fd, cm->fields, cm->numFields, PETSC_INT);
681: PetscBinaryRead(fd, &cm->numNodes, 1, PETSC_INT);
682: PetscBinaryRead(fd, &cm->numGhostNodes, 1, PETSC_INT);
683: PetscBinaryRead(fd, &cm->numOverlapNodes, 1, PETSC_INT);
684: PetscBinaryRead(fd, &cm->numClasses, 1, PETSC_INT);
685: PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);
686: for(f = 0; f < cm->numFields; f++) {
687: PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);
688: PetscBinaryRead(fd, cm->fieldClasses[f], cm->numClasses, PETSC_INT);
689: }
690: PetscMalloc(cm->numOverlapNodes * sizeof(int), &cm->classes);
691: PetscBinaryRead(fd, cm->classes, cm->numOverlapNodes, PETSC_INT);
692: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizes);
693: PetscBinaryRead(fd, cm->classSizes, cm->numClasses, PETSC_INT);
694: PetscBinaryRead(fd, &cm->isReduced, 1, PETSC_INT);
695: if (cm->isReduced == PETSC_TRUE) {
696: PetscMalloc(cm->numFields * sizeof(int *), &cm->reduceFieldClasses);
697: for(f = 0; f < cm->numFields; f++) {
698: PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);
699: PetscBinaryRead(fd, cm->reduceFieldClasses[f], cm->numClasses, PETSC_INT);
700: }
701: }
702: PetscBinaryRead(fd, &cm->isConstrained, 1, PETSC_INT);
703: PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);
704: PetscBinaryRead(fd, cm->isClassConstrained, cm->numClasses, PETSC_INT);
705: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);
706: PetscBinaryRead(fd, cm->classSizeDiffs, cm->numClasses, PETSC_INT);
707: PetscBinaryRead(fd, &cm->mapSize, 1, PETSC_INT);
708: PetscBinaryRead(fd, &cm->numOldClasses, 1, PETSC_INT);
709: PetscBinaryRead(fd, &hasClassMap, 1, PETSC_INT);
710: if (hasClassMap) {
711: PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);
712: for(bd = 0; bd < cm->mapSize; bd++) {
713: PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);
714: PetscBinaryRead(fd, cm->classMap[bd], cm->numOldClasses, PETSC_INT);
715: }
716: }
718: PetscLogObjectMemory(cm, (cm->numFields + cm->mapSize) * sizeof(int *) + (cm->numFields*(cm->numClasses + 1) +
719: cm->numOverlapNodes + cm->numOldClasses*cm->mapSize + cm->numClasses*3) * sizeof(int));
720: *map = cm;
721: }
722: return(0);
723: }
724: EXTERN_C_END
726: EXTERN_C_BEGIN
727: #undef __FUNCT__
729: int FieldClassMapCreate_Triangular_2D(FieldClassMap cm)
730: {
731: ParameterDict dict;
732: Grid grid;
733: Mesh mesh;
734: Partition part;
735: Discretization disc, firstDisc;
736: int *fields;
737: int numOverlapElements, numCorners;
738: int f, field, comp, elem, corner, node, nclass;
739: PetscTruth match, isconst, islin, isquad;
740: int ierr;
743: PetscMemcpy(cm->ops, &cmOps, sizeof(FieldClassMapOps));
744: PetscStrallocpy(CLASS_MAP_SER_TRIANGULAR_2D_BINARY, &cm->serialize_name);
745: /* Get arguments */
746: PetscObjectGetParameterDict((PetscObject) cm, &dict);
747: ParameterDictGetInteger(dict, "numFields", &cm->numFields);
748: ParameterDictGetObject(dict, "fields", (void **) &fields);
749: ParameterDictGetObject(dict, "grid", (void **) &grid);
750: /* Get fields */
751: PetscMalloc(cm->numFields * sizeof(int), &cm->fields);
752: PetscMemcpy(cm->fields, fields, cm->numFields * sizeof(int));
754: /* Get number of classes -- this is very simplistic, needs to be changed for other discretizations */
755: GridGetFieldDisc(grid, 0, &firstDisc);
756: cm->numClasses = 1;
757: for(f = 0; f < cm->numFields; f++) {
758: GridGetFieldDisc(grid, cm->fields[f], &disc);
759: PetscStrcmp(disc->type_name, firstDisc->type_name, &match);
760: if (match == PETSC_FALSE) {
761: cm->numClasses++;
762: break;
763: }
764: }
765: GridGetMesh(grid, &mesh);
766: MeshGetPartition(mesh, &part);
767: PartitionGetNumNodes(part, &cm->numNodes);
768: PartitionGetNumOverlapNodes(part, &cm->numOverlapNodes);
769: cm->numGhostNodes = cm->numOverlapNodes - cm->numNodes;
770: PetscMalloc(cm->numOverlapNodes * sizeof(int), &cm->classes);
772: /* Setup field class membership -- this is very simplistic, needs to be changed for other discretizations */
773: PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);
774: for(f = 0; f < cm->numFields; f++) {
775: field = cm->fields[f];
776: PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);
777: GridGetFieldDisc(grid, field, &disc);
778: for(nclass = 0; nclass < cm->numClasses; nclass++) {
779: if (grid->dim == 1) {
780: PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_1D_CONSTANT, &isconst);
781: PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_1D_LINEAR, &islin);
782: PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_1D_QUADRATIC, &isquad);
783: } else if (grid->dim == 2) {
784: PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_2D_LINEAR, &islin);
785: PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_2D_QUADRATIC, &isquad);
786: }
787: /* These only work for pairs in 2D (all are fine in 1D, how do I generalize this?) */
788: if (isconst == PETSC_TRUE) {
789: cm->fieldClasses[f][nclass] = (nclass == 1 ? 1 : 0);
790: } else if (islin == PETSC_TRUE) {
791: cm->fieldClasses[f][nclass] = (nclass == 0 ? 1 : 0);
792: } else if (isquad == PETSC_TRUE) {
793: cm->fieldClasses[f][nclass] = 1;
794: } else {
795: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid discretization type");
796: }
797: }
798: }
800: /* Setup classes */
801: PetscMemzero(cm->classes, cm->numOverlapNodes * sizeof(int));
802: if (cm->numClasses > 1) {
803: /* Set all midnodes to class 1 */
804: PartitionGetNumOverlapElements(part, &numOverlapElements);
805: MeshGetNumCorners(mesh, &numCorners);
806: for(elem = 0; elem < numOverlapElements; elem++) {
807: for(corner = grid->dim+1; corner < numCorners; corner++) {
808: MeshGetNodeFromElement(mesh, elem, corner, &node);
809: if (cm->classes[node] == 0)
810: cm->classes[node] = 1;
811: }
812: }
813: }
815: /* Communicate ghost classes */
816: PartitionGhostNodeExchange(part, INSERT_VALUES, SCATTER_FORWARD, cm->classes, &cm->classes[cm->numNodes]);
817:
819: /* Setup class sizes */
820: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizes);
821: PetscMemzero(cm->classSizes, cm->numClasses * sizeof(int));
822: for(nclass = 0; nclass < cm->numClasses; nclass++) {
823: for(f = 0; f < cm->numFields; f++) {
824: field = cm->fields[f];
825: if (cm->fieldClasses[f][nclass]) {
826: GridGetFieldComponents(grid, field, &comp);
827: cm->classSizes[nclass] += comp;
828: }
829: }
830: }
832: /* Initialize constraint support */
833: PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);
834: PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);
835: PetscMemzero(cm->isClassConstrained, cm->numClasses * sizeof(int));
836: PetscMemzero(cm->classSizeDiffs, cm->numClasses * sizeof(int));
838: PetscLogObjectMemory(cm, (cm->numFields + cm->numOverlapNodes + (cm->numFields + 4)*cm->numClasses) * sizeof(int) +
839: cm->numFields * sizeof(int *));
840: return(0);
841: }
842: EXTERN_C_END