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: 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]), 39 N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]), 40 N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]))); 41 } 42 PetscCall(PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ)); 43 *geom = g; 44 PetscFunctionReturn(0); 45 } 46 47 /*@C 48 PetscFEGeomDestroy - Destroy a PetscFEGeom object 49 50 Input Parameter: 51 . geom - PetscFEGeom object 52 53 Level: beginner 54 55 .seealso: PetscFEGeomCreate() 56 @*/ 57 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 58 { 59 PetscFunctionBegin; 60 if (!*geom) PetscFunctionReturn(0); 61 PetscCall(PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ)); 62 PetscCall(PetscFree((*geom)->invJ)); 63 PetscCall(PetscFree2((*geom)->face,(*geom)->n)); 64 PetscCall(PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1])); 65 PetscCall(PetscFree(*geom)); 66 PetscFunctionReturn(0); 67 } 68 69 /*@C 70 PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom 71 72 Input Parameters: 73 + geom - PetscFEGeom object 74 . cStart - The first cell in the chunk 75 - cEnd - The first cell not in the chunk 76 77 Output Parameter: 78 . chunkGeom - The chunk of cells 79 80 Level: intermediate 81 82 .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate() 83 @*/ 84 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 85 { 86 PetscInt Nq; 87 PetscInt dE; 88 89 PetscFunctionBegin; 90 PetscValidPointer(geom,1); 91 PetscValidPointer(chunkGeom,4); 92 if (!(*chunkGeom)) { 93 PetscCall(PetscNew(chunkGeom)); 94 } 95 Nq = geom->numPoints; 96 dE= geom->dimEmbed; 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 = &geom->v[Nq*dE*cStart]; 103 (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart]; 104 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL; 105 (*chunkGeom)->detJ = &geom->detJ[Nq*cStart]; 106 (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL; 107 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 108 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL; 109 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL; 110 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL; 111 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL; 112 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL; 113 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL; 114 (*chunkGeom)->isAffine = geom->isAffine; 115 PetscFunctionReturn(0); 116 } 117 118 /*@C 119 PetscFEGeomRestoreChunk - Restore the chunk 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: PetscFEGeomGetChunk(), PetscFEGeomCreate() 130 @*/ 131 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 132 { 133 PetscFunctionBegin; 134 PetscCall(PetscFree(*chunkGeom)); 135 PetscFunctionReturn(0); 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 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 151 nothing needs to be done with it afterwards. 152 153 In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in. 154 155 Level: intermediate 156 157 .seealso: 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(0); 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 . f - The face 194 - p - The point 195 196 Output Parameter: 197 . pgeom - The cell geometry of face f at point p 198 199 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 200 nothing needs to be done with it afterwards. 201 202 Level: intermediate 203 204 .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate() 205 @*/ 206 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom) 207 { 208 const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE; 209 const PetscInt dim = bd ? geom->dimEmbed : geom->dim; 210 const PetscInt dE = geom->dimEmbed; 211 const PetscInt Np = geom->numPoints; 212 213 PetscFunctionBeginHot; 214 pgeom->dim = dim; 215 pgeom->dimEmbed = dE; 216 //pgeom->isAffine = geom->isAffine; 217 if (geom->isAffine) { 218 if (!p) { 219 if (bd) { 220 pgeom->J = &geom->suppJ[0][c*Np*dE*dE]; 221 pgeom->invJ = &geom->suppInvJ[0][c*Np*dE*dE]; 222 pgeom->detJ = &geom->suppDetJ[0][c*Np]; 223 } else { 224 pgeom->J = &geom->J[c*Np*dE*dE]; 225 pgeom->invJ = &geom->invJ[c*Np*dE*dE]; 226 pgeom->detJ = &geom->detJ[c*Np]; 227 } 228 } 229 } else { 230 if (bd) { 231 pgeom->J = &geom->suppJ[0][(c*Np+p)*dE*dE]; 232 pgeom->invJ = &geom->suppInvJ[0][(c*Np+p)*dE*dE]; 233 pgeom->detJ = &geom->suppDetJ[0][c*Np+p]; 234 } else { 235 pgeom->J = &geom->J[(c*Np+p)*dE*dE]; 236 pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE]; 237 pgeom->detJ = &geom->detJ[c*Np+p]; 238 } 239 } 240 PetscFunctionReturn(0); 241 } 242 243 /*@ 244 PetscFEGeomComplete - Calculate derived quntites from base geometry specification 245 246 Input Parameter: 247 . geom - PetscFEGeom object 248 249 Level: intermediate 250 251 .seealso: PetscFEGeomCreate() 252 @*/ 253 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 254 { 255 PetscInt i, j, N, dE; 256 257 PetscFunctionBeginHot; 258 N = geom->numPoints * geom->numCells; 259 dE = geom->dimEmbed; 260 switch (dE) { 261 case 3: 262 for (i = 0; i < N; i++) { 263 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 264 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 265 } 266 break; 267 case 2: 268 for (i = 0; i < N; i++) { 269 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 270 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 271 } 272 break; 273 case 1: 274 for (i = 0; i < N; i++) { 275 geom->detJ[i] = PetscAbsReal(geom->J[i]); 276 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 277 } 278 break; 279 } 280 if (geom->n) { 281 for (i = 0; i < N; i++) { 282 for (j = 0; j < dE; j++) { 283 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 284 } 285 } 286 } 287 PetscFunctionReturn(0); 288 } 289