1*ccb4e88aSToby Isaac const char help[] = "Test basic creation and evaluation of PETSCSPACEPTRIMMED"; 2*ccb4e88aSToby Isaac 3*ccb4e88aSToby Isaac #include <petscfe.h> 4*ccb4e88aSToby Isaac 5*ccb4e88aSToby Isaac static PetscErrorCode test(PetscInt dim, PetscInt formDegree, PetscInt degree, PetscInt nCopies) 6*ccb4e88aSToby Isaac { 7*ccb4e88aSToby Isaac MPI_Comm comm = PETSC_COMM_SELF; 8*ccb4e88aSToby Isaac PetscSpace sp; 9*ccb4e88aSToby Isaac PetscInt Nf, Nb; 10*ccb4e88aSToby Isaac PetscInt maxDexp, maxD, d; 11*ccb4e88aSToby Isaac PetscInt Nbexp, Bsize, Dsize, Hsize; 12*ccb4e88aSToby Isaac PetscReal *B, *D, *H; 13*ccb4e88aSToby Isaac PetscQuadrature quad; 14*ccb4e88aSToby Isaac PetscInt npoints; 15*ccb4e88aSToby Isaac const PetscReal *points; 16*ccb4e88aSToby Isaac PetscErrorCode ierr; 17*ccb4e88aSToby Isaac 18*ccb4e88aSToby Isaac PetscFunctionBegin; 19*ccb4e88aSToby Isaac ierr = PetscSpaceCreate(comm, &sp);CHKERRQ(ierr); 20*ccb4e88aSToby Isaac ierr = PetscObjectSetName((PetscObject)sp, "ptrimmed");CHKERRQ(ierr); 21*ccb4e88aSToby Isaac ierr = PetscSpaceSetType(sp, PETSCSPACEPTRIMMED);CHKERRQ(ierr); 22*ccb4e88aSToby Isaac ierr = PetscSpaceSetNumVariables(sp, dim);CHKERRQ(ierr); 23*ccb4e88aSToby Isaac ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nf);CHKERRQ(ierr); 24*ccb4e88aSToby Isaac ierr = PetscSpaceSetNumComponents(sp, Nf * nCopies);CHKERRQ(ierr); 25*ccb4e88aSToby Isaac ierr = PetscSpaceSetDegree(sp, degree, PETSC_DETERMINE);CHKERRQ(ierr); 26*ccb4e88aSToby Isaac ierr = PetscSpacePTrimmedSetFormDegree(sp, formDegree);CHKERRQ(ierr); 27*ccb4e88aSToby Isaac ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 28*ccb4e88aSToby Isaac ierr = PetscSpaceView(sp, NULL);CHKERRQ(ierr); 29*ccb4e88aSToby Isaac 30*ccb4e88aSToby Isaac ierr = PetscDTPTrimmedSize(dim, formDegree == 0 ? degree : degree + 1, PetscAbsInt(formDegree), &Nbexp);CHKERRQ(ierr); 31*ccb4e88aSToby Isaac Nbexp *= nCopies; 32*ccb4e88aSToby Isaac ierr = PetscSpaceGetDimension(sp, &Nb);CHKERRQ(ierr); 33*ccb4e88aSToby Isaac if (Nb != Nbexp) SETERRQ2(comm, PETSC_ERR_PLIB, "Space dimension mismatch, %D != %D", Nbexp, Nb); 34*ccb4e88aSToby Isaac 35*ccb4e88aSToby Isaac maxDexp = (PetscAbsInt(formDegree) == dim || formDegree == 0) ? degree : degree + 1; 36*ccb4e88aSToby Isaac ierr = PetscSpaceGetDegree(sp, &d, &maxD);CHKERRQ(ierr); 37*ccb4e88aSToby Isaac if (degree != d) SETERRQ2(comm, PETSC_ERR_PLIB, "Space degree mismatch, %D != %D", degree, d); 38*ccb4e88aSToby Isaac if (maxDexp != maxD) SETERRQ2(comm, PETSC_ERR_PLIB, "Space max degree mismatch, %D != %D", maxDexp, maxD); 39*ccb4e88aSToby Isaac 40*ccb4e88aSToby Isaac ierr = PetscDTStroudConicalQuadrature(dim, 1, maxD + 1, -1., 1., &quad);CHKERRQ(ierr); 41*ccb4e88aSToby Isaac ierr = PetscQuadratureGetData(quad, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); 42*ccb4e88aSToby Isaac 43*ccb4e88aSToby Isaac Bsize = npoints * Nb * Nf * nCopies; 44*ccb4e88aSToby Isaac Dsize = dim * Bsize; 45*ccb4e88aSToby Isaac Hsize = dim * Dsize; 46*ccb4e88aSToby Isaac ierr = PetscMalloc3(Bsize, &B, Dsize, &D, Hsize, &H);CHKERRQ(ierr); 47*ccb4e88aSToby Isaac ierr = PetscSpaceEvaluate(sp, npoints, points, B, D, H);CHKERRQ(ierr); 48*ccb4e88aSToby Isaac for (PetscInt i = 0; i < Bsize; i++) { 49*ccb4e88aSToby Isaac if (PetscIsInfOrNanReal(B[i])) SETERRQ1(comm, PETSC_ERR_PLIB, "Bad value B[%D]", i); 50*ccb4e88aSToby Isaac } 51*ccb4e88aSToby Isaac for (PetscInt i = 0; i < Dsize; i++) { 52*ccb4e88aSToby Isaac if (PetscIsInfOrNanReal(D[i])) SETERRQ1(comm, PETSC_ERR_PLIB, "Bad value D[%D]", i); 53*ccb4e88aSToby Isaac } 54*ccb4e88aSToby Isaac for (PetscInt i = 0; i < Hsize; i++) { 55*ccb4e88aSToby Isaac if (PetscIsInfOrNanReal(H[i])) SETERRQ1(comm, PETSC_ERR_PLIB, "Bad value H[%H]", i); 56*ccb4e88aSToby Isaac } 57*ccb4e88aSToby Isaac ierr = PetscFree3(B, D, H);CHKERRQ(ierr); 58*ccb4e88aSToby Isaac ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 59*ccb4e88aSToby Isaac ierr = PetscSpaceDestroy(&sp);CHKERRQ(ierr); 60*ccb4e88aSToby Isaac PetscFunctionReturn(0); 61*ccb4e88aSToby Isaac } 62*ccb4e88aSToby Isaac 63*ccb4e88aSToby Isaac int main(int argc, char *argv[]) 64*ccb4e88aSToby Isaac { 65*ccb4e88aSToby Isaac PetscErrorCode ierr; 66*ccb4e88aSToby Isaac 67*ccb4e88aSToby Isaac ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 68*ccb4e88aSToby Isaac for (PetscInt dim = 0; dim <= 3; dim++) { 69*ccb4e88aSToby Isaac for (PetscInt formDegree = -dim; formDegree <= dim; formDegree++) { 70*ccb4e88aSToby Isaac for (PetscInt degree = 0; degree <= 4; degree++) { 71*ccb4e88aSToby Isaac if (formDegree == 0 && degree == 0) continue; 72*ccb4e88aSToby Isaac for (PetscInt nCopies = 1; nCopies <= PetscMax(2,dim); nCopies++) { 73*ccb4e88aSToby Isaac ierr = test(dim, formDegree, degree, nCopies);CHKERRQ(ierr); 74*ccb4e88aSToby Isaac } 75*ccb4e88aSToby Isaac } 76*ccb4e88aSToby Isaac } 77*ccb4e88aSToby Isaac } 78*ccb4e88aSToby Isaac ierr = PetscFinalize(); 79*ccb4e88aSToby Isaac return ierr; 80*ccb4e88aSToby Isaac } 81*ccb4e88aSToby Isaac 82*ccb4e88aSToby Isaac /*TEST 83*ccb4e88aSToby Isaac 84*ccb4e88aSToby Isaac test: 85*ccb4e88aSToby Isaac 86*ccb4e88aSToby Isaac TEST*/ 87