1 const char help[] = "Test basic creation and evaluation of PETSCSPACEPTRIMMED"; 2 3 #include <petscfe.h> 4 5 static PetscErrorCode test(PetscInt dim, PetscInt formDegree, PetscInt degree, PetscInt nCopies) { 6 MPI_Comm comm = PETSC_COMM_SELF; 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; 15 16 PetscFunctionBegin; 17 PetscCall(PetscSpaceCreate(comm, &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)); 27 28 PetscCall(PetscDTPTrimmedSize(dim, formDegree == 0 ? degree : degree + 1, PetscAbsInt(formDegree), &Nbexp)); 29 Nbexp *= nCopies; 30 PetscCall(PetscSpaceGetDimension(sp, &Nb)); 31 PetscCheck(Nb == Nbexp, comm, PETSC_ERR_PLIB, "Space dimension mismatch, %" PetscInt_FMT " != %" PetscInt_FMT, Nbexp, Nb); 32 33 maxDexp = (PetscAbsInt(formDegree) == dim || formDegree == 0) ? degree : degree + 1; 34 PetscCall(PetscSpaceGetDegree(sp, &d, &maxD)); 35 PetscCheck(degree == d, comm, PETSC_ERR_PLIB, "Space degree mismatch, %" PetscInt_FMT " != %" PetscInt_FMT, degree, d); 36 PetscCheck(maxDexp == maxD, comm, PETSC_ERR_PLIB, "Space max degree mismatch, %" PetscInt_FMT " != %" PetscInt_FMT, maxDexp, maxD); 37 38 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, maxD + 1, -1., 1., &quad)); 39 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &npoints, &points, NULL)); 40 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]), comm, PETSC_ERR_PLIB, "Bad value B[%" PetscInt_FMT "]", i); } 47 for (PetscInt i = 0; i < Dsize; i++) { PetscCheck(!PetscIsInfOrNanReal(D[i]), comm, PETSC_ERR_PLIB, "Bad value D[%" PetscInt_FMT "]", i); } 48 for (PetscInt i = 0; i < Hsize; i++) { PetscCheck(!PetscIsInfOrNanReal(H[i]), comm, 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(0); 53 } 54 55 int main(int argc, char *argv[]) { 56 PetscFunctionBeginUser; 57 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 58 for (PetscInt dim = 0; dim <= 3; dim++) { 59 for (PetscInt formDegree = -dim; formDegree <= dim; formDegree++) { 60 for (PetscInt degree = 0; degree <= 4; degree++) { 61 if (formDegree == 0 && degree == 0) continue; 62 for (PetscInt nCopies = 1; nCopies <= PetscMax(2, dim); nCopies++) { PetscCall(test(dim, formDegree, degree, nCopies)); } 63 } 64 } 65 } 66 PetscCall(PetscFinalize()); 67 return 0; 68 } 69 70 /*TEST 71 72 test: 73 74 TEST*/ 75