1 static const char help[] = "Tests DMCreateMassMatrix and DMCreateMassMatrixLumped"; 2 3 #include <petscdmplex.h> 4 #include <petscfe.h> 5 #include <petscds.h> 6 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 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 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