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*/