Actual source code: ex1.c

  1: const char help[] = "Test basic creation and evaluation of PETSCSPACEPTRIMMED";

  3: #include <petscfe.h>

  5: static PetscErrorCode test(PetscInt dim, PetscInt formDegree, PetscInt degree, PetscInt nCopies)
  6: {
  7:   MPI_Comm         comm = PETSC_COMM_SELF;
  8:   PetscSpace       sp;
  9:   PetscInt         Nf, Nb;
 10:   PetscInt         maxDexp, maxD, d;
 11:   PetscInt         Nbexp, Bsize, Dsize, Hsize;
 12:   PetscReal       *B, *D, *H;
 13:   PetscQuadrature  quad;
 14:   PetscInt         npoints;
 15:   const PetscReal *points;

 17:   PetscSpaceCreate(comm, &sp);
 18:   PetscObjectSetName((PetscObject)sp, "ptrimmed");
 19:   PetscSpaceSetType(sp, PETSCSPACEPTRIMMED);
 20:   PetscSpaceSetNumVariables(sp, dim);
 21:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nf);
 22:   PetscSpaceSetNumComponents(sp, Nf * nCopies);
 23:   PetscSpaceSetDegree(sp, degree, PETSC_DETERMINE);
 24:   PetscSpacePTrimmedSetFormDegree(sp, formDegree);
 25:   PetscSpaceSetUp(sp);
 26:   PetscSpaceView(sp, NULL);

 28:   PetscDTPTrimmedSize(dim, formDegree == 0 ? degree : degree + 1, PetscAbsInt(formDegree), &Nbexp);
 29:   Nbexp *= nCopies;
 30:   PetscSpaceGetDimension(sp, &Nb);

 33:   maxDexp = (PetscAbsInt(formDegree) == dim || formDegree == 0) ? degree : degree + 1;
 34:   PetscSpaceGetDegree(sp, &d, &maxD);

 38:   PetscDTStroudConicalQuadrature(dim, 1, maxD + 1, -1., 1., &quad);
 39:   PetscQuadratureGetData(quad, NULL, NULL, &npoints, &points, NULL);

 41:   Bsize = npoints * Nb * Nf * nCopies;
 42:   Dsize = dim * Bsize;
 43:   Hsize = dim * Dsize;
 44:   PetscMalloc3(Bsize, &B, Dsize, &D, Hsize, &H);
 45:   PetscSpaceEvaluate(sp, npoints, points, B, D, H);
 49:   PetscFree3(B, D, H);
 50:   PetscQuadratureDestroy(&quad);
 51:   PetscSpaceDestroy(&sp);
 52:   return 0;
 53: }

 55: int main(int argc, char *argv[])
 56: {
 58:   PetscInitialize(&argc, &argv, NULL, help);
 59:   for (PetscInt dim = 0; dim <= 3; dim++) {
 60:     for (PetscInt formDegree = -dim; formDegree <= dim; formDegree++) {
 61:       for (PetscInt degree = 0; degree <= 4; degree++) {
 62:         if (formDegree == 0 && degree == 0) continue;
 63:         for (PetscInt nCopies = 1; nCopies <= PetscMax(2, dim); nCopies++) test(dim, formDegree, degree, nCopies);
 64:       }
 65:     }
 66:   }
 67:   PetscFinalize();
 68:   return 0;
 69: }

 71: /*TEST

 73:   test:

 75: TEST*/