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)-formm 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: {
 79:   if (N <= 3) {
 80:     if (!k) {
 81:       *wv = w[0];
 82:     } else {
 83:       if (N == 1) {
 84:         *wv = w[0] * v[0];
 85:       } else if (N == 2) {
 86:         if (k == 1) {
 87:           *wv = w[0] * v[0] + w[1] * v[1];
 88:         } else {
 89:           *wv = w[0] * (v[0] * v[3] - v[1] * v[2]);
 90:         }
 91:       } else {
 92:         if (k == 1) {
 93:           *wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
 94:         } else if (k == 2) {
 95:           *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]);
 96:         } else {
 97:           *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]));
 98:         }
 99:       }
100:     }
101:   } else {
102:     PetscInt  Nk, Nf;
103:     PetscInt *subset, *perm;
104:     PetscInt  i, j, l;
105:     PetscReal sum = 0.;

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

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

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

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

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

141:    Output Parameter:
142: .  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}}),
143:              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}.

145:    Level: intermediate

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

156:   if (N <= 3) {
157:     PetscInt Njk;

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

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

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

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

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

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

218:    Output Parameter:
219: .  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

221:    Level: intermediate

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

232:   if (N <= 3) {
233:     PetscInt Njk;

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

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

277:       PetscDTEnumSubset(N, j + k, i, subset);
278:       for (l = 0; l < JKj; l++) {
279:         PetscBool jkOdd;
280:         PetscInt  m, jInd, kInd;

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

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

298:    Input Parameters:
299: +  N - the dimension of the origin vector space of the linear transformation, M >= 0
300: .  M - the dimension of the image vector space of the linear transformation, N >= 0
301: .  L - a linear transformation, an [M x N] matrix in row-major format
302: .  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
303:        the Hodge star operator (see note).
304: -  w - a |k|-form in the image space, size [M choose |k|]

306:    Output Parameter:
307: .  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).

309:    Level: intermediate

311:    Note: negative form degrees accommodate, 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,
312:    but its normal trace is integrated on faces, like a 2-form.  The correct pullback then is to apply the Hodge star transformation from (M-2)-form to 2-form, pullback as a 2-form,
313:    then the inverse Hodge star transformation.

315: .seealso: `PetscDTAltV`, `PetscDTAltVPullbackMatrix()`, `PetscDTAltVStar()`
316: @*/
317: PetscErrorCode PetscDTAltVPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *Lstarw)
318: {
319:   PetscInt i, j, Nk, Mk;

323:   if (N <= 3 && M <= 3) {
324:     PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
325:     PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
326:     if (!k) {
327:       Lstarw[0] = w[0];
328:     } else if (k == 1) {
329:       for (i = 0; i < Nk; i++) {
330:         PetscReal sum = 0.;

332:         for (j = 0; j < Mk; j++) sum += L[j * Nk + i] * w[j];
333:         Lstarw[i] = sum;
334:       }
335:     } else if (k == -1) {
336:       PetscReal mult[3] = {1., -1., 1.};

338:       for (i = 0; i < Nk; i++) {
339:         PetscReal sum = 0.;

341:         for (j = 0; j < Mk; j++) sum += L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * w[j] * mult[j];
342:         Lstarw[i] = mult[i] * sum;
343:       }
344:     } else if (k == 2) {
345:       PetscInt pairs[3][2] = {
346:         {0, 1},
347:         {0, 2},
348:         {1, 2}
349:       };

351:       for (i = 0; i < Nk; i++) {
352:         PetscReal sum = 0.;
353:         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];
354:         Lstarw[i] = sum;
355:       }
356:     } else if (k == -2) {
357:       PetscInt pairs[3][2] = {
358:         {1, 2},
359:         {2, 0},
360:         {0, 1}
361:       };
362:       PetscInt offi = (N == 2) ? 2 : 0;
363:       PetscInt offj = (M == 2) ? 2 : 0;

365:       for (i = 0; i < Nk; i++) {
366:         PetscReal sum = 0.;

368:         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];
369:         Lstarw[i] = sum;
370:       }
371:     } else {
372:       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]);

