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:   PetscSpaceGetDegree(space, &degree, NULL);
 11:   PetscObjectGetOptionsPrefix((PetscObject)space, &prefix);
 12:   PetscSpaceCreate(PetscObjectComm((PetscObject)space), subspace);
 13:   PetscSpaceSetType(*subspace, PETSCSPACEPOLYNOMIAL);
 14:   PetscSpaceSetNumVariables(*subspace, Nvs);
 15:   PetscSpaceSetNumComponents(*subspace, Ncs);
 16:   PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE);
 17:   PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix);
 18:   PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "tensorcomp_");
 19:   if (((PetscObject)space)->name) {
 20:     PetscObjectGetName((PetscObject)space, &name);
 21:     PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s tensor component", name);
 22:     PetscObjectSetName((PetscObject)*subspace, subname);
 23:   } else PetscObjectSetName((PetscObject)*subspace, "tensor component");
 24:   return 0;
 25: }

 27: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
 28: {
 29:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
 30:   PetscInt           Ns, Nc, i, Nv, deg;
 31:   PetscBool          uniform = PETSC_TRUE;

 33:   PetscSpaceGetNumVariables(sp, &Nv);
 34:   if (!Nv) return 0;
 35:   PetscSpaceGetNumComponents(sp, &Nc);
 36:   PetscSpaceTensorGetNumSubspaces(sp, &Ns);
 37:   PetscSpaceGetDegree(sp, &deg, NULL);
 38:   if (Ns > 1) {
 39:     PetscSpace s0;

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

 45:       PetscSpaceTensorGetSubspace(sp, i, &si);
 46:       if (si != s0) {
 47:         uniform = PETSC_FALSE;
 48:         break;
 49:       }
 50:     }
 51:   }
 52:   Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv, 1) : Ns;
 53:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace tensor options");
 54:   PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL, 0);
 55:   PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL);
 56:   PetscOptionsHeadEnd();
 59:   if (Ns != tens->numTensSpaces) PetscSpaceTensorSetNumSubspaces(sp, Ns);
 60:   if (uniform) {
 61:     PetscInt   Nvs = Nv / Ns;
 62:     PetscInt   Ncs;
 63:     PetscSpace subspace;

 66:     Ncs = (PetscInt)PetscPowReal((PetscReal)Nc, 1. / Ns);
 68:     PetscSpaceTensorGetSubspace(sp, 0, &subspace);
 69:     if (!subspace) PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace);
 70:     else PetscObjectReference((PetscObject)subspace);
 71:     PetscSpaceSetFromOptions(subspace);
 72:     for (i = 0; i < Ns; i++) PetscSpaceTensorSetSubspace(sp, i, subspace);
 73:     PetscSpaceDestroy(&subspace);
 74:   } else {
 75:     for (i = 0; i < Ns; i++) {
 76:       PetscSpace subspace;

 78:       PetscSpaceTensorGetSubspace(sp, i, &subspace);
 79:       if (!subspace) {
 80:         char tprefix[128];

 82:         PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace);
 83:         PetscSNPrintf(tprefix, 128, "%d_", (int)i);
 84:         PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix);
 85:       } else PetscObjectReference((PetscObject)subspace);
 86:       PetscSpaceSetFromOptions(subspace);
 87:       PetscSpaceTensorSetSubspace(sp, i, subspace);
 88:       PetscSpaceDestroy(&subspace);
 89:     }
 90:   }
 91:   return 0;
 92: }

 94: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
 95: {
 96:   PetscSpace_Tensor *tens    = (PetscSpace_Tensor *)sp->data;
 97:   PetscBool          uniform = PETSC_TRUE;
 98:   PetscInt           Ns      = tens->numTensSpaces, i, n;

100:   for (i = 1; i < Ns; i++) {
101:     if (tens->tensspaces[i] != tens->tensspaces[0]) {
102:       uniform = PETSC_FALSE;
103:       break;
104:     }
105:   }
106:   if (uniform) PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces (all identical)\n", Ns);
107:   else PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces\n", Ns);
108:   n = uniform ? 1 : Ns;
109:   for (i = 0; i < n; i++) {
110:     PetscViewerASCIIPushTab(v);
111:     PetscSpaceView(tens->tensspaces[i], v);
112:     PetscViewerASCIIPopTab(v);
113:   }
114:   return 0;
115: }

