Actual source code: fieldClassMap.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: fieldClassMap.c,v 1.6 2000/01/31 17:42:32 knepley Exp $";
  3: #endif

 5:  #include src/grid/gridimpl.h

  7: /* Logging support */
  8: int CLASS_MAP_COOKIE;

 10: #undef  __FUNCT__
 12: /*@
 13:   FieldClassMapDestroy - Destroys a class map object.

 15:   Collective on FieldClassMap

 17:   Input Parameter:
 18: . map - The map

 20:   Level: beginner

 22: .keywords: class, field class, class map, destroy
 23: .seealso FieldClassMapCreate()
 24: @*/
 25: int FieldClassMapDestroy(FieldClassMap map)
 26: {
 27:   int f, bd;

 32:   if (--map->refct > 0)
 33:     return(0);
 34:   PetscFree(map->fields);
 35:   for(f = 0; f < map->numFields; f++) {
 36:     PetscFree(map->fieldClasses[f]);
 37:   }
 38:   PetscFree(map->fieldClasses);
 39:   PetscFree(map->classes);
 40:   PetscFree(map->classSizes);
 41:   if (map->isReduced == PETSC_TRUE) {
 42:     for(f = 0; f < map->numFields; f++) {
 43:       PetscFree(map->reduceFieldClasses[f]);
 44:     }
 45:     PetscFree(map->reduceFieldClasses);
 46:   }
 47:   PetscFree(map->isClassConstrained);
 48:   PetscFree(map->classSizeDiffs);
 49:   if (map->classMap != PETSC_NULL) {
 50:     for(bd = 0; bd < map->mapSize; bd++) {
 51:       PetscFree(map->classMap[bd]);
 52:     }
 53:     PetscFree(map->classMap);
 54:   }
 55:   (*map->ops->destroy)(map);
 56:   PetscLogObjectDestroy(map);
 57:   PetscHeaderDestroy(map);
 58:   return(0);
 59: }

 61: #undef  __FUNCT__
 63: /*@
 64:   FieldClassMapView - Views a class map object.

 66:   Collective on FieldClassMap

 68:   Input Parameters:
 69: + map    - The map
 70: - viewer - The viewer with which to view the map

 72:   Level: beginner

 74: .keywords: class, field class, class map, view
 75: .seealso: FieldClassMapDestroy(), PetscViewerDrawOpen()
 76: @*/
 77: int FieldClassMapView(FieldClassMap map, PetscViewer viewer)
 78: {

 83:   if (viewer == PETSC_NULL)
 84:     viewer = PETSC_VIEWER_STDOUT_SELF;
 85:   else
 87:   (*map->ops->view)(map, viewer);
 88:   return(0);
 89: }

 91: /*MC
 92:   FieldClassMapRegister - Adds a creation method to the class map package.

 94:   Synopsis:

 96:   FieldClassMapRegister(char *name, char *path, char *func_name, int (*create_func)(FieldClassMap))

 98:   Not Collective

100:   Input Parameters:
101: + name        - The name of a new user-defined creation routine
102: . path        - The path (either absolute or relative) of the library containing this routine
103: . func_name   - The name of routine to create method context
104: - create_func - The creation routine itself

106:   Notes:
107:   FieldClassMapRegister() may be called multiple times to add several user-defined mapes.

109:   If dynamic libraries are used, then the fourth input argument (create_func) is ignored.

111:   Sample usage:
112: .vb
113:   FieldClassMapRegister("my_map", /home/username/my_lib/lib/libO/solaris/mylib.a, "MyFunc", MyFunc);
114: .ve

116:   Then, your map type can be chosen with the procedural interface via
117: $     FieldClassMapSetType(vec, "my_map")
118:   or at runtime via the option
119: $     -class_map_type my_map

121:   Level: advanced

123:   $PETSC_ARCH and $BOPT occuring in pathname will be replaced with appropriate values.

125: .keywords: class, field class, class map, register
126: .seealso: FieldClassMapRegisterAll(), FieldClassMapRegisterDestroy()
127: M*/
128: #undef __FUNCT__  
130: int FieldClassMapRegister_Private(const char sname[], const char path[], const char name[],
131:                                   int (*function)(FieldClassMap))
132: {
133:   char fullname[256];
134:   int  ierr;

137:   PetscStrcpy(fullname, path);
138:   PetscStrcat(fullname, ":");
139:   PetscStrcat(fullname, name);
140:   PetscFListAdd(&FieldClassMapList, sname, fullname, (void (*)(void)) function);
141:   return(0);
142: }

