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