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