Actual source code: ex16.c
1: const char help[] = "Test PETSCFEVECTOR";
3: #include <petscfe.h>
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;
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;
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;
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);
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 %d and %d both match to vector basis %d", (int)pre_image[e], (int)b, (int)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 %d", (int)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, " %3d", (int)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;
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: }
112: int main(int argc, char **argv)
113: {
114: PetscFE scalar, vector;
115: PetscViewer viewer;
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: }
147: /*TEST
149: test:
150: suffix: 0
152: TEST*/