Welcome to spl’s documentation!

PSYDAC is a Python/Fortran 2003 library for B-Splines/NURBS and Computer Aided Design Algorithms.

PSYDAC can be used in three different ways:

  1. Fortran 90/95 subroutines through the file fortran/src/bsplines/bspline.F90
  2. Fortran 2003 objects, mainly the mapping and cad objects
  3. The same objects as in 2. but through Python

First Steps with PSYDAC

This document is meant to give a tutorial-like overview of PSYDAC.

The green arrows designate “more info” links leading to advanced sections about the described task.

By reading this tutorial, you’ll be able to:

  • compile a simple PSYDAC file
  • get familiar with parallel programing paradigms
  • create, modify and build a PSYDAC project.

Install PSYDAC

Examples

In this section, we describe some features of PSYDAC on simple examples.

See script.

Dive into PSYDAC

Contents

Introduction

Input and Output

B-Splines and NURBS

We start this section by recalling some basic properies about B-splines curves and surfaces. We also recall some fundamental algorithms (knot insertion and degree elevation).

For a basic introduction to the subject, we refer to the books [LP95] and [Far02].

A B-Splines family, (N_i)_{ 1 \leqslant i \leqslant n} of order k, can be generated using a non-decreasing sequence of knots T=(t_i)_{1\leqslant i \leqslant n + k}.

B-Splines series

The j-th B-Spline of order k is defined by the recurrence relation:

N_j^k = w_j^k N_j^{k-1} + ( 1 - w_{j+1}^k ) N_{j+1}^{k-1}

where,

