1 /* 2 Common tools for constructing discretizations 3 */ 4 #if !defined(PETSCDT_H) 5 #define PETSCDT_H 6 7 #include <petscsys.h> 8 9 /* SUBMANSEC = DT */ 10 11 PETSC_EXTERN PetscClassId PETSCQUADRATURE_CLASSID; 12 13 /*S 14 PetscQuadrature - Quadrature rule for integration. 15 16 Level: beginner 17 18 .seealso: `PetscQuadratureCreate()`, `PetscQuadratureDestroy()` 19 S*/ 20 typedef struct _p_PetscQuadrature *PetscQuadrature; 21 22 /*E 23 PetscGaussLobattoLegendreCreateType - algorithm used to compute the Gauss-Lobatto-Legendre nodes and weights 24 25 Level: intermediate 26 27 $ `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` - compute the nodes via linear algebra 28 $ `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON` - compute the nodes by solving a nonlinear equation with Newton's method 29 30 E*/ 31 typedef enum { 32 PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, 33 PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON 34 } PetscGaussLobattoLegendreCreateType; 35 36 /*E 37 PetscDTNodeType - A description of strategies for generating nodes (both 38 quadrature nodes and nodes for Lagrange polynomials) 39 40 Level: intermediate 41 42 $ `PETSCDTNODES_DEFAULT` - Nodes chosen by PETSc 43 $ `PETSCDTNODES_GAUSSJACOBI` - Nodes at either Gauss-Jacobi or Gauss-Lobatto-Jacobi quadrature points 44 $ `PETSCDTNODES_EQUISPACED` - Nodes equispaced either including the endpoints or excluding them 45 $ `PETSCDTNODES_TANHSINH` - Nodes at Tanh-Sinh quadrature points 46 47 Note: 48 A `PetscDTNodeType` can be paired with a `PetscBool` to indicate whether 49 the nodes include endpoints or not, and in the case of `PETSCDT_GAUSSJACOBI` 50 with exponents for the weight function. 51 52 E*/ 53 typedef enum { 54 PETSCDTNODES_DEFAULT = -1, 55 PETSCDTNODES_GAUSSJACOBI, 56 PETSCDTNODES_EQUISPACED, 57 PETSCDTNODES_TANHSINH 58 } PetscDTNodeType; 59 60 PETSC_EXTERN const char *const *const PetscDTNodeTypes; 61 62 /*E 63 PetscDTSimplexQuadratureType - A description of classes of quadrature rules for simplices 64 65 Level: intermediate 66 67 $ `PETSCDTSIMPLEXQUAD_DEFAULT` - Quadrature rule chosen by PETSc 68 $ `PETSCDTSIMPLEXQUAD_CONIC` - Quadrature rules constructed as 69 conically-warped tensor products of 1D 70 Gauss-Jacobi quadrature rules. These are 71 explicitly computable in any dimension for any 72 degree, and the tensor-product structure can be 73 exploited by sum-factorization methods, but 74 they are not efficient in terms of nodes per 75 polynomial degree. 76 $ `PETSCDTSIMPLEXQUAD_MINSYM` - Quadrature rules that are fully symmetric 77 (symmetries of the simplex preserve the nodes 78 and weights) with minimal (or near minimal) 79 number of nodes. In dimensions higher than 1 80 these are not simple to compute, so lookup 81 tables are used. 82 83 .seealso: `PetscDTSimplexQuadrature()` 84 E*/ 85 typedef enum { 86 PETSCDTSIMPLEXQUAD_DEFAULT = -1, 87 PETSCDTSIMPLEXQUAD_CONIC = 0, 88 PETSCDTSIMPLEXQUAD_MINSYM 89 } PetscDTSimplexQuadratureType; 90 91 PETSC_EXTERN const char *const *const PetscDTSimplexQuadratureTypes; 92 93 PETSC_EXTERN PetscErrorCode PetscQuadratureCreate(MPI_Comm, PetscQuadrature *); 94 PETSC_EXTERN PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature, PetscQuadrature *); 95 PETSC_EXTERN PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature, PetscInt *); 96 PETSC_EXTERN PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature, PetscInt); 97 PETSC_EXTERN PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature, PetscInt *); 98 PETSC_EXTERN PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature, PetscInt); 99 PETSC_EXTERN PetscErrorCode PetscQuadratureEqual(PetscQuadrature, PetscQuadrature, PetscBool *); 100 PETSC_EXTERN PetscErrorCode PetscQuadratureGetData(PetscQuadrature, PetscInt *, PetscInt *, PetscInt *, const PetscReal *[], const PetscReal *[]); 101 PETSC_EXTERN PetscErrorCode PetscQuadratureSetData(PetscQuadrature, PetscInt, PetscInt, PetscInt, const PetscReal[], const PetscReal[]); 102 PETSC_EXTERN PetscErrorCode PetscQuadratureView(PetscQuadrature, PetscViewer); 103 PETSC_EXTERN PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *); 104 105 PETSC_EXTERN PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature, PetscQuadrature, PetscQuadrature *); 106 PETSC_EXTERN PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature, PetscInt, const PetscReal[], const PetscReal[], PetscQuadrature *); 107 108 PETSC_EXTERN PetscErrorCode PetscQuadraturePushForward(PetscQuadrature, PetscInt, const PetscReal[], const PetscReal[], const PetscReal[], PetscInt, PetscQuadrature *); 109 110 PETSC_EXTERN PetscErrorCode PetscDTLegendreEval(PetscInt, const PetscReal *, PetscInt, const PetscInt *, PetscReal *, PetscReal *, PetscReal *); 111 PETSC_EXTERN PetscErrorCode PetscDTJacobiNorm(PetscReal, PetscReal, PetscInt, PetscReal *); 112 PETSC_EXTERN PetscErrorCode PetscDTJacobiEval(PetscInt, PetscReal, PetscReal, const PetscReal *, PetscInt, const PetscInt *, PetscReal *, PetscReal *, PetscReal *); 113 PETSC_EXTERN PetscErrorCode PetscDTJacobiEvalJet(PetscReal, PetscReal, PetscInt, const PetscReal[], PetscInt, PetscInt, PetscReal[]); 114 PETSC_EXTERN PetscErrorCode PetscDTPKDEvalJet(PetscInt, PetscInt, const PetscReal[], PetscInt, PetscInt, PetscReal[]); 115 PETSC_EXTERN PetscErrorCode PetscDTPTrimmedSize(PetscInt, PetscInt, PetscInt, PetscInt *); 116 PETSC_EXTERN PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt, PetscInt, const PetscReal[], PetscInt, PetscInt, PetscInt, PetscReal[]); 117 PETSC_EXTERN PetscErrorCode PetscDTGaussQuadrature(PetscInt, PetscReal, PetscReal, PetscReal *, PetscReal *); 118 PETSC_EXTERN PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *, PetscReal *); 119 PETSC_EXTERN PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *, PetscReal *); 120 PETSC_EXTERN PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt, PetscGaussLobattoLegendreCreateType, PetscReal *, PetscReal *); 121 PETSC_EXTERN PetscErrorCode PetscDTReconstructPoly(PetscInt, PetscInt, const PetscReal *, PetscInt, const PetscReal *, PetscReal *); 122 PETSC_EXTERN PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt, PetscInt, PetscInt, PetscReal, PetscReal, PetscQuadrature *); 123 PETSC_EXTERN PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt, PetscInt, PetscInt, PetscReal, PetscReal, PetscQuadrature *); 124 PETSC_EXTERN PetscErrorCode PetscDTSimplexQuadrature(PetscInt, PetscInt, PetscDTSimplexQuadratureType, PetscQuadrature *); 125 126 PETSC_EXTERN PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt, PetscInt, PetscReal, PetscReal, PetscQuadrature *); 127 PETSC_EXTERN PetscErrorCode PetscDTTanhSinhIntegrate(void (*)(const PetscReal[], void *, PetscReal *), PetscReal, PetscReal, PetscInt, void *, PetscReal *); 128 PETSC_EXTERN PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*)(const PetscReal[], void *, PetscReal *), PetscReal, PetscReal, PetscInt, void *, PetscReal *); 129 130 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt, PetscReal *, PetscReal *, const PetscReal *, PetscReal *); 131 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 132 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 133 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***, PetscReal ***); 134 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***, PetscReal ***); 135 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 136 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 137 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 138 PETSC_EXTERN PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt, PetscReal *, PetscReal *, PetscReal ***); 139 140 PETSC_EXTERN PetscErrorCode PetscDTAltVApply(PetscInt, PetscInt, const PetscReal *, const PetscReal *, PetscReal *); 141 PETSC_EXTERN PetscErrorCode PetscDTAltVWedge(PetscInt, PetscInt, PetscInt, const PetscReal *, const PetscReal *, PetscReal *); 142 PETSC_EXTERN PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt, PetscInt, PetscInt, const PetscReal *, PetscReal *); 143 PETSC_EXTERN PetscErrorCode PetscDTAltVPullback(PetscInt, PetscInt, const PetscReal *, PetscInt, const PetscReal *, PetscReal *); 144 PETSC_EXTERN PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt, PetscInt, const PetscReal *, PetscInt, PetscReal *); 145 PETSC_EXTERN PetscErrorCode PetscDTAltVInterior(PetscInt, PetscInt, const PetscReal *, const PetscReal *, PetscReal *); 146 PETSC_EXTERN PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt, PetscInt, const PetscReal *, PetscReal *); 147 PETSC_EXTERN PetscErrorCode PetscDTAltVInteriorPattern(PetscInt, PetscInt, PetscInt (*)[3]); 148 PETSC_EXTERN PetscErrorCode PetscDTAltVStar(PetscInt, PetscInt, PetscInt, const PetscReal *, PetscReal *); 149 150 PETSC_EXTERN PetscErrorCode PetscDTBaryToIndex(PetscInt, PetscInt, const PetscInt[], PetscInt *); 151 PETSC_EXTERN PetscErrorCode PetscDTIndexToBary(PetscInt, PetscInt, PetscInt, PetscInt[]); 152 PETSC_EXTERN PetscErrorCode PetscDTGradedOrderToIndex(PetscInt, const PetscInt[], PetscInt *); 153 PETSC_EXTERN PetscErrorCode PetscDTIndexToGradedOrder(PetscInt, PetscInt, PetscInt[]); 154 155 #if defined(PETSC_USE_64BIT_INDICES) 156 #define PETSC_FACTORIAL_MAX 20 157 #define PETSC_BINOMIAL_MAX 61 158 #else 159 #define PETSC_FACTORIAL_MAX 12 160 #define PETSC_BINOMIAL_MAX 29 161 #endif 162 163 /*MC 164 PetscDTFactorial - Approximate n! as a real number 165 166 Input Parameter: 167 . n - a non-negative integer 168 169 Output Parameter: 170 . factorial - n! 171 172 Level: beginner 173 M*/ 174 static inline PetscErrorCode PetscDTFactorial(PetscInt n, PetscReal *factorial) { 175 PetscReal f = 1.0; 176 177 PetscFunctionBegin; 178 *factorial = -1.0; 179 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Factorial called with negative number %" PetscInt_FMT, n); 180 for (PetscInt i = 1; i < n + 1; ++i) f *= (PetscReal)i; 181 *factorial = f; 182 PetscFunctionReturn(0); 183 } 184 185 /*MC 186 PetscDTFactorialInt - Compute n! as an integer 187 188 Input Parameter: 189 . n - a non-negative integer 190 191 Output Parameter: 192 . factorial - n! 193 194 Level: beginner 195 196 Note: this is limited to n such that n! can be represented by PetscInt, which is 12 if PetscInt is a signed 32-bit integer and 20 if PetscInt is a signed 64-bit integer. 197 M*/ 198 static inline PetscErrorCode PetscDTFactorialInt(PetscInt n, PetscInt *factorial) { 199 PetscInt facLookup[13] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600}; 200 201 PetscFunctionBegin; 202 *factorial = -1; 203 PetscCheck(n >= 0 && n <= PETSC_FACTORIAL_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of elements %" PetscInt_FMT " is not in supported range [0,%d]", n, PETSC_FACTORIAL_MAX); 204 if (n <= 12) { 205 *factorial = facLookup[n]; 206 } else { 207 PetscInt f = facLookup[12]; 208 PetscInt i; 209 210 for (i = 13; i < n + 1; ++i) f *= i; 211 *factorial = f; 212 } 213 PetscFunctionReturn(0); 214 } 215 216 /*MC 217 PetscDTBinomial - Approximate the binomial coefficient "n choose k" 218 219 Input Parameters: 220 + n - a non-negative integer 221 - k - an integer between 0 and n, inclusive 222 223 Output Parameter: 224 . binomial - approximation of the binomial coefficient n choose k 225 226 Level: beginner 227 M*/ 228 static inline PetscErrorCode PetscDTBinomial(PetscInt n, PetscInt k, PetscReal *binomial) { 229 PetscFunctionBeginHot; 230 *binomial = -1.0; 231 PetscCheck(n >= 0 && k >= 0 && k <= n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Binomial arguments (%" PetscInt_FMT " %" PetscInt_FMT ") must be non-negative, k <= n", n, k); 232 if (n <= 3) { 233 PetscInt binomLookup[4][4] = { 234 {1, 0, 0, 0}, 235 {1, 1, 0, 0}, 236 {1, 2, 1, 0}, 237 {1, 3, 3, 1} 238 }; 239 240 *binomial = (PetscReal)binomLookup[n][k]; 241 } else { 242 PetscReal binom = 1.0; 243 244 k = PetscMin(k, n - k); 245 for (PetscInt i = 0; i < k; i++) binom = (binom * (PetscReal)(n - i)) / (PetscReal)(i + 1); 246 *binomial = binom; 247 } 248 PetscFunctionReturn(0); 249 } 250 251 /*MC 252 PetscDTBinomialInt - Compute the binomial coefficient "n choose k" 253 254 Input Parameters: 255 + n - a non-negative integer 256 - k - an integer between 0 and n, inclusive 257 258 Output Parameter: 259 . binomial - the binomial coefficient n choose k 260 261 Note: this is limited by integers that can be represented by `PetscInt` 262 263 Level: beginner 264 M*/ 265 static inline PetscErrorCode PetscDTBinomialInt(PetscInt n, PetscInt k, PetscInt *binomial) { 266 PetscInt bin; 267 268 PetscFunctionBegin; 269 *binomial = -1; 270 PetscCheck(n >= 0 && k >= 0 && k <= n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Binomial arguments (%" PetscInt_FMT " %" PetscInt_FMT ") must be non-negative, k <= n", n, k); 271 PetscCheck(n <= PETSC_BINOMIAL_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Binomial elements %" PetscInt_FMT " is larger than max for PetscInt, %d", n, PETSC_BINOMIAL_MAX); 272 if (n <= 3) { 273 PetscInt binomLookup[4][4] = { 274 {1, 0, 0, 0}, 275 {1, 1, 0, 0}, 276 {1, 2, 1, 0}, 277 {1, 3, 3, 1} 278 }; 279 280 bin = binomLookup[n][k]; 281 } else { 282 PetscInt binom = 1; 283 284 k = PetscMin(k, n - k); 285 for (PetscInt i = 0; i < k; i++) binom = (binom * (n - i)) / (i + 1); 286 bin = binom; 287 } 288 *binomial = bin; 289 PetscFunctionReturn(0); 290 } 291 292 /*MC 293 PetscDTEnumPerm - Get a permutation of n integers from its encoding into the integers [0, n!) as a sequence of swaps. 294 295 A permutation can be described by the operations that convert the lists [0, 1, ..., n-1] into the permutation, 296 by a sequence of swaps, where the ith step swaps whatever number is in ith position with a number that is in 297 some position j >= i. This swap is encoded as the difference (j - i). The difference d_i at step i is less than 298 (n - i). This sequence of n-1 differences [d_0, ..., d_{n-2}] is encoded as the number 299 (n-1)! * d_0 + (n-2)! * d_1 + ... + 1! * d_{n-2}. 300 301 Input Parameters: 302 + n - a non-negative integer (see note about limits below) 303 - k - an integer in [0, n!) 304 305 Output Parameters: 306 + perm - the permuted list of the integers [0, ..., n-1] 307 - isOdd - if not NULL, returns wether the permutation used an even or odd number of swaps. 308 309 Note: 310 Limited to n such that n! can be represented by `PetscInt`, which is 12 if `PetscInt` is a signed 32-bit integer and 20 if `PetscInt` is a signed 64-bit integer. 311 312 Level: beginner 313 M*/ 314 static inline PetscErrorCode PetscDTEnumPerm(PetscInt n, PetscInt k, PetscInt *perm, PetscBool *isOdd) { 315 PetscInt odd = 0; 316 PetscInt i; 317 PetscInt work[PETSC_FACTORIAL_MAX]; 318 PetscInt *w; 319 320 PetscFunctionBegin; 321 if (isOdd) *isOdd = PETSC_FALSE; 322 PetscCheck(n >= 0 && n <= PETSC_FACTORIAL_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of elements %" PetscInt_FMT " is not in supported range [0,%d]", n, PETSC_FACTORIAL_MAX); 323 w = &work[n - 2]; 324 for (i = 2; i <= n; i++) { 325 *(w--) = k % i; 326 k /= i; 327 } 328 for (i = 0; i < n; i++) perm[i] = i; 329 for (i = 0; i < n - 1; i++) { 330 PetscInt s = work[i]; 331 PetscInt swap = perm[i]; 332 333 perm[i] = perm[i + s]; 334 perm[i + s] = swap; 335 odd ^= (!!s); 336 } 337 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 338 PetscFunctionReturn(0); 339 } 340 341 /*MC 342 PetscDTPermIndex - Encode a permutation of n into an integer in [0, n!). This inverts `PetscDTEnumPerm`. 343 344 Input Parameters: 345 + n - a non-negative integer (see note about limits below) 346 - perm - the permuted list of the integers [0, ..., n-1] 347 348 Output Parameters: 349 + k - an integer in [0, n!) 350 - isOdd - if not NULL, returns wether the permutation used an even or odd number of swaps. 351 352 Note: 353 Limited to n such that n! can be represented by `PetscInt`, which is 12 if `PetscInt` is a signed 32-bit integer and 20 if `PetscInt` is a signed 64-bit integer. 354 355 Level: beginner 356 M*/ 357 static inline PetscErrorCode PetscDTPermIndex(PetscInt n, const PetscInt *perm, PetscInt *k, PetscBool *isOdd) { 358 PetscInt odd = 0; 359 PetscInt i, idx; 360 PetscInt work[PETSC_FACTORIAL_MAX]; 361 PetscInt iwork[PETSC_FACTORIAL_MAX]; 362 363 PetscFunctionBeginHot; 364 *k = -1; 365 if (isOdd) *isOdd = PETSC_FALSE; 366 PetscCheck(n >= 0 && n <= PETSC_FACTORIAL_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of elements %" PetscInt_FMT " is not in supported range [0,%d]", n, PETSC_FACTORIAL_MAX); 367 for (i = 0; i < n; i++) work[i] = i; /* partial permutation */ 368 for (i = 0; i < n; i++) iwork[i] = i; /* partial permutation inverse */ 369 for (idx = 0, i = 0; i < n - 1; i++) { 370 PetscInt j = perm[i]; 371 PetscInt icur = work[i]; 372 PetscInt jloc = iwork[j]; 373 PetscInt diff = jloc - i; 374 375 idx = idx * (n - i) + diff; 376 /* swap (i, jloc) */ 377 work[i] = j; 378 work[jloc] = icur; 379 iwork[j] = i; 380 iwork[icur] = jloc; 381 odd ^= (!!diff); 382 } 383 *k = idx; 384 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 385 PetscFunctionReturn(0); 386 } 387 388 /*MC 389 PetscDTEnumSubset - Get an ordered subset of the integers [0, ..., n - 1] from its encoding as an integers in [0, n choose k). 390 The encoding is in lexicographic order. 391 392 Input Parameters: 393 + n - a non-negative integer (see note about limits below) 394 . k - an integer in [0, n] 395 - j - an index in [0, n choose k) 396 397 Output Parameter: 398 . subset - the jth subset of size k of the integers [0, ..., n - 1] 399 400 Note: 401 Limited by arguments such that n choose k can be represented by `PetscInt` 402 403 Level: beginner 404 405 .seealso: `PetscDTSubsetIndex()` 406 M*/ 407 static inline PetscErrorCode PetscDTEnumSubset(PetscInt n, PetscInt k, PetscInt j, PetscInt *subset) { 408 PetscInt Nk; 409 410 PetscFunctionBeginHot; 411 PetscCall(PetscDTBinomialInt(n, k, &Nk)); 412 for (PetscInt i = 0, l = 0; i < n && l < k; i++) { 413 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 414 PetscInt Nminusk = Nk - Nminuskminus; 415 416 if (j < Nminuskminus) { 417 subset[l++] = i; 418 Nk = Nminuskminus; 419 } else { 420 j -= Nminuskminus; 421 Nk = Nminusk; 422 } 423 } 424 PetscFunctionReturn(0); 425 } 426 427 /*MC 428 PetscDTSubsetIndex - Convert an ordered subset of k integers from the set [0, ..., n - 1] to its encoding as an integers in [0, n choose k) in lexicographic order. 429 This is the inverse of `PetscDTEnumSubset`. 430 431 Input Parameters: 432 + n - a non-negative integer (see note about limits below) 433 . k - an integer in [0, n] 434 - subset - an ordered subset of the integers [0, ..., n - 1] 435 436 Output Parameter: 437 . index - the rank of the subset in lexicographic order 438 439 Note: 440 Limited by arguments such that n choose k can be represented by `PetscInt` 441 442 Level: beginner 443 444 .seealso: `PetscDTEnumSubset()` 445 M*/ 446 static inline PetscErrorCode PetscDTSubsetIndex(PetscInt n, PetscInt k, const PetscInt *subset, PetscInt *index) { 447 PetscInt j = 0, Nk; 448 449 PetscFunctionBegin; 450 *index = -1; 451 PetscCall(PetscDTBinomialInt(n, k, &Nk)); 452 for (PetscInt i = 0, l = 0; i < n && l < k; i++) { 453 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 454 PetscInt Nminusk = Nk - Nminuskminus; 455 456 if (subset[l] == i) { 457 l++; 458 Nk = Nminuskminus; 459 } else { 460 j += Nminuskminus; 461 Nk = Nminusk; 462 } 463 } 464 *index = j; 465 PetscFunctionReturn(0); 466 } 467 468 /*MC 469 PetscDTEnumSubset - Split the integers [0, ..., n - 1] into two complementary ordered subsets, the first subset of size k and being the jth subset of that size in lexicographic order. 470 471 Input Parameters: 472 + n - a non-negative integer (see note about limits below) 473 . k - an integer in [0, n] 474 - j - an index in [0, n choose k) 475 476 Output Parameters: 477 + perm - the jth subset of size k of the integers [0, ..., n - 1], followed by its complementary set. 478 - isOdd - if not NULL, return whether perm is an even or odd permutation. 479 480 Note: 481 Limited by arguments such that n choose k can be represented by `PetscInt` 482 483 Level: beginner 484 485 .seealso: `PetscDTEnumSubset()`, `PetscDTSubsetIndex()` 486 M*/ 487 static inline PetscErrorCode PetscDTEnumSplit(PetscInt n, PetscInt k, PetscInt j, PetscInt *perm, PetscBool *isOdd) { 488 PetscInt i, l, m, Nk, odd = 0; 489 PetscInt *subcomp = perm + k; 490 491 PetscFunctionBegin; 492 if (isOdd) *isOdd = PETSC_FALSE; 493 PetscCall(PetscDTBinomialInt(n, k, &Nk)); 494 for (i = 0, l = 0, m = 0; i < n && l < k; i++) { 495 PetscInt Nminuskminus = (Nk * (k - l)) / (n - i); 496 PetscInt Nminusk = Nk - Nminuskminus; 497 498 if (j < Nminuskminus) { 499 perm[l++] = i; 500 Nk = Nminuskminus; 501 } else { 502 subcomp[m++] = i; 503 j -= Nminuskminus; 504 odd ^= ((k - l) & 1); 505 Nk = Nminusk; 506 } 507 } 508 for (; i < n; i++) subcomp[m++] = i; 509 if (isOdd) *isOdd = odd ? PETSC_TRUE : PETSC_FALSE; 510 PetscFunctionReturn(0); 511 } 512 513 struct _p_PetscTabulation { 514 PetscInt K; /* Indicates a k-jet, namely tabulated derivatives up to order k */ 515 PetscInt Nr; /* The number of tabulation replicas (often 1) */ 516 PetscInt Np; /* The number of tabulation points in a replica */ 517 PetscInt Nb; /* The number of functions tabulated */ 518 PetscInt Nc; /* The number of function components */ 519 PetscInt cdim; /* The coordinate dimension */ 520 PetscReal **T; /* The tabulation T[K] of functions and their derivatives 521 T[0] = B[Nr*Np][Nb][Nc]: The basis function values at quadrature points 522 T[1] = D[Nr*Np][Nb][Nc][cdim]: The basis function derivatives at quadrature points 523 T[2] = H[Nr*Np][Nb][Nc][cdim][cdim]: The basis function second derivatives at quadrature points */ 524 }; 525 typedef struct _p_PetscTabulation *PetscTabulation; 526 527 typedef PetscErrorCode (*PetscProbFunc)(const PetscReal[], const PetscReal[], PetscReal[]); 528 529 typedef enum { 530 DTPROB_DENSITY_CONSTANT, 531 DTPROB_DENSITY_GAUSSIAN, 532 DTPROB_DENSITY_MAXWELL_BOLTZMANN, 533 DTPROB_NUM_DENSITY 534 } DTProbDensityType; 535 PETSC_EXTERN const char *const DTProbDensityTypes[]; 536 537 PETSC_EXTERN PetscErrorCode PetscPDFMaxwellBoltzmann1D(const PetscReal[], const PetscReal[], PetscReal[]); 538 PETSC_EXTERN PetscErrorCode PetscCDFMaxwellBoltzmann1D(const PetscReal[], const PetscReal[], PetscReal[]); 539 PETSC_EXTERN PetscErrorCode PetscPDFMaxwellBoltzmann2D(const PetscReal[], const PetscReal[], PetscReal[]); 540 PETSC_EXTERN PetscErrorCode PetscCDFMaxwellBoltzmann2D(const PetscReal[], const PetscReal[], PetscReal[]); 541 PETSC_EXTERN PetscErrorCode PetscPDFMaxwellBoltzmann3D(const PetscReal[], const PetscReal[], PetscReal[]); 542 PETSC_EXTERN PetscErrorCode PetscCDFMaxwellBoltzmann3D(const PetscReal[], const PetscReal[], PetscReal[]); 543 PETSC_EXTERN PetscErrorCode PetscPDFGaussian1D(const PetscReal[], const PetscReal[], PetscReal[]); 544 PETSC_EXTERN PetscErrorCode PetscCDFGaussian1D(const PetscReal[], const PetscReal[], PetscReal[]); 545 PETSC_EXTERN PetscErrorCode PetscPDFSampleGaussian1D(const PetscReal[], const PetscReal[], PetscReal[]); 546 PETSC_EXTERN PetscErrorCode PetscPDFGaussian2D(const PetscReal[], const PetscReal[], PetscReal[]); 547 PETSC_EXTERN PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal[], const PetscReal[], PetscReal[]); 548 PETSC_EXTERN PetscErrorCode PetscPDFGaussian3D(const PetscReal[], const PetscReal[], PetscReal[]); 549 PETSC_EXTERN PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal[], const PetscReal[], PetscReal[]); 550 PETSC_EXTERN PetscErrorCode PetscPDFConstant1D(const PetscReal[], const PetscReal[], PetscReal[]); 551 PETSC_EXTERN PetscErrorCode PetscCDFConstant1D(const PetscReal[], const PetscReal[], PetscReal[]); 552 PETSC_EXTERN PetscErrorCode PetscPDFSampleConstant1D(const PetscReal[], const PetscReal[], PetscReal[]); 553 PETSC_EXTERN PetscErrorCode PetscPDFConstant2D(const PetscReal[], const PetscReal[], PetscReal[]); 554 PETSC_EXTERN PetscErrorCode PetscCDFConstant2D(const PetscReal[], const PetscReal[], PetscReal[]); 555 PETSC_EXTERN PetscErrorCode PetscPDFSampleConstant2D(const PetscReal[], const PetscReal[], PetscReal[]); 556 PETSC_EXTERN PetscErrorCode PetscPDFConstant3D(const PetscReal[], const PetscReal[], PetscReal[]); 557 PETSC_EXTERN PetscErrorCode PetscCDFConstant3D(const PetscReal[], const PetscReal[], PetscReal[]); 558 PETSC_EXTERN PetscErrorCode PetscPDFSampleConstant3D(const PetscReal[], const PetscReal[], PetscReal[]); 559 PETSC_EXTERN PetscErrorCode PetscProbCreateFromOptions(PetscInt, const char[], const char[], PetscProbFunc *, PetscProbFunc *, PetscProbFunc *); 560 561 #include <petscvec.h> 562 563 PETSC_EXTERN PetscErrorCode PetscProbComputeKSStatistic(Vec, PetscProbFunc, PetscReal *); 564 565 #endif 566