Actual source code: dtaltv.c

petsc-3.14.1 2020-11-03
Report Typos and Errors
  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 Arguments:
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 Arguments:
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: {

78:   if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
79:   if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
80:   if (N <= 3) {
81:     if (!k) {
82:       *wv = w;
83:     } else {
84:       if (N == 1)        {*wv = w * v;}
85:       else if (N == 2) {
86:         if (k == 1)      {*wv = w * v + w * v;}
87:         else             {*wv = w * (v * v - v * v);}
88:       } else {
89:         if (k == 1)      {*wv = w * v + w * v + w * v;}
90:         else if (k == 2) {
91:           *wv = w * (v * v - v * v) +
92:                 w * (v * v - v * v) +
93:                 w * (v * v - v * v);
94:         } else {
95:           *wv = w * (v * (v * v - v * v) +
96:                         v * (v * v - v * v) +
97:                         v * (v * v - v * v));
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++) {
121:           prod *= v[perm[l] * N + subset[l]];
122:         }
123:         subsum += prod;
124:       }
125:       sum += w[i] * subsum;
126:     }
127:     PetscFree2(subset, perm);
128:     *wv = sum;
129:   }
130:   return(0);
131: }

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

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

143:    Output Arguments:
144: .  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}}),
145:              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}.

147:    Level: intermediate

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

157:   if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
158:   if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree");
159:   if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension");
160:   if (N <= 3) {
161:     PetscInt Njk;

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

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

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

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

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

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

223:    Output Arguments:
224: .  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

226:    Level: intermediate

228: .seealso: PetscDTAltV, PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
229: @*/
230: PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, PetscReal *awedgeMat)
231: {
232:   PetscInt       i;

236:   if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
237:   if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree");
238:   if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension");
239:   if (N <= 3) {
240:     PetscInt Njk;

242:     PetscDTBinomialInt(N, j+k, &Njk);
243:     if (!j) {
244:       for (i = 0; i < Njk * Njk; i++) {awedgeMat[i] = 0.;}
245:       for (i = 0; i < Njk; i++) {awedgeMat[i * (Njk + 1)] = a;}
246:     } else if (!k) {
247:       for (i = 0; i < Njk; i++) {awedgeMat[i] = a[i];}
248:     } else {
249:       if (N == 2) {
250:         awedgeMat = -a; awedgeMat =  a;
251:       } else {
252:         if (j+k == 2) {
253:           awedgeMat = -a; awedgeMat =  a; awedgeMat =    0.;
254:           awedgeMat = -a; awedgeMat =    0.; awedgeMat =  a;
255:           awedgeMat =    0.; awedgeMat = -a; awedgeMat =  a;
256:         } else {
257:           awedgeMat =  a; awedgeMat = -a; awedgeMat =  a;
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++) {
282:           subsetj[m] = subset[subsetjk[m]];
283:         }
284:         for (m = 0; m < k; m++) {
285:           subsetk[m] = subset[subsetjk[j+m]];
286:         }
287:         PetscDTSubsetIndex(N, j, subsetj, &jInd);
288:         PetscDTSubsetIndex(N, k, subsetk, &kInd);
289:         awedgeMat[i * Nk + kInd] += jkOdd ? - a[jInd] : a[jInd];
290:       }
291:     }
292:     PetscFree4(subset, subsetjk, subsetj, subsetk);
293:   }
294:   return(0);
295: }

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

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

307:    Output Arguments:
308: .  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).

310:    Level: intermediate

312:    Note: negative form degrees accomodate, e.g., H-div conforming vector fields.  An H-div conforming vector field stores its degrees of freedom as (dx, dy, dz), like a 1-form,
313:    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,
314:    then the inverse Hodge star transformation.

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

