Actual source code: ex23.c

  1: static char help[] = "Test for function and field projection\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscds.h>

  6: typedef struct {
  7:   PetscBool multifield; /* Different numbers of input and output fields */
  8:   PetscBool subdomain;  /* Try with a volumetric submesh */
  9:   PetscBool submesh;    /* Try with a boundary submesh */
 10:   PetscBool auxfield;   /* Try with auxiliary fields */
 11: } AppCtx;

 13: /* (x + y)*dim + d */
 14: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 15: {
 16:   PetscInt c;
 17:   for (c = 0; c < Nc; ++c) u[c] = (x[0] + x[1]) * Nc + c;
 18:   return 0;
 19: }

 21: /* {x, y, z} */
 22: static PetscErrorCode linear2(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 23: {
 24:   PetscInt c;
 25:   for (c = 0; c < Nc; ++c) u[c] = x[c];
 26:   return 0;
 27: }

 29: /* {u_x, u_y, u_z} */
 30: static void linear_vector(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
 31: {
 32:   PetscInt d;
 33:   for (d = 0; d < uOff[1] - uOff[0]; ++d) f[d] = u[d + uOff[0]];
 34: }

 36: /* p */
 37: static void linear_scalar(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
 38: {
 39:   f[0] = u[uOff[1]];
 40: }

 42: /* {div u, p^2} */
 43: static void divergence_sq(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
 44: {
 45:   PetscInt d;
 46:   f[0] = 0.0;
 47:   for (d = 0; d < dim; ++d) f[0] += u_x[uOff_x[0] + d * dim + d];
 48:   f[1] = PetscSqr(u[uOff[1]]);
 49: }

 51: static PetscErrorCode ProcessOptions(AppCtx *options)
 52: {
 53:   options->multifield = PETSC_FALSE;
 54:   options->subdomain  = PETSC_FALSE;
 55:   options->submesh    = PETSC_FALSE;
 56:   options->auxfield   = PETSC_FALSE;

 58:   PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
 59:   PetscOptionsBool("-multifield", "Flag for trying different numbers of input/output fields", "ex23.c", options->multifield, &options->multifield, NULL);
 60:   PetscOptionsBool("-subdomain", "Flag for trying volumetric submesh", "ex23.c", options->subdomain, &options->subdomain, NULL);
 61:   PetscOptionsBool("-submesh", "Flag for trying boundary submesh", "ex23.c", options->submesh, &options->submesh, NULL);
 62:   PetscOptionsBool("-auxfield", "Flag for trying auxiliary fields", "ex23.c", options->auxfield, &options->auxfield, NULL);
 63:   PetscOptionsEnd();
 64:   return 0;
 65: }

 67: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 68: {
 69:   DMCreate(comm, dm);
 70:   DMSetType(*dm, DMPLEX);
 71:   DMSetFromOptions(*dm);
 72:   DMViewFromOptions(*dm, NULL, "-orig_dm_view");
 73:   return 0;
 74: }

 76: static PetscErrorCode SetupDiscretization(DM dm, PetscInt dim, PetscBool simplex, AppCtx *user)
 77: {
 78:   PetscFE  fe;
 79:   MPI_Comm comm;

 82:   PetscObjectGetComm((PetscObject)dm, &comm);
 83:   PetscFECreateDefault(comm, dim, dim, simplex, "velocity_", -1, &fe);
 84:   PetscObjectSetName((PetscObject)fe, "velocity");
 85:   DMSetField(dm, 0, NULL, (PetscObject)fe);
 86:   PetscFEDestroy(&fe);
 87:   PetscFECreateDefault(comm, dim, 1, simplex, "pressure_", -1, &fe);
 88:   PetscObjectSetName((PetscObject)fe, "pressure");
 89:   DMSetField(dm, 1, NULL, (PetscObject)fe);
 90:   PetscFEDestroy(&fe);
 91:   DMCreateDS(dm);
 92:   return 0;
 93: }

 95: static PetscErrorCode SetupOutputDiscretization(DM dm, PetscInt dim, PetscBool simplex, AppCtx *user)
 96: {
 97:   PetscFE  fe;
 98:   MPI_Comm comm;

101:   PetscObjectGetComm((PetscObject)dm, &comm);
102:   PetscFECreateDefault(comm, dim, dim, simplex, "output_", -1, &fe);
103:   PetscObjectSetName((PetscObject)fe, "output");
104:   DMSetField(dm, 0, NULL, (PetscObject)fe);
105:   PetscFEDestroy(&fe);
106:   DMCreateDS(dm);
107:   return 0;
108: }

110: static PetscErrorCode CreateSubdomainMesh(DM dm, DMLabel *domLabel, DM *subdm, AppCtx *user)
111: {
112:   DMLabel   label;
113:   PetscBool simplex;
114:   PetscInt  dim, cStart, cEnd, c;

117:   DMPlexIsSimplex(dm, &simplex);
118:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
119:   DMLabelCreate(PETSC_COMM_SELF, "subdomain", &label);
120:   for (c = cStart + (cEnd - cStart) / 2; c < cEnd; ++c) DMLabelSetValue(label, c, 1);
121:   DMPlexFilter(dm, label, 1, subdm);
122:   DMGetDimension(*subdm, &dim);
123:   SetupDiscretization(*subdm, dim, simplex, user);
124:   PetscObjectSetName((PetscObject)*subdm, "subdomain");
125:   DMViewFromOptions(*subdm, NULL, "-sub_dm_view");
126:   if (domLabel) *domLabel = label;
127:   else DMLabelDestroy(&label);
128:   return 0;
129: }

131: static PetscErrorCode CreateBoundaryMesh(DM dm, DMLabel *bdLabel, DM *subdm, AppCtx *user)
132: {
133:   DMLabel   label;
134:   PetscBool simplex;
135:   PetscInt  dim;

138:   DMPlexIsSimplex(dm, &simplex);
139:   DMLabelCreate(PETSC_COMM_SELF, "sub", &label);
140:   DMPlexMarkBoundaryFaces(dm, 1, label);
141:   DMPlexLabelComplete(dm, label);
142:   DMPlexCreateSubmesh(dm, label, 1, PETSC_TRUE, subdm);
143:   DMGetDimension(*subdm, &dim);
144:   SetupDiscretization(*subdm, dim, simplex, user);
145:   PetscObjectSetName((PetscObject)*subdm, "boundary");
146:   DMViewFromOptions(*subdm, NULL, "-sub_dm_view");
147:   if (bdLabel) *bdLabel = label;
148:   else DMLabelDestroy(&label);
149:   return 0;
150: }

152: static PetscErrorCode CreateAuxiliaryVec(DM dm, DM *auxdm, Vec *la, AppCtx *user)
153: {
154:   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
155:   PetscBool simplex;
156:   PetscInt  dim, Nf, f;

159:   DMGetDimension(dm, &dim);
160:   DMPlexIsSimplex(dm, &simplex);
161:   DMGetNumFields(dm, &Nf);
162:   PetscMalloc1(Nf, &afuncs);
163:   for (f = 0; f < Nf; ++f) afuncs[f] = linear;
164:   DMClone(dm, auxdm);
165:   SetupDiscretization(*auxdm, dim, simplex, user);
166:   DMCreateLocalVector(*auxdm, la);
167:   DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, *la);
168:   VecViewFromOptions(*la, NULL, "-local_aux_view");
169:   PetscFree(afuncs);
170:   return 0;
171: }

173: static PetscErrorCode TestFunctionProjection(DM dm, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
174: {
175:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
176:   Vec      x, lx;
177:   PetscInt Nf, f;
178:   PetscInt val[1] = {1};
179:   char     lname[PETSC_MAX_PATH_LEN];

182:   if (dmAux) DMSetAuxiliaryVec(dm, NULL, 0, 0, la);
183:   DMGetNumFields(dm, &Nf);
184:   PetscMalloc1(Nf, &funcs);
185:   for (f = 0; f < Nf; ++f) funcs[f] = linear;
186:   DMGetGlobalVector(dm, &x);
187:   PetscStrcpy(lname, "Function ");
188:   PetscStrcat(lname, name);
189:   PetscObjectSetName((PetscObject)x, lname);
190:   if (!label) DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, x);
191:   else DMProjectFunctionLabel(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, x);
192:   VecViewFromOptions(x, NULL, "-func_view");
193:   DMRestoreGlobalVector(dm, &x);
194:   DMGetLocalVector(dm, &lx);
195:   PetscStrcpy(lname, "Local Function ");
196:   PetscStrcat(lname, name);
197:   PetscObjectSetName((PetscObject)lx, lname);
198:   if (!label) DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_VALUES, lx);
199:   else DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, lx);
200:   VecViewFromOptions(lx, NULL, "-local_func_view");
201:   DMRestoreLocalVector(dm, &lx);
202:   PetscFree(funcs);
203:   if (dmAux) DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL);
204:   return 0;
205: }

207: static PetscErrorCode TestFieldProjection(DM dm, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
208: {
209:   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
210:   void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
211:   Vec      lx, lu;
212:   PetscInt Nf, f;
213:   PetscInt val[1] = {1};
214:   char     lname[PETSC_MAX_PATH_LEN];

217:   if (dmAux) DMSetAuxiliaryVec(dm, NULL, 0, 0, la);
218:   DMGetNumFields(dm, &Nf);
219:   PetscMalloc2(Nf, &funcs, Nf, &afuncs);
220:   for (f = 0; f < Nf; ++f) afuncs[f] = linear;
221:   funcs[0] = linear_vector;
222:   funcs[1] = linear_scalar;
223:   DMGetLocalVector(dm, &lu);
224:   PetscStrcpy(lname, "Local Field Input ");
225:   PetscStrcat(lname, name);
226:   PetscObjectSetName((PetscObject)lu, lname);
227:   if (!label) DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, lu);
228:   else DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu);
229:   VecViewFromOptions(lu, NULL, "-local_input_view");
230:   DMGetLocalVector(dm, &lx);
231:   PetscStrcpy(lname, "Local Field ");
232:   PetscStrcat(lname, name);
233:   PetscObjectSetName((PetscObject)lx, lname);
234:   if (!label) DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx);
235:   else DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx);
236:   VecViewFromOptions(lx, NULL, "-local_field_view");
237:   DMRestoreLocalVector(dm, &lx);
238:   DMRestoreLocalVector(dm, &lu);
239:   PetscFree2(funcs, afuncs);
240:   if (dmAux) DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL);
241:   return 0;
242: }

