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\n", (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\n", (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\n", (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*/