117: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
118: {
119:   PetscBool iascii;

121:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
122:   if (iascii) PetscSpaceTensorView_Ascii(sp, viewer);
123:   return 0;
124: }

126: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
127: {
128:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
129:   PetscInt           Nc, Nv, Ns;
130:   PetscBool          uniform = PETSC_TRUE;
131:   PetscInt           deg, maxDeg;
132:   PetscInt           Ncprod;

134:   if (tens->setupCalled) return 0;
135:   PetscSpaceGetNumVariables(sp, &Nv);
136:   PetscSpaceGetNumComponents(sp, &Nc);
137:   PetscSpaceTensorGetNumSubspaces(sp, &Ns);
138:   if (Ns == PETSC_DEFAULT) {
139:     Ns = Nv;
140:     PetscSpaceTensorSetNumSubspaces(sp, Ns);
141:   }
142:   if (!Ns) {
143:     SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
144:   } else {
145:     PetscSpace s0;

148:     PetscSpaceTensorGetSubspace(sp, 0, &s0);
149:     for (PetscInt i = 1; i < Ns; i++) {
150:       PetscSpace si;

152:       PetscSpaceTensorGetSubspace(sp, i, &si);
153:       if (si != s0) {
154:         uniform = PETSC_FALSE;
155:         break;
156:       }
157:     }
158:     if (uniform) {
159:       PetscInt Nvs = Nv / Ns;
160:       PetscInt Ncs;

163:       Ncs = (PetscInt)(PetscPowReal((PetscReal)Nc, 1. / Ns));
165:       if (!s0) PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0);
166:       else PetscObjectReference((PetscObject)s0);
167:       PetscSpaceSetUp(s0);
168:       for (PetscInt i = 0; i < Ns; i++) PetscSpaceTensorSetSubspace(sp, i, s0);
169:       PetscSpaceDestroy(&s0);
170:       Ncprod = PetscPowInt(Ncs, Ns);
171:     } else {
172:       PetscInt Nvsum = 0;

174:       Ncprod = 1;
175:       for (PetscInt i = 0; i < Ns; i++) {
176:         PetscInt   Nvs, Ncs;
177:         PetscSpace si;

179:         PetscSpaceTensorGetSubspace(sp, i, &si);
180:         if (!si) PetscSpaceTensorCreateSubspace(sp, 1, 1, &si);
181:         else PetscObjectReference((PetscObject)si);
182:         PetscSpaceSetUp(si);
183:         PetscSpaceTensorSetSubspace(sp, i, si);
184:         PetscSpaceGetNumVariables(si, &Nvs);
185:         PetscSpaceGetNumComponents(si, &Ncs);
186:         PetscSpaceDestroy(&si);

188:         Nvsum += Nvs;
189:         Ncprod *= Ncs;
190:       }

194:     }
195:     if (Ncprod != Nc) {
196:       PetscInt    Ncopies = Nc / Ncprod;
197:       PetscInt    Nv      = sp->Nv;
198:       const char *prefix;
199:       const char *name;
200:       char        subname[PETSC_MAX_PATH_LEN];
201:       PetscSpace  subsp;

203:       PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);
204:       PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
205:       PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);
206:       PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");
207:       if (((PetscObject)sp)->name) {
208:         PetscObjectGetName((PetscObject)sp, &name);
209:         PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name);
210:         PetscObjectSetName((PetscObject)subsp, subname);
211:       } else PetscObjectSetName((PetscObject)subsp, "sum component");
212:       PetscSpaceSetType(subsp, PETSCSPACETENSOR);
213:       PetscSpaceSetNumVariables(subsp, Nv);
214:       PetscSpaceSetNumComponents(subsp, Ncprod);
215:       PetscSpaceTensorSetNumSubspaces(subsp, Ns);
216:       for (PetscInt i = 0; i < Ns; i++) {
217:         PetscSpace ssp;

219:         PetscSpaceTensorGetSubspace(sp, i, &ssp);
220:         PetscSpaceTensorSetSubspace(subsp, i, ssp);
221:       }
222:       PetscSpaceSetUp(subsp);
223:       PetscSpaceSetType(sp, PETSCSPACESUM);
224:       PetscSpaceSumSetNumSubspaces(sp, Ncopies);
225:       for (PetscInt i = 0; i < Ncopies; i++) PetscSpaceSumSetSubspace(sp, i, subsp);
226:       PetscSpaceDestroy(&subsp);
227:       PetscSpaceSetUp(sp);
228:       return 0;
229:     }
230:   }
231:   deg    = PETSC_MAX_INT;
232:   maxDeg = 0;
233:   for (PetscInt i = 0; i < Ns; i++) {
234:     PetscSpace si;
235:     PetscInt   iDeg, iMaxDeg;

237:     PetscSpaceTensorGetSubspace(sp, i, &si);
238:     PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);
239:     deg = PetscMin(deg, iDeg);
240:     maxDeg += iMaxDeg;
241:   }
242:   sp->degree        = deg;
243:   sp->maxDegree     = maxDeg;
244:   tens->uniform     = uniform;
245:   tens->setupCalled = PETSC_TRUE;
246:   return 0;
247: }