244: static PetscErrorCode TestFieldProjectionMultiple(DM dm, DM dmIn, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
245: {
246:   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
247:   void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
248:   Vec      lx, lu;
249:   PetscInt Nf, NfIn;
250:   PetscInt val[1] = {1};
251:   char     lname[PETSC_MAX_PATH_LEN];

254:   if (dmAux) DMSetAuxiliaryVec(dm, NULL, 0, 0, la);
255:   DMGetNumFields(dm, &Nf);
256:   DMGetNumFields(dmIn, &NfIn);
257:   PetscMalloc2(Nf, &funcs, NfIn, &afuncs);
258:   funcs[0]  = divergence_sq;
259:   afuncs[0] = linear2;
260:   afuncs[1] = linear;
261:   DMGetLocalVector(dmIn, &lu);
262:   PetscStrcpy(lname, "Local MultiField Input ");
263:   PetscStrcat(lname, name);
264:   PetscObjectSetName((PetscObject)lu, lname);
265:   if (!label) DMProjectFunctionLocal(dmIn, 0.0, afuncs, NULL, INSERT_VALUES, lu);
266:   else DMProjectFunctionLabelLocal(dmIn, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu);
267:   VecViewFromOptions(lu, NULL, "-local_input_view");
268:   DMGetLocalVector(dm, &lx);
269:   PetscStrcpy(lname, "Local MultiField ");
270:   PetscStrcat(lname, name);
271:   PetscObjectSetName((PetscObject)lx, lname);
272:   if (!label) DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx);
273:   else DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx);
274:   VecViewFromOptions(lx, NULL, "-local_field_view");
275:   DMRestoreLocalVector(dm, &lx);
276:   DMRestoreLocalVector(dmIn, &lu);
277:   PetscFree2(funcs, afuncs);
278:   if (dmAux) DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL);
279:   return 0;
280: }

