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++) { 48 PetscCheckFalse(PetscIsInfOrNanReal(B[i]),comm, PETSC_ERR_PLIB, "Bad value B[%" PetscInt_FMT "]", i); 49 } 50 for (PetscInt i = 0; i < Dsize; i++) { 51 PetscCheckFalse(PetscIsInfOrNanReal(D[i]),comm, PETSC_ERR_PLIB, "Bad value D[%" PetscInt_FMT "]", i); 52 } 53 for (PetscInt i = 0; i < Hsize; i++) { 54 PetscCheckFalse(PetscIsInfOrNanReal(H[i]),comm, PETSC_ERR_PLIB, "Bad value H[%" PetscInt_FMT "]", i); 55 } 56 PetscCall(PetscFree3(B, D, H)); 57 PetscCall(PetscQuadratureDestroy(&quad)); 58 PetscCall(PetscSpaceDestroy(&sp)); 59 PetscFunctionReturn(0); 60 } 61 62 int main(int argc, char *argv[]) 63 { 64 65 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 66 for (PetscInt dim = 0; dim <= 3; dim++) { 67 for (PetscInt formDegree = -dim; formDegree <= dim; formDegree++) { 68 for (PetscInt degree = 0; degree <= 4; degree++) { 69 if (formDegree == 0 && degree == 0) continue; 70 for (PetscInt nCopies = 1; nCopies <= PetscMax(2,dim); nCopies++) { 71 PetscCall(test(dim, formDegree, degree, nCopies)); 72 } 73 } 74 } 75 } 76 PetscCall(PetscFinalize()); 77 return 0; 78 } 79 80 /*TEST 81 82 test: 83 84 TEST*/ 85