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 } else { 23 ierr = PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);CHKERRQ(ierr); 24 } 25 *geom = g; 26 PetscFunctionReturn(0); 27 } 28 29 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 30 { 31 PetscErrorCode ierr; 32 33 PetscFunctionBegin; 34 if (!*geom) PetscFunctionReturn(0); 35 ierr = PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);CHKERRQ(ierr); 36 ierr = PetscFree((*geom)->invJ);CHKERRQ(ierr); 37 ierr = PetscFree4((*geom)->face,(*geom)->n,(*geom)->suppInvJ[0],(*geom)->suppInvJ[1]);CHKERRQ(ierr); 38 ierr = PetscFree(*geom);CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 42 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 43 { 44 PetscInt Nq; 45 PetscInt dE; 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 PetscValidPointer(geom,1); 50 PetscValidPointer(chunkGeom,2); 51 if (!(*chunkGeom)) { 52 ierr = PetscNew(chunkGeom);CHKERRQ(ierr); 53 } 54 Nq = geom->numPoints; 55 dE= geom->dimEmbed; 56 (*chunkGeom)->dim = geom->dim; 57 (*chunkGeom)->dimEmbed = geom->dimEmbed; 58 (*chunkGeom)->numPoints = geom->numPoints; 59 (*chunkGeom)->numCells = cEnd - cStart; 60 (*chunkGeom)->xi = geom->xi; 61 (*chunkGeom)->v = &geom->v[Nq*dE*cStart]; 62 (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart]; 63 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL; 64 (*chunkGeom)->detJ = &geom->detJ[Nq*cStart]; 65 (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL; 66 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 67 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL; 68 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL; 69 (*chunkGeom)->isAffine = geom->isAffine; 70 PetscFunctionReturn(0); 71 } 72 73 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 74 { 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 ierr = PetscFree(*chunkGeom);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 82 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 83 { 84 PetscInt i, j, N, dE; 85 86 PetscFunctionBeginHot; 87 N = geom->numPoints * geom->numCells; 88 dE = geom->dimEmbed; 89 switch (dE) { 90 case 3: 91 for (i = 0; i < N; i++) { 92 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 93 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 94 } 95 break; 96 case 2: 97 for (i = 0; i < N; i++) { 98 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 99 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 100 } 101 break; 102 case 1: 103 for (i = 0; i < N; i++) { 104 geom->detJ[i] = PetscAbsReal(geom->J[i]); 105 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 106 } 107 break; 108 } 109 if (geom->n) { 110 for (i = 0; i < N; i++) { 111 for (j = 0; j < dE; j++) { 112 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 113 } 114 } 115 } 116 PetscFunctionReturn(0); 117 } 118 119