144: /*MC
145:   FieldClassMapSerializeRegister - Adds a serialization method to the class map package.

147:   Synopsis:

149:   FieldClassMapSerializeRegister(char *serialize_name, char *path, char *serialize_func_name,
150:                                  int (*serialize_func)(MPI_Comm, FieldClassMap *, PetscViewer, PetscTruth))

152:   Not Collective

154:   Input Parameters:
155: + serialize_name      - The name of a new user-defined serialization routine
156: . path                - The path (either absolute or relative) of the library containing this routine
157: . serialize_func_name - The name of routine to create method context
158: - serialize_func      - The serialization routine itself

160:   Notes:
161:   FieldClassMapSerializeRegister() may be called multiple times to add several user-defined solvers.

163:   If dynamic libraries are used, then the fourth input argument (serialize_func) is ignored.

165:   Sample usage:
166: .vb
167:   FieldClassMapSerializeRegister("my_store", /home/username/my_lib/lib/libO/solaris/mylib.a, "MyStoreFunc", MyStoreFunc);
168: .ve

170:   Then, your serialization can be chosen with the procedural interface via
171: $     FieldClassMapSetSerializeType(map, "my_store")
172:   or at runtime via the option
173: $     -class_map_serialize_type my_store

175:   Level: advanced

177:   $PETSC_ARCH and $BOPT occuring in pathname will be replaced with appropriate values.

179: .keywords: class, field class, class map, register
180: .seealso: FieldClassMapSerializeRegisterAll(), FieldClassMapSerializeRegisterDestroy()
181: M*/
182: #undef __FUNCT__  
184: int FieldClassMapSerializeRegister_Private(const char sname[], const char path[], const char name[],
185:                                            int (*function)(MPI_Comm, FieldClassMap *, PetscViewer, PetscTruth))
186: {
187:   char fullname[256];
188:   int  ierr;

191:   PetscStrcpy(fullname, path);
192:   PetscStrcat(fullname, ":");
193:   PetscStrcat(fullname, name);
194:   PetscFListAdd(&FieldClassMapSerializeList, sname, fullname, (void (*)(void)) function);
195:   return(0);
196: }