w_j^k (x) = \frac{x-t_j}{t_{j+k-1}-t_{j}} \hspace{2cm} N_j^1(x) = \chi_{ \left[ t_j, t_{j+1} \right[ }(x)

for k \geq 1 and 1 \leq j \leq n.

We note some important properties of a B-splines basis:

  • B-splines are piecewise polynomial of degree p=k-1,
  • Compact support; the support of N_j^k is contained in \left[ t_j, t_{j+k} \right] ,
  • If x \in~ ] t_j,t_{j+1} [, then only the B-splines \{ N_{j-k+1}^k,\cdots,N_{j}^k \} are non vanishing at x,
  • Positivity: \forall j \in \{1,\cdots,n \}~~N_j(x) >0, ~~\forall x \in ] t_j, t_{j+k} [,
  • Partition of unity \sum_{i=1}^n N_i^{k}(x) = 1, \forall x \in \mathbb{R},
  • Local linear independence,
  • If a knot t_i has a multiplicity m_i then the B-spline is \mathcal{C}^{(p-m_i)} at t_i.

Knots vector families

There are two kind of knots vectors, called clamped and unclamped. Both families contains uniform and non-uniform sequences.

The following are examples of such knots vectors

  1. Clamped knots (open knots vector)
  • uniform

T_1 &= \{0, 0, 0, 1, 2, 3, 4, 5, 5, 5 \} \\ T_2 &= \{-0.2, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 0.8 \}

_images/bsplines_t1_p2.png _images/bsplines_t2_p2.png
  • non-uniform

T_3 &= \{0, 0, 0, 1, 3, 4, 5, 5, 5 \} \\ T_4 &= \{-0.2, -0.2, 0.4, 0.6, 0.8, 0.8 \}

_images/bsplines_t3_p2.png _images/bsplines_t4_p2.png
  1. Unclamped knots
  • uniform

T_5 &= \{0, 1, 2, 3, 4, 5, 6, 7 \} \\ T_6 &= \{-0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 \}

_images/bsplines_t5_p2.png _images/bsplines_t6_p2.png
  • non-uniform

T_7 &= \{0, 0, 3, 4, 7, 8, 9 \} \\ T_8 &= \{-0.2, 0.2, 0.4, 0.6, 1.0, 2.0, 2.5 \}

_images/bsplines_t7_p2.png _images/bsplines_t8_p2.png

B-Spline curve

The B-spline curve in \mathbb{R}^d associated to knots vector T=(t_i)_{1\leqslant i \leqslant n + k} and the control polygon (\mathbf{P}_i)_{ 1 \leqslant i \leqslant n} is defined by :

\mathcal{C}(t) = \sum_{i=1}^n N_i^k(t) \textbf{P}_i

In (Fig. ref{figBSplineCurve}), we give an example of a quadratic B-Spline curve, and its corresponding knot vector and control points.

_images/courbe_bsplines.png _images/basis_fct_p2_N5.png

We have the following properties for a B-spline curve:

  • If n=k, then \mathcal{C} is just a B’ezier-curve,
  • \mathcal{C} is a piecewise polynomial curve,
  • The curve interpolates its extremas if the associated multiplicity of the first and the last knot are maximum (i.e. equal to k), i.e. open knot vector,
  • Invariance with respect to affine transformations,
  • Strong convex-hull property:

if t_i \leq t \leq t_{i+1}, then \mathcal{C}(t) is inside the convex-hull associated to the control points \mathbf{P}_{i-p},\cdots,\mathbf{P}_{i},

  • Local modification : moving the i^{th} control point \mathbf{P}_{i} affects \mathcal{C}(t), only in the interval [t_i,t_{i+k}],
  • The control polygon approaches the behavior of the curve.

Note

In order to model a singular curve, we can use multiple control points : \mathbf{P}_{i}=\mathbf{P}_{i+1}.

Multivariate tensor product splines

Let us consider d knot vectors \mathcal{T} = \{T^1,T^2,\cdots,T^d\}. For simplicity, we consider that these knot vectors are open, which means that k knots on each side are duplicated so that the spline is interpolating on the boundary, and of bounds 0 and 1. In the sequel we will use the notation I=[0,1]. Each knot vector T^i, will generate a basis for a Schoenberg space, \mathcal{S}_{k_{i}}(T^i,I). The tensor product of all these spaces is also a Schoenberg space, namely \mathcal{S}_{\mathbf{k}}(\mathcal{T}), where \mathbf{k}=\{k_1,\cdots,k_d\}. The cube \mathcal{P}=I^d=[0,1]^d, will be referred to as a patch.

The basis for \mathcal{S}_{\mathbf{k}}(\mathcal{T}) is defined by a tensor product :

N_{\mathbf{i}}^{\mathbf{k}} := N_{i_1}^{k_1} \otimes N_{i_2}^{k_2} \otimes \cdots \otimes N_{i_d}^{k_d}

where, \mathbf{i}=\{i_1,\cdots , i_d \}.

A typical cell from \mathcal{P} is a cube of the form : Q_{\mathbf{i}}=[\xi_{i_1}, \xi_{i_1+1}] \otimes \cdots \otimes [\xi_{i_d}, \xi_{i_d+1}].

Deriving a B-spline curve

The derivative of a B-spline curve is obtained as:

\mathcal{C}^{\prime}(t) = \sum_{i=1}^{n} {N_{i}^{k}}^{\prime}(t) \mathbf{P}_i = \sum_{i=1}^{n} \left(\frac{p}{t_{i+p}-t_{i}}N_{i}^{k-1}(t) \mathbf{P}_i - \frac{p}{t_{i+1+p}-t_{i+1}}N_{i+1}^{k-1}(t) \mathbf{P}_i \right) = \sum_{i=1}^{n-1} {N_{i}^{k-1}}^{\ast}(t) \mathbf{Q}_i

where \mathbf{Q}_i = p \frac{\mathbf{P}_{i+1} - \mathbf{P}_i}{t_{i+1+p}-t_{i+1}}, and \{{N_{i}^{k-1}}^{\ast},~~1 \leq i \leq n-1\} are generated using the knot vector T^{\ast}, which is obtained from T by reducing by one the multiplicity of the first and the last knot (in the case of open knot vector), i.e. by removing the first and the last knot.

More generally, by introducing the B-splines family \{ {N_{i}^{k-j}}^{\ast}, 1 \leq i \leq n-j \} generated by the knots vector T^{j^{\ast}} obtained from T by removing the first and the last knot j times, we have the following result:

proposition

The j^{th} derivative of the curve \mathcal{C} is given by

\mathcal{C}^{(j)}(t) = \sum_{i=1}^{n-j} {N_{i}^{k-j}}^{\ast}(t) \mathbf{P}_i^{(j)}`

where, for j>0

\mathbf{P}_i^{(j)} = \frac{p-j+1}{t_{i+p+1}-t_{i+j}} \left( \mathbf{P}_{i+1}^{(j-1)} - \mathbf{P}_i^{(j-1)} \right) \\ \mbox{and} ~ ~ ~ \mathbf{P}_i^{(0)} = \mathbf{P}_i.

By denoting \mathcal{C}^{\prime} and \mathcal{C}^{\prime\prime} the first and second derivative of the B-spline curve \mathcal{C}, it is easy to show that:

We have,

  • \mathcal{C}^{\prime}(0) = \frac{p}{t_{p+2}} \left(\mathbf{P}_{2} - \mathbf{P}_1\right),
  • \mathcal{C}^{\prime}(1) = \frac{p}{1-t_{n}} \left(\mathbf{P}_{n} - \mathbf{P}_{n-1}\right),
  • \mathcal{C}^{\prime\prime}(0) = \frac{p(p-1)}{t_{p+2}} \left( \frac{1}{t_{p+2}}\mathbf{P}_{1} - \{ \frac{1}{t_{p+2}} + \frac{1}{t_{p+3}} \} \mathbf{P}_2 + \frac{1}{t_{p+3}}\mathbf{P}_{3} \right),
  • \mathcal{C}^{\prime\prime}(1) = \frac{p(p-1)}{1-t_{n}} \left( \frac{1}{1-t_{n}}\mathbf{P}_{n} - \{ \frac{1}{1-t_{n}} + \frac{1}{1-t_{n-1}} \} \mathbf{P}_{n-1} + \frac{1}{1-t_{n-1}}\mathbf{P}_{n-2} \right).

Example

Let us consider the quadratic B-spline curve associated to the knots vector T=\{000~\frac{2}{5}~\frac{3}{5}~111 \} and the control points \{ P_i, 1 \leq i \leq 5 \}:

\mathcal{C}(t) = \sum_{i=1}^{5} {N_{i}^{3}}^{\prime}(t) \mathbf{P}_i

we have,

\mathcal{C}^{\prime}(t) = \sum_{i=1}^{4} {N_{i}^{2}}^{\ast}(t) \mathbf{Q}_i

where

\mathbf{Q}_1 = 5 \{\mathbf{P}_{2} - \mathbf{P}_1\}, ~~~~\mathbf{Q}_2 = \frac{10}{3} \{ \mathbf{P}_{3} - \mathbf{P}_2\}, \\ \mathbf{Q}_3 = \frac{10}{3} \{ \mathbf{P}_{4} - \mathbf{P}_3\},~~~~\mathbf{Q}_4 = 5 \{\mathbf{P}_{5} - \mathbf{P}_4\}.

The B-splines \{ {N_{i}^{2}}^{\ast},~~1 \leq i \leq 4\} are associated to the knot vector T^{\ast}=\{00~\frac{2}{5}~\frac{3}{5}~11 \}.

Fundamental geometric operations

By inserting new knots into the knot vector, we add new control points without changing the shape of the B-Spline curve. This can be done using the DeBoor algorithm [dB01]. We can also elevate the degree of the B-Spline family and keep unchanged the curve [HHM05]. In (Fig. ref{refinement_curve_B_Spline}), we apply these algorithms on a quadratic B-Spline curve and we show the position of the new control points.

Knot insertion

After modification, we denote by \widetilde{n}, \widetilde{k}, \widetilde{T} the new parameters. (\textbf{Q}_i) are the new control points.

One can insert a new knot t, where t_j \leqslant t < t_{j+1}. For this purpose we use the DeBoor algorithm [dB01]:

\widetilde{n} = n+1 \\ \widetilde{k} = k \\ \widetilde{T} = \{ t_1,.., t_j, t, t_{j+1},.., t_{n+k}\} \\ \alpha_i = \left\{\begin{array}{cc}1 & 1 \leqslant i \leqslant j-k+1 \\\frac{t-t_i}{t_{i+k-1}-t_i} & j-k+2 \leqslant i \leqslant j \\0 & j+1 \leqslant i \end{array}\right. \\ \textbf{Q}_i = \alpha_i \textbf{P}_i + (1-\alpha_i) \textbf{P}_{i-1}

Many other algorithms exist, like blossoming for fast insertion algorithm. For more details about this topic, we refer to [NT93].

Order elevation

We can elevate the order of the basis, without changing the curve. Several algorithms exist for this purpose. We used the one by Huang et al. [PP91], [HHM05].

A quadratic B-spline curve and its control points. The knot vector is T = \{ 000, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1 1 1 \}.

_images/curve.png

The curve after a h-refinement by inserting the knots \{ 0.15, 0.35\} while the degree is kept equal to 2.

_images/curve_p0_n9.png

The curve after a p-refinement, the degree was raised by 1 (using cubic B-splines).

_images/curve_p2_n0.png

The curve after duplicating the multiplicity of the internal knots \{ \frac{1}{4}, \frac{1}{2}, \frac{3}{4} \}, this leads to a B’ezier description. We can then, split the curve into 4 pieces (sub-domains), each one will corresponds to a quadratic B’ezier curve.

_images/curve_p0_n3_bezier.png

Translation

Rotation

Todo

not yet available

Scaling

Todo

not yet available

References

[dB01](1, 2) C. de Boor. A Practical Guide to Splines. Applied Mathematical Sciences. Springer New York, 2001. ISBN 9780387953663. URL: https://books.google.de/books?id=m0QDJvBI_ecC.
[Far02]G. Farin. Curves and surfaces for CAGD: a practical guide. Morgan Kaufmann Pub. Inc., San Francisco, CA, USA, 2002. ISBN 1-55860-737-4.
[HHM05](1, 2) Qi-Xing Huang, Shi-Min Hu, and Ralph R. Martin. Fast degree elevation and knot insertion for b-spline curves. Computer Aided Geometric Design, 22(2):183 – 197, 2005. URL: http://www.sciencedirect.com/science/article/B6TYN-4DXBTHR-2/2/d5b3eec2f4c230c8051623c1c000beae, doi:DOI: 10.1016/j.cagd.2004.11.001.
[LP95]W. Tiller L. Piegl. The NURBS Book. Springer-Verlag, Berlin, Heidelberg, 1995. second ed.
[NT93]Goldman R. N. and Lyche T. Knot Insertion and Deletion Algorithms for B-Spline Curves and Surfaces. SIAM, Philadelphia, USA, 1993. ISBN 9780898713060.
[PP91]Hartmut Prautzsch and Bruce Piper. A fast algorithm to raise the degree of spline curves. Comput. Aided Geom. Des., 8:253–265, October 1991. URL: http://portal.acm.org/citation.cfm?id=124930.124932, doi:10.1016/0167-8396(91)90015-4.

GLT

Where do the GLTs come from?

The main aim of this paragraph is to present a crucial example that highlights the importance of the GLT algebra when dealing with linear systems coming from the discretization of PDEs. Let us start with some preliminaries. In detail, we will recall the notion of symbol of a matrix-sequence and the basic idea behind the GLT theory.

Spectral preliminaries

The following one is a rather informal definition of symbol of a matrix-sequence.

example:

When d_n=n, d=1, D=[0,\pi], \{A_n\}_n\sim_\lambda f means

References

Exterior Algebra

Let V be a real vector space of dimension n.

Definition, Alternating algebraic forms:

For each k, we define \Alt^k V as the space of alternating k-linear maps V \times \cdots \times V \rightarrow \mathbb{R}.

Note

  • \Alt^0 = \mathbb{R},
  • \Alt^1 = V^{*} is the dual space of V (the space of covectors)

Definition, Exterior product:

For \omega \in \Alt^j and \eta \in \Alt^k, their exterior (wedge) product is given by:

(\omega \wedge \eta ) (v_1, \cdots, v_{j+k}) = \sum_{\sigma} (\mathrm{sign}~ \sigma) \omega (v_{\sigma(1)}, \cdots, v_{\sigma(j)}) \eta (v_{\sigma(j+1)}, \cdots, v_{\sigma(j+k)})

for all v_i \in V. Where the sum is over all permutations \sigma of \{ 1,\cdots,j+k \}, for which \sigma(1)< \cdots <\sigma(j) and \sigma(j+1)< \cdots <\sigma(j+k).

Note

  • The exterior product is bilinear, associative,
  • anti-commutative: \eta \wedge \omega = (-1)^{jk} \omega \wedge \eta for all \omega \in \Alt^j and \eta \in \Alt^k.

Definition, Grassmann Algebra:

Grassmann Algebra is defined by:

\Alt V := \bigoplus_k \Alt^k V

This is a anti-commutative graded algebra. Also called Exterior Algebra of V^{*}

In the case of V=\mathbb{R}^n, we have:

  • \Alt V^0 \sim \mathbb{R},
  • \Alt V^1 \sim \mathbb{R}^n,
  • \Alt V^{n-1} \sim \mathbb{R}^n, using Riesz representation theorem,
  • \Alt V^n \sim \mathbb{R}, using the map v \longmapsto \det(v,v_1,\cdots,v_{n-1}).

Basis

Let v_1,\cdots,v_n be a basis of V and \mu_1,\cdots,\mu_n the associated dual basis for V^* (\mu_i(v_j) = \delta_{ij}).

For any increasing permutations \sigma, \rho : \{ 1,\cdots,k \} \longrightarrow \{ 1,\cdots,n \}, we have:

\mu_{\sigma(1)} \wedge \cdots \wedge \mu_{\sigma(k)} (v_{\rho(1)}, \cdots, v_{\rho(k)}) = \chi_{\sigma,\rho}

thus the \binom {n}{k} algebraic k-forms \mu_{\sigma(1)} \wedge \cdots \wedge \mu_{\sigma(k)}, form a basis for \Alt^k V and \dim \Alt^k V = \binom {n}{k}.

Definition, Interior product:

Let \omega be a k-form, and v \in V. The interior product of \omega and v is the (k-1)-form \omega \lrcorner v defined by:

\omega \lrcorner v (v_1,\cdots,v_{k-1}) = \omega (v,v_1,\cdots,v_{k-1})

  • We have for \omega \in \Alt^k V, \eta \in \Alt^l V and v \in V:

(\omega \wedge \eta) \lrcorner v = (\omega \lrcorner v)\wedge \eta + (-1)^k \omega \wedge (\eta \lrcorner v)

Definition, Inner product:

If V is has an inner product, then \Alt^k V is endowed with an inner product given by:

(\omega , \eta) = \sum_{\rho} \omega (v_{\rho(1)}, \cdots, v_{\rho(k)}) \eta (v_{\rho(1)}, \cdots, v_{\rho(k)}), ~~~\forall \omega, \eta \in \Alt^k V.

where the sum is over increasing sequences \rho : \{ 1,\cdots,k \} \longrightarrow \{ 1,\cdots,n \}, and v_1, \cdots,v_n is any orthonormal basis.

Orientation and Volume form

Todo

add Orientation and Volume form

Definition, Pullback:

A linear transformation of vector spaces L: V \rightarrow W induces a transformation L^{*}: \Alt W \rightarrow \Alt V, called the pullback, and given by:

L^{*} \omega (v_{1}, \cdots, v_{k}) = \omega (L v_{1}, \cdots, L v_{k}),~~~~~ \forall \omega \in \Alt^k W,~~~ v_{1}, \cdots, v_{k} \in V

  • The pullback acts contravariantly: if U \xrightarrow{~K~} V \xrightarrow{~L~} W then,

    \Alt W \xrightarrow{~K^{*}~} \Alt V \xrightarrow{~L^{*}~} \Alt U

  • L^{*} (\omega \wedge \eta) = L^{*} \omega \wedge L^{*} \eta

Let V be a subspace of W. For the inclusion \imath_V : V \longrightarrow W, we can define its pullback \imath_V^{*}: this is a surjection of \Alt W onto \Alt V.

If W has an inner product and \pi_V : W \longrightarrow V is the orthogonal projection. We can define its pullback \pi_V^{*} : this an injection of \Alt V onto \Alt W.

Let us consider the composition : W shortstack{\pi_V \ \longrightarrow} V shortstack{\imath_V \ \longrightarrow} W, and its pullback \pi_V^* \imath_V^*.

Definition, The tangential and normal parts:

  • \pi_V^* \imath_V^* associates for each \omega \in \Alt^k its tangential part \omega_{\parallel} with respect to V :

(\pi_V^* \imath_V^* \omega) (v_1,\cdots,v_k) = \omega (\pi_V v_1, \cdots, \pi_V v_k), ~~~~~\forall v_1,\cdots,v_k \in W.

  • \omega - \pi_V^* \imath_V^* \omega associates for each \omega \in \Alt^k its normal part \omega_{\perp} with respect to V.

The tangential part of \omega vanishes if and only if the image of \omega in \Alt^k V vanishes.

Let V be an oriented inner product space, with volume form \mbox{vol}. Let \omega \in \Alt^k V. We can define a new linear map L_{\omega} as the composition of \Alt^{n-k} V \longrightarrow \Alt^n V such as:

\mu \longmapsto \omega \wedge \mu

and the canonical isomorphism of \Alt^n V onto \mathbb{R}, and using the Riesz representation theorem, there exists an element \star \omega \in \Alt^{n-k} V such that : L_{\omega} (\mu) = (\star \omega , \mu), i.e.:

\omega \wedge \mu = (\star \omega , \mu) \mbox{vol}, ~~~\omega \in \Alt^{k}, ~\mu \in \Alt^{n-k}

Definition, The Hodge star operation:

The linear map which maps \Alt^k V onto \Alt^{n-k} V \omega \longmapsto \star \omega is called the Hodge star operator.

  • If e_1,\cdots,e_n is any positively oriented orthonormal basis, and \sigma a permutation, we have

\omega(e_{\sigma(1)}, \cdots, e_{\sigma(k)}) = (\mathrm{sign} \sigma) \star \omega(e_{\sigma(k+1)}, \cdots, e_{\sigma(n)})

  • \star \star \omega = (-1)^{k(n-k)} \omega, ~~~\forall \omega \in \Alt^k V, thus the Hodge star is an isometry.
  • (\star \omega)_{\parallel} = \star (\omega_{\perp}) and (\star \omega)_{\perp} = \star (\omega_{\parallel})
  • the image of \star \omega in \Alt^k V vanishes if and only if \omega_{\perp} vanishes.

\begin{tabular}{|c|l|} \hline $\Alt^0 \mathbb{R}^3 \cong \mathbb{R}$ & $c \leftrightarrow c$ \\ % \hline $\Alt^1 \mathbb{R}^3 \cong \mathbb{R}^3$ & $u_1 \diff x_1 + u_2 \diff x_2 + u_3 \diff x_3 \leftrightarrow u$ \\ % \hline $\Alt^2 \mathbb{R}^3 \cong \mathbb{R}^3$ & $u_3 \diff x_1 \wedge \diff x_2 - u_2 \diff x_1 \wedge \diff x_3 + u_1 \diff x_2 \wedge \diff x_3 + \leftrightarrow u$ \\ % \hline $\Alt^3 \mathbb{R}^3 \cong \mathbb{R}$ & $c \diff x_1 \wedge \diff x_2 \wedge \diff x_3 \leftrightarrow c$ \\ \hline \end{tabular}

\begin{tabular}{|c|l|} \hline $ \wedge : \Alt^1 \mathbb{R}^3 \times \Alt^1 \mathbb{R}^3 \longrightarrow \Alt^2 \mathbb{R}^3$ & $\times : \mathbb{R}^3 \times \mathbb{R}^3 \longrightarrow \mathbb{R}^3$ \\ $ \wedge : \Alt^1 \mathbb{R}^3 \times \Alt^2 \mathbb{R}^3 \longrightarrow \Alt^3 \mathbb{R}^3$ & $\cdot : \mathbb{R}^3 \times \mathbb{R}^3 \longrightarrow \mathbb{R}$ \\ \hline \end{tabular}

\begin{tabular}{|c|l|} \hline $ L^* : \Alt^0 \mathbb{R}^3 \longrightarrow \Alt^0 \mathbb{R}^3 $ & $\id : \mathbb{R} \longrightarrow \mathbb{R}$ \\ $ L^* : \Alt^1 \mathbb{R}^3 \longrightarrow \Alt^1 \mathbb{R}^3 $ & $L^T : \mathbb{R}^3 \longrightarrow \mathbb{R}^3$ \\ $ L^* : \Alt^2 \mathbb{R}^3 \longrightarrow \Alt^2 \mathbb{R}^3 $ & $(\det L )L^{-1} : \mathbb{R}^3 \longrightarrow \mathbb{R}^3$ \\ $ L^* : \Alt^3 \mathbb{R}^3 \longrightarrow \Alt^3 \mathbb{R}^3 $ & $(\det L) : \mathbb{R} \longrightarrow \mathbb{R}$ ~~~($c \longmapsto c \det L$) \\ \hline \end{tabular}

\begin{tabular}{|c|l|} \hline $ \lrcorner v : \Alt^1 \mathbb{R}^3 \longrightarrow \Alt^0 \mathbb{R}^3 $ & $v \cdot : \mathbb{R}^3 \longrightarrow \mathbb{R}$ \\ $ \lrcorner v : \Alt^2 \mathbb{R}^3 \longrightarrow \Alt^1 \mathbb{R}^3 $ & $v \times : \mathbb{R}^3 \longrightarrow \mathbb{R}^3$ \\ $ \lrcorner v : \Alt^3 \mathbb{R}^3 \longrightarrow \Alt^2 \mathbb{R}^3 $ & $v : \mathbb{R} \longrightarrow \mathbb{R}^3$ ~~~($c \longmapsto c v$) \\ \hline \end{tabular}

\begin{tabular}{|c|l|} \hline inner product on $\Alt^k \mathbb{R}^3$ induced & dot product on $\mathbb{R}$ and $\mathbb{R}^3$ \\ by dot product on $\mathbb{R}^3$ & \\ $\volume = \diff x_1 \wedge \diff x_2 \wedge \diff x_3$ & $(v_1,v_2,v_3) \longmapsto \det(v_1|v_2|v_3)$ \\ \hline \end{tabular}

\begin{tabular}{|c|l|} \hline $ \star : \Alt^0 \mathbb{R}^3 \longrightarrow \Alt^3 \mathbb{R}^3 $ & $\id : \mathbb{R} \longrightarrow \mathbb{R}$ \\ $ \star : \Alt^1 \mathbb{R}^3 \longrightarrow \Alt^2 \mathbb{R}^3 $ & $\id : \mathbb{R}^3 \longrightarrow \mathbb{R}^3$ \\ \hline \end{tabular}

Exterior Calculus on manifolds and Differential forms

Let \Omega be a smooth manifold, of dimension n.

  • \forall x \in \Omega we denote by T_x \Omega the tangent space. This is a vector space of dimension n,
  • tangent bundle \{ (x,v), ~~ x \in \Omega, v \in T_x \Omega \},
  • Applying the exterior algebra to the tangent spaces, we obtain the exterior forms bundle, whose elements are pairs (x,\mu) with x \in \Omega and \mu \in \Alt^k T_x \Omega.
  • a differential k-form \omega is a section of this bundle. This is a map which associates to each x \in \Omega an element \omega_x \in \Alt^k T_x \Omega,
  • if the map \mathcal{L}_{\omega}^k : x \longmapsto \omega_x (v_1(x), \cdots, v_k(x)) is smooth (whenever v_i are smooth), we say that \omega is a smooth differential k-form,
  • we define \Lambda^k(\Omega) the space of all smooth k-forms on \Omega,
  • \Lambda^0(\Omega) = \mathcal{C}^{\infty}(\Omega),
  • if the map \mathcal{L}_{\omega}^k is \mathcal{C}^{m}(\Omega), we define differential k-forms with less smoothness \mathcal{C}^{m} \Lambda^k (\Omega).

Let \Omega be a smooth manifold, of dimension n.

Exterior product:

if \omega \in \Lambda^k(\Omega) and \eta \in \Lambda^j(\Omega), we may define \omega \wedge \eta as (\omega \wedge \eta)_x = \omega_x \wedge \eta_x and the Grassmann algebra \Lambda(\Omega) := \bigoplus_k \Lambda^k(\Omega)

Differential forms can be differentiated and integrated, without recourse to any additional structure, such as a metric or a measure.

Exterior differentiation:

For each \omega \in \Lambda^k(\Omega), can define the (k+1)-form \diff \omega \in \Lambda^{k+1}(\Omega), such as:

\diff\omega_x(v_1,\cdots,v_{k+1}) = \sum_{j=1}^{k+1} (-1)^j \partial_{v_j} \omega_x(v_1,\cdots,\hat{v_j},\cdots,v_{k+1})

where the hat is used to indicated a suppressed argument.

This defines a graded linear operator of degree +1, of \Lambda(\Omega) onto \Lambda(\Omega).

We have the following properties:

  • \diff \circ \diff = 0
  • \diff (\omega \wedge \eta) = \diff \omega \wedge \eta + (-1)^k \omega \wedge \diff \eta, ~~\forall \omega \in \Lambda^k(\Omega), \eta \in \Lambda^j(\Omega),
  • (Pullback) let \phi be a smooth map of \Omega onto \Omega^{\prime}. Then \phi^*(\omega \wedge \eta) = \phi^*(\omega) \wedge \phi^*(\eta) and \phi^* (\diff \omega) = \diff (\phi^* \omega),
  • (Interior product) the interior product of a differential k-form \omega with a vector field v,
  • we obtain a (k-1)-form by : (\omega \lrcorner v)_x := \omega_x \lrcorner v_x,
  • (Trace operator) the pullback i_{\partial \Omega}^* of i_{\partial \Omega} is the trace operator \trace

Integration:

  • If f is an oriented, piecewise smooth k-dimensional submanifold of \Omega, and \omega is a continuous k-form, then th integral \int_f \omega is well defined :
    • [0-forms] can be evaluated at points,
    • [1-forms] can be integrated over directed curves,
    • [2-forms] can be integrated over directed surfaces,
  • (Inner product) The L^2-inner product of two differential k-forms on an oriented Riemannian manifold \Omega is defined as :

(\omega,\eta)_{L^2 \Lambda^k} = \int_{\Omega} (\omega_x,\eta_x) \volume = \int \omega \wedge \star \eta

The completion of \Lambda^k(\Omega) in the corresponding norm defines the Hilbert space L^2 \Lambda^k(\Omega).

We have the following results:

  • (Integration) if \phi is an orientation-preserving diffeomorphism, then

\int_{\Omega} \phi^* \omega = \int_{\Omega^{\prime}} \omega, ~~~ \forall \omega \in \Lambda^n(\Omega^{\prime})

Theorem, Stokes theorem:

If \Omega is an oriented n-manifold with boundary \partial \Omega, then

\int_{\Omega} \diff \omega = \int_{\partial \Omega} \trace \omega, ~~~ \forall \omega \in \Lambda^{n-1}(\Omega)

Theorem, Integration by parts:

If \Omega is an oriented n-manifold with boundary \partial \Omega, then

\int_{\Omega} \diff \omega \wedge \eta = (-1)^{k-1} \int_{\Omega} \omega \wedge \diff \eta + \int_{\partial \Omega} \trace \omega \wedge \trace \eta, ~~~ \forall \omega \in \Lambda^{k}(\Omega), \eta \in \Lambda^{n-k-1}(\Omega)

Sobolev spaces of differential forms

As for the classical case, we can define the Sobolev spaces as:

  • H^s \Lambda^k(\Omega) is the space of differential k-forms such that \mathcal{L}_{\omega}^k \in H^s(\Omega).
  • H \Lambda^k(\Omega) = \{ \omega \in L^2 \Lambda^k(\Omega),~~ \diff \omega \in L^2 \Lambda^{k+1}(\Omega) \}. The associated norm is :

\| \omega \|_{H \Lambda^k}^2 = \| \omega \|_{H \Lambda}^2 := \| \omega \|_{L^2 \Lambda^k}^2 + \| \diff \omega \|_{L^2 \Lambda^{k+1}}^2

  • H \Lambda^{0}(\Omega) coincides with H^1 \Lambda^{0}(\Omega),
  • H \Lambda^{n}(\Omega) coincides with L^2 \Lambda^{n}(\Omega),
  • for 0 < k < n, we have H^1 \Lambda^k(\Omega) \subset H \Lambda^k(\Omega) \subset L^2 \Lambda^k(\Omega), strictly.

\begin{tabular}{|c|c c c c c|} \hline $k$ & $\Lambda^k$ & $H \Lambda^k$ & $\diff \omega$ & $\int_f \omega$ & $\kappa \omega$ \\ \hline & & & & & \\ 0 & $\mathcal{C}^{\infty}$ & $H^1$ & $\nabla \omega$ & $\omega(f)$ & $0$ \\ 1 & $\mathcal{C}^{\infty}(\mathbb{R}^3)$ & $H(\rots,\mathbb{R}^3)$ & $\rots \omega$ & $\int_f \omega \cdot t \diff \mathcal{H}_1$ & $x \longmapsto x \cdot \omega(x)$ \\ 2 & $\mathcal{C}^{\infty}(\mathbb{R}^3)$ & $H(\divs, \mathbb{R}^3)$ & $\divs \omega$ & $\int_f \omega \cdot n \diff \mathcal{H}_2$ & $x \longmapsto x \times \omega(x)$ \\ 3 & $\mathcal{C}^{\infty}$ & $L^2$ & $0$ & $\int_f \omega \diff \mathcal{H}_3$ & $x \longmapsto x \omega(x)$ \\ & & & & & \\ \hline \end{tabular}

Cohomology and De Rham Complex

The De Rham complex is the sequence of spaces and mappings

0 \xrightarrow{\quad} \Lambda^0(\Omega) \xrightarrow{~\diff~} \Lambda^1(\Omega) \xrightarrow{~\diff~} \cdots \xrightarrow{~\diff~} \Lambda^n(\Omega) \xrightarrow{\quad} 0

Since, \diff \circ \diff = 0, we have

\mathcal{R}(\diff : \Lambda^{k-1}(\Omega) \longrightarrow \Lambda^k(\Omega)) \subset \mathcal{N}(\diff : \Lambda^{k}(\Omega) \longrightarrow \Lambda^{k+1}(\Omega))

If \Omega is an oriented Riemannian manifold, we have the following cohomology:

0 \xrightarrow{\quad} H \Lambda^0(\Omega) \xrightarrow{~\diff~} H \Lambda^1(\Omega) \xrightarrow{~\diff~} \cdots \xrightarrow{~\diff~} H \Lambda^n(\Omega) \xrightarrow{\quad} 0

The coderivative operator \delta : \Lambda^{k}(\Omega) \longrightarrow \Lambda^{k-1}(\Omega) is defined as:

\star \delta \omega = (-1)^k \diff \star \omega,~~~ \omega \in \Lambda^k(\Omega)

  • we have

(\diff \omega , \eta ) = (\omega , \delta \eta ) + \int_{\partial \Omega} \trace \omega \wedge \trace \eta, ~~~ \forall \omega \in \Lambda^{k}(\Omega), \eta \in \Lambda^{k+1}(\Omega),

  • \delta is a graded linear operator of degree -1.
  • \delta is the formal adjoint of \diff whenever \omega or \eta vanishes near the boundary.
  • we define the spaces

H^* \Lambda^k(\Omega) = \{ \omega \in L^2 \Lambda^k(\Omega),~~ \delta \omega \in L^2 \Lambda^{k-1}(\Omega) \}.

we have H^* \Lambda^k(\Omega) = \star H \Lambda^{n-k}(\Omega).

  • we obtain the dual complex

0 \xleftarrow{\quad} H^* \Lambda^0(\Omega) \xleftarrow{~\delta~} H^* \Lambda^1(\Omega) \xleftarrow{~\delta~} \cdots \xleftarrow{~\delta~} H^* \Lambda^n(\Omega) \xleftarrow{\quad} 0

Cohomology with boundary conditions

Let \Lambda_0^k(\Omega) be the subspace of \Lambda^k(\Omega) of smooth k-forms with compact support. We have \diff \Lambda_0^k \subset \Lambda_0^{k+1}.

The De Rham complex with the compact support is

0 \xrightarrow{\quad} \Lambda^0_0(\Omega) \xrightarrow{~\diff~} \Lambda^1_0(\Omega) \xrightarrow{~\diff~} \cdots \xrightarrow{~\diff~} \Lambda^n_0(\Omega) \xrightarrow{\quad} 0

Recall that the closure of \Lambda_0^k(\Omega) in H \Lambda^k(\Omega) is

H_0 \Lambda^k(\Omega) = \{ \omega \in H \Lambda^k(\Omega),~~ \trace \omega =0\}.

The L^2 version of the last complex is

0 \xrightarrow{\quad} H_0 \Lambda^0(\Omega) \xrightarrow{~\diff~} H_0 \Lambda^1(\Omega) \xrightarrow{~\diff~} \cdots \xrightarrow{~\diff~} H_0 \Lambda^n(\Omega) \xrightarrow{\quad} 0

Definition, Harmonic forms:

The harmonic k-forms are the differential k-forms that verify the differential equations

\left\{ \begin{aligned} \diff \omega &=& 0,\\ \delta \omega &=& 0,\\ \trace \star \omega &=& 0.\\ \end{aligned} \right.

this defines the following space,

\mathfrak{H}^k (\Omega) = \{ \omega \in H \Lambda^k(\Omega) \cap H_0^* \Lambda^k(\Omega),~~\diff \omega = 0, \delta \omega = 0 \}

We can also define the following space,

\mathfrak{H}_0^k (\Omega) = \{ \omega \in H_0 \Lambda^k(\Omega) \cap H^* \Lambda^k(\Omega),~~\diff \omega = 0, \delta \omega = 0 \}

As we can see, \star \mathfrak{H}^k (\Omega) = \mathfrak{H}_0^{n-k} (\Omega).

Proposition, Poincaré duality:

There is an isomorphism between the k th De Rham cohomology space and the (n-k) th cohomology space with boundary conditions.

Homological Algebra and Hilbert complexes

Homological Algebra

  • A cochain complex is a sequence of vector spaces and linear maps
  • k-cocycles \mathfrak{Z}^k := \mathcal{N}(d_k),
  • k-coboundaries \mathfrak{B}^k := \mathcal{R}(d_{k-1}),
  • k-cohomology \mathcal{H}^k(V) := \mathfrak{Z}^k / \mathfrak{B}^k,
  • we say that the sequence is exact, if the cohomology vanishes (i.e. \forall~k,~~ \mathcal{H}^k(V) = \{0\}),
  • Given two cochain complexes V,V^{\prime}, a cochain map f =(f_k) (such as \diff^{\prime}_k f_k = f_{k+1} \diff_k)

\begin{array}{ccccccccc} \cdots & \longrightarrow & V_{k-1} & \mbox{\shortstack{$\diff_{k-1}$ \\ $\longrightarrow$}} & V_{k} & \mbox{\shortstack{$\diff_k$ \\ $\longrightarrow$}} & V_{k+1} & \longrightarrow~\cdots \\ & & \downarrow f_{k-1} & & \downarrow f_{k} & & \downarrow f_{k+1} & & \\ \cdots & \longrightarrow & V_{k-1}^{\prime} & \mbox{\shortstack{$\diff_{k-1}^{\prime}$ \\ $\longrightarrow$}} & V_{k}^{\prime} & \mbox{\shortstack{$\diff_k^{\prime}$ \\ $\longrightarrow$}} & V_{k+1}^{\prime} & \longrightarrow~\cdots \end{array}

  • f_k maps k-cochains to k-cochains and k-coboundaries to k-coboundaries, thus induces a map \mathcal{H}^k(f) : \mathcal{H}^k(V) \longrightarrow \mathcal{H}^k(V^{\prime}).

Let V^{\prime} \subset V be two cochain complexes,

  • The inclusion \imath_V is a cochain map and thus induces a map of cohomology \mathcal{H}^k(V^{\prime}) \longrightarrow \mathcal{H}^k(V),

  • If there exists a cochain projection of V onto V^{\prime}, (this leads to \pi \circ \imath = \id_{V^{\prime}}) so \mathcal{H}^k(\pi) \circ \mathcal{H}^k(\imath) = \id_{\mathcal{H}^k(V^{\prime})}.

    \begin{array}{ccccccc} \cdots & \longrightarrow & V_{k-1} & \mbox{\shortstack{$\diff_{k-1}$ \\ $\longrightarrow$}} & V_{k} & \longrightarrow~\cdots \\ & & \pi_{k-1} \downarrow \uparrow \imath & & \pi_{k} \downarrow \uparrow \imath & & \\ \cdots & \longrightarrow & V_{k-1}^{\prime} & \mbox{\shortstack{$\diff_{k-1}$ \\ $\longrightarrow$}} & V_{k}^{\prime} & \longrightarrow~\cdots \end{array}

Thus, \mathcal{H}^k(\imath) is injective and \mathcal{H}^k(\pi) is surjective. Hence, if one of the cohomology spaces \mathcal{H}^k(V) vanishes, then so does \mathcal{H}^k(V^{\prime})

Cycles and boundaries of the De Rham complex

  • k-cocycles

\mathfrak{Z}^k = \{ \omega \in H\Lambda^k(\Omega),~~ \diff \omega = 0 \}, ~~~ \mathfrak{Z}^{*k} = \{ \omega \in H^*\Lambda^k(\Omega),~~ \delta \omega = 0 \},

\mathfrak{Z}_0^k = \{ \omega \in H_0\Lambda^k(\Omega),~~ \diff \omega = 0 \}, ~~~ \mathfrak{Z}_0^{*k} = \{ \omega \in H_0^*\Lambda^k(\Omega),~~ \delta \omega = 0 \},

  • k-coboundaries

\mathfrak{B}^k = \diff H\Lambda^{k-1}(\Omega), ~~~ \mathfrak{B}^{* k} = \delta \Lambda^{k+1}(\Omega),

\mathfrak{B}_0^k = \diff H_0\Lambda^{k-1}(\Omega), ~~~ \mathfrak{B}_0^{* k} = \delta \Lambda_0^{k+1}(\Omega),

  • each of the spaces of cycles is closed in \mathcal{H} \Lambda^k(\Omega) (\mathcal{H}^* \Lambda^k(\Omega)), as well in L^2 \Lambda^k(\Omega).
  • each of the spaces of boundaries is closed in L^2 \Lambda^k(\Omega).
  • let \perp denotes the orthogonal complement in L^2 \Lambda^k(\Omega),

\mathfrak{Z}^{k \perp} \subset \mathfrak{B}^{k \perp} = \mathfrak{Z}_0^{* k} , ~~~ \mathfrak{Z}^{* k \perp} \subset \mathfrak{B}^{* k \perp} = \mathfrak{Z}_0^{k}

\mathfrak{Z}_0^{k \perp} \subset \mathfrak{B}_0^{k \perp} = \mathfrak{Z}^{* k} , ~~~ \mathfrak{Z}_0^{* k \perp} \subset \mathfrak{B}_0^{* k \perp} = \mathfrak{Z}^{k}

The Hodge decomposition

There are two Hodge decompositions, with different boundary conditions,

  1. L^2 \Lambda^k(\Omega) = \underbrace{\mathfrak{B}^{k}}_{\mathfrak{Z}_0^{* k\perp}} \oplus \underbrace{\mathfrak{H}^{k} \oplus \mathfrak{B}_0^{* k}}_{\mathfrak{Z}_0^{* k}=\mathfrak{B}^{k\perp}} = \overbrace{\mathfrak{B}^{k} \oplus \mathfrak{H}^{k}}^{\mathfrak{Z}^{k}=\mathfrak{B}_0^{* k\perp}} \oplus \overbrace{\mathfrak{B}_0^{* k}}^{\mathfrak{Z}^{k\perp}}

  2. L^2 \Lambda^k(\Omega) = \underbrace{\mathfrak{B}_0^{k}}_{\mathfrak{Z}^{* k\perp}} \oplus \underbrace{\mathfrak{H}_0^{k} \oplus \mathfrak{B}^{* k}}_{\mathfrak{Z}^{* k}=\mathfrak{B}_0^{k\perp}} = \overbrace{\mathfrak{B}_0^{k} \oplus \mathfrak{H}_0^{k}}^{\mathfrak{Z}_0^{k}=\mathfrak{B}^{* k\perp}} \oplus \overbrace{\mathfrak{B}^{* k}}^{\mathfrak{Z}_0^{k\perp}}

Summary

\begin{tabular}{|c||c|c|c|c|} \hline $\omega^k \in \Lambda^k(\Omega)$ & $k=0$ & $k=1$ & $k=2$ & $k=3$ \\ \hline $\diff \omega^k$ & $\Grad u$ & $\Curl \uu$ & $\Div \uu$ & $-$ \\ $\delta \omega^k$ & $-$ & $-\Div \uu$ & $\Curl \uu$ & $-\Grad u$ \\ $\mathfrak{i}_{\boldsymbol{\beta}} \omega^k$ & $-$ & $\boldsymbol{\beta} \cdot \uu$ & $\uu \times \boldsymbol{\beta}$ & $u \boldsymbol{\beta}$ \\ $\mathfrak{j}_{\boldsymbol{\beta}} \omega^k$ & $u \boldsymbol{\beta}$ & $-\uu \times \boldsymbol{\beta}$ & $\boldsymbol{\beta} \cdot \uu$ & $-$ \\ $L_{\boldsymbol{\beta}} \omega^k$ & $\boldsymbol{\beta} \cdot \Grad u$ & $\Grad \left(\boldsymbol{\beta} \cdot \uu \right) + \left(\Curl \uu \right) \times \boldsymbol{\beta}$ & $\Curl \left(\uu \times \boldsymbol{\beta} \right) + \boldsymbol{\beta} \Div \uu$ & $\Div \left( u \boldsymbol{\beta} \right)$ \\ $\mathcal{L}_{\boldsymbol{\beta}} \omega^k$ & $-\Div \left( u \boldsymbol{\beta} \right)$ & $-\Curl \left(\uu \times \boldsymbol{\beta} \right) - \boldsymbol{\beta} \Div \uu$ & $-\Grad \left(\boldsymbol{\beta} \cdot \uu \right) - \left(\Curl \uu \right) \times \boldsymbol{\beta}$ & $-\boldsymbol{\beta} \cdot \Grad u$ \\ \hline $\tr \omega^k$ & $u(\xx)$ & $\uu(\xx) \times \nn(\xx)$ & $\uu(\xx) \cdot \nn(\xx)$ & $-$ \\ \hline \hline $H \Lambda^k(\Omega)$ & $\Hgrad$ & $\Hcurl$ & $\Hdiv$ & $\Ltwo$ \\ $V_k$ & $\Vgrad$ & $\Vcurl$ & $\Vdiv$ & $\Vltwo$ \\ \hline \end{tabular}

References

DeRham sequences

here without boundary conditions

\mathbb{R} \hookrightarrow \Hgrad \xrightarrow{\quad \Grad \quad} \Hcurl \xrightarrow{\quad \Curl \quad} \Hdiv \xrightarrow{\quad \Div \quad} \Ltwo \xrightarrow{\quad} 0

Pullbacks

In the case where the physical domain \Omega := \mathcal{F}(\hat{\Omega}) is the image of a logical domain \hat{\Omega} by a smooth mapping \mathcal{F} (at least \mathcal{C}^1), we have the following parallel diagrams

\begin{array}{ccccccc} \Hgrad & \xrightarrow{\quad \Grad \quad} & \Hcurl & \xrightarrow{\quad \Curl \quad} & \Hdiv & \xrightarrow{\quad \Div \quad} & \Ltwo \\ \igrad \Bigg\uparrow & & \icurl \Bigg\uparrow & & \idiv \Bigg\uparrow & & \iltwo \Bigg\uparrow \\ \HgradLogical & \xrightarrow{\quad \Grad \quad} & \HcurlLogical & \xrightarrow{\quad \Curl \quad} & \HdivLogical & \xrightarrow{\quad \Div \quad} & \LtwoLogical \\ % \end{array}

Where the mappings \igrad, \icurl, \idiv and \iltwo are called pullbacks and are given by

\phi (x) :=& \igrad \hat{\phi} (\hat{x}) = \hat{\phi}(\mathcal{F}^{-1}(x)) \\ \Psi (x) :=& \icurl \hat{\Psi} (\hat{x}) = \left( D \mathcal{F} \right)^{-T} \hat{\Psi}(\mathcal{F}^{-1}(x)) \\ \Phi (x) :=& \idiv \hat{\Phi} (\hat{x}) = \frac{1}{J} D \mathcal{F} \hat{\Phi}(\mathcal{F}^{-1}(x)) \\ \rho (x) :=& \iltwo \hat{\rho} (\hat{x}) = \hat{\rho}(\mathcal{F}^{-1}(x))

where D \mathcal{F} is the jacobian matrix of the mapping \mathcal{F}.

Note

The pullbacks \igrad, \icurl, \idiv and \iltwo are isomorphisms between the corresponding spaces.

Discrete Spaces

Let us suppose that we have a sequence of finite subspaces for each of the spaces involved in the DeRham sequence. The discrete DeRham sequence stands for the following commutative diagram between continuous and discrete spaces

\begin{array}{ccccccc} \Hgrad & \xrightarrow{\quad \Grad \quad} & \Hcurl & \xrightarrow{\quad \Curl \quad} & \Hdiv & \xrightarrow{\quad \Div \quad} & \Ltwo \\ \Pigrad \Bigg\downarrow & & \Picurl \Bigg\downarrow & & \Pidiv \Bigg\downarrow & & \Piltwo \Bigg\downarrow \\ \Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Curl \quad} & \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo \\ % \end{array}

When using a Finite Elements methods, we often deal with a reference element, and thus we need also to apply the pullbacks on the discrete spaces. In fact, we have again the following parallel diagram

\begin{array}{ccccccc} \Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Curl \quad} & \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo \\ \igrad \Bigg\uparrow & & \icurl \Bigg\uparrow & & \idiv \Bigg\uparrow & & \iltwo \Bigg\uparrow \\ \VgradLogical & \xrightarrow{\quad \Grad \quad} & \VcurlLogical & \xrightarrow{\quad \Curl \quad} & \VdivLogical & \xrightarrow{\quad \Div \quad} & \VltwoLogical \\ % \end{array}

Since, the pullbacks are isomorphisms in the previous diagram, we can define a one-to-one correspondance

\phi :=& \igrad \hat{\phi}, \quad \phi \in \Vgrad, \hat{\phi} \in \VgradLogical \\ \Psi :=& \icurl \hat{\Psi}, \quad \Psi \in \Vcurl, \hat{\Psi} \in \VcurlLogical \\ \Phi :=& \idiv \hat{\Phi}, \quad \Phi \in \Vdiv, \hat{\Phi} \in \VdivLogical \\ \rho :=& \iltwo \hat{\rho}, \quad \rho \in \Vltwo, \hat{\rho} \in \VltwoLogical

We have then, the following results

\Grad \phi =& \icurl \Grad \hat{\phi} , \quad \phi \in \Vgrad \\ \Curl \Psi =& \idiv \Curl \hat{\Psi} , \quad \Psi \in \Vcurl \\ \Div \Phi =& \iltwo \Div \hat{\Phi} , \quad \Phi \in \Vdiv

Projectors

In some cases, one may need to define projectors on smooth functions

\begin{array}{ccccccc} \Cinfinity & \xrightarrow{\quad \Grad \quad} & \Cinfinity & \xrightarrow{\quad \Curl \quad} & \Cinfinity & \xrightarrow{\quad \Div \quad} & \Cinfinity \\ \Pigrad \Bigg\downarrow & & \Picurl \Bigg\downarrow & & \Pidiv \Bigg\downarrow & & \Piltwo \Bigg\downarrow \\ \Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Curl \quad} & \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo \\ \end{array}

Discrete DeRham sequence for B-Splines

Buffa et al [BSV09] show the construction of a discrete DeRham sequence using B-Splines, (here without boundary conditions)

\begin{array}{ccccccc} \Hgrad & \xrightarrow{\quad \Grad \quad} & \Hcurl & \xrightarrow{\quad \Curl \quad} & \Hdiv & \xrightarrow{\quad \Div \quad} & \Ltwo \\ \Pigrad \Bigg\downarrow & & \Picurl \Bigg\downarrow & & \Pidiv \Bigg\downarrow & & \Piltwo \Bigg\downarrow \\ \Vgradspline & \xrightarrow{\quad \Grad \quad}& \Vcurlspline & \xrightarrow{\quad \Curl \quad} & \Vdivspline& \xrightarrow{\quad \Div \quad} & \Vltwospline \\ \end{array}

1d case

  1. DeRham sequence is reduced to

\mathbb{R} \hookrightarrow \underbrace{\mathcal{S}^{p}}_{\VgradLogical} \xrightarrow{\quad \Grad \quad} \underbrace{\mathcal{S}^{p-1}}_{\VltwoLogical} \xrightarrow{\quad} 0

  1. The recursion formula for derivative writes

{N_i^p}'(t)=D_i^{p}(t)-D_{i+1}^{p}(t) \quad \mbox{where} \quad D_{i}^{p}(t) = \frac{p}{t_{i+p+1}-t_i}N_i^{p-1}(t)

  1. we have \mathcal{S}^{p-1} = \mathbf{span}\{ N_i^{p-1}, 1 \leq i \leq n-1 \} = \mathbf{span}\{ D_i^p, 1 \leq i \leq n-1 \} which is a change of basis as a diagonal matrix
  2. Now if u \in S^p, with and expansion u = \sum_i u_i N_i^p, we have

\begin{align*} u^{\prime} = \sum_i u_i \left( N_i^p \right)^{\prime} = \sum_i (-u_{i-1} + u_i) D_i^p % \label{} \end{align*}

  1. If we introduce the B-Splines coefficients vector \mathbf{u} := \left( u_i \right)_{1 \leq i \leq n} (and \mathbf{u}^{\star} for the derivative), we have

\mathbf{u}^{\star} = D \mathbf{u}

where D is the incidence matrix (of entries -1 and +1)

Discrete derivatives:

\mathcal{G} = D

2d case

In 2d, the are two De-Rham complexes:

\begin{array}{ccccc} \Hgrad & \xrightarrow{\quad \Grad \quad} & \Hcurl & \xrightarrow{\quad \Rots \quad} & \Ltwo \\ \Pigrad \Bigg\downarrow & & \Picurl \Bigg\downarrow & & \Piltwo \Bigg\downarrow \\ \Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Rots \quad} & \Vltwo \\ \end{array}

and

\begin{array}{ccccc} \Hgrad & \xrightarrow{\quad \Curl \quad} & \Hdiv & \xrightarrow{\quad \Div \quad} & \Ltwo \\ \Pigrad \Bigg\downarrow & & \Pidiv \Bigg\downarrow & & \Piltwo \Bigg\downarrow \\ \Vgrad & \xrightarrow{\quad \Grad \quad} & \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo \\ \end{array}

Let I be the identity matrix, we have

Discrete derivatives:

\mathcal{G} = \begin{pmatrix} D \otimes I \\ I \otimes D \end{pmatrix}

\mathcal{C} = \begin{pmatrix} I \otimes D \\ - D \otimes I \end{pmatrix} \quad \mbox{[scalar curl],} \quad \mathcal{C} = \begin{pmatrix} - I \otimes D & D \otimes I \end{pmatrix} \quad \mbox{[vectorial curl]}

\mathcal{D} = \begin{pmatrix} D \otimes I & I \otimes D \end{pmatrix}

3d case

Discrete derivatives:

\mathcal{G} = \begin{pmatrix} D \otimes I \otimes I \\ I \otimes D \otimes I \\ I \otimes I \otimes D \end{pmatrix}

\mathcal{C} = \begin{pmatrix} 0 & - I \otimes I \otimes D & I \otimes D \otimes I \\ I \otimes I \otimes D & 0 & - D \otimes I \otimes I \\ - I \otimes D \otimes I & D \otimes I \otimes I & 0 \end{pmatrix}

\mathcal{D} = \begin{pmatrix} D \otimes I \otimes I & I \otimes D \otimes I & I \otimes I \otimes D \end{pmatrix}

Note

From now on, we will denote the discrete derivative by \mathbb{D}_k for the one going from V_k to V_{k+1}.

Algebraic identities

Let us consider the discretization of the exterior derivative

\omega^{k+1} = \diff \omega^k

multiplying by a test function \eta^{k+1} and integrating over the whole computation domain, we get

\left( \eta^{k+1}, \omega^{k+1} \right)_{k+1} = \left( \eta^{k+1}, \diff \omega^{k} \right)_{k+1}

let E^{k+1}, W^{k} and W^{k+1} be the vector representation of \eta^{k+1}, \omega^{k} and \omega^{k+1}. We get

{E^{k+1}}^T M_{k+1} W^{k+1} = {E^{k+1}}^T D_{k+1,k} W^{k}

where

D_{k+1,k} = \left( \left( \eta^{k+1}_i, \diff \omega^{k}_j \right)_{k+1} \right)_{i,j}

On the other hand, using the coderivative, we get

\left( \eta^{k+1}, \omega^{k+1} \right)_{k+1} = \left( \delta \eta^{k+1}, \omega^{k} \right)_{k} + BC

Let us now introduce the following matrix

D_{k,k+1} = \left( \left( \delta \eta^{k+1}_i, \omega^{k}_j \right)_{k} \right)_{i,j}

hence,

{E^{k+1}}^T D_{k,k+1} W^{k} = \left( \mathbb{D}^{\star}_{k+1} E^{k+1} \right)^T M_{k} W^{k}

Therefor, we have the following important result

Proposition:

  • D_{k+1,k} = D_{k,k+1} + BC
  • D_{k+1,k} = M_{k+1} \mathbb{D}^T_k
  • D_{k,k+1} = {\mathbb{D}^{\star}_{k+1}}^T M_{k}

References

[BSV09]A. Buffa, G. Sangalli, and R. Vazquez. Isogeometric analysis in electromagnetics: b-splines approximation. Comput. Methods Appl. Mech. Engrg, 199:1143–1152, 2009.

API

You will find here both the Fortran doxygen documentation as well as the Python-API.

Fortran API

Doxygen

Python API

psydac

psydac package
Subpackages
psydac.api package
Subpackages
psydac.api.ast package
Subpackages
psydac.api.ast.tests package
Submodules
psydac.api.ast.tests.boundary module
psydac.api.ast.tests.poisson module
psydac.api.ast.tests.system_1 module
psydac.api.ast.tests.system_2 module
psydac.api.ast.tests.system_3 module
psydac.api.ast.tests.test_nodes module
Module contents
Submodules
psydac.api.ast.basic module
psydac.api.ast.evaluation module
psydac.api.ast.expr module
psydac.api.ast.fem module
psydac.api.ast.glt module
psydac.api.ast.linalg module
psydac.api.ast.linalg_kernels module
psydac.api.ast.nodes module
psydac.api.ast.parser module
psydac.api.ast.utilities module
Module contents
psydac.api.printing package
Submodules
psydac.api.printing.pycode module
Module contents
psydac.api.tests package
Submodules
psydac.api.tests.build_domain module
psydac.api.tests.test_2d_biharmonic module
psydac.api.tests.test_2d_laplace module
psydac.api.tests.test_2d_mapping_biharmonic module
psydac.api.tests.test_2d_mapping_laplace module
psydac.api.tests.test_2d_mapping_poisson module
psydac.api.tests.test_2d_multipatch_mapping_maxwell module
psydac.api.tests.test_2d_multipatch_mapping_poisson module
psydac.api.tests.test_2d_multipatch_poisson module
psydac.api.tests.test_2d_navier_stokes module
psydac.api.tests.test_2d_poisson module
psydac.api.tests.test_api_1d_compatible_spaces module
psydac.api.tests.test_api_2d_compatible_spaces module
psydac.api.tests.test_api_2d_fields module
psydac.api.tests.test_api_2d_scalar_analytical_mapping module
psydac.api.tests.test_api_2d_system module
psydac.api.tests.test_api_2d_vector module
psydac.api.tests.test_api_2d_vector_mapping module
psydac.api.tests.test_api_2d_vector_multipatch_mapping module
psydac.api.tests.test_api_3d_scalar module
psydac.api.tests.test_api_3d_scalar_mapping module
psydac.api.tests.test_api_3d_vector module
psydac.api.tests.test_api_3d_vector_mapping module
psydac.api.tests.test_api_expr_2d_scalar module
psydac.api.tests.test_api_feec_1d module
psydac.api.tests.test_api_feec_2d module
psydac.api.tests.test_api_feec_3d module
psydac.api.tests.test_api_glt_2d_scalar module
psydac.api.tests.test_api_glt_2d_scalar_mapping module
psydac.api.tests.test_api_glt_2d_vector module
psydac.api.tests.test_assembly module
psydac.api.tests.test_equation module
psydac.api.tests.test_postprocessing module
psydac.api.tests.test_quadorder module
psydac.api.tests.utils module
Module contents
Submodules
psydac.api.basic module
psydac.api.discretization module
psydac.api.equation module
psydac.api.essential_bc module
psydac.api.expr module
psydac.api.feec module
psydac.api.fem module
psydac.api.glt module
psydac.api.grid module
psydac.api.postprocessing module
psydac.api.settings module
psydac.api.utilities module
Module contents
psydac.cad package
Subpackages
psydac.cad.tests package
Submodules
psydac.cad.tests.test_geometry module
Module contents
Submodules
psydac.cad.cad module
psydac.cad.geometry module
psydac.cad.multipatch module
psydac.cad.utils module
Module contents
psydac.cmd package
Submodules
psydac.cmd.mesh module
Module contents
psydac.core package
Subpackages
psydac.core.tests package
Submodules
psydac.core.tests.test_bsplines module
psydac.core.tests.test_bsplines_pyccel module
psydac.core.tests.test_kernels module
Module contents
Submodules
psydac.core.bsplines module
psydac.core.bsplines_pyccel module
psydac.core.kernels module
Module contents
psydac.ddm package
Subpackages
psydac.ddm.tests package
Submodules
psydac.ddm.tests.test_cart_1d module
psydac.ddm.tests.test_cart_2d module
psydac.ddm.tests.test_cart_3d module
psydac.ddm.tests.test_multicart_2d module
psydac.ddm.tests.test_partition module
Module contents
Submodules
psydac.ddm.basic module
psydac.ddm.blocking_data_exchanger module
psydac.ddm.cart module
psydac.ddm.interface_data_exchanger module
psydac.ddm.nonblocking_data_exchanger module
psydac.ddm.partition module
psydac.ddm.petsc module
psydac.ddm.utilities module
Module contents
psydac.feec package
Subpackages
psydac.feec.multipatch package
Subpackages
psydac.feec.multipatch.examples package
Submodules
psydac.feec.multipatch.examples.h1_source_pbms_conga_2d module
psydac.feec.multipatch.examples.hcurl_eigen_pbms_conga_2d module
psydac.feec.multipatch.examples.hcurl_source_pbms_conga_2d module
psydac.feec.multipatch.examples.mixed_source_pbms_conga_2d module
psydac.feec.multipatch.examples.ppc_test_cases module
Module contents
Submodules
psydac.feec.multipatch.api module
psydac.feec.multipatch.fem_linear_operators module
psydac.feec.multipatch.multipatch_domain_utilities module
psydac.feec.multipatch.operators module
psydac.feec.multipatch.plotting_utilities module
psydac.feec.multipatch.utilities module
psydac.feec.multipatch.utils_conga_2d module
Module contents
psydac.feec.tests package
Submodules
psydac.feec.tests.test_axis_projection module
psydac.feec.tests.test_commuting_projections module
psydac.feec.tests.test_differentiation_matrices module
psydac.feec.tests.test_global_projectors module
psydac.feec.tests.test_projections_parallel module
Module contents
Submodules
psydac.feec.derivatives module
psydac.feec.global_projectors module
psydac.feec.pull_push module
psydac.feec.pushforward module
Module contents
psydac.fem package
Subpackages
psydac.fem.tests package
Submodules
psydac.fem.tests.analytical_profiles_1d module
psydac.fem.tests.analytical_profiles_base module
psydac.fem.tests.splines_error_bounds module
psydac.fem.tests.test_eval_fields_parallel module
psydac.fem.tests.test_product module
psydac.fem.tests.test_spline_histopolation module
psydac.fem.tests.test_spline_interpolation module
psydac.fem.tests.test_splines module
psydac.fem.tests.test_splines_par module
psydac.fem.tests.utilities module
Module contents
Submodules
psydac.fem.basic module
psydac.fem.context module
psydac.fem.grid module
psydac.fem.partitioning module
psydac.fem.projectors module
psydac.fem.splines module
psydac.fem.tensor module
psydac.fem.vector module
Module contents
psydac.linalg package
Subpackages
psydac.linalg.tests package
Submodules
psydac.linalg.tests.test_block module
psydac.linalg.tests.test_fft module
psydac.linalg.tests.test_kron_direct_solver module
psydac.linalg.tests.test_kron_stencil_matrix module
psydac.linalg.tests.test_linalg module
psydac.linalg.tests.test_solvers module
psydac.linalg.tests.test_stencil_interface_matrix module
psydac.linalg.tests.test_stencil_matrix module
psydac.linalg.tests.test_stencil_vector module
psydac.linalg.tests.test_stencil_vector_space module
Module contents
Submodules
psydac.linalg.basic module
psydac.linalg.block module
psydac.linalg.direct_solvers module
psydac.linalg.fft module
psydac.linalg.kernels module
psydac.linalg.kron module
psydac.linalg.solvers module
psydac.linalg.stencil module
psydac.linalg.topetsc module
psydac.linalg.utilities module
Module contents
psydac.mapping package
Subpackages
psydac.mapping.tests package
Submodules
psydac.mapping.tests.test_discrete_mapping module
psydac.mapping.tests.visual_test_discrete_mapping_2d module
psydac.mapping.tests.visual_test_discrete_mapping_3d_surface module
Module contents
Submodules
psydac.mapping.discrete module
Module contents
psydac.polar package
Subpackages
psydac.polar.tests package
Submodules
psydac.polar.tests.test_c1_linops module
psydac.polar.tests.test_c1_projections module
Module contents
Submodules
psydac.polar.c1_cart module
psydac.polar.c1_linops module
psydac.polar.c1_projections module
psydac.polar.c1_spaces module
psydac.polar.dense module
Module contents
psydac.pyccel package
Subpackages
psydac.pyccel.ast package
Submodules
psydac.pyccel.ast.basic module
psydac.pyccel.ast.builtins module
psydac.pyccel.ast.core module
psydac.pyccel.ast.datatypes module
psydac.pyccel.ast.functionalexpr module
psydac.pyccel.ast.mathext module
psydac.pyccel.ast.numbers module
psydac.pyccel.ast.numpyext module
psydac.pyccel.ast.singleton module
psydac.pyccel.ast.type_inference module
psydac.pyccel.ast.utilities module
Module contents
psydac.pyccel.codegen package
Subpackages
psydac.pyccel.codegen.printing package
Submodules
psydac.pyccel.codegen.printing.pycode module
Module contents
Module contents
psydac.pyccel.errors package
Submodules
psydac.pyccel.errors.errors module
psydac.pyccel.errors.messages module
Module contents
psydac.pyccel.symbolic package
Submodules
psydac.pyccel.symbolic.lambdify module
Module contents
Module contents
psydac.utilities package
Submodules
psydac.utilities.quadratures module
psydac.utilities.utils module
psydac.utilities.vtk module
Module contents
Submodules
psydac.version module
Module contents

Indices and tables