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