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 PetscDTAltVStar(PetscInt, PetscInt, PetscInt, const PetscReal *, PetscReal *); 71 72 PETSC_STATIC_INLINE PetscErrorCode PetscDTBinomial(PetscInt n, PetscInt k, PetscInt *binomial) 73 { 74 PetscFunctionBeginHot; 75 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); 76 if (n <= 3) { 77 PetscInt binomLookup[4][4] = {{1, 0, 0, 0}, {1, 1, 0, 0}, {1, 2, 1, 0}, {1, 3, 3, 1}}; 78 79 *binomial = binomLookup[n][k]; 80 } else { 81 PetscReal binom = 1; 82 PetscInt i; 83 84 k = PetscMin(k, n - k); 85 for (i = 0; i < k; i++) binom = (binom * (n - i)) / (i + 1); 86 *binomial = binom; 87 } 88 PetscFunctionReturn(0); 89 } 90 91 PETSC_STATIC_INLINE PetscErrorCode PetscDTEnumPerm(PetscInt n, PetscInt k, PetscInt *work, PetscInt *perm, PetscBool *isOdd) 92 { 93 PetscInt odd = 0; 94 PetscInt i; 95 PetscInt *w = &work[n - 2]; 96 97 PetscFunctionBeginHot; 98 for (i = 2; i <= n; i++) { 99 *(w--) = k % i; 100 k /= i; 101 } 102 for (i = 0; i < n; i++) perm[i] = i; 103 for (i = 0; i < n - 1; i++) { 104 PetscInt s = work[i]; 105 PetscInt swap = perm[i]; 106 107 perm[i] = perm[i + s]; 108 perm[i + s] = swap; 109 odd ^= (!!s); 110 } 111 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 112 PetscFunctionReturn(0); 113 } 114 115 PETSC_STATIC_INLINE PetscErrorCode PetscDTEnumSubset(PetscInt n, PetscInt k, PetscInt j, PetscInt *subset) 116 { 117 PetscInt Nk, i, l; 118 PetscErrorCode ierr; 119 120 PetscFunctionBeginHot; 121 ierr = PetscDTBinomial(n, k, &Nk);CHKERRQ(ierr); 122 for (i = 0, l = 0; i < n && l < k; i++) { 123 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 124 PetscInt Nminusk = Nk - Nminuskminus; 125 126 if (j < Nminuskminus) { 127 subset[l++] = i; 128 Nk = Nminuskminus; 129 } else { 130 j -= Nminuskminus; 131 Nk = Nminusk; 132 } 133 } 134 PetscFunctionReturn(0); 135 } 136 137 PETSC_STATIC_INLINE PetscErrorCode PetscDTSubsetIndex(PetscInt n, PetscInt k, const PetscInt *subset, PetscInt *index) 138 { 139 PetscInt i, j = 0, l, Nk; 140 PetscErrorCode ierr; 141 142 PetscFunctionBeginHot; 143 ierr = PetscDTBinomial(n, k, &Nk);CHKERRQ(ierr); 144 for (i = 0, l = 0; i < n && l < k; i++) { 145 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 146 PetscInt Nminusk = Nk - Nminuskminus; 147 148 if (subset[l] == i) { 149 l++; 150 Nk = Nminuskminus; 151 } else { 152 j += Nminuskminus; 153 Nk = Nminusk; 154 } 155 } 156 *index = j; 157 PetscFunctionReturn(0); 158 } 159 160 161 PETSC_STATIC_INLINE PetscErrorCode PetscDTEnumSplit(PetscInt n, PetscInt k, PetscInt j, PetscInt *subset, PetscBool *isOdd) 162 { 163 PetscInt i, l, m, *subcomp, Nk; 164 PetscInt odd; 165 PetscErrorCode ierr; 166 167 PetscFunctionBeginHot; 168 ierr = PetscDTBinomial(n, k, &Nk);CHKERRQ(ierr); 169 odd = 0; 170 subcomp = &subset[k]; 171 for (i = 0, l = 0, m = 0; i < n && l < k; i++) { 172 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 173 PetscInt Nminusk = Nk - Nminuskminus; 174 175 if (j < Nminuskminus) { 176 subset[l++] = i; 177 Nk = Nminuskminus; 178 } else { 179 subcomp[m++] = i; 180 j -= Nminuskminus; 181 odd ^= ((k - l) & 1); 182 Nk = Nminusk; 183 } 184 } 185 for (; i < n; i++) { 186 subcomp[m++] = i; 187 } 188 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 189 PetscFunctionReturn(0); 190 } 191 192 #endif 193