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 - mode - Type of geometry data to store 11 12 Output Parameter: 13 . geom - The `PetscFEGeom` object, which is a struct not a `PetscObject` 14 15 Level: beginner 16 17 .seealso: `PetscFEGeom`, `PetscQuadrature`, `PetscFEGeomDestroy()`, `PetscFEGeomComplete()` 18 @*/ 19 PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscFEGeomMode mode, 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->mode = mode; 29 g->xi = p; 30 g->numCells = numCells; 31 g->numPoints = Nq; 32 g->dim = dim; 33 g->dimEmbed = dimEmbed; 34 N = numCells * Nq; 35 PetscCall(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ)); 36 if (mode == PETSC_FEGEOM_BOUNDARY || mode == PETSC_FEGEOM_COHESIVE) { 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 - an array of cells of length `cEnd` - `cStart` 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)->mode = geom->mode; 97 (*chunkGeom)->dim = geom->dim; 98 (*chunkGeom)->dimEmbed = geom->dimEmbed; 99 (*chunkGeom)->numPoints = geom->numPoints; 100 (*chunkGeom)->numCells = cEnd - cStart; 101 (*chunkGeom)->xi = geom->xi; 102 (*chunkGeom)->v = PetscSafePointerPlusOffset(geom->v, Nq * dE * cStart); 103 (*chunkGeom)->J = PetscSafePointerPlusOffset(geom->J, Nq * dE * dE * cStart); 104 (*chunkGeom)->invJ = PetscSafePointerPlusOffset(geom->invJ, Nq * dE * dE * cStart); 105 (*chunkGeom)->detJ = PetscSafePointerPlusOffset(geom->detJ, Nq * cStart); 106 (*chunkGeom)->n = PetscSafePointerPlusOffset(geom->n, Nq * dE * cStart); 107 (*chunkGeom)->face = PetscSafePointerPlusOffset(geom->face, cStart); 108 (*chunkGeom)->suppJ[0] = PetscSafePointerPlusOffset(geom->suppJ[0], Nq * dE * dE * cStart); 109 (*chunkGeom)->suppJ[1] = PetscSafePointerPlusOffset(geom->suppJ[1], Nq * dE * dE * cStart); 110 (*chunkGeom)->suppInvJ[0] = PetscSafePointerPlusOffset(geom->suppInvJ[0], Nq * dE * dE * cStart); 111 (*chunkGeom)->suppInvJ[1] = PetscSafePointerPlusOffset(geom->suppInvJ[1], Nq * dE * dE * cStart); 112 (*chunkGeom)->suppDetJ[0] = PetscSafePointerPlusOffset(geom->suppDetJ[0], Nq * cStart); 113 (*chunkGeom)->suppDetJ[1] = PetscSafePointerPlusOffset(geom->suppDetJ[1], Nq * cStart); 114 (*chunkGeom)->isAffine = geom->isAffine; 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 /*@C 119 PetscFEGeomRestoreChunk - Restore the chunk obtained with `PetscFEGeomCreateChunk()` 120 121 Input Parameters: 122 + geom - `PetscFEGeom` object 123 . cStart - The first cell in the chunk 124 . cEnd - The first cell not in the chunk 125 - chunkGeom - The chunk of cells 126 127 Level: intermediate 128 129 .seealso: `PetscFEGeom`, `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()` 130 @*/ 131 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 132 { 133 PetscFunctionBegin; 134 PetscCall(PetscFree(*chunkGeom)); 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 /*@C 139 PetscFEGeomGetPoint - Get the geometry for cell `c` at point `p` as a `PetscFEGeom` 140 141 Input Parameters: 142 + geom - `PetscFEGeom` object 143 . c - The cell 144 . p - The point 145 - pcoords - The reference coordinates of point `p`, or `NULL` 146 147 Output Parameter: 148 . pgeom - The geometry of cell `c` at point `p` 149 150 Level: intermediate 151 152 Notes: 153 For affine geometries, this only copies to `pgeom` at point 0. Since we copy pointers into `pgeom`, 154 nothing needs to be done with it afterwards. 155 156 In the affine case, `pgeom` must have storage for the integration point coordinates in pgeom->v if `pcoords` is passed in. 157 158 .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 159 @*/ 160 PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom) 161 { 162 const PetscInt dim = geom->dim; 163 const PetscInt dE = geom->dimEmbed; 164 const PetscInt Np = geom->numPoints; 165 166 PetscFunctionBeginHot; 167 pgeom->mode = geom->mode; 168 pgeom->dim = dim; 169 pgeom->dimEmbed = dE; 170 //pgeom->isAffine = geom->isAffine; 171 if (geom->isAffine) { 172 if (!p) { 173 pgeom->xi = geom->xi; 174 pgeom->J = &geom->J[c * Np * dE * dE]; 175 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 176 pgeom->detJ = &geom->detJ[c * Np]; 177 pgeom->n = PetscSafePointerPlusOffset(geom->n, c * Np * dE); 178 } 179 if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v); 180 } else { 181 pgeom->v = &geom->v[(c * Np + p) * dE]; 182 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 183 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 184 pgeom->detJ = &geom->detJ[c * Np + p]; 185 pgeom->n = PetscSafePointerPlusOffset(geom->n, (c * Np + p) * dE); 186 } 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 190 /*@C 191 PetscFEGeomGetCellPoint - Get the cell geometry for cell `c` at point `p` as a `PetscFEGeom` 192 193 Input Parameters: 194 + geom - `PetscFEGeom` object 195 . c - The cell 196 - p - The point 197 198 Output Parameter: 199 . pgeom - The cell geometry of cell `c` at point `p` 200 201 Level: intermediate 202 203 Notes: 204 For PETSC_FEGEOM_BOUNDARY mode, this gives the geometry for supporting cell 0. For PETSC_FEGEOM_COHESIVE mode, 205 this gives the bulk geometry for that internal face. 206 207 For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into `pgeom`, 208 nothing needs to be done with it afterwards. 209 210 .seealso: `PetscFEGeom`, `PetscFEGeomMode`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 211 @*/ 212 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom) 213 { 214 const PetscBool bd = geom->mode == PETSC_FEGEOM_BOUNDARY ? PETSC_TRUE : PETSC_FALSE; 215 const PetscInt dim = bd ? geom->dimEmbed : geom->dim; 216 const PetscInt dE = geom->dimEmbed; 217 const PetscInt Np = geom->numPoints; 218 219 PetscFunctionBeginHot; 220 pgeom->mode = geom->mode; 221 pgeom->dim = dim; 222 pgeom->dimEmbed = dE; 223 //pgeom->isAffine = geom->isAffine; 224 if (geom->isAffine) { 225 if (!p) { 226 if (bd) { 227 pgeom->J = &geom->suppJ[0][c * Np * dE * dE]; 228 pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE]; 229 pgeom->detJ = &geom->suppDetJ[0][c * Np]; 230 } else { 231 pgeom->J = &geom->J[c * Np * dE * dE]; 232 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 233 pgeom->detJ = &geom->detJ[c * Np]; 234 } 235 } 236 } else { 237 if (bd) { 238 pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE]; 239 pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE]; 240 pgeom->detJ = &geom->suppDetJ[0][c * Np + p]; 241 } else { 242 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 243 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 244 pgeom->detJ = &geom->detJ[c * Np + p]; 245 } 246 } 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 /*@C 251 PetscFEGeomComplete - Calculate derived quantities from a base geometry specification 252 253 Input Parameter: 254 . geom - `PetscFEGeom` object 255 256 Level: intermediate 257 258 .seealso: `PetscFEGeom`, `PetscFEGeomCreate()` 259 @*/ 260 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 261 { 262 PetscInt i, j, N, dE; 263 264 PetscFunctionBeginHot; 265 N = geom->numPoints * geom->numCells; 266 dE = geom->dimEmbed; 267 switch (dE) { 268 case 3: 269 for (i = 0; i < N; i++) { 270 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 271 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 272 } 273 break; 274 case 2: 275 for (i = 0; i < N; i++) { 276 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 277 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 278 } 279 break; 280 case 1: 281 for (i = 0; i < N; i++) { 282 geom->detJ[i] = PetscAbsReal(geom->J[i]); 283 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 284 } 285 break; 286 } 287 if (geom->n) { 288 for (i = 0; i < N; i++) { 289 for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.); 290 } 291 } 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294