xref: /petsc/src/dm/impls/plex/tests/ex33.c (revision bcd4bb4a4158aa96f212e9537e87b40407faf83e)
1 static char help[] = "Tests for high order geometry\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petscds.h>
5 
6 typedef struct {
7   PetscReal volume; /* Analytical volume of the mesh */
8   PetscReal tol;    /* Tolerance for volume check */
9 } AppCtx;
10 
11 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12 {
13   PetscFunctionBegin;
14   options->volume = -1.0;
15   options->tol    = PETSC_SMALL;
16 
17   PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
18   PetscCall(PetscOptionsReal("-volume", "The analytical volume of the mesh", "ex33.c", options->volume, &options->volume, NULL));
19   PetscCall(PetscOptionsReal("-tol", "The tolerance for the volume check", "ex33.c", options->tol, &options->tol, NULL));
20   PetscOptionsEnd();
21   PetscFunctionReturn(PETSC_SUCCESS);
22 }
23 
24 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
25 {
26   PetscFunctionBegin;
27   PetscCall(DMCreate(comm, dm));
28   PetscCall(DMSetType(*dm, DMPLEX));
29   PetscCall(DMSetFromOptions(*dm));
30   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
31   PetscFunctionReturn(PETSC_SUCCESS);
32 }
33 
34 static void volume(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 vol[])
35 {
36   vol[0] = 1.;
37 }
38 
39 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx)
40 {
41   PetscDS        ds;
42   PetscFE        fe;
43   DMPolytopeType ct;
44   PetscInt       dim, cStart;
45 
46   PetscFunctionBeginUser;
47   PetscCall(DMGetDimension(dm, &dim));
48   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
49   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
50   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, PETSC_DETERMINE, &fe));
51   PetscCall(PetscFESetName(fe, "scalar"));
52   PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
53   PetscCall(PetscFEDestroy(&fe));
54   PetscCall(DMCreateDS(dm));
55   PetscCall(DMGetDS(dm, &ds));
56   PetscCall(PetscDSSetObjective(ds, 0, volume));
57   PetscFunctionReturn(PETSC_SUCCESS);
58 }
59 
60 static PetscErrorCode CheckVolume(DM dm, AppCtx *ctx)
61 {
62   Vec         u;
63   PetscScalar result;
64   PetscReal   vol, tol = ctx->tol;
65 
66   PetscFunctionBeginUser;
67   PetscCall(DMGetGlobalVector(dm, &u));
68   PetscCall(DMPlexComputeIntegralFEM(dm, u, &result, ctx));
69   vol = PetscRealPart(result);
70   PetscCall(DMRestoreGlobalVector(dm, &u));
71   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Volume: %g\n", (double)vol));
72   PetscCheck(ctx->volume <= 0.0 || PetscAbsReal(ctx->volume - vol) <= tol, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Calculated volume %g != %g actual volume (error %g > %g tol)", (double)vol, (double)ctx->volume, (double)PetscAbsReal(ctx->volume - vol), (double)tol);
73   PetscFunctionReturn(PETSC_SUCCESS);
74 }
75 
76 int main(int argc, char **argv)
77 {
78   DM     dm;
79   AppCtx user;
80 
81   PetscFunctionBeginUser;
82   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
83   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
84   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
85   PetscCall(CreateDiscretization(dm, &user));
86   PetscCall(CheckVolume(dm, &user));
87   PetscCall(DMDestroy(&dm));
88   PetscCall(PetscFinalize());
89   return 0;
90 }
91 
92 /*TEST
93 
94   testset:
95     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -volume 4.
96 
97     test:
98       suffix: square_0
99       args: -dm_coord_petscspace_degree 1
100 
101     test:
102       suffix: square_1
103       args: -dm_coord_petscspace_degree 2
104 
105     test:
106       suffix: square_2
107       args: -dm_refine 1 -dm_coord_petscspace_degree 1
108 
109     test:
110       suffix: square_3
111       args: -dm_refine 1 -dm_coord_petscspace_degree 2
112 
113   testset:
114     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -volume 8.
115 
116     test:
117       suffix: cube_0
118       args: -dm_coord_petscspace_degree 1
119 
120     test:
121       suffix: cube_1
122       args: -dm_coord_petscspace_degree 2
123 
124     test:
125       suffix: cube_2
126       args: -dm_refine 1 -dm_coord_petscspace_degree 1
127 
128     test:
129       suffix: cube_3
130       args: -dm_refine 1 -dm_coord_petscspace_degree 2
131 
132   testset:
133     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. \
134           -volume 4. -dm_coord_remap -dm_coord_map shear -dm_coord_map_params 0.0,0.0,3.0
135 
136     test:
137       suffix: shear_0
138       args: -dm_coord_petscspace_degree 1
139 
140     test:
141       suffix: shear_1
142       args: -dm_coord_petscspace_degree 2
143 
144     test:
145       suffix: shear_2
146       args: -dm_refine 1 -dm_coord_petscspace_degree 1
147 
148     test:
149       suffix: shear_3
150       args: -dm_refine 1 -dm_coord_petscspace_degree 2
151 
152   testset:
153     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. \
154           -volume 8. -dm_coord_remap -dm_coord_map shear -dm_coord_map_params 0.0,0.0,3.0,4.0
155 
156     test:
157       suffix: shear_4
158       args: -dm_coord_petscspace_degree 1
159 
160     test:
161       suffix: shear_5
162       args: -dm_coord_petscspace_degree 2
163 
164     test:
165       suffix: shear_6
166       args: -dm_refine 1 -dm_coord_petscspace_degree 1
167 
168     test:
169       suffix: shear_7
170       args: -dm_refine 1 -dm_coord_petscspace_degree 2
171 
172   testset:
173     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower 1.,-1. -dm_plex_box_upper 3.,1. \
174           -dm_coord_petscspace_degree 1 -volume 8. -dm_coord_remap -dm_coord_map flare
175 
176     test:
177       suffix: flare_0
178       args:
179 
180     test:
181       suffix: flare_1
182       args: -dm_refine 2
183 
184   testset:
185     # Area: 3/4 \pi = 2.3562
186     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -petscfe_default_quadrature_order 2 \
187           -volume 2.35619449019235 -dm_coord_remap -dm_coord_map annulus
188 
189     test:
190       # Area: (a+b)/2 h = 3/\sqrt{2} (sqrt{2} - 1/\sqrt{2}) = 3/2
191       suffix: annulus_0
192       requires: double
193       args: -dm_coord_petscspace_degree 1 -volume 1.5
194 
195     test:
196       suffix: annulus_1
197       requires: double
198       args: -dm_refine 3 -dm_coord_petscspace_degree 1 -tol .016
199 
200     test:
201       suffix: annulus_2
202       requires: double
203       args: -dm_refine 3 -dm_coord_petscspace_degree 2 -tol .0038
204 
205     test:
206       suffix: annulus_3
207       requires: double
208       args: -dm_refine 3 -dm_coord_petscspace_degree 3 -tol 2.2e-6
209 
210     test:
211       suffix: annulus_4
212       requires: double
213       args: -dm_refine 2 -dm_coord_petscspace_degree 2 -tol .00012
214 
215     test:
216       suffix: annulus_5
217       requires: double
218       args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol 1.2e-7
219 
220   testset:
221     # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
222     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. \
223           -volume 14.66076571675238 -dm_coord_remap -dm_coord_map shell
224 
225     test:
226       suffix: shell_0
227       requires: double
228       args: -dm_refine 1 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -volume 5.633164922 -tol 1.0e-7
229 
230     test:
231       suffix: shell_1
232       requires: double
233       args: -dm_refine 2 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -tol 3.1
234 
235     test:
236       suffix: shell_2
237       requires: double
238       args: -dm_refine 2 -dm_coord_petscspace_degree 2 -petscfe_default_quadrature_order 2 -tol .1
239 
240     test:
241       suffix: shell_3
242       requires: double
243       args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol .02
244 
245     test:
246       suffix: shell_4
247       requires: double
248       args: -dm_refine 2 -dm_coord_petscspace_degree 4 -petscfe_default_quadrature_order 4 -tol .006
249 
250   test:
251     # Volume: 1.0
252     suffix: gmsh_q2
253     requires: double
254     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
255 
256   test:
257     # Volume: 1.0
258     suffix: gmsh_q3
259     requires: double
260     nsize: {{1 2}}
261     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
262 
263   test:
264     # Volume: 1.0
265     suffix: gmsh_3d_q2
266     requires: double
267     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/cube_q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
268 
269   test:
270     # Volume: 1.0
271     suffix: gmsh_3d_q3
272     requires: double
273     nsize: {{1 2}}
274     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/cube_q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
275 
276 TEST*/
277