Hodge * Operator

 Basics of wedge products
  In the description below, we are assuming that our underlying space is an orientable $C^\infty$ manifold.

We know vector space has a dual space that consists of functionals. Similarly, if we have a tangent plane there is a dual to it. Typically tangent space is spanned by ${ \frac{\partial }{\partial x^i} }$ where $i=1,2,\cdots,n$ for an n-dimensional basis for tangent space. Since this is a vector space, there is dual to it called ``cotangent space''. Just like vector space duals are functions to real numbers, the corresponding vector space dual is consists of differentials $dx_i$ and these are maps to real. These are called ``covectors''. 

If you consider Riemman integral $\int f dx$, here we consider $dx$ to be limit of $\delta x$ which is considered as small strip in $x-y$ plane. This concept can guide thinking on $dx_i$.
 

The differentials $dx_i$ have an operation called ``wedge'' operation, given by
$$\begin{equation} dx \wedge dy = - dy \wedge dx \end{equation}$$

In regular integration $dxdy$ is considered area element and $dxdydz$ volume element. Similarly, wedge product gives oriented area or more generally oriented volume.

Wedge operation maps to reals. Then $dx \wedge dx = - dx \wedge dx$ implies $dx \wedge dx=0$.
For a function $f(x,y,z)$ in $\mathcal{R}^3$, total differential leads to the following equation.
$$\begin{equation}   df = \frac{\partial{f}}{\partial{x}}dx + \frac{\partial{f}}{\partial{y}}dy + \frac{\partial{f}}{\partial{z}}dz \end{equation}$$
Here $dx,dy, dz$ are covectors or duals. If as in above each expression in addition, $df$ involves single differential or covector $dx,dy,dz$, these are called $1$-form. If you wedge two $1$-forms we get $2$($dx\wedge dy$) form etc.,$0$ forms are simply functions.

Collection of $k$ covectors is denoted by $\Lambda^k_n(M)$, where $M$ is $n$ dimensional manifold or space under consideration.

From the above equation it is evident that the $d$ operator acts on functions - that is $0-form$ and produces $1$ form. Using this as an example, we can write $d:\Lambda^k_n(M) \rightarrow \Lambda^{k+1}_{n+1}(M)$ which means that if we apply $d$ to a $k$ form, we get $k+1$ form.

From equation (2) we can also infer the dual space has basis consists of $dx,dy,dz$ or $3$-dimensional. Similarly for $n=3$ space we take $2$-forms, then the dimension is also $3$. In general, for  $\Lambda^k_n(M)$ the dimension of vector space is given by $^nC_k$. Since $^nC_k=^nC_{n-k}$, the dimensions of $k$ covectors and $n-k$ covector space is same.

Hodge *
Let $M$ be an $n$-dimensional $C^\infty$ Manifold. For a given integer $0 \leq k \leq n$, $\Lambda^k T^*_nM$ and $\Lambda^{n-k}T^{*}_nM$ have same dimension as vector space and they are isomorphic. Here $T^*M$ is dual to Tangent space $TM$. If $M$ has a Riemanian metric and oriented, for each point $p \in M$, there is a natural isomorphism
$$\begin{equation}   *:\Lambda^k T^*_nM \equiv \Lambda^{n-k}T^{*}_nM \end{equation}$$
By varying $p$ we get linear isomorphism
$$\begin{equation}   *:\mathcal{A}^k(M) \rightarrow \mathcal{A}^{n-k}(M) \end{equation}$$
This operator is called ``Hodge'' star operator where $\mathcal{A}^k$ represents vector space of $k$ forms etc.,

A bit of explanation for above definition. $C^\infty$ Manifold means a smooth manifold.
Riemannian metric is a positive definite inner product defined defined on the tangent space at each point of the Manifold $g_p: T_pM \times T_p M \rightarrow \mathcal{R}$ in such a way $g_p$ is $C^\infty$ at point $p$. Very standard way to express this is $ds^2 = \sum_{i,j=1}^n g_{ij}dx_i dx_j$.

