1 static const char help[] = "Tests DMCreateMassMatrix and DMCreateMassMatrixLumped";
2
3 #include <petscdmplex.h>
4 #include <petscfe.h>
5 #include <petscds.h>
6
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 obj[])7 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 obj[])
8 {
9 obj[0] = 1.0;
10 }
11
SetupDiscretization(DM dm)12 static PetscErrorCode SetupDiscretization(DM dm)
13 {
14 PetscFE fe0, fe1, fe2;
15 DM plex;
16 PetscInt dim;
17 PetscBool simplex;
18 PetscDS ds;
19
20 PetscFunctionBeginUser;
21 PetscCall(DMGetDimension(dm, &dim));
22 PetscCall(DMConvert(dm, DMPLEX, &plex));
23 PetscCall(DMPlexIsSimplex(plex, &simplex));
24 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &simplex, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
25 PetscCall(DMDestroy(&plex));
26
27 /* Create finite element */
28 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, "f0_", -1, &fe0));
29 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 2, simplex, "f1_", -1, &fe1));
30 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, simplex, "f2_", -1, &fe2));
31 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe0));
32 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe1));
33 PetscCall(DMSetField(dm, 2, NULL, (PetscObject)fe2));
34 PetscCall(DMCreateDS(dm));
35 PetscCall(DMGetDS(dm, &ds));
36 PetscCall(PetscDSSetObjective(ds, 0, volume));
37 PetscCall(PetscDSSetObjective(ds, 1, volume));
38 PetscCall(PetscDSSetObjective(ds, 2, volume));
39 PetscCall(PetscFEDestroy(&fe0));
40 PetscCall(PetscFEDestroy(&fe1));
41 PetscCall(PetscFEDestroy(&fe2));
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
45 #define CheckVals(a, b, rtol, atol, msg) \
46 do { \
47 if (!PetscIsCloseAtTolScalar(a, b, rtol, atol)) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s: %g (%g - %g)\n", msg, (double)PetscAbsScalar(a - b), (double)PetscAbsScalar(a), (double)PetscAbsScalar(b))); \
48 } while (0)
49
main(int argc,char ** argv)50 int main(int argc, char **argv)
51 {
52 DM dm;
53 Mat M;
54 Vec llM, lM, ones, work;
55 PetscReal rtol = PETSC_SMALL, atol = 0.0;
56 PetscScalar vals[6];
57 PetscInt dim;
58 PetscBool amr = PETSC_FALSE;
59
60 PetscFunctionBeginUser;
61 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
62 PetscCall(PetscOptionsGetBool(NULL, NULL, "-amr", &amr, NULL));
63 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
64 PetscCall(DMSetType(dm, DMPLEX));
65 PetscCall(DMSetFromOptions(dm));
66 PetscCall(DMGetDimension(dm, &dim));
67 if (amr) {
68 DM dmConv;
69
70 PetscCall(DMConvert(dm, dim == 2 ? DMP4EST : DMP8EST, &dmConv));
71 PetscCall(DMDestroy(&dm));
72 dm = dmConv;
73 PetscCall(DMSetFromOptions(dm));
74 PetscCall(DMSetUp(dm));
75 }
76 PetscCall(DMViewFromOptions(dm, NULL, "-mesh_view"));
77 PetscCall(SetupDiscretization(dm));
78 PetscCall(DMGetGlobalVector(dm, &ones));
79 PetscCall(DMGetGlobalVector(dm, &work));
80 PetscCall(VecSet(ones, 1.0));
81 PetscCall(DMCreateMassMatrixLumped(dm, &llM, &lM));
82 PetscCall(VecViewFromOptions(llM, NULL, "-local_lumped_mass_view"));
83 PetscCall(VecViewFromOptions(lM, NULL, "-lumped_mass_view"));
84 PetscCall(DMCreateMassMatrix(dm, NULL, &M));
85 PetscCall(MatViewFromOptions(M, NULL, "-mass_view"));
86 PetscCall(MatMult(M, ones, work));
87 PetscCall(VecViewFromOptions(work, NULL, "-mass_rowsum_view"));
88 PetscCall(VecSum(work, &vals[3]));
89 PetscCall(VecSet(work, 0.0));
90 PetscCall(DMLocalToGlobal(dm, llM, ADD_VALUES, work));
91 PetscCall(VecSum(work, &vals[4]));
92 PetscCall(VecSum(lM, &vals[5]));
93 PetscCall(DMPlexComputeIntegralFEM(dm, ones, vals, NULL));
94 CheckVals(vals[0], vals[1], rtol, atol, "Error volume");
95 CheckVals((3 + dim) * vals[0], vals[3], rtol, atol, "Error mass");
96 CheckVals((3 + dim) * vals[0], vals[4], rtol, atol, "Error local lumped mass");
97 CheckVals((3 + dim) * vals[0], vals[5], rtol, atol, "Error lumped mass");
98 PetscCall(MatDestroy(&M));
99 PetscCall(VecDestroy(&lM));
100 PetscCall(VecDestroy(&llM));
101 PetscCall(DMRestoreGlobalVector(dm, &work));
102 PetscCall(DMRestoreGlobalVector(dm, &ones));
103 PetscCall(DMDestroy(&dm));
104 PetscCall(PetscFinalize());
105 return 0;
106 }
107
108 /*TEST
109
110 test:
111 suffix: 0
112 output_file: output/empty.out
113 nsize: {{1 2}}
114 args: -dm_plex_box_faces 3,3,3 -dm_plex_dim {{2 3}} -dm_plex_simplex 0 -f0_petscspace_degree {{0 1}} -f1_petscspace_degree {{0 1}} -f2_petscspace_degree {{0 1}} -petscpartitioner_type simple
115
116 test:
117 TODO: broken
118 suffix: 0_amr
119 output_file: output/empty.out
120 requires: p4est
121 nsize: {{1 2}}
122 args: -dm_plex_dim {{2 3}} -dm_plex_simplex 0 -f0_petscspace_degree {{0 1}} -f1_petscspace_degree {{0 1}} -f2_petscspace_degree {{0 1}} -amr -dm_p4est_refine_pattern hash -dm_forest_initial_refinement 0 -dm_forest_maximum_refinement 2
123
124 test:
125 suffix: 1
126 output_file: output/empty.out
127 requires: triangle
128 nsize: {{1 2}}
129 args: -dm_plex_box_faces 3,3 -dm_plex_dim 2 -dm_plex_simplex 1 -f0_petscspace_degree {{0 1}} -f1_petscspace_degree {{0 1}} -f2_petscspace_degree {{0 1}} -petscpartitioner_type simple
130
131 test:
132 suffix: 2
133 output_file: output/empty.out
134 requires: ctetgen
135 nsize: {{1 2}}
136 args: -dm_plex_box_faces 3,3,3 -dm_plex_dim 3 -dm_plex_simplex 1 -f0_petscspace_degree {{0 1}} -f1_petscspace_degree {{0 1}} -f2_petscspace_degree {{0 1}} -petscpartitioner_type simple
137
138 TEST*/
139