1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petsc/private/petscimpl.h>
3
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;
10
PetscFEDestroy_Vector(PetscFE fe)11 static PetscErrorCode PetscFEDestroy_Vector(PetscFE fe)
12 {
13 PetscFE_Vec *v;
14
15 PetscFunctionBegin;
16 v = (PetscFE_Vec *)fe->data;
17 PetscCall(PetscFEDestroy(&v->scalar_fe));
18 PetscCall(PetscFree(v));
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
PetscFEView_Vector_Ascii(PetscFE fe,PetscViewer v)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;
30
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 }
58
PetscFEView_Vector(PetscFE fe,PetscViewer v)59 static PetscErrorCode PetscFEView_Vector(PetscFE fe, PetscViewer v)
60 {
61 PetscBool isascii;
62
63 PetscFunctionBegin;
64 PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
65 if (isascii) PetscCall(PetscFEView_Vector_Ascii(fe, v));
66 PetscFunctionReturn(PETSC_SUCCESS);
67 }
68
PetscFESetUp_Vector(PetscFE fe)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;
78
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;
98
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));
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
PetscFEGetDimension_Vector(PetscFE fe,PetscInt * dim)106 static PetscErrorCode PetscFEGetDimension_Vector(PetscFE fe, PetscInt *dim)
107 {
108 PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
109
110 PetscFunctionBegin;
111 PetscCall(PetscFEGetDimension(v->scalar_fe, dim));
112 *dim *= v->num_copies;
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
PetscFEVectorInsertTabulation(PetscFE fe,PetscInt npoints,const PetscReal points[],PetscInt k,PetscInt scalar_Nb,PetscInt scalar_point_stride,const PetscReal scalar_Tk[],PetscReal Tk[])116 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[])
117 {
118 PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
119 PetscInt scalar_Nc;
120 PetscInt Nc;
121 PetscInt cdim;
122 PetscInt dblock;
123 PetscInt block;
124 PetscInt scalar_block;
125
126 PetscFunctionBegin;
127 PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
128 Nc = scalar_Nc * v->num_copies;
129 PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
130 dblock = PetscPowInt(cdim, k);
131 block = Nc * dblock;
132 scalar_block = scalar_Nc * dblock;
133 for (PetscInt p = 0; p < npoints; p++) {
134 const PetscReal *scalar_Tp = &scalar_Tk[p * scalar_point_stride];
135 PetscReal *Tp = &Tk[p * scalar_point_stride * v->num_copies * v->num_copies];
136
137 for (PetscInt j = 0; j < v->num_copies; j++) {
138 for (PetscInt scalar_i = 0; scalar_i < scalar_Nb; scalar_i++) {
139 PetscInt i = v->interleave_basis ? (scalar_i * v->num_copies + j) : (j * scalar_Nb + scalar_i);
140 const PetscReal *scalar_Ti = &scalar_Tp[scalar_i * scalar_block];
141 PetscReal *Ti = &Tp[i * block];
142
143 for (PetscInt scalar_c = 0; scalar_c < scalar_Nc; scalar_c++) {
144 PetscInt c = v->interleave_components ? (scalar_c * v->num_copies + j) : (j * scalar_Nc + scalar_c);
145 const PetscReal *scalar_Tc = &scalar_Ti[scalar_c * dblock];
146 PetscReal *Tc = &Ti[c * dblock];
147
148 PetscCall(PetscArraycpy(Tc, scalar_Tc, dblock));
149 }
150 }
151 }
152 }
153 PetscFunctionReturn(PETSC_SUCCESS);
154 }
155
PetscFEComputeTabulation_Vector(PetscFE fe,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation T)156 static PetscErrorCode PetscFEComputeTabulation_Vector(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
157 {
158 PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
159 PetscInt scalar_Nc;
160 PetscInt scalar_Nb;
161 PetscInt cdim;
162 PetscTabulation scalar_T;
163
164 PetscFunctionBegin;
165 PetscAssert(npoints == T->Nr * T->Np, PetscObjectComm((PetscObject)fe), PETSC_ERR_PLIB, "Expected to be able decode PetscFECreateTabulation() from PetscTabulation fields");
166 PetscCall(PetscFECreateTabulation(v->scalar_fe, T->Nr, T->Np, points, K, &scalar_T));
167 PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
168 PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_Nb));
169 PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
170 for (PetscInt k = 0; k <= T->K; k++) {
171 PetscReal *Tk = T->T[k];
172 const PetscReal *scalar_Tk = scalar_T->T[k];
173 PetscInt dblock = PetscPowInt(cdim, k);
174 PetscInt scalar_block = scalar_Nc * dblock;
175 PetscInt scalar_point_stride = scalar_Nb * scalar_block;
176
177 if (v->interleave_basis) {
178 PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, scalar_Nb, scalar_point_stride, scalar_Tk, Tk));
179 } else {
180 PetscDualSpace scalar_dsp;
181 PetscSection ref_section;
182 PetscInt pStart, pEnd;
183
184 PetscCall(PetscFEGetDualSpace(v->scalar_fe, &scalar_dsp));
185 PetscCall(PetscDualSpaceGetSection(scalar_dsp, &ref_section));
186 PetscCall(PetscSectionGetChart(ref_section, &pStart, &pEnd));
187 for (PetscInt p = pStart; p < pEnd; p++) {
188 PetscInt dof, off;
189 PetscReal *Tp;
190 const PetscReal *scalar_Tp;
191
192 PetscCall(PetscSectionGetDof(ref_section, p, &dof));
193 PetscCall(PetscSectionGetOffset(ref_section, p, &off));
194 scalar_Tp = &scalar_Tk[off * scalar_block];
195 Tp = &Tk[off * scalar_block * v->num_copies * v->num_copies];
196 PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, dof, scalar_point_stride, scalar_Tp, Tp));
197 }
198 }
199 }
200 PetscCall(PetscTabulationDestroy(&scalar_T));
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
PetscFECreatePointTrace_Vector(PetscFE fe,PetscInt refPoint,PetscFE * trFE)204 static PetscErrorCode PetscFECreatePointTrace_Vector(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
205 {
206 PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
207 PetscFE scalar_trFE;
208 const char *name;
209
210 PetscFunctionBegin;
211 PetscCall(PetscFECreatePointTrace(v->scalar_fe, refPoint, &scalar_trFE));
212 PetscCall(PetscFECreateVector(scalar_trFE, v->num_copies, v->interleave_basis, v->interleave_components, trFE));
213 PetscCall(PetscFEDestroy(&scalar_trFE));
214 PetscCall(PetscObjectGetName((PetscObject)fe, &name));
215 if (name) PetscCall(PetscFESetName(*trFE, name));
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
219 PETSC_INTERN PetscErrorCode PetscFEIntegrate_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
220 PETSC_INTERN PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS, PetscInt, PetscBdPointFn *, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
221 PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS, PetscDS, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
222 PETSC_INTERN PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS, PetscWeakForm, PetscFEJacobianType, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
223 PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS, PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
224
PetscFEInitialize_Vector(PetscFE fe)225 static PetscErrorCode PetscFEInitialize_Vector(PetscFE fe)
226 {
227 PetscFunctionBegin;
228 fe->ops->setfromoptions = NULL;
229 fe->ops->setup = PetscFESetUp_Vector;
230 fe->ops->view = PetscFEView_Vector;
231 fe->ops->destroy = PetscFEDestroy_Vector;
232 fe->ops->getdimension = PetscFEGetDimension_Vector;
233 fe->ops->createpointtrace = PetscFECreatePointTrace_Vector;
234 fe->ops->computetabulation = PetscFEComputeTabulation_Vector;
235 fe->ops->integrate = PetscFEIntegrate_Basic;
236 fe->ops->integratebd = PetscFEIntegrateBd_Basic;
237 fe->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
238 fe->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
239 fe->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
240 fe->ops->integratejacobianaction = NULL;
241 fe->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
242 fe->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
243 fe->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
244 PetscFunctionReturn(PETSC_SUCCESS);
245 }
246
247 /*MC
248 PETSCFEVECTOR = "vector" - A vector-valued `PetscFE` object that is repeated copies
249 of the same underlying finite element.
250
251 Level: intermediate
252
253 .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PetscFECreateVector()`
254 M*/
PetscFECreate_Vector(PetscFE fe)255 PETSC_EXTERN PetscErrorCode PetscFECreate_Vector(PetscFE fe)
256 {
257 PetscFE_Vec *v;
258
259 PetscFunctionBegin;
260 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1);
261 PetscCall(PetscNew(&v));
262 fe->data = v;
263
264 PetscCall(PetscFEInitialize_Vector(fe));
265 PetscFunctionReturn(PETSC_SUCCESS);
266 }
267
268 /*@
269 PetscFECreateVector - Create a vector-valued `PetscFE` from multiple copies of an underlying
270 `PetscFE`.
271
272 Collective
273
274 Input Parameters:
275 + scalar_fe - a `PetscFE` finite element
276 . num_copies - a positive integer
277 . interleave_basis - if `PETSC_TRUE`, the first `num_copies` basis vectors
278 of the output finite element will be copies of the
279 first basis vector of `scalar_fe`, and so on for the
280 other basis vectors; otherwise all of the first-copy
281 basis vectors will come first, followed by all of the
282 second-copy, and so on.
283 - interleave_components - if `PETSC_TRUE`, the first `num_copies` components
284 of the output finite element will be copies of the
285 first component of `scalar_fe`, and so on for the
286 other components; otherwise all of the first-copy
287 components will come first, followed by all of the
288 second-copy, and so on.
289
290 Output Parameter:
291 . vector_fe - a `PetscFE` of type `PETSCFEVECTOR` that represent a discretization space with `num_copies` copies of `scalar_fe`
292
293 Level: intermediate
294
295 .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PETSCFEVECTOR`
296 @*/
PetscFECreateVector(PetscFE scalar_fe,PetscInt num_copies,PetscBool interleave_basis,PetscBool interleave_components,PetscFE * vector_fe)297 PetscErrorCode PetscFECreateVector(PetscFE scalar_fe, PetscInt num_copies, PetscBool interleave_basis, PetscBool interleave_components, PetscFE *vector_fe)
298 {
299 MPI_Comm comm;
300 PetscFE fe_vec;
301 PetscFE_Vec *v;
302 PetscInt scalar_Nc;
303 PetscQuadrature quad;
304
305 PetscFunctionBegin;
306 PetscValidHeaderSpecific(scalar_fe, PETSCFE_CLASSID, 1);
307 PetscCall(PetscObjectGetComm((PetscObject)scalar_fe, &comm));
308 PetscCheck(num_copies >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Expected positive number of copies, got %" PetscInt_FMT, num_copies);
309 PetscCall(PetscFECreate(comm, vector_fe));
310 fe_vec = *vector_fe;
311 PetscCall(PetscFESetType(fe_vec, PETSCFEVECTOR));
312 v = (PetscFE_Vec *)fe_vec->data;
313 PetscCall(PetscObjectReference((PetscObject)scalar_fe));
314 v->scalar_fe = scalar_fe;
315 v->num_copies = num_copies;
316 v->interleave_basis = interleave_basis;
317 v->interleave_components = interleave_components;
318 PetscCall(PetscFEGetNumComponents(scalar_fe, &scalar_Nc));
319 PetscCall(PetscFESetNumComponents(fe_vec, scalar_Nc * num_copies));
320 PetscCall(PetscFEGetQuadrature(scalar_fe, &quad));
321 PetscCall(PetscFESetQuadrature(fe_vec, quad));
322 PetscCall(PetscFEGetFaceQuadrature(scalar_fe, &quad));
323 PetscCall(PetscFESetFaceQuadrature(fe_vec, quad));
324 {
325 PetscSpace scalar_sp;
326 PetscSpace *copies;
327 PetscSpace sp;
328
329 PetscCall(PetscFEGetBasisSpace(scalar_fe, &scalar_sp));
330 PetscCall(PetscMalloc1(num_copies, &copies));
331 for (PetscInt i = 0; i < num_copies; i++) {
332 PetscCall(PetscObjectReference((PetscObject)scalar_sp));
333 copies[i] = scalar_sp;
334 }
335 PetscCall(PetscSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
336 PetscCall(PetscSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
337 PetscCall(PetscSpaceSetUp(sp));
338 PetscCall(PetscFESetBasisSpace(fe_vec, sp));
339 PetscCall(PetscSpaceDestroy(&sp));
340 for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscSpaceDestroy(&copies[i]));
341 PetscCall(PetscFree(copies));
342 }
343 {
344 PetscDualSpace scalar_sp;
345 PetscDualSpace *copies;
346 PetscDualSpace sp;
347
348 PetscCall(PetscFEGetDualSpace(scalar_fe, &scalar_sp));
349 PetscCall(PetscMalloc1(num_copies, &copies));
350 for (PetscInt i = 0; i < num_copies; i++) {
351 PetscCall(PetscObjectReference((PetscObject)scalar_sp));
352 copies[i] = scalar_sp;
353 }
354 PetscCall(PetscDualSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
355 PetscCall(PetscDualSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
356 PetscCall(PetscDualSpaceSetUp(sp));
357 PetscCall(PetscFESetDualSpace(fe_vec, sp));
358 PetscCall(PetscDualSpaceDestroy(&sp));
359 for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscDualSpaceDestroy(&copies[i]));
360 PetscCall(PetscFree(copies));
361 }
362 PetscFunctionReturn(PETSC_SUCCESS);
363 }
364