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 CHKERRQ(PetscSpaceCreate(comm, &sp)); 19 CHKERRQ(PetscObjectSetName((PetscObject)sp, "ptrimmed")); 20 CHKERRQ(PetscSpaceSetType(sp, PETSCSPACEPTRIMMED)); 21 CHKERRQ(PetscSpaceSetNumVariables(sp, dim)); 22 CHKERRQ(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nf)); 23 CHKERRQ(PetscSpaceSetNumComponents(sp, Nf * nCopies)); 24 CHKERRQ(PetscSpaceSetDegree(sp, degree, PETSC_DETERMINE)); 25 CHKERRQ(PetscSpacePTrimmedSetFormDegree(sp, formDegree)); 26 CHKERRQ(PetscSpaceSetUp(sp)); 27 CHKERRQ(PetscSpaceView(sp, NULL)); 28 29 CHKERRQ(PetscDTPTrimmedSize(dim, formDegree == 0 ? degree : degree + 1, PetscAbsInt(formDegree), &Nbexp)); 30 Nbexp *= nCopies; 31 CHKERRQ(PetscSpaceGetDimension(sp, &Nb)); 32 PetscCheckFalse(Nb != Nbexp,comm, PETSC_ERR_PLIB, "Space dimension mismatch, %D != %D", Nbexp, Nb); 33 34 maxDexp = (PetscAbsInt(formDegree) == dim || formDegree == 0) ? degree : degree + 1; 35 CHKERRQ(PetscSpaceGetDegree(sp, &d, &maxD)); 36 PetscCheckFalse(degree != d,comm, PETSC_ERR_PLIB, "Space degree mismatch, %D != %D", degree, d); 37 PetscCheckFalse(maxDexp != maxD,comm, PETSC_ERR_PLIB, "Space max degree mismatch, %D != %D", maxDexp, maxD); 38 39 CHKERRQ(PetscDTStroudConicalQuadrature(dim, 1, maxD + 1, -1., 1., &quad)); 40 CHKERRQ(PetscQuadratureGetData(quad, NULL, NULL, &npoints, &points, NULL)); 41 42 Bsize = npoints * Nb * Nf * nCopies; 43 Dsize = dim * Bsize; 44 Hsize = dim * Dsize; 45 CHKERRQ(PetscMalloc3(Bsize, &B, Dsize, &D, Hsize, &H)); 46 CHKERRQ(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[%D]", i); 49 } 50 for (PetscInt i = 0; i < Dsize; i++) { 51 PetscCheckFalse(PetscIsInfOrNanReal(D[i]),comm, PETSC_ERR_PLIB, "Bad value D[%D]", i); 52 } 53 for (PetscInt i = 0; i < Hsize; i++) { 54 PetscCheckFalse(PetscIsInfOrNanReal(H[i]),comm, PETSC_ERR_PLIB, "Bad value H[%H]", i); 55 } 56 CHKERRQ(PetscFree3(B, D, H)); 57 CHKERRQ(PetscQuadratureDestroy(&quad)); 58 CHKERRQ(PetscSpaceDestroy(&sp)); 59 PetscFunctionReturn(0); 60 } 61 62 int main(int argc, char *argv[]) 63 { 64 PetscErrorCode ierr; 65 66 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 67 for (PetscInt dim = 0; dim <= 3; dim++) { 68 for (PetscInt formDegree = -dim; formDegree <= dim; formDegree++) { 69 for (PetscInt degree = 0; degree <= 4; degree++) { 70 if (formDegree == 0 && degree == 0) continue; 71 for (PetscInt nCopies = 1; nCopies <= PetscMax(2,dim); nCopies++) { 72 CHKERRQ(test(dim, formDegree, degree, nCopies)); 73 } 74 } 75 } 76 } 77 ierr = PetscFinalize(); 78 return ierr; 79 } 80 81 /*TEST 82 83 test: 84 85 TEST*/ 86