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