249: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
250: {
251:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
252:   PetscInt           Ns, i;

254:   Ns = tens->numTensSpaces;
255:   if (tens->heightsubspaces) {
256:     PetscInt d;

258:     /* sp->Nv is the spatial dimension, so it is equal to the number
259:      * of subspaces on higher co-dimension points */
260:     for (d = 0; d < sp->Nv; ++d) PetscSpaceDestroy(&tens->heightsubspaces[d]);
261:   }
262:   PetscFree(tens->heightsubspaces);
263:   for (i = 0; i < Ns; i++) PetscSpaceDestroy(&tens->tensspaces[i]);
264:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL);
265:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL);
266:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL);
267:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL);
268:   PetscFree(tens->tensspaces);
269:   PetscFree(tens);
270:   return 0;
271: }

273: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
274: {
275:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
276:   PetscInt           i, Ns, d;

278:   PetscSpaceSetUp(sp);
279:   Ns = tens->numTensSpaces;
280:   d  = 1;
281:   for (i = 0; i < Ns; i++) {
282:     PetscInt id;

284:     PetscSpaceGetDimension(tens->tensspaces[i], &id);
285:     d *= id;
286:   }
287:   *dim = d;
288:   return 0;
289: }

291: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
292: {
293:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
294:   DM                 dm   = sp->dm;
295:   PetscInt           Nc   = sp->Nc;
296:   PetscInt           Nv   = sp->Nv;
297:   PetscInt           Ns;
298:   PetscReal         *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
299:   PetscInt           pdim;

301:   if (!tens->setupCalled) {
302:     PetscSpaceSetUp(sp);
303:     PetscSpaceEvaluate(sp, npoints, points, B, D, H);
304:     return 0;
305:   }
306:   Ns = tens->numTensSpaces;
307:   PetscSpaceGetDimension(sp, &pdim);
308:   DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints);
309:   if (B || D || H) DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB);
310:   if (D || H) DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD);
311:   if (H) DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH);
312:   if (B) {
313:     for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.;
314:   }
315:   if (D) {
316:     for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.;
317:   }
318:   if (H) {
319:     for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.;
320:   }
321:   for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
322:     PetscInt sNv, sNc, spdim;
323:     PetscInt vskip, cskip;

