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 { 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; 16 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)); 28 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); 33 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); 38 39 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, maxD + 1, -1., 1., &quad)); 40 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &npoints, &points, NULL)); 41 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 } 55 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 } 71 72 /*TEST 73 74 test: 75 76 TEST*/ 77