FEM Kernel
The main FEM loop.
This main FEM loop will assample the differenc matrices that are required in a FEM solver. It takes the following steps:
Loop over all elements.
- Compute the element based integrals.
calculate the quantities in the reference element.
integrate using Quadrature rules and include the mapping from reference to global axis system.
Assamble these element based contributions into a global operator.
Bram Lagerweij COHMAS Mechanical Engineering KAUST 2021
- fem.element_mass(phi_xq, wq_detJ)
Compute the elmement mass matrix.
This matrix is defined as:
\[\begin{split}M &= \int_{\Omega} \phi_j(x) \,\, \phi_i(x) \,\, dV \quad \forall \quad \phi_i, \phi_j \in V_h \\ &= \sum_{e=0}^N \int_{\Omega_e} \phi_j(x) \,\, \phi_i(x) \, dV \quad \forall \quad \phi_i, \phi_j \in V^e_h \\ &= \sum_{e=0}^N M_e\end{split}\]but here we only compute the portion contributed by our current element. Hence we only need to consider the trial functions within each elements. when integrated in reference element coordinates, \(\xi\) this is:
\[M = \int_0^1 \phi_j(\xi) \phi_i(\xi) det(J)dV\]To evaluate these integras Gaussian quadrature is used such that thi integal becomes:
\[M = \sum_{q=0}^{N_q} \phi_j(\xi_q) \phi_i(\xi_q) det(J) w_q\]- Parameters
phi_xq (array_like(float), shape((dofe, num_q))) – For each shape function the value at the quadrature points.
wq_detJ (array_like(float), shape((dofe, num_q))) – Integration weight including the mapping from local to global coordinates.
- Returns
me – Element mass matrix.
- Return type
array_like(float), shape((dofe, dofe))
- fem.element_rhs(phi_xq, wq_detJ, f_xq)
Compute the elmement right hand side vector.
- Parameters
phi_xq (array_like(float), shape((dofs, num_q))) – For each shape function the value at the quadrature points.
f_xq (array_like(float), shape(num_q)) – The value of the right hand side equation evaluated at the quadrature points.
wq_detJ (array_like(float), shape((dofs, num_q))) – Integration weight including the mapping from local to global coordinates.
- Returns
fe – Element right hand side in our system of equations.
- Return type
array_like(float), shape(dofe)
- fem.element_stiffness(invJ_dphi_xq, wq_detJ)
Compute the elmement stiffness matrix.
This matrix is defined as:
\[\begin{split}S &= \int_{\Omega} -\partial_x\phi_j(x) \,\, \partial_x\phi_i(x) \,\, dV \quad \forall \quad \phi_i, \phi_j \in V_h \\ &= \sum_{e=0}^N \int_{\Omega_e} -\partial_x\phi_j(x) \,\, \partial_x\phi_i(x) \,\, dV \quad \forall \quad \phi_i, \phi_j \in V^e_h \\ &= \sum_{e=0}^N S_e\end{split}\]but here we only compute the portion contributed by our current element. Hence we only need to consider the trial functions within each elements. When integrated in reference element coordinates, \(\xi\) this is:
\[S_e = \int_0^1 J^{-1}\, \partial_{\xi} \phi_j(\xi) \,\, J^{-1}\,\partial_{\xi}\phi_i(\xi) \,\,det(J)dV\]To evaluate these integras Gaussian quadrature is used such that thi integal becomes:
\[S_e = \sum_{q=0}^{N_q} J^{-1}\, \partial_{\xi} \phi_j(\xi_q) \,\, J^{-1}\,\partial_{\xi}\phi_i(\xi_q) \,\,det(J) w_q\]- Parameters
invJ_dphi_xq (array_like(float), shape((dofs, num_q))) – For each shape function its derivative value at the quadrature points times the inverse Jacobian.
wq_detJ (array_like(float), shape((dofe, num_q))) – Integration weight including the mapping from local to global coordinates.
- Returns
Se – Element mass matrix.
- Return type
array_like(float), shape((dofe, dofe))
- fem.element_transport(phi_xq, invJ_dphi_xq, wq_detJ)
Compute the elmement transport matrix.
This matrix is defined as:
\[\begin{split}T &= \int_{\Omega} \partial_x\phi_j(x) \,\, \phi_i(x) \,\, dV \quad \forall \quad \phi_i, \phi_j \in V_h \\ &= \sum_{e=0}^N \int_{\Omega_e} \partial_x\phi_j(x) \,\, \phi_i(x) \, dV \quad \forall \quad \phi_i, \phi_j \in V^e_h \\ &= \sum_{e=0}^N T_e\end{split}\]but here we only compute the portion contributed by our current element. Hence we only need to consider the trial functions within each elements. when integrated in reference element coordinates, \(\xi\) this is:
\[T_e = \int_0^1 J^{-1}\,\partial_{\xi}\phi_j(\xi) \,\, \phi_i(\xi) \,\, det(J)dV\]To evaluate these integras Gaussian quadrature is used such that thi integal becomes:
\[T_e = \sum_{q=0}^{N_q} J^{-1}\,\partial_{\xi}\phi_j(\xi_q) \,\, \phi_i(\xi_q) \,\,det(J) w_q\]- Parameters
phi_xq (array_like(float), shape((dofe, num_q))) – For each shape function the value at the quadrature points.
invJ_dphi_xq (array_like(float), shape((dofs, num_q))) – For each shape function its derivative value at the quadrature points times the inverse Jacobian.
wq_detJ (array_like(float), shape((dofe, num_q))) – Integration weight including the mapping from local to global coordinates.
- Returns
Te – Element mass matrix.
- Return type
array_like(float), shape((dofe, dofe))
- fem.interpolate(mesh, u, x_inter)
Obtain the field \(u(x)\) any points x_inter following the FE interpolation.
- Parameters
mesh (Mesh) – The mesh class specifying all discretization.
u (array_like(float), shape(dofs)) – The field u at the degrees of freedom.
x_inter (array_like(float)) – The location where we want to obtain our interpolated field.
- Returns
The field u at the interpolation points x_inter.
- Return type
array_like(float)
- fem.kernel1d(mesh, rhs=None, mass=False, transport=False, stiffness=False)
Create the global FEM system by looping over the elements.
- Parameters
mesh (Mesh) – The mesh class specifying all discretization.
rhs (callable) – Function that acts as our right hand side (nonhomogeneous term), set equal to None if the rhs is zero valued.
mass (bool, optional) – Return a mass matrix. Default is False.
transport (bool, optional) – Return the transport matrix. Default is False.
stiffness (bool, optional) – Return the stiffness matrix. Default is False.
- Returns
f (array_like(float), shape(dofe)) – Global right hand side in our system of equations. Only when rhs != None, None otherwise.
M (COO (value, (row, column))) – Global mass matrix, ready to be converted to COO. Repeating indices do exist. Only mass == True, None otherwise.
T (COO (value, (row, column))) – Global transport matrix, ready to be converted to COO. Repeating indices do exist. Only ‘transport == True`, None otherwise.
S (COO (value, (row, column))) – Global stiffness matrix, ready to be converted to COO. Repeating indices do exist. Only ‘stiffness == True`, None otherwise.