1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 /*@C 4 PetscFEGeomCreate - Create a `PetscFEGeom` object to manage geometry for a group of cells 5 6 Input Parameters: 7 + quad - A `PetscQuadrature` determining the tabulation 8 . numCells - The number of cells in the group 9 . dimEmbed - The coordinate dimension 10 - faceData - Flag to construct geometry data for the faces 11 12 Output Parameter: 13 . geom - The `PetscFEGeom` object 14 15 Level: beginner 16 17 .seealso: `PetscFEGeom`, `PetscQuadrature`, `PetscFEGeomDestroy()`, `PetscFEGeomComplete()` 18 @*/ 19 PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom) 20 { 21 PetscFEGeom *g; 22 PetscInt dim, Nq, N; 23 const PetscReal *p; 24 25 PetscFunctionBegin; 26 PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, &p, NULL)); 27 PetscCall(PetscNew(&g)); 28 g->xi = p; 29 g->numCells = numCells; 30 g->numPoints = Nq; 31 g->dim = dim; 32 g->dimEmbed = dimEmbed; 33 g->isCohesive = PETSC_FALSE; 34 N = numCells * Nq; 35 PetscCall(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ)); 36 if (faceData) { 37 PetscCall(PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n)); 38 PetscCall(PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]), N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]), N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]))); 39 } 40 PetscCall(PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ)); 41 *geom = g; 42 PetscFunctionReturn(PETSC_SUCCESS); 43 } 44 45 /*@C 46 PetscFEGeomDestroy - Destroy a `PetscFEGeom` object 47 48 Input Parameter: 49 . geom - `PetscFEGeom` object 50 51 Level: beginner 52 53 .seealso: `PetscFEGeom`, `PetscFEGeomCreate()` 54 @*/ 55 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 56 { 57 PetscFunctionBegin; 58 if (!*geom) PetscFunctionReturn(PETSC_SUCCESS); 59 PetscCall(PetscFree3((*geom)->v, (*geom)->J, (*geom)->detJ)); 60 PetscCall(PetscFree((*geom)->invJ)); 61 PetscCall(PetscFree2((*geom)->face, (*geom)->n)); 62 PetscCall(PetscFree6((*geom)->suppJ[0], (*geom)->suppJ[1], (*geom)->suppInvJ[0], (*geom)->suppInvJ[1], (*geom)->suppDetJ[0], (*geom)->suppDetJ[1])); 63 PetscCall(PetscFree(*geom)); 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 /*@C 68 PetscFEGeomGetChunk - Get a chunk of cells in the group as a `PetscFEGeom` 69 70 Input Parameters: 71 + geom - `PetscFEGeom` object 72 . cStart - The first cell in the chunk 73 - cEnd - The first cell not in the chunk 74 75 Output Parameter: 76 . chunkGeom - The chunk of cells 77 78 Level: intermediate 79 80 Note: 81 Use `PetscFEGeomRestoreChunk()` to return the result 82 83 .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 84 @*/ 85 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 86 { 87 PetscInt Nq; 88 PetscInt dE; 89 90 PetscFunctionBegin; 91 PetscAssertPointer(geom, 1); 92 PetscAssertPointer(chunkGeom, 4); 93 if (!(*chunkGeom)) PetscCall(PetscNew(chunkGeom)); 94 Nq = geom->numPoints; 95 dE = geom->dimEmbed; 96 (*chunkGeom)->dim = geom->dim; 97 (*chunkGeom)->dimEmbed = geom->dimEmbed; 98 (*chunkGeom)->numPoints = geom->numPoints; 99 (*chunkGeom)->numCells = cEnd - cStart; 100 (*chunkGeom)->xi = geom->xi; 101 (*chunkGeom)->v = &geom->v[Nq * dE * cStart]; 102 (*chunkGeom)->J = &geom->J[Nq * dE * dE * cStart]; 103 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq * dE * dE * cStart] : NULL; 104 (*chunkGeom)->detJ = &geom->detJ[Nq * cStart]; 105 (*chunkGeom)->n = geom->n ? &geom->n[Nq * dE * cStart] : NULL; 106 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 107 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq * dE * dE * cStart] : NULL; 108 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq * dE * dE * cStart] : NULL; 109 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq * dE * dE * cStart] : NULL; 110 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq * dE * dE * cStart] : NULL; 111 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq * cStart] : NULL; 112 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq * cStart] : NULL; 113 (*chunkGeom)->isAffine = geom->isAffine; 114 PetscFunctionReturn(PETSC_SUCCESS); 115 } 116 117 /*@C 118 PetscFEGeomRestoreChunk - Restore the chunk obtained with `PetscFEGeomCreateChunk()` 119 120 Input Parameters: 121 + geom - `PetscFEGeom` object 122 . cStart - The first cell in the chunk 123 . cEnd - The first cell not in the chunk 124 - chunkGeom - The chunk of cells 125 126 Level: intermediate 127 128 .seealso: `PetscFEGeom`, `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()` 129 @*/ 130 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 131 { 132 PetscFunctionBegin; 133 PetscCall(PetscFree(*chunkGeom)); 134 PetscFunctionReturn(PETSC_SUCCESS); 135 } 136 137 /*@C 138 PetscFEGeomGetPoint - Get the geometry for cell c at point p as a `PetscFEGeom` 139 140 Input Parameters: 141 + geom - `PetscFEGeom` object 142 . c - The cell 143 . p - The point 144 - pcoords - The reference coordinates of point p, or NULL 145 146 Output Parameter: 147 . pgeom - The geometry of cell c at point p 148 149 Level: intermediate 150 151 Notes: 152 For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 153 nothing needs to be done with it afterwards. 154 155 In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in. 156 157 .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 158 @*/ 159 PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom) 160 { 161 const PetscInt dim = geom->dim; 162 const PetscInt dE = geom->dimEmbed; 163 const PetscInt Np = geom->numPoints; 164 165 PetscFunctionBeginHot; 166 pgeom->dim = dim; 167 pgeom->dimEmbed = dE; 168 //pgeom->isAffine = geom->isAffine; 169 if (geom->isAffine) { 170 if (!p) { 171 pgeom->xi = geom->xi; 172 pgeom->J = &geom->J[c * Np * dE * dE]; 173 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 174 pgeom->detJ = &geom->detJ[c * Np]; 175 pgeom->n = geom->n ? &geom->n[c * Np * dE] : NULL; 176 } 177 if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v); 178 } else { 179 pgeom->v = &geom->v[(c * Np + p) * dE]; 180 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 181 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 182 pgeom->detJ = &geom->detJ[c * Np + p]; 183 pgeom->n = geom->n ? &geom->n[(c * Np + p) * dE] : NULL; 184 } 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187 188 /*@C 189 PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a `PetscFEGeom` 190 191 Input Parameters: 192 + geom - `PetscFEGeom` object 193 . c - The face 194 - p - The point 195 196 Output Parameter: 197 . pgeom - The cell geometry of face f at point p 198 199 Level: intermediate 200 201 Note: 202 For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 203 nothing needs to be done with it afterwards. 204 205 .seealso: `PetscFEGeom()`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 206 @*/ 207 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom) 208 { 209 const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE; 210 const PetscInt dim = bd ? geom->dimEmbed : geom->dim; 211 const PetscInt dE = geom->dimEmbed; 212 const PetscInt Np = geom->numPoints; 213 214 PetscFunctionBeginHot; 215 pgeom->dim = dim; 216 pgeom->dimEmbed = dE; 217 //pgeom->isAffine = geom->isAffine; 218 if (geom->isAffine) { 219 if (!p) { 220 if (bd) { 221 pgeom->J = &geom->suppJ[0][c * Np * dE * dE]; 222 pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE]; 223 pgeom->detJ = &geom->suppDetJ[0][c * Np]; 224 } else { 225 pgeom->J = &geom->J[c * Np * dE * dE]; 226 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 227 pgeom->detJ = &geom->detJ[c * Np]; 228 } 229 } 230 } else { 231 if (bd) { 232 pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE]; 233 pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE]; 234 pgeom->detJ = &geom->suppDetJ[0][c * Np + p]; 235 } else { 236 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 237 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 238 pgeom->detJ = &geom->detJ[c * Np + p]; 239 } 240 } 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 /*@ 245 PetscFEGeomComplete - Calculate derived quantities from base geometry specification 246 247 Input Parameter: 248 . geom - `PetscFEGeom` object 249 250 Level: intermediate 251 252 .seealso: `PetscFEGeom`, `PetscFEGeomCreate()` 253 @*/ 254 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 255 { 256 PetscInt i, j, N, dE; 257 258 PetscFunctionBeginHot; 259 N = geom->numPoints * geom->numCells; 260 dE = geom->dimEmbed; 261 switch (dE) { 262 case 3: 263 for (i = 0; i < N; i++) { 264 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 265 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 266 } 267 break; 268 case 2: 269 for (i = 0; i < N; i++) { 270 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 271 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 272 } 273 break; 274 case 1: 275 for (i = 0; i < N; i++) { 276 geom->detJ[i] = PetscAbsReal(geom->J[i]); 277 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 278 } 279 break; 280 } 281 if (geom->n) { 282 for (i = 0; i < N; i++) { 283 for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.); 284 } 285 } 286 PetscFunctionReturn(PETSC_SUCCESS); 287 } 288