1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom) 4 { 5 PetscFEGeom *g; 6 PetscInt dim, Nq, N; 7 const PetscReal *p; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 ierr = PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);CHKERRQ(ierr); 12 ierr = PetscNew(&g);CHKERRQ(ierr); 13 g->xi = p; 14 g->numCells = numCells; 15 g->numPoints = Nq; 16 g->dim = dim; 17 g->dimEmbed = dimEmbed; 18 N = numCells * Nq; 19 ierr = PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);CHKERRQ(ierr); 20 if (faceData) { 21 ierr = PetscCalloc4(numCells, &g->face, N * dimEmbed, &g->n, N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]));CHKERRQ(ierr); 22 } 23 ierr = PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);CHKERRQ(ierr); 24 *geom = g; 25 PetscFunctionReturn(0); 26 } 27 28 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 29 { 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 if (!*geom) PetscFunctionReturn(0); 34 ierr = PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);CHKERRQ(ierr); 35 ierr = PetscFree((*geom)->invJ);CHKERRQ(ierr); 36 ierr = PetscFree4((*geom)->face,(*geom)->n,(*geom)->suppInvJ[0],(*geom)->suppInvJ[1]);CHKERRQ(ierr); 37 ierr = PetscFree(*geom);CHKERRQ(ierr); 38 PetscFunctionReturn(0); 39 } 40 41 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 42 { 43 PetscInt Nq; 44 PetscInt dE; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 PetscValidPointer(geom,1); 49 PetscValidPointer(chunkGeom,2); 50 if (!(*chunkGeom)) { 51 ierr = PetscNew(chunkGeom);CHKERRQ(ierr); 52 } 53 Nq = geom->numPoints; 54 dE= geom->dimEmbed; 55 (*chunkGeom)->dim = geom->dim; 56 (*chunkGeom)->dimEmbed = geom->dimEmbed; 57 (*chunkGeom)->numPoints = geom->numPoints; 58 (*chunkGeom)->numCells = cEnd - cStart; 59 (*chunkGeom)->xi = geom->xi; 60 (*chunkGeom)->v = &geom->v[Nq*dE*cStart]; 61 (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart]; 62 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL; 63 (*chunkGeom)->detJ = &geom->detJ[Nq*cStart]; 64 (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL; 65 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 66 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL; 67 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL; 68 (*chunkGeom)->isAffine = geom->isAffine; 69 PetscFunctionReturn(0); 70 } 71 72 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 73 { 74 PetscErrorCode ierr; 75 76 PetscFunctionBegin; 77 ierr = PetscFree(*chunkGeom);CHKERRQ(ierr); 78 PetscFunctionReturn(0); 79 } 80 81 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 82 { 83 PetscInt i, j, N, dE; 84 85 PetscFunctionBeginHot; 86 N = geom->numPoints * geom->numCells; 87 dE = geom->dimEmbed; 88 switch (dE) { 89 case 3: 90 for (i = 0; i < N; i++) { 91 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 92 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 93 } 94 break; 95 case 2: 96 for (i = 0; i < N; i++) { 97 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 98 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 99 } 100 break; 101 case 1: 102 for (i = 0; i < N; i++) { 103 geom->detJ[i] = PetscAbsReal(geom->J[i]); 104 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 105 } 106 break; 107 } 108 if (geom->n) { 109 for (i = 0; i < N; i++) { 110 for (j = 0; j < dE; j++) { 111 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 112 } 113 } 114 } 115 PetscFunctionReturn(0); 116 } 117 118