Actual source code: fecomposite.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petsc/private/dtimpl.h>
3: #include <petscblaslapack.h>
4: #include <petscdmplextransform.h>
6: static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
7: {
8: PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
10: PetscFree(cmp->embedding);
11: PetscFree(cmp);
12: return 0;
13: }
15: static PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
16: {
17: PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
18: DM K;
19: DMPolytopeType ct;
20: DMPlexTransform tr;
21: PetscReal *subpoint;
22: PetscBLASInt *pivots;
23: PetscBLASInt n, info;
24: PetscScalar *work, *invVscalar;
25: PetscInt dim, pdim, spdim, j, s;
26: PetscSection section;
28: /* Get affine mapping from reference cell to each subcell */
29: PetscDualSpaceGetDM(fem->dualSpace, &K);
30: DMGetDimension(K, &dim);
31: DMPlexGetCellType(K, 0, &ct);
32: DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
33: DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
34: DMPlexRefineRegularGetAffineTransforms(tr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
35: DMPlexTransformDestroy(&tr);
36: /* Determine dof embedding into subelements */
37: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
38: PetscSpaceGetDimension(fem->basisSpace, &spdim);
39: PetscMalloc1(cmp->numSubelements * spdim, &cmp->embedding);
40: DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);
41: PetscDualSpaceGetSection(fem->dualSpace, §ion);
42: for (s = 0; s < cmp->numSubelements; ++s) {
43: PetscInt sd = 0;
44: PetscInt closureSize;
45: PetscInt *closure = NULL;
47: DMPlexGetTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
48: for (j = 0; j < closureSize; j++) {
49: PetscInt point = closure[2 * j];
50: PetscInt dof, off, k;
52: PetscSectionGetDof(section, point, &dof);
53: PetscSectionGetOffset(section, point, &off);
54: for (k = 0; k < dof; k++) cmp->embedding[s * spdim + sd++] = off + k;
55: }
56: DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
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) invVscalar[(s * spdim + j) * spdim + k] += Bf[q * spdim + k] * weights[q];
83: }
84: PetscFree(Bf);
85: }
86: n = spdim;
87: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s * spdim * spdim], &n, pivots, &info));
88: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s * spdim * spdim], &n, pivots, work, &n, &info));
89: }
90: #if defined(PETSC_USE_COMPLEX)
91: for (s = 0; s < cmp->numSubelements * spdim * spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
92: PetscFree(invVscalar);
93: #endif
94: PetscFree2(pivots, work);
95: return 0;
96: }
98: static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
99: {
100: PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
101: DM dm;
102: DMPolytopeType ct;
103: PetscInt pdim; /* Dimension of FE space P */
104: PetscInt spdim; /* Dimension of subelement FE space P */
105: PetscInt dim; /* Spatial dimension */
106: PetscInt comp; /* Field components */
107: PetscInt *subpoints;
108: PetscReal *B = K >= 0 ? T->T[0] : NULL;
109: PetscReal *D = K >= 1 ? T->T[1] : NULL;
110: PetscReal *H = K >= 2 ? T->T[2] : NULL;
111: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint;
112: PetscInt p, s, d, e, j, k;
114: PetscDualSpaceGetDM(fem->dualSpace, &dm);
115: DMGetDimension(dm, &dim);
116: DMPlexGetCellType(dm, 0, &ct);
117: PetscSpaceGetDimension(fem->basisSpace, &spdim);
118: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
119: PetscFEGetNumComponents(fem, &comp);
120: /* Divide points into subelements */
121: DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);
122: DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);
123: for (p = 0; p < npoints; ++p) {
124: for (s = 0; s < cmp->numSubelements; ++s) {
125: PetscBool inside;
127: /* Apply transform, and check that point is inside cell */
128: for (d = 0; d < dim; ++d) {
129: subpoint[d] = -1.0;
130: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s * dim + d) * dim + e] * (points[p * dim + e] - cmp->v0[s * dim + e]);
131: }
132: DMPolytopeInCellTest(ct, subpoint, &inside);
133: if (inside) {
134: subpoints[p] = s;
135: break;
136: }
137: }
139: }
140: DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
141: /* Evaluate the prime basis functions at all points */
142: if (K >= 0) DMGetWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB);
143: if (K >= 1) DMGetWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD);
144: if (K >= 2) DMGetWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH);
145: PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
146: /* Translate to the nodal basis */
147: if (K >= 0) PetscArrayzero(B, npoints * pdim * comp);
148: if (K >= 1) PetscArrayzero(D, npoints * pdim * comp * dim);
149: if (K >= 2) PetscArrayzero(H, npoints * pdim * comp * dim * dim);
150: for (p = 0; p < npoints; ++p) {
151: const PetscInt s = subpoints[p];
153: if (B) {
154: /* Multiply by V^{-1} (spdim x spdim) */
155: for (j = 0; j < spdim; ++j) {
156: const PetscInt i = (p * pdim + cmp->embedding[s * spdim + j]) * comp;
158: B[i] = 0.0;
159: for (k = 0; k < spdim; ++k) B[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpB[p * spdim + k];
160: }
161: }
162: if (D) {
163: /* Multiply by V^{-1} (spdim x spdim) */
164: for (j = 0; j < spdim; ++j) {
165: for (d = 0; d < dim; ++d) {
166: const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim + d;
168: D[i] = 0.0;
169: for (k = 0; k < spdim; ++k) D[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpD[(p * spdim + k) * dim + d];
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) H[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpH[(p * spdim + k) * dim * dim + d];
181: }
182: }
183: }
184: }
185: DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
186: if (K >= 0) DMRestoreWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB);
187: if (K >= 1) DMRestoreWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD);
188: if (K >= 2) DMRestoreWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH);
189: return 0;
190: }
192: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
193: {
194: fem->ops->setfromoptions = NULL;
195: fem->ops->setup = PetscFESetUp_Composite;
196: fem->ops->view = NULL;
197: fem->ops->destroy = PetscFEDestroy_Composite;
198: fem->ops->getdimension = PetscFEGetDimension_Basic;
199: fem->ops->createtabulation = PetscFECreateTabulation_Composite;
200: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
201: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
202: fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
203: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
204: return 0;
205: }
207: /*MC
208: PETSCFECOMPOSITE = "composite" - A `PetscFEType` that represents a composite element
210: Level: intermediate
212: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
213: M*/
214: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
215: {
216: PetscFE_Composite *cmp;
219: PetscNew(&cmp);
220: fem->data = cmp;
222: cmp->numSubelements = -1;
223: cmp->v0 = NULL;
224: cmp->jac = NULL;
226: PetscFEInitialize_Composite(fem);
227: return 0;
228: }
230: /*@C
231: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
233: Not collective
235: Input Parameter:
236: . fem - The `PetscFE` object
238: Output Parameters:
239: + blockSize - The number of elements in a block
240: . numBlocks - The number of blocks in a batch
241: . batchSize - The number of elements in a batch
242: - numBatches - The number of batches in a chunk
244: Level: intermediate
246: .seealso: `PetscFE`, `PetscFECreate()`
247: @*/
248: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
249: {
250: PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
253: if (numSubelements) {
255: *numSubelements = cmp->numSubelements;
256: }
257: if (v0) {
259: *v0 = cmp->v0;
260: }
261: if (jac) {
263: *jac = cmp->jac;
264: }
265: if (invjac) {
267: *invjac = cmp->invjac;
268: }
269: return 0;
270: }