Actual source code: dspacesimple.c

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

      4: static PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
      5: {
      6:   PetscDualSpace_Simple *s  = (PetscDualSpace_Simple *) sp->data;
      7:   DM                     dm = sp->dm;
      8:   PetscInt               dim, pStart, pEnd;
      9:   PetscSection           section;
     10:   PetscErrorCode         ierr;

     13:   DMGetDimension(dm, &dim);
     14:   DMPlexGetChart(dm, &pStart, &pEnd);
     15:   PetscSectionCreate(PETSC_COMM_SELF, &section);
     16:   PetscSectionSetChart(section, pStart, pEnd);
     17:   PetscSectionSetDof(section, pStart, s->dim);
     18:   PetscSectionSetUp(section);
     19:   sp->pointSection = section;
     20:   return(0);
     21: }

     23: static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
     24: {
     25:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
     26:   PetscErrorCode         ierr;

     29:   PetscFree(s->numDof);
     30:   PetscFree(s);
     31:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
     32:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
     33:   return(0);
     34: }

     36: static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew)
     37: {
     38:   PetscInt       dim, d;

     42:   PetscDualSpaceGetDimension(sp, &dim);
     43:   PetscDualSpaceSimpleSetDimension(spNew, dim);
     44:   for (d = 0; d < dim; ++d) {
     45:     PetscQuadrature q;

     47:     PetscDualSpaceGetFunctional(sp, d, &q);
     48:     PetscDualSpaceSimpleSetFunctional(spNew, d, q);
     49:   }
     50:   return(0);
     51: }

     53: static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
     54: {
     56:   return(0);
     57: }

     59: static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
     60: {
     61:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
     62:   DM                     dm;
     63:   PetscInt               spatialDim, f;
     64:   PetscErrorCode         ierr;

     67:   for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
     68:   PetscFree(sp->functional);
     69:   s->dim = dim;
     70:   PetscCalloc1(s->dim, &sp->functional);
     71:   PetscFree(s->numDof);
     72:   PetscDualSpaceGetDM(sp, &dm);
     73:   DMGetCoordinateDim(dm, &spatialDim);
     74:   PetscCalloc1(spatialDim+1, &s->numDof);
     75:   s->numDof[spatialDim] = dim;
     76:   return(0);
     77: }

     79: static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
     80: {
     81:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
     82:   PetscReal             *weights;
     83:   PetscInt               Nc, c, Nq, p;
     84:   PetscErrorCode         ierr;

     87:   if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
     88:   PetscQuadratureDuplicate(q, &sp->functional[f]);
     89:   /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
     90:   PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **) &weights);
     91:   for (c = 0; c < Nc; ++c) {
     92:     PetscReal vol = 0.0;

     94:     for (p = 0; p < Nq; ++p) vol += weights[p*Nc+c];
     95:     for (p = 0; p < Nq; ++p) weights[p*Nc+c] /= (vol == 0.0 ? 1.0 : vol);
     96:   }
     97:   return(0);
     98: }

    100: /*@
    101:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

    103:   Logically Collective on sp

    105:   Input Parameters:
    106: + sp  - the PetscDualSpace
    107: - dim - the basis dimension

    109:   Level: intermediate

    111: .seealso: PetscDualSpaceSimpleSetFunctional()
    112: @*/
    113: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
    114: {

    120:   if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up");
    121:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
    122:   return(0);
    123: }

    125: /*@
    126:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

    128:   Not Collective

    130:   Input Parameters:
    131: + sp  - the PetscDualSpace
    132: . f - the basis index
    133: - q - the basis functional

    135:   Level: intermediate

    137:   Note: The quadrature will be reweighted so that it has unit volume.

    139: .seealso: PetscDualSpaceSimpleSetDimension()
    140: @*/
    141: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
    142: {

    147:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
    148:   return(0);
    149: }

    151: static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
    152: {
    154:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Simple;
    155:   sp->ops->setup                = PetscDualSpaceSetUp_Simple;
    156:   sp->ops->view                 = NULL;
    157:   sp->ops->destroy              = PetscDualSpaceDestroy_Simple;
    158:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Simple;
    159:   sp->ops->createheightsubspace = NULL;
    160:   sp->ops->createpointsubspace  = NULL;
    161:   sp->ops->getsymmetries        = NULL;
    162:   sp->ops->apply                = PetscDualSpaceApplyDefault;
    163:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
    164:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
    165:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
    166:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
    167:   return(0);
    168: }

    170: /*MC
    171:   PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals

    173:   Level: intermediate

    175: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
    176: M*/

    178: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
    179: {
    180:   PetscDualSpace_Simple *s;
    181:   PetscErrorCode         ierr;

    185:   PetscNewLog(sp,&s);
    186:   sp->data = s;

    188:   s->dim    = 0;
    189:   s->numDof = NULL;

    191:   PetscDualSpaceInitialize_Simple(sp);
    192:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
    193:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
    194:   return(0);
    195: }