198: #undef  __FUNCT__
200: /*@
201:   FieldClassMapCreateSubset - This function creates a new class map from an existing
202:   one on a subset of the fields.
203:   
204:   Collective on FieldClassMap

206:   Input Parameter:
207: . map    - The map

209:   Output Parameter:
210: . newMap - The new map

212:   Level: beginner

214: .keywords: class, field class, class map, subset
215: .seealso FieldClassMapCreate()
216: @*/
217: int FieldClassMapCreateSubset(FieldClassMap map, int numFields, int *fields, FieldClassMap *newMap)
218: {
219:   FieldClassMap cm;
220:   int           f, g, bd;
221:   int           ierr;

226:   /* Check for consistency */
227:   for(f = 0; f < numFields; f++) {
228:     for(g = 0; g < map->numFields; g++) {
229:       if (map->fields[g] == fields[f])
230:         break;
231:     }
232:     if (g == map->numFields) {
233:       SETERRQ1(PETSC_ERR_ARG_INCOMP, "Fields not a subset: field %d not in the existing order", fields[f]);
234:     }
235:   }
236:   /* Create new map */
237:   FieldClassMapCreate(map->comm, &cm);
238:   PetscMemcpy(cm->ops, map->ops, sizeof(FieldClassMapOps));
239:   PetscStrallocpy(map->type_name, &cm->type_name);
240:   cm->data = PETSC_NULL;

242:   cm->numFields        = numFields;
243:   PetscMalloc(cm->numFields * sizeof(int), &cm->fields);
244:   PetscMemcpy(cm->fields, fields, cm->numFields * sizeof(int));
245:   cm->numNodes         = map->numNodes;
246:   cm->numGhostNodes    = map->numGhostNodes;
247:   cm->numOverlapNodes  = map->numOverlapNodes;
248:   cm->numClasses       = map->numClasses;
249:   PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);
250:   for(f = 0; f < cm->numFields; f++) {
251:     for(g = 0; g < map->numFields; g++)
252:       if (cm->fields[f] == map->fields[g])
253:         break;
254:     if (g == map->numFields) {
255:       SETERRQ1(PETSC_ERR_PLIB, "Unable to locate field %d in superset", cm->fields[f]);
256:     }
257:     PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);
258:     PetscMemcpy(cm->fieldClasses[f], map->fieldClasses[g], cm->numClasses * sizeof(int));
259:   }
260:   PetscMalloc(cm->numOverlapNodes  * sizeof(int), &cm->classes);
261:   PetscMemcpy(cm->classes, map->classes, cm->numOverlapNodes * sizeof(int));
262:   PetscMalloc(cm->numClasses       * sizeof(int), &cm->classSizes);
263:   PetscMemcpy(cm->classSizes, map->classSizes, cm->numClasses * sizeof(int));
264: #if 0
265:   for(nclass = 0; nclass < cm->numClasses; nclass++) {
266:     for(g = 0; g < map->numFields; g++) {
267:       for(f = 0; f < cm->numFields; f++)
268:         if (cm->fields[f] == map->fields[g])
269:           break;
270:       if (f < cm->numFields) continue;
271:       /* Subtract size of missing fields */
272:       if (map->fieldClasses[g][nclass]) {
273:         GridGetFieldComponents(grid, map->fields[g], &comp);
274:         cm->classSizes[nclass] -= comp;
275:       }
276:     }
277:   }
278: #endif
279:   cm->isReduced = map->isReduced;
280:   if (map->isReduced == PETSC_TRUE) {
281:     PetscMalloc(cm->numFields * sizeof(int *), &cm->reduceFieldClasses);
282:     for(f = 0; f < cm->numFields; f++) {
283:       for(g = 0; g < map->numFields; g++) {
284:         if (cm->fields[f] == map->fields[g])
285:           break;
286:       }
287:       if (g == map->numFields) SETERRQ1(PETSC_ERR_PLIB, "Unable to locate field %d in superset", cm->fields[f]);
288:       PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);
289:       PetscMemcpy(cm->reduceFieldClasses[f], map->reduceFieldClasses[g], cm->numClasses * sizeof(int));
290:     }
291:   }
292:   cm->isConstrained = map->isConstrained;
293:   PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);
294:   PetscMemcpy(cm->isClassConstrained, map->isClassConstrained, cm->numClasses * sizeof(int));
295:   PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);
296:   PetscMemcpy(cm->classSizeDiffs, map->classSizeDiffs, cm->numClasses * sizeof(int));
297:   cm->mapSize       = map->mapSize;
298:   cm->numOldClasses = map->numOldClasses;
299:   if (map->classMap != PETSC_NULL) {
300:     PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);
301:     for(bd = 0; bd < cm->mapSize; bd++) {
302:       PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);
303:       PetscMemcpy(cm->classMap[bd], map->classMap[bd], cm->numOldClasses * sizeof(int));
304:     }
305:   }

307:   PetscLogObjectMemory(cm, (cm->numFields + cm->mapSize) * sizeof(int *) + (cm->numFields*(cm->numClasses + 1) +
308:                    cm->numOverlapNodes + cm->numClasses*3 + cm->numOldClasses*cm->mapSize) * sizeof(int));
309:   *newMap = cm;
310:   return(0);
311: }

