(c) Juan Gomez 2019. Thanks to Universidad EAFIT for support. This material is part of the course Introduction to Finite Element Analysis

# SOLUTION: Interpolation in Two-Dimensional Domains

## By Juan Gomez

### Class activity
**Problem 1**

Extend the 2D interpolation scheme discussed above to the case of a vector valued function in the context of linear elasticity. For that purpose:

* Assume that at each nodal point of the square domain the displacement vector with horizontal and vertical componentes denoted by $u$ and $v$ respectively is known.
* Using these nodal values compute the horizontal and vertical displacement components over the element. Note that the contribution from the nodal dispalcements from the Q-th node can be written like:


$$
\begin{bmatrix}u\\v\end{bmatrix}=\left[\cdots\begin{array}{cc}N^Q&0\\0&N^Q\end{array}\cdots\right]\begin{bmatrix}\vdots\\u^Q\\v^Q\\\vdots\end{bmatrix}\\
$$

where $N^Q$ is the corresponding shape function.


* Using these same nodal values compute the strain field given by:

$$\varepsilon_{xx}=\frac12\left(\frac{\partial u}{\partial x}\right)$$

$$\varepsilon_{yy}=\frac12\left(\frac{\partial v}{\partial y}\right)$$

$$\gamma_{yy}=\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$

It can be observed that in order to find the strain components we need first derivatives of $u$ and $v$ with respect to $x$ and $y$. These derivatives satisfy:

$$\begin{bmatrix}\frac{\partial u}{\partial x}\\\frac{\partial u}{\partial x}\end{bmatrix} = \left[\cdots\begin{array}{cc}\frac{\partial N^Q}{\partial x}&0\\0&\frac{\partial N^Q}{\partial x}\end{array}\cdots\right] \begin{bmatrix}\vdots\\u^Q\\v^Q\\\vdots\end{bmatrix}$$

$$\begin{bmatrix}\frac{\partial u}{\partial y}\\\frac{\partial u}{\partial y}\end{bmatrix} = \left[\cdots\begin{array}{cc}\frac{\partial N^Q}{\partial y}&0\\0&\frac{\partial N^Q}{\partial y}\end{array}\cdots\right] \begin{bmatrix}\vdots\\u^Q\\v^Q\\\vdots\end{bmatrix}$$

thus the vector storing the strain components satisfy:

$$\begin{bmatrix}\varepsilon_{xx}\\\varepsilon_{yy}\\\gamma_{xy}\end{bmatrix}=\begin{bmatrix}&\frac{\partial N^Q}{\partial x}&0&\\\cdots&0&\frac{\partial N^Q}{\partial y}&\cdots\\&\frac{\partial N^Q}{\partial y}&\frac{\partial N^Q}{\partial x}&\end{bmatrix}\begin{bmatrix}\vdots\\u^Q\\v^Q\\\vdots\end{bmatrix}$$.

The expression above can be written in matrix form after considering the contribution from all nodal points in the element as:

$$\left(\varepsilon\right)=\left[B\right]\left(u\right)$$

where we recognize matrix $\left[B\right]$ to be given by:



$$B = \begin{bmatrix}&\frac{\partial N^Q}{\partial x}&0&\\\cdots&0&\frac{\partial N^Q}{\partial y}&\cdots\\&\frac{\partial N^Q}{\partial y}&\frac{\partial N^Q}{\partial x}&\end{bmatrix}$$


```python
%matplotlib notebook        
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
from scipy import interpolate
import math
```

The first subroutine (taken from SolidsPy) stores the interpolation shape functions for the element in the 2D array $N(x , y)$. In this case the element is assumed to be a size 2.0 perfect square.


```python
def sha4(x, y):
    """Shape functions for a 4-noded quad element

    Parameters
    ----------
    x : float
      x coordinate for a point within the element.
    y : float
      y coordinate for a point within the element.

    Returns
    -------
    N : Numpy array
      Array of interpolation functions.

    Examples
    --------
    We can check evaluating at two different points, namely (0, 0) and
    (1, 1). Thus

    >>> N = sha4(0, 0)
    >>> N_ex = np.array([
    ...    [1/4, 0, 1/4, 0, 1/4, 0, 1/4, 0],
    ...    [0, 1/4, 0, 1/4, 0, 1/4, 0, 1/4]])
    >>> np.allclose(N, N_ex)
    True

    and

    >>> N = sha4(1, 1)
    >>> N_ex = np.array([
    ...    [0, 0, 0, 0, 1, 0, 0, 0],
    ...    [0, 0, 0, 0, 0, 1, 0, 0]])
    >>> np.allclose(N, N_ex)
    True

    """
    N = np.zeros((2, 8))
    H = 0.25*np.array(
        [(1 - x)*(1 - y),
         (1 + x)*(1 - y),
         (1 + x)*(1 + y),
         (1 - x)*(1 + y)])
    N[0, ::2] = H
    N[1, 1::2] = H

    return N
```

This second subroutine (adaptade from SolidsPy) stores the shape functions derivatives in the strain matrix $B$ 


