xref: /petsc/src/dm/impls/plex/tests/ex23.c (revision 44c45e3bf38ab8eca9b0c5f7f428c800d7bb968e) 
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 CreateAuxiliaryVec(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) {ierr = DMSetAuxiliaryVec(dm, NULL, 0, la);CHKERRQ(ierr);}
215   ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
216   ierr = PetscMalloc1(Nf, &funcs);CHKERRQ(ierr);
217   for (f = 0; f < Nf; ++f) funcs[f] = linear;
218   ierr = DMGetGlobalVector(dm, &x);CHKERRQ(ierr);
219   ierr = PetscStrcpy(lname, "Function ");CHKERRQ(ierr);
220   ierr = PetscStrcat(lname, name);CHKERRQ(ierr);
221   ierr = PetscObjectSetName((PetscObject) x, lname);CHKERRQ(ierr);
222   if (!label) {ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, x);CHKERRQ(ierr);}
223   else        {ierr = DMProjectFunctionLabel(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, x);CHKERRQ(ierr);}
224   ierr = VecViewFromOptions(x, NULL, "-func_view");CHKERRQ(ierr);
225   ierr = DMRestoreGlobalVector(dm, &x);CHKERRQ(ierr);
226   ierr = DMGetLocalVector(dm, &lx);CHKERRQ(ierr);
227   ierr = PetscStrcpy(lname, "Local Function ");CHKERRQ(ierr);
228   ierr = PetscStrcat(lname, name);CHKERRQ(ierr);
229   ierr = PetscObjectSetName((PetscObject) lx, lname);CHKERRQ(ierr);
230   if (!label) {ierr = DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_VALUES, lx);CHKERRQ(ierr);}
231   else        {ierr = DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, lx);CHKERRQ(ierr);}
232   ierr = VecViewFromOptions(lx, NULL, "-local_func_view");CHKERRQ(ierr);
233   ierr = DMRestoreLocalVector(dm, &lx);CHKERRQ(ierr);
234   ierr = PetscFree(funcs);CHKERRQ(ierr);
235   if (dmAux) {ierr = DMSetAuxiliaryVec(dm, NULL, 0, NULL);CHKERRQ(ierr);}
236   PetscFunctionReturn(0);
237 }
238 
239 static PetscErrorCode TestFieldProjection(DM dm, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
240 {
241   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
242   void           (**funcs)(PetscInt, PetscInt, PetscInt,
243                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
244                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
245                            PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
246   Vec               lx, lu;
247   PetscInt          Nf, f;
248   PetscInt          val[1] = {1};
249   char              lname[PETSC_MAX_PATH_LEN];
250   PetscErrorCode    ierr;
251 
252   PetscFunctionBeginUser;
253   if (dmAux) {ierr = DMSetAuxiliaryVec(dm, NULL, 0, la);CHKERRQ(ierr);}
254   ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
255   ierr = PetscMalloc2(Nf, &funcs, Nf, &afuncs);CHKERRQ(ierr);
256   for (f = 0; f < Nf; ++f) afuncs[f]  = linear;
257   funcs[0] = linear_vector;
258   funcs[1] = linear_scalar;
259   ierr = DMGetLocalVector(dm, &lu);CHKERRQ(ierr);
260   ierr = PetscStrcpy(lname, "Local Field Input ");CHKERRQ(ierr);
261   ierr = PetscStrcat(lname, name);CHKERRQ(ierr);
262   ierr = PetscObjectSetName((PetscObject) lu, lname);CHKERRQ(ierr);
263   if (!label) {ierr = DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, lu);CHKERRQ(ierr);}
264   else        {ierr = DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu);CHKERRQ(ierr);}
265   ierr = VecViewFromOptions(lu, NULL, "-local_input_view");CHKERRQ(ierr);
266   ierr = DMGetLocalVector(dm, &lx);CHKERRQ(ierr);
267   ierr = PetscStrcpy(lname, "Local Field ");CHKERRQ(ierr);
268   ierr = PetscStrcat(lname, name);CHKERRQ(ierr);
269   ierr = PetscObjectSetName((PetscObject) lx, lname);CHKERRQ(ierr);
270   if (!label) {ierr = DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx);CHKERRQ(ierr);}
271   else        {ierr = DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx);CHKERRQ(ierr);}
272   ierr = VecViewFromOptions(lx, NULL, "-local_field_view");CHKERRQ(ierr);
273   ierr = DMRestoreLocalVector(dm, &lx);CHKERRQ(ierr);
274   ierr = DMRestoreLocalVector(dm, &lu);CHKERRQ(ierr);
275   ierr = PetscFree2(funcs, afuncs);CHKERRQ(ierr);
276   if (dmAux) {ierr = DMSetAuxiliaryVec(dm, NULL, 0, NULL);CHKERRQ(ierr);}
277   PetscFunctionReturn(0);
278 }
279 
280 static PetscErrorCode TestFieldProjectionMultiple(DM dm, DM dmIn, DM dmAux, DMLabel label, Vec la, const char name[], AppCtx *user)
281 {
282   PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
283   void           (**funcs)(PetscInt, PetscInt, PetscInt,
284                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
285                            const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
286                            PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
287   Vec               lx, lu;
288   PetscInt          Nf, NfIn;
289   PetscInt          val[1] = {1};
290   char              lname[PETSC_MAX_PATH_LEN];
291   PetscErrorCode    ierr;
292 
293   PetscFunctionBeginUser;
294   if (dmAux) {ierr = DMSetAuxiliaryVec(dm, NULL, 0, la);CHKERRQ(ierr);}
295   ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
296   ierr = DMGetNumFields(dmIn, &NfIn);CHKERRQ(ierr);
297   ierr = PetscMalloc2(Nf, &funcs, NfIn, &afuncs);CHKERRQ(ierr);
298   funcs[0]  = divergence_sq;
299   afuncs[0] = linear2;
300   afuncs[1] = linear;
301   ierr = DMGetLocalVector(dmIn, &lu);CHKERRQ(ierr);
302   ierr = PetscStrcpy(lname, "Local MultiField Input ");CHKERRQ(ierr);
303   ierr = PetscStrcat(lname, name);CHKERRQ(ierr);
304   ierr = PetscObjectSetName((PetscObject) lu, lname);CHKERRQ(ierr);
305   if (!label) {ierr = DMProjectFunctionLocal(dmIn, 0.0, afuncs, NULL, INSERT_VALUES, lu);CHKERRQ(ierr);}
306   else        {ierr = DMProjectFunctionLabelLocal(dmIn, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu);CHKERRQ(ierr);}
307   ierr = VecViewFromOptions(lu, NULL, "-local_input_view");CHKERRQ(ierr);
308   ierr = DMGetLocalVector(dm, &lx);CHKERRQ(ierr);
309   ierr = PetscStrcpy(lname, "Local MultiField ");CHKERRQ(ierr);
310   ierr = PetscStrcat(lname, name);CHKERRQ(ierr);
311   ierr = PetscObjectSetName((PetscObject) lx, lname);CHKERRQ(ierr);
312   if (!label) {ierr = DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx);CHKERRQ(ierr);}
313   else        {ierr = DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx);CHKERRQ(ierr);}
314   ierr = VecViewFromOptions(lx, NULL, "-local_field_view");CHKERRQ(ierr);
315   ierr = DMRestoreLocalVector(dm, &lx);CHKERRQ(ierr);
316   ierr = DMRestoreLocalVector(dmIn, &lu);CHKERRQ(ierr);
317   ierr = PetscFree2(funcs, afuncs);CHKERRQ(ierr);
318   if (dmAux) {ierr = DMSetAuxiliaryVec(dm, NULL, 0, NULL);CHKERRQ(ierr);}
319   PetscFunctionReturn(0);
320 }
321 
322 int main(int argc, char **argv)
323 {
324   DM             dm, subdm, auxdm;
325   Vec            la;
326   AppCtx         user;
327   PetscErrorCode ierr;
328 
329   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
330   ierr = ProcessOptions(&user);CHKERRQ(ierr);
331   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
332   ierr = SetupDiscretization(dm, user.dim, user.cellSimplex, &user);CHKERRQ(ierr);
333   /* Volumetric Mesh Projection */
334   if (!user.multifield) {
335     ierr = TestFunctionProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user);CHKERRQ(ierr);
336     ierr = TestFieldProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user);CHKERRQ(ierr);
337   } else {
338     DM dmOut;
339 
340     ierr = DMClone(dm, &dmOut);CHKERRQ(ierr);
341     ierr = SetupOutputDiscretization(dmOut, user.dim, user.cellSimplex, &user);CHKERRQ(ierr);
342     ierr = TestFieldProjectionMultiple(dmOut, dm, NULL, NULL, NULL, "Volumetric Primary", &user);CHKERRQ(ierr);
343     ierr = DMDestroy(&dmOut);CHKERRQ(ierr);
344   }
345   if (user.auxfield) {
346     /* Volumetric Mesh Projection with Volumetric Data */
347     ierr = CreateAuxiliaryVec(dm, &auxdm, &la, &user);CHKERRQ(ierr);
348     ierr = TestFunctionProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user);CHKERRQ(ierr);
349     ierr = TestFieldProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user);CHKERRQ(ierr);
350     ierr = VecDestroy(&la);CHKERRQ(ierr);
351     /* Update of Volumetric Auxiliary Data with primary Volumetric Data */
352     ierr = DMGetLocalVector(dm, &la);CHKERRQ(ierr);
353     ierr = VecSet(la, 1.0);CHKERRQ(ierr);
354     ierr = TestFieldProjection(auxdm, dm, NULL, la, "Volumetric Auxiliary Update with Volumetric Primary", &user);CHKERRQ(ierr);
355     ierr = DMRestoreLocalVector(dm, &la);CHKERRQ(ierr);
356     ierr = DMDestroy(&auxdm);CHKERRQ(ierr);
357   }
358   if (user.subdomain) {
359     DMLabel domLabel;
360 
361     /* Subdomain Mesh Projection */
362     ierr = CreateSubdomainMesh(dm, &domLabel, &subdm, &user);CHKERRQ(ierr);
363     ierr = TestFunctionProjection(subdm, NULL, NULL, NULL, "Subdomain Primary", &user);CHKERRQ(ierr);
364     ierr = TestFieldProjection(subdm, NULL, NULL, NULL, "Subdomain Primary", &user);CHKERRQ(ierr);
365     if (user.auxfield) {
366       /* Subdomain Mesh Projection with Subdomain Data */
367       ierr = CreateAuxiliaryVec(subdm, &auxdm, &la, &user);CHKERRQ(ierr);
368       ierr = TestFunctionProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Subdomain Auxiliary", &user);CHKERRQ(ierr);
369       ierr = TestFieldProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Subdomain Auxiliary", &user);CHKERRQ(ierr);
370       ierr = VecDestroy(&la);CHKERRQ(ierr);
371       ierr = DMDestroy(&auxdm);CHKERRQ(ierr);
372       /* Subdomain Mesh Projection with Volumetric Data */
373       ierr = CreateAuxiliaryVec(dm, &auxdm, &la, &user);CHKERRQ(ierr);
374       ierr = TestFunctionProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Volumetric Auxiliary", &user);CHKERRQ(ierr);
375       ierr = TestFieldProjection(subdm, auxdm, NULL, la, "Subdomain Primary and Volumetric Auxiliary", &user);CHKERRQ(ierr);
376       ierr = VecDestroy(&la);CHKERRQ(ierr);
377       ierr = DMDestroy(&auxdm);CHKERRQ(ierr);
378       /* Volumetric Mesh Projection with Subdomain Data */
379       ierr = CreateAuxiliaryVec(subdm, &auxdm, &la, &user);CHKERRQ(ierr);
380       ierr = TestFunctionProjection(subdm, auxdm, domLabel, la, "Volumetric Primary and Subdomain Auxiliary", &user);CHKERRQ(ierr);
381       ierr = TestFieldProjection(subdm, auxdm, domLabel, la, "Volumetric Primary and Subdomain Auxiliary", &user);CHKERRQ(ierr);
382       ierr = VecDestroy(&la);CHKERRQ(ierr);
383       ierr = DMDestroy(&auxdm);CHKERRQ(ierr);
384     }
385     ierr = DMDestroy(&subdm);CHKERRQ(ierr);
386     ierr = DMLabelDestroy(&domLabel);CHKERRQ(ierr);
387   }
388   if (user.submesh) {
389     DMLabel bdLabel;
390 
391     /* Boundary Mesh Projection */
392     ierr = CreateBoundaryMesh(dm, &bdLabel, &subdm, &user);CHKERRQ(ierr);
393     ierr = TestFunctionProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user);CHKERRQ(ierr);
394     ierr = TestFieldProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user);CHKERRQ(ierr);
395     if (user.auxfield) {
396       /* Boundary Mesh Projection with Boundary Data */
397       ierr = CreateAuxiliaryVec(subdm, &auxdm, &la, &user);CHKERRQ(ierr);
398       ierr = TestFunctionProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);
399       ierr = TestFieldProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);
400       ierr = VecDestroy(&la);CHKERRQ(ierr);
401       ierr = DMDestroy(&auxdm);CHKERRQ(ierr);
402       /* Volumetric Mesh Projection with Boundary Data */
403       ierr = CreateAuxiliaryVec(subdm, &auxdm, &la, &user);CHKERRQ(ierr);
404       ierr = TestFunctionProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);
405       ierr = TestFieldProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);
406       ierr = VecDestroy(&la);CHKERRQ(ierr);
407       ierr = DMDestroy(&auxdm);CHKERRQ(ierr);
408     }
409     ierr = DMLabelDestroy(&bdLabel);CHKERRQ(ierr);
410     ierr = DMDestroy(&subdm);CHKERRQ(ierr);
411   }
412   ierr = DMDestroy(&dm);CHKERRQ(ierr);
413   ierr = PetscFinalize();
414   return ierr;
415 }
416 
417 /*TEST
418 
419   test:
420     suffix: 0
421     requires: triangle
422     args: -dim 2 -func_view -local_func_view -local_input_view -local_field_view
423   test:
424     suffix: mf_0
425     requires: triangle
426     args: -dim 2 -velocity_petscspace_degree 1 -velocity_petscfe_default_quadrature_order 2 \
427          -pressure_petscspace_degree 2 -pressure_petscfe_default_quadrature_order 2 \
428          -multifield -output_petscspace_degree 1 -output_petscfe_default_quadrature_order 2 \
429          -local_input_view -local_field_view
430   test:
431     suffix: 1
432     requires: triangle
433     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
434   test:
435     suffix: 2
436     requires: triangle
437     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
438 
439 TEST*/
440 
441 /*
442   Post-processing wants to project a function of the fields into some FE space
443   - This is DMProjectField()
444   - What about changing the number of components of the output, like displacement to stress? Aux vars
445 
446   Update of state variables
447   - This is DMProjectField(), but solution must be the aux var
448 */
449