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: }