Actual source code: spacetensor.c

petsc-master 2019-08-18
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>

  3: static PetscErrorCode PetscSpaceTensorCreateSubspace(PetscSpace space, PetscInt Nvs, PetscSpace *subspace)
  4: {
  5:   PetscInt    degree;
  6:   const char *prefix;

 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, 1);
 16:   PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE);
 17:   PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix);
 18:   PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "subspace_");
 19:   return(0);
 20: }

 22: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
 23: {
 24:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
 25:   PetscInt           Ns, Nc, i, Nv, deg;
 26:   PetscBool          uniform = PETSC_TRUE;
 27:   PetscErrorCode     ierr;

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

 38:     PetscSpaceTensorGetSubspace(sp, 0, &s0);
 39:     for (i = 1; i < Ns; i++) {
 40:       PetscSpace si;

 42:       PetscSpaceTensorGetSubspace(sp, i, &si);
 43:       if (si != s0) {uniform = PETSC_FALSE; break;}
 44:     }
 45:   }
 46:   Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv,1) : Ns;
 47:   PetscOptionsHead(PetscOptionsObject,"PetscSpace tensor options");
 48:   PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL,0);
 49:   PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL);
 50:   PetscOptionsTail();
 51:   if (Ns < 0 || (Nv > 0 && Ns == 0)) SETERRQ1(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a tensor space made up of %D spaces\n",Ns);
 52:   if (Nv > 0 && Ns > Nv) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a tensor space with %D subspaces over %D variables\n", Ns, Nv);
 53:   if (Ns != tens->numTensSpaces) {PetscSpaceTensorSetNumSubspaces(sp, Ns);}
 54:   if (uniform) {
 55:     PetscInt   Nvs = Nv / Ns;
 56:     PetscSpace subspace;

 58:     if (Nv % Ns) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONG,"Cannot use %D uniform subspaces for %D variable space\n", Ns, Nv);
 59:     PetscSpaceTensorGetSubspace(sp, 0, &subspace);
 60:     if (!subspace) {PetscSpaceTensorCreateSubspace(sp, Nvs, &subspace);}
 61:     else           {PetscObjectReference((PetscObject)subspace);}
 62:     PetscSpaceSetFromOptions(subspace);
 63:     for (i = 0; i < Ns; i++) {PetscSpaceTensorSetSubspace(sp, i, subspace);}
 64:     PetscSpaceDestroy(&subspace);
 65:   } else {
 66:     for (i = 0; i < Ns; i++) {
 67:       PetscSpace subspace;

 69:       PetscSpaceTensorGetSubspace(sp, i, &subspace);
 70:       if (!subspace) {
 71:         char tprefix[128];

 73:         PetscSpaceTensorCreateSubspace(sp, 1, &subspace);
 74:         PetscSNPrintf(tprefix, 128, "%d_",(int)i);
 75:         PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix);
 76:       } else {
 77:         PetscObjectReference((PetscObject)subspace);
 78:       }
 79:       PetscSpaceSetFromOptions(subspace);
 80:       PetscSpaceTensorSetSubspace(sp, i, subspace);
 81:       PetscSpaceDestroy(&subspace);
 82:     }
 83:   }
 84:   return(0);
 85: }

 87: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
 88: {
 89:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
 90:   PetscBool          uniform = PETSC_TRUE;
 91:   PetscInt           Ns = tens->numTensSpaces, i, n;
 92:   PetscErrorCode     ierr;

 95:   for (i = 1; i < Ns; i++) {
 96:     if (tens->tensspaces[i] != tens->tensspaces[0]) {uniform = PETSC_FALSE; break;}
 97:   }
 98:   if (uniform) {PetscViewerASCIIPrintf(v, "Tensor space of %D subspaces (all identical)\n", Ns);}
 99:   else         {PetscViewerASCIIPrintf(v, "Tensor space of %D subspaces\n", Ns);}
100:   n = uniform ? 1 : Ns;
101:   for (i = 0; i < n; i++) {
102:     PetscViewerASCIIPushTab(v);
103:     PetscSpaceView(tens->tensspaces[i], v);
104:     PetscViewerASCIIPopTab(v);
105:   }
106:   return(0);
107: }

109: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
110: {
111:   PetscBool      iascii;

115:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
116:   if (iascii) {PetscSpaceTensorView_Ascii(sp, viewer);}
117:   return(0);
118: }

120: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
121: {
122:   PetscSpace_Tensor *tens    = (PetscSpace_Tensor *) sp->data;
123:   PetscInt           Nv, Ns, i;
124:   PetscBool          uniform = PETSC_TRUE;
125:   PetscInt           deg, maxDeg;
126:   PetscErrorCode     ierr;

129:   if (tens->setupCalled) return(0);
130:   PetscSpaceGetNumVariables(sp, &Nv);
131:   PetscSpaceTensorGetNumSubspaces(sp, &Ns);
132:   if (Ns == PETSC_DEFAULT) {
133:     Ns = Nv;
134:     PetscSpaceTensorSetNumSubspaces(sp, Ns);
135:   }
136:   if (!Ns) {
137:     if (Nv) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
138:   } else {
139:     PetscSpace s0;

141:     if (Nv > 0 && Ns > Nv) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a tensor space with %D subspaces over %D variables\n", Ns, Nv);
142:     PetscSpaceTensorGetSubspace(sp, 0, &s0);
143:     for (i = 1; i < Ns; i++) {
144:       PetscSpace si;

146:       PetscSpaceTensorGetSubspace(sp, i, &si);
147:       if (si != s0) {uniform = PETSC_FALSE; break;}
148:     }
149:     if (uniform) {
150:       PetscInt   Nvs = Nv / Ns;

152:       if (Nv % Ns) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONG,"Cannot use %D uniform subspaces for %D variable space\n", Ns, Nv);
153:       if (!s0) {PetscSpaceTensorCreateSubspace(sp, Nvs, &s0);}
154:       else     {PetscObjectReference((PetscObject) s0);}
155:       PetscSpaceSetUp(s0);
156:       for (i = 0; i < Ns; i++) {PetscSpaceTensorSetSubspace(sp, i, s0);}
157:       PetscSpaceDestroy(&s0);
158:     } else {
159:       for (i = 0 ; i < Ns; i++) {
160:         PetscSpace si;

162:         PetscSpaceTensorGetSubspace(sp, i, &si);
163:         if (!si) {PetscSpaceTensorCreateSubspace(sp, 1, &si);}
164:         else     {PetscObjectReference((PetscObject) si);}
165:         PetscSpaceSetUp(si);
166:         PetscSpaceTensorSetSubspace(sp, i, si);
167:         PetscSpaceDestroy(&si);
168:       }
169:     }
170:   }
171:   deg = PETSC_MAX_INT;
172:   maxDeg = 0;
173:   for (i = 0; i < Ns; i++) {
174:     PetscSpace si;
175:     PetscInt   iDeg, iMaxDeg;

177:     PetscSpaceTensorGetSubspace(sp, i, &si);
178:     PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);
179:     deg    = PetscMin(deg, iDeg);
180:     maxDeg += iMaxDeg;
181:   }
182:   sp->degree    = deg;
183:   sp->maxDegree = maxDeg;
184:   tens->uniform = uniform;
185:   tens->setupCalled = PETSC_TRUE;
186:   return(0);
187: }

189: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
190: {
191:   PetscSpace_Tensor *tens    = (PetscSpace_Tensor *) sp->data;
192:   PetscInt           Ns, i;
193:   PetscErrorCode     ierr;

196:   Ns = tens->numTensSpaces;
197:   if (tens->heightsubspaces) {
198:     PetscInt d;

200:     /* sp->Nv is the spatial dimension, so it is equal to the number
201:      * of subspaces on higher co-dimension points */
202:     for (d = 0; d < sp->Nv; ++d) {
203:       PetscSpaceDestroy(&tens->heightsubspaces[d]);
204:     }
205:   }
206:   PetscFree(tens->heightsubspaces);
207:   for (i = 0; i < Ns; i++) {PetscSpaceDestroy(&tens->tensspaces[i]);}
208:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", NULL);
209:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", NULL);
210:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", NULL);
211:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", NULL);
212:   PetscFree(tens->tensspaces);
213:   PetscFree(tens);
214:   return(0);
215: }

217: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
218: {
219:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
220:   PetscInt           i, Ns, Nc, d;
221:   PetscErrorCode     ierr;

224:   PetscSpaceSetUp(sp);
225:   Ns = tens->numTensSpaces;
226:   Nc = sp->Nc;
227:   d  = 1;
228:   for (i = 0; i < Ns; i++) {
229:     PetscInt id;

231:     PetscSpaceGetDimension(tens->tensspaces[i], &id);
232:     d *= id;
233:   }
234:   d *= Nc;
235:   *dim = d;
236:   return(0);
237: }

239: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
240: {
241:   PetscSpace_Tensor *tens  = (PetscSpace_Tensor *) sp->data;
242:   DM               dm      = sp->dm;
243:   PetscInt         Nc      = sp->Nc;
244:   PetscInt         Nv      = sp->Nv;
245:   PetscInt         Ns;
246:   PetscReal       *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
247:   PetscInt         c, pdim, d, e, der, der2, i, l, si, p, s, step;
248:   PetscErrorCode   ierr;

251:   if (!tens->setupCalled) {PetscSpaceSetUp(sp);}
252:   Ns = tens->numTensSpaces;
253:   PetscSpaceGetDimension(sp,&pdim);
254:   pdim /= Nc;
255:   DMGetWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);
256:   if (B || D || H) {DMGetWorkArray(dm, npoints*pdim,       MPIU_REAL, &sB);}
257:   if (D || H)      {DMGetWorkArray(dm, npoints*pdim*Nv,    MPIU_REAL, &sD);}
258:   if (H)           {DMGetWorkArray(dm, npoints*pdim*Nv*Nv, MPIU_REAL, &sH);}
259:   if (B) {
260:     for (i = 0; i < npoints*pdim*Nc*Nc; i++) B[i] = 0.;
261:     for (i = 0; i < npoints*pdim; i++) B[i * Nc*Nc] = 1.;
262:   }
263:   if (D) {
264:     for (i = 0; i < npoints*pdim*Nc*Nc*Nv; i++) D[i] = 0.;
265:     for (i = 0; i < npoints*pdim; i++) {
266:       for (l = 0; l < Nv; l++) {
267:         D[i * Nc*Nc*Nv + l] = 1.;
268:       }
269:     }
270:   }
271:   if (H) {
272:     for (i = 0; i < npoints*pdim*Nc*Nc*Nv*Nv; i++) D[i] = 0.;
273:     for (i = 0; i < npoints*pdim; i++) {
274:       for (l = 0; l < Nv*Nv; l++) {
275:         H[i * Nc*Nc*Nv*Nv + l] = 1.;
276:       }
277:     }
278:   }
279:   for (s = 0, d = 0, step = 1; s < Ns; s++) {
280:     PetscInt sNv, spdim;
281:     PetscInt skip, j, k;

283:     PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv);
284:     PetscSpaceGetDimension(tens->tensspaces[s], &spdim);
285:     if ((pdim % step) || (pdim % spdim))  SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %d, Ns %D, pdim %D, s %D, step %D, spdim %D", Nv, Ns, pdim, s, step, spdim);
286:     skip = pdim / (step * spdim);
287:     for (p = 0; p < npoints; ++p) {
288:       for (i = 0; i < sNv; i++) {
289:         lpoints[p * sNv + i] = points[p*Nv + d + i];
290:       }
291:     }
292:     PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH);
293:     if (B) {
294:       for (p = 0; p < npoints; p++) {
295:         for (k = 0; k < skip; k++) {
296:           for (si = 0; si < spdim; si++) {
297:             for (j = 0; j < step; j++) {
298:               i = (k * spdim + si) * step + j;
299:               B[(pdim * p + i) * Nc * Nc] *= sB[spdim * p + si];
300:             }
301:           }
302:         }
303:       }
304:     }
305:     if (D) {
306:       for (p = 0; p < npoints; p++) {
307:         for (k = 0; k < skip; k++) {
308:           for (si = 0; si < spdim; si++) {
309:             for (j = 0; j < step; j++) {
310:               i = (k * spdim + si) * step + j;
311:               for (der = 0; der < Nv; der++) {
312:                 if (der >= d && der < d + sNv) {
313:                   D[(pdim * p + i) * Nc*Nc*Nv + der] *= sD[(spdim * p + si) * sNv + der - d];
314:                 } else {
315:                   D[(pdim * p + i) * Nc*Nc*Nv + der] *= sB[spdim * p + si];
316:                 }
317:               }
318:             }
319:           }
320:         }
321:       }
322:     }
323:     if (H) {
324:       for (p = 0; p < npoints; p++) {
325:         for (k = 0; k < skip; k++) {
326:           for (si = 0; si < spdim; si++) {
327:             for (j = 0; j < step; j++) {
328:               i = (k * spdim + si) * step + j;
329:               for (der = 0; der < Nv; der++) {
330:                 for (der2 = 0; der2 < Nv; der2++) {
331:                   if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
332:                     H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sH[((spdim * p + si) * sNv + der - d) * sNv + der2 - d];
333:                   } else if (der >= d && der < d + sNv) {
334:                     H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sD[(spdim * p + si) * sNv + der - d];
335:                   } else if (der2 >= d && der2 < d + sNv) {
336:                     H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sD[(spdim * p + si) * sNv + der2 - d];
337:                   } else {
338:                     H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sB[spdim * p + si];
339:                   }
340:                 }
341:               }
342:             }
343:           }
344:         }
345:       }
346:     }
347:     d += sNv;
348:     step *= spdim;
349:   }
350:   if (B && Nc > 1) {
351:     /* Make direct sum basis for multicomponent space */
352:     for (p = 0; p < npoints; ++p) {
353:       for (i = 0; i < pdim; ++i) {
354:         for (c = 1; c < Nc; ++c) {
355:           B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
356:         }
357:       }
358:     }
359:   }
360:   if (D && Nc > 1) {
361:     /* Make direct sum basis for multicomponent space */
362:     for (p = 0; p < npoints; ++p) {
363:       for (i = 0; i < pdim; ++i) {
364:         for (c = 1; c < Nc; ++c) {
365:           for (d = 0; d < Nv; ++d) {
366:             D[((p*pdim*Nc + i*Nc + c)*Nc + c)*Nv + d] = D[(p*pdim + i)*Nc*Nc*Nv + d];
367:           }
368:         }
369:       }
370:     }
371:   }
372:   if (H && Nc > 1) {
373:     /* Make direct sum basis for multicomponent space */
374:     for (p = 0; p < npoints; ++p) {
375:       for (i = 0; i < pdim; ++i) {
376:         for (c = 1; c < Nc; ++c) {
377:           for (d = 0; d < Nv; ++d) {
378:             for (e = 0; e < Nv; ++e) {
379:               H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*Nv + d)*Nv + e] = H[((p*pdim + i)*Nc*Nc*Nv + d)*Nv + e];
380:             }
381:           }
382:         }
383:       }
384:     }
385:   }
386:   if (H)           {DMRestoreWorkArray(dm, npoints*pdim*Nv*Nv, MPIU_REAL, &sH);}
387:   if (D || H)      {DMRestoreWorkArray(dm, npoints*pdim*Nv,    MPIU_REAL, &sD);}
388:   if (B || D || H) {DMRestoreWorkArray(dm, npoints*pdim,       MPIU_REAL, &sB);}
389:   DMRestoreWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);
390:   return(0);
391: }

