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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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:   PetscFunctionBegin;
 54:   options->multifield = PETSC_FALSE;
 55:   options->subdomain  = PETSC_FALSE;
 56:   options->submesh    = PETSC_FALSE;
 57:   options->auxfield   = PETSC_FALSE;

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

 68: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 69: {
 70:   PetscFunctionBegin;
 71:   PetscCall(DMCreate(comm, dm));
 72:   PetscCall(DMSetType(*dm, DMPLEX));
 73:   PetscCall(DMSetFromOptions(*dm));
 74:   PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

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

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

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

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

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

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

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

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

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

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

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

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

209: static PetscErrorCode TestFieldProjection(DM dm, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
210: {
211:   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
212:   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[]);
213:   Vec      lx, lu;
214:   PetscInt Nf, f;
215:   PetscInt val[1] = {1};
216:   char     lname[PETSC_MAX_PATH_LEN];

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

246: static PetscErrorCode TestFieldProjectionMultiple(DM dm, DM dmIn, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
247: {
248:   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
249:   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[]);
250:   Vec      lx, lu;
251:   PetscInt Nf, NfIn;
252:   PetscInt val[1] = {1};
253:   char     lname[PETSC_MAX_PATH_LEN];

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

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

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

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

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

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

383: /*TEST

385:   test:
386:     suffix: 0
387:     requires: triangle
388:     args: -dm_plex_box_faces 1,1 -func_view -local_func_view -local_input_view -local_field_view
389:   test:
390:     suffix: mf_0
391:     requires: triangle
392:     args: -dm_plex_box_faces 1,1 -velocity_petscspace_degree 1 -velocity_petscfe_default_quadrature_order 2 \
393:          -pressure_petscspace_degree 2 -pressure_petscfe_default_quadrature_order 2 \
394:          -multifield -output_petscspace_degree 1 -output_petscfe_default_quadrature_order 2 \
395:          -local_input_view -local_field_view
396:   test:
397:     suffix: 1
398:     requires: triangle
399:     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
400:   test:
401:     suffix: 2
402:     requires: triangle
403:     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

405: TEST*/

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

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