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