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