1 /// @file 2 /// DM setup for solid mechanics example using PETSc 3 4 #include "../include/boundary.h" 5 #include "../include/setup-dm.h" 6 #include "../include/petsc-macros.h" 7 8 // ----------------------------------------------------------------------------- 9 // Setup DM 10 // ----------------------------------------------------------------------------- 11 PetscErrorCode CreateBCLabel(DM dm, const char name[]) { 12 PetscErrorCode ierr; 13 DMLabel label; 14 15 PetscFunctionBeginUser; 16 17 ierr = DMCreateLabel(dm, name); CHKERRQ(ierr); 18 ierr = DMGetLabel(dm, name, &label); CHKERRQ(ierr); 19 ierr = DMPlexMarkBoundaryFaces(dm, 1, label); CHKERRQ(ierr); 20 ierr = DMPlexLabelComplete(dm, label); CHKERRQ(ierr); 21 22 PetscFunctionReturn(0); 23 }; 24 25 // Read mesh and distribute DM in parallel 26 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx app_ctx, DM *dm) { 27 PetscErrorCode ierr; 28 const char *filename = app_ctx->mesh_file; 29 PetscBool interpolate = PETSC_TRUE; 30 DM distributed_mesh = NULL; 31 PetscPartitioner part; 32 33 PetscFunctionBeginUser; 34 35 if (!*filename) { 36 PetscInt dim = 3, faces[3] = {3, 3, 3}; 37 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", 38 faces, &dim, NULL); CHKERRQ(ierr); 39 if (!dim) dim = 3; 40 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, faces, NULL, 41 NULL, NULL, interpolate, dm); CHKERRQ(ierr); 42 } else { 43 ierr = DMPlexCreateFromFile(comm, filename, interpolate, dm); CHKERRQ(ierr); 44 } 45 46 // Distribute DM in parallel 47 ierr = DMPlexGetPartitioner(*dm, &part); CHKERRQ(ierr); 48 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 49 ierr = DMPlexDistribute(*dm, 0, NULL, &distributed_mesh); CHKERRQ(ierr); 50 if (distributed_mesh) { 51 ierr = DMDestroy(dm); CHKERRQ(ierr); 52 *dm = distributed_mesh; 53 } 54 ierr = DMViewFromOptions(*dm, NULL, "-dm_view"); CHKERRQ(ierr); 55 56 PetscFunctionReturn(0); 57 }; 58 59 // Create FE by degree 60 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 61 PetscBool isSimplex, const char prefix[], 62 PetscInt order, PetscFE *fem) { 63 PetscQuadrature q, fq; 64 DM K; 65 PetscSpace P; 66 PetscDualSpace Q; 67 PetscInt quadPointsPerEdge; 68 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 69 PetscErrorCode ierr; 70 71 PetscFunctionBeginUser; 72 /* Create space */ 73 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P); CHKERRQ(ierr); 74 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix); CHKERRQ(ierr); 75 ierr = PetscSpacePolynomialSetTensor(P, tensor); CHKERRQ(ierr); 76 ierr = PetscSpaceSetFromOptions(P); CHKERRQ(ierr); 77 ierr = PetscSpaceSetNumComponents(P, Nc); CHKERRQ(ierr); 78 ierr = PetscSpaceSetNumVariables(P, dim); CHKERRQ(ierr); 79 ierr = PetscSpaceSetDegree(P, order, order); CHKERRQ(ierr); 80 ierr = PetscSpaceSetUp(P); CHKERRQ(ierr); 81 ierr = PetscSpacePolynomialGetTensor(P, &tensor); CHKERRQ(ierr); 82 /* Create dual space */ 83 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q); 84 CHKERRQ(ierr); 85 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE); CHKERRQ(ierr); 86 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix); CHKERRQ(ierr); 87 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K); CHKERRQ(ierr); 88 ierr = PetscDualSpaceSetDM(Q, K); CHKERRQ(ierr); 89 ierr = DMDestroy(&K); CHKERRQ(ierr); 90 ierr = PetscDualSpaceSetNumComponents(Q, Nc); CHKERRQ(ierr); 91 ierr = PetscDualSpaceSetOrder(Q, order); CHKERRQ(ierr); 92 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor); CHKERRQ(ierr); 93 ierr = PetscDualSpaceSetFromOptions(Q); CHKERRQ(ierr); 94 ierr = PetscDualSpaceSetUp(Q); CHKERRQ(ierr); 95 /* Create element */ 96 ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), fem); CHKERRQ(ierr); 97 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix); CHKERRQ(ierr); 98 ierr = PetscFESetFromOptions(*fem); CHKERRQ(ierr); 99 ierr = PetscFESetBasisSpace(*fem, P); CHKERRQ(ierr); 100 ierr = PetscFESetDualSpace(*fem, Q); CHKERRQ(ierr); 101 ierr = PetscFESetNumComponents(*fem, Nc); CHKERRQ(ierr); 102 ierr = PetscFESetUp(*fem); CHKERRQ(ierr); 103 ierr = PetscSpaceDestroy(&P); CHKERRQ(ierr); 104 ierr = PetscDualSpaceDestroy(&Q); CHKERRQ(ierr); 105 /* Create quadrature */ 106 quadPointsPerEdge = PetscMax(order + 1,1); 107 if (isSimplex) { 108 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, 109 &q); CHKERRQ(ierr); 110 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, 111 &fq); CHKERRQ(ierr); 112 } else { 113 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, 114 &q); CHKERRQ(ierr); 115 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, 116 &fq); CHKERRQ(ierr); 117 } 118 ierr = PetscFESetQuadrature(*fem, q); CHKERRQ(ierr); 119 ierr = PetscFESetFaceQuadrature(*fem, fq); CHKERRQ(ierr); 120 ierr = PetscQuadratureDestroy(&q); CHKERRQ(ierr); 121 ierr = PetscQuadratureDestroy(&fq); CHKERRQ(ierr); 122 123 PetscFunctionReturn(0); 124 } 125 126 // Setup DM with FE space of appropriate degree 127 PetscErrorCode SetupDMByDegree(DM dm, AppCtx app_ctx, PetscInt order, 128 PetscBool boundary, PetscInt num_comp_u) { 129 PetscErrorCode ierr; 130 PetscInt dim; 131 PetscFE fe; 132 IS face_set_is; // Index Set for Face Sets 133 const char *name = "Face Sets"; // PETSc internal requirement 134 PetscInt num_face_sets; // Number of FaceSets in face_set_is 135 const PetscInt *face_set_ids; // id of each FaceSet 136 137 PetscFunctionBeginUser; 138 139 // Setup DM 140 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 141 ierr = PetscFECreateByDegree(dm, dim, num_comp_u, PETSC_FALSE, NULL, order, 142 &fe); 143 CHKERRQ(ierr); 144 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 145 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 146 ierr = DMCreateDS(dm); CHKERRQ(ierr); 147 148 // Add Dirichlet (Essential) boundary 149 if (boundary) { 150 if (app_ctx->forcing_choice == FORCE_MMS) { 151 if (app_ctx->test_mode) { 152 // -- Test mode - box mesh 153 PetscBool has_label; 154 PetscInt marker_ids[1] = {1}; 155 ierr = DMHasLabel(dm, "marker", &has_label); CHKERRQ(ierr); 156 if (!has_label) { 157 ierr = CreateBCLabel(dm, "marker"); CHKERRQ(ierr); 158 } 159 DMLabel label; 160 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 161 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, "marker", 1, marker_ids, 162 0, 0, NULL, (void(*)(void))BCMMS, NULL, NULL, NULL); 163 CHKERRQ(ierr); 164 } else { 165 // -- ExodusII mesh with MMS 166 ierr = DMGetLabelIdIS(dm, name, &face_set_is); CHKERRQ(ierr); 167 ierr = ISGetSize(face_set_is,&num_face_sets); CHKERRQ(ierr); 168 ierr = ISGetIndices(face_set_is, &face_set_ids); CHKERRQ(ierr); 169 DMLabel label; 170 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 171 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, "Face Sets", 172 num_face_sets, face_set_ids, 0, 0, NULL, 173 (void(*)(void))BCMMS, NULL, NULL, NULL); 174 CHKERRQ(ierr); 175 ierr = ISRestoreIndices(face_set_is, &face_set_ids); CHKERRQ(ierr); 176 ierr = ISDestroy(&face_set_is); CHKERRQ(ierr); 177 } 178 } else { 179 // -- ExodusII mesh with user specified BCs 180 // -- Clamp BCs 181 DMLabel label; 182 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 183 for (PetscInt i = 0; i < app_ctx->bc_clamp_count; i++) { 184 char bcName[25]; 185 snprintf(bcName, sizeof bcName, "clamp_%d", app_ctx->bc_clamp_faces[i]); 186 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, bcName, label, "Face Sets", 1, 187 &app_ctx->bc_clamp_faces[i], 0, 0, 188 NULL, (void(*)(void))BCClamp, NULL, 189 (void *)&app_ctx->bc_clamp_max[i], NULL); 190 CHKERRQ(ierr); 191 } 192 } 193 } 194 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 195 CHKERRQ(ierr); 196 197 // Cleanup 198 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 199 200 PetscFunctionReturn(0); 201 }; 202