1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2
3 /*@
4 PetscSpaceSumGetNumSubspaces - Get the number of spaces in the sum space
5
6 Input Parameter:
7 . sp - the function 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: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
18 @*/
PetscSpaceSumGetNumSubspaces(PetscSpace sp,PetscInt * numSumSpaces)19 PetscErrorCode PetscSpaceSumGetNumSubspaces(PetscSpace sp, PetscInt *numSumSpaces)
20 {
21 PetscFunctionBegin;
22 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
23 PetscAssertPointer(numSumSpaces, 2);
24 PetscTryMethod(sp, "PetscSpaceSumGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numSumSpaces));
25 PetscFunctionReturn(PETSC_SUCCESS);
26 }
27
28 /*@
29 PetscSpaceSumSetNumSubspaces - Set the number of spaces in the sum space
30
31 Input Parameters:
32 + sp - the function 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: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
41 @*/
PetscSpaceSumSetNumSubspaces(PetscSpace sp,PetscInt numSumSpaces)42 PetscErrorCode PetscSpaceSumSetNumSubspaces(PetscSpace sp, PetscInt numSumSpaces)
43 {
44 PetscFunctionBegin;
45 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
46 PetscTryMethod(sp, "PetscSpaceSumSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numSumSpaces));
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
50 /*@
51 PetscSpaceSumGetConcatenate - Get the concatenate flag for this space.
52
53 Input Parameter:
54 . sp - the function space object
55
56 Output Parameter:
57 . concatenate - flag indicating whether subspaces are concatenated.
58
59 Level: intermediate
60
61 Notes:
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: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetConcatenate()`
67 @*/
PetscSpaceSumGetConcatenate(PetscSpace sp,PetscBool * concatenate)68 PetscErrorCode PetscSpaceSumGetConcatenate(PetscSpace sp, PetscBool *concatenate)
69 {
70 PetscFunctionBegin;
71 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
72 PetscTryMethod(sp, "PetscSpaceSumGetConcatenate_C", (PetscSpace, PetscBool *), (sp, concatenate));
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
76 /*@
77 PetscSpaceSumSetConcatenate - Sets the concatenate flag for this space.
78
79 Input Parameters:
80 + sp - the function 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: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetConcatenate()`
91 @*/
PetscSpaceSumSetConcatenate(PetscSpace sp,PetscBool concatenate)92 PetscErrorCode PetscSpaceSumSetConcatenate(PetscSpace sp, PetscBool concatenate)
93 {
94 PetscFunctionBegin;
95 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
96 PetscTryMethod(sp, "PetscSpaceSumSetConcatenate_C", (PetscSpace, PetscBool), (sp, concatenate));
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99
100 /*@
101 PetscSpaceSumGetSubspace - Get a space in the sum space
102
103 Input Parameters:
104 + sp - the function space object
105 - s - The space number
106
107 Output Parameter:
108 . subsp - the `PetscSpace`
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: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
116 @*/
PetscSpaceSumGetSubspace(PetscSpace sp,PetscInt s,PetscSpace * subsp)117 PetscErrorCode PetscSpaceSumGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
118 {
119 PetscFunctionBegin;
120 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
121 PetscAssertPointer(subsp, 3);
122 PetscTryMethod(sp, "PetscSpaceSumGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
126 /*@
127 PetscSpaceSumSetSubspace - Set a space in the sum space
128
129 Input Parameters:
130 + sp - the function 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: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
140 @*/
PetscSpaceSumSetSubspace(PetscSpace sp,PetscInt s,PetscSpace subsp)141 PetscErrorCode PetscSpaceSumSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
142 {
143 PetscFunctionBegin;
144 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
145 if (subsp) PetscValidHeaderSpecific(subsp, PETSCSPACE_CLASSID, 3);
146 PetscTryMethod(sp, "PetscSpaceSumSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149
PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space,PetscInt * numSumSpaces)150 static PetscErrorCode PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space, PetscInt *numSumSpaces)
151 {
152 PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
153
154 PetscFunctionBegin;
155 *numSumSpaces = sum->numSumSpaces;
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space,PetscInt numSumSpaces)159 static PetscErrorCode PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space, PetscInt numSumSpaces)
160 {
161 PetscSpace_Sum *sum = (PetscSpace_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 if (Ns >= 0) {
168 PetscInt s;
169 for (s = 0; s < Ns; ++s) PetscCall(PetscSpaceDestroy(&sum->sumspaces[s]));
170 PetscCall(PetscFree(sum->sumspaces));
171 }
172
173 Ns = sum->numSumSpaces = numSumSpaces;
174 PetscCall(PetscCalloc1(Ns, &sum->sumspaces));
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
PetscSpaceSumGetConcatenate_Sum(PetscSpace sp,PetscBool * concatenate)178 static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp, PetscBool *concatenate)
179 {
180 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
181
182 PetscFunctionBegin;
183 *concatenate = sum->concatenate;
184 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
PetscSpaceSumSetConcatenate_Sum(PetscSpace sp,PetscBool concatenate)187 static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp, PetscBool concatenate)
188 {
189 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
190
191 PetscFunctionBegin;
192 PetscCheck(!sum->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change space concatenation after setup called.");
193
194 sum->concatenate = concatenate;
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
PetscSpaceSumGetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace * subspace)198 static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace *subspace)
199 {
200 PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
201 PetscInt Ns = sum->numSumSpaces;
202
203 PetscFunctionBegin;
204 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceSumSetNumSubspaces() first");
205 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
206
207 *subspace = sum->sumspaces[s];
208 PetscFunctionReturn(PETSC_SUCCESS);
209 }
210
PetscSpaceSumSetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace subspace)211 static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace subspace)
212 {
213 PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
214 PetscInt Ns = sum->numSumSpaces;
215
216 PetscFunctionBegin;
217 PetscCheck(!sum->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
218 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceSumSetNumSubspaces() first");
219 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
220
221 PetscCall(PetscObjectReference((PetscObject)subspace));
222 PetscCall(PetscSpaceDestroy(&sum->sumspaces[s]));
223 sum->sumspaces[s] = subspace;
224 PetscFunctionReturn(PETSC_SUCCESS);
225 }
226
PetscSpaceSetFromOptions_Sum(PetscSpace sp,PetscOptionItems PetscOptionsObject)227 static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscSpace sp, PetscOptionItems PetscOptionsObject)
228 {
229 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
230 PetscInt Ns, Nc, Nv, deg, i;
231 PetscBool concatenate = PETSC_TRUE;
232 const char *prefix;
233
234 PetscFunctionBegin;
235 PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
236 if (!Nv) PetscFunctionReturn(PETSC_SUCCESS);
237 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
238 PetscCall(PetscSpaceSumGetNumSubspaces(sp, &Ns));
239 PetscCall(PetscSpaceGetDegree(sp, °, NULL));
240 Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
241
242 PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace sum options");
243 PetscCall(PetscOptionsBoundedInt("-petscspace_sum_spaces", "The number of subspaces", "PetscSpaceSumSetNumSubspaces", Ns, &Ns, NULL, 0));
244 PetscCall(PetscOptionsBool("-petscspace_sum_concatenate", "Subspaces are concatenated components of the final space", "PetscSpaceSumSetFromOptions", concatenate, &concatenate, NULL));
245 PetscOptionsHeadEnd();
246
247 PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a sum space of %" PetscInt_FMT " spaces", Ns);
248 if (Ns != sum->numSumSpaces) PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ns));
249 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
250 for (i = 0; i < Ns; ++i) {
251 PetscInt sNv;
252 PetscSpace subspace;
253
254 PetscCall(PetscSpaceSumGetSubspace(sp, i, &subspace));
255 if (!subspace) {
256 char subspacePrefix[256];
257
258 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subspace));
259 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subspace, prefix));
260 PetscCall(PetscSNPrintf(subspacePrefix, 256, "sumcomp_%" PetscInt_FMT "_", i));
261 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, subspacePrefix));
262 } else PetscCall(PetscObjectReference((PetscObject)subspace));
263 PetscCall(PetscSpaceSetFromOptions(subspace));
264 PetscCall(PetscSpaceGetNumVariables(subspace, &sNv));
265 PetscCheck(sNv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " has not been set properly, number of variables is 0.", i);
266 PetscCall(PetscSpaceSumSetSubspace(sp, i, subspace));
267 PetscCall(PetscSpaceDestroy(&subspace));
268 }
269 PetscFunctionReturn(PETSC_SUCCESS);
270 }
271
PetscSpaceSetUp_Sum(PetscSpace sp)272 static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
273 {
274 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
275 PetscBool concatenate = PETSC_TRUE;
276 PetscBool uniform;
277 PetscInt Nv, Ns, Nc, i, sum_Nc = 0, deg = PETSC_INT_MAX, maxDeg = PETSC_INT_MIN;
278 PetscInt minNc, maxNc;
279
280 PetscFunctionBegin;
281 if (sum->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
282
283 PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
284 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
285 PetscCall(PetscSpaceSumGetNumSubspaces(sp, &Ns));
286 if (Ns == PETSC_DEFAULT) {
287 Ns = 1;
288 PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ns));
289 }
290 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);
291 uniform = PETSC_TRUE;
292 if (Ns) {
293 PetscSpace s0;
294
295 PetscCall(PetscSpaceSumGetSubspace(sp, 0, &s0));
296 for (PetscInt i = 1; i < Ns; i++) {
297 PetscSpace si;
298
299 PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
300 if (si != s0) {
301 uniform = PETSC_FALSE;
302 break;
303 }
304 }
305 }
306
307 minNc = Nc;
308 maxNc = Nc;
309 for (i = 0; i < Ns; ++i) {
310 PetscInt sNv, sNc, iDeg, iMaxDeg;
311 PetscSpace si;
312
313 PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
314 PetscCall(PetscSpaceSetUp(si));
315 PetscCall(PetscSpaceGetNumVariables(si, &sNv));
316 PetscCheck(sNv == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " has %" PetscInt_FMT " variables, space has %" PetscInt_FMT ".", i, sNv, Nv);
317 PetscCall(PetscSpaceGetNumComponents(si, &sNc));
318 if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
319 minNc = PetscMin(minNc, sNc);
320 maxNc = PetscMax(maxNc, sNc);
321 sum_Nc += sNc;
322 PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
323 PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg));
324 deg = PetscMin(deg, iDeg);
325 maxDeg = PetscMax(maxDeg, iMaxDeg);
326 }
327
328 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);
329 else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");
330
331 sp->degree = deg;
332 sp->maxDegree = maxDeg;
333 sum->concatenate = concatenate;
334 sum->uniform = uniform;
335 sum->setupcalled = PETSC_TRUE;
336 PetscFunctionReturn(PETSC_SUCCESS);
337 }
338
PetscSpaceSumView_Ascii(PetscSpace sp,PetscViewer v)339 static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp, PetscViewer v)
340 {
341 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
342 PetscBool concatenate = sum->concatenate;
343 PetscInt i, Ns = sum->numSumSpaces;
344
345 PetscFunctionBegin;
346 if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
347 else PetscCall(PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
348 for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
349 PetscCall(PetscViewerASCIIPushTab(v));
350 PetscCall(PetscSpaceView(sum->sumspaces[i], v));
351 PetscCall(PetscViewerASCIIPopTab(v));
352 }
353 PetscFunctionReturn(PETSC_SUCCESS);
354 }
355
PetscSpaceView_Sum(PetscSpace sp,PetscViewer viewer)356 static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp, PetscViewer viewer)
357 {
358 PetscBool isascii;
359
360 PetscFunctionBegin;
361 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
362 if (isascii) PetscCall(PetscSpaceSumView_Ascii(sp, viewer));
363 PetscFunctionReturn(PETSC_SUCCESS);
364 }
365
PetscSpaceDestroy_Sum(PetscSpace sp)366 static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
367 {
368 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
369 PetscInt i, Ns = sum->numSumSpaces;
370
371 PetscFunctionBegin;
372 for (i = 0; i < Ns; ++i) PetscCall(PetscSpaceDestroy(&sum->sumspaces[i]));
373 PetscCall(PetscFree(sum->sumspaces));
374 if (sum->heightsubspaces) {
375 PetscInt d;
376
377 /* sp->Nv is the spatial dimension, so it is equal to the number
378 * of subspaces on higher co-dimension points */
379 for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&sum->heightsubspaces[d]));
380 }
381 PetscCall(PetscFree(sum->heightsubspaces));
382 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", NULL));
383 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", NULL));
384 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", NULL));
385 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", NULL));
386 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", NULL));
387 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", NULL));
388 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetInterleave_C", NULL));
389 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetInterleave_C", NULL));
390 PetscCall(PetscFree(sum));
391 PetscFunctionReturn(PETSC_SUCCESS);
392 }
393
PetscSpaceGetDimension_Sum(PetscSpace sp,PetscInt * dim)394 static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp, PetscInt *dim)
395 {
396 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
397 PetscInt i, d = 0, Ns = sum->numSumSpaces;
398
399 PetscFunctionBegin;
400 if (!sum->setupcalled) {
401 PetscCall(PetscSpaceSetUp(sp));
402 PetscCall(PetscSpaceGetDimension(sp, dim));
403 PetscFunctionReturn(PETSC_SUCCESS);
404 }
405
406 for (i = 0; i < Ns; ++i) {
407 PetscInt id;
408
409 PetscCall(PetscSpaceGetDimension(sum->sumspaces[i], &id));
410 d += id;
411 }
412
413 *dim = d;
414 PetscFunctionReturn(PETSC_SUCCESS);
415 }
416
PetscSpaceEvaluate_Sum(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])417 static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
418 {
419 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
420 PetscBool concatenate = sum->concatenate;
421 DM dm = sp->dm;
422 PetscInt Nc = sp->Nc, Nv = sp->Nv, Ns = sum->numSumSpaces;
423 PetscInt i, s, offset, ncoffset, pdimfull, numelB, numelD, numelH;
424 PetscReal *sB = NULL, *sD = NULL, *sH = NULL;
425
426 PetscFunctionBegin;
427 if (!sum->setupcalled) {
428 PetscCall(PetscSpaceSetUp(sp));
429 PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
430 PetscFunctionReturn(PETSC_SUCCESS);
431 }
432 PetscCall(PetscSpaceGetDimension(sp, &pdimfull));
433 numelB = npoints * pdimfull * Nc;
434 numelD = numelB * Nv;
435 numelH = numelD * Nv;
436 if (B || D || H) PetscCall(DMGetWorkArray(dm, numelB, MPIU_REAL, &sB));
437 if (D || H) PetscCall(DMGetWorkArray(dm, numelD, MPIU_REAL, &sD));
438 if (H) PetscCall(DMGetWorkArray(dm, numelH, MPIU_REAL, &sH));
439 if (B)
440 for (i = 0; i < numelB; ++i) B[i] = 0.;
441 if (D)
442 for (i = 0; i < numelD; ++i) D[i] = 0.;
443 if (H)
444 for (i = 0; i < numelH; ++i) H[i] = 0.;
445
446 for (s = 0, offset = 0, ncoffset = 0; s < Ns; ++s) {
447 PetscInt sNv, spdim, sNc, p;
448
449 PetscCall(PetscSpaceGetNumVariables(sum->sumspaces[s], &sNv));
450 PetscCall(PetscSpaceGetNumComponents(sum->sumspaces[s], &sNc));
451 PetscCall(PetscSpaceGetDimension(sum->sumspaces[s], &spdim));
452 PetscCheck(offset + spdim <= pdimfull, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspace dimensions exceed target space dimension.");
453 if (s == 0 || !sum->uniform) PetscCall(PetscSpaceEvaluate(sum->sumspaces[s], npoints, points, sB, sD, sH));
454 if (B || D || H) {
455 for (p = 0; p < npoints; ++p) {
456 PetscInt j;
457
458 for (j = 0; j < spdim; ++j) {
459 PetscInt c;
460 PetscInt b = sum->interleave_basis ? (j * Ns + s) : (j + offset);
461
462 for (c = 0; c < sNc; ++c) {
463 PetscInt compoffset, BInd, sBInd;
464
465 compoffset = concatenate ? (sum->interleave_components ? (c * Ns + s) : (c + ncoffset)) : c;
466 BInd = (p * pdimfull + b) * Nc + compoffset;
467 sBInd = (p * spdim + j) * sNc + c;
468 if (B) B[BInd] = sB[sBInd];
469 if (D || H) {
470 PetscInt v;
471
472 for (v = 0; v < Nv; ++v) {
473 PetscInt DInd, sDInd;
474
475 DInd = BInd * Nv + v;
476 sDInd = sBInd * Nv + v;
477 if (D) D[DInd] = sD[sDInd];
478 if (H) {
479 PetscInt v2;
480
481 for (v2 = 0; v2 < Nv; ++v2) {
482 PetscInt HInd, sHInd;
483
484 HInd = DInd * Nv + v2;
485 sHInd = sDInd * Nv + v2;
486 H[HInd] = sH[sHInd];
487 }
488 }
489 }
490 }
491 }
492 }
493 }
494 }
495 offset += spdim;
496 ncoffset += sNc;
497 }
498
499 if (H) PetscCall(DMRestoreWorkArray(dm, numelH, MPIU_REAL, &sH));
500 if (D || H) PetscCall(DMRestoreWorkArray(dm, numelD, MPIU_REAL, &sD));
501 if (B || D || H) PetscCall(DMRestoreWorkArray(dm, numelB, MPIU_REAL, &sB));
502 PetscFunctionReturn(PETSC_SUCCESS);
503 }
504
PetscSpaceGetHeightSubspace_Sum(PetscSpace sp,PetscInt height,PetscSpace * subsp)505 static PetscErrorCode PetscSpaceGetHeightSubspace_Sum(PetscSpace sp, PetscInt height, PetscSpace *subsp)
506 {
507 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
508 PetscInt Nc, dim, order;
509 PetscBool tensor;
510
511 PetscFunctionBegin;
512 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
513 PetscCall(PetscSpaceGetNumVariables(sp, &dim));
514 PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
515 PetscCall(PetscSpacePolynomialGetTensor(sp, &tensor));
516 PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
517 if (!sum->heightsubspaces) PetscCall(PetscCalloc1(dim, &sum->heightsubspaces));
518 if (height <= dim) {
519 if (!sum->heightsubspaces[height - 1]) {
520 PetscSpace sub;
521 const char *name;
522
523 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
524 PetscCall(PetscObjectGetName((PetscObject)sp, &name));
525 PetscCall(PetscObjectSetName((PetscObject)sub, name));
526 PetscCall(PetscSpaceSetType(sub, PETSCSPACESUM));
527 PetscCall(PetscSpaceSumSetNumSubspaces(sub, sum->numSumSpaces));
528 PetscCall(PetscSpaceSumSetConcatenate(sub, sum->concatenate));
529 PetscCall(PetscSpaceSetNumComponents(sub, Nc));
530 PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
531 for (PetscInt i = 0; i < sum->numSumSpaces; i++) {
532 PetscSpace subh;
533
534 PetscCall(PetscSpaceGetHeightSubspace(sum->sumspaces[i], height, &subh));
535 PetscCall(PetscSpaceSumSetSubspace(sub, i, subh));
536 }
537 PetscCall(PetscSpaceSetUp(sub));
538 sum->heightsubspaces[height - 1] = sub;
539 }
540 *subsp = sum->heightsubspaces[height - 1];
541 } else {
542 *subsp = NULL;
543 }
544 PetscFunctionReturn(PETSC_SUCCESS);
545 }
546
547 /*@
548 PetscSpaceSumSetInterleave - Set whether the basis functions and components of a uniform sum are interleaved
549
550 Logically collective
551
552 Input Parameters:
553 + sp - a `PetscSpace` of type `PETSCSPACESUM`
554 . interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
555 - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscSpaceSumGetConcatenate()`),
556 interleave the concatenated components
557
558 Level: developer
559
560 .seealso: `PetscSpace`, `PETSCSPACESUM`, `PETSCFEVECTOR`, `PetscSpaceSumGetInterleave()`
561 @*/
PetscSpaceSumSetInterleave(PetscSpace sp,PetscBool interleave_basis,PetscBool interleave_components)562 PetscErrorCode PetscSpaceSumSetInterleave(PetscSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
563 {
564 PetscFunctionBegin;
565 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
566 PetscTryMethod(sp, "PetscSpaceSumSetInterleave_C", (PetscSpace, PetscBool, PetscBool), (sp, interleave_basis, interleave_components));
567 PetscFunctionReturn(PETSC_SUCCESS);
568 }
569
PetscSpaceSumSetInterleave_Sum(PetscSpace sp,PetscBool interleave_basis,PetscBool interleave_components)570 static PetscErrorCode PetscSpaceSumSetInterleave_Sum(PetscSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
571 {
572 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
573
574 PetscFunctionBegin;
575 sum->interleave_basis = interleave_basis;
576 sum->interleave_components = interleave_components;
577 PetscFunctionReturn(PETSC_SUCCESS);
578 }
579
580 /*@
581 PetscSpaceSumGetInterleave - Get whether the basis functions and components of a uniform sum are interleaved
582
583 Logically collective
584
585 Input Parameter:
586 . sp - a `PetscSpace` of type `PETSCSPACESUM`
587
588 Output Parameters:
589 + interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
590 - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscSpaceSumGetConcatenate()`),
591 interleave the concatenated components
592
593 Level: developer
594
595 .seealso: `PetscSpace`, `PETSCSPACESUM`, `PETSCFEVECTOR`, `PetscSpaceSumSetInterleave()`
596 @*/
PetscSpaceSumGetInterleave(PetscSpace sp,PeOp PetscBool * interleave_basis,PeOp PetscBool * interleave_components)597 PetscErrorCode PetscSpaceSumGetInterleave(PetscSpace sp, PeOp PetscBool *interleave_basis, PeOp PetscBool *interleave_components)
598 {
599 PetscFunctionBegin;
600 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
601 if (interleave_basis) PetscAssertPointer(interleave_basis, 2);
602 if (interleave_components) PetscAssertPointer(interleave_components, 3);
603 PetscTryMethod(sp, "PetscSpaceSumGetInterleave_C", (PetscSpace, PetscBool *, PetscBool *), (sp, interleave_basis, interleave_components));
604 PetscFunctionReturn(PETSC_SUCCESS);
605 }
606
PetscSpaceSumGetInterleave_Sum(PetscSpace sp,PetscBool * interleave_basis,PetscBool * interleave_components)607 static PetscErrorCode PetscSpaceSumGetInterleave_Sum(PetscSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
608 {
609 PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
610
611 PetscFunctionBegin;
612 if (interleave_basis) *interleave_basis = sum->interleave_basis;
613 if (interleave_components) *interleave_components = sum->interleave_components;
614 PetscFunctionReturn(PETSC_SUCCESS);
615 }
616
PetscSpaceInitialize_Sum(PetscSpace sp)617 static PetscErrorCode PetscSpaceInitialize_Sum(PetscSpace sp)
618 {
619 PetscFunctionBegin;
620 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
621 sp->ops->setup = PetscSpaceSetUp_Sum;
622 sp->ops->view = PetscSpaceView_Sum;
623 sp->ops->destroy = PetscSpaceDestroy_Sum;
624 sp->ops->getdimension = PetscSpaceGetDimension_Sum;
625 sp->ops->evaluate = PetscSpaceEvaluate_Sum;
626 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Sum;
627
628 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", PetscSpaceSumGetNumSubspaces_Sum));
629 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", PetscSpaceSumSetNumSubspaces_Sum));
630 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", PetscSpaceSumGetSubspace_Sum));
631 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", PetscSpaceSumSetSubspace_Sum));
632 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", PetscSpaceSumGetConcatenate_Sum));
633 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", PetscSpaceSumSetConcatenate_Sum));
634 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetInterleave_C", PetscSpaceSumGetInterleave_Sum));
635 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetInterleave_C", PetscSpaceSumSetInterleave_Sum));
636 PetscFunctionReturn(PETSC_SUCCESS);
637 }
638
639 /*MC
640 PETSCSPACESUM = "sum" - A `PetscSpace` object that encapsulates a sum of subspaces.
641
642 Level: intermediate
643
644 Note:
645 That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
646 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
647 same number of variables.
648
649 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceSumGetNumSubspaces()`, `PetscSpaceSumSetNumSubspaces()`,
650 `PetscSpaceSumGetConcatenate()`, `PetscSpaceSumSetConcatenate()`
651 M*/
PetscSpaceCreate_Sum(PetscSpace sp)652 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
653 {
654 PetscSpace_Sum *sum;
655
656 PetscFunctionBegin;
657 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
658 PetscCall(PetscNew(&sum));
659 sum->numSumSpaces = PETSC_DEFAULT;
660 sp->data = sum;
661 PetscCall(PetscSpaceInitialize_Sum(sp));
662 PetscFunctionReturn(PETSC_SUCCESS);
663 }
664
PetscSpaceCreateSum(PetscInt numSubspaces,const PetscSpace subspaces[],PetscBool concatenate,PetscSpace * sumSpace)665 PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces, const PetscSpace subspaces[], PetscBool concatenate, PetscSpace *sumSpace)
666 {
667 PetscInt i, Nv, Nc = 0;
668
669 PetscFunctionBegin;
670 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
671 PetscCall(PetscSpaceSetType(*sumSpace, PETSCSPACESUM));
672 PetscCall(PetscSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
673 PetscCall(PetscSpaceSumSetConcatenate(*sumSpace, concatenate));
674 for (i = 0; i < numSubspaces; ++i) {
675 PetscInt sNc;
676
677 PetscCall(PetscSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
678 PetscCall(PetscSpaceGetNumComponents(subspaces[i], &sNc));
679 if (concatenate) Nc += sNc;
680 else Nc = sNc;
681 }
682 PetscCall(PetscSpaceGetNumVariables(subspaces[0], &Nv));
683 PetscCall(PetscSpaceSetNumComponents(*sumSpace, Nc));
684 PetscCall(PetscSpaceSetNumVariables(*sumSpace, Nv));
685 PetscCall(PetscSpaceSetUp(*sumSpace));
686 PetscFunctionReturn(PETSC_SUCCESS);
687 }
688