324:   if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
325:   if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
326:   if (N <= 3 && M <= 3) {

328:     PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
329:     PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
330:     if (!k) {
331:       Lstarw = w;
332:     } else if (k == 1) {
333:       for (i = 0; i < Nk; i++) {
334:         PetscReal sum = 0.;

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

342:       for (i = 0; i < Nk; i++) {
343:         PetscReal sum = 0.;

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

353:       for (i = 0; i < Nk; i++) {
354:         PetscReal sum = 0.;
355:         for (j = 0; j < Mk; j++) {
356:           sum += (L[pairs[j] * N + pairs[i]] * L[pairs[j] * N + pairs[i]] -
357:                   L[pairs[j] * N + pairs[i]] * L[pairs[j] * N + pairs[i]]) * w[j];
358:         }
359:         Lstarw[i] = sum;
360:       }
361:     } else if (k == -2) {
362:       PetscInt  pairs = {{1,2},{2,0},{0,1}};
363:       PetscInt  offi = (N == 2) ? 2 : 0;
364:       PetscInt  offj = (M == 2) ? 2 : 0;

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

369:         for (j = 0; j < Mk; j++) {
370:           sum += (L[pairs[offj + j] * N + pairs[offi + i]] *
371:                   L[pairs[offj + j] * N + pairs[offi + i]] -
372:                   L[pairs[offj + j] * N + pairs[offi + i]] *
373:                   L[pairs[offj + j] * N + pairs[offi + i]]) * w[j];

375:         }
376:         Lstarw[i] = sum;
377:       }
378:     } else {
379:       PetscReal detL = L * (L * L - L * L) +
380:                        L * (L * L - L * L) +
381:                        L * (L * L - L * L);

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

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

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

428:       PetscMalloc1(Nk, &sLsw);
429:       PetscDTAltVStar(N, N - k, -1,  Lstarw, sLsw);
430:       for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i];
431:       PetscFree(sLsw);
432:     }
433:     PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
434:     PetscFree(walloc);
435:   }
436:   return(0);
437: }

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

442:    Input Arguments:
443: +  N - the dimension of the origin vector space of the linear transformation, N >= 0
444: .  M - the dimension of the image vector space of the linear transformation, M >= 0
445: .  L - a linear transformation, an [M x N] matrix in row-major format
446: -  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())

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

451:    Level: intermediate

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

465:   if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
466:   if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
467:   if (N <= 3 && M <= 3) {
468:     PetscReal mult = {1., -1., 1.};

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

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

498:       for (i = 0; i < Nk; i++) {
499:         for (j = 0; j < Mk; j++) {
500:           Lstar[i * Mk + j] = L[pairs[offj + j] * N + pairs[offi + i]] *
501:                               L[pairs[offj + j] * N + pairs[offi + i]] -
502:                               L[pairs[offj + j] * N + pairs[offi + i]] *
503:                               L[pairs[offj + j] * N + pairs[offi + i]];
504:         }
505:       }
506:     } else {
507:       PetscReal detL = L * (L * L - L * L) +
508:                        L * (L * L - L * L) +
509:                        L * (L * L - L * L);

511:       for (i = 0; i < Nk; i++) {Lstar[i] = detL;}
512:     }
513:   } else {
514:     if (k < 0) {
515:       negative = PETSC_TRUE;
516:       k = -k;
517:     }
518:     PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
519:     PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
520:     PetscDTFactorialInt(PetscAbsInt(k), &Nf);
521:     PetscMalloc5(M, &subsetw, N, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);
522:     for (i = 0; i < Nk * Mk; i++) Lstar[i] = 0.;
523:     for (i = 0; i < Mk; i++) {
524:       PetscBool iOdd;
525:       PetscInt  iidx, jidx;

527:       PetscDTEnumSplit(M, k, i, subsetw, &iOdd);
528:       iidx = negative ? Mk - 1 - i : i;
529:       iOdd = negative ? (PetscBool) (iOdd ^ ((k * (M-k)) & 1)) : PETSC_FALSE;
530:       for (j = 0; j < Nk; j++) {
531:         PetscBool jOdd;

533:         PetscDTEnumSplit(N, k, j, subsetv, &jOdd);
534:         jidx = negative ? Nk - 1 - j : j;
535:         jOdd = negative ? (PetscBool) (iOdd ^ jOdd ^ ((k * (N-k)) & 1)) : PETSC_FALSE;
536:         for (p = 0; p < Nf; p++) {
537:           PetscReal prod;
538:           PetscBool isOdd;

540:           PetscDTEnumPerm(k, p, perm, &isOdd);
541:           isOdd = (PetscBool) (isOdd ^ jOdd);
542:           prod = isOdd ? -1. : 1.;
543:           for (l = 0; l < k; l++) {
544:             prod *= L[subsetw[perm[l]] * N + subsetv[l]];
545:           }
546:           Lstar[jidx * Mk + iidx] += prod;
547:         }
548:       }
549:     }
550:     PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
551:   }
552:   return(0);
553: }

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

