1 #include <petsc/private/petscimpl.h> 2 #include <petsc/private/dtimpl.h> 3 4 /*@ 5 PetscDTAltVApply - Apply a k-form to a set of k N-dimensional vectors 6 7 Input Arguments: 8 + N - the dimension of the space 9 . k - the index of the alternating form 10 . w - the alternating form 11 - v - the set of k vectors of size N, size [k x N], row-major 12 13 Output Arguments: 14 . wv - w[v[0],...,v[k-1]] 15 16 Level: intermediate 17 18 .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 19 @*/ 20 PetscErrorCode PetscDTAltVApply(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wv) 21 { 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension"); 26 if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 27 if (N <= 3) { 28 if (!k) { 29 *wv = w[0]; 30 } else { 31 if (N == 1) {*wv = w[0] * v[0];} 32 else if (N == 2) { 33 if (k == 1) {*wv = w[0] * v[0] + w[1] * v[1];} 34 else {*wv = w[0] * (v[0] * v[3] - v[1] * v[2]);} 35 } else { 36 if (k == 1) {*wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];} 37 else if (k == 2) { 38 *wv = w[0] * (v[0] * v[4] - v[1] * v[3]) + 39 w[1] * (v[0] * v[5] - v[2] * v[3]) + 40 w[2] * (v[1] * v[5] - v[2] * v[4]); 41 } else { 42 *wv = w[0] * (v[0] * (v[4] * v[8] - v[5] * v[7]) + 43 v[1] * (v[5] * v[6] - v[3] * v[8]) + 44 v[2] * (v[3] * v[7] - v[4] * v[6])); 45 } 46 } 47 } 48 } else { 49 PetscInt Nk, Nf; 50 PetscInt *subset, *perm; 51 PetscInt i, j, l; 52 PetscReal sum = 0.; 53 54 ierr = PetscDTFactorialInt(k, &Nf);CHKERRQ(ierr); 55 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 56 ierr = PetscMalloc2(k, &subset, k, &perm);CHKERRQ(ierr); 57 for (i = 0; i < Nk; i++) { 58 PetscReal subsum = 0.; 59 60 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 61 for (j = 0; j < Nf; j++) { 62 PetscBool permOdd; 63 PetscReal prod; 64 65 ierr = PetscDTEnumPerm(k, j, perm, &permOdd);CHKERRQ(ierr); 66 prod = permOdd ? -1. : 1.; 67 for (l = 0; l < k; l++) { 68 prod *= v[perm[l] * N + subset[l]]; 69 } 70 subsum += prod; 71 } 72 sum += w[i] * subsum; 73 } 74 ierr = PetscFree2(subset, perm);CHKERRQ(ierr); 75 *wv = sum; 76 } 77 PetscFunctionReturn(0); 78 } 79 80 /*@ 81 PetscDTAltVWedge - Compute the wedge product of two alternating forms 82 83 Input Arguments: 84 + N - the dimension of the space 85 . j - the index of the form a 86 . k - the index of the form b 87 . a - the j-form 88 - b - the k-form 89 90 Output Arguments: 91 . awedgeb - a wedge b, a (j+k)-form 92 93 Level: intermediate 94 95 .seealso: PetscDTAltVWedgeMatrix(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 96 @*/ 97 PetscErrorCode PetscDTAltVWedge(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, const PetscReal *b, PetscReal *awedgeb) 98 { 99 PetscInt i; 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension"); 104 if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree"); 105 if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension"); 106 if (N <= 3) { 107 PetscInt Njk; 108 109 ierr = PetscDTBinomialInt(N, j+k, &Njk);CHKERRQ(ierr); 110 if (!j) {for (i = 0; i < Njk; i++) {awedgeb[i] = a[0] * b[i];}} 111 else if (!k) {for (i = 0; i < Njk; i++) {awedgeb[i] = a[i] * b[0];}} 112 else { 113 if (N == 2) {awedgeb[0] = a[0] * b[1] - a[1] * b[0];} 114 else { 115 if (j+k == 2) { 116 awedgeb[0] = a[0] * b[1] - a[1] * b[0]; 117 awedgeb[1] = a[0] * b[2] - a[2] * b[0]; 118 awedgeb[2] = a[1] * b[2] - a[2] * b[1]; 119 } else { 120 awedgeb[0] = a[0] * b[2] - a[1] * b[1] + a[2] * b[0]; 121 } 122 } 123 } 124 } else { 125 PetscInt Njk; 126 PetscInt JKj; 127 PetscInt *subset, *subsetjk, *subsetj, *subsetk; 128 PetscInt i; 129 130 ierr = PetscDTBinomialInt(N, j+k, &Njk);CHKERRQ(ierr); 131 ierr = PetscDTBinomialInt(j+k, j, &JKj);CHKERRQ(ierr); 132 ierr = PetscMalloc4(j+k, &subset, j+k, &subsetjk, j, &subsetj, k, &subsetk);CHKERRQ(ierr); 133 for (i = 0; i < Njk; i++) { 134 PetscReal sum = 0.; 135 PetscInt l; 136 137 ierr = PetscDTEnumSubset(N, j+k, i, subset);CHKERRQ(ierr); 138 for (l = 0; l < JKj; l++) { 139 PetscBool jkOdd; 140 PetscInt m, jInd, kInd; 141 142 ierr = PetscDTEnumSplit(j+k, j, l, subsetjk, &jkOdd);CHKERRQ(ierr); 143 for (m = 0; m < j; m++) { 144 subsetj[m] = subset[subsetjk[m]]; 145 } 146 for (m = 0; m < k; m++) { 147 subsetk[m] = subset[subsetjk[j+m]]; 148 } 149 ierr = PetscDTSubsetIndex(N, j, subsetj, &jInd);CHKERRQ(ierr); 150 ierr = PetscDTSubsetIndex(N, k, subsetk, &kInd);CHKERRQ(ierr); 151 sum += jkOdd ? -(a[jInd] * b[kInd]) : (a[jInd] * b[kInd]); 152 } 153 awedgeb[i] = sum; 154 } 155 ierr = PetscFree4(subset, subsetjk, subsetj, subsetk);CHKERRQ(ierr); 156 } 157 PetscFunctionReturn(0); 158 } 159 160 /*@ 161 PetscDTAltVWedgeMatrix - Compute the matrix defined by the wedge product with an alternating form 162 163 Input Arguments: 164 + N - the dimension of the space 165 . j - the index of the form a 166 . k - the index of the form that (a wedge) will be applied to 167 - a - the j-form 168 169 Output Arguments: 170 . awedge - a wedge, an [(N choose j+k) x (N choose k)] matrix in row-major order, such that (a wedge) * b = a wedge b 171 172 Level: intermediate 173 174 .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 175 @*/ 176 PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, PetscReal *awedgeMat) 177 { 178 PetscInt i; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension"); 183 if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree"); 184 if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension"); 185 if (N <= 3) { 186 PetscInt Njk; 187 188 ierr = PetscDTBinomialInt(N, j+k, &Njk);CHKERRQ(ierr); 189 if (!j) { 190 for (i = 0; i < Njk * Njk; i++) {awedgeMat[i] = 0.;} 191 for (i = 0; i < Njk; i++) {awedgeMat[i * (Njk + 1)] = a[0];} 192 } else if (!k) { 193 for (i = 0; i < Njk; i++) {awedgeMat[i] = a[i];} 194 } else { 195 if (N == 2) { 196 awedgeMat[0] = -a[1]; awedgeMat[1] = a[0]; 197 } else { 198 if (j+k == 2) { 199 awedgeMat[0] = -a[1]; awedgeMat[1] = a[0]; awedgeMat[2] = 0.; 200 awedgeMat[3] = -a[2]; awedgeMat[4] = 0.; awedgeMat[5] = a[0]; 201 awedgeMat[6] = 0.; awedgeMat[7] = -a[2]; awedgeMat[8] = a[1]; 202 } else { 203 awedgeMat[0] = a[2]; awedgeMat[1] = -a[1]; awedgeMat[2] = a[0]; 204 } 205 } 206 } 207 } else { 208 PetscInt Njk; 209 PetscInt Nk; 210 PetscInt JKj, i; 211 PetscInt *subset, *subsetjk, *subsetj, *subsetk; 212 213 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 214 ierr = PetscDTBinomialInt(N, j+k, &Njk);CHKERRQ(ierr); 215 ierr = PetscDTBinomialInt(j+k, j, &JKj);CHKERRQ(ierr); 216 ierr = PetscMalloc4(j+k, &subset, j+k, &subsetjk, j, &subsetj, k, &subsetk);CHKERRQ(ierr); 217 for (i = 0; i < Njk * Nk; i++) awedgeMat[i] = 0.; 218 for (i = 0; i < Njk; i++) { 219 PetscInt l; 220 221 ierr = PetscDTEnumSubset(N, j+k, i, subset);CHKERRQ(ierr); 222 for (l = 0; l < JKj; l++) { 223 PetscBool jkOdd; 224 PetscInt m, jInd, kInd; 225 226 ierr = PetscDTEnumSplit(j+k, j, l, subsetjk, &jkOdd);CHKERRQ(ierr); 227 for (m = 0; m < j; m++) { 228 subsetj[m] = subset[subsetjk[m]]; 229 } 230 for (m = 0; m < k; m++) { 231 subsetk[m] = subset[subsetjk[j+m]]; 232 } 233 ierr = PetscDTSubsetIndex(N, j, subsetj, &jInd);CHKERRQ(ierr); 234 ierr = PetscDTSubsetIndex(N, k, subsetk, &kInd);CHKERRQ(ierr); 235 awedgeMat[i * Nk + kInd] += jkOdd ? - a[jInd] : a[jInd]; 236 } 237 } 238 ierr = PetscFree4(subset, subsetjk, subsetj, subsetk);CHKERRQ(ierr); 239 } 240 PetscFunctionReturn(0); 241 } 242 243 /*@ 244 PetscDTAltVPullback - Compute the pullback of an alternating form under a linear transformation 245 246 Input Arguments: 247 + N - the dimension of the origin space 248 . M - the dimension of the image space 249 . L - the linear transformation, an [M x N] matrix in row-major format 250 . k - the index of the form. A negative form degree indicates that Pullback should be conjugated by the Hodge star operator (see note) 251 - w - the k-form in the image space 252 253 Output Arguments: 254 . Lstarw - the pullback of w to a k-form in the origin space 255 256 Level: intermediate 257 258 Note: negative form degrees accomodate, e.g., H-div conforming vector fields. An H-div conforming vector field stores its degrees of freedom as (dx, dy, dz), like a 1-form, 259 but its normal trace is integrated on faces, like a 2-form. The correct pullback then is to apply the Hodge star transformation from 1-form to 2-form, pullback as a 2-form, 260 then the inverse Hodge star transformation. 261 262 .seealso: PetscDTAltVPullbackMatrix(), PetscDTAltVStar() 263 @*/ 264 PetscErrorCode PetscDTAltVPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *Lstarw) 265 { 266 PetscInt i, j, Nk, Mk; 267 PetscErrorCode ierr; 268 269 PetscFunctionBegin; 270 if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions"); 271 if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 272 if (N <= 3 && M <= 3) { 273 274 ierr = PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 275 ierr = PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 276 if (!k) { 277 Lstarw[0] = w[0]; 278 } else if (k == 1) { 279 for (i = 0; i < Nk; i++) { 280 PetscReal sum = 0.; 281 282 for (j = 0; j < Mk; j++) {sum += L[j * Nk + i] * w[j];} 283 Lstarw[i] = sum; 284 } 285 } else if (k == -1) { 286 PetscReal mult[3] = {1., -1., 1.}; 287 288 for (i = 0; i < Nk; i++) { 289 PetscReal sum = 0.; 290 291 for (j = 0; j < Mk; j++) { 292 sum += L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * w[j] * mult[j]; 293 } 294 Lstarw[i] = mult[i] * sum; 295 } 296 } else if (k == 2) { 297 PetscInt pairs[3][2] = {{0,1},{0,2},{1,2}}; 298 299 for (i = 0; i < Nk; i++) { 300 PetscReal sum = 0.; 301 for (j = 0; j < Mk; j++) { 302 sum += (L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] - 303 L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]]) * w[j]; 304 } 305 Lstarw[i] = sum; 306 } 307 } else if (k == -2) { 308 PetscInt pairs[3][2] = {{1,2},{2,0},{0,1}}; 309 PetscInt offi = (N == 2) ? 2 : 0; 310 PetscInt offj = (M == 2) ? 2 : 0; 311 312 for (i = 0; i < Nk; i++) { 313 PetscReal sum = 0.; 314 315 for (j = 0; j < Mk; j++) { 316 sum += (L[pairs[offj + j][0] * N + pairs[offi + i][0]] * 317 L[pairs[offj + j][1] * N + pairs[offi + i][1]] - 318 L[pairs[offj + j][1] * N + pairs[offi + i][0]] * 319 L[pairs[offj + j][0] * N + pairs[offi + i][1]]) * w[j]; 320 321 } 322 Lstarw[i] = sum; 323 } 324 } else { 325 PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + 326 L[1] * (L[5] * L[6] - L[3] * L[8]) + 327 L[2] * (L[3] * L[7] - L[4] * L[6]); 328 329 for (i = 0; i < Nk; i++) {Lstarw[i] = detL * w[i];} 330 } 331 } else { 332 PetscInt Nf, l, p; 333 PetscReal *Lw, *Lwv; 334 PetscInt *subsetw, *subsetv; 335 PetscInt *perm; 336 PetscReal *walloc = NULL; 337 const PetscReal *ww = NULL; 338 PetscBool negative = PETSC_FALSE; 339 340 ierr = PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 341 ierr = PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 342 ierr = PetscDTFactorialInt(PetscAbsInt(k), &Nf);CHKERRQ(ierr); 343 if (k < 0) { 344 negative = PETSC_TRUE; 345 k = -k; 346 ierr = PetscMalloc1(Mk, &walloc);CHKERRQ(ierr); 347 ierr = PetscDTAltVStar(M, M - k, 1, w, walloc);CHKERRQ(ierr); 348 ww = walloc; 349 } else { 350 ww = w; 351 } 352 ierr = PetscMalloc5(k, &subsetw, k, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);CHKERRQ(ierr); 353 for (i = 0; i < Nk; i++) Lstarw[i] = 0.; 354 for (i = 0; i < Mk; i++) { 355 ierr = PetscDTEnumSubset(M, k, i, subsetw);CHKERRQ(ierr); 356 for (j = 0; j < Nk; j++) { 357 ierr = PetscDTEnumSubset(N, k, j, subsetv);CHKERRQ(ierr); 358 for (p = 0; p < Nf; p++) { 359 PetscReal prod; 360 PetscBool isOdd; 361 362 ierr = PetscDTEnumPerm(k, p, perm, &isOdd);CHKERRQ(ierr); 363 prod = isOdd ? -ww[i] : ww[i]; 364 for (l = 0; l < k; l++) { 365 prod *= L[subsetw[perm[l]] * N + subsetv[l]]; 366 } 367 Lstarw[j] += prod; 368 } 369 } 370 } 371 if (negative) { 372 PetscReal *sLsw; 373 374 ierr = PetscMalloc1(Nk, &sLsw);CHKERRQ(ierr); 375 ierr = PetscDTAltVStar(N, N - k, -1, Lstarw, sLsw);CHKERRQ(ierr); 376 for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i]; 377 ierr = PetscFree(sLsw);CHKERRQ(ierr); 378 } 379 ierr = PetscFree5(subsetw, subsetv, perm, Lw, Lwv);CHKERRQ(ierr); 380 ierr = PetscFree(walloc);CHKERRQ(ierr); 381 } 382 PetscFunctionReturn(0); 383 } 384 385 /*@ 386 PetscDTAltVPullbackMatrix - Compute the pullback matrix for k-forms under a linear transformation 387 388 Input Arguments: 389 + N - the dimension of the origin space 390 . M - the dimension of the image space 391 . L - the linear transformation, an [M x N] matrix in row-major format 392 - k - the index of the alternating forms. A negative form degree indicates that the pullback should be conjugated by the Hodge star operator (see note in PetscDTAltvPullback()) 393 394 Output Arguments: 395 . Lstar - the pullback matrix, an [(N choose k) x (M choose k)] matrix in row-major format such that Lstar * w = L^* w 396 397 Level: intermediate 398 399 .seealso: PetscDTAltVPullback(), PetscDTAltVStar() 400 @*/ 401 PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, PetscReal *Lstar) 402 { 403 PetscInt Nk, Mk, Nf, i, j, l, p; 404 PetscReal *Lw, *Lwv; 405 PetscInt *subsetw, *subsetv; 406 PetscInt *perm; 407 PetscBool negative = PETSC_FALSE; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions"); 412 if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 413 if (N <= 3 && M <= 3) { 414 PetscReal mult[3] = {1., -1., 1.}; 415 416 ierr = PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 417 ierr = PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 418 if (!k) { 419 Lstar[0] = 1.; 420 } else if (k == 1) { 421 for (i = 0; i < Nk; i++) {for (j = 0; j < Mk; j++) {Lstar[i * Mk + j] = L[j * Nk + i];}} 422 } else if (k == -1) { 423 for (i = 0; i < Nk; i++) { 424 for (j = 0; j < Mk; j++) { 425 Lstar[i * Mk + j] = L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * mult[i] * mult[j]; 426 } 427 } 428 } else if (k == 2) { 429 PetscInt pairs[3][2] = {{0,1},{0,2},{1,2}}; 430 431 for (i = 0; i < Nk; i++) { 432 for (j = 0; j < Mk; j++) { 433 Lstar[i * Mk + j] = L[pairs[j][0] * N + pairs[i][0]] * 434 L[pairs[j][1] * N + pairs[i][1]] - 435 L[pairs[j][1] * N + pairs[i][0]] * 436 L[pairs[j][0] * N + pairs[i][1]]; 437 } 438 } 439 } else if (k == -2) { 440 PetscInt pairs[3][2] = {{1,2},{2,0},{0,1}}; 441 PetscInt offi = (N == 2) ? 2 : 0; 442 PetscInt offj = (M == 2) ? 2 : 0; 443 444 for (i = 0; i < Nk; i++) { 445 for (j = 0; j < Mk; j++) { 446 Lstar[i * Mk + j] = L[pairs[offj + j][0] * N + pairs[offi + i][0]] * 447 L[pairs[offj + j][1] * N + pairs[offi + i][1]] - 448 L[pairs[offj + j][1] * N + pairs[offi + i][0]] * 449 L[pairs[offj + j][0] * N + pairs[offi + i][1]]; 450 } 451 } 452 } else { 453 PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + 454 L[1] * (L[5] * L[6] - L[3] * L[8]) + 455 L[2] * (L[3] * L[7] - L[4] * L[6]); 456 457 for (i = 0; i < Nk; i++) {Lstar[i] = detL;} 458 } 459 } else { 460 if (k < 0) { 461 negative = PETSC_TRUE; 462 k = -k; 463 } 464 ierr = PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 465 ierr = PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 466 ierr = PetscDTFactorialInt(PetscAbsInt(k), &Nf);CHKERRQ(ierr); 467 ierr = PetscMalloc5(M, &subsetw, N, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);CHKERRQ(ierr); 468 for (i = 0; i < Nk * Mk; i++) Lstar[i] = 0.; 469 for (i = 0; i < Mk; i++) { 470 PetscBool iOdd; 471 PetscInt iidx, jidx; 472 473 ierr = PetscDTEnumSplit(M, k, i, subsetw, &iOdd);CHKERRQ(ierr); 474 iidx = negative ? Mk - 1 - i : i; 475 iOdd = negative ? iOdd ^ ((k * (M-k)) & 1) : PETSC_FALSE; 476 for (j = 0; j < Nk; j++) { 477 PetscBool jOdd; 478 479 ierr = PetscDTEnumSplit(N, k, j, subsetv, &jOdd);CHKERRQ(ierr); 480 jidx = negative ? Nk - 1 - j : j; 481 jOdd = negative ? iOdd ^ jOdd ^ ((k * (N-k)) & 1) : PETSC_FALSE; 482 for (p = 0; p < Nf; p++) { 483 PetscReal prod; 484 PetscBool isOdd; 485 486 ierr = PetscDTEnumPerm(k, p, perm, &isOdd);CHKERRQ(ierr); 487 isOdd ^= jOdd; 488 prod = isOdd ? -1. : 1.; 489 for (l = 0; l < k; l++) { 490 prod *= L[subsetw[perm[l]] * N + subsetv[l]]; 491 } 492 Lstar[jidx * Mk + iidx] += prod; 493 } 494 } 495 } 496 ierr = PetscFree5(subsetw, subsetv, perm, Lw, Lwv);CHKERRQ(ierr); 497 } 498 PetscFunctionReturn(0); 499 } 500 501 /*@ 502 PetscDTAltVInterior - Compute the interior product of an alternating form with a vector 503 504 Input Arguments: 505 + N - the dimension of the origin space 506 . k - the index of the alternating forms 507 . w - the k-form 508 - v - the N dimensional vector 509 510 Output Arguments: 511 . wIntv - the (k-1) form (w int v). 512 513 Level: intermediate 514 515 .seealso: PetscDTAltVInteriorMatrix(), PetscDTAltVInteriorPattern(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 516 @*/ 517 PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv) 518 { 519 PetscInt i, Nk, Nkm; 520 PetscErrorCode ierr; 521 522 PetscFunctionBegin; 523 if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 524 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 525 ierr = PetscDTBinomialInt(N, k-1, &Nkm);CHKERRQ(ierr); 526 if (N <= 3) { 527 if (k == 1) { 528 PetscReal sum = 0.; 529 530 for (i = 0; i < N; i++) { 531 sum += w[i] * v[i]; 532 } 533 wIntv[0] = sum; 534 } else if (k == N) { 535 PetscReal mult[3] = {1., -1., 1.}; 536 537 for (i = 0; i < N; i++) { 538 wIntv[N - 1 - i] = w[0] * v[i] * mult[i]; 539 } 540 } else { 541 wIntv[0] = - w[0]*v[1] - w[1]*v[2]; 542 wIntv[1] = w[0]*v[0] - w[2]*v[2]; 543 wIntv[2] = w[1]*v[0] + w[2]*v[1]; 544 } 545 } else { 546 PetscInt *subset, *work; 547 548 ierr = PetscMalloc2(k, &subset, k, &work);CHKERRQ(ierr); 549 for (i = 0; i < Nkm; i++) wIntv[i] = 0.; 550 for (i = 0; i < Nk; i++) { 551 PetscInt j, l, m; 552 553 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 554 for (j = 0; j < k; j++) { 555 PetscInt idx; 556 PetscBool flip = (j & 1); 557 558 for (l = 0, m = 0; l < k; l++) { 559 if (l != j) work[m++] = subset[l]; 560 } 561 ierr = PetscDTSubsetIndex(N, k - 1, work, &idx);CHKERRQ(ierr); 562 wIntv[idx] += flip ? -(w[i] * v[subset[j]]) : (w[i] * v[subset[j]]); 563 } 564 } 565 ierr = PetscFree2(subset, work);CHKERRQ(ierr); 566 } 567 PetscFunctionReturn(0); 568 } 569 570 /*@ 571 PetscDTAltVInteriorMatrix - Compute the matrix of the linear transformation induced on an alternating form by the interior product with a vector 572 573 Input Arguments: 574 + N - the dimension of the origin space 575 . k - the index of the alternating forms 576 - v - the N dimensional vector 577 578 Output Arguments: 579 . intvMat - an [(N choose (k-1)) x (N choose k)] matrix, row-major: (intvMat) * w = (w int v) 580 581 Level: intermediate 582 583 .seealso: PetscDTAltVInterior(), PetscDTAltVInteriorPattern(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 584 @*/ 585 PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat) 586 { 587 PetscInt i, Nk, Nkm; 588 PetscErrorCode ierr; 589 590 PetscFunctionBegin; 591 if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 592 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 593 ierr = PetscDTBinomialInt(N, k-1, &Nkm);CHKERRQ(ierr); 594 if (N <= 3) { 595 if (k == 1) { 596 for (i = 0; i < N; i++) intvMat[i] = v[i]; 597 } else if (k == N) { 598 PetscReal mult[3] = {1., -1., 1.}; 599 600 for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i]; 601 } else { 602 intvMat[0] = -v[1]; intvMat[1] = -v[2]; intvMat[2] = 0.; 603 intvMat[3] = v[0]; intvMat[4] = 0.; intvMat[5] = -v[2]; 604 intvMat[6] = 0.; intvMat[7] = v[0]; intvMat[8] = v[1]; 605 } 606 } else { 607 PetscInt *subset, *work; 608 609 ierr = PetscMalloc2(k, &subset, k, &work);CHKERRQ(ierr); 610 for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.; 611 for (i = 0; i < Nk; i++) { 612 PetscInt j, l, m; 613 614 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 615 for (j = 0; j < k; j++) { 616 PetscInt idx; 617 PetscBool flip = (j & 1); 618 619 for (l = 0, m = 0; l < k; l++) { 620 if (l != j) work[m++] = subset[l]; 621 } 622 ierr = PetscDTSubsetIndex(N, k - 1, work, &idx);CHKERRQ(ierr); 623 intvMat[idx * Nk + i] += flip ? -v[subset[j]] : v[subset[j]]; 624 } 625 } 626 ierr = PetscFree2(subset, work);CHKERRQ(ierr); 627 } 628 PetscFunctionReturn(0); 629 } 630 631 /*@ 632 PetscDTAltVInteriorPattern - compute the sparsity and sign pattern of the interior product matrix computed in PetscDTAltVInteriorMatrix() 633 634 Input Arguments: 635 + N - the dimension of the origin space 636 - k - the index of the alternating forms 637 638 Output Arguments: 639 . indices - The interior product matrix has (N choose k) * k non-zeros. indices[i][0] and indices[i][1] are the row and column of a non-zero, 640 and its value is equal to the vector coordinate v[j] if indices[i][2] = j, or -v[j] if indices[i][2] = -(j+1) 641 642 Level: intermediate 643 644 Note: this form is useful when the interior product needs to be computed at multiple locations, as when computing the Koszul differential 645 646 .seealso: PetscDTAltVInterior(), PetscDTAltVInteriorMatrix(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 647 @*/ 648 PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3]) 649 { 650 PetscInt i, Nk, Nkm; 651 PetscErrorCode ierr; 652 653 PetscFunctionBegin; 654 if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 655 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 656 ierr = PetscDTBinomialInt(N, k-1, &Nkm);CHKERRQ(ierr); 657 if (N <= 3) { 658 if (k == 1) { 659 for (i = 0; i < N; i++) { 660 indices[i][0] = 0; 661 indices[i][1] = i; 662 indices[i][2] = i; 663 } 664 } else if (k == N) { 665 PetscInt val[3] = {0, -2, 2}; 666 667 for (i = 0; i < N; i++) { 668 indices[i][0] = N - 1 - i; 669 indices[i][1] = 0; 670 indices[i][2] = val[i]; 671 } 672 } else { 673 indices[0][0] = 0; indices[0][1] = 0; indices[0][2] = -(1 + 1); 674 indices[1][0] = 0; indices[1][1] = 1; indices[1][2] = -(2 + 1); 675 indices[2][0] = 1; indices[2][1] = 0; indices[2][2] = 0; 676 indices[3][0] = 1; indices[3][1] = 2; indices[3][2] = -(2 + 1); 677 indices[4][0] = 2; indices[4][1] = 1; indices[4][2] = 0; 678 indices[5][0] = 2; indices[5][1] = 2; indices[5][2] = 1; 679 } 680 } else { 681 PetscInt *subset, *work; 682 683 ierr = PetscMalloc2(k, &subset, k, &work);CHKERRQ(ierr); 684 for (i = 0; i < Nk; i++) { 685 PetscInt j, l, m; 686 687 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 688 for (j = 0; j < k; j++) { 689 PetscInt idx; 690 PetscBool flip = (j & 1); 691 692 for (l = 0, m = 0; l < k; l++) { 693 if (l != j) work[m++] = subset[l]; 694 } 695 ierr = PetscDTSubsetIndex(N, k - 1, work, &idx);CHKERRQ(ierr); 696 indices[i * k + j][0] = idx; 697 indices[i * k + j][1] = i; 698 indices[i * k + j][2] = flip ? -(subset[j] + 1) : subset[j]; 699 } 700 } 701 ierr = PetscFree2(subset, work);CHKERRQ(ierr); 702 } 703 PetscFunctionReturn(0); 704 } 705 706 /*@ 707 PetscDTAltVStar - Apply a power of the Hodge star transformation to an alternating form 708 709 Input Arguments: 710 + N - the dimension of the space 711 . k - the index of the alternating form 712 . pow - the number of times to apply the Hodge star operator 713 - w - the alternating form 714 715 Output Arguments: 716 . starw = (star)^pow w 717 718 Level: intermediate 719 720 .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 721 @*/ 722 PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw) 723 { 724 PetscInt Nk, i; 725 PetscErrorCode ierr; 726 727 PetscFunctionBegin; 728 if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 729 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 730 pow = pow % 4; 731 pow = (pow + 4) % 4; /* make non-negative */ 732 /* pow is now 0, 1, 2, 3 */ 733 if (N <= 3) { 734 if (pow & 1) { 735 PetscReal mult[3] = {1., -1., 1.}; 736 737 for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i]; 738 } else { 739 for (i = 0; i < Nk; i++) starw[i] = w[i]; 740 } 741 if (pow > 1 && ((k * (N - k)) & 1)) { 742 for (i = 0; i < Nk; i++) starw[i] = -starw[i]; 743 } 744 } else { 745 PetscInt *subset; 746 747 ierr = PetscMalloc1(N, &subset);CHKERRQ(ierr); 748 if (pow % 2) { 749 PetscInt l = (pow == 1) ? k : N - k; 750 for (i = 0; i < Nk; i++) { 751 PetscBool sOdd; 752 PetscInt j, idx; 753 754 ierr = PetscDTEnumSplit(N, l, i, subset, &sOdd);CHKERRQ(ierr); 755 ierr = PetscDTSubsetIndex(N, l, subset, &idx);CHKERRQ(ierr); 756 ierr = PetscDTSubsetIndex(N, N-l, &subset[l], &j);CHKERRQ(ierr); 757 starw[j] = sOdd ? -w[idx] : w[idx]; 758 } 759 } else { 760 for (i = 0; i < Nk; i++) starw[i] = w[i]; 761 } 762 /* star^2 = -1^(k * (N - k)) */ 763 if (pow > 1 && (k * (N - k)) % 2) { 764 for (i = 0; i < Nk; i++) starw[i] = -starw[i]; 765 } 766 ierr = PetscFree(subset);CHKERRQ(ierr); 767 } 768 PetscFunctionReturn(0); 769 } 770