Actual source code: fevector.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petsc/private/petscimpl.h>
4: typedef struct _n_PetscFE_Vec {
5: PetscFE scalar_fe;
6: PetscInt num_copies;
7: PetscBool interleave_basis;
8: PetscBool interleave_components;
9: } PetscFE_Vec;
11: static PetscErrorCode PetscFEDestroy_Vector(PetscFE fe)
12: {
13: PetscFE_Vec *v;
15: PetscFunctionBegin;
16: v = (PetscFE_Vec *)fe->data;
17: PetscCall(PetscFEDestroy(&v->scalar_fe));
18: PetscCall(PetscFree(v));
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: static PetscErrorCode PetscFEView_Vector_Ascii(PetscFE fe, PetscViewer v)
23: {
24: PetscInt dim, Nc, scalar_Nc;
25: PetscSpace basis = NULL;
26: PetscDualSpace dual = NULL;
27: PetscQuadrature quad = NULL;
28: PetscFE_Vec *vec;
29: PetscViewerFormat fmt;
31: PetscFunctionBegin;
32: vec = (PetscFE_Vec *)fe->data;
33: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
34: PetscCall(PetscFEGetNumComponents(fe, &Nc));
35: PetscCall(PetscFEGetNumComponents(vec->scalar_fe, &scalar_Nc));
36: PetscCall(PetscFEGetBasisSpace(fe, &basis));
37: PetscCall(PetscFEGetDualSpace(fe, &dual));
38: PetscCall(PetscFEGetQuadrature(fe, &quad));
39: PetscCall(PetscViewerGetFormat(v, &fmt));
40: PetscCall(PetscViewerASCIIPushTab(v));
41: if (scalar_Nc == 1) {
42: PetscCall(PetscViewerASCIIPrintf(v, "Vector Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components\n", dim, Nc));
43: } else {
44: PetscCall(PetscViewerASCIIPrintf(v, "Vector Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components (%" PetscInt_FMT " copies of finite element with %" PetscInt_FMT " components)\n", dim, Nc, vec->num_copies, scalar_Nc));
45: }
46: if (fmt == PETSC_VIEWER_ASCII_INFO_DETAIL) {
47: PetscCall(PetscViewerASCIIPrintf(v, "Original finite element:\n"));
48: PetscCall(PetscViewerASCIIPushTab(v));
49: PetscCall(PetscFEView(vec->scalar_fe, v));
50: PetscCall(PetscViewerASCIIPopTab(v));
51: }
52: if (basis) PetscCall(PetscSpaceView(basis, v));
53: if (dual) PetscCall(PetscDualSpaceView(dual, v));
54: if (quad) PetscCall(PetscQuadratureView(quad, v));
55: PetscCall(PetscViewerASCIIPopTab(v));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode PetscFEView_Vector(PetscFE fe, PetscViewer v)
60: {
61: PetscBool iascii;
63: PetscFunctionBegin;
64: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
65: if (iascii) PetscCall(PetscFEView_Vector_Ascii(fe, v));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode PetscFESetUp_Vector(PetscFE fe)
70: {
71: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
72: PetscDualSpace dsp;
73: PetscInt n, Ncopies = v->num_copies;
74: PetscInt scalar_n;
75: PetscInt *d, *d_mapped;
76: PetscDualSpace_Sum *sum;
77: PetscBool is_sum;
79: PetscFunctionBegin;
80: PetscCall(PetscFESetUp(v->scalar_fe));
81: PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_n));
82: PetscCall(PetscFEGetDualSpace(fe, &dsp));
83: PetscCall(PetscObjectTypeCompare((PetscObject)dsp, PETSCDUALSPACESUM, &is_sum));
84: PetscCheck(is_sum, PetscObjectComm((PetscObject)fe), PETSC_ERR_ARG_INCOMP, "Expected PETSCDUALSPACESUM dual space");
85: sum = (PetscDualSpace_Sum *)dsp->data;
86: n = Ncopies * scalar_n;
87: PetscCall(PetscCalloc1(n * n, &fe->invV));
88: PetscCall(PetscMalloc2(scalar_n, &d, scalar_n, &d_mapped));
89: for (PetscInt i = 0; i < scalar_n; i++) d[i] = i;
90: for (PetscInt c = 0; c < Ncopies; c++) {
91: PetscCall(ISLocalToGlobalMappingApply(sum->all_rows[c], scalar_n, d, d_mapped));
92: for (PetscInt i = 0; i < scalar_n; i++) {
93: PetscInt iw = d_mapped[i];
94: PetscReal *row_w = &fe->invV[iw * n];
95: const PetscReal *row_r = &v->scalar_fe->invV[i * scalar_n];
96: PetscInt j0 = v->interleave_basis ? c : c * scalar_n;
97: PetscInt jstride = v->interleave_basis ? Ncopies : 1;
99: for (PetscInt j = 0; j < scalar_n; j++) row_w[j0 + j * jstride] = row_r[j];
100: }
101: }
102: PetscCall(PetscFree2(d, d_mapped));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static PetscErrorCode PetscFEGetDimension_Vector(PetscFE fe, PetscInt *dim)
108: {
109: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
111: PetscFunctionBegin;
112: PetscCall(PetscFEGetDimension(v->scalar_fe, dim));
113: *dim *= v->num_copies;
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode PetscFEVectorInsertTabulation(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt k, PetscInt scalar_Nb, PetscInt scalar_point_stride, const PetscReal scalar_Tk[], PetscReal Tk[])
118: {
119: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
120: PetscInt scalar_Nc;
121: PetscInt Nc;
122: PetscInt cdim;
123: PetscInt dblock;
124: PetscInt block;
125: PetscInt scalar_block;
127: PetscFunctionBegin;
128: PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
129: Nc = scalar_Nc * v->num_copies;
130: PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
131: dblock = PetscPowInt(cdim, k);
132: block = Nc * dblock;
133: scalar_block = scalar_Nc * dblock;
134: for (PetscInt p = 0; p < npoints; p++) {
135: const PetscReal *scalar_Tp = &scalar_Tk[p * scalar_point_stride];
136: PetscReal *Tp = &Tk[p * scalar_point_stride * v->num_copies * v->num_copies];
138: for (PetscInt j = 0; j < v->num_copies; j++) {
139: for (PetscInt scalar_i = 0; scalar_i < scalar_Nb; scalar_i++) {
140: PetscInt i = v->interleave_basis ? (scalar_i * v->num_copies + j) : (j * scalar_Nb + scalar_i);
141: const PetscReal *scalar_Ti = &scalar_Tp[scalar_i * scalar_block];
142: PetscReal *Ti = &Tp[i * block];
144: for (PetscInt scalar_c = 0; scalar_c < scalar_Nc; scalar_c++) {
145: PetscInt c = v->interleave_components ? (scalar_c * v->num_copies + j) : (j * scalar_Nc + scalar_c);
146: const PetscReal *scalar_Tc = &scalar_Ti[scalar_c * dblock];
147: PetscReal *Tc = &Ti[c * dblock];
149: PetscCall(PetscArraycpy(Tc, scalar_Tc, dblock));
150: }
151: }
152: }
153: }
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode PetscFECreateTabulation_Vector(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
158: {
159: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
160: PetscInt scalar_Nc;
161: PetscInt scalar_Nb;
162: PetscInt cdim;
163: PetscTabulation scalar_T;
165: PetscFunctionBegin;
166: PetscAssert(npoints == T->Nr * T->Np, PetscObjectComm((PetscObject)fe), PETSC_ERR_PLIB, "Expected to be able decode PetscFECreateTabulation() from PetscTabulation fields");
167: PetscCall(PetscFECreateTabulation(v->scalar_fe, T->Nr, T->Np, points, K, &scalar_T));
168: PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
169: PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_Nb));
170: PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
171: for (PetscInt k = 0; k <= T->K; k++) {
172: PetscReal *Tk = T->T[k];
173: const PetscReal *scalar_Tk = scalar_T->T[k];
174: PetscInt dblock = PetscPowInt(cdim, k);
175: PetscInt scalar_block = scalar_Nc * dblock;
176: PetscInt scalar_point_stride = scalar_Nb * scalar_block;
178: if (v->interleave_basis) {
179: PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, scalar_Nb, scalar_point_stride, scalar_Tk, Tk));
180: } else {
181: PetscDualSpace scalar_dsp;
182: PetscSection ref_section;
183: PetscInt pStart, pEnd;
185: PetscCall(PetscFEGetDualSpace(v->scalar_fe, &scalar_dsp));
186: PetscCall(PetscDualSpaceGetSection(scalar_dsp, &ref_section));
187: PetscCall(PetscSectionGetChart(ref_section, &pStart, &pEnd));
188: for (PetscInt p = pStart; p < pEnd; p++) {
189: PetscInt dof, off;
190: PetscReal *Tp;
191: const PetscReal *scalar_Tp;
193: PetscCall(PetscSectionGetDof(ref_section, p, &dof));
194: PetscCall(PetscSectionGetOffset(ref_section, p, &off));
195: scalar_Tp = &scalar_Tk[off * scalar_block];
196: Tp = &Tk[off * scalar_block * v->num_copies * v->num_copies];
197: PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, dof, scalar_point_stride, scalar_Tp, Tp));
198: }
199: }
200: }
201: PetscCall(PetscTabulationDestroy(&scalar_T));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode PetscFECreatePointTrace_Vector(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
206: {
207: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
208: PetscFE scalar_trFE;
209: const char *name;
211: PetscFunctionBegin;
212: PetscCall(PetscFECreatePointTrace(v->scalar_fe, refPoint, &scalar_trFE));
213: PetscCall(PetscFECreateVector(scalar_trFE, v->num_copies, v->interleave_basis, v->interleave_components, trFE));
214: PetscCall(PetscFEDestroy(&scalar_trFE));
215: PetscCall(PetscObjectGetName((PetscObject)fe, &name));
216: if (name) PetscCall(PetscFESetName(*trFE, name));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: PETSC_INTERN PetscErrorCode PetscFEIntegrate_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
221: PETSC_INTERN PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
222: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS, PetscDS, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
223: PETSC_INTERN PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS, PetscWeakForm, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
224: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS, PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
226: static PetscErrorCode PetscFEInitialize_Vector(PetscFE fe)
227: {
228: PetscFunctionBegin;
229: fe->ops->setfromoptions = NULL;
230: fe->ops->setup = PetscFESetUp_Vector;
231: fe->ops->view = PetscFEView_Vector;
232: fe->ops->destroy = PetscFEDestroy_Vector;
233: fe->ops->getdimension = PetscFEGetDimension_Vector;
234: fe->ops->createpointtrace = PetscFECreatePointTrace_Vector;
235: fe->ops->createtabulation = PetscFECreateTabulation_Vector;
236: fe->ops->integrate = PetscFEIntegrate_Basic;
237: fe->ops->integratebd = PetscFEIntegrateBd_Basic;
238: fe->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
239: fe->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
240: fe->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
241: fe->ops->integratejacobianaction = NULL;
242: fe->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
243: fe->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
244: fe->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*MC
249: PETSCFEVECTOR = "vector" - A vector-valued `PetscFE` object that is repeated copies
250: of the same underlying finite element.
252: Level: intermediate
254: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PetscFECreateVector()`
255: M*/
256: PETSC_EXTERN PetscErrorCode PetscFECreate_Vector(PetscFE fe)
257: {
258: PetscFE_Vec *v;
260: PetscFunctionBegin;
262: PetscCall(PetscNew(&v));
263: fe->data = v;
265: PetscCall(PetscFEInitialize_Vector(fe));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: PetscFECreateVector - Create a vector-valued `PetscFE` from multiple copies of an underlying
271: `PetscFE`.
273: Collective
275: Input Parameters:
276: + scalar_fe - a `PetscFE` finite element
277: . num_copies - a positive integer
278: . interleave_basis - if `PETSC_TRUE`, the first `num_copies` basis vectors
279: of the output finite element will be copies of the
280: first basis vector of `scalar_fe`, and so on for the
281: other basis vectors; otherwise all of the first-copy
282: basis vectors will come first, followed by all of the
283: second-copy, and so on.
284: - interleave_components - if `PETSC_TRUE`, the first `num_copies` components
285: of the output finite element will be copies of the
286: first component of `scalar_fe`, and so on for the
287: other components; otherwise all of the first-copy
288: components will come first, followed by all of the
289: second-copy, and so on.
291: Output Parameter:
292: . vector_fe - a `PetscFE` of type `PETSCFEVECTOR` that represent a discretization space with `num_copies` copies of `scalar_fe`
294: Level: intermediate
296: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PETSCFEVECTOR`
297: @*/
298: PetscErrorCode PetscFECreateVector(PetscFE scalar_fe, PetscInt num_copies, PetscBool interleave_basis, PetscBool interleave_components, PetscFE *vector_fe)
299: {
300: MPI_Comm comm;
301: PetscFE fe_vec;
302: PetscFE_Vec *v;
303: PetscInt scalar_Nc;
304: PetscQuadrature quad;
306: PetscFunctionBegin;
308: PetscCall(PetscObjectGetComm((PetscObject)scalar_fe, &comm));
309: PetscCheck(num_copies >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Expected positive number of copies, got %" PetscInt_FMT, num_copies);
310: PetscCall(PetscFECreate(comm, vector_fe));
311: fe_vec = *vector_fe;
312: PetscCall(PetscFESetType(fe_vec, PETSCFEVECTOR));
313: v = (PetscFE_Vec *)fe_vec->data;
314: PetscCall(PetscObjectReference((PetscObject)scalar_fe));
315: v->scalar_fe = scalar_fe;
316: v->num_copies = num_copies;
317: v->interleave_basis = interleave_basis;
318: v->interleave_components = interleave_components;
319: PetscCall(PetscFEGetNumComponents(scalar_fe, &scalar_Nc));
320: PetscCall(PetscFESetNumComponents(fe_vec, scalar_Nc * num_copies));
321: PetscCall(PetscFEGetQuadrature(scalar_fe, &quad));
322: PetscCall(PetscFESetQuadrature(fe_vec, quad));
323: PetscCall(PetscFEGetFaceQuadrature(scalar_fe, &quad));
324: PetscCall(PetscFESetFaceQuadrature(fe_vec, quad));
325: {
326: PetscSpace scalar_sp;
327: PetscSpace *copies;
328: PetscSpace sp;
330: PetscCall(PetscFEGetBasisSpace(scalar_fe, &scalar_sp));
331: PetscCall(PetscMalloc1(num_copies, &copies));
332: for (PetscInt i = 0; i < num_copies; i++) {
333: PetscCall(PetscObjectReference((PetscObject)scalar_sp));
334: copies[i] = scalar_sp;
335: }
336: PetscCall(PetscSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
337: PetscCall(PetscSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
338: PetscCall(PetscSpaceSetUp(sp));
339: PetscCall(PetscFESetBasisSpace(fe_vec, sp));
340: PetscCall(PetscSpaceDestroy(&sp));
341: for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscSpaceDestroy(&copies[i]));
342: PetscCall(PetscFree(copies));
343: }
344: {
345: PetscDualSpace scalar_sp;
346: PetscDualSpace *copies;
347: PetscDualSpace sp;
349: PetscCall(PetscFEGetDualSpace(scalar_fe, &scalar_sp));
350: PetscCall(PetscMalloc1(num_copies, &copies));
351: for (PetscInt i = 0; i < num_copies; i++) {
352: PetscCall(PetscObjectReference((PetscObject)scalar_sp));
353: copies[i] = scalar_sp;
354: }
355: PetscCall(PetscDualSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
356: PetscCall(PetscDualSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
357: PetscCall(PetscDualSpaceSetUp(sp));
358: PetscCall(PetscFESetDualSpace(fe_vec, sp));
359: PetscCall(PetscDualSpaceDestroy(&sp));
360: for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscDualSpaceDestroy(&copies[i]));
361: PetscCall(PetscFree(copies));
362: }
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }