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 // PETSc FE Boilerplate 107 // ----------------------------------------------------------------------------- 108 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 109 PetscBool isSimplex, const char prefix[], 110 PetscInt order, PetscFE *fem) { 111 PetscQuadrature q, fq; 112 DM K; 113 PetscSpace P; 114 PetscDualSpace Q; 115 PetscInt quadPointsPerEdge; 116 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 117 PetscErrorCode ierr; 118 119 PetscFunctionBeginUser; 120 /* Create space */ 121 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P); CHKERRQ(ierr); 122 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix); CHKERRQ(ierr); 123 ierr = PetscSpacePolynomialSetTensor(P, tensor); CHKERRQ(ierr); 124 ierr = PetscSpaceSetFromOptions(P); CHKERRQ(ierr); 125 ierr = PetscSpaceSetNumComponents(P, Nc); CHKERRQ(ierr); 126 ierr = PetscSpaceSetNumVariables(P, dim); CHKERRQ(ierr); 127 ierr = PetscSpaceSetDegree(P, order, order); CHKERRQ(ierr); 128 ierr = PetscSpaceSetUp(P); CHKERRQ(ierr); 129 ierr = PetscSpacePolynomialGetTensor(P, &tensor); CHKERRQ(ierr); 130 /* Create dual space */ 131 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q); 132 CHKERRQ(ierr); 133 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE); CHKERRQ(ierr); 134 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix); CHKERRQ(ierr); 135 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K); CHKERRQ(ierr); 136 ierr = PetscDualSpaceSetDM(Q, K); CHKERRQ(ierr); 137 ierr = DMDestroy(&K); CHKERRQ(ierr); 138 ierr = PetscDualSpaceSetNumComponents(Q, Nc); CHKERRQ(ierr); 139 ierr = PetscDualSpaceSetOrder(Q, order); CHKERRQ(ierr); 140 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor); CHKERRQ(ierr); 141 ierr = PetscDualSpaceSetFromOptions(Q); CHKERRQ(ierr); 142 ierr = PetscDualSpaceSetUp(Q); CHKERRQ(ierr); 143 /* Create element */ 144 ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), fem); CHKERRQ(ierr); 145 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix); CHKERRQ(ierr); 146 ierr = PetscFESetFromOptions(*fem); CHKERRQ(ierr); 147 ierr = PetscFESetBasisSpace(*fem, P); CHKERRQ(ierr); 148 ierr = PetscFESetDualSpace(*fem, Q); CHKERRQ(ierr); 149 ierr = PetscFESetNumComponents(*fem, Nc); CHKERRQ(ierr); 150 ierr = PetscFESetUp(*fem); CHKERRQ(ierr); 151 ierr = PetscSpaceDestroy(&P); CHKERRQ(ierr); 152 ierr = PetscDualSpaceDestroy(&Q); CHKERRQ(ierr); 153 /* Create quadrature */ 154 quadPointsPerEdge = PetscMax(order + 1,1); 155 if (isSimplex) { 156 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, 157 &q); CHKERRQ(ierr); 158 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, 159 &fq); CHKERRQ(ierr); 160 } else { 161 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, 162 &q); CHKERRQ(ierr); 163 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, 164 &fq); CHKERRQ(ierr); 165 } 166 ierr = PetscFESetQuadrature(*fem, q); CHKERRQ(ierr); 167 ierr = PetscFESetFaceQuadrature(*fem, fq); CHKERRQ(ierr); 168 ierr = PetscQuadratureDestroy(&q); CHKERRQ(ierr); 169 ierr = PetscQuadratureDestroy(&fq); CHKERRQ(ierr); 170 171 PetscFunctionReturn(0); 172 }; 173 174 // ----------------------------------------------------------------------------- 175 // Create BC label 176 // ----------------------------------------------------------------------------- 177 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) { 178 int ierr; 179 DMLabel label; 180 181 PetscFunctionBeginUser; 182 183 ierr = DMCreateLabel(dm, name); CHKERRQ(ierr); 184 ierr = DMGetLabel(dm, name, &label); CHKERRQ(ierr); 185 ierr = DMPlexMarkBoundaryFaces(dm, 1, label); CHKERRQ(ierr); 186 ierr = DMPlexLabelComplete(dm, label); CHKERRQ(ierr); 187 188 PetscFunctionReturn(0); 189 }; 190 191 // ----------------------------------------------------------------------------- 192 // This function sets up a DM for a given degree 193 // ----------------------------------------------------------------------------- 194 PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt num_comp_u, 195 PetscInt dim, bool enforce_bc, BCFunction bc_func) { 196 PetscInt ierr, marker_ids[1] = {1}; 197 PetscFE fe; 198 199 PetscFunctionBeginUser; 200 201 // Setup FE 202 ierr = PetscFECreateByDegree(dm, dim, num_comp_u, PETSC_FALSE, NULL, degree, 203 &fe); 204 CHKERRQ(ierr); 205 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 206 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 207 208 // Setup DM 209 ierr = DMCreateDS(dm); CHKERRQ(ierr); 210 if (enforce_bc) { 211 PetscBool has_label; 212 DMHasLabel(dm, "marker", &has_label); 213 if (!has_label) {CreateBCLabel(dm, "marker");} 214 DMLabel label; 215 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 216 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, "marker", 1, 217 marker_ids, 0, 0, NULL, 218 (void(*)(void))bc_func, NULL, NULL, NULL); 219 CHKERRQ(ierr); 220 } 221 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 222 CHKERRQ(ierr); 223 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 224 225 PetscFunctionReturn(0); 226 }; 227 228 // ----------------------------------------------------------------------------- 229 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 230 // ----------------------------------------------------------------------------- 231 PetscInt Involute(PetscInt i) { 232 return i >= 0 ? i : -(i + 1); 233 }; 234 235 // ----------------------------------------------------------------------------- 236 // Get CEED restriction data from DMPlex 237 // ----------------------------------------------------------------------------- 238 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 239 CeedInt topo_dim, CeedInt height, DMLabel domain_label, CeedInt value, 240 CeedElemRestriction *elem_restr) { 241 PetscSection section; 242 PetscInt p, num_elem, num_dof, *elem_restr_offsets, e_offset, num_fields, dim, 243 depth; 244 DMLabel depth_label; 245 IS depth_is, iter_is; 246 Vec U_loc; 247 const PetscInt *iter_indices; 248 PetscErrorCode ierr; 249 250 PetscFunctionBeginUser; 251 252 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 253 dim -= height; 254 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 255 ierr = PetscSectionGetNumFields(section, &num_fields); CHKERRQ(ierr); 256 PetscInt num_comp[num_fields], field_off[num_fields+1]; 257 field_off[0] = 0; 258 for (PetscInt f = 0; f < num_fields; f++) { 259 ierr = PetscSectionGetFieldComponents(section, f, &num_comp[f]); CHKERRQ(ierr); 260 field_off[f+1] = field_off[f] + num_comp[f]; 261 } 262 263 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 264 ierr = DMPlexGetDepthLabel(dm, &depth_label); CHKERRQ(ierr); 265 ierr = DMLabelGetStratumIS(depth_label, depth - height, &depth_is); 266 CHKERRQ(ierr); 267 if (domain_label) { 268 IS domain_is; 269 ierr = DMLabelGetStratumIS(domain_label, value, &domain_is); CHKERRQ(ierr); 270 if (domain_is) { // domain_is is non-empty 271 ierr = ISIntersect(depth_is, domain_is, &iter_is); CHKERRQ(ierr); 272 ierr = ISDestroy(&domain_is); CHKERRQ(ierr); 273 } else { // domain_is is NULL (empty) 274 iter_is = NULL; 275 } 276 ierr = ISDestroy(&depth_is); CHKERRQ(ierr); 277 } else { 278 iter_is = depth_is; 279 } 280 if (iter_is) { 281 ierr = ISGetLocalSize(iter_is, &num_elem); CHKERRQ(ierr); 282 ierr = ISGetIndices(iter_is, &iter_indices); CHKERRQ(ierr); 283 } else { 284 num_elem = 0; 285 iter_indices = NULL; 286 } 287 ierr = PetscMalloc1(num_elem*PetscPowInt(P, topo_dim), &elem_restr_offsets); 288 CHKERRQ(ierr); 289 for (p = 0, e_offset = 0; p < num_elem; p++) { 290 PetscInt c = iter_indices[p]; 291 PetscInt num_indices, *indices, num_nodes; 292 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 293 &num_indices, &indices, NULL, NULL); 294 CHKERRQ(ierr); 295 bool flip = false; 296 if (height > 0) { 297 PetscInt num_cells, num_faces, start = -1; 298 const PetscInt *orients, *faces, *cells; 299 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 300 ierr = DMPlexGetSupportSize(dm, c, &num_cells); CHKERRQ(ierr); 301 if (num_cells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 302 "Expected one cell in support of exterior face, but got %D cells", 303 num_cells); 304 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 305 ierr = DMPlexGetConeSize(dm, cells[0], &num_faces); CHKERRQ(ierr); 306 for (PetscInt i=0; i<num_faces; i++) {if (faces[i] == c) start = i;} 307 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, 308 "Could not find face %D in cone of its support", 309 c); 310 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 311 if (orients[start] < 0) flip = true; 312 } 313 if (num_indices % field_off[num_fields]) SETERRQ1(PETSC_COMM_SELF, 314 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 315 c); 316 num_nodes = num_indices / field_off[num_fields]; 317 for (PetscInt i = 0; i < num_nodes; i++) { 318 PetscInt ii = i; 319 if (flip) { 320 if (P == num_nodes) ii = num_nodes - 1 - i; 321 else if (P*P == num_nodes) { 322 PetscInt row = i / P, col = i % P; 323 ii = row + col * P; 324 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, 325 "No support for flipping point with %D nodes != P (%D) or P^2", 326 num_nodes, P); 327 } 328 // Check that indices are blocked by node and thus can be coalesced as a single field with 329 // field_off[num_fields] = sum(num_comp) components. 330 for (PetscInt f = 0; f < num_fields; f++) { 331 for (PetscInt j = 0; j < num_comp[f]; j++) { 332 if (Involute(indices[field_off[f]*num_nodes + ii*num_comp[f] + j]) 333 != Involute(indices[ii*num_comp[0]]) + field_off[f] + j) 334 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 335 "Cell %D closure indices not interlaced for node %D field %D component %D", 336 c, ii, f, j); 337 } 338 } 339 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 340 PetscInt loc = Involute(indices[ii*num_comp[0]]); 341 elem_restr_offsets[e_offset++] = loc; 342 } 343 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 344 &num_indices, &indices, NULL, NULL); 345 CHKERRQ(ierr); 346 } 347 if (e_offset != num_elem*PetscPowInt(P, topo_dim)) 348 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 349 "ElemRestriction of size (%D,%D) initialized %D nodes", num_elem, 350 PetscPowInt(P, topo_dim),e_offset); 351 if (iter_is) { 352 ierr = ISRestoreIndices(iter_is, &iter_indices); CHKERRQ(ierr); 353 } 354 ierr = ISDestroy(&iter_is); CHKERRQ(ierr); 355 356 ierr = DMGetLocalVector(dm, &U_loc); CHKERRQ(ierr); 357 ierr = VecGetLocalSize(U_loc, &num_dof); CHKERRQ(ierr); 358 ierr = DMRestoreLocalVector(dm, &U_loc); CHKERRQ(ierr); 359 CeedElemRestrictionCreate(ceed, num_elem, PetscPowInt(P, topo_dim), 360 field_off[num_fields], 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 361 elem_restr_offsets, elem_restr); 362 ierr = PetscFree(elem_restr_offsets); CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 }; 365 366 // ----------------------------------------------------------------------------- 367