.. !split
.. _fem:approx:2D:
Approximation of functions in 2D
================================
All the concepts and algorithms developed for approximation of 1D functions
:math:`f(x)` can readily be extended to 2D functions :math:`f(x,y)` and 3D functions
:math:`f(x,y,z)`. Basically, the extensions consists of defining basis functions
:math:`{\psi}_i(x,y)` or :math:`{\psi}_i(x,y,z)` over some domain :math:`\Omega`, and
for the least squares and Galerkin methods, the integration is done over
:math:`\Omega`.
As in 1D, the least squares and projection/Galerkin methods
two lead to linear systems
.. math::
\sum_{j\in{\mathcal{I}_s}} A_{i,j}c_j &= b_i,\quad i\in{\mathcal{I}_s},\\
A_{i,j} &= ({\psi}_i,{\psi}_j),\\
b_i &= (f,{\psi}_i),
where the inner product of two functions :math:`f(x,y)` and :math:`g(x,y)` is defined
completely analogously to the 1D case :ref:`(3.2) <Eq:fem:approx:LS:innerprod>`:
.. math::
(f,g) = \int_\Omega f(x,y)g(x,y) dx dy
.. _fem:approx:2D:global:
2D basis functions as tensor products of 1D functions
-----------------------------------------------------
.. index:: tensor product
One straightforward
way to construct a basis in 2D is to combine
1D basis functions. Say we have the 1D vector space
.. _Eq:fem:approx:2D:Vx:
.. math::
:label: fem:approx:2D:Vx
V_x = \mbox{span}\{ \hat{\psi}_0(x),\ldots,\hat{\psi}_{N_x}(x)\}
{\thinspace .}
A similar space for variation in :math:`y` can be defined,
.. _Eq:fem:approx:2D:Vy:
.. math::
:label: fem:approx:2D:Vy
V_y = \mbox{span}\{ \hat{\psi}_0(y),\ldots,\hat{\psi}_{N_y}(y)\}
{\thinspace .}
We can then form 2D basis functions as *tensor products* of 1D basis functions.
.. admonition:: Tensor products
Given two vectors :math:`a=(a_0,\ldots,a_M)` and :math:`b=(b_0,\ldots,b_N)`,
their *outer tensor product*, also called the *dyadic product*,
is :math:`p=a\otimes b`, defined through
.. math::
p_{i,j}=a_ib_j,\quad i=0,\ldots,M,\ j=0,\ldots,N{\thinspace .}
In the tensor terminology,
:math:`a` and :math:`b` are first-order tensors (vectors with one index, also termed
rank-1 tensors), and then their outer
tensor product is a second-order tensor (matrix with two indices, also
termed rank-2 tensor). The
corresponding *inner tensor product* is the well-known scalar or dot
product of two vectors: :math:`p=a\cdot b = \sum_{j=0}^N a_jb_j`. Now,
:math:`p` is a rank-0 tensor.
Tensors are typically represented by arrays in computer code.
In the above example, :math:`a` and :math:`b` are represented by
one-dimensional arrays of length
:math:`M` and :math:`N`, respectively, while :math:`p=a\otimes b` must be represented
by a two-dimensional array of size :math:`M\times N`.
`Tensor products <http://en.wikipedia.org/wiki/Tensor_product>`_ can
be used in a variety of context.
.. The following is from `<http://en.wikipedia.org/wiki/Tensor_product>`_,
.. Notation and examples
Given the vector spaces :math:`V_x` and :math:`V_y` as defined
in :eq:`fem:approx:2D:Vx` and :eq:`fem:approx:2D:Vy`, the
tensor product space :math:`V=V_x\otimes V_y` has a basis formed
as the tensor product of the basis for :math:`V_x` and :math:`V_y`.
That is, if :math:`\left\{ {\varphi}_i(x) \right\}_{i\in{\mathcal{I}_x}}`
and :math:`\left\{ {\varphi}_i(y) \right\}_{i\in {\mathcal{I}_y}}` are basis for
:math:`V_x` and :math:`V_y`, respectively, the elements in the basis for :math:`V` arise
from the tensor product:
:math:`\left\{ {\varphi}_i(x){\varphi}_j(y) \right\}_{i\in {\mathcal{I}_x},j\in {\mathcal{I}_y}}`.
The index sets are :math:`I_x=\{0,\ldots,N_x\}` and :math:`I_y=\{0,\ldots,N_y\}`.
The notation for a basis function in 2D can employ a double index as in
.. math::
{\psi}_{p,q}(x,y) = \hat{\psi}_p(x)\hat{\psi}_q(y),
\quad p\in{\mathcal{I}_x},q\in{\mathcal{I}_y}{\thinspace .}
The expansion for :math:`u` is then written as a double sum
.. math::
u = \sum_{p\in{\mathcal{I}_x}}\sum_{q\in{\mathcal{I}_y}} c_{p,q}{\psi}_{p,q}(x,y){\thinspace .}
Alternatively, we may employ a single index,
.. math::
{\psi}_i(x,y) = \hat{\psi}_p(x)\hat{\psi}_q(y),
and use the standard form for :math:`u`,
.. math::
u = \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x,y){\thinspace .}
The single index is related to the double index through
:math:`i=p N_y + q` or :math:`i=q N_x + p`.
Example: Polynomial basis in 2D
-------------------------------
Suppose we choose :math:`\hat{\psi}_p(x)=x^p`, and try an approximation with
:math:`N_x=N_y=1`:
.. math::
{\psi}_{0,0}=1,\quad {\psi}_{1,0}=x, \quad {\psi}_{0,1}=y,
\quad {\psi}_{1,1}=xy
{\thinspace .}
Using a mapping to one index like :math:`i=q N_x + p`, we get
.. math::
{\psi}_0=1,\quad {\psi}_1=x, \quad {\psi}_2=y,\quad{\psi}_3 =xy
{\thinspace .}
With the specific choice :math:`f(x,y) = (1+x^2)(1+2y^2)` on
:math:`\Omega = [0,L_x]\times [0,L_y]`, we can perform actual calculations:
.. math::
A_{0,0} &= ({\psi}_0,{\psi}_0) = \int_0^{L_y}\int_{0}^{L_x}
{\psi}_0(x,y)^2 dx dy = \int_0^{L_y}\int_{0}^{L_x}dx dy = L_xL_y,\\
A_{1,0} &= ({\psi}_1,{\psi}_0) = \int_0^{L_y}\int_{0}^{L_x} x dxdy =
{\frac{1}{2}}L_x^2L_y,\\
A_{0,1} &= ({\psi}_0,{\psi}_1) = \int_0^{L_y}\int_{0}^{L_x} y dxdy =
{\frac{1}{2}}L_y^2L_x,\\
A_{0,1} &= ({\psi}_0,{\psi}_1) = \int_0^{L_y}\int_{0}^{L_x} xy dxdy =
\int_0^{L_y}ydy \int_{0}^{L_x} xdx =
\frac{1}{4}L_y^2L_x^2
{\thinspace .}
The right-hand side vector has the entries
.. math::
b_{0} &= ({\psi}_0,f) = \int_0^{L_y}\int_{0}^{L_x}1\cdot (1+x^2)(1+2y^2) dxdy\\
&= \int_0^{L_y}(1+2y^2)dy \int_{0}^{L_x} (1+x^2)dx
= (L_y + \frac{2}{3}L_y^3)(L_x + \frac{1}{3}L_x^3)\\
b_{1} &= ({\psi}_1,f) = \int_0^{L_y}\int_{0}^{L_x} x(1+x^2)(1+2y^2) dxdy\\
&=\int_0^{L_y}(1+2y^2)dy \int_{0}^{L_x} x(1+x^2)dx
= (L_y + \frac{2}{3}L_y^3)({\frac{1}{2}}L_x^2 + \frac{1}{4}L_x^4)\\
b_{2} &= ({\psi}_2,f) = \int_0^{L_y}\int_{0}^{L_x} y(1+x^2)(1+2y^2) dxdy\\
&= \int_0^{L_y}y(1+2y^2)dy \int_{0}^{L_x} (1+x^2)dx
= ({\frac{1}{2}}L_y + {\frac{1}{2}}L_y^4)(L_x + \frac{1}{3}L_x^3)\\
b_{3} &= ({\psi}_2,f) = \int_0^{L_y}\int_{0}^{L_x} xy(1+x^2)(1+2y^2) dxdy\\
&= \int_0^{L_y}y(1+2y^2)dy \int_{0}^{L_x} x(1+x^2)dx
= ({\frac{1}{2}}L_y^2 + {\frac{1}{2}}L_y^4)({\frac{1}{2}}L_x^2 + \frac{1}{4}L_x^4)
{\thinspace .}
There is a general pattern in these calculations that we can explore.
An arbitrary matrix entry has the formula
.. math::
A_{i,j} &= ({\psi}_i,{\psi}_j) = \int_0^{L_y}\int_{0}^{L_x}
{\psi}_i{\psi}_j dx dy \\
&= \int_0^{L_y}\int_{0}^{L_x}
{\psi}_{p,q}{\psi}_{r,s} dx dy
= \int_0^{L_y}\int_{0}^{L_x}
\hat{\psi}_p(x)\hat{\psi}_q(y)\hat{\psi}_r(x)\hat{\psi}_s(y) dx dy\\
&= \int_0^{L_y} \hat{\psi}_q(y)\hat{\psi}_s(y)dy
\int_{0}^{L_x} \hat{\psi}_p(x) \hat{\psi}_r(x) dx\\
&= \hat A^{(x)}_{p,r}\hat A^{(y)}_{q,s},
where
.. math::
\hat A^{(x)}_{p,r} = \int_{0}^{L_x} \hat{\psi}_p(x) \hat{\psi}_r(x) dx,
\quad
\hat A^{(y)}_{q,s} = \int_0^{L_y} \hat{\psi}_q(y)\hat{\psi}_s(y)dy,
are matrix entries for one-dimensional approximations. Moreover,
:math:`i=q N_y+q` and :math:`j=s N_y+r`.
With :math:`\hat{\psi}_p(x)=x^p` we have
.. math::
\hat A^{(x)}_{p,r} = \frac{1}{p+r+1}L_x^{p+r+1},\quad
\hat A^{(y)}_{q,s} = \frac{1}{q+s+1}L_y^{q+s+1},
and
.. math::
A_{i,j} = \hat A^{(x)}_{p,r} \hat A^{(y)}_{q,s} =
\frac{1}{p+r+1}L_x^{p+r+1} \frac{1}{q+s+1}L_y^{q+s+1},
for :math:`p,r\in{\mathcal{I}_x}` and :math:`q,s\in{\mathcal{I}_y}`.
Corresponding reasoning for the right-hand side leads to
.. math::
b_i &= ({\psi}_i,f) = \int_0^{L_y}\int_{0}^{L_x}{\psi}_i f\,dxdx\\
&= \int_0^{L_y}\int_{0}^{L_x}\hat{\psi}_p(x)\hat{\psi}_q(y) f\,dxdx\\
&= \int_0^{L_y}\hat{\psi}_q(y) (1+2y^2)dy
\int_0^{L_y}\hat{\psi}_p(x) x^p (1+x^2)dx\\
&= \int_0^{L_y} y^q (1+2y^2)dy
\int_0^{L_y}x^p (1+x^2)dx\\
&= (\frac{1}{q+1} L_y^{q+1} + \frac{2}{q+3}L_y^{q+3})
(\frac{1}{p+1} L_x^{p+1} + \frac{2}{q+3}L_x^{p+3})
Choosing :math:`L_x=L_y=2`, we have
.. math::
A =
\left[\begin{array}{cccc}
4 & 4 & 4 & 4\\
4 & \frac{16}{3} & 4 & \frac{16}{3}\\
4 & 4 & \frac{16}{3} & \frac{16}{3}\\
4 & \frac{16}{3} & \frac{16}{3} & \frac{64}{9}
\end{array}\right],\quad
b = \left[\begin{array}{c}
\frac{308}{9}\\\frac{140}{3}\\44\\60\end{array}\right],
\quad c = \left[
\begin{array}{r}
-\frac{1}{9} \\
\frac{4}{3} \\
- \frac{2}{3} \\
8
\end{array}\right]
{\thinspace .}
Figure :ref:`fem:approx:fe:2D:fig:ubilinear` illustrates the result.
.. _fem:approx:fe:2D:fig:ubilinear:
.. figure:: approx2D_bilinear.png
:width: 800
*Approximation of a 2D quadratic function (left) by a 2D bilinear function (right) using the Galerkin or least squares method*
.. _fem:approx:2D:global:code:
Implementation (3)
-------------------
The ``least_squares`` function from
the section :ref:`fem:approx:global:orth` and/or the
file `approx1D.py <http://tinyurl.com/jvzzcfn/fem/fe_approx1D.py>`_
can with very small modifications solve 2D approximation problems.
First, let ``Omega`` now be a list of the intervals in :math:`x` and :math:`y` direction.
For example, :math:`\Omega = [0,L_x]\times [0,L_y]` can be represented
by ``Omega = [[0, L_x], [0, L_y]]``.
Second, the symbolic integration must be extended to 2D:
.. code-block:: python
import sympy as sp
integrand = psi[i]*psi[j]
I = sp.integrate(integrand,
(x, Omega[0][0], Omega[0][1]),
(y, Omega[1][0], Omega[1][1]))
provided ``integrand`` is an expression involving the ``sympy`` symbols ``x``
and ``y``.
The 2D version of numerical integration becomes
.. code-block:: python
if isinstance(I, sp.Integral):
integrand = sp.lambdify([x,y], integrand)
I = sp.mpmath.quad(integrand,
[Omega[0][0], Omega[0][1]],
[Omega[1][0], Omega[1][1]])
The right-hand side integrals are modified in a similar way.
Third, we must construct a list of 2D basis functions. Here are two
examples based on tensor products of 1D "Taylor-style" polynomials :math:`x^i`
and 1D sine functions :math:`\sin((i+1)\pi x)`:
.. code-block:: python
def taylor(x, y, Nx, Ny):
return [x**i*y**j for i in range(Nx+1) for j in range(Ny+1)]
def sines(x, y, Nx, Ny):
return [sp.sin(sp.pi*(i+1)*x)*sp.sin(sp.pi*(j+1)*y)
for i in range(Nx+1) for j in range(Ny+1)]
The complete code appears in
`approx2D.py <http://tinyurl.com/jvzzcfn/fem/fe_approx2D.py>`_.
The previous hand calculation where a quadratic :math:`f` was approximated by
a bilinear function can be computed symbolically by
>>> from approx2D import *
>>> f = (1+x**2)*(1+2*y**2)
>>> psi = taylor(x, y, 1, 1)
>>> Omega = [[0, 2], [0, 2]]
>>> u = least_squares(f, psi, Omega)
>>> print u
8*x*y - 2*x/3 + 4*y/3 - 1/9
>>> print sp.expand(f)
2*x**2*y**2 + x**2 + 2*y**2 + 1
We may continue with adding higher powers to the basis:
>>> psi = taylor(x, y, 2, 2)
>>> u = least_squares(f, psi, Omega)
>>> print u
2*x**2*y**2 + x**2 + 2*y**2 + 1
>>> print u-f
0
For :math:`N_x\geq 2` and :math:`N_y\geq 2` we recover the exact function :math:`f`, as
expected, since in that case :math:`f\in V` (see
the section :ref:`fem:approx:global:exact`).
.. _fem:approx:3D:global:
Extension to 3D
---------------
Extension to 3D is in principle straightforward once the 2D extension
is understood. The only major difference is that we need the
repeated outer tensor product,
.. math::
V = V_x\otimes V_y\otimes V_z{\thinspace .}
In general, given vectors (first-order tensors)
:math:`a^{(q)} = (a^{(q)}_0,\ldots,a^{(q)}_{N_q}`, :math:`q=0,\ldots,m`,
the tensor product :math:`p=a^{(0)}\otimes\cdots\otimes a^{m}` has
elements
.. math::
p_{i_0,i_1,\ldots,i_m} = a^{(0)}_{i_1}a^{(1)}_{i_1}\cdots a^{(m)}_{i_m}{\thinspace .}
The basis functions in 3D are then
.. math::
{\psi}_{p,q,r}(x,y,z) = \hat{\psi}_p(x)\hat{\psi}_q(y)\hat{\psi}_r(z),
with :math:`p\in{\mathcal{I}_x}`, :math:`q\in{\mathcal{I}_y}`, :math:`r\in{\mathcal{I}_z}`. The expansion of :math:`u` becomes
.. math::
u(x,y,z) = \sum_{p\in{\mathcal{I}_x}}\sum_{q\in{\mathcal{I}_y}}\sum_{r\in{\mathcal{I}_z}} c_{p,q,r}
{\psi}_{p,q,r}(x,y,z){\thinspace .}
A single index can be introduced also here, e.g., :math:`i=N_xN_yr + q_Nx + p`,
:math:`u=\sum_i c_i{\psi}_i(x,y,z)`.
.. admonition:: Use of tensor product spaces
Constructing a multi-dimensional space and basis from tensor products
of 1D spaces is a standard technique when working with global basis
functions. In the world of finite elements, constructing basis functions
by tensor products is much used on quadrilateral and hexahedra cell
shapes, but not on triangles and tetrahedra. Also, the global
finite element basis functions are almost exclusively denoted by a single
index and not by the natural tuple of indices that arises from
tensor products.