313: #undef  __FUNCT__
315: /*@
316:   FieldClassMapDuplicate - Duplicates a class map object.
317:   
318:   Collective on FieldClassMap

320:   Input Parameter:
321: . map    - The map

323:   Output Parameter:
324: . newMap - The new map

326:   Level: beginner

328: .keywords: class, field class, class map, duplicate
329: .seealso FieldClassMapCreate()
330: @*/
331: int FieldClassMapDuplicate(FieldClassMap map, FieldClassMap *newMap)
332: {
333:   FieldClassMap cm;
334:   int           f, bd;
335:   int           ierr;

340:   FieldClassMapCreate(map->comm, &cm);
341:   PetscMemcpy(cm->ops, map->ops, sizeof(FieldClassMapOps));
342:   PetscStrallocpy(map->type_name, &cm->type_name);
343:   cm->data = PETSC_NULL;

345:   cm->numFields        = map->numFields;
346:   PetscMalloc(cm->numFields * sizeof(int), &cm->fields);
347:   PetscMemcpy(cm->fields, map->fields, cm->numFields * sizeof(int));
348:   cm->numNodes         = map->numNodes;
349:   cm->numGhostNodes    = map->numGhostNodes;
350:   cm->numOverlapNodes  = map->numOverlapNodes;
351:   cm->numClasses       = map->numClasses;
352:   PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);
353:   for(f = 0; f < cm->numFields; f++) {
354:     PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);
355:     PetscMemcpy(cm->fieldClasses[f], map->fieldClasses[f], cm->numClasses * sizeof(int));
356:   }
357:   PetscMalloc(cm->numOverlapNodes * sizeof(int), &cm->classes);
358:   PetscMemcpy(cm->classes, map->classes, cm->numOverlapNodes * sizeof(int));
359:   PetscMalloc(cm->numClasses      * sizeof(int), &cm->classSizes);
360:   PetscMemcpy(cm->classSizes, map->classSizes, cm->numClasses * sizeof(int));
361:   cm->isReduced        = map->isReduced;
362:   if (map->isReduced == PETSC_TRUE) {
363:     PetscMalloc(cm->numFields * sizeof(int *), &cm->reduceFieldClasses);
364:     for(f = 0; f < cm->numFields; f++) {
365:       PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);
366:       PetscMemcpy(cm->reduceFieldClasses[f], map->reduceFieldClasses[f], cm->numClasses * sizeof(int));
367:     }
368:   }
369:   cm->isConstrained = map->isConstrained;
370:   PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);
371:   PetscMemcpy(cm->isClassConstrained, map->isClassConstrained, cm->numClasses * sizeof(int));
372:   PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);
373:   PetscMemcpy(cm->classSizeDiffs, map->classSizeDiffs, cm->numClasses * sizeof(int));
374:   cm->mapSize          = map->mapSize;
375:   cm->numOldClasses    = map->numOldClasses;
376:   if (map->classMap != PETSC_NULL) {
377:     PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);
378:     for(bd = 0; bd < cm->mapSize; bd++) {
379:       PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);
380:       PetscMemcpy(cm->classMap[bd], map->classMap[bd], cm->numOldClasses * sizeof(int));
381:     }
382:   }

384:   PetscLogObjectMemory(cm, (cm->numFields + cm->mapSize) * sizeof(int *) + (cm->numFields*(cm->numClasses + 1) +
385:                        cm->numOverlapNodes + cm->numClasses*3 + cm->numOldClasses*cm->mapSize) * sizeof(int));
386:   *newMap = cm;
387:   return(0);
388: }

390: #undef  __FUNCT__
392: /*@
393:   FieldClassMapConstrain - This function creates a new class map from an existing one by implementing
394:   constraints and boundary conditions registered with the grid.

396:   Collective on FieldClassMap

398:   Input Parameters:
399: + map      - The map
400: . grid     - The grid
401: . useBC    - The flag for reducing boundary conditions
402: - useConst - The flag for using constraints

404:   Output Parameter:
405: . newMap   - The constrained map

407:   Level: advanced

409: .keywords: class, field class, class map, constraint
410: .seealso FieldClassMapCreate()
411: @*/
412: int FieldClassMapConstrain(FieldClassMap map, Grid grid, PetscTruth useBC, PetscTruth useConst, FieldClassMap *newMap)
413: {

419:   (*map->ops->constrain)(map, grid, useBC, useConst, newMap);
420:   return(0);
421: }

423: #undef  __FUNCT__
425: /*@
426:   FieldClassMapReduce - This function creates a new class map from an existing one which
427:   lies in the space of variables reduced by boundary conditions.

429:   Collective on FieldClassMap

431:   Input Parameters:
432: + map    - The map
433: - grid   - The grid

435:   Output Parameter:
436: . newMap - The reduction map

438:   Note:
439:   This function is normally used to create operators which map elements from the space of
440:   boundary values to the calculation space. Thus you can construct a matrix to add boundary
441:   conditions to the rhs.

443:   Level: advanced

445: .keywords: class, field class, class map, reduction
446: .seealso: FieldClassMapCreate(), FieldClassMapConstrain()
447: @*/
448: int FieldClassMapReduce(FieldClassMap map, Grid grid, FieldClassMap *newMap)
449: {

455:   (*map->ops->reduce)(map, grid, newMap);
456:   return(0);
457: }