325:     PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv);
326:     PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc);
327:     PetscSpaceGetDimension(tens->tensspaces[s], &spdim);
330:     vskip = pdim / (vstep * spdim);
331:     cskip = Nc / (cstep * sNc);
332:     for (PetscInt p = 0; p < npoints; ++p) {
333:       for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i];
334:     }
335:     PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH);
336:     if (B) {
337:       for (PetscInt k = 0; k < vskip; k++) {
338:         for (PetscInt si = 0; si < spdim; si++) {
339:           for (PetscInt j = 0; j < vstep; j++) {
340:             PetscInt i = (k * spdim + si) * vstep + j;

342:             for (PetscInt l = 0; l < cskip; l++) {
343:               for (PetscInt sc = 0; sc < sNc; sc++) {
344:                 for (PetscInt m = 0; m < cstep; m++) {
345:                   PetscInt c = (l * sNc + sc) * cstep + m;

347:                   for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
348:                 }
349:               }
350:             }
351:           }
352:         }
353:       }
354:     }
355:     if (D) {
356:       for (PetscInt k = 0; k < vskip; k++) {
357:         for (PetscInt si = 0; si < spdim; si++) {
358:           for (PetscInt j = 0; j < vstep; j++) {
359:             PetscInt i = (k * spdim + si) * vstep + j;

361:             for (PetscInt l = 0; l < cskip; l++) {
362:               for (PetscInt sc = 0; sc < sNc; sc++) {
363:                 for (PetscInt m = 0; m < cstep; m++) {
364:                   PetscInt c = (l * sNc + sc) * cstep + m;

366:                   for (PetscInt der = 0; der < Nv; der++) {
367:                     if (der >= d && der < d + sNv) {
368:                       for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
369:                     } else {
370:                       for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc];
371:                     }
372:                   }
373:                 }
374:               }
375:             }
376:           }
377:         }
378:       }
379:     }
380:     if (H) {
381:       for (PetscInt k = 0; k < vskip; k++) {
382:         for (PetscInt si = 0; si < spdim; si++) {
383:           for (PetscInt j = 0; j < vstep; j++) {
384:             PetscInt i = (k * spdim + si) * vstep + j;

386:             for (PetscInt l = 0; l < cskip; l++) {
387:               for (PetscInt sc = 0; sc < sNc; sc++) {
388:                 for (PetscInt m = 0; m < cstep; m++) {
389:                   PetscInt c = (l * sNc + sc) * cstep + m;

391:                   for (PetscInt der = 0; der < Nv; der++) {
392:                     for (PetscInt der2 = 0; der2 < Nv; der2++) {
393:                       if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
394:                         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];
395:                       } else if (der >= d && der < d + sNv) {
396:                         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];
397:                       } else if (der2 >= d && der2 < d + sNv) {
398:                         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];
399:                       } else {
400:                         for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
401:                       }
402:                     }
403:                   }
404:                 }
405:               }
406:             }
407:           }
408:         }
409:       }
410:     }
411:     d += sNv;
412:     vstep *= spdim;
413:     cstep *= sNc;
414:   }
415:   if (H) DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH);
416:   if (D || H) DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD);
417:   if (B || D || H) DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB);
418:   DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints);
419:   return 0;
420: }

422: /*@
423:   PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product

425:   Input Parameters:
426: + sp  - the function space object
427: - numTensSpaces - the number of spaces

429:   Level: intermediate

431: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
432: @*/
433: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
434: {
436:   PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
437:   return 0;
438: }

440: /*@
441:   PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product

443:   Input Parameter:
444: . sp  - the function space object

446:   Output Parameter:
447: . numTensSpaces - the number of spaces

449:   Level: intermediate

451: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
452: @*/
453: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
454: {
457:   PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
458:   return 0;
459: }

461: /*@
462:   PetscSpaceTensorSetSubspace - Set a space in the tensor product

464:   Input Parameters:
465: + sp    - the function space object
466: . s     - The space number
467: - subsp - the number of spaces

469:   Level: intermediate

471: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
472: @*/
473: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
474: {
477:   PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
478:   return 0;
479: }