374:       for (i = 0; i < Nk; i++) Lstarw[i] = detL * w[i];
375:     }
376:   } else {
377:     PetscInt         Nf, l, p;
378:     PetscReal       *Lw, *Lwv;
379:     PetscInt        *subsetw, *subsetv;
380:     PetscInt        *perm;
381:     PetscReal       *walloc   = NULL;
382:     const PetscReal *ww       = NULL;
383:     PetscBool        negative = PETSC_FALSE;

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

407:           PetscDTEnumPerm(k, p, perm, &isOdd);
408:           prod = isOdd ? -ww[i] : ww[i];
409:           for (l = 0; l < k; l++) prod *= L[subsetw[perm[l]] * N + subsetv[l]];
410:           Lstarw[j] += prod;
411:         }
412:       }
413:     }
414:     if (negative) {
415:       PetscReal *sLsw;

417:       PetscMalloc1(Nk, &sLsw);
418:       PetscDTAltVStar(N, N - k, -1, Lstarw, sLsw);
419:       for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i];
420:       PetscFree(sLsw);
421:     }
422:     PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
423:     PetscFree(walloc);
424:   }
425:   return 0;
426: }

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

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

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

441:    Level: intermediate

443: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVStar()`
444: @*/
445: PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, PetscReal *Lstar)
446: {
447:   PetscInt   Nk, Mk, Nf, i, j, l, p;
448:   PetscReal *Lw, *Lwv;
449:   PetscInt  *subsetw, *subsetv;
450:   PetscInt  *perm;
451:   PetscBool  negative = PETSC_FALSE;

455:   if (N <= 3 && M <= 3) {
456:     PetscReal mult[3] = {1., -1., 1.};

458:     PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
459:     PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
460:     if (!k) {
461:       Lstar[0] = 1.;
462:     } else if (k == 1) {
463:       for (i = 0; i < Nk; i++) {
464:         for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[j * Nk + i];
465:       }
466:     } else if (k == -1) {
467:       for (i = 0; i < Nk; i++) {
468:         for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * mult[i] * mult[j];
469:       }
470:     } else if (k == 2) {
471:       PetscInt pairs[3][2] = {
472:         {0, 1},
473:         {0, 2},
474:         {1, 2}
475:       };

477:       for (i = 0; i < Nk; i++) {
478:         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]];
479:       }
480:     } else if (k == -2) {
481:       PetscInt pairs[3][2] = {
482:         {1, 2},
483:         {2, 0},
484:         {0, 1}
485:       };
486:       PetscInt offi = (N == 2) ? 2 : 0;
487:       PetscInt offj = (M == 2) ? 2 : 0;

489:       for (i = 0; i < Nk; i++) {
490:         for (j = 0; j < Mk; j++) {
491:           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]];
492:         }
493:       }
494:     } else {
495:       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]);

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

513:       PetscDTEnumSplit(M, k, i, subsetw, &iOdd);
514:       iidx = negative ? Mk - 1 - i : i;
515:       iOdd = negative ? (PetscBool)(iOdd ^ ((k * (M - k)) & 1)) : PETSC_FALSE;
516:       for (j = 0; j < Nk; j++) {
517:         PetscBool jOdd;

519:         PetscDTEnumSplit(N, k, j, subsetv, &jOdd);
520:         jidx = negative ? Nk - 1 - j : j;
521:         jOdd = negative ? (PetscBool)(iOdd ^ jOdd ^ ((k * (N - k)) & 1)) : PETSC_FALSE;
522:         for (p = 0; p < Nf; p++) {
523:           PetscReal prod;
524:           PetscBool isOdd;

526:           PetscDTEnumPerm(k, p, perm, &isOdd);
527:           isOdd = (PetscBool)(isOdd ^ jOdd);
528:           prod  = isOdd ? -1. : 1.;
529:           for (l = 0; l < k; l++) prod *= L[subsetw[perm[l]] * N + subsetv[l]];
530:           Lstar[jidx * Mk + iidx] += prod;
531:         }
532:       }
533:     }
534:     PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
535:   }
536:   return 0;
537: }

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

542:    Input Parameters:
543: +  N - the dimension of the vector space, N >= 0
544: .  k - the degree k of the k-form w, 0 <= k <= N
545: .  w - a k-form, size [N choose k]
546: -  v - an N dimensional vector

548:    Output Parameter:
549: .  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}).

551:    Level: intermediate

553: .seealso: `PetscDTAltV`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
554: @*/
555: PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv)
556: {
557:   PetscInt i, Nk, Nkm;

560:   PetscDTBinomialInt(N, k, &Nk);
561:   PetscDTBinomialInt(N, k - 1, &Nkm);
562:   if (N <= 3) {
563:     if (k == 1) {
564:       PetscReal sum = 0.;

566:       for (i = 0; i < N; i++) sum += w[i] * v[i];
567:       wIntv[0] = sum;
568:     } else if (k == N) {
569:       PetscReal mult[3] = {1., -1., 1.};

571:       for (i = 0; i < N; i++) wIntv[N - 1 - i] = w[0] * v[i] * mult[i];
572:     } else {
573:       wIntv[0] = -w[0] * v[1] - w[1] * v[2];
574:       wIntv[1] = w[0] * v[0] - w[2] * v[2];
575:       wIntv[2] = w[1] * v[0] + w[2] * v[1];
576:     }
577:   } else {
578:     PetscInt *subset, *work;

580:     PetscMalloc2(k, &subset, k, &work);
581:     for (i = 0; i < Nkm; i++) wIntv[i] = 0.;
582:     for (i = 0; i < Nk; i++) {
583:       PetscInt j, l, m;

585:       PetscDTEnumSubset(N, k, i, subset);
586:       for (j = 0; j < k; j++) {
587:         PetscInt  idx;
588:         PetscBool flip = (PetscBool)(j & 1);

590:         for (l = 0, m = 0; l < k; l++) {
591:           if (l != j) work[m++] = subset[l];
592:         }
593:         PetscDTSubsetIndex(N, k - 1, work, &idx);
594:         wIntv[idx] += flip ? -(w[i] * v[subset[j]]) : (w[i] * v[subset[j]]);
595:       }
596:     }
597:     PetscFree2(subset, work);
598:   }
599:   return 0;
600: }

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

605:    Input Parameters:
606: +  N - the dimension of the vector space, N >= 0
607: .  k - the degree k of the k-forms on which intvMat acts, 0 <= k <= N
608: -  v - an N dimensional vector

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

613:    Level: intermediate

615: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
616: @*/
617: PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat)
618: {
619:   PetscInt i, Nk, Nkm;

622:   PetscDTBinomialInt(N, k, &Nk);
623:   PetscDTBinomialInt(N, k - 1, &Nkm);
624:   if (N <= 3) {
625:     if (k == 1) {
626:       for (i = 0; i < N; i++) intvMat[i] = v[i];
627:     } else if (k == N) {
628:       PetscReal mult[3] = {1., -1., 1.};

630:       for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i];
631:     } else {
632:       intvMat[0] = -v[1];
633:       intvMat[1] = -v[2];
634:       intvMat[2] = 0.;
635:       intvMat[3] = v[0];
636:       intvMat[4] = 0.;
637:       intvMat[5] = -v[2];
638:       intvMat[6] = 0.;
639:       intvMat[7] = v[0];
640:       intvMat[8] = v[1];
641:     }
642:   } else {
643:     PetscInt *subset, *work;

645:     PetscMalloc2(k, &subset, k, &work);
646:     for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.;
647:     for (i = 0; i < Nk; i++) {
648:       PetscInt j, l, m;

650:       PetscDTEnumSubset(N, k, i, subset);
651:       for (j = 0; j < k; j++) {
652:         PetscInt  idx;
653:         PetscBool flip = (PetscBool)(j & 1);

655:         for (l = 0, m = 0; l < k; l++) {
656:           if (l != j) work[m++] = subset[l];
657:         }
658:         PetscDTSubsetIndex(N, k - 1, work, &idx);
659:         intvMat[idx * Nk + i] += flip ? -v[subset[j]] : v[subset[j]];
660:       }
661:     }
662:     PetscFree2(subset, work);
663:   }
664:   return 0;
665: }

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

670:    Input Parameters:
671: +  N - the dimension of the vector space, N >= 0
672: -  k - the degree of the k-forms on which intvMat from `PetscDTAltVInteriorMatrix()` acts, 0 <= k <= N.

674:    Output Parameter:
675: .  indices - The interior product matrix intvMat has size [(N choose (k-1)) x (N choose k)] and has (N choose k) * k
676:              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
677:              coordinate v[j] if indices[i][2] = j, or -v[j] if indices[i][2] = -(j+1)

679:    Level: intermediate

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

684: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
685: @*/
686: PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3])
687: {
688:   PetscInt i, Nk, Nkm;

691:   PetscDTBinomialInt(N, k, &Nk);
692:   PetscDTBinomialInt(N, k - 1, &Nkm);
693:   if (N <= 3) {
694:     if (k == 1) {
695:       for (i = 0; i < N; i++) {
696:         indices[i][0] = 0;
697:         indices[i][1] = i;
698:         indices[i][2] = i;
699:       }
700:     } else if (k == N) {
701:       PetscInt val[3] = {0, -2, 2};

703:       for (i = 0; i < N; i++) {
704:         indices[i][0] = N - 1 - i;
705:         indices[i][1] = 0;
706:         indices[i][2] = val[i];
707:       }
708:     } else {
709:       indices[0][0] = 0;
710:       indices[0][1] = 0;
711:       indices[0][2] = -(1 + 1);
712:       indices[1][0] = 0;
713:       indices[1][1] = 1;
714:       indices[1][2] = -(2 + 1);
715:       indices[2][0] = 1;
716:       indices[2][1] = 0;
717:       indices[2][2] = 0;
718:       indices[3][0] = 1;
719:       indices[3][1] = 2;
720:       indices[3][2] = -(2 + 1);
721:       indices[4][0] = 2;
722:       indices[4][1] = 1;
723:       indices[4][2] = 0;
724:       indices[5][0] = 2;
725:       indices[5][1] = 2;
726:       indices[5][2] = 1;
727:     }
728:   } else {
729:     PetscInt *subset, *work;

731:     PetscMalloc2(k, &subset, k, &work);
732:     for (i = 0; i < Nk; i++) {
733:       PetscInt j, l, m;

735:       PetscDTEnumSubset(N, k, i, subset);
736:       for (j = 0; j < k; j++) {
737:         PetscInt  idx;
738:         PetscBool flip = (PetscBool)(j & 1);

740:         for (l = 0, m = 0; l < k; l++) {
741:           if (l != j) work[m++] = subset[l];
742:         }
743:         PetscDTSubsetIndex(N, k - 1, work, &idx);
744:         indices[i * k + j][0] = idx;
745:         indices[i * k + j][1] = i;
746:         indices[i * k + j][2] = flip ? -(subset[j] + 1) : subset[j];
747:       }
748:     }
749:     PetscFree2(subset, work);
750:   }
751:   return 0;
752: }

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

757:    Input Parameters:
758: +  N - the dimension of the vector space, N >= 0
759: .  k - the degree k of the k-form w, 0 <= k <= N
760: .  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.
761: -  w - a k-form, size [N choose k]

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

766:    Level: intermediate

768: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
769: @*/
770: PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw)
771: {
772:   PetscInt Nk, i;

775:   PetscDTBinomialInt(N, k, &Nk);
776:   pow = pow % 4;
777:   pow = (pow + 4) % 4; /* make non-negative */
778:   /* pow is now 0, 1, 2, 3 */
779:   if (N <= 3) {
780:     if (pow & 1) {
781:       PetscReal mult[3] = {1., -1., 1.};

783:       for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i];
784:     } else {
785:       for (i = 0; i < Nk; i++) starw[i] = w[i];
786:     }
787:     if (pow > 1 && ((k * (N - k)) & 1)) {
788:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
789:     }
790:   } else {
791:     PetscInt *subset;

793:     PetscMalloc1(N, &subset);
794:     if (pow % 2) {
795:       PetscInt l = (pow == 1) ? k : N - k;
796:       for (i = 0; i < Nk; i++) {
797:         PetscBool sOdd;
798:         PetscInt  j, idx;

800:         PetscDTEnumSplit(N, l, i, subset, &sOdd);
801:         PetscDTSubsetIndex(N, l, subset, &idx);
802:         PetscDTSubsetIndex(N, N - l, &subset[l], &j);
803:         starw[j] = sOdd ? -w[idx] : w[idx];
804:       }
805:     } else {
806:       for (i = 0; i < Nk; i++) starw[i] = w[i];
807:     }
808:     /* star^2 = -1^(k * (N - k)) */
809:     if (pow > 1 && (k * (N - k)) % 2) {
810:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
811:     }
812:     PetscFree(subset);
813:   }
814:   return 0;
815: }