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: the degrees of freedom are ordered lexicographically by their associated subsets)
 64: -  v - a set of k vectors of size N, size [k x N], each vector stored contiguously

 66:    Output Parameter:
 67: .  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 entry (j,k) given by the s_{i,k}'th coordinate of v_j

 69:    Level: intermediate

 71: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
 72: @*/
 73: PetscErrorCode PetscDTAltVApply(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wv)
 74: {
 77:   if (N <= 3) {
 78:     if (!k) {
 79:       *wv = w[0];
 80:     } else {
 81:       if (N == 1) {
 82:         *wv = w[0] * v[0];
 83:       } else if (N == 2) {
 84:         if (k == 1) {
 85:           *wv = w[0] * v[0] + w[1] * v[1];
 86:         } else {
 87:           *wv = w[0] * (v[0] * v[3] - v[1] * v[2]);
 88:         }
 89:       } else {
 90:         if (k == 1) {
 91:           *wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
 92:         } else if (k == 2) {
 93:           *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]);
 94:         } else {
 95:           *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]));
 96:         }
 97:       }
 98:     }
 99:   } else {
100:     PetscInt  Nk, Nf;
101:     PetscInt *subset, *perm;
102:     PetscInt  i, j, l;
103:     PetscReal sum = 0.;

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

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

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

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

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

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

143:    Level: intermediate

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

154:   if (N <= 3) {
155:     PetscInt Njk;

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

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

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

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

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

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

216:    Output Parameter:
217: .  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

219:    Level: intermediate

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

230:   if (N <= 3) {
231:     PetscInt Njk;

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

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

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

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

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

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

303:    Output Parameter:
304: .  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).

306:    Level: intermediate

308:    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,
309:    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,
310:    then the inverse Hodge star transformation.

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

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

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

335:       for (i = 0; i < Nk; i++) {
336:         PetscReal sum = 0.;

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

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

362:       for (i = 0; i < Nk; i++) {
363:         PetscReal sum = 0.;

365:         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];
366:         Lstarw[i] = sum;
367:       }
368:     } else {
369:       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]);

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

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

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

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

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

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

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

437:    Level: intermediate

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

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

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

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

485:       for (i = 0; i < Nk; i++) {
486:         for (j = 0; j < Mk; j++) {
487:           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]];
488:         }
489:       }
490:     } else {
491:       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]);

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

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

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

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

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

538:    Input Parameters:
539: +  N - the dimension of the vector space, N >= 0
540: .  k - the degree k of the k-form w, 0 <= k <= N
541: .  w - a k-form, size [N choose k]
542: -  v - an N dimensional vector

544:    Output Parameter:
545: .  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}).

547:    Level: intermediate

549: .seealso: `PetscDTAltV`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
550: @*/
551: PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv)
552: {
553:   PetscInt i, Nk, Nkm;

556:   PetscDTBinomialInt(N, k, &Nk);
557:   PetscDTBinomialInt(N, k - 1, &Nkm);
558:   if (N <= 3) {
559:     if (k == 1) {
560:       PetscReal sum = 0.;

562:       for (i = 0; i < N; i++) sum += w[i] * v[i];
563:       wIntv[0] = sum;
564:     } else if (k == N) {
565:       PetscReal mult[3] = {1., -1., 1.};

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

576:     PetscMalloc2(k, &subset, k, &work);
577:     for (i = 0; i < Nkm; i++) wIntv[i] = 0.;
578:     for (i = 0; i < Nk; i++) {
579:       PetscInt j, l, m;

581:       PetscDTEnumSubset(N, k, i, subset);
582:       for (j = 0; j < k; j++) {
583:         PetscInt  idx;
584:         PetscBool flip = (PetscBool)(j & 1);

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

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

601:    Input Parameters:
602: +  N - the dimension of the vector space, N >= 0
603: .  k - the degree k of the k-forms on which intvMat acts, 0 <= k <= N
604: -  v - an N dimensional vector

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

609:    Level: intermediate

611: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
612: @*/
613: PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat)
614: {
615:   PetscInt i, Nk, Nkm;

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

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

641:     PetscMalloc2(k, &subset, k, &work);
642:     for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.;
643:     for (i = 0; i < Nk; i++) {
644:       PetscInt j, l, m;

646:       PetscDTEnumSubset(N, k, i, subset);
647:       for (j = 0; j < k; j++) {
648:         PetscInt  idx;
649:         PetscBool flip = (PetscBool)(j & 1);

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

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

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

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

675:    Level: intermediate

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

679: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
680: @*/
681: PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3])
682: {
683:   PetscInt i, Nk, Nkm;

686:   PetscDTBinomialInt(N, k, &Nk);
687:   PetscDTBinomialInt(N, k - 1, &Nkm);
688:   if (N <= 3) {
689:     if (k == 1) {
690:       for (i = 0; i < N; i++) {
691:         indices[i][0] = 0;
692:         indices[i][1] = i;
693:         indices[i][2] = i;
694:       }
695:     } else if (k == N) {
696:       PetscInt val[3] = {0, -2, 2};

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

726:     PetscMalloc2(k, &subset, k, &work);
727:     for (i = 0; i < Nk; i++) {
728:       PetscInt j, l, m;

730:       PetscDTEnumSubset(N, k, i, subset);
731:       for (j = 0; j < k; j++) {
732:         PetscInt  idx;
733:         PetscBool flip = (PetscBool)(j & 1);

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

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

752:    Input Parameters:
753: +  N - the dimension of the vector space, N >= 0
754: .  k - the degree k of the k-form w, 0 <= k <= N
755: .  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.
756: -  w - a k-form, size [N choose k]

758:    Output Parameter:
759: .  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.

761:    Level: intermediate

763: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
764: @*/
765: PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw)
766: {
767:   PetscInt Nk, i;

770:   PetscDTBinomialInt(N, k, &Nk);
771:   pow = pow % 4;
772:   pow = (pow + 4) % 4; /* make non-negative */
773:   /* pow is now 0, 1, 2, 3 */
774:   if (N <= 3) {
775:     if (pow & 1) {
776:       PetscReal mult[3] = {1., -1., 1.};

778:       for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i];
779:     } else {
780:       for (i = 0; i < Nk; i++) starw[i] = w[i];
781:     }
782:     if (pow > 1 && ((k * (N - k)) & 1)) {
783:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
784:     }
785:   } else {
786:     PetscInt *subset;

788:     PetscMalloc1(N, &subset);
789:     if (pow % 2) {
790:       PetscInt l = (pow == 1) ? k : N - k;
791:       for (i = 0; i < Nk; i++) {
792:         PetscBool sOdd;
793:         PetscInt  j, idx;

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