1 const char help[] = "Test PETSCFEVECTOR";
2
3 #include <petscfe.h>
4
PetscFEVectorTest(PetscFE orig_fe,PetscInt n_copies,PetscBool interleave_basis,PetscBool interleave_components)5 static PetscErrorCode PetscFEVectorTest(PetscFE orig_fe, PetscInt n_copies, PetscBool interleave_basis, PetscBool interleave_components)
6 {
7 PetscFE vec_fe, dup_fe;
8 PetscQuadrature quad;
9 PetscInt num_points;
10 const PetscReal *points;
11 PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)orig_fe));
12 PetscTabulation orig_T, vec_T, dup_T;
13 PetscSpace space;
14 PetscInt Nb, vNb, vNb_s, vNb_d, Nc, vNc, cdim;
15 PetscDualSpace dual_space, dup_dual_space;
16 PetscBool ib_s, ic_s, ib_d, ic_d;
17
18 PetscFunctionBegin;
19 PetscCall(PetscFEGetQuadrature(orig_fe, &quad));
20 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &num_points, &points, NULL));
21 PetscCall(PetscFECreateVector(orig_fe, n_copies, interleave_basis, interleave_components, &vec_fe));
22 PetscCall(PetscFEGetBasisSpace(vec_fe, &space));
23 PetscCall(PetscFEGetDualSpace(vec_fe, &dual_space));
24 PetscCall(PetscObjectSetName((PetscObject)vec_fe, "vector fe"));
25 PetscCall(PetscObjectSetName((PetscObject)space, "vector basis space"));
26 PetscCall(PetscObjectSetName((PetscObject)dual_space, "vector dual space"));
27 PetscCall(PetscFEView(vec_fe, viewer));
28 PetscCall(PetscFECreateTabulation(orig_fe, 1, num_points, points, 1, &orig_T));
29 PetscCall(PetscFECreateTabulation(vec_fe, 1, num_points, points, 1, &vec_T));
30 PetscCall(PetscFEGetDimension(orig_fe, &Nb));
31 PetscCall(PetscFEGetDimension(vec_fe, &vNb));
32 PetscCall(PetscFEGetNumComponents(orig_fe, &Nc));
33 PetscCall(PetscFEGetNumComponents(vec_fe, &vNc));
34 PetscCall(PetscFEGetSpatialDimension(orig_fe, &cdim));
35 {
36 PetscInt *pre_image;
37 PetscInt c_stride = interleave_components ? n_copies : 1;
38 PetscInt c_incr = interleave_components ? 1 : Nc;
39
40 PetscCall(PetscMalloc1(vNb, &pre_image));
41 for (PetscInt e = 0; e < vNb; e++) pre_image[e] = -1;
42 for (PetscInt copy = 0, coffset = 0; copy < n_copies; copy++, coffset += c_incr) {
43 for (PetscInt b = 0; b < Nb; b++) {
44 for (PetscInt e = 0; e < vNb; e++) {
45 PetscReal err = 0.0;
46
47 for (PetscInt k = 0; k <= orig_T->K; k++) {
48 const PetscReal *s_Tk = orig_T->T[k];
49 const PetscReal *v_Tk = vec_T->T[k];
50 PetscInt dblock = PetscPowInt(cdim, k);
51
52 for (PetscInt p = 0; p < num_points; p++) {
53 const PetscReal *s_Tp = &s_Tk[(p * Nb + b) * Nc * dblock];
54 const PetscReal *v_Tp = &v_Tk[(p * vNb + e) * vNc * dblock];
55 for (PetscInt c = 0; c < Nc; c++) {
56 PetscInt vc = coffset + c * c_stride;
57 const PetscReal *s_Tc = &s_Tp[c * dblock];
58 const PetscReal *v_Tc = &v_Tp[vc * dblock];
59 for (PetscInt d = 0; d < PetscPowInt(cdim, k); d++) err = PetscMax(err, PetscAbsReal(s_Tc[d] - v_Tc[d]));
60 }
61 }
62 }
63 if (err < PETSC_SMALL) {
64 PetscCheck(pre_image[e] == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Original basis %" PetscInt_FMT " and %" PetscInt_FMT " both match to vector basis %" PetscInt_FMT, pre_image[e], b, e);
65 pre_image[e] = b;
66 }
67 }
68 }
69 }
70 for (PetscInt e = 0; e < vNb; e++) PetscCheck(pre_image[e] >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No original basis matched to %" PetscInt_FMT, e);
71 PetscCall(PetscViewerASCIIPrintf(viewer, "Vector basis to original basis:"));
72 for (PetscInt e = 0; e < vNb; e++) {
73 if (!(e % 16)) PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
74 PetscCall(PetscViewerASCIIPrintf(viewer, " %3" PetscInt_FMT, pre_image[e]));
75 }
76 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
77 PetscCall(PetscFree(pre_image));
78 }
79 PetscCall(PetscSpaceSumGetInterleave(space, &ib_s, &ic_s));
80 PetscCall(PetscDualSpaceSumGetInterleave(dual_space, &ib_d, &ic_d));
81 PetscCheck(ib_s == interleave_basis, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interleave basis of space does not match");
82 PetscCheck(ic_s == interleave_components, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interleave components of space does not match");
83 PetscCheck(ib_d == interleave_basis, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interleave basis of dual space does not match");
84 PetscCheck(ic_d == interleave_components, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interleave components of dual space does not match");
85 PetscCall(PetscSpaceGetDimension(space, &vNb_s));
86 PetscCall(PetscDualSpaceGetDimension(dual_space, &vNb_d));
87 PetscCheck(vNb_s == vNb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dimension of space does not match");
88 PetscCheck(vNb_s == vNb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dimension of dual space does not match");
89 PetscCall(PetscObjectReference((PetscObject)space));
90 PetscCall(PetscDualSpaceDuplicate(dual_space, &dup_dual_space)); // not necessary just testing interface
91 PetscCall(PetscDualSpaceSetUp(dup_dual_space));
92 PetscCall(PetscFECreateFromSpaces(space, dup_dual_space, NULL, NULL, &dup_fe));
93 PetscCall(PetscFECreateTabulation(dup_fe, 1, num_points, points, 1, &dup_T));
94 {
95 PetscReal err = 0.0;
96
97 for (PetscInt k = 0; k <= vec_T->K; k++) {
98 PetscInt dblock = PetscPowInt(cdim, k);
99 PetscInt size = num_points * vNb * vNc * dblock;
100 for (PetscInt i = 0; i < size; i++) err = PetscMax(err, PetscAbsReal(vec_T->T[k][i] - dup_T->T[k][i]));
101 }
102 PetscCheck(err < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error between direct tabulation and indirect tabulation: %g", (double)err);
103 }
104 PetscCall(PetscTabulationDestroy(&dup_T));
105 PetscCall(PetscTabulationDestroy(&vec_T));
106 PetscCall(PetscTabulationDestroy(&orig_T));
107 PetscCall(PetscFEDestroy(&dup_fe));
108 PetscCall(PetscFEDestroy(&vec_fe));
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
main(int argc,char ** argv)112 int main(int argc, char **argv)
113 {
114 PetscFE scalar, vector;
115 PetscViewer viewer;
116
117 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
118 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, 3, 1, PETSC_TRUE, 3, PETSC_DETERMINE, &scalar));
119 viewer = PETSC_VIEWER_STDOUT_SELF;
120 PetscCall(PetscObjectSetName((PetscObject)scalar, "base FE (scalar)"));
121 PetscCall(PetscFEView(scalar, viewer));
122 PetscCall(PetscViewerASCIIPushTab(viewer));
123 for (PetscInt n_copies = 1; n_copies <= 3; n_copies++) {
124 PetscCall(PetscFEVectorTest(scalar, n_copies, PETSC_FALSE, PETSC_FALSE));
125 PetscCall(PetscFEVectorTest(scalar, n_copies, PETSC_FALSE, PETSC_TRUE));
126 PetscCall(PetscFEVectorTest(scalar, n_copies, PETSC_TRUE, PETSC_FALSE));
127 PetscCall(PetscFEVectorTest(scalar, n_copies, PETSC_TRUE, PETSC_TRUE));
128 }
129 PetscCall(PetscViewerASCIIPopTab(viewer));
130 PetscCall(PetscFEDestroy(&scalar));
131 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, 3, 3, PETSC_TRUE, 3, PETSC_DETERMINE, &vector));
132 PetscCall(PetscObjectSetName((PetscObject)vector, "base FE (vector)"));
133 PetscCall(PetscFEView(vector, viewer));
134 PetscCall(PetscViewerASCIIPushTab(viewer));
135 for (PetscInt n_copies = 1; n_copies <= 3; n_copies++) {
136 PetscCall(PetscFEVectorTest(vector, n_copies, PETSC_FALSE, PETSC_FALSE));
137 PetscCall(PetscFEVectorTest(vector, n_copies, PETSC_FALSE, PETSC_TRUE));
138 PetscCall(PetscFEVectorTest(vector, n_copies, PETSC_TRUE, PETSC_FALSE));
139 PetscCall(PetscFEVectorTest(vector, n_copies, PETSC_TRUE, PETSC_TRUE));
140 }
141 PetscCall(PetscViewerASCIIPopTab(viewer));
142 PetscCall(PetscFEDestroy(&vector));
143 PetscCall(PetscFinalize());
144 return 0;
145 }
146
147 /*TEST
148
149 test:
150 suffix: 0
151
152 TEST*/
153