1 const char help[] = "Test PETSCFEVECTOR"; 2 3 #include <petscfe.h> 4 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 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