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