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 ierr = PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);CHKERRQ(ierr); 28 ierr = PetscNew(&g);CHKERRQ(ierr); 29 g->xi = p; 30 g->numCells = numCells; 31 g->numPoints = Nq; 32 g->dim = dim; 33 g->dimEmbed = dimEmbed; 34 g->isHybrid = PETSC_FALSE; 35 N = numCells * Nq; 36 ierr = PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);CHKERRQ(ierr); 37 if (faceData) { 38 ierr = PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);CHKERRQ(ierr); 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 ierr = PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);CHKERRQ(ierr); 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 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 if (!*geom) PetscFunctionReturn(0); 64 ierr = PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);CHKERRQ(ierr); 65 ierr = PetscFree((*geom)->invJ);CHKERRQ(ierr); 66 ierr = PetscFree2((*geom)->face,(*geom)->n);CHKERRQ(ierr); 67 ierr = PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);CHKERRQ(ierr); 68 ierr = PetscFree(*geom);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 /*@C 73 PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom 74 75 Input Parameters: 76 + geom - PetscFEGeom object 77 . cStart - The first cell in the chunk 78 - cEnd - The first cell not in the chunk 79 80 Output Parameter: 81 . chunkGeom - The chunk of cells 82 83 Level: intermediate 84 85 .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate() 86 @*/ 87 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 88 { 89 PetscInt Nq; 90 PetscInt dE; 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 PetscValidPointer(geom,1); 95 PetscValidPointer(chunkGeom,4); 96 if (!(*chunkGeom)) { 97 ierr = PetscNew(chunkGeom);CHKERRQ(ierr); 98 } 99 Nq = geom->numPoints; 100 dE= geom->dimEmbed; 101 (*chunkGeom)->dim = geom->dim; 102 (*chunkGeom)->dimEmbed = geom->dimEmbed; 103 (*chunkGeom)->numPoints = geom->numPoints; 104 (*chunkGeom)->numCells = cEnd - cStart; 105 (*chunkGeom)->xi = geom->xi; 106 (*chunkGeom)->v = &geom->v[Nq*dE*cStart]; 107 (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart]; 108 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL; 109 (*chunkGeom)->detJ = &geom->detJ[Nq*cStart]; 110 (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL; 111 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 112 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL; 113 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL; 114 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL; 115 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL; 116 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL; 117 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL; 118 (*chunkGeom)->isAffine = geom->isAffine; 119 PetscFunctionReturn(0); 120 } 121 122 /*@C 123 PetscFEGeomRestoreChunk - Restore the chunk 124 125 Input Parameters: 126 + geom - PetscFEGeom object 127 . cStart - The first cell in the chunk 128 . cEnd - The first cell not in the chunk 129 - chunkGeom - The chunk of cells 130 131 Level: intermediate 132 133 .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate() 134 @*/ 135 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 136 { 137 PetscErrorCode ierr; 138 139 PetscFunctionBegin; 140 ierr = PetscFree(*chunkGeom);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 /*@C 145 PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom 146 147 Input Parameters: 148 + geom - PetscFEGeom object 149 . c - The cell 150 . p - The point 151 - pcoords - The reference coordinates of point p, or NULL 152 153 Output Parameter: 154 . pgeom - The geometry of cell c at point p 155 156 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 157 nothing needs to be done with it afterwards. 158 159 In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in. 160 161 Level: intermediate 162 163 .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate() 164 @*/ 165 PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom) 166 { 167 const PetscInt dim = geom->dim; 168 const PetscInt dE = geom->dimEmbed; 169 const PetscInt Np = geom->numPoints; 170 171 PetscFunctionBeginHot; 172 pgeom->dim = dim; 173 pgeom->dimEmbed = dE; 174 //pgeom->isAffine = geom->isAffine; 175 if (geom->isAffine) { 176 if (!p) { 177 pgeom->xi = geom->xi; 178 pgeom->J = &geom->J[c*Np*dE*dE]; 179 pgeom->invJ = &geom->invJ[c*Np*dE*dE]; 180 pgeom->detJ = &geom->detJ[c*Np]; 181 pgeom->n = geom->n ? &geom->n[c*Np*dE] : NULL; 182 } 183 if (pcoords) {CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c*Np*dE], pgeom->J, pcoords, pgeom->v);} 184 } else { 185 pgeom->v = &geom->v[(c*Np+p)*dE]; 186 pgeom->J = &geom->J[(c*Np+p)*dE*dE]; 187 pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE]; 188 pgeom->detJ = &geom->detJ[c*Np+p]; 189 pgeom->n = geom->n ? &geom->n[(c*Np+p)*dE] : NULL; 190 } 191 PetscFunctionReturn(0); 192 } 193 194 /*@C 195 PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom 196 197 Input Parameters: 198 + geom - PetscFEGeom object 199 . f - The face 200 - p - The point 201 202 Output Parameter: 203 . pgeom - The cell geometry of face f at point p 204 205 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 206 nothing needs to be done with it afterwards. 207 208 Level: intermediate 209 210 .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate() 211 @*/ 212 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom) 213 { 214 const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isHybrid ? 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->dim = dim; 221 pgeom->dimEmbed = dE; 222 //pgeom->isAffine = geom->isAffine; 223 if (geom->isAffine) { 224 if (!p) { 225 if (bd) { 226 pgeom->J = &geom->suppJ[0][c*Np*dE*dE]; 227 pgeom->invJ = &geom->suppInvJ[0][c*Np*dE*dE]; 228 pgeom->detJ = &geom->suppDetJ[0][c*Np]; 229 } else { 230 pgeom->J = &geom->J[c*Np*dE*dE]; 231 pgeom->invJ = &geom->invJ[c*Np*dE*dE]; 232 pgeom->detJ = &geom->detJ[c*Np]; 233 } 234 } 235 } else { 236 if (bd) { 237 pgeom->J = &geom->suppJ[0][(c*Np+p)*dE*dE]; 238 pgeom->invJ = &geom->suppInvJ[0][(c*Np+p)*dE*dE]; 239 pgeom->detJ = &geom->suppDetJ[0][c*Np+p]; 240 } else { 241 pgeom->J = &geom->J[(c*Np+p)*dE*dE]; 242 pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE]; 243 pgeom->detJ = &geom->detJ[c*Np+p]; 244 } 245 } 246 PetscFunctionReturn(0); 247 } 248 249 /*@ 250 PetscFEGeomComplete - Calculate derived quntites from base geometry specification 251 252 Input Parameter: 253 . geom - PetscFEGeom object 254 255 Level: intermediate 256 257 .seealso: PetscFEGeomCreate() 258 @*/ 259 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 260 { 261 PetscInt i, j, N, dE; 262 263 PetscFunctionBeginHot; 264 N = geom->numPoints * geom->numCells; 265 dE = geom->dimEmbed; 266 switch (dE) { 267 case 3: 268 for (i = 0; i < N; i++) { 269 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 270 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 271 } 272 break; 273 case 2: 274 for (i = 0; i < N; i++) { 275 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 276 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 277 } 278 break; 279 case 1: 280 for (i = 0; i < N; i++) { 281 geom->detJ[i] = PetscAbsReal(geom->J[i]); 282 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 283 } 284 break; 285 } 286 if (geom->n) { 287 for (i = 0; i < N; i++) { 288 for (j = 0; j < dE; j++) { 289 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 290 } 291 } 292 } 293 PetscFunctionReturn(0); 294 } 295