xref: /petsc/src/dm/impls/plex/tests/ex23.c (revision f97672e55eacc8688507b9471cd7ec2664d7f203)
1 static char help[] = "Test for function and field projection\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petscds.h>
5 
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;
12 
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 }
20 
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 }
28 
29 /* {u_x, u_y, u_z} */
30 static void linear_vector(PetscInt dim, PetscInt Nf, PetscInt NfAux,
31                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
32                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
33                           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
34 {
35   PetscInt d;
36   for (d = 0; d < uOff[1]-uOff[0]; ++d) f[d] = u[d+uOff[0]];
37 }
38 
39 /* p */
40 static void linear_scalar(PetscInt dim, PetscInt Nf, PetscInt NfAux,
41                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
42                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
43                           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
44 {
45   f[0] = u[uOff[1]];
46 }
47 
48 /* {div u, p^2} */
49 static void divergence_sq(PetscInt dim, PetscInt Nf, PetscInt NfAux,
50                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
51                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
52                           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
53 {
54   PetscInt d;
55   f[0] = 0.0;
56   for (d = 0; d < dim; ++d) f[0] += u_x[uOff_x[0]+d*dim+d];
57   f[1] = PetscSqr(u[uOff[1]]);
58 }
59 
60 static PetscErrorCode ProcessOptions(AppCtx *options)
61 {
62   PetscFunctionBegin;
63   options->multifield  = PETSC_FALSE;
64   options->subdomain   = PETSC_FALSE;
65   options->submesh     = PETSC_FALSE;
66   options->auxfield    = PETSC_FALSE;
67 
68   PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
69   PetscCall(PetscOptionsBool("-multifield", "Flag for trying different numbers of input/output fields", "ex23.c", options->multifield, &options->multifield, NULL));
70   PetscCall(PetscOptionsBool("-subdomain", "Flag for trying volumetric submesh", "ex23.c", options->subdomain, &options->subdomain, NULL));
71   PetscCall(PetscOptionsBool("-submesh", "Flag for trying boundary submesh", "ex23.c", options->submesh, &options->submesh, NULL));
72   PetscCall(PetscOptionsBool("-auxfield", "Flag for trying auxiliary fields", "ex23.c", options->auxfield, &options->auxfield, NULL));
73   PetscOptionsEnd();
74   PetscFunctionReturn(0);
75 }
76 
77 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
78 {
79   PetscFunctionBegin;
80   PetscCall(DMCreate(comm, dm));
81   PetscCall(DMSetType(*dm, DMPLEX));
82   PetscCall(DMSetFromOptions(*dm));
83   PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
84   PetscFunctionReturn(0);
85 }
86 
87 static PetscErrorCode SetupDiscretization(DM dm, PetscInt dim, PetscBool simplex, AppCtx *user)
88 {
89   PetscFE        fe;
90   MPI_Comm       comm;
91 
92   PetscFunctionBeginUser;
93   PetscCall(PetscObjectGetComm((PetscObject) dm, &comm));
94   PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "velocity_", -1, &fe));
95   PetscCall(PetscObjectSetName((PetscObject) fe, "velocity"));
96   PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe));
97   PetscCall(PetscFEDestroy(&fe));
98   PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pressure_", -1, &fe));
99   PetscCall(PetscObjectSetName((PetscObject) fe, "pressure"));
100   PetscCall(DMSetField(dm, 1, NULL, (PetscObject) fe));
101   PetscCall(PetscFEDestroy(&fe));
102   PetscCall(DMCreateDS(dm));
103   PetscFunctionReturn(0);
104 }
105 
106 static PetscErrorCode SetupOutputDiscretization(DM dm, PetscInt dim, PetscBool simplex, AppCtx *user)
107 {
108   PetscFE        fe;
109   MPI_Comm       comm;
110 
111   PetscFunctionBeginUser;
112   PetscCall(PetscObjectGetComm((PetscObject) dm, &comm));
113   PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "output_", -1, &fe));
114   PetscCall(PetscObjectSetName((PetscObject) fe, "output"));
115   PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe));
116   PetscCall(PetscFEDestroy(&fe));
117   PetscCall(DMCreateDS(dm));
118   PetscFunctionReturn(0);
119 }
120 
121 static PetscErrorCode CreateSubdomainMesh(DM dm, DMLabel *domLabel, DM *subdm, AppCtx *user)
122 {
123   DMLabel        label;
124   PetscBool      simplex;
125   PetscInt       dim, cStart, cEnd, c;
126 
127   PetscFunctionBeginUser;
128   PetscCall(DMPlexIsSimplex(dm, &simplex));
129   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
130   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subdomain", &label));
131   for (c = cStart + (cEnd-cStart)/2; c < cEnd; ++c) PetscCall(DMLabelSetValue(label, c, 1));
132   PetscCall(DMPlexFilter(dm, label, 1, subdm));
133   PetscCall(DMGetDimension(*subdm, &dim));
134   PetscCall(SetupDiscretization(*subdm, dim, simplex, user));
135   PetscCall(PetscObjectSetName((PetscObject) *subdm, "subdomain"));
136   PetscCall(DMViewFromOptions(*subdm, NULL, "-sub_dm_view"));
137   if (domLabel) *domLabel = label;
138   else          PetscCall(DMLabelDestroy(&label));
139   PetscFunctionReturn(0);
140 }
141 
142 static PetscErrorCode CreateBoundaryMesh(DM dm, DMLabel *bdLabel, DM *subdm, AppCtx *user)
143 {
144   DMLabel        label;
145   PetscBool      simplex;
146   PetscInt       dim;
147 
148   PetscFunctionBeginUser;
149   PetscCall(DMPlexIsSimplex(dm, &simplex));
150   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "sub", &label));
151   PetscCall(DMPlexMarkBoundaryFaces(dm, 1, label));
152   PetscCall(DMPlexLabelComplete(dm, label));
153   PetscCall(DMPlexCreateSubmesh(dm, label, 1, PETSC_TRUE, subdm));
154   PetscCall(DMGetDimension(*subdm, &dim));
155   PetscCall(SetupDiscretization(*subdm, dim, simplex, user));
156   PetscCall(PetscObjectSetName((PetscObject) *subdm, "boundary"));
157   PetscCall(DMViewFromOptions(*subdm, NULL, "-sub_dm_view"));
158   if (bdLabel) *bdLabel = label;
159   else         PetscCall(DMLabelDestroy(&label));
160   PetscFunctionReturn(0);
161 }
162 
163 static PetscErrorCode CreateAuxiliaryVec(DM dm, DM *auxdm, Vec *la, AppCtx *user)
164 {
165   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
166   PetscBool         simplex;
167   PetscInt          dim, Nf, f;
168 
169   PetscFunctionBeginUser;
170   PetscCall(DMGetDimension(dm, &dim));
171   PetscCall(DMPlexIsSimplex(dm, &simplex));
172   PetscCall(DMGetNumFields(dm, &Nf));
173   PetscCall(PetscMalloc1(Nf, &afuncs));
174   for (f = 0; f < Nf; ++f) afuncs[f]  = linear;
175   PetscCall(DMClone(dm, auxdm));
176   PetscCall(SetupDiscretization(*auxdm, dim, simplex, user));
177   PetscCall(DMCreateLocalVector(*auxdm, la));
178   PetscCall(DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, *la));
179   PetscCall(VecViewFromOptions(*la, NULL, "-local_aux_view"));
180   PetscCall(PetscFree(afuncs));
181   PetscFunctionReturn(0);
182 }
183 
184 static PetscErrorCode TestFunctionProjection(DM dm, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
185 {
186   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
187   Vec               x, lx;
188   PetscInt          Nf, f;
189   PetscInt          val[1] = {1};
190   char              lname[PETSC_MAX_PATH_LEN];
191 
192   PetscFunctionBeginUser;
193   if (dmAux) PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, la));
194   PetscCall(DMGetNumFields(dm, &Nf));
195   PetscCall(PetscMalloc1(Nf, &funcs));
196   for (f = 0; f < Nf; ++f) funcs[f] = linear;
197   PetscCall(DMGetGlobalVector(dm, &x));
198   PetscCall(PetscStrcpy(lname, "Function "));
199   PetscCall(PetscStrcat(lname, name));
200   PetscCall(PetscObjectSetName((PetscObject) x, lname));
201   if (!label) PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, x));
202   else        PetscCall(DMProjectFunctionLabel(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, x));
203   PetscCall(VecViewFromOptions(x, NULL, "-func_view"));
204   PetscCall(DMRestoreGlobalVector(dm, &x));
205   PetscCall(DMGetLocalVector(dm, &lx));
206   PetscCall(PetscStrcpy(lname, "Local Function "));
207   PetscCall(PetscStrcat(lname, name));
208   PetscCall(PetscObjectSetName((PetscObject) lx, lname));
209   if (!label) PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_VALUES, lx));
210   else        PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, lx));
211   PetscCall(VecViewFromOptions(lx, NULL, "-local_func_view"));
212   PetscCall(DMRestoreLocalVector(dm, &lx));
213   PetscCall(PetscFree(funcs));
214   if (dmAux) PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL));
215   PetscFunctionReturn(0);
216 }
217 
218 static PetscErrorCode TestFieldProjection(DM dm, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
219 {
220   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
221   void           (**funcs)(PetscInt, PetscInt, PetscInt,
222                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
223                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
224                            PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
225   Vec               lx, lu;
226   PetscInt          Nf, f;
227   PetscInt          val[1] = {1};
228   char              lname[PETSC_MAX_PATH_LEN];
229 
230   PetscFunctionBeginUser;
231   if (dmAux) PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, la));
232   PetscCall(DMGetNumFields(dm, &Nf));
233   PetscCall(PetscMalloc2(Nf, &funcs, Nf, &afuncs));
234   for (f = 0; f < Nf; ++f) afuncs[f]  = linear;
235   funcs[0] = linear_vector;
236   funcs[1] = linear_scalar;
237   PetscCall(DMGetLocalVector(dm, &lu));
238   PetscCall(PetscStrcpy(lname, "Local Field Input "));
239   PetscCall(PetscStrcat(lname, name));
240   PetscCall(PetscObjectSetName((PetscObject) lu, lname));
241   if (!label) PetscCall(DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, lu));
242   else        PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu));
243   PetscCall(VecViewFromOptions(lu, NULL, "-local_input_view"));
244   PetscCall(DMGetLocalVector(dm, &lx));
245   PetscCall(PetscStrcpy(lname, "Local Field "));
246   PetscCall(PetscStrcat(lname, name));
247   PetscCall(PetscObjectSetName((PetscObject) lx, lname));
248   if (!label) PetscCall(DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx));
249   else        PetscCall(DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx));
250   PetscCall(VecViewFromOptions(lx, NULL, "-local_field_view"));
251   PetscCall(DMRestoreLocalVector(dm, &lx));
252   PetscCall(DMRestoreLocalVector(dm, &lu));
253   PetscCall(PetscFree2(funcs, afuncs));
254   if (dmAux) PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL));
255   PetscFunctionReturn(0);
256 }
257 
258 static PetscErrorCode TestFieldProjectionMultiple(DM dm, DM dmIn, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
259 {
260   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
261   void           (**funcs)(PetscInt, PetscInt, PetscInt,
262                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
263                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
264                            PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
265   Vec               lx, lu;
266   PetscInt          Nf, NfIn;
267   PetscInt          val[1] = {1};
268   char              lname[PETSC_MAX_PATH_LEN];
269 
270   PetscFunctionBeginUser;
271   if (dmAux) PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, la));
272   PetscCall(DMGetNumFields(dm, &Nf));
273   PetscCall(DMGetNumFields(dmIn, &NfIn));
274   PetscCall(PetscMalloc2(Nf, &funcs, NfIn, &afuncs));
275   funcs[0]  = divergence_sq;
276   afuncs[0] = linear2;
277   afuncs[1] = linear;
278   PetscCall(DMGetLocalVector(dmIn, &lu));
279   PetscCall(PetscStrcpy(lname, "Local MultiField Input "));
280   PetscCall(PetscStrcat(lname, name));
281   PetscCall(PetscObjectSetName((PetscObject) lu, lname));
282   if (!label) PetscCall(DMProjectFunctionLocal(dmIn, 0.0, afuncs, NULL, INSERT_VALUES, lu));
283   else        PetscCall(DMProjectFunctionLabelLocal(dmIn, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu));
284   PetscCall(VecViewFromOptions(lu, NULL, "-local_input_view"));
285   PetscCall(DMGetLocalVector(dm, &lx));
286   PetscCall(PetscStrcpy(lname, "Local MultiField "));
287   PetscCall(PetscStrcat(lname, name));
288   PetscCall(PetscObjectSetName((PetscObject) lx, lname));
289   if (!label) PetscCall(DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx));
290   else        PetscCall(DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx));
291   PetscCall(VecViewFromOptions(lx, NULL, "-local_field_view"));
292   PetscCall(DMRestoreLocalVector(dm, &lx));
293   PetscCall(DMRestoreLocalVector(dmIn, &lu));
294   PetscCall(PetscFree2(funcs, afuncs));
295   if (dmAux) PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL));
296   PetscFunctionReturn(0);
297 }
298 
299 int main(int argc, char **argv)
300 {
301   DM             dm, subdm, auxdm;
302   Vec            la;
303   PetscInt       dim;
304   PetscBool      simplex;
305   AppCtx         user;
306 
307   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
308   PetscCall(ProcessOptions(&user));
309   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
310   PetscCall(DMGetDimension(dm, &dim));
311   PetscCall(DMPlexIsSimplex(dm, &simplex));
312   PetscCall(SetupDiscretization(dm, dim, simplex, &user));
313   /* Volumetric Mesh Projection */
314   if (!user.multifield) {
315     PetscCall(TestFunctionProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user));
316     PetscCall(TestFieldProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user));
317   } else {
318     DM dmOut;
319 
320     PetscCall(DMClone(dm, &dmOut));
321     PetscCall(SetupOutputDiscretization(dmOut, dim, simplex, &user));
322     PetscCall(TestFieldProjectionMultiple(dmOut, dm, NULL, NULL, NULL, "Volumetric Primary", &user));
323     PetscCall(DMDestroy(&dmOut));
324   }
325   if (user.auxfield) {
326     /* Volumetric Mesh Projection with Volumetric Data */
327     PetscCall(CreateAuxiliaryVec(dm, &auxdm, &la, &user));
328     PetscCall(TestFunctionProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user));
329     PetscCall(TestFieldProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user));
330     PetscCall(VecDestroy(&la));
331     /* Update of Volumetric Auxiliary Data with primary Volumetric Data */
332     PetscCall(DMGetLocalVector(dm, &la));
333     PetscCall(VecSet(la, 1.0));
334     PetscCall(TestFieldProjection(auxdm, dm, NULL, la, "Volumetric Auxiliary Update with Volumetric Primary", &user));
335     PetscCall(DMRestoreLocalVector(dm, &la));
336     PetscCall(DMDestroy(&auxdm));
337   }
338   if (user.subdomain) {
339     DMLabel domLabel;
340 
341     /* Subdomain Mesh Projection */
342     PetscCall(CreateSubdomainMesh(dm, &domLabel, &subdm, &user));
343     PetscCall(TestFunctionProjection(subdm, NULL, NULL, NULL, "Subdomain Primary", &user));
344     PetscCall(TestFieldProjection(subdm, NULL, NULL, NULL, "Subdomain Primary", &user));
345     if (user.auxfield) {
346       /* Subdomain Mesh Projection with Subdomain Data */
347       PetscCall(CreateAuxiliaryVec(subdm, &auxdm, &la, &user));
348       PetscCall(TestFunctionProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Subdomain Auxiliary", &user));
349       PetscCall(TestFieldProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Subdomain Auxiliary", &user));
350       PetscCall(VecDestroy(&la));
351       PetscCall(DMDestroy(&auxdm));
352       /* Subdomain Mesh Projection with Volumetric Data */
353       PetscCall(CreateAuxiliaryVec(dm, &auxdm, &la, &user));
354       PetscCall(TestFunctionProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Volumetric Auxiliary", &user));
355       PetscCall(TestFieldProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Volumetric Auxiliary", &user));
356       PetscCall(VecDestroy(&la));
357       PetscCall(DMDestroy(&auxdm));
358       /* Volumetric Mesh Projection with Subdomain Data */
359       PetscCall(CreateAuxiliaryVec(subdm, &auxdm, &la, &user));
360       PetscCall(TestFunctionProjection(subdm, auxdm, domLabel, la, "Volumetric Primary and Subdomain Auxiliary", &user));
361       PetscCall(TestFieldProjection(subdm, auxdm, domLabel, la, "Volumetric Primary and Subdomain Auxiliary", &user));
362       PetscCall(VecDestroy(&la));
363       PetscCall(DMDestroy(&auxdm));
364     }
365     PetscCall(DMDestroy(&subdm));
366     PetscCall(DMLabelDestroy(&domLabel));
367   }
368   if (user.submesh) {
369     DMLabel bdLabel;
370 
371     /* Boundary Mesh Projection */
372     PetscCall(CreateBoundaryMesh(dm, &bdLabel, &subdm, &user));
373     PetscCall(TestFunctionProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user));
374     PetscCall(TestFieldProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user));
375     if (user.auxfield) {
376       /* Boundary Mesh Projection with Boundary Data */
377       PetscCall(CreateAuxiliaryVec(subdm, &auxdm, &la, &user));
378       PetscCall(TestFunctionProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user));
379       PetscCall(TestFieldProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user));
380       PetscCall(VecDestroy(&la));
381       PetscCall(DMDestroy(&auxdm));
382       /* Volumetric Mesh Projection with Boundary Data */
383       PetscCall(CreateAuxiliaryVec(subdm, &auxdm, &la, &user));
384       PetscCall(TestFunctionProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user));
385       PetscCall(TestFieldProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user));
386       PetscCall(VecDestroy(&la));
387       PetscCall(DMDestroy(&auxdm));
388     }
389     PetscCall(DMLabelDestroy(&bdLabel));
390     PetscCall(DMDestroy(&subdm));
391   }
392   PetscCall(DMDestroy(&dm));
393   PetscCall(PetscFinalize());
394   return 0;
395 }
396 
397 /*TEST
398 
399   test:
400     suffix: 0
401     requires: triangle
402     args: -dm_plex_box_faces 1,1 -func_view -local_func_view -local_input_view -local_field_view
403   test:
404     suffix: mf_0
405     requires: triangle
406     args: -dm_plex_box_faces 1,1 -velocity_petscspace_degree 1 -velocity_petscfe_default_quadrature_order 2 \
407          -pressure_petscspace_degree 2 -pressure_petscfe_default_quadrature_order 2 \
408          -multifield -output_petscspace_degree 1 -output_petscfe_default_quadrature_order 2 \
409          -local_input_view -local_field_view
410   test:
411     suffix: 1
412     requires: triangle
413     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
414   test:
415     suffix: 2
416     requires: triangle
417     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
418 
419 TEST*/
420 
421 /*
422   Post-processing wants to project a function of the fields into some FE space
423   - This is DMProjectField()
424   - What about changing the number of components of the output, like displacement to stress? Aux vars
425 
426   Update of state variables
427   - This is DMProjectField(), but solution must be the aux var
428 */
429