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:   PetscFunctionBegin;
 18:   PetscCall(PetscSpaceCreate(comm, &sp));
 19:   PetscCall(PetscObjectSetName((PetscObject)sp, "ptrimmed"));
 20:   PetscCall(PetscSpaceSetType(sp, PETSCSPACEPTRIMMED));
 21:   PetscCall(PetscSpaceSetNumVariables(sp, dim));
 22:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nf));
 23:   PetscCall(PetscSpaceSetNumComponents(sp, Nf * nCopies));
 24:   PetscCall(PetscSpaceSetDegree(sp, degree, PETSC_DETERMINE));
 25:   PetscCall(PetscSpacePTrimmedSetFormDegree(sp, formDegree));
 26:   PetscCall(PetscSpaceSetUp(sp));
 27:   PetscCall(PetscSpaceView(sp, NULL));

 29:   PetscCall(PetscDTPTrimmedSize(dim, formDegree == 0 ? degree : degree + 1, PetscAbsInt(formDegree), &Nbexp));
 30:   Nbexp *= nCopies;
 31:   PetscCall(PetscSpaceGetDimension(sp, &Nb));
 32:   PetscCheck(Nb == Nbexp, comm, PETSC_ERR_PLIB, "Space dimension mismatch, %" PetscInt_FMT " != %" PetscInt_FMT, Nbexp, Nb);

 34:   maxDexp = (PetscAbsInt(formDegree) == dim || formDegree == 0) ? degree : degree + 1;
 35:   PetscCall(PetscSpaceGetDegree(sp, &d, &maxD));
 36:   PetscCheck(degree == d, comm, PETSC_ERR_PLIB, "Space degree mismatch, %" PetscInt_FMT " != %" PetscInt_FMT, degree, d);
 37:   PetscCheck(maxDexp == maxD, comm, PETSC_ERR_PLIB, "Space max degree mismatch, %" PetscInt_FMT " != %" PetscInt_FMT, maxDexp, maxD);

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

 42:   Bsize = npoints * Nb * Nf * nCopies;
 43:   Dsize = dim * Bsize;
 44:   Hsize = dim * Dsize;
 45:   PetscCall(PetscMalloc3(Bsize, &B, Dsize, &D, Hsize, &H));
 46:   PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
 47:   for (PetscInt i = 0; i < Bsize; i++) PetscCheck(!PetscIsInfOrNanReal(B[i]), comm, PETSC_ERR_PLIB, "Bad value B[%" PetscInt_FMT "]", i);
 48:   for (PetscInt i = 0; i < Dsize; i++) PetscCheck(!PetscIsInfOrNanReal(D[i]), comm, PETSC_ERR_PLIB, "Bad value D[%" PetscInt_FMT "]", i);
 49:   for (PetscInt i = 0; i < Hsize; i++) PetscCheck(!PetscIsInfOrNanReal(H[i]), comm, PETSC_ERR_PLIB, "Bad value H[%" PetscInt_FMT "]", i);
 50:   PetscCall(PetscFree3(B, D, H));
 51:   PetscCall(PetscQuadratureDestroy(&quad));
 52:   PetscCall(PetscSpaceDestroy(&sp));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

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

 72: /*TEST

 74:   test:

 76: TEST*/