282: int main(int argc, char **argv)
283: {
284:   DM        dm, subdm, auxdm;
285:   Vec       la;
286:   PetscInt  dim;
287:   PetscBool simplex;
288:   AppCtx    user;

291:   PetscInitialize(&argc, &argv, NULL, help);
292:   ProcessOptions(&user);
293:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
294:   DMGetDimension(dm, &dim);
295:   DMPlexIsSimplex(dm, &simplex);
296:   SetupDiscretization(dm, dim, simplex, &user);
297:   /* Volumetric Mesh Projection */
298:   if (!user.multifield) {
299:     TestFunctionProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user);
300:     TestFieldProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user);
301:   } else {
302:     DM dmOut;

304:     DMClone(dm, &dmOut);
305:     SetupOutputDiscretization(dmOut, dim, simplex, &user);
306:     TestFieldProjectionMultiple(dmOut, dm, NULL, NULL, NULL, "Volumetric Primary", &user);
307:     DMDestroy(&dmOut);
308:   }
309:   if (user.auxfield) {
310:     /* Volumetric Mesh Projection with Volumetric Data */
311:     CreateAuxiliaryVec(dm, &auxdm, &la, &user);
312:     TestFunctionProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user);
313:     TestFieldProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user);
314:     VecDestroy(&la);
315:     /* Update of Volumetric Auxiliary Data with primary Volumetric Data */
316:     DMGetLocalVector(dm, &la);
317:     VecSet(la, 1.0);
318:     TestFieldProjection(auxdm, dm, NULL, la, "Volumetric Auxiliary Update with Volumetric Primary", &user);
319:     DMRestoreLocalVector(dm, &la);
320:     DMDestroy(&auxdm);
321:   }
322:   if (user.subdomain) {
323:     DMLabel domLabel;

325:     /* Subdomain Mesh Projection */
326:     CreateSubdomainMesh(dm, &domLabel, &subdm, &user);
327:     TestFunctionProjection(subdm, NULL, NULL, NULL, "Subdomain Primary", &user);
328:     TestFieldProjection(subdm, NULL, NULL, NULL, "Subdomain Primary", &user);
329:     if (user.auxfield) {
330:       /* Subdomain Mesh Projection with Subdomain Data */
331:       CreateAuxiliaryVec(subdm, &auxdm, &la, &user);
332:       TestFunctionProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Subdomain Auxiliary", &user);
333:       TestFieldProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Subdomain Auxiliary", &user);
334:       VecDestroy(&la);
335:       DMDestroy(&auxdm);
336:       /* Subdomain Mesh Projection with Volumetric Data */
337:       CreateAuxiliaryVec(dm, &auxdm, &la, &user);
338:       TestFunctionProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Volumetric Auxiliary", &user);
339:       TestFieldProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Volumetric Auxiliary", &user);
340:       VecDestroy(&la);
341:       DMDestroy(&auxdm);
342:       /* Volumetric Mesh Projection with Subdomain Data */
343:       CreateAuxiliaryVec(subdm, &auxdm, &la, &user);
344:       TestFunctionProjection(subdm, auxdm, domLabel, la, "Volumetric Primary and Subdomain Auxiliary", &user);
345:       TestFieldProjection(subdm, auxdm, domLabel, la, "Volumetric Primary and Subdomain Auxiliary", &user);
346:       VecDestroy(&la);
347:       DMDestroy(&auxdm);
348:     }
349:     DMDestroy(&subdm);
350:     DMLabelDestroy(&domLabel);
351:   }
352:   if (user.submesh) {
353:     DMLabel bdLabel;

355:     /* Boundary Mesh Projection */
356:     CreateBoundaryMesh(dm, &bdLabel, &subdm, &user);
357:     TestFunctionProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user);
358:     TestFieldProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user);
359:     if (user.auxfield) {
360:       /* Boundary Mesh Projection with Boundary Data */
361:       CreateAuxiliaryVec(subdm, &auxdm, &la, &user);
362:       TestFunctionProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user);
363:       TestFieldProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user);
364:       VecDestroy(&la);
365:       DMDestroy(&auxdm);
366:       /* Volumetric Mesh Projection with Boundary Data */
367:       CreateAuxiliaryVec(subdm, &auxdm, &la, &user);
368:       TestFunctionProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user);
369:       TestFieldProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user);
370:       VecDestroy(&la);
371:       DMDestroy(&auxdm);
372:     }
373:     DMLabelDestroy(&bdLabel);
374:     DMDestroy(&subdm);
375:   }
376:   DMDestroy(&dm);
377:   PetscFinalize();
378:   return 0;
379: }

381: /*TEST

383:   test:
384:     suffix: 0
385:     requires: triangle
386:     args: -dm_plex_box_faces 1,1 -func_view -local_func_view -local_input_view -local_field_view
387:   test:
388:     suffix: mf_0
389:     requires: triangle
390:     args: -dm_plex_box_faces 1,1 -velocity_petscspace_degree 1 -velocity_petscfe_default_quadrature_order 2 \
391:          -pressure_petscspace_degree 2 -pressure_petscfe_default_quadrature_order 2 \
392:          -multifield -output_petscspace_degree 1 -output_petscfe_default_quadrature_order 2 \
393:          -local_input_view -local_field_view
394:   test:
395:     suffix: 1
396:     requires: triangle
397:     args: -dm_plex_box_faces 1,1 -velocity_petscspace_degree 1 -velocity_petscfe_default_quadrature_order 2 -pressure_petscspace_degree 2 -pressure_petscfe_default_quadrature_order 2 -func_view -local_func_view -local_input_view -local_field_view -submesh -auxfield
398:   test:
399:     suffix: 2
400:     requires: triangle
401:     args: -dm_plex_box_faces 1,1 -velocity_petscspace_degree 1 -velocity_petscfe_default_quadrature_order 2 -pressure_petscspace_degree 2 -pressure_petscfe_default_quadrature_order 2 -func_view -local_func_view -local_input_view -local_field_view -subdomain -auxfield

403: TEST*/

405: /*
406:   Post-processing wants to project a function of the fields into some FE space
407:   - This is DMProjectField()
408:   - What about changing the number of components of the output, like displacement to stress? Aux vars

410:   Update of state variables
411:   - This is DMProjectField(), but solution must be the aux var
412: */