558:    Input Arguments:
559: +  N - the dimension of the vector space, N >= 0
560: .  k - the degree k of the k-form w, 0 <= k <= N
561: .  w - a k-form, size [N choose k]
562: -  v - an N dimensional vector

564:    Output Arguments:
565: .  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}).

567:    Level: intermediate

569: .seealso: PetscDTAltV, PetscDTAltVInteriorMatrix(), PetscDTAltVInteriorPattern(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
570: @*/
571: PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv)
572: {
573:   PetscInt        i, Nk, Nkm;
574:   PetscErrorCode  ierr;

577:   if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
578:   PetscDTBinomialInt(N, k,   &Nk);
579:   PetscDTBinomialInt(N, k-1, &Nkm);
580:   if (N <= 3) {
581:     if (k == 1) {
582:       PetscReal sum = 0.;

584:       for (i = 0; i < N; i++) {
585:         sum += w[i] * v[i];
586:       }
587:       wIntv = sum;
588:     } else if (k == N) {
589:       PetscReal mult = {1., -1., 1.};

591:       for (i = 0; i < N; i++) {
592:         wIntv[N - 1 - i] = w * v[i] * mult[i];
593:       }
594:     } else {
595:       wIntv = - w*v - w*v;
596:       wIntv =   w*v - w*v;
597:       wIntv =   w*v + w*v;
598:     }
599:   } else {
600:     PetscInt       *subset, *work;

602:     PetscMalloc2(k, &subset, k, &work);
603:     for (i = 0; i < Nkm; i++) wIntv[i] = 0.;
604:     for (i = 0; i < Nk; i++) {
605:       PetscInt  j, l, m;

607:       PetscDTEnumSubset(N, k, i, subset);
608:       for (j = 0; j < k; j++) {
609:         PetscInt  idx;
610:         PetscBool flip = (PetscBool) (j & 1);

612:         for (l = 0, m = 0; l < k; l++) {
613:           if (l != j) work[m++] = subset[l];
614:         }
615:         PetscDTSubsetIndex(N, k - 1, work, &idx);
616:         wIntv[idx] += flip ? -(w[i] * v[subset[j]]) :  (w[i] * v[subset[j]]);
617:       }
618:     }
619:     PetscFree2(subset, work);
620:   }
621:   return(0);
622: }

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

627:    Input Arguments:
628: +  N - the dimension of the vector space, N >= 0
629: .  k - the degree k of the k-forms on which intvMat acts, 0 <= k <= N
630: -  v - an N dimensional vector

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

635:    Level: intermediate

637: .seealso: PetscDTAltV, PetscDTAltVInterior(), PetscDTAltVInteriorPattern(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
638: @*/
639: PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat)
640: {
641:   PetscInt        i, Nk, Nkm;
642:   PetscErrorCode  ierr;

645:   if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
646:   PetscDTBinomialInt(N, k,   &Nk);
647:   PetscDTBinomialInt(N, k-1, &Nkm);
648:   if (N <= 3) {
649:     if (k == 1) {
650:       for (i = 0; i < N; i++) intvMat[i] = v[i];
651:     } else if (k == N) {
652:       PetscReal mult = {1., -1., 1.};

654:       for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i];
655:     } else {
656:       intvMat = -v; intvMat = -v; intvMat =    0.;
657:       intvMat =  v; intvMat =    0.; intvMat = -v;
658:       intvMat =    0.; intvMat =  v; intvMat =  v;
659:     }
660:   } else {
661:     PetscInt       *subset, *work;

663:     PetscMalloc2(k, &subset, k, &work);
664:     for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.;
665:     for (i = 0; i < Nk; i++) {
666:       PetscInt  j, l, m;

668:       PetscDTEnumSubset(N, k, i, subset);
669:       for (j = 0; j < k; j++) {
670:         PetscInt  idx;
671:         PetscBool flip = (PetscBool) (j & 1);

673:         for (l = 0, m = 0; l < k; l++) {
674:           if (l != j) work[m++] = subset[l];
675:         }
676:         PetscDTSubsetIndex(N, k - 1, work, &idx);
677:         intvMat[idx * Nk + i] += flip ? -v[subset[j]] :  v[subset[j]];
678:       }
679:     }
680:     PetscFree2(subset, work);
681:   }
682:   return(0);
683: }

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

