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