393: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
394: {

399:   PetscTryMethod(sp,"PetscSpaceTensorSetNumSubspaces_C",(PetscSpace,PetscInt),(sp,numTensSpaces));
400:   return(0);
401: }

403: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
404: {

410:   PetscTryMethod(sp,"PetscSpaceTensorGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numTensSpaces));
411:   return(0);
412: }

414: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
415: {

421:   PetscTryMethod(sp,"PetscSpaceTensorSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));
422:   return(0);
423: }

425: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
426: {

432:   PetscTryMethod(sp,"PetscSpaceTensorGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));
433:   return(0);
434: }

436: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
437: {
438:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
439:   PetscInt           Ns;
440:   PetscErrorCode     ierr;

443:   if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change number of subspaces after setup called\n");
444:   Ns = tens->numTensSpaces;
445:   if (numTensSpaces == Ns) return(0);
446:   if (Ns >= 0) {
447:     PetscInt s;

449:     for (s = 0; s < Ns; s++) {PetscSpaceDestroy(&tens->tensspaces[s]);}
450:     PetscFree(tens->tensspaces);
451:   }
452:   Ns = tens->numTensSpaces = numTensSpaces;
453:   PetscCalloc1(Ns, &tens->tensspaces);
454:   return(0);
455: }

457: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
458: {
459:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;

462:   *numTensSpaces = tens->numTensSpaces;
463:   return(0);
464: }

466: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
467: {
468:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
469:   PetscInt           Ns;
470:   PetscErrorCode     ierr;

473:   if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change subspace after setup called\n");
474:   Ns = tens->numTensSpaces;
475:   if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n");
476:   if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
477:   PetscObjectReference((PetscObject)subspace);
478:   PetscSpaceDestroy(&tens->tensspaces[s]);
479:   tens->tensspaces[s] = subspace;
480:   return(0);
481: }

483: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
484: {
485:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
486:   PetscInt         Nc, dim, order, i;
487:   PetscSpace       bsp;

491:   if (!tens->uniform) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Can only get a generic height subspace of a uniform tensor space: this tensor space is not uniform.\n");
492:   PetscSpaceGetNumComponents(sp, &Nc);
493:   PetscSpaceGetNumVariables(sp, &dim);
494:   PetscSpaceGetDegree(sp, &order, NULL);
495:   if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);}
496:   if (!tens->heightsubspaces) {PetscCalloc1(dim, &tens->heightsubspaces);}
497:   if (height <= dim) {
498:     if (!tens->heightsubspaces[height-1]) {
499:       PetscSpace  sub;
500:       const char *name;

502:       PetscSpaceTensorGetSubspace(sp, 0, &bsp);
503:       PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
504:       PetscObjectGetName((PetscObject) sp,  &name);
505:       PetscObjectSetName((PetscObject) sub,  name);
506:       PetscSpaceSetType(sub, PETSCSPACETENSOR);
507:       PetscSpaceSetNumComponents(sub, Nc);
508:       PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
509:       PetscSpaceSetNumVariables(sub, dim-height);
510:       PetscSpaceTensorSetNumSubspaces(sub, dim-height);
511:       for (i = 0; i < dim - height; i++) {
512:         PetscSpaceTensorSetSubspace(sub, i, bsp);
513:       }
514:       PetscSpaceSetUp(sub);
515:       tens->heightsubspaces[height-1] = sub;
516:     }
517:     *subsp = tens->heightsubspaces[height-1];
518:   } else {
519:     *subsp = NULL;
520:   }
521:   return(0);
522: }

524: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
525: {
526:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
527:   PetscInt           Ns;

530:   Ns = tens->numTensSpaces;
531:   if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n");
532:   if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
533:   *subspace = tens->tensspaces[s];
534:   return(0);
535: }

537: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
538: {

542:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Tensor;
543:   sp->ops->setup             = PetscSpaceSetUp_Tensor;
544:   sp->ops->view              = PetscSpaceView_Tensor;
545:   sp->ops->destroy           = PetscSpaceDestroy_Tensor;
546:   sp->ops->getdimension      = PetscSpaceGetDimension_Tensor;
547:   sp->ops->evaluate          = PetscSpaceEvaluate_Tensor;
548:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
549:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor);
550:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor);
551:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor);
552:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor);
553:   return(0);
554: }

556: /*MC
557:   PETSCSPACETENSOR = "tensor" - A PetscSpace object that encapsulates a tensor product space.
558:                      Subspaces are scalar spaces (num of componenents = 1), so the components
559:                      of a vector-valued tensor space are assumed to be identical.

561:   Level: intermediate

563: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
564: M*/

566: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
567: {
568:   PetscSpace_Tensor *tens;
569:   PetscErrorCode   ierr;

573:   PetscNewLog(sp,&tens);
574:   sp->data = tens;

576:   tens->numTensSpaces = PETSC_DEFAULT;

578:   PetscSpaceInitialize_Tensor(sp);
579:   return(0);
580: }