Going back to Hodge star, if we write $V$ for $T_p M$, and use $V^*$ for $T_p^* M$, we get,
$$\begin{equation}   *: \Lambda^k V^* \rightarrow \Lambda^{n-k} V^* \end{equation}$$

This linear map can be defined by setting
$$\begin{equation}   *(\theta_1 \wedge \theta_2 \wedge \cdots \wedge \theta_k) = \theta_{k+1}\wedge \cdots \theta_n \end{equation}$$
In particular
$$\begin{equation}   *1 = \theta_1 \wedge \theta_2 \wedge \cdots \wedge \theta_n = 1 \end{equation}$$

Let $\omega$ be a one form for $X \in \Xi(M)$. For example, $\omega = f_x dx + f_y dy + f_z dz$ where $f_x = \frac{\partial f}{\partial x}$ etc., Then,
$$\begin{equation}   div X = *d*\omega \end{equation}$$

For example, for $\omega$ given in above example, we can work this out as follows:
$$\begin{align*}   *\omega = *(f_x dx + f_y dy + f_z dz) = f_x*dx+f_y*dy+f_z*dz \\   = f_x dy \wedge dz + f_y dz \wedge dx + f_z dx \wedge dy \\   d*\omega = f_{xx} dx \wedge dy \wedge dz + f_{yy}  dx \wedge dy \wedge dz + f_{zz}  dx \wedge dy \wedge dz \\   *d*\omega = f_{xx}+f_{yy}+f_{zz} = \frac{\partial^2f}{\partial x^2}+ \frac{\partial^2f}{\partial y^2}+ \frac{\partial^2f}{\partial z^2} = div X \end{align*}$$

Chain complex

 A cochain complex $\mathcal{C}$ is a collection of vector spaces ${C^k}_{k\in\mathcal{Z}}$ together with sequence of linear maps $d_k:C^k \rightarrow C^{k+1}$
$$\begin{equation}   \cdots \rightarrow C^{-1} \xrightarrow{d_{-1}} C^{0} \xrightarrow{d_{0}} C^1 \xrightarrow{d_{1}} C^2\xrightarrow{d_{2}}\cdots \end{equation}$$
with
$$\begin{equation}   d_k \circ d_{k-1} = 0 \end{equation}$$
${d_k}$ are collection of linear maps known as ``differentials'' of the cochain complex.
One relevant example of Cochain complex is the vector space $\Omega^{*}(M)$ of differential forms on Manifold together with exterior derivative.
$$\begin{equation}   \cdots \rightarrow \Omega^{-1}(M) \xrightarrow{d_{-1}} \Omega^{0}(M) \xrightarrow{d_{0}} \Omega^1(M) \xrightarrow{d_{1}} \Omega^2(M)\xrightarrow{d_{2}}\cdots,\;d\circ d = 0 \end{equation}$$
Above cochain complex is known as deRahm complex.

Mixed formulation of Laplacian

The following deals with mixed problem of Dirichelet kind for Poisson equation. To get started, we cast Poisson equation into first order equation.

