1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2
PetscSpaceTensorCreateSubspace(PetscSpace space,PetscInt Nvs,PetscInt Ncs,PetscSpace * subspace)3 static PetscErrorCode PetscSpaceTensorCreateSubspace(PetscSpace space, PetscInt Nvs, PetscInt Ncs, PetscSpace *subspace)
4 {
5 PetscInt degree;
6 const char *prefix;
7 const char *name;
8 char subname[PETSC_MAX_PATH_LEN];
9
10 PetscFunctionBegin;
11 PetscCall(PetscSpaceGetDegree(space, °ree, NULL));
12 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)space, &prefix));
13 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)space), subspace));
14 PetscCall(PetscSpaceSetType(*subspace, PETSCSPACEPOLYNOMIAL));
15 PetscCall(PetscSpaceSetNumVariables(*subspace, Nvs));
16 PetscCall(PetscSpaceSetNumComponents(*subspace, Ncs));
17 PetscCall(PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE));
18 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix));
19 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "tensorcomp_"));
20 if (((PetscObject)space)->name) {
21 PetscCall(PetscObjectGetName((PetscObject)space, &name));
22 PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s tensor component", name));
23 PetscCall(PetscObjectSetName((PetscObject)*subspace, subname));
24 } else PetscCall(PetscObjectSetName((PetscObject)*subspace, "tensor component"));
25 PetscFunctionReturn(PETSC_SUCCESS);
26 }
27
PetscSpaceSetFromOptions_Tensor(PetscSpace sp,PetscOptionItems PetscOptionsObject)28 static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscSpace sp, PetscOptionItems PetscOptionsObject)
29 {
30 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
31 PetscInt Ns, Nc, i, Nv, deg;
32 PetscBool uniform = PETSC_TRUE;
33
34 PetscFunctionBegin;
35 PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
36 if (!Nv) PetscFunctionReturn(PETSC_SUCCESS);
37 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
38 PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns));
39 PetscCall(PetscSpaceGetDegree(sp, °, NULL));
40 if (Ns > 1) {
41 PetscSpace s0;
42
43 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
44 for (i = 1; i < Ns; i++) {
45 PetscSpace si;
46
47 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
48 if (si != s0) {
49 uniform = PETSC_FALSE;
50 break;
51 }
52 }
53 }
54 Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv, 1) : Ns;
55 PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace tensor options");
56 PetscCall(PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL, 0));
57 PetscCall(PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL));
58 PetscOptionsHeadEnd();
59 PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space made up of %" PetscInt_FMT " spaces", Ns);
60 PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv);
61 if (Ns != tens->numTensSpaces) PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns));
62 if (uniform) {
63 PetscInt Nvs = Nv / Ns;
64 PetscInt Ncs;
65 PetscSpace subspace;
66
67 PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv);
68 Ncs = (PetscInt)PetscPowReal((PetscReal)Nc, 1. / Ns);
69 PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc);
70 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &subspace));
71 if (!subspace) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace));
72 else PetscCall(PetscObjectReference((PetscObject)subspace));
73 PetscCall(PetscSpaceSetFromOptions(subspace));
74 for (i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace));
75 PetscCall(PetscSpaceDestroy(&subspace));
76 } else {
77 for (i = 0; i < Ns; i++) {
78 PetscSpace subspace;
79
80 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &subspace));
81 if (!subspace) {
82 char tprefix[128];
83
84 PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace));
85 PetscCall(PetscSNPrintf(tprefix, 128, "%" PetscInt_FMT "_", i));
86 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix));
87 } else PetscCall(PetscObjectReference((PetscObject)subspace));
88 PetscCall(PetscSpaceSetFromOptions(subspace));
89 PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace));
90 PetscCall(PetscSpaceDestroy(&subspace));
91 }
92 }
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
PetscSpaceTensorView_Ascii(PetscSpace sp,PetscViewer v)96 static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
97 {
98 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
99 PetscBool uniform = PETSC_TRUE;
100 PetscInt Ns = tens->numTensSpaces, i, n;
101
102 PetscFunctionBegin;
103 for (i = 1; i < Ns; i++) {
104 if (tens->tensspaces[i] != tens->tensspaces[0]) {
105 uniform = PETSC_FALSE;
106 break;
107 }
108 }
109 if (uniform) PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces (all identical)\n", Ns));
110 else PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces\n", Ns));
111 n = uniform ? 1 : Ns;
112 for (i = 0; i < n; i++) {
113 PetscCall(PetscViewerASCIIPushTab(v));
114 PetscCall(PetscSpaceView(tens->tensspaces[i], v));
115 PetscCall(PetscViewerASCIIPopTab(v));
116 }
117 PetscFunctionReturn(PETSC_SUCCESS);
118 }
119
PetscSpaceView_Tensor(PetscSpace sp,PetscViewer viewer)120 static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
121 {
122 PetscBool isascii;
123
124 PetscFunctionBegin;
125 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
126 if (isascii) PetscCall(PetscSpaceTensorView_Ascii(sp, viewer));
127 PetscFunctionReturn(PETSC_SUCCESS);
128 }
129
PetscSpaceSetUp_Tensor(PetscSpace sp)130 static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
131 {
132 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
133 PetscInt Nc, Nv, Ns;
134 PetscBool uniform = PETSC_TRUE;
135 PetscInt deg, maxDeg;
136 PetscInt Ncprod;
137
138 PetscFunctionBegin;
139 if (tens->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
140 PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
141 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
142 PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns));
143 if (Ns == PETSC_DEFAULT) {
144 Ns = Nv;
145 PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns));
146 }
147 if (!Ns) {
148 SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
149 } else {
150 PetscSpace s0 = NULL;
151
152 PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv);
153 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
154 for (PetscInt i = 1; i < Ns; i++) {
155 PetscSpace si;
156
157 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
158 if (si != s0) {
159 uniform = PETSC_FALSE;
160 break;
161 }
162 }
163 if (uniform) {
164 PetscInt Nvs = Nv / Ns;
165 PetscInt Ncs;
166
167 PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv);
168 Ncs = (PetscInt)(PetscPowReal((PetscReal)Nc, 1. / Ns));
169 PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc);
170 if (!s0) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0));
171 else PetscCall(PetscObjectReference((PetscObject)s0));
172 PetscCall(PetscSpaceSetUp(s0));
173 for (PetscInt i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, s0));
174 PetscCall(PetscSpaceDestroy(&s0));
175 Ncprod = PetscPowInt(Ncs, Ns);
176 } else {
177 PetscInt Nvsum = 0;
178
179 Ncprod = 1;
180 for (PetscInt i = 0; i < Ns; i++) {
181 PetscInt Nvs, Ncs;
182 PetscSpace si = NULL;
183
184 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
185 if (!si) PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &si));
186 else PetscCall(PetscObjectReference((PetscObject)si));
187 PetscCall(PetscSpaceSetUp(si));
188 PetscCall(PetscSpaceTensorSetSubspace(sp, i, si));
189 PetscCall(PetscSpaceGetNumVariables(si, &Nvs));
190 PetscCall(PetscSpaceGetNumComponents(si, &Ncs));
191 PetscCall(PetscSpaceDestroy(&si));
192
193 Nvsum += Nvs;
194 Ncprod *= Ncs;
195 }
196
197 PetscCheck(Nvsum == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Sum of subspace variables %" PetscInt_FMT " does not equal the number of variables %" PetscInt_FMT, Nvsum, Nv);
198 PetscCheck(Nc % Ncprod == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Product of subspace components %" PetscInt_FMT " does not divide the number of components %" PetscInt_FMT, Ncprod, Nc);
199 }
200 if (Ncprod != Nc) {
201 PetscInt Ncopies = Nc / Ncprod;
202 PetscInt Nv = sp->Nv;
203 const char *prefix;
204 const char *name;
205 char subname[PETSC_MAX_PATH_LEN];
206 PetscSpace subsp;
207
208 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
209 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
210 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
211 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
212 if (((PetscObject)sp)->name) {
213 PetscCall(PetscObjectGetName((PetscObject)sp, &name));
214 PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
215 PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
216 } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
217 PetscCall(PetscSpaceSetType(subsp, PETSCSPACETENSOR));
218 PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
219 PetscCall(PetscSpaceSetNumComponents(subsp, Ncprod));
220 PetscCall(PetscSpaceTensorSetNumSubspaces(subsp, Ns));
221 for (PetscInt i = 0; i < Ns; i++) {
222 PetscSpace ssp;
223
224 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &ssp));
225 PetscCall(PetscSpaceTensorSetSubspace(subsp, i, ssp));
226 }
227 PetscCall(PetscSpaceSetUp(subsp));
228 PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
229 PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ncopies));
230 for (PetscInt i = 0; i < Ncopies; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
231 PetscCall(PetscSpaceDestroy(&subsp));
232 PetscCall(PetscSpaceSetUp(sp));
233 PetscFunctionReturn(PETSC_SUCCESS);
234 }
235 }
236 deg = PETSC_INT_MAX;
237 maxDeg = 0;
238 for (PetscInt i = 0; i < Ns; i++) {
239 PetscSpace si;
240 PetscInt iDeg, iMaxDeg;
241
242 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
243 PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg));
244 deg = PetscMin(deg, iDeg);
245 maxDeg += iMaxDeg;
246 }
247 sp->degree = deg;
248 sp->maxDegree = maxDeg;
249 tens->uniform = uniform;
250 tens->setupcalled = PETSC_TRUE;
251 PetscFunctionReturn(PETSC_SUCCESS);
252 }
253
PetscSpaceDestroy_Tensor(PetscSpace sp)254 static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
255 {
256 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
257 PetscInt Ns, i;
258
259 PetscFunctionBegin;
260 Ns = tens->numTensSpaces;
261 if (tens->heightsubspaces) {
262 PetscInt d;
263
264 /* sp->Nv is the spatial dimension, so it is equal to the number
265 * of subspaces on higher co-dimension points */
266 for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&tens->heightsubspaces[d]));
267 }
268 PetscCall(PetscFree(tens->heightsubspaces));
269 for (i = 0; i < Ns; i++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[i]));
270 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL));
271 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL));
272 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL));
273 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL));
274 PetscCall(PetscFree(tens->tensspaces));
275 PetscCall(PetscFree(tens));
276 PetscFunctionReturn(PETSC_SUCCESS);
277 }
278
PetscSpaceGetDimension_Tensor(PetscSpace sp,PetscInt * dim)279 static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
280 {
281 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
282 PetscInt i, Ns, d;
283
284 PetscFunctionBegin;
285 PetscCall(PetscSpaceSetUp(sp));
286 Ns = tens->numTensSpaces;
287 d = 1;
288 for (i = 0; i < Ns; i++) {
289 PetscInt id;
290
291 PetscCall(PetscSpaceGetDimension(tens->tensspaces[i], &id));
292 d *= id;
293 }
294 *dim = d;
295 PetscFunctionReturn(PETSC_SUCCESS);
296 }
297
PetscSpaceEvaluate_Tensor(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])298 static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
299 {
300 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
301 DM dm = sp->dm;
302 PetscInt Nc = sp->Nc;
303 PetscInt Nv = sp->Nv;
304 PetscInt Ns;
305 PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
306 PetscInt pdim;
307
308 PetscFunctionBegin;
309 if (!tens->setupcalled) {
310 PetscCall(PetscSpaceSetUp(sp));
311 PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
312 PetscFunctionReturn(PETSC_SUCCESS);
313 }
314 Ns = tens->numTensSpaces;
315 PetscCall(PetscSpaceGetDimension(sp, &pdim));
316 PetscCall(DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
317 if (B || D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
318 if (D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
319 if (H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
320 if (B) {
321 for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.;
322 }
323 if (D) {
324 for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.;
325 }
326 if (H) {
327 for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.;
328 }
329 for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
330 PetscInt sNv, sNc, spdim;
331 PetscInt vskip, cskip;
332
333 PetscCall(PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv));
334 PetscCall(PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc));
335 PetscCall(PetscSpaceGetDimension(tens->tensspaces[s], &spdim));
336 PetscCheck(!(pdim % vstep) && !(pdim % spdim), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", pdim %" PetscInt_FMT ", s %" PetscInt_FMT ", vstep %" PetscInt_FMT ", spdim %" PetscInt_FMT, Nv, Ns, pdim, s, vstep, spdim);
337 PetscCheck(!(Nc % cstep) && !(Nc % sNc), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", Nc %" PetscInt_FMT ", s %" PetscInt_FMT ", cstep %" PetscInt_FMT ", sNc %" PetscInt_FMT, Nv, Ns, Nc, s, cstep, spdim);
338 vskip = pdim / (vstep * spdim);
339 cskip = Nc / (cstep * sNc);
340 for (PetscInt p = 0; p < npoints; ++p) {
341 for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i];
342 }
343 PetscCall(PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH));
344 if (B) {
345 for (PetscInt k = 0; k < vskip; k++) {
346 for (PetscInt si = 0; si < spdim; si++) {
347 for (PetscInt j = 0; j < vstep; j++) {
348 PetscInt i = (k * spdim + si) * vstep + j;
349
350 for (PetscInt l = 0; l < cskip; l++) {
351 for (PetscInt sc = 0; sc < sNc; sc++) {
352 for (PetscInt m = 0; m < cstep; m++) {
353 PetscInt c = (l * sNc + sc) * cstep + m;
354
355 for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
356 }
357 }
358 }
359 }
360 }
361 }
362 }
363 if (D) {
364 for (PetscInt k = 0; k < vskip; k++) {
365 for (PetscInt si = 0; si < spdim; si++) {
366 for (PetscInt j = 0; j < vstep; j++) {
367 PetscInt i = (k * spdim + si) * vstep + j;
368
369 for (PetscInt l = 0; l < cskip; l++) {
370 for (PetscInt sc = 0; sc < sNc; sc++) {
371 for (PetscInt m = 0; m < cstep; m++) {
372 PetscInt c = (l * sNc + sc) * cstep + m;
373
374 for (PetscInt der = 0; der < Nv; der++) {
375 if (der >= d && der < d + sNv) {
376 for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
377 } else {
378 for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc];
379 }
380 }
381 }
382 }
383 }
384 }
385 }
386 }
387 }
388 if (H) {
389 for (PetscInt k = 0; k < vskip; k++) {
390 for (PetscInt si = 0; si < spdim; si++) {
391 for (PetscInt j = 0; j < vstep; j++) {
392 PetscInt i = (k * spdim + si) * vstep + j;
393
394 for (PetscInt l = 0; l < cskip; l++) {
395 for (PetscInt sc = 0; sc < sNc; sc++) {
396 for (PetscInt m = 0; m < cstep; m++) {
397 PetscInt c = (l * sNc + sc) * cstep + m;
398
399 for (PetscInt der = 0; der < Nv; der++) {
400 for (PetscInt der2 = 0; der2 < Nv; der2++) {
401 if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
402 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sH[(((spdim * p + si) * sNc + sc) * sNv + der - d) * sNv + der2 - d];
403 } else if (der >= d && der < d + sNv) {
404 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
405 } else if (der2 >= d && der2 < d + sNv) {
406 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der2 - d];
407 } else {
408 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
409 }
410 }
411 }
412 }
413 }
414 }
415 }
416 }
417 }
418 }
419 d += sNv;
420 vstep *= spdim;
421 cstep *= sNc;
422 }
423 if (H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
424 if (D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
425 if (B || D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
426 PetscCall(DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
427 PetscFunctionReturn(PETSC_SUCCESS);
428 }
429
430 /*@
431 PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product space
432
433 Input Parameters:
434 + sp - the function space object
435 - numTensSpaces - the number of spaces
436
437 Level: intermediate
438
439 Note:
440 The name NumSubspaces is misleading because it is actually setting the number of defining spaces of the tensor product space, not a number of Subspaces of it
441
442 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
443 @*/
PetscSpaceTensorSetNumSubspaces(PetscSpace sp,PetscInt numTensSpaces)444 PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
445 {
446 PetscFunctionBegin;
447 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
448 PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
449 PetscFunctionReturn(PETSC_SUCCESS);
450 }
451
452 /*@
453 PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product space
454
455 Input Parameter:
456 . sp - the function space object
457
458 Output Parameter:
459 . numTensSpaces - the number of spaces
460
461 Level: intermediate
462
463 Note:
464 The name NumSubspaces is misleading because it is actually getting the number of defining spaces of the tensor product space, not a number of Subspaces of it
465
466 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
467 @*/
PetscSpaceTensorGetNumSubspaces(PetscSpace sp,PetscInt * numTensSpaces)468 PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
469 {
470 PetscFunctionBegin;
471 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
472 PetscAssertPointer(numTensSpaces, 2);
473 PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
477 /*@
478 PetscSpaceTensorSetSubspace - Set a space in the tensor product space
479
480 Input Parameters:
481 + sp - the function space object
482 . s - The space number
483 - subsp - the number of spaces
484
485 Level: intermediate
486
487 Note:
488 The name SetSubspace is misleading because it is actually setting one of the defining spaces of the tensor product space, not a Subspace of it
489
490 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
491 @*/
PetscSpaceTensorSetSubspace(PetscSpace sp,PetscInt s,PetscSpace subsp)492 PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
493 {
494 PetscFunctionBegin;
495 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
496 if (subsp) PetscValidHeaderSpecific(subsp, PETSCSPACE_CLASSID, 3);
497 PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
498 PetscFunctionReturn(PETSC_SUCCESS);
499 }
500
501 /*@
502 PetscSpaceTensorGetSubspace - Get a space in the tensor product space
503
504 Input Parameters:
505 + sp - the function space object
506 - s - The space number
507
508 Output Parameter:
509 . subsp - the `PetscSpace`
510
511 Level: intermediate
512
513 Note:
514 The name GetSubspace is misleading because it is actually getting one of the defining spaces of the tensor product space, not a Subspace of it
515
516 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
517 @*/
PetscSpaceTensorGetSubspace(PetscSpace sp,PetscInt s,PetscSpace * subsp)518 PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
519 {
520 PetscFunctionBegin;
521 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
522 PetscAssertPointer(subsp, 3);
523 PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
524 PetscFunctionReturn(PETSC_SUCCESS);
525 }
526
PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space,PetscInt numTensSpaces)527 static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
528 {
529 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
530 PetscInt Ns;
531
532 PetscFunctionBegin;
533 PetscCheck(!tens->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
534 Ns = tens->numTensSpaces;
535 if (numTensSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
536 if (Ns >= 0) {
537 PetscInt s;
538
539 for (s = 0; s < Ns; s++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
540 PetscCall(PetscFree(tens->tensspaces));
541 }
542 Ns = tens->numTensSpaces = numTensSpaces;
543 PetscCall(PetscCalloc1(Ns, &tens->tensspaces));
544 PetscFunctionReturn(PETSC_SUCCESS);
545 }
546
PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space,PetscInt * numTensSpaces)547 static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
548 {
549 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
550
551 PetscFunctionBegin;
552 *numTensSpaces = tens->numTensSpaces;
553 PetscFunctionReturn(PETSC_SUCCESS);
554 }
555
PetscSpaceTensorSetSubspace_Tensor(PetscSpace space,PetscInt s,PetscSpace subspace)556 static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
557 {
558 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
559 PetscInt Ns;
560
561 PetscFunctionBegin;
562 PetscCheck(!tens->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
563 Ns = tens->numTensSpaces;
564 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
565 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
566 PetscCall(PetscObjectReference((PetscObject)subspace));
567 PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
568 tens->tensspaces[s] = subspace;
569 PetscFunctionReturn(PETSC_SUCCESS);
570 }
571
PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp,PetscInt height,PetscSpace * subsp)572 static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
573 {
574 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
575 PetscInt Nc, dim, order, i;
576 PetscSpace bsp;
577
578 PetscFunctionBegin;
579 PetscCall(PetscSpaceGetNumVariables(sp, &dim));
580 PetscCheck(tens->uniform && tens->numTensSpaces == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can only get a generic height subspace of a uniform tensor space of 1d spaces.");
581 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
582 PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
583 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);
584 if (!tens->heightsubspaces) PetscCall(PetscCalloc1(dim, &tens->heightsubspaces));
585 if (height <= dim) {
586 if (!tens->heightsubspaces[height - 1]) {
587 PetscSpace sub;
588 const char *name;
589
590 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &bsp));
591 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
592 PetscCall(PetscObjectGetName((PetscObject)sp, &name));
593 PetscCall(PetscObjectSetName((PetscObject)sub, name));
594 PetscCall(PetscSpaceSetType(sub, PETSCSPACETENSOR));
595 PetscCall(PetscSpaceSetNumComponents(sub, Nc));
596 PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE));
597 PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
598 PetscCall(PetscSpaceTensorSetNumSubspaces(sub, dim - height));
599 for (i = 0; i < dim - height; i++) PetscCall(PetscSpaceTensorSetSubspace(sub, i, bsp));
600 PetscCall(PetscSpaceSetUp(sub));
601 tens->heightsubspaces[height - 1] = sub;
602 }
603 *subsp = tens->heightsubspaces[height - 1];
604 } else {
605 *subsp = NULL;
606 }
607 PetscFunctionReturn(PETSC_SUCCESS);
608 }
609
PetscSpaceTensorGetSubspace_Tensor(PetscSpace space,PetscInt s,PetscSpace * subspace)610 static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
611 {
612 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
613 PetscInt Ns;
614
615 PetscFunctionBegin;
616 Ns = tens->numTensSpaces;
617 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
618 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
619 *subspace = tens->tensspaces[s];
620 PetscFunctionReturn(PETSC_SUCCESS);
621 }
622
PetscSpaceInitialize_Tensor(PetscSpace sp)623 static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
624 {
625 PetscFunctionBegin;
626 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor;
627 sp->ops->setup = PetscSpaceSetUp_Tensor;
628 sp->ops->view = PetscSpaceView_Tensor;
629 sp->ops->destroy = PetscSpaceDestroy_Tensor;
630 sp->ops->getdimension = PetscSpaceGetDimension_Tensor;
631 sp->ops->evaluate = PetscSpaceEvaluate_Tensor;
632 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
633 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor));
634 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor));
635 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor));
636 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor));
637 PetscFunctionReturn(PETSC_SUCCESS);
638 }
639
640 /*MC
641 PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space.
642 A tensor product is created of the components of the subspaces as well.
643
644 Level: intermediate
645
646 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceTensorSetSubspace()`,
647 `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceTensorSetNumSubspaces()`
648 M*/
649
PetscSpaceCreate_Tensor(PetscSpace sp)650 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
651 {
652 PetscSpace_Tensor *tens;
653
654 PetscFunctionBegin;
655 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
656 PetscCall(PetscNew(&tens));
657 sp->data = tens;
658
659 tens->numTensSpaces = PETSC_DEFAULT;
660
661 PetscCall(PetscSpaceInitialize_Tensor(sp));
662 PetscFunctionReturn(PETSC_SUCCESS);
663 }
664