1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2
3 /*@
4 PetscDualSpaceSumGetNumSubspaces - Get the number of spaces in the sum space
5
6 Input Parameter:
7 . sp - the dual space object
8
9 Output Parameter:
10 . numSumSpaces - the number of spaces
11
12 Level: intermediate
13
14 Note:
15 The name NumSubspaces is slightly misleading because it is actually getting the number of defining spaces of the sum, not a number of Subspaces of it
16
17 .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetNumSubspaces()`
18 @*/
PetscDualSpaceSumGetNumSubspaces(PetscDualSpace sp,PetscInt * numSumSpaces)19 PetscErrorCode PetscDualSpaceSumGetNumSubspaces(PetscDualSpace sp, PetscInt *numSumSpaces)
20 {
21 PetscFunctionBegin;
22 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
23 PetscAssertPointer(numSumSpaces, 2);
24 PetscTryMethod(sp, "PetscDualSpaceSumGetNumSubspaces_C", (PetscDualSpace, PetscInt *), (sp, numSumSpaces));
25 PetscFunctionReturn(PETSC_SUCCESS);
26 }
27
28 /*@
29 PetscDualSpaceSumSetNumSubspaces - Set the number of spaces in the sum space
30
31 Input Parameters:
32 + sp - the dual space object
33 - numSumSpaces - the number of spaces
34
35 Level: intermediate
36
37 Note:
38 The name NumSubspaces is slightly misleading because it is actually setting the number of defining spaces of the sum, not a number of Subspaces of it
39
40 .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetNumSubspaces()`
41 @*/
PetscDualSpaceSumSetNumSubspaces(PetscDualSpace sp,PetscInt numSumSpaces)42 PetscErrorCode PetscDualSpaceSumSetNumSubspaces(PetscDualSpace sp, PetscInt numSumSpaces)
43 {
44 PetscFunctionBegin;
45 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
46 PetscTryMethod(sp, "PetscDualSpaceSumSetNumSubspaces_C", (PetscDualSpace, PetscInt), (sp, numSumSpaces));
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
50 /*@
51 PetscDualSpaceSumGetConcatenate - Get the concatenate flag for this space.
52
53 Input Parameter:
54 . sp - the dual space object
55
56 Output Parameter:
57 . concatenate - flag indicating whether subspaces are concatenated.
58
59 Level: intermediate
60
61 Note:
62 A concatenated sum space will have the number of components equal to the sum of the number of
63 components of all subspaces. A non-concatenated, or direct sum space will have the same
64 number of components as its subspaces.
65
66 .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetConcatenate()`
67 @*/
PetscDualSpaceSumGetConcatenate(PetscDualSpace sp,PetscBool * concatenate)68 PetscErrorCode PetscDualSpaceSumGetConcatenate(PetscDualSpace sp, PetscBool *concatenate)
69 {
70 PetscFunctionBegin;
71 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
72 PetscTryMethod(sp, "PetscDualSpaceSumGetConcatenate_C", (PetscDualSpace, PetscBool *), (sp, concatenate));
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
76 /*@
77 PetscDualSpaceSumSetConcatenate - Sets the concatenate flag for this space.
78
79 Input Parameters:
80 + sp - the dual space object
81 - concatenate - are subspaces concatenated components (true) or direct summands (false)
82
83 Level: intermediate
84
85 Notes:
86 A concatenated sum space will have the number of components equal to the sum of the number of
87 components of all subspaces. A non-concatenated, or direct sum space will have the same
88 number of components as its subspaces .
89
90 .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetConcatenate()`
91 @*/
PetscDualSpaceSumSetConcatenate(PetscDualSpace sp,PetscBool concatenate)92 PetscErrorCode PetscDualSpaceSumSetConcatenate(PetscDualSpace sp, PetscBool concatenate)
93 {
94 PetscFunctionBegin;
95 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
96 PetscTryMethod(sp, "PetscDualSpaceSumSetConcatenate_C", (PetscDualSpace, PetscBool), (sp, concatenate));
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99
100 /*@
101 PetscDualSpaceSumGetSubspace - Get a space in the sum space
102
103 Input Parameters:
104 + sp - the dual space object
105 - s - The space number
106
107 Output Parameter:
108 . subsp - the `PetscDualSpace`
109
110 Level: intermediate
111
112 Note:
113 The name GetSubspace is slightly misleading because it is actually getting one of the defining spaces of the sum, not a Subspace of it
114
115 .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetSubspace()`
116 @*/
PetscDualSpaceSumGetSubspace(PetscDualSpace sp,PetscInt s,PetscDualSpace * subsp)117 PetscErrorCode PetscDualSpaceSumGetSubspace(PetscDualSpace sp, PetscInt s, PetscDualSpace *subsp)
118 {
119 PetscFunctionBegin;
120 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
121 PetscAssertPointer(subsp, 3);
122 PetscTryMethod(sp, "PetscDualSpaceSumGetSubspace_C", (PetscDualSpace, PetscInt, PetscDualSpace *), (sp, s, subsp));
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
126 /*@
127 PetscDualSpaceSumSetSubspace - Set a space in the sum space
128
129 Input Parameters:
130 + sp - the dual space object
131 . s - The space number
132 - subsp - the number of spaces
133
134 Level: intermediate
135
136 Note:
137 The name SetSubspace is slightly misleading because it is actually setting one of the defining spaces of the sum, not a Subspace of it
138
139 .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetSubspace()`
140 @*/
PetscDualSpaceSumSetSubspace(PetscDualSpace sp,PetscInt s,PetscDualSpace subsp)141 PetscErrorCode PetscDualSpaceSumSetSubspace(PetscDualSpace sp, PetscInt s, PetscDualSpace subsp)
142 {
143 PetscFunctionBegin;
144 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
145 if (subsp) PetscValidHeaderSpecific(subsp, PETSCDUALSPACE_CLASSID, 3);
146 PetscTryMethod(sp, "PetscDualSpaceSumSetSubspace_C", (PetscDualSpace, PetscInt, PetscDualSpace), (sp, s, subsp));
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149
PetscDualSpaceSumGetNumSubspaces_Sum(PetscDualSpace space,PetscInt * numSumSpaces)150 static PetscErrorCode PetscDualSpaceSumGetNumSubspaces_Sum(PetscDualSpace space, PetscInt *numSumSpaces)
151 {
152 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
153
154 PetscFunctionBegin;
155 *numSumSpaces = sum->numSumSpaces;
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
PetscDualSpaceSumSetNumSubspaces_Sum(PetscDualSpace space,PetscInt numSumSpaces)159 static PetscErrorCode PetscDualSpaceSumSetNumSubspaces_Sum(PetscDualSpace space, PetscInt numSumSpaces)
160 {
161 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
162 PetscInt Ns = sum->numSumSpaces;
163
164 PetscFunctionBegin;
165 PetscCheck(!sum->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
166 if (numSumSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
167 for (PetscInt s = 0; s < Ns; ++s) PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[s]));
168 PetscCall(PetscFree(sum->sumspaces));
169
170 Ns = sum->numSumSpaces = numSumSpaces;
171 PetscCall(PetscCalloc1(Ns, &sum->sumspaces));
172 PetscFunctionReturn(PETSC_SUCCESS);
173 }
174
PetscDualSpaceSumGetConcatenate_Sum(PetscDualSpace sp,PetscBool * concatenate)175 static PetscErrorCode PetscDualSpaceSumGetConcatenate_Sum(PetscDualSpace sp, PetscBool *concatenate)
176 {
177 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
178
179 PetscFunctionBegin;
180 *concatenate = sum->concatenate;
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
PetscDualSpaceSumSetConcatenate_Sum(PetscDualSpace sp,PetscBool concatenate)184 static PetscErrorCode PetscDualSpaceSumSetConcatenate_Sum(PetscDualSpace sp, PetscBool concatenate)
185 {
186 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
187
188 PetscFunctionBegin;
189 PetscCheck(!sum->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change space concatenation after setup called.");
190
191 sum->concatenate = concatenate;
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
PetscDualSpaceSumGetSubspace_Sum(PetscDualSpace space,PetscInt s,PetscDualSpace * subspace)195 static PetscErrorCode PetscDualSpaceSumGetSubspace_Sum(PetscDualSpace space, PetscInt s, PetscDualSpace *subspace)
196 {
197 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
198 PetscInt Ns = sum->numSumSpaces;
199
200 PetscFunctionBegin;
201 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscDualSpaceSumSetNumSubspaces() first");
202 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
203
204 *subspace = sum->sumspaces[s];
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
PetscDualSpaceSumSetSubspace_Sum(PetscDualSpace space,PetscInt s,PetscDualSpace subspace)208 static PetscErrorCode PetscDualSpaceSumSetSubspace_Sum(PetscDualSpace space, PetscInt s, PetscDualSpace subspace)
209 {
210 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
211 PetscInt Ns = sum->numSumSpaces;
212
213 PetscFunctionBegin;
214 PetscCheck(!sum->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
215 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscDualSpaceSumSetNumSubspaces() first");
216 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
217
218 PetscCall(PetscObjectReference((PetscObject)subspace));
219 PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[s]));
220 sum->sumspaces[s] = subspace;
221 PetscFunctionReturn(PETSC_SUCCESS);
222 }
223
PetscDualSpaceDuplicate_Sum(PetscDualSpace sp,PetscDualSpace spNew)224 static PetscErrorCode PetscDualSpaceDuplicate_Sum(PetscDualSpace sp, PetscDualSpace spNew)
225 {
226 PetscInt num_subspaces, Nc;
227 PetscBool concatenate, interleave_basis, interleave_components;
228 PetscDualSpace subsp_first = NULL;
229 PetscDualSpace subsp_dup_first = NULL;
230 DM K;
231
232 PetscFunctionBegin;
233 PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
234 PetscCall(PetscDualSpaceSumSetNumSubspaces(spNew, num_subspaces));
235 PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
236 PetscCall(PetscDualSpaceSumSetConcatenate(spNew, concatenate));
237 PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
238 PetscCall(PetscDualSpaceSumSetInterleave(spNew, interleave_basis, interleave_components));
239 PetscCall(PetscDualSpaceGetDM(sp, &K));
240 PetscCall(PetscDualSpaceSetDM(spNew, K));
241 PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
242 PetscCall(PetscDualSpaceSetNumComponents(spNew, Nc));
243 for (PetscInt s = 0; s < num_subspaces; s++) {
244 PetscDualSpace subsp, subspNew;
245
246 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
247 if (s == 0) {
248 subsp_first = subsp;
249 PetscCall(PetscDualSpaceDuplicate(subsp, &subsp_dup_first));
250 subspNew = subsp_dup_first;
251 } else if (subsp == subsp_first) {
252 PetscCall(PetscObjectReference((PetscObject)subsp_dup_first));
253 subspNew = subsp_dup_first;
254 } else {
255 PetscCall(PetscDualSpaceDuplicate(subsp, &subspNew));
256 }
257 PetscCall(PetscDualSpaceSumSetSubspace(spNew, s, subspNew));
258 PetscCall(PetscDualSpaceDestroy(&subspNew));
259 }
260 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262
PetscDualSpaceSumCreateQuadrature(PetscInt Ns,PetscInt cdim,PetscBool uniform_points,PetscQuadrature subquads[],PetscQuadrature * fullquad)263 static PetscErrorCode PetscDualSpaceSumCreateQuadrature(PetscInt Ns, PetscInt cdim, PetscBool uniform_points, PetscQuadrature subquads[], PetscQuadrature *fullquad)
264 {
265 PetscReal *points;
266 PetscInt Npoints;
267
268 PetscFunctionBegin;
269 if (uniform_points) {
270 PetscCall(PetscObjectReference((PetscObject)subquads[0]));
271 *fullquad = subquads[0];
272 PetscFunctionReturn(PETSC_SUCCESS);
273 }
274 Npoints = 0;
275 for (PetscInt s = 0; s < Ns; s++) {
276 PetscInt subNpoints;
277
278 if (!subquads[s]) continue;
279 PetscCall(PetscQuadratureGetData(subquads[s], NULL, NULL, &subNpoints, NULL, NULL));
280 Npoints += subNpoints;
281 }
282 *fullquad = NULL;
283 if (!Npoints) PetscFunctionReturn(PETSC_SUCCESS);
284 PetscCall(PetscMalloc1(Npoints * cdim, &points));
285 for (PetscInt s = 0, offset = 0; s < Ns; s++) {
286 PetscInt subNpoints;
287 const PetscReal *subpoints;
288
289 PetscCall(PetscQuadratureGetData(subquads[s], NULL, NULL, &subNpoints, &subpoints, NULL));
290 PetscCall(PetscArraycpy(&points[offset], subpoints, cdim * subNpoints));
291 offset += cdim * subNpoints;
292 }
293 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, fullquad));
294 PetscCall(PetscQuadratureSetData(*fullquad, cdim, 0, Npoints, points, NULL));
295 PetscFunctionReturn(PETSC_SUCCESS);
296 }
297
PetscDualSpaceSumCreateMatrix(PetscDualSpace sp,Mat submats[],ISLocalToGlobalMapping map_rows[],ISLocalToGlobalMapping map_cols[],PetscQuadrature fullquad,Mat * fullmat)298 static PetscErrorCode PetscDualSpaceSumCreateMatrix(PetscDualSpace sp, Mat submats[], ISLocalToGlobalMapping map_rows[], ISLocalToGlobalMapping map_cols[], PetscQuadrature fullquad, Mat *fullmat)
299 {
300 Mat mat;
301 PetscInt *i = NULL, *j = NULL;
302 PetscScalar *v = NULL;
303 PetscInt npoints;
304 PetscInt nrows = 0, ncols, nnz = 0;
305 PetscInt Nc, num_subspaces;
306
307 PetscFunctionBegin;
308 PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
309 PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
310 PetscCall(PetscQuadratureGetData(fullquad, NULL, NULL, &npoints, NULL, NULL));
311 ncols = Nc * npoints;
312 for (PetscInt s = 0; s < num_subspaces; s++) {
313 // Get the COO data for each matrix, map the is and js, and append to growing COO data
314 PetscInt sNrows, sNcols;
315 Mat smat;
316 const PetscInt *si;
317 const PetscInt *sj;
318 PetscScalar *sv;
319 PetscMemType memtype;
320 PetscInt snz;
321 PetscInt snz_actual;
322 PetscInt *cooi;
323 PetscInt *cooj;
324 PetscScalar *coov;
325 PetscScalar *v_new;
326 PetscInt *i_new;
327 PetscInt *j_new;
328
329 if (!submats[s]) continue;
330 PetscCall(MatGetSize(submats[s], &sNrows, &sNcols));
331 nrows += sNrows;
332 PetscCall(MatConvert(submats[s], MATSEQAIJ, MAT_INITIAL_MATRIX, &smat));
333 PetscCall(MatBindToCPU(smat, PETSC_TRUE));
334 PetscCall(MatSeqAIJGetCSRAndMemType(smat, &si, &sj, &sv, &memtype));
335 PetscCheck(memtype == PETSC_MEMTYPE_HOST, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not convert matrix to host memory");
336 snz = si[sNrows];
337
338 PetscCall(PetscMalloc1(nnz + snz, &v_new));
339 PetscCall(PetscArraycpy(v_new, v, nnz));
340 PetscCall(PetscFree(v));
341 v = v_new;
342
343 PetscCall(PetscMalloc1(nnz + snz, &i_new));
344 PetscCall(PetscArraycpy(i_new, i, nnz));
345 PetscCall(PetscFree(i));
346 i = i_new;
347
348 PetscCall(PetscMalloc1(nnz + snz, &j_new));
349 PetscCall(PetscArraycpy(j_new, j, nnz));
350 PetscCall(PetscFree(j));
351 j = j_new;
352
353 PetscCall(PetscMalloc2(snz, &cooi, snz, &cooj));
354
355 coov = &v[nnz];
356
357 snz_actual = 0;
358 for (PetscInt row = 0; row < sNrows; row++) {
359 for (PetscInt k = si[row]; k < si[row + 1]; k++, snz_actual++) {
360 cooi[snz_actual] = row;
361 cooj[snz_actual] = sj[k];
362 coov[snz_actual] = sv[k];
363 }
364 }
365 PetscCall(MatDestroy(&smat));
366
367 if (snz_actual > 0) {
368 PetscCall(ISLocalToGlobalMappingApply(map_rows[s], snz_actual, cooi, &i[nnz]));
369 PetscCall(ISLocalToGlobalMappingApply(map_cols[s], snz_actual, cooj, &j[nnz]));
370 nnz += snz_actual;
371 }
372 PetscCall(PetscFree2(cooi, cooj));
373 }
374 PetscCall(MatCreate(PETSC_COMM_SELF, fullmat));
375 mat = *fullmat;
376 PetscCall(MatSetSizes(mat, nrows, ncols, nrows, ncols));
377 PetscCall(MatSetType(mat, MATSEQAIJ));
378 PetscCall(MatSetPreallocationCOO(mat, nnz, i, j));
379 PetscCall(MatSetValuesCOO(mat, v, INSERT_VALUES));
380 PetscCall(PetscFree(i));
381 PetscCall(PetscFree(j));
382 PetscCall(PetscFree(v));
383 PetscFunctionReturn(PETSC_SUCCESS);
384 }
385
PetscDualSpaceSumCreateMappings(PetscDualSpace sp,PetscBool interior,PetscBool uniform_points,ISLocalToGlobalMapping map_row[],ISLocalToGlobalMapping map_col[])386 static PetscErrorCode PetscDualSpaceSumCreateMappings(PetscDualSpace sp, PetscBool interior, PetscBool uniform_points, ISLocalToGlobalMapping map_row[], ISLocalToGlobalMapping map_col[])
387 {
388 PetscSection section;
389 PetscInt pStart, pEnd, Ns, Nc;
390 PetscInt num_points = 0;
391 PetscBool interleave_basis, interleave_components, concatenate;
392 PetscInt *roffset;
393
394 PetscFunctionBegin;
395 if (interior) {
396 PetscCall(PetscDualSpaceGetInteriorSection(sp, §ion));
397 } else {
398 PetscCall(PetscDualSpaceGetSection(sp, §ion));
399 }
400 PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
401 PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
402 PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
403 PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
404 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
405 for (PetscInt p = pStart; p < pEnd; p++) {
406 PetscInt dof;
407
408 PetscCall(PetscSectionGetDof(section, p, &dof));
409 num_points += (dof > 0);
410 }
411 PetscCall(PetscCalloc1(pEnd - pStart, &roffset));
412 for (PetscInt s = 0, coffset = 0; s < Ns; s++) {
413 PetscDualSpace subsp;
414 PetscSection subsection;
415 IS is_row, is_col;
416 PetscInt subNb;
417 PetscInt sNc, sNpoints, sNcols;
418 PetscQuadrature quad;
419
420 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
421 PetscCall(PetscDualSpaceGetNumComponents(subsp, &sNc));
422 if (interior) {
423 PetscCall(PetscDualSpaceGetInteriorSection(subsp, &subsection));
424 PetscCall(PetscDualSpaceGetInteriorData(subsp, &quad, NULL));
425 } else {
426 PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
427 PetscCall(PetscDualSpaceGetAllData(subsp, &quad, NULL));
428 }
429 PetscCall(PetscSectionGetStorageSize(subsection, &subNb));
430 if (num_points == 1) {
431 PetscInt rstride;
432
433 rstride = interleave_basis ? Ns : 1;
434 PetscCall(ISCreateStride(PETSC_COMM_SELF, subNb, roffset[0], rstride, &is_row));
435 roffset[0] += interleave_basis ? 1 : subNb;
436 } else {
437 PetscInt *rows;
438
439 PetscCall(PetscMalloc1(subNb, &rows));
440 for (PetscInt p = pStart; p < pEnd; p++) {
441 PetscInt subdof, dof, off, suboff, stride;
442
443 PetscCall(PetscSectionGetOffset(subsection, p, &suboff));
444 PetscCall(PetscSectionGetDof(subsection, p, &subdof));
445 PetscCall(PetscSectionGetOffset(section, p, &off));
446 PetscCall(PetscSectionGetDof(section, p, &dof));
447 PetscCheck(subdof * Ns == dof || !interleave_basis, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Basis cannot be interleaved");
448 stride = interleave_basis ? Ns : 1;
449 for (PetscInt k = 0; k < subdof; k++) rows[suboff + k] = off + roffset[p - pStart] + k * stride;
450 roffset[p - pStart] += interleave_basis ? 1 : subdof;
451 }
452 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, subNb, rows, PETSC_OWN_POINTER, &is_row));
453 }
454
455 sNpoints = 0;
456 if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &sNpoints, NULL, NULL));
457 sNcols = sNpoints * sNc;
458
459 if (!concatenate) {
460 PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, 1, &is_col));
461 coffset += uniform_points ? 0 : sNcols;
462 } else {
463 if (uniform_points && interleave_components) {
464 PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, Ns, &is_col));
465 coffset += 1;
466 } else {
467 PetscInt *cols;
468
469 PetscCall(PetscMalloc1(sNcols, &cols));
470 for (PetscInt p = 0, r = 0; p < sNpoints; p++) {
471 for (PetscInt c = 0; c < sNc; c++, r++) cols[r] = coffset + p * Nc + c;
472 }
473 coffset += uniform_points ? sNc : Nc * sNpoints + sNc;
474 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sNcols, cols, PETSC_OWN_POINTER, &is_col));
475 }
476 }
477 PetscCall(ISLocalToGlobalMappingCreateIS(is_row, &map_row[s]));
478 PetscCall(ISLocalToGlobalMappingCreateIS(is_col, &map_col[s]));
479 PetscCall(ISDestroy(&is_row));
480 PetscCall(ISDestroy(&is_col));
481 }
482 PetscCall(PetscFree(roffset));
483 PetscFunctionReturn(PETSC_SUCCESS);
484 }
485
PetscDualSpaceCreateFacetSubspace_Sum(PetscDualSpace sp,PetscInt f,PetscDualSpace * bdsp)486 static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Sum(PetscDualSpace sp, PetscInt f, PetscDualSpace *bdsp)
487 {
488 PetscInt k, Nc, Nk, Nknew, Ncnew, Ns;
489 PetscInt dim, pointDim = -1;
490 PetscInt depth, Ncopies;
491 PetscBool any;
492 DM dm, K;
493 DMPolytopeType ct;
494
495 PetscFunctionBegin;
496 PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
497 any = PETSC_FALSE;
498 for (PetscInt s = 0; s < Ns; s++) {
499 PetscDualSpace subsp, subpointsp;
500
501 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
502 PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
503 if (subpointsp) any = PETSC_TRUE;
504 }
505 if (!any) {
506 *bdsp = NULL;
507 PetscFunctionReturn(PETSC_SUCCESS);
508 }
509 PetscCall(PetscDualSpaceGetDM(sp, &dm));
510 PetscCall(DMGetDimension(dm, &dim));
511 PetscCall(DMPlexGetDepth(dm, &depth));
512 PetscCall(PetscDualSpaceDuplicate(sp, bdsp));
513 PetscCheck((depth == dim) || (depth == 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
514 pointDim = (depth == dim) ? (dim - 1) : 0;
515 PetscCall(DMPlexGetCellType(dm, f, &ct));
516 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
517 PetscCall(PetscDualSpaceSetDM(*bdsp, K));
518 PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
519 PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
520 PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
521 Ncopies = Nc / Nk;
522 PetscCall(PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew));
523 Ncnew = Nknew * Ncopies;
524 PetscCall(PetscDualSpaceSetNumComponents(*bdsp, Ncnew));
525 for (PetscInt s = 0; s < Ns; s++) {
526 PetscDualSpace subsp, subpointsp;
527
528 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
529 PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
530 if (subpointsp) {
531 PetscCall(PetscObjectReference((PetscObject)subpointsp));
532 } else {
533 // make an empty dual space
534 PetscInt subNc, subNcopies, subpointNc;
535
536 PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subsp), &subpointsp));
537 PetscCall(PetscDualSpaceGetNumComponents(subsp, &subNc));
538 subNcopies = subNc / Nk;
539 subpointNc = subNcopies * Nknew;
540 PetscCall(PetscDualSpaceSetType(subpointsp, PETSCDUALSPACESIMPLE));
541 PetscCall(PetscDualSpaceSimpleSetDimension(subpointsp, 0));
542 PetscCall(PetscDualSpaceSetFormDegree(subpointsp, k));
543 PetscCall(PetscDualSpaceSetNumComponents(subpointsp, subpointNc));
544 }
545 // K should be equal to subpointsp->dm, but we want it to be exactly the same
546 PetscCall(PetscObjectReference((PetscObject)K));
547 PetscCall(DMDestroy(&subpointsp->dm));
548 subpointsp->dm = K;
549 PetscCall(PetscDualSpaceSetUp(subpointsp));
550 PetscCall(PetscDualSpaceSumSetSubspace(*bdsp, s, subpointsp));
551 PetscCall(PetscDualSpaceDestroy(&subpointsp));
552 }
553 PetscCall(DMDestroy(&K));
554 PetscCall(PetscDualSpaceSetUp(*bdsp));
555 PetscFunctionReturn(PETSC_SUCCESS);
556 }
557
PetscDualSpaceSumIsUniform(PetscDualSpace sp,PetscBool * is_uniform)558 static PetscErrorCode PetscDualSpaceSumIsUniform(PetscDualSpace sp, PetscBool *is_uniform)
559 {
560 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
561 PetscBool uniform = PETSC_TRUE;
562
563 PetscFunctionBegin;
564 for (PetscInt s = 1; s < sum->numSumSpaces; s++) {
565 if (sum->sumspaces[s] != sum->sumspaces[0]) {
566 uniform = PETSC_FALSE;
567 break;
568 }
569 }
570 *is_uniform = uniform;
571 PetscFunctionReturn(PETSC_SUCCESS);
572 }
573
PetscDualSpaceGetSymmetries_Sum(PetscDualSpace sp,const PetscInt **** perms,const PetscScalar **** flips)574 static PetscErrorCode PetscDualSpaceGetSymmetries_Sum(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
575 {
576 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
577
578 PetscFunctionBegin;
579 if (!sum->symComputed) {
580 PetscInt Ns;
581 PetscBool any_perms = PETSC_FALSE;
582 PetscBool any_flips = PETSC_FALSE;
583 PetscInt ***symperms = NULL;
584 PetscScalar ***symflips = NULL;
585
586 sum->symComputed = PETSC_TRUE;
587 PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
588 for (PetscInt s = 0; s < Ns; s++) {
589 PetscDualSpace subsp;
590 const PetscInt ***sub_perms;
591 const PetscScalar ***sub_flips;
592
593 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
594 PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
595 if (sub_perms) any_perms = PETSC_TRUE;
596 if (sub_flips) any_flips = PETSC_TRUE;
597 }
598 if (any_perms || any_flips) {
599 DM K;
600 PetscInt pStart, pEnd, numPoints;
601 PetscInt spintdim;
602
603 PetscCall(PetscDualSpaceGetDM(sp, &K));
604 PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
605 numPoints = pEnd - pStart;
606 PetscCall(PetscCalloc1(numPoints, &symperms));
607 PetscCall(PetscCalloc1(numPoints, &symflips));
608 PetscCall(PetscDualSpaceGetBoundarySymmetries_Internal(sp, symperms, symflips));
609 // get interior symmetries
610 PetscCall(PetscDualSpaceGetInteriorDimension(sp, &spintdim));
611 if (spintdim) {
612 PetscInt groupSize;
613 PetscInt **cellPerms;
614 PetscScalar **cellFlips;
615 DMPolytopeType ct;
616
617 PetscCall(DMPlexGetCellType(K, 0, &ct));
618 groupSize = DMPolytopeTypeGetNumArrangements(ct);
619 sum->numSelfSym = groupSize;
620 sum->selfSymOff = groupSize / 2;
621 PetscCall(PetscCalloc1(groupSize, &cellPerms));
622 PetscCall(PetscCalloc1(groupSize, &cellFlips));
623 symperms[0] = &cellPerms[groupSize / 2];
624 symflips[0] = &cellFlips[groupSize / 2];
625 for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
626 PetscBool any_o_perms = PETSC_FALSE;
627 PetscBool any_o_flips = PETSC_FALSE;
628
629 for (PetscInt s = 0; s < Ns; s++) {
630 PetscDualSpace subsp;
631 const PetscInt ***sub_perms;
632 const PetscScalar ***sub_flips;
633
634 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
635 PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
636 if (sub_perms && sub_perms[0] && sub_perms[0][o]) any_o_perms = PETSC_TRUE;
637 if (sub_flips && sub_flips[0] && sub_flips[0][o]) any_o_flips = PETSC_TRUE;
638 }
639 if (any_o_perms) {
640 PetscInt *o_perm;
641
642 PetscCall(PetscMalloc1(spintdim, &o_perm));
643 for (PetscInt i = 0; i < spintdim; i++) o_perm[i] = i;
644 for (PetscInt s = 0; s < Ns; s++) {
645 PetscDualSpace subsp;
646 const PetscInt ***sub_perms;
647 const PetscScalar ***sub_flips;
648
649 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
650 PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
651 if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
652 PetscInt subspdim;
653 PetscInt *range, *domain;
654 PetscInt *range_mapped, *domain_mapped;
655
656 PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
657 PetscCall(PetscMalloc4(subspdim, &range, subspdim, &range_mapped, subspdim, &domain, subspdim, &domain_mapped));
658 for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
659 PetscCall(PetscArraycpy(range, sub_perms[0][o], subspdim));
660 PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
661 PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, range, range_mapped));
662 for (PetscInt i = 0; i < subspdim; i++) o_perm[domain_mapped[i]] = range_mapped[i];
663 PetscCall(PetscFree4(range, range_mapped, domain, domain_mapped));
664 }
665 }
666 symperms[0][o] = o_perm;
667 }
668 if (any_o_flips) {
669 PetscScalar *o_flip;
670
671 PetscCall(PetscMalloc1(spintdim, &o_flip));
672 for (PetscInt i = 0; i < spintdim; i++) o_flip[i] = 1.0;
673 for (PetscInt s = 0; s < Ns; s++) {
674 PetscDualSpace subsp;
675 const PetscInt ***sub_perms;
676 const PetscScalar ***sub_flips;
677
678 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
679 PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
680 if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
681 PetscInt subspdim;
682 PetscInt *domain;
683 PetscInt *domain_mapped;
684
685 PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
686 PetscCall(PetscMalloc2(subspdim, &domain, subspdim, &domain_mapped));
687 for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
688 PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
689 for (PetscInt i = 0; i < subspdim; i++) o_flip[domain_mapped[i]] = sub_perms[0][o][i];
690 PetscCall(PetscFree2(domain, domain_mapped));
691 }
692 }
693 symflips[0][o] = o_flip;
694 }
695 }
696 {
697 PetscBool any_perms = PETSC_FALSE;
698 PetscBool any_flips = PETSC_FALSE;
699 for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
700 if (symperms[0][o]) any_perms = PETSC_TRUE;
701 if (symflips[0][o]) any_flips = PETSC_TRUE;
702 }
703 if (!any_perms) {
704 PetscCall(PetscFree(cellPerms));
705 symperms[0] = NULL;
706 }
707 if (!any_flips) {
708 PetscCall(PetscFree(cellFlips));
709 symflips[0] = NULL;
710 }
711 }
712 }
713 if (!any_perms) {
714 PetscCall(PetscFree(symperms));
715 symperms = NULL;
716 }
717 if (!any_flips) {
718 PetscCall(PetscFree(symflips));
719 symflips = NULL;
720 }
721 }
722 sum->symperms = symperms;
723 sum->symflips = symflips;
724 }
725 if (perms) *perms = (const PetscInt ***)sum->symperms;
726 if (flips) *flips = (const PetscScalar ***)sum->symflips;
727 PetscFunctionReturn(PETSC_SUCCESS);
728 }
729
PetscDualSpaceSetUp_Sum(PetscDualSpace sp)730 static PetscErrorCode PetscDualSpaceSetUp_Sum(PetscDualSpace sp)
731 {
732 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
733 PetscBool concatenate = PETSC_TRUE;
734 PetscBool uniform;
735 PetscInt Ns, Nc, i, sum_Nc = 0;
736 PetscInt minNc, maxNc;
737 PetscInt minForm, maxForm, cdim, depth;
738 DM K;
739 PetscQuadrature *all_quads = NULL;
740 PetscQuadrature *int_quads = NULL;
741 Mat *all_mats = NULL;
742 Mat *int_mats = NULL;
743
744 PetscFunctionBegin;
745 if (sum->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
746 sum->setupcalled = PETSC_TRUE;
747
748 PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
749 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);
750
751 // step 1: make sure they share a DM
752 PetscCall(PetscDualSpaceGetDM(sp, &K));
753 PetscCall(DMGetCoordinateDim(K, &cdim));
754 PetscCall(DMPlexGetDepth(K, &depth));
755 PetscCall(PetscDualSpaceSumIsUniform(sp, &sp->uniform));
756 uniform = sp->uniform;
757 {
758 for (PetscInt s = 0; s < Ns; s++) {
759 PetscDualSpace subsp;
760 DM sub_K;
761
762 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
763 PetscCall(PetscDualSpaceSetUp(subsp));
764 PetscCall(PetscDualSpaceGetDM(subsp, &sub_K));
765 if (s == 0 && K == NULL) {
766 PetscCall(PetscDualSpaceSetDM(sp, sub_K));
767 K = sub_K;
768 }
769 PetscCheck(sub_K == K, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " does not have the same DM as the sum space", s);
770 }
771 }
772
773 // step 2: count components
774 PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
775 PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
776 minNc = Nc;
777 maxNc = Nc;
778 minForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MAX : sp->k;
779 maxForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MIN : sp->k;
780 for (i = 0; i < Ns; ++i) {
781 PetscInt sNc, formDegree;
782 PetscDualSpace si;
783
784 PetscCall(PetscDualSpaceSumGetSubspace(sp, i, &si));
785 PetscCall(PetscDualSpaceSetUp(si));
786 PetscCall(PetscDualSpaceGetNumComponents(si, &sNc));
787 if (sNc != Nc) PetscCheck(concatenate, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace as a different number of components but space does not concatenate components");
788 minNc = PetscMin(minNc, sNc);
789 maxNc = PetscMax(maxNc, sNc);
790 sum_Nc += sNc;
791 PetscCall(PetscDualSpaceGetFormDegree(si, &formDegree));
792 minForm = PetscMin(minForm, formDegree);
793 maxForm = PetscMax(maxForm, formDegree);
794 }
795 sp->k = (minForm != maxForm) ? PETSC_FORM_DEGREE_UNDEFINED : minForm;
796
797 if (concatenate) PetscCheck(sum_Nc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Total number of subspace components (%" PetscInt_FMT ") does not match number of target space components (%" PetscInt_FMT ").", sum_Nc, Nc);
798 else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");
799
800 /* A PetscDualSpace should have a fixed number of components, but
801 if the spaces we are combining have different form degrees, they will not
802 have the same number of components on subcomponents of the boundary,
803 so we do not try to create boundary dual spaces in this case */
804 if (sp->k != PETSC_FORM_DEGREE_UNDEFINED && depth > 0) {
805 PetscInt pStart, pEnd;
806 PetscInt *pStratStart, *pStratEnd;
807
808 PetscCall(DMPlexGetDepth(K, &depth));
809 PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
810 PetscCall(PetscCalloc1(pEnd, &sp->pointSpaces));
811 PetscCall(PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd));
812 for (PetscInt d = 0; d <= depth; ++d) PetscCall(DMPlexGetDepthStratum(K, d, &pStratStart[d], &pStratEnd[d]));
813
814 for (PetscInt p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
815 PetscReal v0[3];
816 DMPolytopeType ptype;
817 PetscReal J[9], detJ;
818 PetscInt q;
819
820 PetscCall(DMPlexComputeCellGeometryAffineFEM(K, p, v0, J, NULL, &detJ));
821 PetscCall(DMPlexGetCellType(K, p, &ptype));
822
823 /* compare to previous facets: if computed, reference that dualspace */
824 for (q = pStratStart[depth - 1]; q < p; q++) {
825 DMPolytopeType qtype;
826
827 PetscCall(DMPlexGetCellType(K, q, &qtype));
828 if (qtype == ptype) break;
829 }
830 if (q < p) { /* this facet has the same dual space as that one */
831 PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[q]));
832 sp->pointSpaces[p] = sp->pointSpaces[q];
833 continue;
834 }
835 /* if not, recursively compute this dual space */
836 PetscCall(PetscDualSpaceCreateFacetSubspace_Sum(sp, p, &sp->pointSpaces[p]));
837 }
838 for (PetscInt h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
839 PetscInt hd = depth - h;
840
841 for (PetscInt p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
842 PetscInt suppSize;
843 const PetscInt *supp;
844 DM qdm;
845 PetscDualSpace qsp, psp;
846 PetscInt c, coneSize, q;
847 const PetscInt *cone;
848 const PetscInt *refCone;
849
850 PetscCall(DMPlexGetSupportSize(K, p, &suppSize));
851 PetscCall(DMPlexGetSupport(K, p, &supp));
852 q = supp[0];
853 qsp = sp->pointSpaces[q];
854 PetscCall(DMPlexGetConeSize(K, q, &coneSize));
855 PetscCall(DMPlexGetCone(K, q, &cone));
856 for (c = 0; c < coneSize; c++)
857 if (cone[c] == p) break;
858 PetscCheck(c != coneSize, PetscObjectComm((PetscObject)K), PETSC_ERR_PLIB, "cone/support mismatch");
859 if (!qsp) {
860 sp->pointSpaces[p] = NULL;
861 continue;
862 }
863 PetscCall(PetscDualSpaceGetDM(qsp, &qdm));
864 PetscCall(DMPlexGetCone(qdm, 0, &refCone));
865 /* get the equivalent dual space from the support dual space */
866 PetscCall(PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp));
867 PetscCall(PetscObjectReference((PetscObject)psp));
868 sp->pointSpaces[p] = psp;
869 }
870 }
871 PetscCall(PetscFree2(pStratStart, pStratEnd));
872 }
873
874 sum->uniform = uniform;
875 PetscCall(PetscCalloc1(Ns, &sum->all_rows));
876 PetscCall(PetscCalloc1(Ns, &sum->all_cols));
877 PetscCall(PetscCalloc1(Ns, &sum->int_rows));
878 PetscCall(PetscCalloc1(Ns, &sum->int_cols));
879 PetscCall(PetscMalloc4(Ns, &all_quads, Ns, &all_mats, Ns, &int_quads, Ns, &int_mats));
880 {
881 // test for uniform all points and uniform interior points
882 PetscBool uniform_all = PETSC_TRUE;
883 PetscBool uniform_interior = PETSC_TRUE;
884 PetscQuadrature quad_all_first = NULL;
885 PetscQuadrature quad_interior_first = NULL;
886 PetscInt pStart, pEnd;
887
888 PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
889 PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &sp->pointSection));
890
891 for (PetscInt p = pStart; p < pEnd; p++) {
892 PetscInt full_dof = 0;
893
894 for (PetscInt s = 0; s < Ns; s++) {
895 PetscDualSpace subsp;
896 PetscSection subsection;
897 PetscInt dof;
898
899 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
900 PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
901 PetscCall(PetscSectionGetDof(subsection, p, &dof));
902 full_dof += dof;
903 }
904 PetscCall(PetscSectionSetDof(sp->pointSection, p, full_dof));
905 }
906 PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection));
907
908 for (PetscInt s = 0; s < Ns; s++) {
909 PetscDualSpace subsp;
910 PetscQuadrature subquad_all;
911 PetscQuadrature subquad_interior;
912 Mat submat_all;
913 Mat submat_interior;
914
915 PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
916 PetscCall(PetscDualSpaceGetAllData(subsp, &subquad_all, &submat_all));
917 PetscCall(PetscDualSpaceGetInteriorData(subsp, &subquad_interior, &submat_interior));
918 if (!s) {
919 quad_all_first = subquad_all;
920 quad_interior_first = subquad_interior;
921 } else {
922 if (subquad_all != quad_all_first) uniform_all = PETSC_FALSE;
923 if (subquad_interior != quad_interior_first) uniform_interior = PETSC_FALSE;
924 }
925 all_quads[s] = subquad_all;
926 int_quads[s] = subquad_interior;
927 all_mats[s] = submat_all;
928 int_mats[s] = submat_interior;
929 }
930 sum->uniform_all_points = uniform_all;
931 sum->uniform_interior_points = uniform_interior;
932
933 PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_TRUE, uniform_interior, sum->int_rows, sum->int_cols));
934 PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_interior, int_quads, &sp->intNodes));
935 if (sp->intNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, int_mats, sum->int_rows, sum->int_cols, sp->intNodes, &sp->intMat));
936
937 PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_FALSE, uniform_all, sum->all_rows, sum->all_cols));
938 PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_all, all_quads, &sp->allNodes));
939 if (sp->allNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, all_mats, sum->all_rows, sum->all_cols, sp->allNodes, &sp->allMat));
940 }
941 PetscCall(PetscFree4(all_quads, all_mats, int_quads, int_mats));
942 PetscCall(PetscDualSpaceComputeFunctionalsFromAllData(sp));
943 PetscFunctionReturn(PETSC_SUCCESS);
944 }
945
PetscDualSpaceSumView_Ascii(PetscDualSpace sp,PetscViewer v)946 static PetscErrorCode PetscDualSpaceSumView_Ascii(PetscDualSpace sp, PetscViewer v)
947 {
948 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
949 PetscBool concatenate = sum->concatenate;
950 PetscInt i, Ns = sum->numSumSpaces;
951
952 PetscFunctionBegin;
953 if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
954 else PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
955 for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
956 PetscCall(PetscViewerASCIIPushTab(v));
957 PetscCall(PetscDualSpaceView(sum->sumspaces[i], v));
958 PetscCall(PetscViewerASCIIPopTab(v));
959 }
960 PetscFunctionReturn(PETSC_SUCCESS);
961 }
962
PetscDualSpaceView_Sum(PetscDualSpace sp,PetscViewer viewer)963 static PetscErrorCode PetscDualSpaceView_Sum(PetscDualSpace sp, PetscViewer viewer)
964 {
965 PetscBool isascii;
966
967 PetscFunctionBegin;
968 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
969 if (isascii) PetscCall(PetscDualSpaceSumView_Ascii(sp, viewer));
970 PetscFunctionReturn(PETSC_SUCCESS);
971 }
972
PetscDualSpaceDestroy_Sum(PetscDualSpace sp)973 static PetscErrorCode PetscDualSpaceDestroy_Sum(PetscDualSpace sp)
974 {
975 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
976 PetscInt i, Ns = sum->numSumSpaces;
977
978 PetscFunctionBegin;
979 if (sum->symperms) {
980 PetscInt **selfSyms = sum->symperms[0];
981
982 if (selfSyms) {
983 PetscInt i, **allocated = &selfSyms[-sum->selfSymOff];
984
985 for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
986 PetscCall(PetscFree(allocated));
987 }
988 PetscCall(PetscFree(sum->symperms));
989 }
990 if (sum->symflips) {
991 PetscScalar **selfSyms = sum->symflips[0];
992
993 if (selfSyms) {
994 PetscInt i;
995 PetscScalar **allocated = &selfSyms[-sum->selfSymOff];
996
997 for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
998 PetscCall(PetscFree(allocated));
999 }
1000 PetscCall(PetscFree(sum->symflips));
1001 }
1002 for (i = 0; i < Ns; ++i) {
1003 PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[i]));
1004 if (sum->all_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_rows[i]));
1005 if (sum->all_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_cols[i]));
1006 if (sum->int_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_rows[i]));
1007 if (sum->int_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_cols[i]));
1008 }
1009 PetscCall(PetscFree(sum->sumspaces));
1010 PetscCall(PetscFree(sum->all_rows));
1011 PetscCall(PetscFree(sum->all_cols));
1012 PetscCall(PetscFree(sum->int_rows));
1013 PetscCall(PetscFree(sum->int_cols));
1014 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", NULL));
1015 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", NULL));
1016 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", NULL));
1017 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", NULL));
1018 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", NULL));
1019 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", NULL));
1020 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", NULL));
1021 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", NULL));
1022 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL));
1023 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL));
1024 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL));
1025 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL));
1026 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL));
1027 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL));
1028 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL));
1029 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL));
1030 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL));
1031 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL));
1032 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL));
1033 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL));
1034 PetscCall(PetscFree(sum));
1035 PetscFunctionReturn(PETSC_SUCCESS);
1036 }
1037
1038 /*@
1039 PetscDualSpaceSumSetInterleave - Set whether the basis functions and components of a uniform sum are interleaved
1040
1041 Logically collective
1042
1043 Input Parameters:
1044 + sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1045 . interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1046 - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1047 interleave the concatenated components
1048
1049 Level: developer
1050
1051 .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumGetInterleave()`
1052 @*/
PetscDualSpaceSumSetInterleave(PetscDualSpace sp,PetscBool interleave_basis,PetscBool interleave_components)1053 PetscErrorCode PetscDualSpaceSumSetInterleave(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1054 {
1055 PetscFunctionBegin;
1056 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
1057 PetscTryMethod(sp, "PetscDualSpaceSumSetInterleave_C", (PetscDualSpace, PetscBool, PetscBool), (sp, interleave_basis, interleave_components));
1058 PetscFunctionReturn(PETSC_SUCCESS);
1059 }
1060
PetscDualSpaceSumSetInterleave_Sum(PetscDualSpace sp,PetscBool interleave_basis,PetscBool interleave_components)1061 static PetscErrorCode PetscDualSpaceSumSetInterleave_Sum(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1062 {
1063 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
1064
1065 PetscFunctionBegin;
1066 sum->interleave_basis = interleave_basis;
1067 sum->interleave_components = interleave_components;
1068 PetscFunctionReturn(PETSC_SUCCESS);
1069 }
1070
1071 /*@
1072 PetscDualSpaceSumGetInterleave - Get whether the basis functions and components of a uniform sum are interleaved
1073
1074 Logically collective
1075
1076 Input Parameter:
1077 . sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1078
1079 Output Parameters:
1080 + interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1081 - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1082 interleave the concatenated components
1083
1084 Level: developer
1085
1086 .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumSetInterleave()`
1087 @*/
PetscDualSpaceSumGetInterleave(PetscDualSpace sp,PeOp PetscBool * interleave_basis,PeOp PetscBool * interleave_components)1088 PetscErrorCode PetscDualSpaceSumGetInterleave(PetscDualSpace sp, PeOp PetscBool *interleave_basis, PeOp PetscBool *interleave_components)
1089 {
1090 PetscFunctionBegin;
1091 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
1092 if (interleave_basis) PetscAssertPointer(interleave_basis, 2);
1093 if (interleave_components) PetscAssertPointer(interleave_components, 3);
1094 PetscTryMethod(sp, "PetscDualSpaceSumGetInterleave_C", (PetscDualSpace, PetscBool *, PetscBool *), (sp, interleave_basis, interleave_components));
1095 PetscFunctionReturn(PETSC_SUCCESS);
1096 }
1097
PetscDualSpaceSumGetInterleave_Sum(PetscDualSpace sp,PetscBool * interleave_basis,PetscBool * interleave_components)1098 static PetscErrorCode PetscDualSpaceSumGetInterleave_Sum(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1099 {
1100 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
1101
1102 PetscFunctionBegin;
1103 if (interleave_basis) *interleave_basis = sum->interleave_basis;
1104 if (interleave_components) *interleave_components = sum->interleave_components;
1105 PetscFunctionReturn(PETSC_SUCCESS);
1106 }
1107
1108 #define PetscDualSpaceSumPassthrough(sp, func, ...) \
1109 do { \
1110 PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data; \
1111 PetscBool is_uniform; \
1112 PetscCall(PetscDualSpaceSumIsUniform(sp, &is_uniform)); \
1113 if (is_uniform && sum->numSumSpaces > 0) { \
1114 PetscDualSpace subsp; \
1115 PetscCall(PetscDualSpaceSumGetSubspace(sp, 0, &subsp)); \
1116 PetscCall(func(subsp, __VA_ARGS__)); \
1117 } \
1118 } while (0)
1119
PetscDualSpaceLagrangeGetContinuity_Sum(PetscDualSpace sp,PetscBool * value)1120 static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Sum(PetscDualSpace sp, PetscBool *value)
1121 {
1122 PetscFunctionBegin;
1123 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetContinuity, value);
1124 PetscFunctionReturn(PETSC_SUCCESS);
1125 }
1126
PetscDualSpaceLagrangeSetContinuity_Sum(PetscDualSpace sp,PetscBool value)1127 static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Sum(PetscDualSpace sp, PetscBool value)
1128 {
1129 PetscFunctionBegin;
1130 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetContinuity, value);
1131 PetscFunctionReturn(PETSC_SUCCESS);
1132 }
1133
PetscDualSpaceLagrangeGetTensor_Sum(PetscDualSpace sp,PetscBool * value)1134 static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Sum(PetscDualSpace sp, PetscBool *value)
1135 {
1136 PetscFunctionBegin;
1137 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTensor, value);
1138 PetscFunctionReturn(PETSC_SUCCESS);
1139 }
1140
PetscDualSpaceLagrangeSetTensor_Sum(PetscDualSpace sp,PetscBool value)1141 static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Sum(PetscDualSpace sp, PetscBool value)
1142 {
1143 PetscFunctionBegin;
1144 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTensor, value);
1145 PetscFunctionReturn(PETSC_SUCCESS);
1146 }
1147
PetscDualSpaceLagrangeGetTrimmed_Sum(PetscDualSpace sp,PetscBool * value)1148 static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Sum(PetscDualSpace sp, PetscBool *value)
1149 {
1150 PetscFunctionBegin;
1151 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTrimmed, value);
1152 PetscFunctionReturn(PETSC_SUCCESS);
1153 }
1154
PetscDualSpaceLagrangeSetTrimmed_Sum(PetscDualSpace sp,PetscBool value)1155 static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Sum(PetscDualSpace sp, PetscBool value)
1156 {
1157 PetscFunctionBegin;
1158 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTrimmed, value);
1159 PetscFunctionReturn(PETSC_SUCCESS);
1160 }
1161
PetscDualSpaceLagrangeGetUseMoments_Sum(PetscDualSpace sp,PetscBool * value)1162 static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Sum(PetscDualSpace sp, PetscBool *value)
1163 {
1164 PetscFunctionBegin;
1165 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetUseMoments, value);
1166 PetscFunctionReturn(PETSC_SUCCESS);
1167 }
1168
PetscDualSpaceLagrangeSetUseMoments_Sum(PetscDualSpace sp,PetscBool value)1169 static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Sum(PetscDualSpace sp, PetscBool value)
1170 {
1171 PetscFunctionBegin;
1172 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetUseMoments, value);
1173 PetscFunctionReturn(PETSC_SUCCESS);
1174 }
1175
PetscDualSpaceLagrangeGetMomentOrder_Sum(PetscDualSpace sp,PetscInt * value)1176 static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Sum(PetscDualSpace sp, PetscInt *value)
1177 {
1178 PetscFunctionBegin;
1179 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetMomentOrder, value);
1180 PetscFunctionReturn(PETSC_SUCCESS);
1181 }
1182
PetscDualSpaceLagrangeSetMomentOrder_Sum(PetscDualSpace sp,PetscInt value)1183 static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Sum(PetscDualSpace sp, PetscInt value)
1184 {
1185 PetscFunctionBegin;
1186 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetMomentOrder, value);
1187 PetscFunctionReturn(PETSC_SUCCESS);
1188 }
1189
PetscDualSpaceLagrangeGetNodeType_Sum(PetscDualSpace sp,PetscDTNodeType * node_type,PetscBool * include_endpoints,PetscReal * exponent)1190 static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType *node_type, PetscBool *include_endpoints, PetscReal *exponent)
1191 {
1192 PetscFunctionBegin;
1193 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetNodeType, node_type, include_endpoints, exponent);
1194 PetscFunctionReturn(PETSC_SUCCESS);
1195 }
1196
PetscDualSpaceLagrangeSetNodeType_Sum(PetscDualSpace sp,PetscDTNodeType node_type,PetscBool include_endpoints,PetscReal exponent)1197 static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType node_type, PetscBool include_endpoints, PetscReal exponent)
1198 {
1199 PetscFunctionBegin;
1200 PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetNodeType, node_type, include_endpoints, exponent);
1201 PetscFunctionReturn(PETSC_SUCCESS);
1202 }
1203
PetscDualSpaceInitialize_Sum(PetscDualSpace sp)1204 static PetscErrorCode PetscDualSpaceInitialize_Sum(PetscDualSpace sp)
1205 {
1206 PetscFunctionBegin;
1207 sp->ops->destroy = PetscDualSpaceDestroy_Sum;
1208 sp->ops->view = PetscDualSpaceView_Sum;
1209 sp->ops->setfromoptions = NULL;
1210 sp->ops->duplicate = PetscDualSpaceDuplicate_Sum;
1211 sp->ops->setup = PetscDualSpaceSetUp_Sum;
1212 sp->ops->createheightsubspace = NULL;
1213 sp->ops->createpointsubspace = NULL;
1214 sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Sum;
1215 sp->ops->apply = PetscDualSpaceApplyDefault;
1216 sp->ops->applyall = PetscDualSpaceApplyAllDefault;
1217 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
1218 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
1219 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
1220
1221 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", PetscDualSpaceSumGetNumSubspaces_Sum));
1222 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", PetscDualSpaceSumSetNumSubspaces_Sum));
1223 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", PetscDualSpaceSumGetSubspace_Sum));
1224 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", PetscDualSpaceSumSetSubspace_Sum));
1225 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", PetscDualSpaceSumGetConcatenate_Sum));
1226 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", PetscDualSpaceSumSetConcatenate_Sum));
1227 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", PetscDualSpaceSumGetInterleave_Sum));
1228 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", PetscDualSpaceSumSetInterleave_Sum));
1229 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Sum));
1230 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Sum));
1231 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Sum));
1232 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Sum));
1233 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Sum));
1234 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Sum));
1235 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Sum));
1236 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Sum));
1237 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Sum));
1238 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Sum));
1239 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Sum));
1240 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Sum));
1241 PetscFunctionReturn(PETSC_SUCCESS);
1242 }
1243
1244 /*MC
1245 PETSCDUALSPACESUM = "sum" - A `PetscDualSpace` object that encapsulates a sum of subspaces.
1246
1247 Level: intermediate
1248
1249 Note:
1250 That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
1251 the direct sum of A and B will also have 2 components while the concatenated sum will have 4 components. In both cases A and B must be defined over the
1252 same reference element.
1253
1254 .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`, `PetscDualSpaceSumGetNumSubspaces()`, `PetscDualSpaceSumSetNumSubspaces()`,
1255 `PetscDualSpaceSumGetConcatenate()`, `PetscDualSpaceSumSetConcatenate()`, `PetscDualSpaceSumSetInterleave()`, `PetscDualSpaceSumGetInterleave()`
1256 M*/
PetscDualSpaceCreate_Sum(PetscDualSpace sp)1257 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Sum(PetscDualSpace sp)
1258 {
1259 PetscDualSpace_Sum *sum;
1260
1261 PetscFunctionBegin;
1262 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
1263 PetscCall(PetscNew(&sum));
1264 sum->numSumSpaces = PETSC_DEFAULT;
1265 sp->data = sum;
1266 sp->k = PETSC_FORM_DEGREE_UNDEFINED;
1267 PetscCall(PetscDualSpaceInitialize_Sum(sp));
1268 PetscFunctionReturn(PETSC_SUCCESS);
1269 }
1270
1271 /*@
1272 PetscDualSpaceCreateSum - Create a finite element dual basis that is the sum of other dual bases
1273
1274 Collective
1275
1276 Input Parameters:
1277 + numSubspaces - the number of spaces that will be added together
1278 . subspaces - an array of length `numSubspaces` of spaces
1279 - concatenate - if `PETSC_FALSE`, the sum-space has the same components as the individual dual spaces (`PetscDualSpaceGetNumComponents()`); if `PETSC_TRUE`, the individual components are concatenated to create a dual space with more components
1280
1281 Output Parameter:
1282 . sumSpace - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1283
1284 Level: advanced
1285
1286 .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCSPACESUM`
1287 @*/
PetscDualSpaceCreateSum(PetscInt numSubspaces,const PetscDualSpace subspaces[],PetscBool concatenate,PetscDualSpace * sumSpace)1288 PetscErrorCode PetscDualSpaceCreateSum(PetscInt numSubspaces, const PetscDualSpace subspaces[], PetscBool concatenate, PetscDualSpace *sumSpace)
1289 {
1290 PetscInt i, Nc = 0;
1291
1292 PetscFunctionBegin;
1293 PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
1294 PetscCall(PetscDualSpaceSetType(*sumSpace, PETSCDUALSPACESUM));
1295 PetscCall(PetscDualSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
1296 PetscCall(PetscDualSpaceSumSetConcatenate(*sumSpace, concatenate));
1297 for (i = 0; i < numSubspaces; ++i) {
1298 PetscInt sNc;
1299
1300 PetscCall(PetscDualSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
1301 PetscCall(PetscDualSpaceGetNumComponents(subspaces[i], &sNc));
1302 if (concatenate) Nc += sNc;
1303 else Nc = sNc;
1304
1305 if (i == 0) {
1306 DM dm;
1307
1308 PetscCall(PetscDualSpaceGetDM(subspaces[i], &dm));
1309 PetscCall(PetscDualSpaceSetDM(*sumSpace, dm));
1310 }
1311 }
1312 PetscCall(PetscDualSpaceSetNumComponents(*sumSpace, Nc));
1313 PetscFunctionReturn(PETSC_SUCCESS);
1314 }
1315