Actual source code: dtaltv.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/petscimpl.h>
2: #include <petsc/private/dtimpl.h>
4: /*MC
5: PetscDTAltV - An interface for common operations on k-forms, also known as alternating algebraic forms or alternating k-linear maps.
6: The name of the interface comes from the notation "Alt V" for the algebra of all k-forms acting vectors in the space V, also known as the exterior algebra of V*.
8: A recommended reference for this material is Section 2 "Exterior algebra and exterior calculus" in "Finite element
9: exterior calculus, homological techniques, and applications", by Arnold, Falk, & Winther (2006, doi:10.1017/S0962492906210018).
11: A k-form w (k is called the "form degree" of w) is an alternating k-linear map acting on tuples (v_1, ..., v_k) of
12: vectors from a vector space V and producing a real number:
13: - alternating: swapping any two vectors in a tuple reverses the sign of the result, e.g. w(v_1, v_2, ..., v_k) = -w(v_2, v_1, ..., v_k)
14: - k-linear: w acts linear in each vector separately, e.g. w(a*v + b*y, v_2, ..., v_k) = a*w(v,v_2,...,v_k) + b*w(y,v_2,...,v_k)
15: This action is implemented as PetscDTAltVApply.
17: The k-forms on a vector space form a vector space themselves, Alt^k V. The dimension of Alt^k V, if V is N dimensional, is N choose k. (This
18: shows that for an N dimensional space, only 0 <= k <= N are valid form degrees.)
19: The standard basis for Alt^k V, used in PetscDTAltV, has one basis k-form for each ordered subset of k coordinates of the N dimensional space:
20: For example, if the coordinate directions of a four dimensional space are (t, x, y, z), then there are 4 choose 2 = 6 ordered subsets of two coordinates.
21: They are, in lexicographic order, (t, x), (t, y), (t, z), (x, y), (x, z) and (y, z). PetscDTAltV also orders the basis of Alt^k V lexicographically
22: by the associated subsets.
24: The unit basis k-form associated with coordinates (c_1, ..., c_k) acts on a set of k vectors (v_1, ..., v_k) by creating a square matrix V where
25: V[i,j] = v_i[c_j] and taking the determinant of V.
27: If j + k <= N, then a j-form f and a k-form g can be multiplied to create a (j+k)-form using the wedge or exterior product, (f wedge g).
28: This is an anticommutative product, (f wedge g) = -(g wedge f). It is sufficient to describe the wedge product of two basis forms.
29: Let f be the basis j-form associated with coordinates (f_1,...,f_j) and g be the basis k-form associated with coordinates (g_1,...,g_k):
30: - If there is any coordinate in both sets, then (f wedge g) = 0.
31: - Otherwise, (f wedge g) is a multiple of the basis (j+k)-form h associated with (f_1,...,f_j,g_1,...,g_k).
32: - In fact it is equal to either h or -h depending on how (f_1,...,f_j,g_1,...,g_k) compares to the same list of coordinates given in ascending order: if it is an even permutation of that list, then (f wedge g) = h, otherwise (f wedge g) = -h.
33: The wedge product is implemented for either two inputs (f and g) in PetscDTAltVWedge, or for one (just f, giving a
34: matrix to multiply against multiple choices of g) in PetscDTAltVWedgeMatrix.
36: If k > 0, a k-form w and a vector v can combine to make a (k-1)-formm through the interior product, (w int v),
37: defined by (w int v)(v_1,...,v_{k-1}) = w(v,v_1,...,v_{k-1}).
39: The interior product is implemented for either two inputs (w and v) in PetscDTAltVInterior, for one (just v, giving a
40: matrix to multiply against multiple choices of w) in PetscDTAltVInteriorMatrix,
41: or for no inputs (giving the sparsity pattern of PetscDTAltVInteriorMatrix) in PetscDTAltVInteriorPattern.
43: When there is a linear map L: V -> W from an N dimensional vector space to an M dimensional vector space,
44: it induces the linear pullback map L^* : Alt^k W -> Alt^k V, defined by L^* w(v_1,...,v_k) = w(L v_1, ..., L v_k).
45: The pullback is implemented as PetscDTAltVPullback (acting on a known w) and PetscDTAltVPullbackMatrix (creating a matrix that computes the actin of L^*).
47: Alt^k V and Alt^(N-k) V have the same dimension, and the Hodge star operator maps between them. We note that Alt^N V is a one dimensional space, and its
48: basis vector is sometime called vol. The Hodge star operator has the property that (f wedge (star g)) = (f,g) vol, where (f,g) is the simple inner product
49: of the basis coefficients of f and g.
50: Powers of the Hodge star operator can be applied with PetscDTAltVStar
52: Level: intermediate
54: .seealso: PetscDTAltVApply(), PetscDTAltVWedge(), PetscDTAltVInterior(), PetscDTAltVPullback(), PetscDTAltVStar()
55: M*/
57: /*@
58: PetscDTAltVApply - Apply an a k-form (an alternating k-linear map) to a set of k N-dimensional vectors
60: Input Arguments:
61: + N - the dimension of the vector space, N >= 0
62: . k - the degree k of the k-form w, 0 <= k <= N
63: . w - a k-form, size [N choose k] (each degree of freedom of a k-form is associated with a subset of k coordinates of the N-dimensional vectors: the degrees of freedom are ordered lexicographically by their associated subsets)
64: - v - a set of k vectors of size N, size [k x N], each vector stored contiguously
66: Output Arguments:
67: . wv - w(v_1,...,v_k) = \sum_i w_i * det(V_i): the degree of freedom w_i is associated with coordinates [s_{i,1},...,s_{i,k}], and the square matrix V_i has entry (j,k) given by the s_{i,k}'th coordinate of v_j
69: Level: intermediate
71: .seealso: PetscDTAltV, PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
72: @*/
73: PetscErrorCode PetscDTAltVApply(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wv)
74: {
78: if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
79: if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
80: if (N <= 3) {
81: if (!k) {
82: *wv = w[0];
83: } else {
84: if (N == 1) {*wv = w[0] * v[0];}
85: else if (N == 2) {
86: if (k == 1) {*wv = w[0] * v[0] + w[1] * v[1];}
87: else {*wv = w[0] * (v[0] * v[3] - v[1] * v[2]);}
88: } else {
89: if (k == 1) {*wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];}
90: else if (k == 2) {
91: *wv = w[0] * (v[0] * v[4] - v[1] * v[3]) +
92: w[1] * (v[0] * v[5] - v[2] * v[3]) +
93: 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]) +
96: v[1] * (v[5] * v[6] - v[3] * v[8]) +
97: 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++) {
121: prod *= v[perm[l] * N + subset[l]];
122: }
123: subsum += prod;
124: }
125: sum += w[i] * subsum;
126: }
127: PetscFree2(subset, perm);
128: *wv = sum;
129: }
130: return(0);
131: }
133: /*@
134: PetscDTAltVWedge - Compute the wedge product of a j-form and a k-form, giving a (j+k) form
136: Input Arguments:
137: + N - the dimension of the vector space, N >= 0
138: . j - the degree j of the j-form a, 0 <= j <= N
139: . k - the degree k of the k-form b, 0 <= k <= N and 0 <= j+k <= N
140: . a - a j-form, size [N choose j]
141: - b - a k-form, size [N choose k]
143: Output Arguments:
144: . awedgeb - the (j+k)-form a wedge b, size [N choose (j+k)]: (a wedge b)(v_1,...,v_{j+k}) = \sum_{s} sign(s) a(v_{s_1},...,v_{s_j}) b(v_{s_{j+1}},...,v_{s_{j+k}}),
145: where the sum is over permutations s such that s_1 < s_2 < ... < s_j and s_{j+1} < s_{j+2} < ... < s_{j+k}.
147: Level: intermediate
149: .seealso: PetscDTAltV, PetscDTAltVWedgeMatrix(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
150: @*/
151: PetscErrorCode PetscDTAltVWedge(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, const PetscReal *b, PetscReal *awedgeb)
152: {
153: PetscInt i;
157: if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
158: if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree");
159: if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension");
160: if (N <= 3) {
161: PetscInt Njk;
163: PetscDTBinomialInt(N, j+k, &Njk);
164: if (!j) {for (i = 0; i < Njk; i++) {awedgeb[i] = a[0] * b[i];}}
165: else if (!k) {for (i = 0; i < Njk; i++) {awedgeb[i] = a[i] * b[0];}}
166: else {
167: if (N == 2) {awedgeb[0] = a[0] * b[1] - a[1] * b[0];}
168: else {
169: if (j+k == 2) {
170: awedgeb[0] = a[0] * b[1] - a[1] * b[0];
171: awedgeb[1] = a[0] * b[2] - a[2] * b[0];
172: awedgeb[2] = a[1] * b[2] - a[2] * b[1];
173: } else {
174: awedgeb[0] = a[0] * b[2] - a[1] * b[1] + a[2] * b[0];
175: }
176: }
177: }
178: } else {
179: PetscInt Njk;
180: PetscInt JKj;
181: PetscInt *subset, *subsetjk, *subsetj, *subsetk;
182: PetscInt i;
184: PetscDTBinomialInt(N, j+k, &Njk);
185: PetscDTBinomialInt(j+k, j, &JKj);
186: PetscMalloc4(j+k, &subset, j+k, &subsetjk, j, &subsetj, k, &subsetk);
187: for (i = 0; i < Njk; i++) {
188: PetscReal sum = 0.;
189: PetscInt l;
191: PetscDTEnumSubset(N, j+k, i, subset);
192: for (l = 0; l < JKj; l++) {
193: PetscBool jkOdd;
194: PetscInt m, jInd, kInd;
196: PetscDTEnumSplit(j+k, j, l, subsetjk, &jkOdd);
197: for (m = 0; m < j; m++) {
198: subsetj[m] = subset[subsetjk[m]];
199: }
200: for (m = 0; m < k; m++) {
201: subsetk[m] = subset[subsetjk[j+m]];
202: }
203: PetscDTSubsetIndex(N, j, subsetj, &jInd);
204: PetscDTSubsetIndex(N, k, subsetk, &kInd);
205: sum += jkOdd ? -(a[jInd] * b[kInd]) : (a[jInd] * b[kInd]);
206: }
207: awedgeb[i] = sum;
208: }
209: PetscFree4(subset, subsetjk, subsetj, subsetk);
210: }
211: return(0);
212: }
214: /*@
215: PetscDTAltVWedgeMatrix - Compute the matrix defined by the wedge product with a given j-form that maps k-forms to (j+k)-forms
217: Input Arguments:
218: + N - the dimension of the vector space, N >= 0
219: . j - the degree j of the j-form a, 0 <= j <= N
220: . k - the degree k of the k-forms that (a wedge) will be applied to, 0 <= k <= N and 0 <= j+k <= N
221: - a - a j-form, size [N choose j]
223: Output Arguments:
224: . awedge - (a wedge), an [(N choose j+k) x (N choose k)] matrix in row-major order, such that (a wedge) * b = a wedge b
226: Level: intermediate
228: .seealso: PetscDTAltV, PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
229: @*/
230: PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, PetscReal *awedgeMat)
231: {
232: PetscInt i;
236: if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension");
237: if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree");
238: if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension");
239: if (N <= 3) {
240: PetscInt Njk;
242: PetscDTBinomialInt(N, j+k, &Njk);
243: if (!j) {
244: for (i = 0; i < Njk * Njk; i++) {awedgeMat[i] = 0.;}
245: for (i = 0; i < Njk; i++) {awedgeMat[i * (Njk + 1)] = a[0];}
246: } else if (!k) {
247: for (i = 0; i < Njk; i++) {awedgeMat[i] = a[i];}
248: } else {
249: if (N == 2) {
250: awedgeMat[0] = -a[1]; awedgeMat[1] = a[0];
251: } else {
252: if (j+k == 2) {
253: awedgeMat[0] = -a[1]; awedgeMat[1] = a[0]; awedgeMat[2] = 0.;
254: awedgeMat[3] = -a[2]; awedgeMat[4] = 0.; awedgeMat[5] = a[0];
255: awedgeMat[6] = 0.; awedgeMat[7] = -a[2]; awedgeMat[8] = a[1];
256: } else {
257: awedgeMat[0] = a[2]; awedgeMat[1] = -a[1]; 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++) {
282: subsetj[m] = subset[subsetjk[m]];
283: }
284: for (m = 0; m < k; m++) {
285: subsetk[m] = subset[subsetjk[j+m]];
286: }
287: PetscDTSubsetIndex(N, j, subsetj, &jInd);
288: PetscDTSubsetIndex(N, k, subsetk, &kInd);
289: awedgeMat[i * Nk + kInd] += jkOdd ? - a[jInd] : a[jInd];
290: }
291: }
292: PetscFree4(subset, subsetjk, subsetj, subsetk);
293: }
294: return(0);
295: }
297: /*@
298: PetscDTAltVPullback - Compute the pullback of a k-form under a linear transformation of the coordinate space
300: Input Arguments:
301: + N - the dimension of the origin vector space of the linear transformation, M >= 0
302: . M - the dimension of the image vector space of the linear transformation, N >= 0
303: . L - a linear transformation, an [M x N] matrix in row-major format
304: . k - the *signed* degree k of the |k|-form w, -(min(M,N)) <= k <= min(M,N). A negative form degree indicates that the pullback should be conjugated by the Hodge star operator (see note).
305: - w - a |k|-form in the image space, size [M choose |k|]
307: Output Arguments:
308: . Lstarw - the pullback of w to a |k|-form in the origin space, size [N choose |k|]: (Lstarw)(v_1,...v_k) = w(L*v_1,...,L*v_k).
310: Level: intermediate
312: Note: negative form degrees accomodate, e.g., H-div conforming vector fields. An H-div conforming vector field stores its degrees of freedom as (dx, dy, dz), like a 1-form,
313: but its normal trace is integrated on faces, like a 2-form. The correct pullback then is to apply the Hodge star transformation from (M-2)-form to 2-form, pullback as a 2-form,
314: then the inverse Hodge star transformation.
316: .seealso: PetscDTAltV, PetscDTAltVPullbackMatrix(), PetscDTAltVStar()
317: @*/
318: PetscErrorCode PetscDTAltVPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *Lstarw)
319: {
320: PetscInt i, j, Nk, Mk;
321: PetscErrorCode ierr;
324: if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
325: if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
326: if (N <= 3 && M <= 3) {
328: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
329: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
330: if (!k) {
331: Lstarw[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++) {
346: sum += L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * w[j] * mult[j];
347: }
348: Lstarw[i] = mult[i] * sum;
349: }
350: } else if (k == 2) {
351: PetscInt pairs[3][2] = {{0,1},{0,2},{1,2}};
353: for (i = 0; i < Nk; i++) {
354: PetscReal sum = 0.;
355: for (j = 0; j < Mk; j++) {
356: sum += (L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] -
357: L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]]) * w[j];
358: }
359: Lstarw[i] = sum;
360: }
361: } else if (k == -2) {
362: PetscInt pairs[3][2] = {{1,2},{2,0},{0,1}};
363: PetscInt offi = (N == 2) ? 2 : 0;
364: PetscInt offj = (M == 2) ? 2 : 0;
366: for (i = 0; i < Nk; i++) {
367: PetscReal sum = 0.;
369: for (j = 0; j < Mk; j++) {
370: sum += (L[pairs[offj + j][0] * N + pairs[offi + i][0]] *
371: L[pairs[offj + j][1] * N + pairs[offi + i][1]] -
372: L[pairs[offj + j][1] * N + pairs[offi + i][0]] *
373: L[pairs[offj + j][0] * N + pairs[offi + i][1]]) * w[j];
375: }
376: Lstarw[i] = sum;
377: }
378: } else {
379: PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) +
380: L[1] * (L[5] * L[6] - L[3] * L[8]) +
381: L[2] * (L[3] * L[7] - L[4] * L[6]);
383: for (i = 0; i < Nk; i++) {Lstarw[i] = detL * w[i];}
384: }
385: } else {
386: PetscInt Nf, l, p;
387: PetscReal *Lw, *Lwv;
388: PetscInt *subsetw, *subsetv;
389: PetscInt *perm;
390: PetscReal *walloc = NULL;
391: const PetscReal *ww = NULL;
392: PetscBool negative = PETSC_FALSE;
394: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
395: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
396: PetscDTFactorialInt(PetscAbsInt(k), &Nf);
397: if (k < 0) {
398: negative = PETSC_TRUE;
399: k = -k;
400: PetscMalloc1(Mk, &walloc);
401: PetscDTAltVStar(M, M - k, 1, w, walloc);
402: ww = walloc;
403: } else {
404: ww = w;
405: }
406: PetscMalloc5(k, &subsetw, k, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);
407: for (i = 0; i < Nk; i++) Lstarw[i] = 0.;
408: for (i = 0; i < Mk; i++) {
409: PetscDTEnumSubset(M, k, i, subsetw);
410: for (j = 0; j < Nk; j++) {
411: PetscDTEnumSubset(N, k, j, subsetv);
412: for (p = 0; p < Nf; p++) {
413: PetscReal prod;
414: PetscBool isOdd;
416: PetscDTEnumPerm(k, p, perm, &isOdd);
417: prod = isOdd ? -ww[i] : ww[i];
418: for (l = 0; l < k; l++) {
419: prod *= L[subsetw[perm[l]] * N + subsetv[l]];
420: }
421: Lstarw[j] += prod;
422: }
423: }
424: }
425: if (negative) {
426: PetscReal *sLsw;
428: PetscMalloc1(Nk, &sLsw);
429: PetscDTAltVStar(N, N - k, -1, Lstarw, sLsw);
430: for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i];
431: PetscFree(sLsw);
432: }
433: PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
434: PetscFree(walloc);
435: }
436: return(0);
437: }
439: /*@
440: PetscDTAltVPullbackMatrix - Compute the pullback matrix for k-forms under a linear transformation
442: Input Arguments:
443: + N - the dimension of the origin vector space of the linear transformation, N >= 0
444: . M - the dimension of the image vector space of the linear transformation, M >= 0
445: . L - a linear transformation, an [M x N] matrix in row-major format
446: - k - the *signed* degree k of the |k|-forms on which Lstar acts, -(min(M,N)) <= k <= min(M,N). A negative form degree indicates that the pullback should be conjugated by the Hodge star operator (see note in PetscDTAltvPullback())
448: Output Arguments:
449: . Lstar - the pullback matrix, an [(N choose |k|) x (M choose |k|)] matrix in row-major format such that Lstar * w = L^* w
451: Level: intermediate
453: .seealso: PetscDTAltV, PetscDTAltVPullback(), PetscDTAltVStar()
454: @*/
455: PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, PetscReal *Lstar)
456: {
457: PetscInt Nk, Mk, Nf, i, j, l, p;
458: PetscReal *Lw, *Lwv;
459: PetscInt *subsetw, *subsetv;
460: PetscInt *perm;
461: PetscBool negative = PETSC_FALSE;
462: PetscErrorCode ierr;
465: if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions");
466: if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
467: if (N <= 3 && M <= 3) {
468: PetscReal mult[3] = {1., -1., 1.};
470: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
471: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
472: if (!k) {
473: Lstar[0] = 1.;
474: } else if (k == 1) {
475: for (i = 0; i < Nk; i++) {for (j = 0; j < Mk; j++) {Lstar[i * Mk + j] = L[j * Nk + i];}}
476: } else if (k == -1) {
477: for (i = 0; i < Nk; i++) {
478: for (j = 0; j < Mk; j++) {
479: Lstar[i * Mk + j] = L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * mult[i] * mult[j];
480: }
481: }
482: } else if (k == 2) {
483: PetscInt pairs[3][2] = {{0,1},{0,2},{1,2}};
485: for (i = 0; i < Nk; i++) {
486: for (j = 0; j < Mk; j++) {
487: Lstar[i * Mk + j] = L[pairs[j][0] * N + pairs[i][0]] *
488: L[pairs[j][1] * N + pairs[i][1]] -
489: L[pairs[j][1] * N + pairs[i][0]] *
490: L[pairs[j][0] * N + pairs[i][1]];
491: }
492: }
493: } else if (k == -2) {
494: PetscInt pairs[3][2] = {{1,2},{2,0},{0,1}};
495: PetscInt offi = (N == 2) ? 2 : 0;
496: PetscInt offj = (M == 2) ? 2 : 0;
498: for (i = 0; i < Nk; i++) {
499: for (j = 0; j < Mk; j++) {
500: Lstar[i * Mk + j] = L[pairs[offj + j][0] * N + pairs[offi + i][0]] *
501: L[pairs[offj + j][1] * N + pairs[offi + i][1]] -
502: L[pairs[offj + j][1] * N + pairs[offi + i][0]] *
503: L[pairs[offj + j][0] * N + pairs[offi + i][1]];
504: }
505: }
506: } else {
507: PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) +
508: L[1] * (L[5] * L[6] - L[3] * L[8]) +
509: L[2] * (L[3] * L[7] - L[4] * L[6]);
511: for (i = 0; i < Nk; i++) {Lstar[i] = detL;}
512: }
513: } else {
514: if (k < 0) {
515: negative = PETSC_TRUE;
516: k = -k;
517: }
518: PetscDTBinomialInt(M, PetscAbsInt(k), &Mk);
519: PetscDTBinomialInt(N, PetscAbsInt(k), &Nk);
520: PetscDTFactorialInt(PetscAbsInt(k), &Nf);
521: PetscMalloc5(M, &subsetw, N, &subsetv, k, &perm, N * k, &Lw, k * k, &Lwv);
522: for (i = 0; i < Nk * Mk; i++) Lstar[i] = 0.;
523: for (i = 0; i < Mk; i++) {
524: PetscBool iOdd;
525: PetscInt iidx, jidx;
527: PetscDTEnumSplit(M, k, i, subsetw, &iOdd);
528: iidx = negative ? Mk - 1 - i : i;
529: iOdd = negative ? (PetscBool) (iOdd ^ ((k * (M-k)) & 1)) : PETSC_FALSE;
530: for (j = 0; j < Nk; j++) {
531: PetscBool jOdd;
533: PetscDTEnumSplit(N, k, j, subsetv, &jOdd);
534: jidx = negative ? Nk - 1 - j : j;
535: jOdd = negative ? (PetscBool) (iOdd ^ jOdd ^ ((k * (N-k)) & 1)) : PETSC_FALSE;
536: for (p = 0; p < Nf; p++) {
537: PetscReal prod;
538: PetscBool isOdd;
540: PetscDTEnumPerm(k, p, perm, &isOdd);
541: isOdd = (PetscBool) (isOdd ^ jOdd);
542: prod = isOdd ? -1. : 1.;
543: for (l = 0; l < k; l++) {
544: prod *= L[subsetw[perm[l]] * N + subsetv[l]];
545: }
546: Lstar[jidx * Mk + iidx] += prod;
547: }
548: }
549: }
550: PetscFree5(subsetw, subsetv, perm, Lw, Lwv);
551: }
552: return(0);
553: }
555: /*@
556: PetscDTAltVInterior - Compute the interior product of a k-form with a vector
558: Input Arguments:
559: + N - the dimension of the vector space, N >= 0
560: . k - the degree k of the k-form w, 0 <= k <= N
561: . w - a k-form, size [N choose k]
562: - v - an N dimensional vector
564: Output Arguments:
565: . wIntv - the (k-1)-form (w int v), size [N choose (k-1)]: (w int v) is defined by its action on (k-1) vectors {v_1, ..., v_{k-1}} as (w inv v)(v_1, ..., v_{k-1}) = w(v, v_1, ..., v_{k-1}).
567: Level: intermediate
569: .seealso: PetscDTAltV, PetscDTAltVInteriorMatrix(), PetscDTAltVInteriorPattern(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
570: @*/
571: PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv)
572: {
573: PetscInt i, Nk, Nkm;
574: PetscErrorCode ierr;
577: if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
578: PetscDTBinomialInt(N, k, &Nk);
579: PetscDTBinomialInt(N, k-1, &Nkm);
580: if (N <= 3) {
581: if (k == 1) {
582: PetscReal sum = 0.;
584: for (i = 0; i < N; i++) {
585: sum += w[i] * v[i];
586: }
587: wIntv[0] = sum;
588: } else if (k == N) {
589: PetscReal mult[3] = {1., -1., 1.};
591: for (i = 0; i < N; i++) {
592: wIntv[N - 1 - i] = w[0] * v[i] * mult[i];
593: }
594: } else {
595: wIntv[0] = - w[0]*v[1] - w[1]*v[2];
596: wIntv[1] = w[0]*v[0] - w[2]*v[2];
597: wIntv[2] = w[1]*v[0] + w[2]*v[1];
598: }
599: } else {
600: PetscInt *subset, *work;
602: PetscMalloc2(k, &subset, k, &work);
603: for (i = 0; i < Nkm; i++) wIntv[i] = 0.;
604: for (i = 0; i < Nk; i++) {
605: PetscInt j, l, m;
607: PetscDTEnumSubset(N, k, i, subset);
608: for (j = 0; j < k; j++) {
609: PetscInt idx;
610: PetscBool flip = (PetscBool) (j & 1);
612: for (l = 0, m = 0; l < k; l++) {
613: if (l != j) work[m++] = subset[l];
614: }
615: PetscDTSubsetIndex(N, k - 1, work, &idx);
616: wIntv[idx] += flip ? -(w[i] * v[subset[j]]) : (w[i] * v[subset[j]]);
617: }
618: }
619: PetscFree2(subset, work);
620: }
621: return(0);
622: }
624: /*@
625: PetscDTAltVInteriorMatrix - Compute the matrix of the linear transformation induced on a k-form by the interior product with a vector
627: Input Arguments:
628: + N - the dimension of the vector space, N >= 0
629: . k - the degree k of the k-forms on which intvMat acts, 0 <= k <= N
630: - v - an N dimensional vector
632: Output Arguments:
633: . intvMat - an [(N choose (k-1)) x (N choose k)] matrix, row-major: (intvMat) * w = (w int v)
635: Level: intermediate
637: .seealso: PetscDTAltV, PetscDTAltVInterior(), PetscDTAltVInteriorPattern(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
638: @*/
639: PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat)
640: {
641: PetscInt i, Nk, Nkm;
642: PetscErrorCode ierr;
645: if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
646: PetscDTBinomialInt(N, k, &Nk);
647: PetscDTBinomialInt(N, k-1, &Nkm);
648: if (N <= 3) {
649: if (k == 1) {
650: for (i = 0; i < N; i++) intvMat[i] = v[i];
651: } else if (k == N) {
652: PetscReal mult[3] = {1., -1., 1.};
654: for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i];
655: } else {
656: intvMat[0] = -v[1]; intvMat[1] = -v[2]; intvMat[2] = 0.;
657: intvMat[3] = v[0]; intvMat[4] = 0.; intvMat[5] = -v[2];
658: intvMat[6] = 0.; intvMat[7] = v[0]; intvMat[8] = v[1];
659: }
660: } else {
661: PetscInt *subset, *work;
663: PetscMalloc2(k, &subset, k, &work);
664: for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.;
665: for (i = 0; i < Nk; i++) {
666: PetscInt j, l, m;
668: PetscDTEnumSubset(N, k, i, subset);
669: for (j = 0; j < k; j++) {
670: PetscInt idx;
671: PetscBool flip = (PetscBool) (j & 1);
673: for (l = 0, m = 0; l < k; l++) {
674: if (l != j) work[m++] = subset[l];
675: }
676: PetscDTSubsetIndex(N, k - 1, work, &idx);
677: intvMat[idx * Nk + i] += flip ? -v[subset[j]] : v[subset[j]];
678: }
679: }
680: PetscFree2(subset, work);
681: }
682: return(0);
683: }
685: /*@
686: PetscDTAltVInteriorPattern - compute the sparsity and sign pattern of the interior product matrix computed in PetscDTAltVInteriorMatrix()
688: Input Arguments:
689: + N - the dimension of the vector space, N >= 0
690: - k - the degree of the k-forms on which intvMat from PetscDTAltVInteriorMatrix() acts, 0 <= k <= N.
692: Output Arguments:
693: . indices - The interior product matrix intvMat has size [(N choose (k-1)) x (N choose k)] and has (N choose k) * k
694: non-zeros. indices[i][0] and indices[i][1] are the row and column of a non-zero, and its value is equal to the vector
695: coordinate v[j] if indices[i][2] = j, or -v[j] if indices[i][2] = -(j+1)
697: Level: intermediate
699: Note: this function is useful when the interior product needs to be computed at multiple locations, as when computing the Koszul differential
701: .seealso: PetscDTAltV, PetscDTAltVInterior(), PetscDTAltVInteriorMatrix(), PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
702: @*/
703: PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3])
704: {
705: PetscInt i, Nk, Nkm;
706: PetscErrorCode ierr;
709: if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
710: PetscDTBinomialInt(N, k, &Nk);
711: PetscDTBinomialInt(N, k-1, &Nkm);
712: if (N <= 3) {
713: if (k == 1) {
714: for (i = 0; i < N; i++) {
715: indices[i][0] = 0;
716: indices[i][1] = i;
717: indices[i][2] = i;
718: }
719: } else if (k == N) {
720: PetscInt val[3] = {0, -2, 2};
722: for (i = 0; i < N; i++) {
723: indices[i][0] = N - 1 - i;
724: indices[i][1] = 0;
725: indices[i][2] = val[i];
726: }
727: } else {
728: indices[0][0] = 0; indices[0][1] = 0; indices[0][2] = -(1 + 1);
729: indices[1][0] = 0; indices[1][1] = 1; indices[1][2] = -(2 + 1);
730: indices[2][0] = 1; indices[2][1] = 0; indices[2][2] = 0;
731: indices[3][0] = 1; indices[3][1] = 2; indices[3][2] = -(2 + 1);
732: indices[4][0] = 2; indices[4][1] = 1; indices[4][2] = 0;
733: indices[5][0] = 2; indices[5][1] = 2; indices[5][2] = 1;
734: }
735: } else {
736: PetscInt *subset, *work;
738: PetscMalloc2(k, &subset, k, &work);
739: for (i = 0; i < Nk; i++) {
740: PetscInt j, l, m;
742: PetscDTEnumSubset(N, k, i, subset);
743: for (j = 0; j < k; j++) {
744: PetscInt idx;
745: PetscBool flip = (PetscBool) (j & 1);
747: for (l = 0, m = 0; l < k; l++) {
748: if (l != j) work[m++] = subset[l];
749: }
750: PetscDTSubsetIndex(N, k - 1, work, &idx);
751: indices[i * k + j][0] = idx;
752: indices[i * k + j][1] = i;
753: indices[i * k + j][2] = flip ? -(subset[j] + 1) : subset[j];
754: }
755: }
756: PetscFree2(subset, work);
757: }
758: return(0);
759: }
761: /*@
762: PetscDTAltVStar - Apply a power of the Hodge star operator, which maps k-forms to (N-k) forms, to a k-form
764: Input Arguments:
765: + N - the dimension of the vector space, N >= 0
766: . k - the degree k of the k-form w, 0 <= k <= N
767: . pow - the number of times to apply the Hodge star operator: pow < 0 indicates that the inverse of the Hodge star operator should be applied |pow| times.
768: - w - a k-form, size [N choose k]
770: Output Arguments:
771: . starw = (star)^pow w. Each degree of freedom of a k-form is associated with a subset S of k coordinates of the N dimensional vector space: the Hodge start operator (star) maps that degree of freedom to the degree of freedom associated with S', the complement of S, with a sign change if the permutation of coordinates {S[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.
773: Level: intermediate
775: .seealso: PetscDTAltV, PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
776: @*/
777: PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw)
778: {
779: PetscInt Nk, i;
780: PetscErrorCode ierr;
783: if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree");
784: PetscDTBinomialInt(N, k, &Nk);
785: pow = pow % 4;
786: pow = (pow + 4) % 4; /* make non-negative */
787: /* pow is now 0, 1, 2, 3 */
788: if (N <= 3) {
789: if (pow & 1) {
790: PetscReal mult[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: PetscMalloc1(N, &subset);
803: if (pow % 2) {
804: PetscInt l = (pow == 1) ? k : N - k;
805: for (i = 0; i < Nk; i++) {
806: PetscBool sOdd;
807: PetscInt j, idx;
809: PetscDTEnumSplit(N, l, i, subset, &sOdd);
810: PetscDTSubsetIndex(N, l, subset, &idx);
811: PetscDTSubsetIndex(N, N-l, &subset[l], &j);
812: starw[j] = sOdd ? -w[idx] : w[idx];
813: }
814: } else {
815: for (i = 0; i < Nk; i++) starw[i] = w[i];
816: }
817: /* star^2 = -1^(k * (N - k)) */
818: if (pow > 1 && (k * (N - k)) % 2) {
819: for (i = 0; i < Nk; i++) starw[i] = -starw[i];
820: }
821: PetscFree(subset);
822: }
823: return(0);
824: }