1 const char help[] = "Test basic creation and evaluation of PETSCSPACEPTRIMMED";
2
3 #include <petscfe.h>
4
test(PetscInt dim,PetscInt formDegree,PetscInt degree,PetscInt nCopies)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
main(int argc,char * argv[])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