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