Actual source code: dtaltv.c
1: #include <petsc/private/petscimpl.h>
2: #include <petsc/private/dtimpl.h>
4: /*MC
5: PetscDTAltV - An interface for common operations on k-forms, also known as alternating algebraic forms or alternating k-linear maps.
6: The name of the interface comes from the notation "Alt V" for the algebra of all k-forms acting vectors in the space V, also known as the exterior algebra of V*.
8: A recommended reference for this material is Section 2 "Exterior algebra and exterior calculus" in "Finite element
9: exterior calculus, homological techniques, and applications", by Arnold, Falk, & Winther (2006, doi:10.1017/S0962492906210018).
11: A k-form w (k is called the "form degree" of w) is an alternating k-linear map acting on tuples (v_1, ..., v_k) of
12: vectors from a vector space V and producing a real number:
13: - alternating: swapping any two vectors in a tuple reverses the sign of the result, e.g. w(v_1, v_2, ..., v_k) = -w(v_2, v_1, ..., v_k)
14: - k-linear: w acts linear in each vector separately, e.g. w(a*v + b*y, v_2, ..., v_k) = a*w(v,v_2,...,v_k) + b*w(y,v_2,...,v_k)
15: This action is implemented as `PetscDTAltVApply()`.
17: The k-forms on a vector space form a vector space themselves, Alt^k V. The dimension of Alt^k V, if V is N dimensional, is N choose k. (This
18: shows that for an N dimensional space, only 0 <= k <= N are valid form degrees.)
19: The standard basis for Alt^k V, used in PetscDTAltV, has one basis k-form for each ordered subset of k coordinates of the N dimensional space:
20: For example, if the coordinate directions of a four dimensional space are (t, x, y, z), then there are 4 choose 2 = 6 ordered subsets of two coordinates.
21: They are, in lexicographic order, (t, x), (t, y), (t, z), (x, y), (x, z) and (y, z). PetscDTAltV also orders the basis of Alt^k V lexicographically
22: by the associated subsets.
24: The unit basis k-form associated with coordinates (c_1, ..., c_k) acts on a set of k vectors (v_1, ..., v_k) by creating a square matrix V where
25: V[i,j] = v_i[c_j] and taking the determinant of V.
27: If j + k <= N, then a j-form f and a k-form g can be multiplied to create a (j+k)-form using the wedge or exterior product, (f wedge g).
28: This is an anticommutative product, (f wedge g) = -(g wedge f). It is sufficient to describe the wedge product of two basis forms.
29: Let f be the basis j-form associated with coordinates (f_1,...,f_j) and g be the basis k-form associated with coordinates (g_1,...,g_k):
30: - If there is any coordinate in both sets, then (f wedge g) = 0.
31: - Otherwise, (f wedge g) is a multiple of the basis (j+k)-form h associated with (f_1,...,f_j,g_1,...,g_k).
32: - In fact it is equal to either h or -h depending on how (f_1,...,f_j,g_1,...,g_k) compares to the same list of coordinates given in ascending order: if it is an even permutation of that list, then (f wedge g) = h, otherwise (f wedge g) = -h.
33: The wedge product is implemented for either two inputs (f and g) in `PetscDTAltVWedge()`, or for one (just f, giving a
34: matrix to multiply against multiple choices of g) in `PetscDTAltVWedgeMatrix()`.
36: If k > 0, a k-form w and a vector v can combine to make a (k-1)-formm through the interior product, (w int v),
37: defined by (w int v)(v_1,...,v_{k-1}) = w(v,v_1,...,v_{k-1}).
39: The interior product is implemented for either two inputs (w and v) in PetscDTAltVInterior, for one (just v, giving a
40: matrix to multiply against multiple choices of w) in `PetscDTAltVInteriorMatrix()`,
41: or for no inputs (giving the sparsity pattern of `PetscDTAltVInteriorMatrix()`) in `PetscDTAltVInteriorPattern()`.
43: When there is a linear map L: V -> W from an N dimensional vector space to an M dimensional vector space,
44: it induces the linear pullback map L^* : Alt^k W -> Alt^k V, defined by L^* w(v_1,...,v_k) = w(L v_1, ..., L v_k).
45: The pullback is implemented as `PetscDTAltVPullback()` (acting on a known w) and `PetscDTAltVPullbackMatrix()` (creating a matrix that computes the actin of L^*).
47: Alt^k V and Alt^(N-k) V have the same dimension, and the Hodge star operator maps between them. We note that Alt^N V is a one dimensional space, and its
48: basis vector is sometime called vol. The Hodge star operator has the property that (f wedge (star g)) = (f,g) vol, where (f,g) is the simple inner product
49: of the basis coefficients of f and g.
50: Powers of the Hodge star operator can be applied with PetscDTAltVStar
52: Level: intermediate
54: .seealso: `PetscDTAltVApply()`, `PetscDTAltVWedge()`, `PetscDTAltVInterior()`, `PetscDTAltVPullback()`, `PetscDTAltVStar()`
55: M*/
57: /*@
58: PetscDTAltVApply - Apply an a k-form (an alternating k-linear map) to a set of k N-dimensional vectors
60: Input Parameters:
61: + N - the dimension of the vector space, N >= 0
62: . k - the degree k of the k-form w, 0 <= k <= N
63: . w - a k-form, size [N choose k] (each degree of freedom of a k-form is associated with a subset of k coordinates of the N-dimensional vectors.
64: The degrees of freedom are ordered lexicographically by their associated subsets)
65: - v - a set of k vectors of size N, size [k x N], each vector stored contiguously
67: Output Parameter:
68: . wv - w(v_1,...,v_k) = \sum_i w_i * det(V_i): the degree of freedom w_i is associated with coordinates [s_{i,1},...,s_{i,k}], and the square matrix V_i has
69: entry (j,k) given by the s_{i,k}'th coordinate of v_j
71: Level: intermediate
73: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
74: @*/
75: PetscErrorCode PetscDTAltVApply(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wv)
76: {
79: if (N <= 3) {
80: if (!k) {
81: *wv = w[0];
82: } else {
83: if (N == 1) {
84: *wv = w[0] * v[0];
85: } else if (N == 2) {
86: if (k == 1) {
87: *wv = w[0] * v[0] + w[1] * v[1];
88: } else {
89: *wv = w[0] * (v[0] * v[3] - v[1] * v[2]);
90: }
91: } else {
92: if (k == 1) {
93: *wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
94: } else if (k == 2) {
95: *wv = w[0] * (v[0] * v[4] - v[1] * v[3]) + w[1] * (v[0] * v[5] - v[2] * v[3]) + w[2] * (v[1] * v[5] - v[2] * v[4]);
96: } else {
97: *wv = w[0] * (v[0] * (v[4] * v[8] - v[5] * v[7]) + v[1] * (v[5] * v[6] - v[3] * v[8]) + v[2] * (v[3] * v[7] - v[4] * v[6]));
98: }
99: }
100: }
101: } else {
102: PetscInt Nk, Nf;
103: PetscInt *subset, *perm;
104: PetscInt i, j, l;
105: PetscReal sum = 0.;
107: PetscDTFactorialInt(k, &Nf);
108: PetscDTBinomialInt(N, k, &Nk);
109: PetscMalloc2(k, &subset, k, &perm);
110: for (i = 0; i < Nk; i++) {
111: PetscReal subsum = 0.;
113: PetscDTEnumSubset(N, k, i, subset);
114: for (j = 0; j < Nf; j++) {
115: PetscBool permOdd;
116: PetscReal prod;
118: PetscDTEnumPerm(k, j, perm, &permOdd);
119: prod = permOdd ? -1. : 1.;
120: for (l = 0; l < k; l++) prod *= v[perm[l] * N + subset[l]];
121: subsum += prod;
122: }
123: sum += w[i] * subsum;
124: }
125: PetscFree2(subset, perm);
126: *wv = sum;
127: }
128: return 0;
129: }
131: /*@
132: PetscDTAltVWedge - Compute the wedge product of a j-form and a k-form, giving a (j+k) form
134: Input Parameters:
135: + N - the dimension of the vector space, N >= 0
136: . j - the degree j of the j-form a, 0 <= j <= N
137: . k - the degree k of the k-form b, 0 <= k <= N and 0 <= j+k <= N
138: . a - a j-form, size [N choose j]
139: - b - a k-form, size [N choose k]
141: Output Parameter:
142: . awedgeb - the (j+k)-form a wedge b, size [N choose (j+k)]: (a wedge b)(v_1,...,v_{j+k}) = \sum_{s} sign(s) a(v_{s_1},...,v_{s_j}) b(v_{s_{j+1}},...,v_{s_{j+k}}),
143: where the sum is over permutations s such that s_1 < s_2 < ... < s_j and s_{j+1} < s_{j+2} < ... < s_{j+k}.
145: Level: intermediate
147: .seealso: `PetscDTAltV`, `PetscDTAltVWedgeMatrix()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
148: @*/
149: PetscErrorCode PetscDTAltVWedge(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, const PetscReal *b, PetscReal *awedgeb)
150: {
151: PetscInt i;
156: if (N <= 3) {
157: PetscInt Njk;
159: PetscDTBinomialInt(N, j + k, &Njk);
160: if (!j) {
161: for (i = 0; i < Njk; i++) awedgeb[i] = a[0] * b[i];
162: } else if (!k) {
163: for (i = 0; i < Njk; i++) awedgeb[i] = a[i] * b[0];
164: } else {
165: if (N == 2) {
166: awedgeb[0] = a[0] * b[1] - a[1] * b[0];
167: } else {
168: if (j + k == 2) {
169: awedgeb[0] = a[0] * b[1] - a[1] * b[0];
170: awedgeb[1] = a[0] * b[2] - a[2] * b[0];
171: awedgeb[2] = a[1] * b[2] - a[2] * b[1];
172: } else {
173: awedgeb[0] = a[0] * b[2] - a[1] * b[1] + a[2] * b[0];
174: }
175: }
176: }
177: } else {
178: PetscInt Njk;
179: PetscInt JKj;
180: PetscInt *subset, *subsetjk, *subsetj, *subsetk;
181: PetscInt i;
183: PetscDTBinomialInt(N, j + k, &Njk);
184: PetscDTBinomialInt(j + k, j, &JKj);
185: PetscMalloc4(j + k, &subset, j + k, &subsetjk, j, &subsetj, k, &subsetk);
186: for (i = 0; i < Njk; i++) {
187: PetscReal sum = 0.;
188: PetscInt l;
190: PetscDTEnumSubset(N, j + k, i, subset);
191: for (l = 0; l < JKj; l++) {
192: PetscBool jkOdd;
193: PetscInt m, jInd, kInd;
195: PetscDTEnumSplit(j + k, j, l, subsetjk, &jkOdd);
196: for (m = 0; m < j; m++) subsetj[m] = subset[subsetjk[m]];
197: for (m = 0; m < k; m++) subsetk[m] = subset[subsetjk[j + m]];
198: PetscDTSubsetIndex(N, j, subsetj, &jInd);
199: PetscDTSubsetIndex(N, k, subsetk, &kInd);
200: sum += jkOdd ? -(a[jInd] * b[kInd]) : (a[jInd] * b[kInd]);
201: }
202: awedgeb[i] = sum;
203: }
204: PetscFree4(subset, subsetjk, subsetj, subsetk);
205: }
206: return 0;
207: }
209: /*@
210: PetscDTAltVWedgeMatrix - Compute the matrix defined by the wedge product with a given j-form that maps k-forms to (j+k)-forms
212: Input Parameters:
213: + N - the dimension of the vector space, N >= 0
214: . j - the degree j of the j-form a, 0 <= j <= N
215: . k - the degree k of the k-forms that (a wedge) will be applied to, 0 <= k <= N and 0 <= j+k <= N
216: - a - a j-form, size [N choose j]
218: Output Parameter:
219: . awedge - (a wedge), an [(N choose j+k) x (N choose k)] matrix in row-major order, such that (a wedge) * b = a wedge b
221: Level: intermediate
223: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
224: @*/
225: PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, PetscReal *awedgeMat)
226: {
227: PetscInt i;
232: if (N <= 3) {
233: PetscInt Njk;
235: PetscDTBinomialInt(N, j + k, &Njk);
236: if (!j) {
237: for (i = 0; i < Njk * Njk; i++) awedgeMat[i] = 0.;
238: for (i = 0; i < Njk; i++) awedgeMat[i * (Njk + 1)] = a[0];
239: } else if (!k) {
240: for (i = 0; i < Njk; i++) awedgeMat[i] = a[i];
241: } else {
242: if (N == 2) {
243: awedgeMat[0] = -a[1];
244: awedgeMat[1] = a[0];
245: } else {
246: if (j + k == 2) {
247: awedgeMat[0] = -a[1];
248: awedgeMat[1] = a[0];
249: awedgeMat[2] = 0.;
250: awedgeMat[3] = -a[2];
251: awedgeMat[4] = 0.;
252: awedgeMat[5] = a[0];
253: awedgeMat[6] = 0.;
254: awedgeMat[7] = -a[2];
255: awedgeMat[8] = a[1];
256: } else {
257: awedgeMat[0] = a[2];
258: awedgeMat[1] = -a[1];
259: awedgeMat[2] = a[0];
260: }
261: }
262: }
263: } else {
264: PetscInt Njk;
265: PetscInt Nk;
266: PetscInt JKj, i;
267: PetscInt *subset, *subsetjk, *subsetj, *subsetk;
269: PetscDTBinomialInt(N, k, &Nk);
270: PetscDTBinomialInt(N, j + k, &Njk);
271: PetscDTBinomialInt(j + k, j, &JKj);
272: PetscMalloc4(j + k, &subset, j + k, &subsetjk, j, &subsetj, k, &subsetk);
273: for (i = 0; i < Njk * Nk; i++) awedgeMat[i] = 0.;
274: for (i = 0; i < Njk; i++) {
275: PetscInt l;
277: PetscDTEnumSubset(N, j + k, i, subset);
278: for (l = 0; l < JKj; l++) {
279: PetscBool jkOdd;
280: PetscInt m, jInd, kInd;
282: PetscDTEnumSplit(j + k, j, l, subsetjk, &jkOdd);
283: for (m = 0; m < j; m++) subsetj[m] = subset[subsetjk[m]];
284: for (m = 0; m < k; m++) subsetk[m] = subset[subsetjk[j + m]];
285: PetscDTSubsetIndex(N, j, subsetj, &jInd);
286: PetscDTSubsetIndex(N, k, subsetk, &kInd);
287: awedgeMat[i * Nk + kInd] += jkOdd ? -a[jInd] : a[jInd];
288: }
289: }
290: PetscFree4(subset, subsetjk, subsetj, subsetk);
291: }
292: return 0;
293: }
295: /*@
296: PetscDTAltVPullback - Compute the pullback of a k-form under a linear transformation of the coordinate space
298: Input Parameters:
299: + N - the dimension of the origin vector space of the linear transformation, M >= 0
300: . M - the dimension of the image vector space of the linear transformation, N >= 0
301: . L - a linear transformation, an [M x N] matrix in row-major format
302: . k - the *signed* degree k of the |k|-form w, -(min(M,N)) <= k <= min(M,N). A negative form degree indicates that the pullback should be conjugated by
303: the Hodge star operator (see note).
304: - w - a |k|-form in the image space, size [M choose |k|]
306: Output Parameter:
307: . Lstarw - the pullback of w to a |k|-form in the origin space, size [N choose |k|]: (Lstarw)(v_1,...v_k) = w(L*v_1,...,L*v_k).
309: Level: intermediate
311: Note: negative form degrees accommodate, e.g., H-div conforming vector fields. An H-div conforming vector field stores its degrees of freedom as (dx, dy, dz), like a 1-form,
312: but its normal trace is integrated on faces, like a 2-form. The correct pullback then is to apply the Hodge star transformation from (M-2)-form to 2-form, pullback as a 2-form,
313: then the inverse Hodge star transformation.
315: .seealso: `PetscDTAltV`, `PetscDTAltVPullbackMatrix()`, `PetscDTAltVStar()`
316: @*/
317: PetscErrorCode PetscDTAltVPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *Lstarw)
318: {
319: PetscInt i, j, Nk, Mk;
323: if (N <= 3 && M <= 3) {
324: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
325: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
326: if (!k) {
327: Lstarw[0] = w[0];
328: } else if (k == 1) {
329: for (i = 0; i < Nk; i++) {
330: PetscReal sum = 0.;
332: for (j = 0; j < Mk; j++) sum += L[j * Nk + i] * w[j];
333: Lstarw[i] = sum;
334: }
335: } else if (k == -1) {
336: PetscReal mult[3] = {1., -1., 1.};
338: for (i = 0; i < Nk; i++) {
339: PetscReal sum = 0.;
341: for (j = 0; j < Mk; j++) sum += L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * w[j] * mult[j];
342: Lstarw[i] = mult[i] * sum;
343: }
344: } else if (k == 2) {
345: PetscInt pairs[3][2] = {
346: {0, 1},
347: {0, 2},
348: {1, 2}
349: };
351: for (i = 0; i < Nk; i++) {
352: PetscReal sum = 0.;
353: for (j = 0; j < Mk; j++) sum += (L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] - L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]]) * w[j];
354: Lstarw[i] = sum;
355: }
356: } else if (k == -2) {
357: PetscInt pairs[3][2] = {
358: {1, 2},
359: {2, 0},
360: {0, 1}
361: };
362: PetscInt offi = (N == 2) ? 2 : 0;
363: PetscInt offj = (M == 2) ? 2 : 0;
365: for (i = 0; i < Nk; i++) {
366: PetscReal sum = 0.;
368: for (j = 0; j < Mk; j++) sum += (L[pairs[offj + j][0] * N + pairs[offi + i][0]] * L[pairs[offj + j][1] * N + pairs[offi + i][1]] - L[pairs[offj + j][1] * N + pairs[offi + i][0]] * L[pairs[offj + j][0] * N + pairs[offi + i][1]]) * w[j];
369: Lstarw[i] = sum;
370: }
371: } else {
372: PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + L[1] * (L[5] * L[6] - L[3] * L[8]) + L[2] * (L[3] * L[7] - L[4] * L[6]);
374: for (i = 0; i < Nk; i++) Lstarw[i] = detL * w[i];
375: }
376: } else {
377: PetscInt Nf, l, p;
378: PetscReal *Lw, *Lwv;
379: PetscInt *subsetw, *subsetv;
380: PetscInt *perm;
381: PetscReal *walloc = NULL;
382: const PetscReal *ww = NULL;
383: PetscBool negative = PETSC_FALSE;
385: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
386: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
387: PetscDTFactorialInt(PetscAbsInt(k), &Nf);
388: if (k < 0) {
389: negative = PETSC_TRUE;
390: k = -k;
391: PetscMalloc1(Mk, &walloc);
392: PetscDTAltVStar(M, M - k, 1, w, walloc);
393: ww = walloc;
394: } else {
395: ww = w;
396: }
397: PetscMalloc5(k, &subsetw, k, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);
398: for (i = 0; i < Nk; i++) Lstarw[i] = 0.;
399: for (i = 0; i < Mk; i++) {
400: PetscDTEnumSubset(M, k, i, subsetw);
401: for (j = 0; j < Nk; j++) {
402: PetscDTEnumSubset(N, k, j, subsetv);
403: for (p = 0; p < Nf; p++) {
404: PetscReal prod;
405: PetscBool isOdd;
407: PetscDTEnumPerm(k, p, perm, &isOdd);
408: prod = isOdd ? -ww[i] : ww[i];
409: for (l = 0; l < k; l++) prod *= L[subsetw[perm[l]] * N + subsetv[l]];
410: Lstarw[j] += prod;
411: }
412: }
413: }
414: if (negative) {
415: PetscReal *sLsw;
417: PetscMalloc1(Nk, &sLsw);
418: PetscDTAltVStar(N, N - k, -1, Lstarw, sLsw);
419: for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i];
420: PetscFree(sLsw);
421: }
422: PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
423: PetscFree(walloc);
424: }
425: return 0;
426: }
428: /*@
429: PetscDTAltVPullbackMatrix - Compute the pullback matrix for k-forms under a linear transformation
431: Input Parameters:
432: + N - the dimension of the origin vector space of the linear transformation, N >= 0
433: . M - the dimension of the image vector space of the linear transformation, M >= 0
434: . L - a linear transformation, an [M x N] matrix in row-major format
435: - k - the *signed* degree k of the |k|-forms on which Lstar acts, -(min(M,N)) <= k <= min(M,N).
436: A negative form degree indicates that the pullback should be conjugated by the Hodge star operator (see note in `PetscDTAltvPullback()`)
438: Output Parameter:
439: . Lstar - the pullback matrix, an [(N choose |k|) x (M choose |k|)] matrix in row-major format such that Lstar * w = L^* w
441: Level: intermediate
443: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVStar()`
444: @*/
445: PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, PetscReal *Lstar)
446: {
447: PetscInt Nk, Mk, Nf, i, j, l, p;
448: PetscReal *Lw, *Lwv;
449: PetscInt *subsetw, *subsetv;
450: PetscInt *perm;
451: PetscBool negative = PETSC_FALSE;
455: if (N <= 3 && M <= 3) {
456: PetscReal mult[3] = {1., -1., 1.};
458: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
459: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
460: if (!k) {
461: Lstar[0] = 1.;
462: } else if (k == 1) {
463: for (i = 0; i < Nk; i++) {
464: for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[j * Nk + i];
465: }
466: } else if (k == -1) {
467: for (i = 0; i < Nk; i++) {
468: for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * mult[i] * mult[j];
469: }
470: } else if (k == 2) {
471: PetscInt pairs[3][2] = {
472: {0, 1},
473: {0, 2},
474: {1, 2}
475: };
477: for (i = 0; i < Nk; i++) {
478: for (j = 0; j < Mk; j++) Lstar[i * Mk + j] = L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] - L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]];
479: }
480: } else if (k == -2) {
481: PetscInt pairs[3][2] = {
482: {1, 2},
483: {2, 0},
484: {0, 1}
485: };
486: PetscInt offi = (N == 2) ? 2 : 0;
487: PetscInt offj = (M == 2) ? 2 : 0;
489: for (i = 0; i < Nk; i++) {
490: for (j = 0; j < Mk; j++) {
491: Lstar[i * Mk + j] = L[pairs[offj + j][0] * N + pairs[offi + i][0]] * L[pairs[offj + j][1] * N + pairs[offi + i][1]] - L[pairs[offj + j][1] * N + pairs[offi + i][0]] * L[pairs[offj + j][0] * N + pairs[offi + i][1]];
492: }
493: }
494: } else {
495: PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + L[1] * (L[5] * L[6] - L[3] * L[8]) + L[2] * (L[3] * L[7] - L[4] * L[6]);
497: for (i = 0; i < Nk; i++) Lstar[i] = detL;
498: }
499: } else {
500: if (k < 0) {
501: negative = PETSC_TRUE;
502: k = -k;
503: }
504: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
505: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
506: PetscDTFactorialInt(PetscAbsInt(k), &Nf);
507: PetscMalloc5(M, &subsetw, N, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);
508: for (i = 0; i < Nk * Mk; i++) Lstar[i] = 0.;
509: for (i = 0; i < Mk; i++) {
510: PetscBool iOdd;
511: PetscInt iidx, jidx;
513: PetscDTEnumSplit(M, k, i, subsetw, &iOdd);
514: iidx = negative ? Mk - 1 - i : i;
515: iOdd = negative ? (PetscBool)(iOdd ^ ((k * (M - k)) & 1)) : PETSC_FALSE;
516: for (j = 0; j < Nk; j++) {
517: PetscBool jOdd;
519: PetscDTEnumSplit(N, k, j, subsetv, &jOdd);
520: jidx = negative ? Nk - 1 - j : j;
521: jOdd = negative ? (PetscBool)(iOdd ^ jOdd ^ ((k * (N - k)) & 1)) : PETSC_FALSE;
522: for (p = 0; p < Nf; p++) {
523: PetscReal prod;
524: PetscBool isOdd;
526: PetscDTEnumPerm(k, p, perm, &isOdd);
527: isOdd = (PetscBool)(isOdd ^ jOdd);
528: prod = isOdd ? -1. : 1.;
529: for (l = 0; l < k; l++) prod *= L[subsetw[perm[l]] * N + subsetv[l]];
530: Lstar[jidx * Mk + iidx] += prod;
531: }
532: }
533: }
534: PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
535: }
536: return 0;
537: }
539: /*@
540: PetscDTAltVInterior - Compute the interior product of a k-form with a vector
542: Input Parameters:
543: + N - the dimension of the vector space, N >= 0
544: . k - the degree k of the k-form w, 0 <= k <= N
545: . w - a k-form, size [N choose k]
546: - v - an N dimensional vector
548: Output Parameter:
549: . wIntv - the (k-1)-form (w int v), size [N choose (k-1)]: (w int v) is defined by its action on (k-1) vectors {v_1, ..., v_{k-1}} as (w inv v)(v_1, ..., v_{k-1}) = w(v, v_1, ..., v_{k-1}).
551: Level: intermediate
553: .seealso: `PetscDTAltV`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
554: @*/
555: PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv)
556: {
557: PetscInt i, Nk, Nkm;
560: PetscDTBinomialInt(N, k, &Nk);
561: PetscDTBinomialInt(N, k - 1, &Nkm);
562: if (N <= 3) {
563: if (k == 1) {
564: PetscReal sum = 0.;
566: for (i = 0; i < N; i++) sum += w[i] * v[i];
567: wIntv[0] = sum;
568: } else if (k == N) {
569: PetscReal mult[3] = {1., -1., 1.};
571: for (i = 0; i < N; i++) wIntv[N - 1 - i] = w[0] * v[i] * mult[i];
572: } else {
573: wIntv[0] = -w[0] * v[1] - w[1] * v[2];
574: wIntv[1] = w[0] * v[0] - w[2] * v[2];
575: wIntv[2] = w[1] * v[0] + w[2] * v[1];
576: }
577: } else {
578: PetscInt *subset, *work;
580: PetscMalloc2(k, &subset, k, &work);
581: for (i = 0; i < Nkm; i++) wIntv[i] = 0.;
582: for (i = 0; i < Nk; i++) {
583: PetscInt j, l, m;
585: PetscDTEnumSubset(N, k, i, subset);
586: for (j = 0; j < k; j++) {
587: PetscInt idx;
588: PetscBool flip = (PetscBool)(j & 1);
590: for (l = 0, m = 0; l < k; l++) {
591: if (l != j) work[m++] = subset[l];
592: }
593: PetscDTSubsetIndex(N, k - 1, work, &idx);
594: wIntv[idx] += flip ? -(w[i] * v[subset[j]]) : (w[i] * v[subset[j]]);
595: }
596: }
597: PetscFree2(subset, work);
598: }
599: return 0;
600: }
602: /*@
603: PetscDTAltVInteriorMatrix - Compute the matrix of the linear transformation induced on a k-form by the interior product with a vector
605: Input Parameters:
606: + N - the dimension of the vector space, N >= 0
607: . k - the degree k of the k-forms on which intvMat acts, 0 <= k <= N
608: - v - an N dimensional vector
610: Output Parameter:
611: . intvMat - an [(N choose (k-1)) x (N choose k)] matrix, row-major: (intvMat) * w = (w int v)
613: Level: intermediate
615: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorPattern()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
616: @*/
617: PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat)
618: {
619: PetscInt i, Nk, Nkm;
622: PetscDTBinomialInt(N, k, &Nk);
623: PetscDTBinomialInt(N, k - 1, &Nkm);
624: if (N <= 3) {
625: if (k == 1) {
626: for (i = 0; i < N; i++) intvMat[i] = v[i];
627: } else if (k == N) {
628: PetscReal mult[3] = {1., -1., 1.};
630: for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i];
631: } else {
632: intvMat[0] = -v[1];
633: intvMat[1] = -v[2];
634: intvMat[2] = 0.;
635: intvMat[3] = v[0];
636: intvMat[4] = 0.;
637: intvMat[5] = -v[2];
638: intvMat[6] = 0.;
639: intvMat[7] = v[0];
640: intvMat[8] = v[1];
641: }
642: } else {
643: PetscInt *subset, *work;
645: PetscMalloc2(k, &subset, k, &work);
646: for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.;
647: for (i = 0; i < Nk; i++) {
648: PetscInt j, l, m;
650: PetscDTEnumSubset(N, k, i, subset);
651: for (j = 0; j < k; j++) {
652: PetscInt idx;
653: PetscBool flip = (PetscBool)(j & 1);
655: for (l = 0, m = 0; l < k; l++) {
656: if (l != j) work[m++] = subset[l];
657: }
658: PetscDTSubsetIndex(N, k - 1, work, &idx);
659: intvMat[idx * Nk + i] += flip ? -v[subset[j]] : v[subset[j]];
660: }
661: }
662: PetscFree2(subset, work);
663: }
664: return 0;
665: }
667: /*@
668: PetscDTAltVInteriorPattern - compute the sparsity and sign pattern of the interior product matrix computed in `PetscDTAltVInteriorMatrix()`
670: Input Parameters:
671: + N - the dimension of the vector space, N >= 0
672: - k - the degree of the k-forms on which intvMat from `PetscDTAltVInteriorMatrix()` acts, 0 <= k <= N.
674: Output Parameter:
675: . indices - The interior product matrix intvMat has size [(N choose (k-1)) x (N choose k)] and has (N choose k) * k
676: non-zeros. indices[i][0] and indices[i][1] are the row and column of a non-zero, and its value is equal to the vector
677: coordinate v[j] if indices[i][2] = j, or -v[j] if indices[i][2] = -(j+1)
679: Level: intermediate
681: Note:
682: This function is useful when the interior product needs to be computed at multiple locations, as when computing the Koszul differential
684: .seealso: `PetscDTAltV`, `PetscDTAltVInterior()`, `PetscDTAltVInteriorMatrix()`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
685: @*/
686: PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3])
687: {
688: PetscInt i, Nk, Nkm;
691: PetscDTBinomialInt(N, k, &Nk);
692: PetscDTBinomialInt(N, k - 1, &Nkm);
693: if (N <= 3) {
694: if (k == 1) {
695: for (i = 0; i < N; i++) {
696: indices[i][0] = 0;
697: indices[i][1] = i;
698: indices[i][2] = i;
699: }
700: } else if (k == N) {
701: PetscInt val[3] = {0, -2, 2};
703: for (i = 0; i < N; i++) {
704: indices[i][0] = N - 1 - i;
705: indices[i][1] = 0;
706: indices[i][2] = val[i];
707: }
708: } else {
709: indices[0][0] = 0;
710: indices[0][1] = 0;
711: indices[0][2] = -(1 + 1);
712: indices[1][0] = 0;
713: indices[1][1] = 1;
714: indices[1][2] = -(2 + 1);
715: indices[2][0] = 1;
716: indices[2][1] = 0;
717: indices[2][2] = 0;
718: indices[3][0] = 1;
719: indices[3][1] = 2;
720: indices[3][2] = -(2 + 1);
721: indices[4][0] = 2;
722: indices[4][1] = 1;
723: indices[4][2] = 0;
724: indices[5][0] = 2;
725: indices[5][1] = 2;
726: indices[5][2] = 1;
727: }
728: } else {
729: PetscInt *subset, *work;
731: PetscMalloc2(k, &subset, k, &work);
732: for (i = 0; i < Nk; i++) {
733: PetscInt j, l, m;
735: PetscDTEnumSubset(N, k, i, subset);
736: for (j = 0; j < k; j++) {
737: PetscInt idx;
738: PetscBool flip = (PetscBool)(j & 1);
740: for (l = 0, m = 0; l < k; l++) {
741: if (l != j) work[m++] = subset[l];
742: }
743: PetscDTSubsetIndex(N, k - 1, work, &idx);
744: indices[i * k + j][0] = idx;
745: indices[i * k + j][1] = i;
746: indices[i * k + j][2] = flip ? -(subset[j] + 1) : subset[j];
747: }
748: }
749: PetscFree2(subset, work);
750: }
751: return 0;
752: }
754: /*@
755: PetscDTAltVStar - Apply a power of the Hodge star operator, which maps k-forms to (N-k) forms, to a k-form
757: Input Parameters:
758: + N - the dimension of the vector space, N >= 0
759: . k - the degree k of the k-form w, 0 <= k <= N
760: . pow - the number of times to apply the Hodge star operator: pow < 0 indicates that the inverse of the Hodge star operator should be applied |pow| times.
761: - w - a k-form, size [N choose k]
763: Output Parameter:
764: . starw = (star)^pow w. Each degree of freedom of a k-form is associated with a subset S of k coordinates of the N dimensional vector space: the Hodge start operator (star) maps that degree of freedom to the degree of freedom associated with S', the complement of S, with a sign change if the permutation of coordinates {S[0], ... S[k-1], S'[0], ... S'[N-k- 1]} is an odd permutation. This implies (star)^2 w = (-1)^{k(N-k)} w, and (star)^4 w = w.
766: Level: intermediate
768: .seealso: `PetscDTAltV`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
769: @*/
770: PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw)
771: {
772: PetscInt Nk, i;
775: PetscDTBinomialInt(N, k, &Nk);
776: pow = pow % 4;
777: pow = (pow + 4) % 4; /* make non-negative */
778: /* pow is now 0, 1, 2, 3 */
779: if (N <= 3) {
780: if (pow & 1) {
781: PetscReal mult[3] = {1., -1., 1.};
783: for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i];
784: } else {
785: for (i = 0; i < Nk; i++) starw[i] = w[i];
786: }
787: if (pow > 1 && ((k * (N - k)) & 1)) {
788: for (i = 0; i < Nk; i++) starw[i] = -starw[i];
789: }
790: } else {
791: PetscInt *subset;
793: PetscMalloc1(N, &subset);
794: if (pow % 2) {
795: PetscInt l = (pow == 1) ? k : N - k;
796: for (i = 0; i < Nk; i++) {
797: PetscBool sOdd;
798: PetscInt j, idx;
800: PetscDTEnumSplit(N, l, i, subset, &sOdd);
801: PetscDTSubsetIndex(N, l, subset, &idx);
802: PetscDTSubsetIndex(N, N - l, &subset[l], &j);
803: starw[j] = sOdd ? -w[idx] : w[idx];
804: }
805: } else {
806: for (i = 0; i < Nk; i++) starw[i] = w[i];
807: }
808: /* star^2 = -1^(k * (N - k)) */
809: if (pow > 1 && (k * (N - k)) % 2) {
810: for (i = 0; i < Nk; i++) starw[i] = -starw[i];
811: }
812: PetscFree(subset);
813: }
814: return 0;
815: }