1 /* 2 Common tools for constructing discretizations 3 */ 4 #if !defined(PETSCDT_H) 5 #define PETSCDT_H 6 7 #include <petscsys.h> 8 9 /*S 10 PetscQuadrature - Quadrature rule for integration. 11 12 Level: beginner 13 14 .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy() 15 S*/ 16 typedef struct _p_PetscQuadrature *PetscQuadrature; 17 18 /*E 19 PetscGaussLobattoLegendreCreateType - algorithm used to compute the Gauss-Lobatto-Legendre nodes and weights 20 21 Level: intermediate 22 23 $ PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA - compute the nodes via linear algebra 24 $ PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON - compute the nodes by solving a nonlinear equation with Newton's method 25 26 E*/ 27 typedef enum {PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON} PetscGaussLobattoLegendreCreateType; 28 29 PETSC_EXTERN PetscErrorCode PetscQuadratureCreate(MPI_Comm, PetscQuadrature *); 30 PETSC_EXTERN PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature, PetscQuadrature *); 31 PETSC_EXTERN PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature, PetscInt*); 32 PETSC_EXTERN PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature, PetscInt); 33 PETSC_EXTERN PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature, PetscInt*); 34 PETSC_EXTERN PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature, PetscInt); 35 PETSC_EXTERN PetscErrorCode PetscQuadratureGetData(PetscQuadrature, PetscInt*, PetscInt*, PetscInt*, const PetscReal *[], const PetscReal *[]); 36 PETSC_EXTERN PetscErrorCode PetscQuadratureSetData(PetscQuadrature, PetscInt, PetscInt, PetscInt, const PetscReal [], const PetscReal []); 37 PETSC_EXTERN PetscErrorCode PetscQuadratureView(PetscQuadrature, PetscViewer); 38 PETSC_EXTERN PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *); 39 40 PETSC_EXTERN PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature, PetscInt, const PetscReal[], const PetscReal[], PetscQuadrature *); 41 42 PETSC_EXTERN PetscErrorCode PetscDTLegendreEval(PetscInt,const PetscReal*,PetscInt,const PetscInt*,PetscReal*,PetscReal*,PetscReal*); 43 PETSC_EXTERN PetscErrorCode PetscDTGaussQuadrature(PetscInt,PetscReal,PetscReal,PetscReal*,PetscReal*); 44 PETSC_EXTERN PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt,PetscGaussLobattoLegendreCreateType,PetscReal*,PetscReal*); 45 PETSC_EXTERN PetscErrorCode PetscDTReconstructPoly(PetscInt,PetscInt,const PetscReal*,PetscInt,const PetscReal*,PetscReal*); 46 PETSC_EXTERN PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt,PetscInt,PetscInt,PetscReal,PetscReal,PetscQuadrature*); 47 PETSC_EXTERN PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt,PetscInt,PetscInt,PetscReal,PetscReal,PetscQuadrature*); 48 49 PETSC_EXTERN PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt, PetscInt, PetscReal, PetscReal, PetscQuadrature *); 50 PETSC_EXTERN PetscErrorCode PetscDTTanhSinhIntegrate(void (*)(PetscReal, PetscReal *), PetscReal, PetscReal, PetscInt, PetscReal *); 51 PETSC_EXTERN PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*)(PetscReal, PetscReal *), PetscReal, PetscReal, PetscInt, PetscReal *); 52 53 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt, PetscReal *, PetscReal *, const PetscReal *, PetscReal *); 54 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 55 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 56 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***, PetscReal ***); 57 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***, PetscReal ***); 58 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 59 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 60 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 61 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 62 63 PETSC_EXTERN PetscErrorCode PetscDTAltVApply(PetscInt, PetscInt, const PetscReal *, const PetscReal *, PetscReal *); 64 PETSC_EXTERN PetscErrorCode PetscDTAltVWedge(PetscInt, PetscInt, PetscInt, const PetscReal *, const PetscReal *, PetscReal *); 65 PETSC_EXTERN PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt, PetscInt, PetscInt, const PetscReal *, PetscReal *); 66 PETSC_EXTERN PetscErrorCode PetscDTAltVPullback(PetscInt, PetscInt, const PetscReal *, PetscInt, const PetscReal *, PetscReal *); 67 PETSC_EXTERN PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt, PetscInt, const PetscReal *, PetscInt, PetscReal *); 68 PETSC_EXTERN PetscErrorCode PetscDTAltVInterior(PetscInt, PetscInt, const PetscReal *, const PetscReal *, PetscReal *); 69 PETSC_EXTERN PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt, PetscInt, const PetscReal *, PetscReal *); 70 PETSC_EXTERN PetscErrorCode PetscDTAltVInteriorPattern(PetscInt, PetscInt, PetscInt (*)[3]); 71 PETSC_EXTERN PetscErrorCode PetscDTAltVStar(PetscInt, PetscInt, PetscInt, const PetscReal *, PetscReal *); 72 73 PETSC_STATIC_INLINE PetscErrorCode PetscDTBinomial(PetscInt n, PetscInt k, PetscInt *binomial) 74 { 75 PetscFunctionBeginHot; 76 if (n < 0 || k < 0 || k > n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Binomal arguments (%D %D) must be non-negative, k <= n\n", n, k); 77 if (n <= 3) { 78 PetscInt binomLookup[4][4] = {{1, 0, 0, 0}, {1, 1, 0, 0}, {1, 2, 1, 0}, {1, 3, 3, 1}}; 79 80 *binomial = binomLookup[n][k]; 81 } else { 82 PetscReal binom = 1; 83 PetscInt i; 84 85 k = PetscMin(k, n - k); 86 for (i = 0; i < k; i++) binom = (binom * (n - i)) / (i + 1); 87 *binomial = binom; 88 } 89 PetscFunctionReturn(0); 90 } 91 92 PETSC_STATIC_INLINE PetscErrorCode PetscDTEnumPerm(PetscInt n, PetscInt k, PetscInt *work, PetscInt *perm, PetscBool *isOdd) 93 { 94 PetscInt odd = 0; 95 PetscInt i; 96 PetscInt *w = &work[n - 2]; 97 98 PetscFunctionBeginHot; 99 for (i = 2; i <= n; i++) { 100 *(w--) = k % i; 101 k /= i; 102 } 103 for (i = 0; i < n; i++) perm[i] = i; 104 for (i = 0; i < n - 1; i++) { 105 PetscInt s = work[i]; 106 PetscInt swap = perm[i]; 107 108 perm[i] = perm[i + s]; 109 perm[i + s] = swap; 110 odd ^= (!!s); 111 } 112 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 113 PetscFunctionReturn(0); 114 } 115 116 PETSC_STATIC_INLINE PetscErrorCode PetscDTEnumSubset(PetscInt n, PetscInt k, PetscInt j, PetscInt *subset) 117 { 118 PetscInt Nk, i, l; 119 PetscErrorCode ierr; 120 121 PetscFunctionBeginHot; 122 ierr = PetscDTBinomial(n, k, &Nk);CHKERRQ(ierr); 123 for (i = 0, l = 0; i < n && l < k; i++) { 124 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 125 PetscInt Nminusk = Nk - Nminuskminus; 126 127 if (j < Nminuskminus) { 128 subset[l++] = i; 129 Nk = Nminuskminus; 130 } else { 131 j -= Nminuskminus; 132 Nk = Nminusk; 133 } 134 } 135 PetscFunctionReturn(0); 136 } 137 138 PETSC_STATIC_INLINE PetscErrorCode PetscDTSubsetIndex(PetscInt n, PetscInt k, const PetscInt *subset, PetscInt *index) 139 { 140 PetscInt i, j = 0, l, Nk; 141 PetscErrorCode ierr; 142 143 PetscFunctionBeginHot; 144 ierr = PetscDTBinomial(n, k, &Nk);CHKERRQ(ierr); 145 for (i = 0, l = 0; i < n && l < k; i++) { 146 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 147 PetscInt Nminusk = Nk - Nminuskminus; 148 149 if (subset[l] == i) { 150 l++; 151 Nk = Nminuskminus; 152 } else { 153 j += Nminuskminus; 154 Nk = Nminusk; 155 } 156 } 157 *index = j; 158 PetscFunctionReturn(0); 159 } 160 161 162 PETSC_STATIC_INLINE PetscErrorCode PetscDTEnumSplit(PetscInt n, PetscInt k, PetscInt j, PetscInt *subset, PetscBool *isOdd) 163 { 164 PetscInt i, l, m, *subcomp, Nk; 165 PetscInt odd; 166 PetscErrorCode ierr; 167 168 PetscFunctionBeginHot; 169 ierr = PetscDTBinomial(n, k, &Nk);CHKERRQ(ierr); 170 odd = 0; 171 subcomp = &subset[k]; 172 for (i = 0, l = 0, m = 0; i < n && l < k; i++) { 173 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 174 PetscInt Nminusk = Nk - Nminuskminus; 175 176 if (j < Nminuskminus) { 177 subset[l++] = i; 178 Nk = Nminuskminus; 179 } else { 180 subcomp[m++] = i; 181 j -= Nminuskminus; 182 odd ^= ((k - l) & 1); 183 Nk = Nminusk; 184 } 185 } 186 for (; i < n; i++) { 187 subcomp[m++] = i; 188 } 189 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 190 PetscFunctionReturn(0); 191 } 192 193 #endif 194