```python
def stdm4NQ(x, y):
    """Strain-displacement interpolator B for a 4-noded quad element

    Parameters
    ----------
    r : float
      r component in the natural space.
    s : float
      s component in the natural space.
    coord : ndarray
      Coordinates of the nodes of the element (4, 2).

    Returns
    -------
    ddet : float
      Determinant evaluated at `(r, s)`.
    B : ndarray
      Strain-displacement interpolator evaluated at `(r, s)`.

    """
    nn = 4
    B = np.zeros((3, 2*nn))
    dhdx = 0.25*np.array([
            [y- 1, -y + 1, y + 1, -y - 1],
            [x - 1, -x - 1, x + 1, -x + 1]])

    B[0,  ::2] = 0.5*dhdx[0, :]
    B[1, 1::2] = 0.5* dhdx[1, :]
    B[2,  ::2] = dhdx[1, :]
    B[2, 1::2] = dhdx[0, :]
    return B
```

We are now ready to conduct the interpolation assuming the vector of nodal displacements for the element is known. The response at a single point is computed like:


```python
v = sym.Matrix(8, 1, [-1.0 , -0.5 , +0.5 , +0.4 , -1.0 , -0.5 ,+0.5 ,+0.4])
x = 0.25
y = 0.25
N = sha4(x, y)
B = stdm4NQ(x, y)
u_sol = N*v
e_sol = B*v
print(u_sol)
print(e_sol)
```

    Matrix([[-0.296875000000000], [-0.0781250000000000]])
    Matrix([[-0.0937500000000000], [-0.0562500000000000], [-0.300000000000000]])


**Problem 2**

In a realistic finite element problem and due to geometric irregularities in the computational domain all the elements of the mesh would have different geometric paramters. These differences not only complicate the interpolation process but even create serious problems on coding a systematic approach. In practice every real (distorted) element is transformed to a canonical element in such a way that the interpolation process is conducted in the canonical space. The figure below shows the relation between these two spaces emphasizing the fact that there is a one to one connecion between a point in both spaces. The mathematical details of the transformation are provided in the **Class Notes.**


<center></center>


* Compute the transformation between the distorted and canonical element spaces required to conduct two-dimensional interpolation in a realistic case.

To conduct the transformation we assume that the relation between the physical space $\overrightarrow x$ and canonical (or mathematical space) $\overrightarrow\eta$ is given by:

$$\overrightarrow x=f(\overrightarrow\eta)$$

In the finite element method this transformation is established using also an interpolation scheme. It then follows that the physical space coordinates can be obtained in terms of the mathematical space coordinates like:

$$\begin{bmatrix}x\\y\end{bmatrix}=\left[\cdots\begin{array}{cc}N^Q&0\\0&N^Q\end{array}\cdots\right]\begin{bmatrix}\vdots\\x^Q\\y^Q\\\vdots\end{bmatrix}$$

and where $x^Q$ and $y^Q$ are the coordinates of the Q-th nodal point in the physical space. The corresponding function for the 4-noded element is then


```python
def sha4nat(r,s):
    """
    Compute the shape functions for bi-linear
    square element of size 2.0.
    """
    sh=sym.zeros(4)
    N=sym.zeros(2 , 8)
    sh[0] =(1.0/4.0)*(1 - r)*(1 - s)
    sh[1] =(1.0/4.0)*(1 + r)*(1 - s)
    sh[2] =(1.0/4.0)*(1 + r)*(1 + s)
    sh[3] =(1.0/4.0)*(1 - r)*(1 + s)
    
    N[0,0]=sh[0]
    N[0,2]=sh[1]
    N[0,4]=sh[2]
    N[0,6]=sh[3]
    
    N[1,1]=sh[0]
    N[1,3]=sh[1]
    N[1,5]=sh[2]
    N[1,7]=sh[3]
#
    return N
```

Let us define the nodal points coordinates


```python
coord = sym.Matrix(8, 1, [0.0 , 0.0 , 1.0 , 0.0 , 1.0 , 1.0 , 0.0 , 1.0])
```

Now conduct the interpolation:


```python
r=0.0
s=0.0
N = sha4nat(r,s)
X = N*coord
print(X)
```

    Matrix([[0.500000000000000], [0.500000000000000]])



```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./nb_style.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>

/*
Template for Notebooks for Modelación computacional.

Based on Lorena Barba template available at:

    https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css
*/

/* Fonts */
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

/* Text */
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>









```python

```

ab060c562075da1f6f4e6e9e75b081a5bef56239
ipynb
Jupyter Notebook
solutions/SOL_03_lagrange_2d.ipynb
AppliedMechanics-EAFIT/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2019-11-26T13:28:30.000Z
2022-02-16T17:57:11.000Z
solutions/SOL_03_lagrange_2d.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
solutions/SOL_03_lagrange_2d.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2020-02-17T07:24:59.000Z
2022-03-02T07:54:28.000Z
Qwen/Qwen-72B
1. YES
2. YES