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, PetscCtx 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, PetscCtx 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:   for (PetscInt d = 0; d < uOff[1] - uOff[0]; ++d) f[d] = u[d + uOff[0]];
 33: }

 35: /* p */
 36: 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[])
 37: {
 38:   f[0] = u[uOff[1]];
 39: }

 41: /* {div u, p^2} */
 42: 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[])
 43: {
 44:   f[0] = 0.0;
 45:   for (PetscInt d = 0; d < dim; ++d) f[0] += u_x[uOff_x[0] + d * dim + d];
 46:   f[1] = PetscSqr(u[uOff[1]]);
 47: }

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

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

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

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

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

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

100:   PetscFunctionBeginUser;
101:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
102:   PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "output_", -1, &fe));
103:   PetscCall(PetscObjectSetName((PetscObject)fe, "output"));
104:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
105:   PetscCall(PetscFEDestroy(&fe));
106:   PetscCall(DMCreateDS(dm));
107:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

158:   PetscFunctionBeginUser;
159:   PetscCall(DMGetDimension(dm, &dim));
160:   PetscCall(DMPlexIsSimplex(dm, &simplex));
161:   PetscCall(DMGetNumFields(dm, &Nf));
162:   PetscCall(PetscMalloc1(Nf, &afuncs));
163:   for (f = 0; f < Nf; ++f) afuncs[f] = linear;
164:   PetscCall(DMClone(dm, auxdm));
165:   PetscCall(SetupDiscretization(*auxdm, dim, simplex, user));
166:   PetscCall(DMCreateLocalVector(*auxdm, la));
167:   PetscCall(DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, *la));
168:   PetscCall(VecViewFromOptions(*la, NULL, "-local_aux_view"));
169:   PetscCall(PetscFree(afuncs));
170:   PetscFunctionReturn(PETSC_SUCCESS);
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;
178:   PetscInt val[1] = {1};
179:   char     lname[PETSC_MAX_PATH_LEN];

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

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

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

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

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

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

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