459: #undef  __FUNCT__
461: /*@
462:   FieldClassMapIsConstrained - This function determines whether or not the map is constrained.

464:   Not collective

466:   Input Parameter:
467: . cm            - The class map

469:   Output Parameter:
470: . isConstrained - The constraint flag

472:   Level: intermediate

474: .keywords: class, field class, class map
475: .seealso: FieldClassMapCreate(), FieldClassMapConstrain()
476: @*/
477: int FieldClassMapIsConstrained(FieldClassMap cm, PetscTruth *isConstrained) {
481:   *isConstrained = cm->isConstrained;
482:   return(0);
483: }

485: #undef  __FUNCT__
487: /*@
488:   FieldClassMapIsDefined - This function determines whether a field is defined on
489:   nodes of a given class, or equivalently whether the field is a member of the class.

491:   Not collective

493:   Input Parameters:
494: + map     - The map
495: . field   - The field
496: - nclass  - The node class

498:   Output Parameter:
499: . defined - The flag for class membership

501:   Level: intermediate

503: .keywords: class, field class, class map
504: .seealso: FieldClassMapCreate(), FieldClassMapConstrain()
505: @*/
506: int FieldClassMapIsDefined(FieldClassMap map, int field, int nclass, PetscTruth *defined)
507: {
508:   int   numFields    = map->numFields;
509:   int   numClasses   = map->numClasses;
510:   int  *fields       = map->fields;
511:   int **fieldClasses = map->fieldClasses;
512:   int   f;

517:   for(f = 0; f < numFields; f++) {
518:     if (field == fields[f]) break;
519:   }
520:   if (f == numFields) SETERRQ1(PETSC_ERR_ARG_WRONG, "Field %d not present in grid", field);
521:   if ((nclass < 0) || (nclass >= numClasses)) {
522:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid class %d should be in [0,%d)", nclass, numClasses);
523:   }
524:   if (fieldClasses[f][nclass]) {
525:     *defined = PETSC_TRUE;
526:   } else {
527:     *defined = PETSC_FALSE;
528:   }
529:   return(0);
530: }

532: #undef  __FUNCT__
534: /*@
535:   FieldClassMapGetNumFields - This function returns the number of fields in the map.

537:   Not collective

539:   Input Parameter:
540: . cm        - The class map

542:   Output Parameter:
543: . numFields - The number of fields in the map

545:   Level: intermediate

547: .keywords: class, class map, field
548: .seealso: FieldClassMapGetField(), FieldClassMapCreate()
549: @*/
550: int FieldClassMapGetNumFields(FieldClassMap cm, int *numFields) {
554:   *numFields = cm->numFields;
555:   return(0);
556: }

558: #undef  __FUNCT__
560: /*@
561:   FieldClassMapGetField - This function returns the canonical field number of a field index in the map.

563:   Not collective

565:   Input Parameters:
566: + cm    - The class map
567: - f     - The field index in the map

569:   Output Parameter:
570: . field - The canonical field number

572:   Level: intermediate

574: .keywords: class, class map, field
575: .seealso: FieldClassMapGetNumFields(), FieldClassMapCreate()
576: @*/
577: int FieldClassMapGetField(FieldClassMap cm, int f, int *field) {
581:   if ((f < 0) || (f >= cm->numFields)) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Field index %d must be in [0,%d)", f, cm->numFields);
582:   *field = cm->fields[f];
583:   return(0);
584: }

586: #undef  __FUNCT__
588: /*@
589:   FieldClassMapGetNodeClass - This function returns the class for the given node.

591:   Not collective

593:   Input Parameters:
594: + map    - The map
595: - node   - The node

597:   Output Parameter:
598: . nclass - The node class

600:   Level: intermediate

602: .keywords: class, class map
603: .seealso: FieldClassMapCreate(), FieldClassMapConstrain()
604: @*/
605: int FieldClassMapGetNodeClass(FieldClassMap map, int node, int *nclass)
606: {
607:   int  numOverlapNodes = map->numOverlapNodes;
608:   int *classes         = map->classes;

613:   if ((node < 0) || (node >= numOverlapNodes)) {
614:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", node, numOverlapNodes);
615:   }
616:   *nclass = classes[node];
617:   return(0);
618: }