688:    Input Arguments:
689: +  N - the dimension of the vector space, N >= 0
690: -  k - the degree of the k-forms on which intvMat from PetscDTAltVInteriorMatrix() acts, 0 <= k <= N.

692:    Output Arguments:
693: .  indices - The interior product matrix intvMat has size [(N choose (k-1)) x (N choose k)] and has (N choose k) * k
694:              non-zeros.  indices[i] and indices[i] are the row and column of a non-zero, and its value is equal to the vector
695:              coordinate v[j] if indices[i] = j, or -v[j] if indices[i] = -(j+1)

697:    Level: intermediate

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

701: .seealso: PetscDTAltV, PetscDTAltVInterior(), PetscDTAltVInteriorMatrix(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
702: @*/
703: PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices))
704: {
705:   PetscInt        i, Nk, Nkm;
706:   PetscErrorCode  ierr;

709:   if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
710:   PetscDTBinomialInt(N, k,   &Nk);
711:   PetscDTBinomialInt(N, k-1, &Nkm);
712:   if (N <= 3) {
713:     if (k == 1) {
714:       for (i = 0; i < N; i++) {
715:         indices[i] = 0;
716:         indices[i] = i;
717:         indices[i] = i;
718:       }
719:     } else if (k == N) {
720:       PetscInt val = {0, -2, 2};

722:       for (i = 0; i < N; i++) {
723:         indices[i] = N - 1 - i;
724:         indices[i] = 0;
725:         indices[i] = val[i];
726:       }
727:     } else {
728:       indices = 0; indices = 0; indices = -(1 + 1);
729:       indices = 0; indices = 1; indices = -(2 + 1);
730:       indices = 1; indices = 0; indices = 0;
731:       indices = 1; indices = 2; indices = -(2 + 1);
732:       indices = 2; indices = 1; indices = 0;
733:       indices = 2; indices = 2; indices = 1;
734:     }
735:   } else {
736:     PetscInt       *subset, *work;

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

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

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

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

764:    Input Arguments:
765: +  N - the dimension of the vector space, N >= 0
766: .  k - the degree k of the k-form w, 0 <= k <= N
767: .  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.
768: -  w - a k-form, size [N choose k]

770:    Output Arguments:
771: .  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, ... S[k-1], S', ... S'[N-k- 1]} is an odd permutation.  This implies (star)^2 w = (-1)^{k(N-k)} w, and (star)^4 w = w.

773:    Level: intermediate

775: .seealso: PetscDTAltV, PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
776: @*/
777: PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw)
778: {
779:   PetscInt        Nk, i;
780:   PetscErrorCode  ierr;

783:   if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
784:   PetscDTBinomialInt(N, k, &Nk);
785:   pow = pow % 4;
786:   pow = (pow + 4) % 4; /* make non-negative */
787:   /* pow is now 0, 1, 2, 3 */
788:   if (N <= 3) {
789:     if (pow & 1) {
790:       PetscReal mult = {1., -1., 1.};

792:       for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i];
793:     } else {
794:       for (i = 0; i < Nk; i++) starw[i] = w[i];
795:     }
796:     if (pow > 1 && ((k * (N - k)) & 1)) {
797:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
798:     }
799:   } else {
800:     PetscInt       *subset;

802:     PetscMalloc1(N, &subset);
803:     if (pow % 2) {
804:       PetscInt l = (pow == 1) ? k : N - k;
805:       for (i = 0; i < Nk; i++) {
806:         PetscBool sOdd;
807:         PetscInt  j, idx;

809:         PetscDTEnumSplit(N, l, i, subset, &sOdd);
810:         PetscDTSubsetIndex(N, l, subset, &idx);
811:         PetscDTSubsetIndex(N, N-l, &subset[l], &j);
812:         starw[j] = sOdd ? -w[idx] : w[idx];
813:       }
814:     } else {
815:       for (i = 0; i < Nk; i++) starw[i] = w[i];
816:     }
817:     /* star^2 = -1^(k * (N - k)) */
818:     if (pow > 1 && (k * (N - k)) % 2) {
819:       for (i = 0; i < Nk; i++) starw[i] = -starw[i];
820:     }
821:     PetscFree(subset);
822:   }
823:   return(0);
824: }