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