481: /*@
482:   PetscSpaceTensorGetSubspace - Get a space in the tensor product

484:   Input Parameters:
485: + sp - the function space object
486: - s  - The space number

488:   Output Parameter:
489: . subsp - the PetscSpace

491:   Level: intermediate

493: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
494: @*/
495: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
496: {
499:   PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
500:   return 0;
501: }

503: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
504: {
505:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
506:   PetscInt           Ns;

509:   Ns = tens->numTensSpaces;
510:   if (numTensSpaces == Ns) return 0;
511:   if (Ns >= 0) {
512:     PetscInt s;

514:     for (s = 0; s < Ns; s++) PetscSpaceDestroy(&tens->tensspaces[s]);
515:     PetscFree(tens->tensspaces);
516:   }
517:   Ns = tens->numTensSpaces = numTensSpaces;
518:   PetscCalloc1(Ns, &tens->tensspaces);
519:   return 0;
520: }

522: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
523: {
524:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;

526:   *numTensSpaces = tens->numTensSpaces;
527:   return 0;
528: }

530: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
531: {
532:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
533:   PetscInt           Ns;

536:   Ns = tens->numTensSpaces;
539:   PetscObjectReference((PetscObject)subspace);
540:   PetscSpaceDestroy(&tens->tensspaces[s]);
541:   tens->tensspaces[s] = subspace;
542:   return 0;
543: }

545: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
546: {
547:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
548:   PetscInt           Nc, dim, order, i;
549:   PetscSpace         bsp;

551:   PetscSpaceGetNumVariables(sp, &dim);
553:   PetscSpaceGetNumComponents(sp, &Nc);
554:   PetscSpaceGetDegree(sp, &order, NULL);
556:   if (!tens->heightsubspaces) PetscCalloc1(dim, &tens->heightsubspaces);
557:   if (height <= dim) {
558:     if (!tens->heightsubspaces[height - 1]) {
559:       PetscSpace  sub;
560:       const char *name;

562:       PetscSpaceTensorGetSubspace(sp, 0, &bsp);
563:       PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub);
564:       PetscObjectGetName((PetscObject)sp, &name);
565:       PetscObjectSetName((PetscObject)sub, name);
566:       PetscSpaceSetType(sub, PETSCSPACETENSOR);
567:       PetscSpaceSetNumComponents(sub, Nc);
568:       PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
569:       PetscSpaceSetNumVariables(sub, dim - height);
570:       PetscSpaceTensorSetNumSubspaces(sub, dim - height);
571:       for (i = 0; i < dim - height; i++) PetscSpaceTensorSetSubspace(sub, i, bsp);
572:       PetscSpaceSetUp(sub);
573:       tens->heightsubspaces[height - 1] = sub;
574:     }
575:     *subsp = tens->heightsubspaces[height - 1];
576:   } else {
577:     *subsp = NULL;
578:   }
579:   return 0;
580: }

582: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
583: {
584:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
585:   PetscInt           Ns;

587:   Ns = tens->numTensSpaces;
590:   *subspace = tens->tensspaces[s];
591:   return 0;
592: }

594: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
595: {
596:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Tensor;
597:   sp->ops->setup             = PetscSpaceSetUp_Tensor;
598:   sp->ops->view              = PetscSpaceView_Tensor;
599:   sp->ops->destroy           = PetscSpaceDestroy_Tensor;
600:   sp->ops->getdimension      = PetscSpaceGetDimension_Tensor;
601:   sp->ops->evaluate          = PetscSpaceEvaluate_Tensor;
602:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
603:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor);
604:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor);
605:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor);
606:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor);
607:   return 0;
608: }

610: /*MC
611:   PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space.
612:                      A tensor product is created of the components of the subspaces as well.

614:   Level: intermediate

616: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
617: M*/

619: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
620: {
621:   PetscSpace_Tensor *tens;

624:   PetscNew(&tens);
625:   sp->data = tens;

627:   tens->numTensSpaces = PETSC_DEFAULT;

629:   PetscSpaceInitialize_Tensor(sp);
630:   return 0;
631: }