Actual source code: ex1.c
1: const char help[] = "Test basic creation and evaluation of PETSCSPACEPTRIMMED";
3: #include <petscfe.h>
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;
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));
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);
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);
39: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, maxD + 1, -1., 1., &quad));
40: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &npoints, &points, NULL));
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: }
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: }
72: /*TEST
74: test:
76: TEST*/