1 #include "../include/petscutils.h" 2 3 // ----------------------------------------------------------------------------- 4 // Convert PETSc MemType to libCEED MemType 5 // ----------------------------------------------------------------------------- 6 CeedMemType MemTypeP2C(PetscMemType mem_type) { 7 return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 8 } 9 10 // ----------------------------------------------------------------------------- 11 // Apply 3D Kershaw mesh transformation 12 // ----------------------------------------------------------------------------- 13 // Transition from a value of "a" for x=0, to a value of "b" for x=1. Optionally 14 // smooth -- see the commented versions at the end. 15 static double step(const double a, const double b, double x) { 16 if (x <= 0) return a; 17 if (x >= 1) return b; 18 return a + (b-a) * (x); 19 } 20 21 // 1D transformation at the right boundary 22 static double right(const double eps, const double x) { 23 return (x <= 0.5) ? (2-eps) * x : 1 + eps*(x-1); 24 } 25 26 // 1D transformation at the left boundary 27 static double left(const double eps, const double x) { 28 return 1-right(eps,1-x); 29 } 30 31 // Apply 3D Kershaw mesh transformation 32 // The eps parameters are in (0, 1] 33 // Uniform mesh is recovered for eps=1 34 PetscErrorCode Kershaw(DM dm_orig, PetscScalar eps) { 35 PetscErrorCode ierr; 36 Vec coord; 37 PetscInt ncoord; 38 PetscScalar *c; 39 40 PetscFunctionBeginUser; 41 ierr = DMGetCoordinatesLocal(dm_orig, &coord); CHKERRQ(ierr); 42 ierr = VecGetLocalSize(coord, &ncoord); CHKERRQ(ierr); 43 ierr = VecGetArray(coord, &c); CHKERRQ(ierr); 44 45 for (PetscInt i = 0; i < ncoord; i += 3) { 46 PetscScalar x = c[i], y = c[i+1], z = c[i+2]; 47 PetscInt layer = x*6; 48 PetscScalar lambda = (x-layer/6.0)*6; 49 c[i] = x; 50 51 switch (layer) { 52 case 0: 53 c[i+1] = left(eps, y); 54 c[i+2] = left(eps, z); 55 break; 56 case 1: 57 case 4: 58 c[i+1] = step(left(eps, y), right(eps, y), lambda); 59 c[i+2] = step(left(eps, z), right(eps, z), lambda); 60 break; 61 case 2: 62 c[i+1] = step(right(eps, y), left(eps, y), lambda/2); 63 c[i+2] = step(right(eps, z), left(eps, z), lambda/2); 64 break; 65 case 3: 66 c[i+1] = step(right(eps, y), left(eps, y), (1+lambda)/2); 67 c[i+2] = step(right(eps, z), left(eps, z), (1+lambda)/2); 68 break; 69 default: 70 c[i+1] = right(eps, y); 71 c[i+2] = right(eps, z); 72 } 73 } 74 ierr = VecRestoreArray(coord, &c); CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 // ----------------------------------------------------------------------------- 79 // Create BC label 80 // ----------------------------------------------------------------------------- 81 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) { 82 int ierr; 83 DMLabel label; 84 85 PetscFunctionBeginUser; 86 87 ierr = DMCreateLabel(dm, name); CHKERRQ(ierr); 88 ierr = DMGetLabel(dm, name, &label); CHKERRQ(ierr); 89 ierr = DMPlexMarkBoundaryFaces(dm, 1, label); CHKERRQ(ierr); 90 91 PetscFunctionReturn(0); 92 }; 93 94 // ----------------------------------------------------------------------------- 95 // This function sets up a DM for a given degree 96 // ----------------------------------------------------------------------------- 97 PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt num_comp_u, 98 PetscInt dim, bool enforce_bc, BCFunction bc_func) { 99 PetscInt ierr, marker_ids[1] = {1}; 100 PetscFE fe; 101 MPI_Comm comm; 102 103 PetscFunctionBeginUser; 104 105 // Setup FE 106 ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); 107 ierr = PetscFECreateLagrange(comm, dim, num_comp_u, PETSC_FALSE, degree, degree, 108 &fe); CHKERRQ(ierr); 109 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 110 { 111 /* create FE field for coordinates */ 112 PetscFE fe_coords; 113 PetscInt num_comp_coord; 114 ierr = DMGetCoordinateDim(dm, &num_comp_coord); CHKERRQ(ierr); 115 ierr = PetscFECreateLagrange(comm, dim, num_comp_coord, PETSC_FALSE, 1, 1, 116 &fe_coords); CHKERRQ(ierr); 117 ierr = DMProjectCoordinates(dm, fe_coords); CHKERRQ(ierr); 118 ierr = PetscFEDestroy(&fe_coords); CHKERRQ(ierr); 119 } 120 121 // Setup DM 122 ierr = DMCreateDS(dm); CHKERRQ(ierr); 123 if (enforce_bc) { 124 PetscBool has_label; 125 DMHasLabel(dm, "marker", &has_label); 126 if (!has_label) {CreateBCLabel(dm, "marker");} 127 DMLabel label; 128 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 129 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, 130 marker_ids, 0, 0, NULL, (void(*)(void))bc_func, 131 NULL, NULL, NULL); CHKERRQ(ierr); 132 } 133 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 134 CHKERRQ(ierr); 135 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 136 137 PetscFunctionReturn(0); 138 }; 139 140 // ----------------------------------------------------------------------------- 141 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 142 // ----------------------------------------------------------------------------- 143 PetscInt Involute(PetscInt i) { 144 return i >= 0 ? i : -(i + 1); 145 }; 146 147 // ----------------------------------------------------------------------------- 148 // Get CEED restriction data from DMPlex 149 // ----------------------------------------------------------------------------- 150 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, 151 DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { 152 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; 153 PetscErrorCode ierr; 154 155 PetscFunctionBeginUser; 156 157 ierr = DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, 158 &elem_size, &num_comp, &num_dof, &elem_restr_offsets); 159 CHKERRQ(ierr); 160 161 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 162 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 163 elem_restr_offsets, elem_restr); 164 ierr = PetscFree(elem_restr_offsets); CHKERRQ(ierr); 165 166 PetscFunctionReturn(0); 167 }; 168 169 // ----------------------------------------------------------------------------- 170