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

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

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

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

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

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

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

146:    Level: intermediate

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

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

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

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

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

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

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

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

220:    Output Parameter:
221: .  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

223:    Level: intermediate

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

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

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

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

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

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

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

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

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

312:    Level: intermediate

314:    Note: 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,
315:    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,
316:    then the inverse Hodge star transformation.

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

324:   PetscFunctionBegin;
325:   PetscCheck(N >= 0 && M >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
326:   PetscCheck(PetscAbsInt(k) <= N && PetscAbsInt(k) <= M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
327:   if (N <= 3 && M <= 3) {
328:     PetscCall(PetscDTBinomialInt(M, PetscAbsInt(k), &Mk));
329:     PetscCall(PetscDTBinomialInt(N, PetscAbsInt(k), &Nk));
330:     if (!k) {
331:       Lstarw[0] = w[0];
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[3] = {1., -1., 1.};

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

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

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

369:       for (i = 0; i < Nk; i++) {
370:         PetscReal sum = 0.;

372:         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];
373:         Lstarw[i] = sum;
374:       }
375:     } else {
376:       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]);

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

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

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

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

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

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

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

445:    Level: intermediate

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

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

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

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

494:       for (i = 0; i < Nk; i++) {
495:         for (j = 0; j < Mk; j++) {
496:           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]];
497:         }
498:       }
499:     } else {
500:       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]);

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

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

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

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

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

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

553:    Output Parameter:
554: .  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}).

556:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

619:    Level: intermediate

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

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

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

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

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

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

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

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

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

686:    Level: intermediate

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

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

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

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

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

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

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

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

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

771:    Output Parameter:
772: .  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.

774:    Level: intermediate

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

782:   PetscFunctionBegin;
783:   PetscCheck(k >= 0 && k <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
784:   PetscCall(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[3] = {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:     PetscCall(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:         PetscCall(PetscDTEnumSplit(N, l, i, subset, &sOdd));
810:         PetscCall(PetscDTSubsetIndex(N, l, subset, &idx));
811:         PetscCall(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:     PetscCall(PetscFree(subset));
822:   }
823:   PetscFunctionReturn(PETSC_SUCCESS);
824: }