Actual source code: ex4.c
1: static const char help[] = "Tests DMCreateMassMatrix and DMCreateMassMatrixLumped";
3: #include <petscdmplex.h>
4: #include <petscfe.h>
5: #include <petscds.h>
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: }
12: static PetscErrorCode SetupDiscretization(DM dm)
13: {
14: PetscFE fe0, fe1, fe2;
15: DM plex;
16: PetscInt dim;
17: PetscBool simplex;
18: PetscDS ds;
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, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
25: PetscCall(DMDestroy(&plex));
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: }
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)
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;
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;
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: }
108: /*TEST
110: test:
111: suffix: 0
112: output_file: output/ex4_ok.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
116: test:
117: TODO: broken
118: suffix: 0_amr
119: output_file: output/ex4_ok.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
124: test:
125: suffix: 1
126: output_file: output/ex4_ok.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
131: test:
132: suffix: 2
133: output_file: output/ex4_ok.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
138: TEST*/