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) {
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: static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
101: {
102: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
103: DM dm;
104: DMPolytopeType ct;
105: PetscInt pdim; /* Dimension of FE space P */
106: PetscInt spdim; /* Dimension of subelement FE space P */
107: PetscInt dim; /* Spatial dimension */
108: PetscInt comp; /* Field components */
109: PetscInt *subpoints;
110: PetscReal *B = K >= 0 ? T->T[0] : NULL;
111: PetscReal *D = K >= 1 ? T->T[1] : NULL;
112: PetscReal *H = K >= 2 ? T->T[2] : NULL;
113: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint;
114: PetscInt p, s, d, e, j, k;
116: PetscDualSpaceGetDM(fem->dualSpace, &dm);
117: DMGetDimension(dm, &dim);
118: DMPlexGetCellType(dm, 0, &ct);
119: PetscSpaceGetDimension(fem->basisSpace, &spdim);
120: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
121: PetscFEGetNumComponents(fem, &comp);
122: /* Divide points into subelements */
123: DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);
124: DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);
125: for (p = 0; p < npoints; ++p) {
126: for (s = 0; s < cmp->numSubelements; ++s) {
127: PetscBool inside;
129: /* Apply transform, and check that point is inside cell */
130: for (d = 0; d < dim; ++d) {
131: subpoint[d] = -1.0;
132: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
133: }
134: DMPolytopeInCellTest(ct, subpoint, &inside);
135: if (inside) {subpoints[p] = s; break;}
136: }
138: }
139: DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
140: /* Evaluate the prime basis functions at all points */
141: if (K >= 0) DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);
142: if (K >= 1) DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);
143: if (K >= 2) DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);
144: PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
145: /* Translate to the nodal basis */
146: if (K >= 0) PetscArrayzero(B, npoints*pdim*comp);
147: if (K >= 1) PetscArrayzero(D, npoints*pdim*comp*dim);
148: if (K >= 2) PetscArrayzero(H, npoints*pdim*comp*dim*dim);
149: for (p = 0; p < npoints; ++p) {
150: const PetscInt s = subpoints[p];
152: if (B) {
153: /* Multiply by V^{-1} (spdim x spdim) */
154: for (j = 0; j < spdim; ++j) {
155: const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
157: B[i] = 0.0;
158: for (k = 0; k < spdim; ++k) {
159: B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
160: }
161: }
162: }
163: if (D) {
164: /* Multiply by V^{-1} (spdim x spdim) */
165: for (j = 0; j < spdim; ++j) {
166: for (d = 0; d < dim; ++d) {
167: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
169: D[i] = 0.0;
170: for (k = 0; k < spdim; ++k) {
171: D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
172: }
173: }
174: }
175: }
176: if (H) {
177: /* Multiply by V^{-1} (pdim x pdim) */
178: for (j = 0; j < spdim; ++j) {
179: for (d = 0; d < dim*dim; ++d) {
180: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
182: H[i] = 0.0;
183: for (k = 0; k < spdim; ++k) {
184: H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
185: }
186: }
187: }
188: }
189: }
190: DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
191: if (K >= 0) DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);
192: if (K >= 1) DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);
193: if (K >= 2) DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);
194: return 0;
195: }
197: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
198: {
199: fem->ops->setfromoptions = NULL;
200: fem->ops->setup = PetscFESetUp_Composite;
201: fem->ops->view = NULL;
202: fem->ops->destroy = PetscFEDestroy_Composite;
203: fem->ops->getdimension = PetscFEGetDimension_Basic;
204: fem->ops->createtabulation = PetscFECreateTabulation_Composite;
205: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
206: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
207: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
208: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
209: return 0;
210: }
212: /*MC
213: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
215: Level: intermediate
217: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
218: M*/
219: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
220: {
221: PetscFE_Composite *cmp;
224: PetscNewLog(fem, &cmp);
225: fem->data = cmp;
227: cmp->numSubelements = -1;
228: cmp->v0 = NULL;
229: cmp->jac = NULL;
231: PetscFEInitialize_Composite(fem);
232: return 0;
233: }
235: /*@C
236: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
238: Not collective
240: Input Parameter:
241: . fem - The PetscFE object
243: Output Parameters:
244: + blockSize - The number of elements in a block
245: . numBlocks - The number of blocks in a batch
246: . batchSize - The number of elements in a batch
247: - numBatches - The number of batches in a chunk
249: Level: intermediate
251: .seealso: PetscFECreate()
252: @*/
253: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
254: {
255: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
262: return 0;
263: }