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