Actual source code: fecomposite.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/petscfeimpl.h>
2: #include <petsc/private/dtimpl.h>
3: #include <petsc/private/dmpleximpl.h>
4: #include <petscblaslapack.h>
6: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
7: {
8: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
9: PetscErrorCode ierr;
12: CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
13: PetscFree(cmp->embedding);
14: PetscFree(cmp);
15: return(0);
16: }
18: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
19: {
20: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
21: DM K;
22: PetscReal *subpoint;
23: PetscBLASInt *pivots;
24: PetscBLASInt n, info;
25: PetscScalar *work, *invVscalar;
26: PetscInt dim, pdim, spdim, j, s;
27: PetscErrorCode ierr;
30: /* Get affine mapping from reference cell to each subcell */
31: PetscDualSpaceGetDM(fem->dualSpace, &K);
32: DMGetDimension(K, &dim);
33: DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
34: CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
35: /* Determine dof embedding into subelements */
36: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
37: PetscSpaceGetDimension(fem->basisSpace, &spdim);
38: PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
39: DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);
40: for (s = 0; s < cmp->numSubelements; ++s) {
41: PetscInt sd = 0;
43: for (j = 0; j < pdim; ++j) {
44: PetscBool inside;
45: PetscQuadrature f;
46: PetscInt d, e;
48: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
49: /* Apply transform to first point, and check that point is inside subcell */
50: for (d = 0; d < dim; ++d) {
51: subpoint[d] = -1.0;
52: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
53: }
54: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
55: if (inside) {cmp->embedding[s*spdim+sd++] = j;}
56: }
57: if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
58: }
59: DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint);
60: /* Construct the change of basis from prime basis to nodal basis for each subelement */
61: PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
62: PetscMalloc2(spdim,&pivots,spdim,&work);
63: #if defined(PETSC_USE_COMPLEX)
64: PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);
65: #else
66: invVscalar = fem->invV;
67: #endif
68: for (s = 0; s < cmp->numSubelements; ++s) {
69: for (j = 0; j < spdim; ++j) {
70: PetscReal *Bf;
71: PetscQuadrature f;
72: const PetscReal *points, *weights;
73: PetscInt Nc, Nq, q, k;
75: PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
76: PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
77: PetscMalloc1(f->numPoints*spdim*Nc,&Bf);
78: PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
79: for (k = 0; k < spdim; ++k) {
80: /* n_j \cdot \phi_k */
81: invVscalar[(s*spdim + j)*spdim+k] = 0.0;
82: for (q = 0; q < Nq; ++q) {
83: invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*weights[q];
84: }
85: }
86: PetscFree(Bf);
87: }
88: n = spdim;
89: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
90: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info));
91: }
92: #if defined(PETSC_USE_COMPLEX)
93: for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
94: PetscFree(invVscalar);
95: #endif
96: PetscFree2(pivots,work);
97: return(0);
98: }
100: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
101: {
102: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
103: DM dm;
104: PetscInt pdim; /* Dimension of FE space P */
105: PetscInt spdim; /* Dimension of subelement FE space P */
106: PetscInt dim; /* Spatial dimension */
107: PetscInt comp; /* Field components */
108: PetscInt *subpoints;
109: PetscReal *tmpB, *tmpD, *tmpH, *subpoint;
110: PetscInt p, s, d, e, j, k;
111: PetscErrorCode ierr;
114: PetscDualSpaceGetDM(fem->dualSpace, &dm);
115: DMGetDimension(dm, &dim);
116: PetscSpaceGetDimension(fem->basisSpace, &spdim);
117: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
118: PetscFEGetNumComponents(fem, &comp);
119: /* Divide points into subelements */
120: DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);
121: DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);
122: for (p = 0; p < npoints; ++p) {
123: for (s = 0; s < cmp->numSubelements; ++s) {
124: PetscBool inside;
126: /* Apply transform, and check that point is inside cell */
127: for (d = 0; d < dim; ++d) {
128: subpoint[d] = -1.0;
129: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
130: }
131: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
132: if (inside) {subpoints[p] = s; break;}
133: }
134: if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
135: }
136: DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
137: /* Evaluate the prime basis functions at all points */
138: if (B) {DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
139: if (D) {DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
140: if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
141: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
142: /* Translate to the nodal basis */
143: if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
144: if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
145: if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
146: for (p = 0; p < npoints; ++p) {
147: const PetscInt s = subpoints[p];
149: if (B) {
150: /* Multiply by V^{-1} (spdim x spdim) */
151: for (j = 0; j < spdim; ++j) {
152: const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
154: B[i] = 0.0;
155: for (k = 0; k < spdim; ++k) {
156: B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
157: }
158: }
159: }
160: if (D) {
161: /* Multiply by V^{-1} (spdim x spdim) */
162: for (j = 0; j < spdim; ++j) {
163: for (d = 0; d < dim; ++d) {
164: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
166: D[i] = 0.0;
167: for (k = 0; k < spdim; ++k) {
168: D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
169: }
170: }
171: }
172: }
173: if (H) {
174: /* Multiply by V^{-1} (pdim x pdim) */
175: for (j = 0; j < spdim; ++j) {
176: for (d = 0; d < dim*dim; ++d) {
177: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
179: H[i] = 0.0;
180: for (k = 0; k < spdim; ++k) {
181: H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
182: }
183: }
184: }
185: }
186: }
187: DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
188: if (B) {DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
189: if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
190: if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
191: return(0);
192: }
194: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
195: {
197: fem->ops->setfromoptions = NULL;
198: fem->ops->setup = PetscFESetUp_Composite;
199: fem->ops->view = NULL;
200: fem->ops->destroy = PetscFEDestroy_Composite;
201: fem->ops->getdimension = PetscFEGetDimension_Basic;
202: fem->ops->gettabulation = PetscFEGetTabulation_Composite;
203: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
204: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
205: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
206: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
207: return(0);
208: }
210: /*MC
211: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
213: Level: intermediate
215: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
216: M*/
217: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
218: {
219: PetscFE_Composite *cmp;
220: PetscErrorCode ierr;
224: PetscNewLog(fem, &cmp);
225: fem->data = cmp;
227: cmp->cellRefiner = REFINER_NOOP;
228: cmp->numSubelements = -1;
229: cmp->v0 = NULL;
230: cmp->jac = NULL;
232: PetscFEInitialize_Composite(fem);
233: return(0);
234: }
236: /*@C
237: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
239: Not collective
241: Input Parameter:
242: . fem - The PetscFE object
244: Output Parameters:
245: + blockSize - The number of elements in a block
246: . numBlocks - The number of blocks in a batch
247: . batchSize - The number of elements in a batch
248: - numBatches - The number of batches in a chunk
250: Level: intermediate
252: .seealso: PetscFECreate()
253: @*/
254: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
255: {
256: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
264: return(0);
265: }