Actual source code: dtaltv.c

      1: #include <petsc/private/petscimpl.h>
      2: #include <petsc/private/dtimpl.h>

      4: /*MC
      5:    PetscDTAltV - An interface for common operations on k-forms, also known as alternating algebraic forms or alternating k-linear maps.
      6:    The name of the interface comes from the notation "Alt V" for the algebra of all k-forms acting vectors in the space V, also known as the exterior algebra of V*.

      8:    A recommended reference for this material is Section 2 "Exterior algebra and exterior calculus" in "Finite element
      9:    exterior calculus, homological techniques, and applications", by Arnold, Falk, & Winther (2006, doi:10.1017/S0962492906210018).

     11:    A k-form w (k is called the "form degree" of w) is an alternating k-linear map acting on tuples (v_1, ..., v_k) of
     12:    vectors from a vector space V and producing a real number:
     13:    - alternating: swapping any two vectors in a tuple reverses the sign of the result, e.g. w(v_1, v_2, ..., v_k) = -w(v_2, v_1, ..., v_k)
     14:    - k-linear: w acts linear in each vector separately, e.g. w(a*v + b*y, v_2, ..., v_k) = a*w(v,v_2,...,v_k) + b*w(y,v_2,...,v_k)
     15:    This action is implemented as `PetscDTAltVApply()`.

     17:    The k-forms on a vector space form a vector space themselves, Alt^k V.  The dimension of Alt^k V, if V is N dimensional, is N choose k.  (This
     18:    shows that for an N dimensional space, only 0 <= k <= N are valid form degrees.)
     19:    The standard basis for Alt^k V, used in PetscDTAltV, has one basis k-form for each ordered subset of k coordinates of the N dimensional space:
     20:    For example, if the coordinate directions of a four dimensional space are (t, x, y, z), then there are 4 choose 2 = 6 ordered subsets of two coordinates.
     21:    They are, in lexicographic order, (t, x), (t, y), (t, z), (x, y), (x, z) and (y, z).  PetscDTAltV also orders the basis of Alt^k V lexicographically
     22:    by the associated subsets.

     24:    The unit basis k-form associated with coordinates (c_1, ..., c_k) acts on a set of k vectors (v_1, ..., v_k) by creating a square matrix V where
     25:    V[i,j] = v_i[c_j] and taking the determinant of V.

     27:    If j + k <= N, then a j-form f and a k-form g can be multiplied to create a (j+k)-form using the wedge or exterior product, (f wedge g).
     28:    This is an anticommutative product, (f wedge g) = -(g wedge f).  It is sufficient to describe the wedge product of two basis forms.
     29:    Let f be the basis j-form associated with coordinates (f_1,...,f_j) and g be the basis k-form associated with coordinates (g_1,...,g_k):
     30:    - If there is any coordinate in both sets, then (f wedge g) = 0.
     31:    - Otherwise, (f wedge g) is a multiple of the basis (j+k)-form h associated with (f_1,...,f_j,g_1,...,g_k).
     32:    - In fact it is equal to either h or -h depending on how (f_1,...,f_j,g_1,...,g_k) compares to the same list of coordinates given in ascending order: if it is an even permutation of that list, then (f wedge g) = h, otherwise (f wedge g) = -h.
     33:    The wedge product is implemented for either two inputs (f and g) in `PetscDTAltVWedge()`, or for one (just f, giving a
     34:    matrix to multiply against multiple choices of g) in `PetscDTAltVWedgeMatrix()`.

     36:    If k > 0, a k-form w and a vector v can combine to make a (k-1)-form through the interior product, (w int v),
     37:    defined by (w int v)(v_1,...,v_{k-1}) = w(v,v_1,...,v_{k-1}).

     39:    The interior product is implemented for either two inputs (w and v) in PetscDTAltVInterior, for one (just v, giving a
     40:    matrix to multiply against multiple choices of w) in `PetscDTAltVInteriorMatrix()`,
     41:    or for no inputs (giving the sparsity pattern of `PetscDTAltVInteriorMatrix()`) in `PetscDTAltVInteriorPattern()`.

     43:    When there is a linear map L: V -> W from an N dimensional vector space to an M dimensional vector space,
     44:    it induces the linear pullback map L^* : Alt^k W -> Alt^k V, defined by L^* w(v_1,...,v_k) = w(L v_1, ..., L v_k).
     45:    The pullback is implemented as `PetscDTAltVPullback()` (acting on a known w) and `PetscDTAltVPullbackMatrix()` (creating a matrix that computes the actin of L^*).

     47:    Alt^k V and Alt^(N-k) V have the same dimension, and the Hodge star operator maps between them.  We note that Alt^N V is a one dimensional space, and its
     48:    basis vector is sometime called vol.  The Hodge star operator has the property that (f wedge (star g)) = (f,g) vol, where (f,g) is the simple inner product
     49:    of the basis coefficients of f and g.
     50:    Powers of the Hodge star operator can be applied with PetscDTAltVStar

     52:    Level: intermediate

     54: .seealso: `PetscDTAltVApply()`, `PetscDTAltVWedge()`, `PetscDTAltVInterior()`, `PetscDTAltVPullback()`, `PetscDTAltVStar()`
     55: M*/

     57: /*@
     58:   PetscDTAltVApply - Apply an a k-form (an alternating k-linear map) to a set of k N-dimensional vectors

     60:   Input Parameters:
     61: + N - the dimension of the vector space, N >= 0
     62: . k - the degree k of the k-form w, 0 <= k <= N
     63: . w - a k-form, size [N choose k] (each degree of freedom of a k-form is associated with a subset of k coordinates of the N-dimensional vectors.
     64:        The degrees of freedom are ordered lexicographically by their associated subsets)
     65: - v - a set of k vectors of size N, size [k x N], each vector stored contiguously

     67:   Output Parameter:
     68: . wv - w(v_1,...,v_k) = \sum_i w_i * det(V_i): the degree of freedom w_i is associated with coordinates [s_{i,1},...,s_{i,k}], and the square matrix V_i has
     69:    entry (j,k) given by the s_{i,k}'th coordinate of v_j

     71:   Level: intermediate

     73: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
     74: @*/
     75: PetscErrorCode PetscDTAltVApply(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wv)
     76: {
     77:   PetscFunctionBegin;
     78:   PetscCheck(N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
     79:   PetscCheck(k >= 0 && k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
     80:   if (N <= 3) {
     81:     if (!k) {
     82:       *wv = w[0];
     83:     } else {
     84:       if (N == 1) {
     85:         *wv = w[0] * v[0];
     86:       } else if (N == 2) {
     87:         if (k == 1) {
     88:           *wv = w[0] * v[0] + w[1] * v[1];
     89:         } else {
     90:           *wv = w[0] * (v[0] * v[3] - v[1] * v[2]);
     91:         }
     92:       } else {
     93:         if (k == 1) {
     94:           *wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
     95:         } else if (k == 2) {
     96:           *wv = w[0] * (v[0] * v[4] - v[1] * v[3]) + w[1] * (v[0] * v[5] - v[2] * v[3]) + w[2] * (v[1] * v[5] - v[2] * v[4]);
     97:         } else {
     98:           *wv = w[0] * (v[0] * (v[4] * v[8] - v[5] * v[7]) + v[1] * (v[5] * v[6] - v[3] * v[8]) + v[2] * (v[3] * v[7] - v[4] * v[6]));
     99:         }
    100:       }
    101:     }
    102:   } else {
    103:     PetscInt  Nk, Nf;
    104:     PetscInt *subset, *perm;
    105:     PetscInt  i, j, l;
    106:     PetscReal sum = 0.;

    108:     PetscCall(PetscDTFactorialInt(k, &Nf));
    109:     PetscCall(PetscDTBinomialInt(N, k, &Nk));
    110:     PetscCall(PetscMalloc2(k, &subset, k, &perm));
    111:     for (i = 0; i < Nk; i++) {
    112:       PetscReal subsum = 0.;

    114:       PetscCall(PetscDTEnumSubset(N, k, i, subset));
    115:       for (j = 0; j < Nf; j++) {
    116:         PetscBool permOdd;
    117:         PetscReal prod;

    119:         PetscCall(PetscDTEnumPerm(k, j, perm, &permOdd));
    120:         prod = permOdd ? -1. : 1.;
    121:         for (l = 0; l < k; l++) prod *= v[perm[l] * N + subset[l]];
    122:         subsum += prod;
    123:       }
    124:       sum += w[i] * subsum;
    125:     }
    126:     PetscCall(PetscFree2(subset, perm));
    127:     *wv = sum;
    128:   }
    129:   PetscFunctionReturn(PETSC_SUCCESS);
    130: }

    132: /*@
    133:   PetscDTAltVWedge - Compute the wedge product of a j-form and a k-form, giving a (j+k) form

    135:   Input Parameters:
    136: + N - the dimension of the vector space, N >= 0
    137: . j - the degree j of the j-form a, 0 <= j <= N
    138: . k - the degree k of the k-form b, 0 <= k <= N and 0 <= j+k <= N
    139: . a - a j-form, size [N choose j]
    140: - b - a k-form, size [N choose k]

    142:   Output Parameter:
    143: . awedgeb - the (j+k)-form a wedge b, size [N choose (j+k)]: (a wedge b)(v_1,...,v_{j+k}) = \sum_{s} sign(s) a(v_{s_1},...,v_{s_j}) b(v_{s_{j+1}},...,v_{s_{j+k}}),
    144:              where the sum is over permutations s such that s_1 < s_2 < ... < s_j and s_{j+1} < s_{j+2} < ... < s_{j+k}.

    146:   Level: intermediate

    148: .seealso: `PetscDTAltV`, `PetscDTAltVWedgeMatrix()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
    149: @*/
    150: PetscErrorCode PetscDTAltVWedge(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, const PetscReal *b, PetscReal *awedgeb)
    151: {
    152:   PetscInt i;

    154:   PetscFunctionBegin;
    155:   PetscCheck(N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
    156:   PetscCheck(j >= 0 && k >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree");
    157:   PetscCheck(j + k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension");
    158:   if (N <= 3) {
    159:     PetscInt Njk;

    161:     PetscCall(PetscDTBinomialInt(N, j + k, &Njk));
    162:     if (!j) {
    163:       for (i = 0; i < Njk; i++) awedgeb[i] = a[0] * b[i];
    164:     } else if (!k) {
    165:       for (i = 0; i < Njk; i++) awedgeb[i] = a[i] * b[0];
    166:     } else {
    167:       if (N == 2) {
    168:         awedgeb[0] = a[0] * b[1] - a[1] * b[0];
    169:       } else {
    170:         if (j + k == 2) {
    171:           awedgeb[0] = a[0] * b[1] - a[1] * b[0];
    172:           awedgeb[1] = a[0] * b[2] - a[2] * b[0];
    173:           awedgeb[2] = a[1] * b[2] - a[2] * b[1];
    174:         } else {
    175:           awedgeb[0] = a[0] * b[2] - a[1] * b[1] + a[2] * b[0];
    176:         }
    177:       }
    178:     }
    179:   } else {
    180:     PetscInt  Njk;
    181:     PetscInt  JKj;
    182:     PetscInt *subset, *subsetjk, *subsetj, *subsetk;
    183:     PetscInt  i;

    185:     PetscCall(PetscDTBinomialInt(N, j + k, &Njk));
    186:     PetscCall(PetscDTBinomialInt(j + k, j, &JKj));
    187:     PetscCall(PetscMalloc4(j + k, &subset, j + k, &subsetjk, j, &subsetj, k, &subsetk));
    188:     for (i = 0; i < Njk; i++) {
    189:       PetscReal sum = 0.;
    190:       PetscInt  l;

    192:       PetscCall(PetscDTEnumSubset(N, j + k, i, subset));
    193:       for (l = 0; l < JKj; l++) {
    194:         PetscBool jkOdd;
    195:         PetscInt  m, jInd, kInd;

    197:         PetscCall(PetscDTEnumSplit(j + k, j, l, subsetjk, &jkOdd));
    198:         for (m = 0; m < j; m++) subsetj[m] = subset[subsetjk[m]];
    199:         for (m = 0; m < k; m++) subsetk[m] = subset[subsetjk[j + m]];
    200:         PetscCall(PetscDTSubsetIndex(N, j, subsetj, &jInd));
    201:         PetscCall(PetscDTSubsetIndex(N, k, subsetk, &kInd));
    202:         sum += jkOdd ? -(a[jInd] * b[kInd]) : (a[jInd] * b[kInd]);
    203:       }
    204:       awedgeb[i] = sum;
    205:     }
    206:     PetscCall(PetscFree4(subset, subsetjk, subsetj, subsetk));
    207:   }
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: /*@
    212:   PetscDTAltVWedgeMatrix - Compute the matrix defined by the wedge product with a given j-form that maps k-forms to (j+k)-forms

    214:   Input Parameters:
    215: + N - the dimension of the vector space, N >= 0
    216: . j - the degree j of the j-form a, 0 <= j <= N
    217: . k - the degree k of the k-forms that (a wedge) will be applied to, 0 <= k <= N and 0 <= j+k <= N
    218: - a - a j-form, size [N choose j]

    220:   Output Parameter:
    221: . awedgeMat - (a wedge), an [(N choose j+k) x (N choose k)] matrix in row-major order, such that (a wedge) * b = a wedge b

    223:   Level: intermediate

    225: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
    226: @*/
    227: PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, PetscReal *awedgeMat)
    228: {
    229:   PetscInt i;

    231:   PetscFunctionBegin;
    232:   PetscCheck(N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
    233:   PetscCheck(j >= 0 && k >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree");
    234:   PetscCheck(j + k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension");
    235:   if (N <= 3) {
    236:     PetscInt Njk;

    238:     PetscCall(PetscDTBinomialInt(N, j + k, &Njk));
    239:     if (!j) {
    240:       for (i = 0; i < Njk * Njk; i++) awedgeMat[i] = 0.;
    241:       for (i = 0; i < Njk; i++) awedgeMat[i * (Njk + 1)] = a[0];
    242:     } else if (!k) {
    243:       for (i = 0; i < Njk; i++) awedgeMat[i] = a[i];
    244:     } else {
    245:       if (N == 2) {
    246:         awedgeMat[0] = -a[1];
    247:         awedgeMat[1] = a[0];
    248:       } else {
    249:         if (j + k == 2) {
    250:           awedgeMat[0] = -a[1];
    251:           awedgeMat[1] = a[0];
    252:           awedgeMat[2] = 0.;
    253:           awedgeMat[3] = -a[2];
    254:           awedgeMat[4] = 0.;
    255:           awedgeMat[5] = a[0];
    256:           awedgeMat[6] = 0.;
    257:           awedgeMat[7] = -a[2];
    258:           awedgeMat[8] = a[1];
    259:         } else {
    260:           awedgeMat[0] = a[2];
    261:           awedgeMat[1] = -a[1];
    262:           awedgeMat[2] = a[0];
    263:         }
    264:       }
    265:     }
    266:   } else {
    267:     PetscInt  Njk;
    268:     PetscInt  Nk;
    269:     PetscInt  JKj, i;
    270:     PetscInt *subset, *subsetjk, *subsetj, *subsetk;

    272:     PetscCall(PetscDTBinomialInt(N, k, &Nk));
    273:     PetscCall(PetscDTBinomialInt(N, j + k, &Njk));
    274:     PetscCall(PetscDTBinomialInt(j + k, j, &JKj));
    275:     PetscCall(PetscMalloc4(j + k, &subset, j + k, &subsetjk, j, &subsetj, k, &subsetk));
    276:     for (i = 0; i < Njk * Nk; i++) awedgeMat[i] = 0.;
    277:     for (i = 0; i < Njk; i++) {
    278:       PetscInt l;

    280:       PetscCall(PetscDTEnumSubset(N, j + k, i, subset));
    281:       for (l = 0; l < JKj; l++) {
    282:         PetscBool jkOdd;
    283:         PetscInt  m, jInd, kInd;

    285:         PetscCall(PetscDTEnumSplit(j + k, j, l, subsetjk, &jkOdd));
    286:         for (m = 0; m < j; m++) subsetj[m] = subset[subsetjk[m]];
    287:         for (m = 0; m < k; m++) subsetk[m] = subset[subsetjk[j + m]];
    288:         PetscCall(PetscDTSubsetIndex(N, j, subsetj, &jInd));
    289:         PetscCall(PetscDTSubsetIndex(N, k, subsetk, &kInd));
    290:         awedgeMat[i * Nk + kInd] += jkOdd ? -a[jInd] : a[jInd];
    291:       }
    292:     }
    293:     PetscCall(PetscFree4(subset, subsetjk, subsetj, subsetk));
    294:   }
    295:   PetscFunctionReturn(PETSC_SUCCESS);
    296: }

    298: /*@
    299:   PetscDTAltVPullback - Compute the pullback of a k-form under a linear transformation of the coordinate space

    301:   Input Parameters:
    302: + N - the dimension of the origin vector space of the linear transformation, M >= 0
    303: . M - the dimension of the image vector space of the linear transformation, N >= 0
    304: . L - a linear transformation, an [M x N] matrix in row-major format
    305: . k - the *signed* degree k of the |k|-form w, -(min(M,N)) <= k <= min(M,N).  A negative form degree indicates that the pullback should be conjugated by
    306:        the Hodge star operator (see note).
    307: - w - a |k|-form in the image space, size [M choose |k|]

    309:   Output Parameter:
    310: . Lstarw - the pullback of w to a |k|-form in the origin space, size [N choose |k|]: (Lstarw)(v_1,...v_k) = w(L*v_1,...,L*v_k).

    312:   Level: intermediate

    314:   Note:
    315:   Negative form degrees accommodate, e.g., H-div conforming vector fields.  An H-div conforming
    316:   vector field stores its degrees of freedom as (dx, dy, dz), like a 1-form, but its normal
    317:   trace is integrated on faces, like a 2-form.  The correct pullback then is to apply the Hodge
    318:   star transformation from (M-2)-form to 2-form, pullback as a 2-form, then invert the Hodge
    319:   star transformation.

    321: .seealso: `PetscDTAltV`, `PetscDTAltVPullbackMatrix()`, `PetscDTAltVStar()`
    322: @*/
    323: PetscErrorCode PetscDTAltVPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *Lstarw)
    324: {
    325:   PetscInt i, j, Nk, Mk;

    327:   PetscFunctionBegin;
    328:   PetscCheck(N >= 0 && M >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
    329:   PetscCheck(PetscAbsInt(k) <= N && PetscAbsInt(k) <= M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
    330:   if (N <= 3 && M <= 3) {
    331:     PetscCall(PetscDTBinomialInt(M, PetscAbsInt(k), &Mk));
    332:     PetscCall(PetscDTBinomialInt(N, PetscAbsInt(k), &Nk));
    333:     if (!k) {
    334:       Lstarw[0] = w[0];
    335:     } else if (k == 1) {
    336:       for (i = 0; i < Nk; i++) {
    337:         PetscReal sum = 0.;

    339:         for (j = 0; j < Mk; j++) sum += L[j * Nk + i] * w[j];
    340:         Lstarw[i] = sum;
    341:       }
    342:     } else if (k == -1) {
    343:       PetscReal mult[3] = {1., -1., 1.};

    345:       for (i = 0; i < Nk; i++) {
    346:         PetscReal sum = 0.;

    348:         for (j = 0; j < Mk; j++) sum += L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * w[j] * mult[j];
    349:         Lstarw[i] = mult[i] * sum;
    350:       }
    351:     } else if (k == 2) {
    352:       PetscInt pairs[3][2] = {
    353:         {0, 1},
    354:         {0, 2},
    355:         {1, 2}
    356:       };

    358:       for (i = 0; i < Nk; i++) {
    359:         PetscReal sum = 0.;
    360:         for (j = 0; j < Mk; j++) sum += (L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] - L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]]) * w[j];
    361:         Lstarw[i] = sum;
    362:       }
    363:     } else if (k == -2) {
    364:       PetscInt pairs[3][2] = {
    365:         {1, 2},
    366:         {2, 0},
    367:         {0, 1}
    368:       };
    369:       PetscInt offi = (N == 2) ? 2 : 0;
    370:       PetscInt offj = (M == 2) ? 2 : 0;

    372:       for (i = 0; i < Nk; i++) {
    373:         PetscReal sum = 0.;

    375:         for (j = 0; j < Mk; j++) sum += (L[pairs[offj + j][0] * N + pairs[offi + i][0]] * L[pairs[offj + j][1] * N + pairs[offi + i][1]] - L[pairs[offj + j][1] * N + pairs[offi + i][0]] * L[pairs[offj + j][0] * N + pairs[offi + i][1]]) * w[j];
    376:         Lstarw[i] = sum;
    377:       }
    378:     } else {
    379:       PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + L[1] * (L[5] * L[6] - L[3] * L[8]) + L[2] * (L[3] * L[7] - L[4] * L[6]);

    381:       for (i = 0; i < Nk; i++) Lstarw[i] = detL * w[i];
    382:     }
    383:   } else {
    384:     PetscInt         Nf, l, p;
    385:     PetscReal       *Lw, *Lwv;
    386:     PetscInt        *subsetw, *subsetv;
    387:     PetscInt        *perm;
    388:     PetscReal       *walloc   = NULL;
    389:     const PetscReal *ww       = NULL;
    390:     PetscBool        negative = PETSC_FALSE;

    392:     PetscCall(PetscDTBinomialInt(M, PetscAbsInt(k), &Mk));
    393:     PetscCall(PetscDTBinomialInt(N, PetscAbsInt(k), &Nk));
    394:     PetscCall(PetscDTFactorialInt(PetscAbsInt(k), &Nf));
    395:     if (k < 0) {
    396:       negative = PETSC_TRUE;
    397:       k        = -k;
    398:       PetscCall(PetscMalloc1(Mk, &walloc));
    399:       PetscCall(PetscDTAltVStar(M, M - k, 1, w, walloc));
    400:       ww = walloc;
    401:     } else {
    402:       ww = w;
    403:     }
    404:     PetscCall(PetscMalloc5(k, &subsetw, k, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv));
    405:     for (i = 0; i < Nk; i++) Lstarw[i] = 0.;
    406:     for (i = 0; i < Mk; i++) {
    407:       PetscCall(PetscDTEnumSubset(M, k, i, subsetw));
    408:       for (j = 0; j < Nk; j++) {
    409:         PetscCall(PetscDTEnumSubset(N, k, j, subsetv));
    410:         for (p = 0; p < Nf; p++) {
    411:           PetscReal prod;
    412:           PetscBool isOdd;

    414:           PetscCall(PetscDTEnumPerm(k, p, perm, &isOdd));
    415:           prod = isOdd ? -ww[i] : ww[i];
    416:           for (l = 0; l < k; l++) prod *= L[subsetw[perm[l]] * N + subsetv[l]];
    417:           Lstarw[j] += prod;
    418:         }
    419:       }
    420:     }
    421:     if (negative) {
    422:       PetscReal *sLsw;

    424:       PetscCall(PetscMalloc1(Nk, &sLsw));
    425:       PetscCall(PetscDTAltVStar(N, N - k, -1, Lstarw, sLsw));
    426:       for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i];
    427:       PetscCall(PetscFree(sLsw));
    428:     }
    429:     PetscCall(PetscFree5(subsetw, subsetv, perm, Lw, Lwv));
    430:     PetscCall(PetscFree(walloc));
    431:   }
    432:   PetscFunctionReturn(PETSC_SUCCESS);
    433: }

    435: /*@
    436:   PetscDTAltVPullbackMatrix - Compute the pullback matrix for k-forms under a linear transformation

    438:   Input Parameters:
    439: + N - the dimension of the origin vector space of the linear transformation, N >= 0
    440: . M - the dimension of the image vector space of the linear transformation, M >= 0
    441: . L - a linear transformation, an [M x N] matrix in row-major format
    442: - k - the *signed* degree k of the |k|-forms on which Lstar acts, -(min(M,N)) <= k <= min(M,N).
    443:        A negative form degree indicates that the pullback should be conjugated by the Hodge star operator (see note in `PetscDTAltvPullback()`)

    445:   Output Parameter:
    446: . Lstar - the pullback matrix, an [(N choose |k|) x (M choose |k|)] matrix in row-major format such that Lstar * w = L^* w

    448:   Level: intermediate

    450: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVStar()`
    451: @*/
    452: PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, PetscReal *Lstar)
    453: {
    454:   PetscInt   Nk, Mk, Nf, i, j, l, p;
    455:   PetscReal *Lw, *Lwv;
    456:   PetscInt  *subsetw, *subsetv;
    457:   PetscInt  *perm;
    458:   PetscBool  negative = PETSC_FALSE;

    460:   PetscFunctionBegin;
    461:   PetscCheck(N >= 0 && M >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
    462:   PetscCheck(PetscAbsInt(k) <= N && PetscAbsInt(k) <= M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
    463:   if (N <= 3 && M <= 3) {
    464:     PetscReal mult[3] = {1., -1., 1.};

    466:     PetscCall(PetscDTBinomialInt(M, PetscAbsInt(k), &Mk));
    467:     PetscCall(PetscDTBinomialInt(N, PetscAbsInt(k), &Nk));
    468:     if (!k) {
    469:       Lstar[0] = 1.;
    470:     } else if (k == 1) {
    471:       for (i = 0; i < Nk; i++) {
    472:         for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[j * Nk + i];
    473:       }
    474:     } else if (k == -1) {
    475:       for (i = 0; i < Nk; i++) {
    476:         for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * mult[i] * mult[j];
    477:       }
    478:     } else if (k == 2) {
    479:       PetscInt pairs[3][2] = {
    480:         {0, 1},
    481:         {0, 2},
    482:         {1, 2}
    483:       };

    485:       for (i = 0; i < Nk; i++) {
    486:         for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] - L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]];
    487:       }
    488:     } else if (k == -2) {
    489:       PetscInt pairs[3][2] = {
    490:         {1, 2},
    491:         {2, 0},
    492:         {0, 1}
    493:       };
    494:       PetscInt offi = (N == 2) ? 2 : 0;
    495:       PetscInt offj = (M == 2) ? 2 : 0;

    497:       for (i = 0; i < Nk; i++) {
    498:         for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[pairs[offj + j][0] * N + pairs[offi + i][0]] * L[pairs[offj + j][1] * N + pairs[offi + i][1]] - L[pairs[offj + j][1] * N + pairs[offi + i][0]] * L[pairs[offj + j][0] * N + pairs[offi + i][1]];
    499:       }
    500:     } else {
    501:       PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + L[1] * (L[5] * L[6] - L[3] * L[8]) + L[2] * (L[3] * L[7] - L[4] * L[6]);

    503:       for (i = 0; i < Nk; i++) Lstar[i] = detL;
    504:     }
    505:   } else {
    506:     if (k < 0) {
    507:       negative = PETSC_TRUE;
    508:       k        = -k;
    509:     }
    510:     PetscCall(PetscDTBinomialInt(M, PetscAbsInt(k), &Mk));
    511:     PetscCall(PetscDTBinomialInt(N, PetscAbsInt(k), &Nk));
    512:     PetscCall(PetscDTFactorialInt(PetscAbsInt(k), &Nf));
    513:     PetscCall(PetscMalloc5(M, &subsetw, N, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv));
    514:     for (i = 0; i < Nk * Mk; i++) Lstar[i] = 0.;
    515:     for (i = 0; i < Mk; i++) {
    516:       PetscBool iOdd;
    517:       PetscInt  iidx, jidx;

    519:       PetscCall(PetscDTEnumSplit(M, k, i, subsetw, &iOdd));
    520:       iidx = negative ? Mk - 1 - i : i;
    521:       iOdd = negative ? (PetscBool)(iOdd ^ ((k * (M - k)) & 1)) : PETSC_FALSE;
    522:       for (j = 0; j < Nk; j++) {
    523:         PetscBool jOdd;

    525:         PetscCall(PetscDTEnumSplit(N, k, j, subsetv, &jOdd));
    526:         jidx = negative ? Nk - 1 - j : j;
    527:         jOdd = negative ? (PetscBool)(iOdd ^ jOdd ^ ((k * (N - k)) & 1)) : PETSC_FALSE;
    528:         for (p = 0; p < Nf; p++) {
    529:           PetscReal prod;
    530:           PetscBool isOdd;

    532:           PetscCall(PetscDTEnumPerm(k, p, perm, &isOdd));
    533:           isOdd = (PetscBool)(isOdd ^ jOdd);
    534:           prod  = isOdd ? -1. : 1.;
    535:           for (l = 0; l < k; l++) prod *= L[subsetw[perm[l]] * N + subsetv[l]];
    536:           Lstar[jidx * Mk + iidx] += prod;
    537:         }
    538:       }
    539:     }
    540:     PetscCall(PetscFree5(subsetw, subsetv, perm, Lw, Lwv));
    541:   }
    542:   PetscFunctionReturn(PETSC_SUCCESS);
    543: }

    545: /*@
    546:   PetscDTAltVInterior - Compute the interior product of a k-form with a vector

    548:   Input Parameters:
    549: + N - the dimension of the vector space, N >= 0
    550: . k - the degree k of the k-form w, 0 <= k <= N
    551: . w - a k-form, size [N choose k]
    552: - v - an N dimensional vector

    554:   Output Parameter:
    555: . wIntv - the (k-1)-form (w int v), size [N choose (k-1)]: (w int v) is defined by its action on (k-1) vectors {v_1, ..., v_{k-1}} as (w inv v)(v_1, ..., v_{k-1}) = w(v, v_1, ..., v_{k-1}).

    557:   Level: intermediate

    559: .seealso: `PetscDTAltV`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
    560: @*/
    561: PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv)
    562: {
    563:   PetscInt i, Nk, Nkm;

    565:   PetscFunctionBegin;
    566:   PetscCheck(k > 0 && k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
    567:   PetscCall(PetscDTBinomialInt(N, k, &Nk));
    568:   PetscCall(PetscDTBinomialInt(N, k - 1, &Nkm));
    569:   if (N <= 3) {
    570:     if (k == 1) {
    571:       PetscReal sum = 0.;

    573:       for (i = 0; i < N; i++) sum += w[i] * v[i];
    574:       wIntv[0] = sum;
    575:     } else if (k == N) {
    576:       PetscReal mult[3] = {1., -1., 1.};

    578:       for (i = 0; i < N; i++) wIntv[N - 1 - i] = w[0] * v[i] * mult[i];
    579:     } else {
    580:       wIntv[0] = -w[0] * v[1] - w[1] * v[2];
    581:       wIntv[1] = w[0] * v[0] - w[2] * v[2];
    582:       wIntv[2] = w[1] * v[0] + w[2] * v[1];
    583:     }
    584:   } else {
    585:     PetscInt *subset, *work;

    587:     PetscCall(PetscMalloc2(k, &subset, k, &work));
    588:     for (i = 0; i < Nkm; i++) wIntv[i] = 0.;
    589:     for (i = 0; i < Nk; i++) {
    590:       PetscInt j, l, m;

    592:       PetscCall(PetscDTEnumSubset(N, k, i, subset));
    593:       for (j = 0; j < k; j++) {
    594:         PetscInt  idx;
    595:         PetscBool flip = (PetscBool)(j & 1);

    597:         for (l = 0, m = 0; l < k; l++) {
    598:           if (l != j) work[m++] = subset[l];
    599:         }
    600:         PetscCall(PetscDTSubsetIndex(N, k - 1, work, &idx));
    601:         wIntv[idx] += flip ? -(w[i] * v[subset[j]]) : (w[i] * v[subset[j]]);
    602:       }
    603:     }
    604:     PetscCall(PetscFree2(subset, work));
    605:   }
    606:   PetscFunctionReturn(PETSC_SUCCESS);
    607: }

    609: /*@
    610:   PetscDTAltVInteriorMatrix - Compute the matrix of the linear transformation induced on a k-form by the interior product with a vector

    612:   Input Parameters:
    613: + N - the dimension of the vector space, N >= 0
    614: . k - the degree k of the k-forms on which intvMat acts, 0 <= k <= N
    615: - v - an N dimensional vector

    617:   Output Parameter:
    618: . intvMat - an [(N choose (k-1)) x (N choose k)] matrix, row-major: (intvMat) * w = (w int v)

    620:   Level: intermediate

    622: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
    623: @*/
    624: PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat)
    625: {
    626:   PetscInt i, Nk, Nkm;

    628:   PetscFunctionBegin;
    629:   PetscCheck(k > 0 && k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
    630:   PetscCall(PetscDTBinomialInt(N, k, &Nk));
    631:   PetscCall(PetscDTBinomialInt(N, k - 1, &Nkm));
    632:   if (N <= 3) {
    633:     if (k == 1) {
    634:       for (i = 0; i < N; i++) intvMat[i] = v[i];
    635:     } else if (k == N) {
    636:       PetscReal mult[3] = {1., -1., 1.};

    638:       for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i];
    639:     } else {
    640:       intvMat[0] = -v[1];
    641:       intvMat[1] = -v[2];
    642:       intvMat[2] = 0.;
    643:       intvMat[3] = v[0];
    644:       intvMat[4] = 0.;
    645:       intvMat[5] = -v[2];
    646:       intvMat[6] = 0.;
    647:       intvMat[7] = v[0];
    648:       intvMat[8] = v[1];
    649:     }
    650:   } else {
    651:     PetscInt *subset, *work;

    653:     PetscCall(PetscMalloc2(k, &subset, k, &work));
    654:     for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.;
    655:     for (i = 0; i < Nk; i++) {
    656:       PetscInt j, l, m;

    658:       PetscCall(PetscDTEnumSubset(N, k, i, subset));
    659:       for (j = 0; j < k; j++) {
    660:         PetscInt  idx;
    661:         PetscBool flip = (PetscBool)(j & 1);

    663:         for (l = 0, m = 0; l < k; l++) {
    664:           if (l != j) work[m++] = subset[l];
    665:         }
    666:         PetscCall(PetscDTSubsetIndex(N, k - 1, work, &idx));
    667:         intvMat[idx * Nk + i] += flip ? -v[subset[j]] : v[subset[j]];
    668:       }
    669:     }
    670:     PetscCall(PetscFree2(subset, work));
    671:   }
    672:   PetscFunctionReturn(PETSC_SUCCESS);
    673: }

    675: /*@
    676:   PetscDTAltVInteriorPattern - compute the sparsity and sign pattern of the interior product matrix computed in `PetscDTAltVInteriorMatrix()`

    678:   Input Parameters:
    679: + N - the dimension of the vector space, $N \ge 0$
    680: - k - the degree of the k-forms on which `intvMat` from `PetscDTAltVInteriorMatrix()` acts, $ 0 le k le N $.

    682:   Output Parameter:
    683: . indices - The interior product matrix `intvMat` has dimensions [(N choose (k-1)) x (N choose k)] and has (N choose k) * k
    684:              non-zeros.  indices[i][0] and indices[i][1] are the row and column of a non-zero, and its value is equal to the vector
    685:              coordinate v[j] if indices[i][2] = j, or -v[j] if indices[i][2] = -(j+1)

    687:   Level: intermediate

    689:   Note:
    690:   This function is useful when the interior product needs to be computed at multiple locations, as when computing the Koszul differential

    692: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
    693: @*/
    694: PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3])
    695: {
    696:   PetscInt i, Nk, Nkm;

    698:   PetscFunctionBegin;
    699:   PetscCheck(k > 0 && k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
    700:   PetscCall(PetscDTBinomialInt(N, k, &Nk));
    701:   PetscCall(PetscDTBinomialInt(N, k - 1, &Nkm));
    702:   if (N <= 3) {
    703:     if (k == 1) {
    704:       for (i = 0; i < N; i++) {
    705:         indices[i][0] = 0;
    706:         indices[i][1] = i;
    707:         indices[i][2] = i;
    708:       }
    709:     } else if (k == N) {
    710:       PetscInt val[3] = {0, -2, 2};

    712:       for (i = 0; i < N; i++) {
    713:         indices[i][0] = N - 1 - i;
    714:         indices[i][1] = 0;
    715:         indices[i][2] = val[i];
    716:       }
    717:     } else {
    718:       indices[0][0] = 0;
    719:       indices[0][1] = 0;
    720:       indices[0][2] = -(1 + 1);
    721:       indices[1][0] = 0;
    722:       indices[1][1] = 1;
    723:       indices[1][2] = -(2 + 1);
    724:       indices[2][0] = 1;
    725:       indices[2][1] = 0;
    726:       indices[2][2] = 0;
    727:       indices[3][0] = 1;
    728:       indices[3][1] = 2;
    729:       indices[3][2] = -(2 + 1);
    730:       indices[4][0] = 2;
    731:       indices[4][1] = 1;
    732:       indices[4][2] = 0;
    733:       indices[5][0] = 2;
    734:       indices[5][1] = 2;
    735:       indices[5][2] = 1;
    736:     }
    737:   } else {
    738:     PetscInt *subset, *work;

    740:     PetscCall(PetscMalloc2(k, &subset, k, &work));
    741:     for (i = 0; i < Nk; i++) {
    742:       PetscInt j, l, m;

    744:       PetscCall(PetscDTEnumSubset(N, k, i, subset));
    745:       for (j = 0; j < k; j++) {
    746:         PetscInt  idx;
    747:         PetscBool flip = (PetscBool)(j & 1);

    749:         for (l = 0, m = 0; l < k; l++) {
    750:           if (l != j) work[m++] = subset[l];
    751:         }
    752:         PetscCall(PetscDTSubsetIndex(N, k - 1, work, &idx));
    753:         indices[i * k + j][0] = idx;
    754:         indices[i * k + j][1] = i;
    755:         indices[i * k + j][2] = flip ? -(subset[j] + 1) : subset[j];
    756:       }
    757:     }
    758:     PetscCall(PetscFree2(subset, work));
    759:   }
    760:   PetscFunctionReturn(PETSC_SUCCESS);
    761: }

    763: /*@
    764:   PetscDTAltVStar - Apply a power of the Hodge star operator, which maps k-forms to (N-k) forms, to a k-form

    766:   Input Parameters:
    767: + N   - the dimension of the vector space, N >= 0
    768: . k   - the degree k of the k-form w, 0 <= k <= N
    769: . pow - the number of times to apply the Hodge star operator: pow < 0 indicates that the inverse of the Hodge star operator should be applied |pow| times.
    770: - w   - a k-form, size [N choose k]

    772:   Output Parameter:
    773: . starw - (star)^pow w

    775:   Level: intermediate

    777:   Notes:
    778:   Each degree of freedom of a k-form is associated with a subset S of k coordinates of the N
    779:   dimensional vector space: the Hodge start operator (star) maps that degree of freedom to the
    780:   degree of freedom associated with S', the complement of S, with a sign change if the
    781:   permutation of coordinates {S[0], ... S[k-1], S'[0], starw- 1]} is an odd permutation.  This
    782:   implies (star)^2 w = (-1)^{k(N-k)} w, and (star)^4 w = w.

    784: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
    785: @*/
    786: PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw)
    787: {
    788:   PetscInt Nk, i;

    790:   PetscFunctionBegin;
    791:   PetscCheck(k >= 0 && k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
    792:   PetscCall(PetscDTBinomialInt(N, k, &Nk));
    793:   pow = pow % 4;
    794:   pow = (pow + 4) % 4; /* make non-negative */
    795:   /* pow is now 0, 1, 2, 3 */
    796:   if (N <= 3) {
    797:     if (pow & 1) {
    798:       PetscReal mult[3] = {1., -1., 1.};

    800:       for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i];
    801:     } else {
    802:       for (i = 0; i < Nk; i++) starw[i] = w[i];
    803:     }
    804:     if (pow > 1 && ((k * (N - k)) & 1)) {
    805:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
    806:     }
    807:   } else {
    808:     PetscInt *subset;

    810:     PetscCall(PetscMalloc1(N, &subset));
    811:     if (pow % 2) {
    812:       PetscInt l = (pow == 1) ? k : N - k;
    813:       for (i = 0; i < Nk; i++) {
    814:         PetscBool sOdd;
    815:         PetscInt  j, idx;

    817:         PetscCall(PetscDTEnumSplit(N, l, i, subset, &sOdd));
    818:         PetscCall(PetscDTSubsetIndex(N, l, subset, &idx));
    819:         PetscCall(PetscDTSubsetIndex(N, N - l, &subset[l], &j));
    820:         starw[j] = sOdd ? -w[idx] : w[idx];
    821:       }
    822:     } else {
    823:       for (i = 0; i < Nk; i++) starw[i] = w[i];
    824:     }
    825:     /* star^2 = -1^(k * (N - k)) */
    826:     if (pow > 1 && (k * (N - k)) % 2) {
    827:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
    828:     }
    829:     PetscCall(PetscFree(subset));
    830:   }
    831:   PetscFunctionReturn(PETSC_SUCCESS);
    832: }