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:   PetscSpace       sp;
  8:   PetscInt         Nf, Nb;
  9:   PetscInt         maxDexp, maxD, d;
 10:   PetscInt         Nbexp, Bsize, Dsize, Hsize;
 11:   PetscReal       *B, *D, *H;
 12:   PetscQuadrature  quad;
 13:   PetscInt         npoints;
 14:   const PetscReal *points;

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

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

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

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

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

 55: int main(int argc, char *argv[])
 56: {
 57:   PetscFunctionBeginUser;
 58:   PetscCall(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++) PetscCall(test(dim, formDegree, degree, nCopies));
 64:       }
 65:     }
 66:   }
 67:   PetscCall(PetscFinalize());
 68:   return 0;
 69: }

 71: /*TEST

 73:   test:

 75: TEST*/