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