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 PetscErrorCode ierr; 25 26 PetscFunctionBegin; 27 CHKERRQ(PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL)); 28 CHKERRQ(PetscNew(&g)); 29 g->xi = p; 30 g->numCells = numCells; 31 g->numPoints = Nq; 32 g->dim = dim; 33 g->dimEmbed = dimEmbed; 34 g->isCohesive = PETSC_FALSE; 35 N = numCells * Nq; 36 CHKERRQ(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ)); 37 if (faceData) { 38 CHKERRQ(PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n)); 39 ierr = PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]), 40 N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]), 41 N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]));CHKERRQ(ierr); 42 } 43 CHKERRQ(PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ)); 44 *geom = g; 45 PetscFunctionReturn(0); 46 } 47 48 /*@C 49 PetscFEGeomDestroy - Destroy a PetscFEGeom object 50 51 Input Parameter: 52 . geom - PetscFEGeom object 53 54 Level: beginner 55 56 .seealso: PetscFEGeomCreate() 57 @*/ 58 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 59 { 60 PetscFunctionBegin; 61 if (!*geom) PetscFunctionReturn(0); 62 CHKERRQ(PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ)); 63 CHKERRQ(PetscFree((*geom)->invJ)); 64 CHKERRQ(PetscFree2((*geom)->face,(*geom)->n)); 65 CHKERRQ(PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1])); 66 CHKERRQ(PetscFree(*geom)); 67 PetscFunctionReturn(0); 68 } 69 70 /*@C 71 PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom 72 73 Input Parameters: 74 + geom - PetscFEGeom object 75 . cStart - The first cell in the chunk 76 - cEnd - The first cell not in the chunk 77 78 Output Parameter: 79 . chunkGeom - The chunk of cells 80 81 Level: intermediate 82 83 .seealso: 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 PetscValidPointer(geom,1); 92 PetscValidPointer(chunkGeom,4); 93 if (!(*chunkGeom)) { 94 CHKERRQ(PetscNew(chunkGeom)); 95 } 96 Nq = geom->numPoints; 97 dE= geom->dimEmbed; 98 (*chunkGeom)->dim = geom->dim; 99 (*chunkGeom)->dimEmbed = geom->dimEmbed; 100 (*chunkGeom)->numPoints = geom->numPoints; 101 (*chunkGeom)->numCells = cEnd - cStart; 102 (*chunkGeom)->xi = geom->xi; 103 (*chunkGeom)->v = &geom->v[Nq*dE*cStart]; 104 (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart]; 105 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL; 106 (*chunkGeom)->detJ = &geom->detJ[Nq*cStart]; 107 (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL; 108 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 109 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL; 110 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL; 111 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL; 112 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL; 113 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL; 114 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL; 115 (*chunkGeom)->isAffine = geom->isAffine; 116 PetscFunctionReturn(0); 117 } 118 119 /*@C 120 PetscFEGeomRestoreChunk - Restore the chunk 121 122 Input Parameters: 123 + geom - PetscFEGeom object 124 . cStart - The first cell in the chunk 125 . cEnd - The first cell not in the chunk 126 - chunkGeom - The chunk of cells 127 128 Level: intermediate 129 130 .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate() 131 @*/ 132 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 133 { 134 PetscFunctionBegin; 135 CHKERRQ(PetscFree(*chunkGeom)); 136 PetscFunctionReturn(0); 137 } 138 139 /*@C 140 PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom 141 142 Input Parameters: 143 + geom - PetscFEGeom object 144 . c - The cell 145 . p - The point 146 - pcoords - The reference coordinates of point p, or NULL 147 148 Output Parameter: 149 . pgeom - The geometry of cell c at point p 150 151 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 152 nothing needs to be done with it afterwards. 153 154 In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in. 155 156 Level: intermediate 157 158 .seealso: 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->dim = dim; 168 pgeom->dimEmbed = dE; 169 //pgeom->isAffine = geom->isAffine; 170 if (geom->isAffine) { 171 if (!p) { 172 pgeom->xi = geom->xi; 173 pgeom->J = &geom->J[c*Np*dE*dE]; 174 pgeom->invJ = &geom->invJ[c*Np*dE*dE]; 175 pgeom->detJ = &geom->detJ[c*Np]; 176 pgeom->n = geom->n ? &geom->n[c*Np*dE] : NULL; 177 } 178 if (pcoords) {CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c*Np*dE], pgeom->J, pcoords, pgeom->v);} 179 } else { 180 pgeom->v = &geom->v[(c*Np+p)*dE]; 181 pgeom->J = &geom->J[(c*Np+p)*dE*dE]; 182 pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE]; 183 pgeom->detJ = &geom->detJ[c*Np+p]; 184 pgeom->n = geom->n ? &geom->n[(c*Np+p)*dE] : NULL; 185 } 186 PetscFunctionReturn(0); 187 } 188 189 /*@C 190 PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom 191 192 Input Parameters: 193 + geom - PetscFEGeom object 194 . f - The face 195 - p - The point 196 197 Output Parameter: 198 . pgeom - The cell geometry of face f at point p 199 200 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 201 nothing needs to be done with it afterwards. 202 203 Level: intermediate 204 205 .seealso: 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(0); 242 } 243 244 /*@ 245 PetscFEGeomComplete - Calculate derived quntites from base geometry specification 246 247 Input Parameter: 248 . geom - PetscFEGeom object 249 250 Level: intermediate 251 252 .seealso: 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++) { 284 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 285 } 286 } 287 } 288 PetscFunctionReturn(0); 289 } 290