Actual source code: spacetensor.c

      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++) H[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: /*@
    394:   PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product

    396:   Input Parameters:
    397: + sp  - the function space object
    398: - numTensSpaces - the number of spaces

    400:   Level: intermediate

    402: .seealso: PetscSpaceTensorGetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
    403: @*/
    404: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
    405: {

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

    414: /*@
    415:   PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product

    417:   Input Parameter:
    418: . sp  - the function space object

    420:   Output Parameter:
    421: . numTensSpaces - the number of spaces

    423:   Level: intermediate

    425: .seealso: PetscSpaceTensorSetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
    426: @*/
    427: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
    428: {

    434:   PetscTryMethod(sp,"PetscSpaceTensorGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numTensSpaces));
    435:   return(0);
    436: }

    438: /*@
    439:   PetscSpaceTensorSetSubspace - Set a space in the tensor product

    441:   Input Parameters:
    442: + sp    - the function space object
    443: . s     - The space number
    444: - subsp - the number of spaces

    446:   Level: intermediate

    448: .seealso: PetscSpaceTensorGetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
    449: @*/
    450: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
    451: {

    457:   PetscTryMethod(sp,"PetscSpaceTensorSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));
    458:   return(0);
    459: }

    461: /*@
    462:   PetscSpaceTensorGetSubspace - Get a space in the tensor product

    464:   Input Parameters:
    465: + sp - the function space object
    466: - s  - The space number

    468:   Output Parameter:
    469: . subsp - the PetscSpace

    471:   Level: intermediate

    473: .seealso: PetscSpaceTensorSetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
    474: @*/
    475: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
    476: {

    482:   PetscTryMethod(sp,"PetscSpaceTensorGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));
    483:   return(0);
    484: }

    486: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
    487: {
    488:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
    489:   PetscInt           Ns;
    490:   PetscErrorCode     ierr;

    493:   if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change number of subspaces after setup called\n");
    494:   Ns = tens->numTensSpaces;
    495:   if (numTensSpaces == Ns) return(0);
    496:   if (Ns >= 0) {
    497:     PetscInt s;

    499:     for (s = 0; s < Ns; s++) {PetscSpaceDestroy(&tens->tensspaces[s]);}
    500:     PetscFree(tens->tensspaces);
    501:   }
    502:   Ns = tens->numTensSpaces = numTensSpaces;
    503:   PetscCalloc1(Ns, &tens->tensspaces);
    504:   return(0);
    505: }

    507: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
    508: {
    509:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;

    512:   *numTensSpaces = tens->numTensSpaces;
    513:   return(0);
    514: }

    516: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
    517: {
    518:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
    519:   PetscInt           Ns;
    520:   PetscErrorCode     ierr;

    523:   if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change subspace after setup called\n");
    524:   Ns = tens->numTensSpaces;
    525:   if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n");
    526:   if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
    527:   PetscObjectReference((PetscObject)subspace);
    528:   PetscSpaceDestroy(&tens->tensspaces[s]);
    529:   tens->tensspaces[s] = subspace;
    530:   return(0);
    531: }

    533: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
    534: {
    535:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
    536:   PetscInt         Nc, dim, order, i;
    537:   PetscSpace       bsp;

    541:   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");
    542:   PetscSpaceGetNumComponents(sp, &Nc);
    543:   PetscSpaceGetNumVariables(sp, &dim);
    544:   PetscSpaceGetDegree(sp, &order, NULL);
    545:   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);
    546:   if (!tens->heightsubspaces) {PetscCalloc1(dim, &tens->heightsubspaces);}
    547:   if (height <= dim) {
    548:     if (!tens->heightsubspaces[height-1]) {
    549:       PetscSpace  sub;
    550:       const char *name;

    552:       PetscSpaceTensorGetSubspace(sp, 0, &bsp);
    553:       PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
    554:       PetscObjectGetName((PetscObject) sp,  &name);
    555:       PetscObjectSetName((PetscObject) sub,  name);
    556:       PetscSpaceSetType(sub, PETSCSPACETENSOR);
    557:       PetscSpaceSetNumComponents(sub, Nc);
    558:       PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
    559:       PetscSpaceSetNumVariables(sub, dim-height);
    560:       PetscSpaceTensorSetNumSubspaces(sub, dim-height);
    561:       for (i = 0; i < dim - height; i++) {
    562:         PetscSpaceTensorSetSubspace(sub, i, bsp);
    563:       }
    564:       PetscSpaceSetUp(sub);
    565:       tens->heightsubspaces[height-1] = sub;
    566:     }
    567:     *subsp = tens->heightsubspaces[height-1];
    568:   } else {
    569:     *subsp = NULL;
    570:   }
    571:   return(0);
    572: }

    574: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
    575: {
    576:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
    577:   PetscInt           Ns;

    580:   Ns = tens->numTensSpaces;
    581:   if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n");
    582:   if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
    583:   *subspace = tens->tensspaces[s];
    584:   return(0);
    585: }

    587: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
    588: {

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

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

    611:   Level: intermediate

    613: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
    614: M*/

    616: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
    617: {
    618:   PetscSpace_Tensor *tens;
    619:   PetscErrorCode   ierr;

    623:   PetscNewLog(sp,&tens);
    624:   sp->data = tens;

    626:   tens->numTensSpaces = PETSC_DEFAULT;

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