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 // Setup DM with FE space of appropriate degree 60 PetscErrorCode SetupDMByDegree(DM dm, AppCtx app_ctx, PetscInt order, 61 PetscBool boundary, PetscInt num_comp_u) { 62 PetscErrorCode ierr; 63 MPI_Comm comm; 64 PetscInt dim; 65 PetscFE fe; 66 IS face_set_is; // Index Set for Face Sets 67 const char *name = "Face Sets"; // PETSc internal requirement 68 PetscInt num_face_sets; // Number of FaceSets in face_set_is 69 const PetscInt *face_set_ids; // id of each FaceSet 70 71 PetscFunctionBeginUser; 72 73 // Setup DM 74 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 75 ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); 76 ierr = PetscFECreateLagrange(comm, dim, num_comp_u, PETSC_FALSE, order, order, 77 &fe); CHKERRQ(ierr); 78 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 79 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 80 ierr = DMCreateDS(dm); CHKERRQ(ierr); 81 82 // Add Dirichlet (Essential) boundary 83 if (boundary) { 84 if (app_ctx->forcing_choice == FORCE_MMS) { 85 if (app_ctx->test_mode) { 86 // -- Test mode - box mesh 87 PetscBool has_label; 88 PetscInt marker_ids[1] = {1}; 89 ierr = DMHasLabel(dm, "marker", &has_label); CHKERRQ(ierr); 90 if (!has_label) { 91 ierr = CreateBCLabel(dm, "marker"); CHKERRQ(ierr); 92 } 93 DMLabel label; 94 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 95 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, "marker", 1, marker_ids, 96 0, 0, NULL, (void(*)(void))BCMMS, NULL, NULL, NULL); 97 CHKERRQ(ierr); 98 } else { 99 // -- ExodusII mesh with MMS 100 ierr = DMGetLabelIdIS(dm, name, &face_set_is); CHKERRQ(ierr); 101 ierr = ISGetSize(face_set_is,&num_face_sets); CHKERRQ(ierr); 102 ierr = ISGetIndices(face_set_is, &face_set_ids); CHKERRQ(ierr); 103 DMLabel label; 104 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 105 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, "Face Sets", 106 num_face_sets, face_set_ids, 0, 0, NULL, 107 (void(*)(void))BCMMS, NULL, NULL, NULL); 108 CHKERRQ(ierr); 109 ierr = ISRestoreIndices(face_set_is, &face_set_ids); CHKERRQ(ierr); 110 ierr = ISDestroy(&face_set_is); CHKERRQ(ierr); 111 } 112 } else { 113 // -- ExodusII mesh with user specified BCs 114 // -- Clamp BCs 115 DMLabel label; 116 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 117 for (PetscInt i = 0; i < app_ctx->bc_clamp_count; i++) { 118 char bcName[25]; 119 snprintf(bcName, sizeof bcName, "clamp_%d", app_ctx->bc_clamp_faces[i]); 120 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, bcName, label, "Face Sets", 1, 121 &app_ctx->bc_clamp_faces[i], 0, 0, 122 NULL, (void(*)(void))BCClamp, NULL, 123 (void *)&app_ctx->bc_clamp_max[i], NULL); 124 CHKERRQ(ierr); 125 } 126 } 127 } 128 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 129 CHKERRQ(ierr); 130 131 // Cleanup 132 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 133 134 PetscFunctionReturn(0); 135 }; 136