Consider the BVP from linear elasticity:
\begin{align*}
-\mu \Delta\textbf{u} – (\lambda + \mu) \nabla(\nabla \cdot \textbf{u})) &= \textbf{f}, \text{ in } \Omega \subset\mathbb{R}^2, \\
\textbf{u} &= \textbf{g}_D, \text{ on } \partial \Omega.
\end{align*}
First, we multiply by a test function and integrate both sides,
\begin{align*}
\int_{\Omega}( -\mu \textbf{v} \Delta \textbf{u} – (\lambda + \mu) \textbf{v} \nabla (\nabla\cdot \textbf{u})))dX = \int_{\Omega} \textbf{v}\textbf{f} dX.
\end{align*}
Now I am not sure what the best way to proceed is. I think I should break everything up into components, but I am not confident. I think I will end up using the divergence theorem (twice?). Would anyone be able to help me work out the details? I think the $\Delta \textbf{u}$ term should be simple, but I have no clue what identities to use for $\nabla (\nabla \cdot \textbf{u})$. I appreciate any help.
[Math] Deriving the weak form for linear elasticity equation
calculus-of-variationsfinite element methodmultivariable-calculuspartial differential equations
Related Solutions
Here is a standard proof by contradiction. Assume that the inequality does not hold. Then for all $n$ there is $u_n\in H^1(\Omega)$ such that $$ \|\nabla u_n\|_{L^2(\Omega)}^2 + \|u_n\|_{L^2(\partial\Omega)}^2 < n \|u_n\|_{H^1(\Omega)}^2. $$ This implies $u_n\ne0$, and we can assume w.l.o.g. $\|u_n\|_{H^1(\Omega)}=1$. Then $\nabla u_n\to0$ in $L^2(\Omega)$ and $u\to0$ in $L^2(\partial\Omega)$ follows immediately.
Since $(u_n)$ is bounded in $H^1$, after extracting a subsequence (denoted the same for simplicity), we have $u_n\rightharpoonup u$ in $H^1(\Omega)$. Then $\nabla u=0$, and $u$ is a constant. Since $\|u\|_{L^2(\partial\Omega)}=0$, it follows $u=0$. By compact embedding, $u_n\to0$ in $L^2(\Omega)$. This implies $u\to0$ strongly in $H^1(\Omega)$, which is a contradiction to $\|u_n\|_{H^1}=1$.
The OP's equation is: $$ a(x,y)\frac{\partial^2 u}{\partial x^2}+b(x,y) \frac{\partial^2 u}{\partial y^2} +c(x,y)\frac{\partial^2 u}{\partial x \partial y}\\+d(x,y)\frac{\partial u}{\partial x} +e(x,y)\frac{\partial u}{\partial y}+f(x,y)u=g(x,y) $$ But, for some good reasons, we shall consider instead: $$ \begin{bmatrix} \partial / \partial x & \partial / \partial y \end{bmatrix} \begin{bmatrix} A(x,y) & C(x,y) \\ C(x,y) & B(x,y) \end{bmatrix} \begin{bmatrix} \partial u / \partial x \\ \partial u/ \partial y \end{bmatrix}\\ +D(x,y)\frac{\partial u}{\partial x}+E(x,y)\frac{\partial u}{\partial y}+F(x,y)u=G(x,y) $$ Seeking for similarities: $$ \frac{\partial}{\partial x}\left[A(x,y)\frac{\partial u}{\partial x}+C(x,y)\frac{\partial u}{\partial y}\right]+ \frac{\partial}{\partial y}\left[C(x,y)\frac{\partial u}{\partial x}+B(x,y)\frac{\partial u}{\partial y}\right]\\ +D(x,y)\frac{\partial u}{\partial x}+E(x,y)\frac{\partial u}{\partial y}+F(x,y)u=G(x,y) \quad \Longleftrightarrow \\ A\frac{\partial^2 u}{\partial x^2}+B\frac{\partial^2 u}{\partial y^2}+2C\frac{\partial^2 u}{\partial x \partial y}\\ +\left[\frac{\partial A}{\partial x}+\frac{\partial C}{\partial y}+D\right]\frac{\partial u}{\partial x} +\left[\frac{\partial C}{\partial x}+\frac{\partial B}{\partial y}+E\right]\frac{\partial u}{\partial y}+Fu=G $$ we conclude that the OP's equation can be rewritten as: $$ \begin{bmatrix} \partial / \partial x & \partial / \partial y \end{bmatrix} \begin{bmatrix} a & c/2 \\ c/2 & b \end{bmatrix} \begin{bmatrix} \partial u / \partial x \\ \partial u/ \partial y \end{bmatrix} +D\frac{\partial u}{\partial x}+E\frac{\partial u}{\partial y}+fu=g $$ provided that: $$ D = d(x,y)-\left[\frac{\partial a}{\partial x}+\frac{\partial c}{\partial y}\right] \quad ; \quad E = e(x,y)-\left[\frac{\partial c}{\partial x}+\frac{\partial b}{\partial y}\right] $$ With these modifications, the equation is suitable for numerical treatment. We only have to "downsize" the numerical method in the accompanying document from 3-D to 2-D.
The terms $\,+fu=g\,$ are easy, so we shall do them first.
It is argued in this answer @ MSE
that integration (points) at the vertices of a finite element are often the best. Accompanying pictures are inserted here too for convenience:
An interesting consequence is splitting of the quadrilateral in four linear triangles:
This makes the discretization of the terms $\,+fu=g\,$ extremely simple:
$$
+ \frac{1}{4} \begin{bmatrix} f_1\cdot\Delta_1 & 0 & 0 & 0 \\
0 & f_2\cdot\Delta_2 & 0 & 0 \\
0 & 0 & f_3\cdot\Delta_3 & 0 \\
0 & 0 & 0 & f_4\cdot\Delta_4 \end{bmatrix}
\begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} =
\frac{1}{4} \begin{bmatrix} g_1\cdot\Delta_1 \\ g_2\cdot\Delta_2 \\ g_3\cdot\Delta_3 \\ g_4\cdot\Delta_4 \end{bmatrix}
$$
Here $\Delta_k$ is twice the area of the triangle numbered as $(k)$.
The diffusion term has the form:
$$
\frac{\partial Q_x}{\partial x} + \frac{\partial Q_y}{\partial y}
$$
with
$$
\begin{bmatrix} Q_x \\ Q_y \end{bmatrix} =
\begin{bmatrix} a & c/2 \\ c/2 & b \end{bmatrix}
\begin{bmatrix} \partial u / \partial x \\ \partial u/ \partial y \end{bmatrix}
$$
If $\,u\,$ is interpreted as a temperature then this may be considered as a heat flow $\,\vec{Q}\,$
in a medium with anisotropic conductivity.
In this way, the diffusion term can be treated with the standard Galerkin method, exactly as in the abovementioned
reference, or according to an answer
at MSE
with very much the same content.
With help of the differentiation matrix, for each (of the four !) triangles in our basic quadrilateral, the $3 \times 3$ Finite Element matrix for diffusion alone
is like this, with $\Delta/2 = $ area of a triangle : $\Delta/4 \times$
$$ -
\begin{bmatrix} (y_2 - y_3) & -(x_2 - x_3) \\
(y_3 - y_1) & -(x_3 - x_1) \\
(y_1 - y_2) & -(x_1 - x_2) \end{bmatrix} / \Delta
\begin{bmatrix} a & c/2 \\ c/2 & b \end{bmatrix}
\begin{bmatrix} +(y_2 - y_3) & +(y_3 - y_1) & +(y_1 - y_2) \\
-(x_2 - x_3) & -(x_3 - x_1) & -(x_1 - x_2) \end{bmatrix} /\Delta
$$
And take notice of the minus sign!
The Finite Element
assembly scheme
is thus employed at an elementary level, using the topology:
1 2 3 2 4 1 3 1 4 4 3 2
Our Ansatz for the advection element matrix resembles the one for diffusion, but without the OP's $(a,b,c)$ tensor:
$$
M = - \frac{\Delta}{4} \times
\begin{bmatrix} +(y_2 - y_3) & -(x_2 - x_3) \\
+(y_3 - y_1) & -(x_3 - x_1) \\
+(y_1 - y_2) & -(x_1 - x_2) \end{bmatrix} / \Delta
\begin{bmatrix} +(y_2 - y_3) & +(y_3 - y_1) & +(y_1 - y_2) \\
-(x_2 - x_3) & -(x_3 - x_1) & -(x_1 - x_2) \end{bmatrix} /\Delta
$$
Now determine the values of $D(x,y)$ and $E(x,y)$ at the midpoints $(x,y)$ of each of the triangle edges
$(i,j) = (1,2) \to (2,3) \to (3,1)$ and form the inner products:
$$
P_{ij} = D(x,y)(x_j-x_i)+E(x,y)(y_j-y_i)
$$
Then multiply these contributions with the Ansatz, while using an upwind scheme , for $i \ne j$ :
$$
M_{ij} := M_{ij}\times\max(0,-P_{ij}) \quad ; \quad M_{ji} := M_{ji}\times\max(0,-P_{ji})
$$
The main diagonal terms must be made equal to minus the sum of the off-diagonal terms in order to finish the advection matrix.
Does the above sound a bit improbable? The secret behind it is in the section Voronoi Regions
of the 2-D reference.
There we find the following formula for the resistor ($R_3$) equivalent of diffusion:
$$
R_3 = \frac{ \mbox{"length" of } R_3 }{ \mbox{conductivity} \, \times \, \mbox{"diameter" of } R_3 }
$$
Turning it upside down and leaving out the conductivity - as has been done - we have for the Ansatz:
$$
\mbox{matrix entry} = \frac{\mbox{"diameter" of edge}}{\mbox{"length" of edge}}
$$
This is multiplied with the inner product of a "velocity" and an edge,
thus resulting in a "flux", being the projection of a "velocity" times the diameter ("area") of the edge.
At last all elementary parts must be assembled together, giving the completed finite element matrix
for the bulk of the problem.
Hopefully the OP can take care of the boundary conditions and take it from here.
Warning: anisotropy can make the latter exercise a bit tricky.
Best Answer
Let us start with some background in tensor calculus.
For arbitrary domain $D$, tensor $A_{ij}$ and vector $b_i$ the following Green's theorem holds: $$\int_D \mathrm{div} A \cdot b \,\mathrm{d}x = \int_{\partial D} (An)\cdot b \,\mathrm{d}s - \int_D A : \nabla v\,\mathrm{d}x,$$ where $n$ is the normal vector to the boundary $\partial D$ and $:$ denotes the double dot product, i.e. for tensors $C_{ij}$ and $D_{ij}$ we have $$C : D = \sum_{i,j=1}^3 C_{ij} D_{ij}.$$ The divergence of a tensor field $A_{ij}$ is defined as the vector with components $$(\mathrm{div} A)_i = \sum_{j=1}^3 \frac{\partial A_{ij}}{\partial x_j}$$ The gradient of a vector field $b_i$ is the tensor with components $$(\nabla b)_{ij} = \frac{\partial b_i}{\partial x_j}.$$
Back to linear elasticity. People usually write the governing equations in the following way: $$-\mathrm{div}\,\sigma(u) = f,$$ where the stress tensor is $$\sigma(u) = 2 \mu \,\epsilon(u) + \lambda \,\mathrm{tr}\,\epsilon(u) I,$$ and the strain tensor is $$\epsilon(u) = \frac12(\nabla u + \nabla u^T).$$ Then you multiply the first equation with a vectorial test function and integrate over the domain $\Omega$ to get $$-\int_\Omega \mathrm{div}\, \sigma(u) \cdot v \,\mathrm{d}x = \int_\Omega f \cdot v \,\mathrm{d}x.$$ Next we use the Green's theorem to get $$-\int_{\partial \Omega} \sigma(u)n \cdot v \,\mathrm{d}s + \int_\Omega \sigma(u) : \nabla v \,\mathrm{d}x = \int_\Omega f \cdot v \,\mathrm{d}x.$$ Since $\sigma$ happens to be symmetric we can equivalently write $$-\int_{\partial \Omega} \sigma(u)n \cdot v \,\mathrm{d}s + \int_\Omega \sigma(u) : (\frac12(\nabla v + \nabla v^T)) \,\mathrm{d}x = \int_\Omega f \cdot v \,\mathrm{d}x,$$ or $$-\int_{\partial \Omega} \sigma(u)n \cdot v \,\mathrm{d}s + \int_\Omega \sigma(u) : \epsilon(v) \,\mathrm{d}x = \int_\Omega f \cdot v \,\mathrm{d}x,$$ which is the weak formulation of linear elasticity.