Actual source code: spacetensor.c

  1: #include <petsc/private/petscfeimpl.h>

  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];

 10:   PetscFunctionBegin;
 11:   PetscCall(PetscSpaceGetDegree(space, &degree, 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: }

 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;

 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, &deg, NULL));
 40:   if (Ns > 1) {
 41:     PetscSpace s0;

 43:     PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
 44:     for (i = 1; i < Ns; i++) {
 45:       PetscSpace si;

 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;

 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;

 80:       PetscCall(PetscSpaceTensorGetSubspace(sp, i, &subspace));
 81:       if (!subspace) {
 82:         char tprefix[128];

 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: }

 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;

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: }

120: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
121: {
122:   PetscBool iascii;

124:   PetscFunctionBegin;
125:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
126:   if (iascii) PetscCall(PetscSpaceTensorView_Ascii(sp, viewer));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

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;

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;

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;

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;

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;

179:       Ncprod = 1;
180:       for (PetscInt i = 0; i < Ns; i++) {
181:         PetscInt   Nvs, Ncs;
182:         PetscSpace si = NULL;

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));

193:         Nvsum += Nvs;
194:         Ncprod *= Ncs;
195:       }

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;

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;

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;

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: }

254: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
255: {
256:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
257:   PetscInt           Ns, i;

259:   PetscFunctionBegin;
260:   Ns = tens->numTensSpaces;
261:   if (tens->heightsubspaces) {
262:     PetscInt d;

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: }

279: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
280: {
281:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
282:   PetscInt           i, Ns, d;

284:   PetscFunctionBegin;
285:   PetscCall(PetscSpaceSetUp(sp));
286:   Ns = tens->numTensSpaces;
287:   d  = 1;
288:   for (i = 0; i < Ns; i++) {
289:     PetscInt id;

291:     PetscCall(PetscSpaceGetDimension(tens->tensspaces[i], &id));
292:     d *= id;
293:   }
294:   *dim = d;
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

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;

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;

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;

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;

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;

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;

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;

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;

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: }

430: /*@
431:   PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product space

433:   Input Parameters:
434: + sp            - the function space object
435: - numTensSpaces - the number of spaces

437:   Level: intermediate

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

442: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
443: @*/
444: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
445: {
446:   PetscFunctionBegin;
448:   PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
449:   PetscFunctionReturn(PETSC_SUCCESS);
450: }

452: /*@
453:   PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product space

455:   Input Parameter:
456: . sp - the function space object

458:   Output Parameter:
459: . numTensSpaces - the number of spaces

461:   Level: intermediate

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

466: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
467: @*/
468: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
469: {
470:   PetscFunctionBegin;
472:   PetscAssertPointer(numTensSpaces, 2);
473:   PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
474:   PetscFunctionReturn(PETSC_SUCCESS);
475: }

477: /*@
478:   PetscSpaceTensorSetSubspace - Set a space in the tensor product space

480:   Input Parameters:
481: + sp    - the function space object
482: . s     - The space number
483: - subsp - the number of spaces

485:   Level: intermediate

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

490: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
491: @*/
492: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
493: {
494:   PetscFunctionBegin;
497:   PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: /*@
502:   PetscSpaceTensorGetSubspace - Get a space in the tensor product space

504:   Input Parameters:
505: + sp - the function space object
506: - s  - The space number

508:   Output Parameter:
509: . subsp - the `PetscSpace`

511:   Level: intermediate

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

516: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
517: @*/
518: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
519: {
520:   PetscFunctionBegin;
522:   PetscAssertPointer(subsp, 3);
523:   PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
524:   PetscFunctionReturn(PETSC_SUCCESS);
525: }

527: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
528: {
529:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
530:   PetscInt           Ns;

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;

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: }

547: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
548: {
549:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;

551:   PetscFunctionBegin;
552:   *numTensSpaces = tens->numTensSpaces;
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
557: {
558:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
559:   PetscInt           Ns;

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: }

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;

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;

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: }

610: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
611: {
612:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
613:   PetscInt           Ns;

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: }

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: }

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.

644:   Level: intermediate

646: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceTensorSetSubspace()`,
647:           `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceTensorSetNumSubspaces()`
648: M*/

650: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
651: {
652:   PetscSpace_Tensor *tens;

654:   PetscFunctionBegin;
656:   PetscCall(PetscNew(&tens));
657:   sp->data = tens;

659:   tens->numTensSpaces = PETSC_DEFAULT;

661:   PetscCall(PetscSpaceInitialize_Tensor(sp));
662:   PetscFunctionReturn(PETSC_SUCCESS);
663: }