xref: /petsc/src/dm/dt/fe/tests/ex2.c (revision b45e3bf4ff73d80a20c3202b6cd9f79d2f2d3efe) 
1 static const char help[] = "Tests for injecting basis functions";
2 
3 #include <petscdmplex.h>
4 #include <petscfe.h>
5 #include <petscds.h>
6 
7 typedef struct {
8   PetscInt its; /* Number of replications for timing */
9 } AppCtx;
10 
11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12 {
13   PetscErrorCode ierr;
14 
15   PetscFunctionBeginUser;
16   options->its = 1;
17 
18   ierr = PetscOptionsBegin(comm, "", "FE Injection Options", "PETSCFE");CHKERRQ(ierr);
19   ierr = PetscOptionsInt("-its", "The number of replications for timing", "ex1.c", options->its, &options->its, NULL);CHKERRQ(ierr);
20   ierr = PetscOptionsEnd();CHKERRQ(ierr);
21   PetscFunctionReturn(0);
22 }
23 
24 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
25 {
26   PetscInt d;
27   *u = 0.0;
28   for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0*PETSC_PI*x[d]);
29   return 0;
30 }
31 
32 static void f0_trig_u(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 f0[])
36 {
37   PetscInt d;
38   for (d = 0; d < dim; ++d) f0[0] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]);
39 }
40 
41 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
42                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
43                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
44                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
45 {
46   PetscInt d;
47   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
48 }
49 
50 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
51                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
52                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
53                   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
54 {
55   PetscInt d;
56   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
57 }
58 
59 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
60 {
61   PetscDS        ds;
62   DMLabel        label;
63   const PetscInt id = 1;
64   PetscErrorCode ierr;
65 
66   PetscFunctionBeginUser;
67   ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
68   ierr = PetscDSSetResidual(ds, 0, f0_trig_u, f1_u);CHKERRQ(ierr);
69   ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
70   ierr = PetscDSSetExactSolution(ds, 0, trig_u, user);CHKERRQ(ierr);
71   ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr);
72   if (label) {ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) trig_u, NULL, user, NULL);CHKERRQ(ierr);}
73   PetscFunctionReturn(0);
74 }
75 
76 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
77 {
78   DM             cdm = dm;
79   PetscFE        fe;
80   char           prefix[PETSC_MAX_PATH_LEN];
81   PetscInt       dim;
82   PetscErrorCode ierr;
83 
84   PetscFunctionBeginUser;
85   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
86   ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr);
87   ierr = DMCreateFEDefault(dm, dim, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr);
88   ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr);
89   /* Set discretization and boundary conditions for each mesh */
90   ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
91   ierr = DMCreateDS(dm);CHKERRQ(ierr);
92   ierr = (*setup)(dm, user);CHKERRQ(ierr);
93   while (cdm) {
94     ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr);
95     ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
96   }
97   ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
98   PetscFunctionReturn(0);
99 }
100 
101 static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom(void *ctx)
102 {
103   PetscFEGeom   *geom = (PetscFEGeom *) ctx;
104   PetscErrorCode ierr;
105 
106   PetscFunctionBegin;
107   ierr = PetscFEGeomDestroy(&geom);CHKERRQ(ierr);
108   PetscFunctionReturn(0);
109 }
110 
111 PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
112 {
113   char            composeStr[33] = {0};
114   PetscObjectId   id;
115   PetscContainer  container;
116   PetscErrorCode  ierr;
117 
118   PetscFunctionBegin;
119   ierr = PetscObjectGetId((PetscObject) quad, &id);CHKERRQ(ierr);
120   ierr = PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%x\n", id);CHKERRQ(ierr);
121   ierr = PetscObjectQuery((PetscObject) cellIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr);
122   if (container) {
123     ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr);
124   } else {
125     ierr = DMFieldCreateFEGeom(coordField, cellIS, quad, faceData, geom);CHKERRQ(ierr);
126     ierr = PetscContainerCreate(PETSC_COMM_SELF, &container);CHKERRQ(ierr);
127     ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr);
128     ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr);
129     ierr = PetscObjectCompose((PetscObject) cellIS, composeStr, (PetscObject) container);CHKERRQ(ierr);
130     ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
131   }
132   PetscFunctionReturn(0);
133 }
134 
135 PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
136 {
137   PetscFunctionBegin;
138   *geom = NULL;
139   PetscFunctionReturn(0);
140 }
141 
142 static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
143 {
144   DMField        coordField;
145   PetscInt       Nf, f, maxDegree;
146   PetscErrorCode ierr;
147 
148   PetscFunctionBeginUser;
149   *affineQuad = NULL;
150   *affineGeom = NULL;
151   *quads      = NULL;
152   *geoms      = NULL;
153   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
154   ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr);
155   ierr = DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);CHKERRQ(ierr);
156   if (maxDegree <= 1) {
157     ierr = DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad);CHKERRQ(ierr);
158     if (*affineQuad) {ierr = CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom);CHKERRQ(ierr);}
159   } else {
160     ierr = PetscCalloc2(Nf, quads, Nf, geoms);CHKERRQ(ierr);
161     for (f = 0; f < Nf; ++f) {
162       PetscFE fe;
163 
164       ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr);
165       ierr = PetscFEGetQuadrature(fe, &(*quads)[f]);CHKERRQ(ierr);
166       ierr = PetscObjectReference((PetscObject) (*quads)[f]);CHKERRQ(ierr);
167       ierr = CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]);CHKERRQ(ierr);
168     }
169   }
170   PetscFunctionReturn(0);
171 }
172 
173 static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
174 {
175   DMField        coordField;
176   PetscInt       Nf, f;
177   PetscErrorCode ierr;
178 
179   PetscFunctionBeginUser;
180   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
181   ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr);
182   if (*affineQuad) {
183     ierr = CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom);CHKERRQ(ierr);
184     ierr = PetscQuadratureDestroy(affineQuad);CHKERRQ(ierr);
185   } else {
186     for (f = 0; f < Nf; ++f) {
187       ierr = CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]);CHKERRQ(ierr);
188       ierr = PetscQuadratureDestroy(&(*quads)[f]);CHKERRQ(ierr);
189     }
190     ierr = PetscFree2(*quads, *geoms);CHKERRQ(ierr);
191   }
192   PetscFunctionReturn(0);
193 }
194 
195 static PetscErrorCode TestEvaluation(DM dm)
196 {
197   PetscFE        fe;
198   PetscSpace     sp;
199   PetscReal     *points;
200   PetscReal     *B, *D, *H;
201   PetscInt       dim, Nb, b, Nc, c, Np, p;
202   PetscErrorCode ierr;
203 
204   PetscFunctionBeginUser;
205   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
206   ierr = DMGetField(dm, 0, NULL, (PetscObject *) &fe);CHKERRQ(ierr);
207   Np = 6;
208   ierr = PetscMalloc1(Np*dim, &points);CHKERRQ(ierr);
209   if (dim == 3) {
210     points[0]  = -1.0; points[1]  = -1.0; points[2]  = -1.0;
211     points[3]  =  1.0; points[4]  = -1.0; points[5]  = -1.0;
212     points[6]  = -1.0; points[7]  =  1.0; points[8]  = -1.0;
213     points[9]  = -1.0; points[10] = -1.0; points[11] =  1.0;
214     points[12] =  1.0; points[13] = -1.0; points[14] =  1.0;
215     points[15] = -1.0; points[16] =  1.0; points[17] =  1.0;
216   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only for 3D right now");
217   ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr);
218   ierr = PetscSpaceGetDimension(sp, &Nb);CHKERRQ(ierr);
219   ierr = PetscSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr);
220   ierr = DMGetWorkArray(dm, Np*Nb*Nc, MPIU_REAL, &B);CHKERRQ(ierr);
221   ierr = DMGetWorkArray(dm, Np*Nb*Nc*dim, MPIU_REAL, &D);CHKERRQ(ierr);
222   ierr = DMGetWorkArray(dm, Np*Nb*Nc*dim*dim, MPIU_REAL, &H);CHKERRQ(ierr);
223   ierr = PetscSpaceEvaluate(sp, Np, points, B, NULL, NULL /*D, H*/);CHKERRQ(ierr);
224   for (p = 0; p < Np; ++p) {
225     ierr = PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT "\n", p);CHKERRQ(ierr);
226     for (b = 0; b < Nb; ++b) {
227       ierr = PetscPrintf(PETSC_COMM_SELF, "B[%" PetscInt_FMT "]:", b);CHKERRQ(ierr);
228       for (c = 0; c < Nc; ++c) {ierr = PetscPrintf(PETSC_COMM_SELF, " %g", (double) B[(p*Nb+b)*Nc+c]);CHKERRQ(ierr);}
229       ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
230 #if 0
231       for (c = 0; c < Nc; ++c) {
232         ierr = PetscPrintf(PETSC_COMM_SELF, " D[%" PetscInt_FMT ",%" PetscInt_FMT "]:", b, c);CHKERRQ(ierr);
233         for (d = 0; d < dim; ++d) {ierr = PetscPrintf(PETSC_COMM_SELF, " %g", (double) B[((p*Nb+b)*Nc+c)*dim+d)]);CHKERRQ(ierr);}
234         ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
235       }
236 #endif
237     }
238   }
239   ierr = DMRestoreWorkArray(dm, Np*Nb, MPIU_REAL, &B);CHKERRQ(ierr);
240   ierr = DMRestoreWorkArray(dm, Np*Nb*dim, MPIU_REAL, &D);CHKERRQ(ierr);
241   ierr = DMRestoreWorkArray(dm, Np*Nb*dim*dim, MPIU_REAL, &H);CHKERRQ(ierr);
242   ierr = PetscFree(points);CHKERRQ(ierr);
243   PetscFunctionReturn(0);
244 }
245 
246 static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its)
247 {
248   PetscDS         ds;
249   PetscFEGeom    *chunkGeom = NULL;
250   PetscQuadrature affineQuad,  *quads = NULL;
251   PetscFEGeom    *affineGeom, **geoms = NULL;
252   PetscScalar    *u, *elemVec;
253   IS              cellIS;
254   PetscInt        depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k;
255   PetscErrorCode  ierr;
256 
257   PetscFunctionBeginUser;
258   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
259   ierr = DMGetStratumIS(dm, "depth", depth, &cellIS);CHKERRQ(ierr);
260   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
261   ierr = DMGetCellDS(dm, cStart, &ds);CHKERRQ(ierr);
262   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
263   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
264   ierr = CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms);CHKERRQ(ierr);
265   ierr = PetscMalloc2(chunkSize*totDim, &u, chunkSize*totDim, &elemVec);CHKERRQ(ierr);
266   /* Assumptions:
267     - Single field
268     - No input data
269     - No auxiliary data
270     - No time-dependence
271   */
272   for (i = 0; i < its; ++i) {
273     for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) {
274       const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS;
275 
276       ierr = PetscArrayzero(elemVec, chunkSize*totDim);CHKERRQ(ierr);
277       /* TODO Replace with DMPlexGetCellFields() */
278       for (k = 0; k < chunkSize*totDim; ++k) u[k] = 1.0;
279       for (f = 0; f < Nf; ++f) {
280         PetscFormKey key;
281         PetscFEGeom     *geom = affineGeom ? affineGeom : geoms[f];
282         /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */
283 
284         key.label = NULL; key.value = 0; key.field = f;
285         ierr = PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom);CHKERRQ(ierr);
286         ierr = PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec);CHKERRQ(ierr);
287       }
288     }
289   }
290   ierr = PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom);CHKERRQ(ierr);
291   ierr = DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms);CHKERRQ(ierr);
292   ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
293   ierr = PetscFree2(u, elemVec);CHKERRQ(ierr);
294   PetscFunctionReturn(0);
295 }
296 
297 static PetscErrorCode TestUnisolvence(DM dm)
298 {
299   Mat M;
300   Vec v;
301   PetscErrorCode ierr;
302 
303   PetscFunctionBeginUser;
304   ierr = DMGetLocalVector(dm, &v);CHKERRQ(ierr);
305   ierr = DMRestoreLocalVector(dm, &v);CHKERRQ(ierr);
306   ierr = DMCreateMassMatrix(dm, dm, &M);CHKERRQ(ierr);
307   ierr = MatViewFromOptions(M, NULL, "-mass_view");CHKERRQ(ierr);
308   ierr = MatDestroy(&M);CHKERRQ(ierr);
309   PetscFunctionReturn(0);
310 }
311 
312 int main(int argc, char **argv)
313 {
314   DM             dm;
315   AppCtx         ctx;
316   PetscMPIInt    size;
317   PetscErrorCode ierr;
318 
319   ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr;
320   ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr);
321   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only.");
322   ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
323   ierr = DMCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr);
324   ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
325   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
326   ierr = PetscObjectSetName((PetscObject) dm, "Mesh");CHKERRQ(ierr);
327   ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_view");CHKERRQ(ierr);
328   ierr = SetupDiscretization(dm, "field", SetupPrimalProblem, &ctx);CHKERRQ(ierr);
329   ierr = TestEvaluation(dm);CHKERRQ(ierr);
330   ierr = TestIntegration(dm, 1, ctx.its);CHKERRQ(ierr);
331   ierr = TestUnisolvence(dm);CHKERRQ(ierr);
332   ierr = DMDestroy(&dm);CHKERRQ(ierr);
333   ierr = PetscFinalize();
334   return ierr;
335 }
336 
337 /*TEST
338   test:
339     suffix: 0
340     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -field_petscspace_degree 0
341 
342   test:
343     suffix: 2
344     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
345           -field_petscspace_type sum \
346           -field_petscspace_variables 3 \
347           -field_petscspace_components 3 \
348           -field_petscspace_sum_spaces 2 \
349           -field_petscspace_sum_concatenate false \
350           -field_sumcomp_0_petscspace_variables 3 \
351           -field_sumcomp_0_petscspace_components 3 \
352           -field_sumcomp_0_petscspace_degree 1 \
353           -field_sumcomp_1_petscspace_variables 3 \
354           -field_sumcomp_1_petscspace_components 3 \
355           -field_sumcomp_1_petscspace_type wxy \
356           -field_petscdualspace_form_degree 0 \
357           -field_petscdualspace_order 1 \
358           -field_petscdualspace_components 3
359 
360 TEST*/
361