Recall, original Poisson equation is
$$\begin{equation}   -u^{''}(x) = f(x) \end{equation}$$
Let's set $\sigma = -\text{grad}\; u$. $grad$ - aka Gradient is a vector field. In its full glory written as
$$\begin{equation}   \nabla u = (\frac{\partial u}{\partial x},\frac{\partial u}{\partial y}, \frac{\partial u}{\partial z} ) \end{equation}$$
clearly $\text{grad}\; u$ is a vector field. For our case, we simply set this to $\sigma= -\text{grad}\; u$.
Our original equation is $-u^{''}(x) = f(x)$. To express this in terms of $\sigma$, we use another operator - div.
Let's write down what div is.
The divergence of a vector field, say $F=(F_x,F_y,F_z)$, is a scalar function that represents the net rate of outward flux per unit volume at each point in the field. It gives a measure of how much the vector field is spreading out or compressing at a given point. The divergence is calculated as follows:
$$\begin{equation}   \text{div} F = \frac{\partial F_x}{\partial x}+ \frac{\partial F_x}{\partial y} + \frac{\partial F_x}{\partial z} \end{equation}$$
Clear from above definition if we $\text{div}\; \text{grad}\; u$, we get
$$\begin{equation}   \text{ div } \text{grad}\; u = \frac{\partial^2 u}{\partial x^2}+ \frac{\partial^2 u}{\partial y^2}+ \frac{\partial^2 u}{\partial z^2} \end{equation}$$
which in our original equation is same as $u^{''}$.
Hence, the original equation becomes,
$$\begin{equation}   \text{div}\; \text{grad}\; u(x) = f(x) \end{equation}$$
This in turn is shortened as,

$$\begin{equation}   \sigma = - \text{grad}\; u, \text{div}\; \sigma = f \end{equation}$$
The pair $(\sigma,u)$ can be characterized as critical point (unique) of the functional
$$\begin{equation}   I(\sigma,u) = \int_\Omega (\frac{1}{2} \sigma.\sigma - u \text{div}  \sigma)dx + \int_\Omega fu dx \end{equation}$$
over $H(div:\Omega) \times L^2(\Omega)$ where $H(div:\Omega) = {\sigma \in L^2; div \sigma \in L^2}$.
Equivalently one can solve weak problem
$$\begin{equation}   \int_{\Omega}\tau.\sigma dx - \int_{\Omega}u \text{div} \tau dx = 0, \tau \in H(div:\Omega) \end{equation}$$
$$\begin{equation}   \int_{\Omega} \text{div} \tau v dx = \int_{\Omega} f v dx \end{equation}$$
This fits into abstract framework if we define, $V=H(div:\Omega) \times L^(\Omega)$.
$$\begin{equation}   B(\sigma,u;\tau,v) = \int_{\Omega}\sigma.\tau dx - \int_{\Omega}u \text{div} \tau dx + \int_{\Omega}\text{div} \sigma v dx, F(\tau,v) = \int_{\Omega}fvdx \end{equation}$$


In this case the bilinear form $B$ is not coercive, and so the choice of subspaces an the analysis is not so simple
as for the standard finite element method for Poisson’s equation. Finite element discretizations based on such saddle point variational principles are called mixed finite element methods. Thus a mixed finite element for Poisson’s equation is obtained by choosing subspaces $\Sigma_h \subset H(div;\Omega)$ and $V_h \subset L_2(\Omega)$ and seeking a critical point of $I$ over $\Sigma_h \times V_h$. The resulting Galerkin method has the form: Find $\sigma_h \in \Sigma_h,u_h \in V_h$ satisfying
$$\begin{equation}       \int_{\Omega} \sigma_h \cdot \tau \, dx - \int_{\Omega} u_h \, \text{div} \tau \, dx = 0, \quad \forall \tau \in \Sigma_h,     \int_{\Omega} \text{div} \sigma_h v \, dx = \int_{\Omega} fv \, dx, \quad \forall v \in V_h.   \end{equation}$$
Since the bilinear form is not coercive, it is not automatic that the linear system is nonsingular.
If $f=0$, then $\int_{\Omega} fv \, dx=0$. This means $\int_{\Omega} \text{div} \sigma_h v \, dx = 0$. This foces $\sigma_h=0$ as $v$ is trail function. This, in turn, via $\int_{\Omega} \sigma_h \cdot \tau \, dx - \int_{\Omega} u_h \, \text{div} \tau \, dx = 0$ forces $\int_{\Omega} u_h \, \text{div} \tau \, dx = 0$. Then $u_h=0$.
Choosing $\tau = \sigma h$ and $v = uh$ and adding the discretized variational equations, it follows
immediately that when $f = 0, \sigma_h = 0$. However, $u_h$ need not vanish.