v0.14.0
cell_forces.cpp

Implementation with forces potential

In the following code, we solve the problem of identification of forces which result in given displacements on some surface inside the body.

In the experiment, a body is covered on top by transparent gel layer, on the interface between movements of markers (bids) is observed. Strain in the body is generated by cells living on the top layer. Here we looking for forces generated by those cells.

Following implementing when needed could be easily extended to curved surfaces arbitrary located in the body.

How to cite us?

DOI DOI

Local curl-free formulation

This problem can be mathematically described by minimisation problem with the constraints, as follows

\[ \left\{ \begin{array}{l} \textrm{min}\,J(u,\rho) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} \\ \textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\ n\cdot\sigma(u) = \rho \quad \textrm{on}\,S_\rho \\ \end{array} \right. \]

Applying standard procedure, i.e. Lagrange multiplier method, we get

\[ \left\{ \begin{array}{l} \textrm{min}\,J(u,\rho,\Upsilon) = \frac{1}{2} (u,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} + (\Upsilon,\textrm{div}[\sigma])_V \\ n\cdot\sigma = \rho \quad \textrm{on}\,S_\rho \end{array} \right. \]

and applying differentiation by parts

\[ J(u,\rho,\Upsilon) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,n\cdot\sigma)_{S_\rho} \]

and using second constraint, i.e. static constrain

\[ J(u,\rho,\Upsilon) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,\rho)_{S_\rho}. \]

Note that force is enforced in week sense.

Notting that

\[ \sigma = \mathcal{A} \left( \textrm{grad}[u] \right)^{s} \]

and using symmetry of operator \(\mathcal{A}\),

\[ (\textrm{grad}[\Upsilon],\sigma)_V = (\textrm{grad}[\Upsilon],\mathcal{A} \textrm{grad}[u] )_V = (\mathcal{A}^{*}\textrm{grad}[\Upsilon],\textrm{grad}[u] )_V = (\Sigma,\textrm{grad}[u] )_V \]

we finally get

\[ J(u,\rho,\Upsilon) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +(\Sigma,\textrm{grad}[u])_V - (\Upsilon,\rho)_{S_\rho} \]

Calculating derivatives, we can get Euler equations in following form

\[ \left\{ \begin{array}{ll} (Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V &= 0 \\ (\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} &= 0 \\ (\epsilon_\rho \rho,\delta \rho)_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} &= 0 \end{array} \right. \]

and applying finite element discretization, above equations are expressed by linear algebraic equations suitable for computer calculations

\[ \left[ \begin{array}{ccc} S & K_{u \Upsilon} & 0 \\ K_{\Upsilon u} & 0 & -B^\textrm{T} \\ 0 & -B & D \end{array} \right] \left[ \begin{array}{c} u \\ \Upsilon \\ \rho \end{array} \right] = \left[ \begin{array}{c} S u_d \\ 0 \\ 0 \end{array} \right] \]

and finally swapping first two rows, we finally get

\[ \left[ \begin{array}{ccc} K & 0 & -B^\textrm{T} \\ S & K & 0 \\ 0 & -B & D \end{array} \right] \left[ \begin{array}{c} u \\ \Upsilon \\ \rho \end{array} \right] = \left[ \begin{array}{c} 0 \\ S u_d \\ 0 \end{array} \right] \]

where

\[ S=\epsilon_u^{-1} I\\ \epsilon_\rho^\textrm{abs} = \epsilon_\rho^\text{ref}\epsilon_u^{-1} \]

Displacement field is \(\mathbf{u}\) and \(\Upsilon\), force field on top surface is \(\boldsymbol\rho\).

We not yet explored nature of cell forces. Since the cells are attached to the body surface and anything else, the body as results those forces can not be subjected to rigid body motion.

We assume that forces are generated by 2d objects, as results only tangential forces can be produced by those forces if the surface is planar. For non-planar surfaces addition force normal to the surface is present, although the magnitude of this force could be calculated purely from geometrical considerations, thus additional physical equations are not needed.

Assuming that straight and very short pre-stressed fibres generate cell forces; a force field is curl-free. Utilising that forces can be expressed by potential field as follows

\[ \rho = \frac{\partial \Phi_\rho}{\partial \mathbf{x}} \]

where \(\Phi\) is force potential field.

Approximation

For this problem \(H^(\Omega)\) function space is required for \(u\), \(\Upsilon\) and \(\rho\). Note that \(\rho\) is given by derivative of potential field \(\Phi\) which derivatives gives potential forces.

Nonlocal, i.e. weak curl-free

In this variant generalisation of the local model to a weakly enforced force curl-free cell force field is considered. This generalisation is a consequence of the observation that cell has some small but finite size. Moreover, a cell has a complex pre-stressed structure which can transfer shear forces.

We define model like before, however this time we add additional term,

\[ \left\{ \begin{array}{l} \textrm{min}\,J(u,\rho) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho-l^2\textrm{curl}^s [\textrm{curl}^s[\rho]])_{S_\rho} \\ \textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\ n\cdot\sigma = \rho \end{array} \right. \]

and reformulating this we have

\[ \left\{ \begin{array}{l} \textrm{min}\,J(u,\rho) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +\frac{\epsilon_l^{-1}}{2}(\textrm{curl}^s[\rho],\textrm{curl}^s[\rho])_{S_\rho} \\ \textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\ n\cdot\sigma = \rho \end{array} \right. \]

where parameter \(\epsilon_l\) controls curl-free term, if this parameter approach zero, this formulation converge to local variant presented above.

Applying the same procedure like before, set of Euler equations is derived

\[ \left\{ \begin{array}{l} (Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V= 0 \\ (\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} = 0 \\ (\epsilon_\rho \rho,\delta \rho)_{S_\rho}+\epsilon_l^{-1}(\textrm{curl}^s[\rho],\textrm{curl}^s[\delta \rho])_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} = 0 \end{array} \right. \]

That give following system of equations

\[ \left[ \begin{array}{ccc} S & K_{u \Upsilon} & 0 \\ K_{\Upsilon u} & 0 & -B^\textrm{T} \\ 0 & -B & D \end{array} \right] \left[ \begin{array}{c} u \\ \Upsilon \\ \rho \end{array} \right] = \left[ \begin{array}{c} S u_d \\ 0 \\ 0 \end{array} \right] \]

and swapping the first two rows we get

\[ \left[ \begin{array}{ccc} K & 0 & -B^\textrm{T} \\ S & K & 0 \\ 0 & -B & D \end{array} \right] \left[ \begin{array}{c} u \\ \Upsilon \\ \rho \end{array} \right] = \left[ \begin{array}{c} 0 \\ S u_d \\ 0 \end{array} \right] \]

Approximation

For this problem \(H^(\Omega)\) function space is required for \(u\), \(\Upsilon\). Note that force field need to be in \(H(\textrm{curl},S_\rho)\) space.

Field split pre-conditioner

In this implementation, the block structure of the matrix is utilised. The PCFIELDSPLIT (see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html) is utilised. Matrix is stored using nested matrix format MATNEST (see )http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html).

The sub-netsetd matrix is created containing two upper-left blocks,

\[ X = \left[ \begin{array}{cc} K & 0 \\ S & K \end{array} \right] \]

this netsed matrix is blocked in

\[ A= \left[ \begin{array}{cc} X & V \\ H & D \end{array} \right] \]

where

\[ V= \left[ \begin{array}{c} -B^\textrm{T} \\ 0 \end{array} \right] \]

and

\[ H= \left[ \begin{array}{cc} 0 & -B^\textrm{T} \end{array} \right] \]

The system associated with X matrix is solved using PCFIELDSPLIT and multiplicative relaxation scheme. In following implementation sub-matrix, K is factorised only once. The system related to matrix \(A\) composed form matrices \(X,D,V,H\) is solved using Schur complement. Note that here we exploit MoFEM capability to create sub-problems in easy way.

Running code

Approximations of displacements

./map_disp \
-my_data_x ./examples/data_x.csv \
-my_data_y ./examples/data_y.csv \
-my_file ./examples/mesh_cell_force_map.cub \
-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
-lambda 0.01 -my_order 3 -scale 0.05 -cube_size 1 -my_nparts 4

If only one layer is specified, another thin polymer layer could be created automatically using prism elements,

./map_disp_prism \
-my_thinckness 0.01 \
-my_data_x ./examples/data_x.csv \
-my_data_y ./examples/data_y.csv \
-my_file ./examples/mesh_for_prisms.cub \
-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
-lambda 0.01 -my_order 3 -scale 0.05 -cube_size 1 -my_nparts 4

Config file

eps_u = 1e-4
eps_rho = 1e2
eps_l = 1e-2
[block_1]
young_modulus=3.6e-3
poisson_ratio=0.49
[block_2]
young_modulus=0.35
poisson_ratio=0.49
displacement_order=3

Calculating forces

mpirun -np 4 ./cell_forces \
-my_file analysis_mesh.h5m \
-my_order 1 -my_order_force 2 -my_max_post_proc_ref_level 0 \
-my_block_config ./examples/block_config.in \
-ksp_type fgmres -ksp_monitor \
-fieldsplit_1_ksp_type fgmres \
-fieldsplit_1_pc_type lu \
-fieldsplit_1_pc_factor_mat_solver_package mumps \
-fieldsplit_1_ksp_max_it 100 \
-fieldsplit_1_ksp_monitor \
-fieldsplit_0_ksp_type gmres \
-fieldsplit_0_ksp_max_it 25 \
-fieldsplit_0_fieldsplit_0_ksp_type preonly \
-fieldsplit_0_fieldsplit_0_pc_type lu \
-fieldsplit_0_fieldsplit_0_pc_factor_mat_solver_package mumps \
-fieldsplit_0_fieldsplit_1_ksp_type preonly \
-fieldsplit_0_fieldsplit_1_pc_type lu \
-fieldsplit_0_fieldsplit_1_pc_factor_mat_solver_package mumps \
-ksp_atol 1e-6 -ksp_rtol 0 -my_eps_u 1e-4 -my_curl 1

As results of above we get

Example: results

Installation with docker

/**
* \file cell_forces.cpp
* \brief Calculate cell forces
* \example cell_forces.cpp
\tableofcontents
\section cell_potential Implementation with forces potential
In the following code, we solve the problem of identification of forces which
result in given displacements on some surface inside the body.
In the experiment, a body is covered on top by transparent gel layer, on the
interface between movements of markers (bids) is observed. Strain in the body is
generated by cells living on the top layer. Here we looking for forces generated
by those cells.
Following implementing when needed could be easily extended to curved surfaces
arbitrary located in the body.
\ref mofem_citation
\htmlonly
<a href="https://doi.org/10.5281/zenodo.439392"><img
src="https://zenodo.org/badge/DOI/10.5281/zenodo.439392.svg" alt="DOI"></a> <a
href="https://doi.org/10.5281/zenodo.439395"><img
src="https://zenodo.org/badge/DOI/10.5281/zenodo.439395.svg" alt="DOI"></a>
\endhtmlonly
\subsection cell_local Local curl-free formulation
This problem can be mathematically described by minimisation problem with the
constraints, as follows
\f[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +
\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} \\
\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
n\cdot\sigma(u) = \rho \quad \textrm{on}\,S_\rho \\
\end{array}
\right.
\f]
Applying standard procedure, i.e. Lagrange multiplier method, we get
\f[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho,\Upsilon) =
\frac{1}{2} (u,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +
(\Upsilon,\textrm{div}[\sigma])_V \\ n\cdot\sigma = \rho \quad
\textrm{on}\,S_\rho \end{array} \right. \f] and applying differentiation by
parts \f[ J(u,\rho,\Upsilon) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,n\cdot\sigma)_{S_\rho}
\f]
and using second constraint, i.e. static constrain
\f[
J(u,\rho,\Upsilon) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,\rho)_{S_\rho}.
\f]
Note that force is enforced in week sense.
Notting that
\f[
\sigma = \mathcal{A} \left( \textrm{grad}[u] \right)^{s}
\f]
and using symmetry of operator \f$\mathcal{A}\f$,
\f[
(\textrm{grad}[\Upsilon],\sigma)_V =
(\textrm{grad}[\Upsilon],\mathcal{A} \textrm{grad}[u] )_V =
(\mathcal{A}^{*}\textrm{grad}[\Upsilon],\textrm{grad}[u] )_V =
(\Sigma,\textrm{grad}[u] )_V
\f]
we finally get
\f[
J(u,\rho,\Upsilon) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+(\Sigma,\textrm{grad}[u])_V - (\Upsilon,\rho)_{S_\rho}
\f]
Calculating derivatives, we can get Euler equations in following form
\f[
\left\{
\begin{array}{ll}
(Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V &= 0 \\
(\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} &= 0
\\
(\epsilon_\rho \rho,\delta \rho)_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} &= 0
\end{array}
\right.
\f]
and applying finite element discretization, above equations are expressed
by linear algebraic equations suitable for computer calculations
\f[
\left[
\begin{array}{ccc}
S & K_{u \Upsilon} & 0 \\
K_{\Upsilon u} & 0 & -B^\textrm{T} \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
S u_d \\ 0 \\ 0
\end{array}
\right]
\f]
and finally swapping first two rows, we finally get
\f[
\left[
\begin{array}{ccc}
K & 0 & -B^\textrm{T} \\
S & K & 0 \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
0 \\ S u_d \\ 0
\end{array}
\right]
\f]
where
\f[
S=\epsilon_u^{-1} I\\
\epsilon_\rho^\textrm{abs} = \epsilon_\rho^\text{ref}\epsilon_u^{-1}
\f]
Displacement field is \f$\mathbf{u}\f$ and \f$\Upsilon\f$, force field on top
surface is \f$\boldsymbol\rho\f$.
We not yet explored nature of cell forces. Since the cells are attached to the
body surface and anything else, the body as results those forces can not be
subjected to rigid body motion.
We assume that forces are generated by 2d objects, as results only tangential
forces can be produced by those forces if the surface is planar. For non-planar
surfaces addition force normal to the surface is present, although the magnitude
of this force could be calculated purely from geometrical considerations, thus
additional physical equations are not needed.
Assuming that straight and very short pre-stressed fibres generate cell forces;
a force field is curl-free. Utilising that forces can be expressed by potential
field as follows
\f[
\rho = \frac{\partial \Phi_\rho}{\partial \mathbf{x}}
\f]
where \f$\Phi\f$ is force potential field.
\subsubsection cell_local_approx Approximation
For this problem \f$H^(\Omega)\f$ function space is required for \f$u\f$,
\f$\Upsilon\f$ and \f$\rho\f$. Note that \f$\rho\f$ is given by derivative of
potential field \f$\Phi\f$ which derivatives gives potential forces.
\subsection cell_nonlocal Nonlocal, i.e. weak curl-free
In this variant generalisation of the local model to a weakly enforced force
curl-free cell force field is considered. This generalisation is a consequence
of the observation that cell has some small but finite size. Moreover, a cell
has a complex pre-stressed structure which can transfer shear forces.
We define model like before, however this time we add additional term,
\f[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +
\frac{\epsilon_\rho}{2}(\rho,\rho-l^2\textrm{curl}^s
[\textrm{curl}^s[\rho]])_{S_\rho} \\
\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
n\cdot\sigma = \rho
\end{array}
\right.
\f]
and reformulating this we have
\f[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+\frac{\epsilon_l^{-1}}{2}(\textrm{curl}^s[\rho],\textrm{curl}^s[\rho])_{S_\rho}
\\
\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
n\cdot\sigma = \rho
\end{array}
\right.
\f]
where parameter \f$\epsilon_l\f$ controls curl-free term, if this parameter
approach zero, this formulation converge to local variant presented above.
Applying the same procedure like before, set of Euler equations is derived
\f[
\left\{
\begin{array}{l}
(Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V= 0 \\
(\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} = 0 \\
(\epsilon_\rho \rho,\delta
\rho)_{S_\rho}+\epsilon_l^{-1}(\textrm{curl}^s[\rho],\textrm{curl}^s[\delta
\rho])_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} = 0 \end{array} \right. \f] That
give following system of equations \f[ \left[ \begin{array}{ccc}
S & K_{u \Upsilon} & 0 \\
K_{\Upsilon u} & 0 & -B^\textrm{T} \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
S u_d \\ 0 \\ 0
\end{array}
\right]
\f]
and swapping the first two rows we get
\f[
\left[
\begin{array}{ccc}
K & 0 & -B^\textrm{T} \\
S & K & 0 \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
0 \\ S u_d \\ 0
\end{array}
\right]
\f]
\subsubsection cell_nonlocal_approx Approximation
For this problem \f$H^(\Omega)\f$ function space is required for \f$u\f$,
\f$\Upsilon\f$. Note that force field need to be in
\f$H(\textrm{curl},S_\rho)\f$ space.
\section cell_solution Field split pre-conditioner
In this implementation, the block structure of the matrix is utilised. The \e
PCFIELDSPLIT
(see
<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html>)
is utilised. Matrix is stored using nested matrix format \e MATNEST
(see
)<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html>).
The sub-netsetd matrix is created containing two upper-left blocks,
\f[
X = \left[
\begin{array}{cc}
K & 0 \\
S & K
\end{array}
\right]
\f]
this netsed matrix is blocked in
\f[
A= \left[
\begin{array}{cc}
X & V \\
H & D
\end{array}
\right]
\f]
where
\f[
V= \left[
\begin{array}{c}
-B^\textrm{T} \\
0
\end{array}
\right]
\f]
and
\f[
H= \left[
\begin{array}{cc}
0 & -B^\textrm{T}
\end{array}
\right]
\f]
The system associated with X matrix is solved using \e PCFIELDSPLIT and
multiplicative relaxation scheme. In following implementation sub-matrix, K is
factorised only once. The system related to matrix \f$A\f$ composed form
matrices \f$X,D,V,H\f$ is solved using Schur complement. Note that here we
exploit MoFEM capability to create sub-problems in easy way.
\section cell_running_code Running code
Approximations of displacements
\code
./map_disp \
-my_data_x ./examples/data_x.csv \
-my_data_y ./examples/data_y.csv \
-my_file ./examples/mesh_cell_force_map.cub \
-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
-lambda 0.01 -my_order 3 -scale 0.05 -cube_size 1 -my_nparts 4
\endcode
If only one layer is specified, another thin polymer layer could be created
automatically using prism elements,
\code
./map_disp_prism \
-my_thinckness 0.01 \
-my_data_x ./examples/data_x.csv \
-my_data_y ./examples/data_y.csv \
-my_file ./examples/mesh_for_prisms.cub \
-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
-lambda 0.01 -my_order 3 -scale 0.05 -cube_size 1 -my_nparts 4
\endcode
Config file
\include users_modules/cell_engineering/examples/block_config.in
Calculating forces
\code
mpirun -np 4 ./cell_forces \
-my_file analysis_mesh.h5m \
-my_order 1 -my_order_force 2 -my_max_post_proc_ref_level 0 \
-my_block_config ./examples/block_config.in \
-ksp_type fgmres -ksp_monitor \
-fieldsplit_1_ksp_type fgmres \
-fieldsplit_1_pc_type lu \
-fieldsplit_1_pc_factor_mat_solver_package mumps \
-fieldsplit_1_ksp_max_it 100 \
-fieldsplit_1_ksp_monitor \
-fieldsplit_0_ksp_type gmres \
-fieldsplit_0_ksp_max_it 25 \
-fieldsplit_0_fieldsplit_0_ksp_type preonly \
-fieldsplit_0_fieldsplit_0_pc_type lu \
-fieldsplit_0_fieldsplit_0_pc_factor_mat_solver_package mumps \
-fieldsplit_0_fieldsplit_1_ksp_type preonly \
-fieldsplit_0_fieldsplit_1_pc_type lu \
-fieldsplit_0_fieldsplit_1_pc_factor_mat_solver_package mumps \
-ksp_atol 1e-6 -ksp_rtol 0 -my_eps_u 1e-4 -my_curl 1
\endcode
As results of above we get
\image html cell_engineering_forces_example.gif "Example: results" width=600px
\section cell_install Installation with docker
- First install docker as per the instructions here:
<https://docs.docker.com/installation/#installation>
- Install cell user module: \e docker \e pull \e likask/cell_engineering
\todo Improve documentation
\todo Generalization to curved surfaces. That should not be difficult adding
metric tensors and taking into account curvature.
\todo Instead of elastic material we can use GEL model developed in other
module. That will allow to take into account drying and other rheological
effects.
\todo For nonlinear material instead of KSP solver we can add SNES solver
\todo Physical equations for cell mechanics could be added directly to minimized
functional or by constrains.
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
using namespace MoFEM;
#include <Hooke.hpp>
#include <CellForces.hpp>
static char help[] = "-my_block_config set block data\n"
"\n";
namespace CellEngineering {
int oRder;
double yOung;
double pOisson;
BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2) {}
};
} // namespace CellEngineering
using namespace boost::numeric;
using namespace CellEngineering;
#include <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
static int debug = 0;
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
PetscBool flg_block_config, flg_file;
char mesh_file_name[255];
char block_config_file[255];
PetscBool flg_order_force;
PetscInt order = 2;
PetscInt order_force = 2;
PetscBool flg_eps_u, flg_eps_rho, flg_eps_l;
double eps_u = 1e-6;
double eps_rho = 1e-3;
double eps_l = 0;
PetscBool is_curl = PETSC_TRUE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
order, &order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-my_order_force",
"default approximation order for force approximation", "", order_force,
&order_force, &flg_order_force);
CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
"", "block_conf.in", block_config_file, 255,
&flg_block_config);
CHKERR PetscOptionsReal("-my_eps_rho", "foce optimality parameter ", "",
eps_rho, &eps_rho, &flg_eps_rho);
CHKERR PetscOptionsReal("-my_eps_u", "displacement optimality parameter ",
"", eps_u, &eps_u, &flg_eps_u);
CHKERR PetscOptionsReal("-my_eps_l", "displacement optimality parameter ",
"", eps_l, &eps_l, &flg_eps_l);
CHKERR PetscOptionsBool("-my_curl", "use Hcurl space to approximate forces",
"", is_curl, &is_curl, NULL);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
// Reade parameters from line command
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"*** ERROR -my_file (MESH FILE NEEDED)");
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
// Read mesh to MOAB
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
// option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
// Create MoFEM (Joseph) database
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
MeshsetsManager *mmanager_ptr;
CHKERR m_field.query_interface(mmanager_ptr);
// print bcs
CHKERR mmanager_ptr->printDisplacementSet();
CHKERR mmanager_ptr->printForceSet();
// print block sets with materials
CHKERR mmanager_ptr->printMaterialsSet();
// stl::bitset see for more details
BitRefLevel bit_level0;
bit_level0.set(0);
{
Range ents3d;
CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
CHKERR m_field.seed_ref_level(ents3d, bit_level0, false);
}
// set app. order
std::vector<Range> setOrderToEnts(10);
// configure blocks by parsing config file
// it allow to set approximation order for each block independently
Range set_order_ents;
std::map<int, BlockOptionData> block_data;
if (flg_block_config) {
double read_eps_u, read_eps_rho, read_eps_l;
try {
ifstream ini_file(block_config_file);
if (!ini_file.is_open()) {
SETERRQ(PETSC_COMM_SELF, 1,
"*** -my_block_config does not exist ***");
}
// std::cerr << block_config_file << std::endl;
po::variables_map vm;
po::options_description config_file_options;
config_file_options.add_options()(
"eps_u", po::value<double>(&read_eps_u)->default_value(-1))(
"eps_rho", po::value<double>(&read_eps_rho)->default_value(-1))(
"eps_l", po::value<double>(&read_eps_l)->default_value(-1));
std::ostringstream str_order;
str_order << "block_" << it->getMeshsetId() << ".displacement_order";
config_file_options.add_options()(
str_order.str().c_str(),
po::value<int>(&block_data[it->getMeshsetId()].oRder)
->default_value(order));
std::ostringstream str_cond;
str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
config_file_options.add_options()(
str_cond.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].yOung)
->default_value(-1));
std::ostringstream str_capa;
str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
config_file_options.add_options()(
str_capa.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].pOisson)
->default_value(-2));
}
po::parsed_options parsed =
parse_config_file(ini_file, config_file_options, true);
store(parsed, vm);
po::notify(vm);
if (block_data[it->getMeshsetId()].oRder == -1)
continue;
if (block_data[it->getMeshsetId()].oRder == order)
continue;
PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
Range block_ents;
CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
true);
// block_ents = block_ents.subset_by_type(MBTET);
Range nodes;
CHKERR moab.get_connectivity(block_ents, nodes, true);
Range ents_to_set_order, ents3d;
CHKERR moab.get_adjacencies(nodes, 3, false, ents3d,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(ents3d, 2, false, ents_to_set_order,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(ents3d, 1, false, ents_to_set_order,
moab::Interface::UNION);
ents_to_set_order = subtract(
ents_to_set_order, ents_to_set_order.subset_by_type(MBQUAD));
ents_to_set_order = subtract(
ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
set_order_ents.merge(ents3d);
set_order_ents.merge(ents_to_set_order);
setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
set_order_ents);
}
CHKERR m_field.synchronise_entities(set_order_ents, 0);
std::vector<std::string> additional_parameters;
additional_parameters =
collect_unrecognized(parsed.options, po::include_positional);
for (std::vector<std::string>::iterator vit =
additional_parameters.begin();
vit != additional_parameters.end(); vit++) {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"** WARNING Unrecognized option %s\n",
vit->c_str());
}
} catch (const std::exception &ex) {
std::ostringstream ss;
ss << ex.what() << std::endl;
SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
}
if (read_eps_u > 0) {
eps_u = read_eps_u;
};
if (read_eps_rho > 0) {
eps_rho = read_eps_rho;
}
if (read_eps_l > 0) {
eps_l = read_eps_l;
}
}
PetscPrintf(PETSC_COMM_WORLD, "epsU = %6.4e epsRho = %6.4e\n", eps_u,
eps_rho);
// Fields
CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE,
CHKERR m_field.add_field("UPSILON", H1, AINSWORTH_LEGENDRE_BASE, 3,
MB_TAG_SPARSE, MF_ZERO);
if (is_curl) {
} else {
}
// add entitities (by tets) to the field
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "U");
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "UPSILON");
CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "U");
CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "UPSILON");
CHKERR m_field.synchronise_field_entities("UPSILON");
Range vertex_to_fix;
Range edges_to_fix;
Range ents_1st_layer;
// If problem is partitioned meshset culd not exist on this part
if (mmanager_ptr->checkMeshset(202, SIDESET)) {
CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
ents_1st_layer, true);
CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 0,
vertex_to_fix, false);
CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 1, edges_to_fix,
false);
if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Should be one vertex only, but is %d", vertex_to_fix.size());
}
}
CHKERR m_field.synchronise_entities(ents_1st_layer, 0);
ents_1st_layer.subset_by_type(MBTRI), MBTRI, "RHO");
Range ents_2nd_layer;
// If problem is partitioned meshset culd not exist on this part
if (mmanager_ptr->checkMeshset(101, SIDESET)) {
CHKERR mmanager_ptr->getEntitiesByDimension(101, SIDESET, 2,
ents_2nd_layer, true);
}
CHKERR m_field.synchronise_entities(ents_2nd_layer, 0);
for (int oo = 2; oo != setOrderToEnts.size(); oo++) {
if (setOrderToEnts[oo].size() > 0) {
CHKERR m_field.synchronise_entities(setOrderToEnts[oo], 0);
CHKERR m_field.set_field_order(setOrderToEnts[oo], "U", oo);
CHKERR m_field.set_field_order(setOrderToEnts[oo], "UPSILON", oo);
}
}
const int through_thickness_order = 2;
{
Range ents3d;
CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
Range ents;
CHKERR moab.get_adjacencies(ents3d, 2, false, ents,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(ents3d, 1, false, ents,
moab::Interface::UNION);
Range prisms;
CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
{
Range quads;
CHKERR moab.get_adjacencies(prisms, 2, false, quads,
moab::Interface::UNION);
Range prism_tris;
prism_tris = quads.subset_by_type(MBTRI);
quads = subtract(quads, prism_tris);
Range quads_edges;
CHKERR moab.get_adjacencies(quads, 1, false, quads_edges,
moab::Interface::UNION);
Range prism_tris_edges;
CHKERR moab.get_adjacencies(prism_tris, 1, false, prism_tris_edges,
moab::Interface::UNION);
quads_edges = subtract(quads_edges, prism_tris_edges);
prisms.merge(quads);
prisms.merge(quads_edges);
}
ents.merge(ents3d);
ents = subtract(ents, set_order_ents);
ents = subtract(ents, prisms);
CHKERR m_field.synchronise_entities(ents, 0);
CHKERR m_field.synchronise_entities(prisms, 0);
CHKERR m_field.set_field_order(ents, "U", order);
CHKERR m_field.set_field_order(ents, "UPSILON", order);
// approx. order through thickness to 2
CHKERR m_field.set_field_order(prisms, "U", through_thickness_order);
CHKERR m_field.set_field_order(prisms, "UPSILON",
through_thickness_order);
}
CHKERR m_field.set_field_order(0, MBVERTEX, "U", 1);
CHKERR m_field.set_field_order(0, MBVERTEX, "UPSILON", 1);
if (is_curl) {
CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
} else {
CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
CHKERR m_field.set_field_order(0, MBVERTEX, "RHO", 1);
}
// Add elastic element
boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(new Hooke<adouble>());
boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>());
NonlinearElasticElement elastic(m_field, 2);
CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
CHKERR elastic.addElement("ELASTIC", "U");
CHKERR elastic.setOperators("U", "MESH_NODE_POSITIONS", false, true);
// Add prisms
CellEngineering::FatPrism fat_prism_rhs(m_field);
CellEngineering::FatPrism fat_prism_lhs(m_field);
{
CHKERR m_field.add_finite_element("ELASTIC_PRISM", MF_ZERO);
CHKERR m_field.modify_finite_element_add_field_row("ELASTIC_PRISM", "U");
CHKERR m_field.modify_finite_element_add_field_col("ELASTIC_PRISM", "U");
CHKERR m_field.modify_finite_element_add_field_data("ELASTIC_PRISM", "U");
Range ents_2nd_layer;
CHKERR mmanager_ptr->getEntitiesByDimension(2, BLOCKSET, 3,
elastic.setOfBlocks[2].tEts);
elastic.setOfBlocks[2].tEts, MBPRISM, "ELASTIC_PRISM");
MatrixDouble inv_jac;
// right hand side operators
fat_prism_rhs.getOpPtrVector().push_back(
fat_prism_rhs.getOpPtrVector().push_back(
fat_prism_rhs.getOpPtrVector().push_back(
"U", elastic.commonData));
fat_prism_rhs.getOpPtrVector().push_back(
"U", elastic.setOfBlocks[2], elastic.commonData, 2, false, false,
true));
fat_prism_rhs.getOpPtrVector().push_back(
"U", elastic.setOfBlocks[2], elastic.commonData));
// Left hand side operators
fat_prism_lhs.getOpPtrVector().push_back(
fat_prism_lhs.getOpPtrVector().push_back(
fat_prism_lhs.getOpPtrVector().push_back(
"U", elastic.commonData));
fat_prism_lhs.getOpPtrVector().push_back(
"U", elastic.setOfBlocks[2], elastic.commonData, 2, true, false,
true));
fat_prism_lhs.getOpPtrVector().push_back(
"U", "U", elastic.setOfBlocks[2], elastic.commonData));
}
// Update material parameters
int id = it->getMeshsetId();
if (block_data[id].yOung > 0) {
elastic.setOfBlocks[id].E = block_data[id].yOung;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Block %d set Young modulus %3.4g\n", id,
elastic.setOfBlocks[id].E);
}
if (block_data[id].pOisson >= -1) {
elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Block %d set Poisson ratio %3.4g\n", id,
elastic.setOfBlocks[id].PoissonRatio);
}
}
// build field
CHKERR m_field.build_fields();
// Control elements
CHKERR m_field.add_finite_element("KUPSUPS");
CHKERR m_field.modify_finite_element_add_field_row("KUPSUPS", "UPSILON");
CHKERR m_field.modify_finite_element_add_field_col("KUPSUPS", "UPSILON");
CHKERR m_field.modify_finite_element_add_field_data("KUPSUPS", "UPSILON");
CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTET, "KUPSUPS");
CHKERR m_field.add_ents_to_finite_element_by_type(0, MBPRISM, "KUPSUPS");
CHKERR m_field.add_finite_element("DISPLACEMENTS_PENALTY");
CHKERR m_field.modify_finite_element_add_field_row("DISPLACEMENTS_PENALTY",
"UPSILON");
CHKERR m_field.modify_finite_element_add_field_col("DISPLACEMENTS_PENALTY",
"U");
CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
"UPSILON");
CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
"U");
CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
"DISP_X");
CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
"DISP_Y");
CHKERR m_field.add_ents_to_finite_element_by_type(ents_2nd_layer, MBTRI,
"DISPLACEMENTS_PENALTY");
// Add element to calculate residual on 1st layer
CHKERR m_field.add_finite_element("BT");
CHKERR m_field.modify_finite_element_add_field_row("BT", "UPSILON");
CHKERR m_field.modify_finite_element_add_field_data("BT", "UPSILON");
CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
"BT");
CHKERR m_field.add_finite_element("B");
CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
"B");
// Add element to calculate residual on 1st layer
CHKERR m_field.add_finite_element("D");
CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
"D");
// build finite elemnts
// build adjacencies
CHKERR m_field.build_adjacencies(bit_level0);
// Register MOFEM DM
DMType dm_name = "MOFEM";
DM dm_control;
{
// Craete dm_control instance
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
CHKERR DMSetType(dm_control, dm_name);
// set dm_control data structure which created mofem datastructures
CHKERR DMMoFEMCreateMoFEM(dm_control, &m_field, "CONTROL_PROB",
bit_level0);
CHKERR DMSetFromOptions(dm_control);
CHKERR DMMoFEMSetSquareProblem(dm_control, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dm_control, PETSC_TRUE);
// add elements to dm_control
CHKERR DMMoFEMAddElement(dm_control, "ELASTIC");
CHKERR DMMoFEMAddElement(dm_control, "ELASTIC_PRISM");
CHKERR DMMoFEMAddElement(dm_control, "KUPSUPS");
CHKERR DMMoFEMAddElement(dm_control, "DISPLACEMENTS_PENALTY");
CHKERR DMMoFEMAddElement(dm_control, "B");
CHKERR DMMoFEMAddElement(dm_control, "BT");
CHKERR DMMoFEMAddElement(dm_control, "D");
CHKERR DMSetUp(dm_control);
}
MatrixDouble inv_jac;
ublas::matrix<Mat> nested_matrices(2, 2);
ublas::vector<IS> nested_is_rows(2);
ublas::vector<IS> nested_is_cols(2);
for (int i = 0; i != 2; i++) {
nested_is_rows[i] = PETSC_NULL;
nested_is_cols[i] = PETSC_NULL;
for (int j = 0; j != 2; j++) {
nested_matrices(i, j) = PETSC_NULL;
}
}
ublas::matrix<Mat> sub_nested_matrices(2, 2);
ublas::vector<IS> sub_nested_is_rows(2);
ublas::vector<IS> sub_nested_is_cols(2);
for (int i = 0; i != 2; i++) {
sub_nested_is_rows[i] = PETSC_NULL;
sub_nested_is_cols[i] = PETSC_NULL;
for (int j = 0; j != 2; j++) {
sub_nested_matrices(i, j) = PETSC_NULL;
}
}
DM dm_sub_volume_control;
{
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
CHKERR DMSetType(dm_sub_volume_control, dm_name);
// set dm_sub_volume_control data structure which created mofem data
// structures
CHKERR DMMoFEMCreateSubDM(dm_sub_volume_control, dm_control,
"SUB_CONTROL_PROB");
CHKERR DMMoFEMSetSquareProblem(dm_sub_volume_control, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dm_sub_volume_control, "U");
CHKERR DMMoFEMAddSubFieldRow(dm_sub_volume_control, "UPSILON");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_volume_control, "U");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_volume_control, "UPSILON");
// add elements to dm_sub_volume_control
CHKERR DMSetUp(dm_sub_volume_control);
const Problem *prb_ptr;
CHKERR m_field.get_problem("SUB_CONTROL_PROB", &prb_ptr);
boost::shared_ptr<Problem::SubProblemData> sub_data =
prb_ptr->getSubData();
CHKERR sub_data->getRowIs(&nested_is_rows[0]);
CHKERR sub_data->getColIs(&nested_is_cols[0]);
// That will be filled at the end
nested_matrices(0, 0) = PETSC_NULL;
}
{
DM dm_sub_sub_elastic;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
// set dm_sub_sub_elastic data structure which created mofem data
// structures
CHKERR DMMoFEMCreateSubDM(dm_sub_sub_elastic, dm_sub_volume_control,
"ELASTIC_PROB");
CHKERR DMMoFEMSetSquareProblem(dm_sub_sub_elastic, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dm_sub_sub_elastic, "ELASTIC");
CHKERR DMMoFEMAddElement(dm_sub_sub_elastic, "ELASTIC_PRISM");
CHKERR DMMoFEMAddSubFieldRow(dm_sub_sub_elastic, "U");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_sub_elastic, "U");
// add elements to dm_sub_sub_elastic
CHKERR DMSetUp(dm_sub_sub_elastic);
Mat Kuu;
Vec Du, Fu;
CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
CHKERR MatZeroEntries(Kuu);
CHKERR VecZeroEntries(Du);
CHKERR VecZeroEntries(Fu);
CHKERR DMoFEMMeshToLocalVector(dm_sub_sub_elastic, Du, INSERT_VALUES,
SCATTER_REVERSE);
// Manage Dirichlet bC
boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
new DirichletDisplacementBc(m_field, "U", Kuu, Du, Fu));
dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
// preproc
dirihlet_bc_ptr.get());
// internal force vector (to take into account Dirichlet boundary
// conditions)
elastic.getLoopFeRhs().snes_f = Fu;
fat_prism_rhs.snes_f = Fu;
CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
&elastic.getLoopFeRhs());
CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
&fat_prism_rhs);
// elastic element matrix
elastic.getLoopFeLhs().snes_B = Kuu;
fat_prism_lhs.snes_B = Kuu;
CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
&elastic.getLoopFeLhs());
CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
&fat_prism_lhs);
// postproc
dirihlet_bc_ptr.get());
CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(Fu);
CHKERR VecAssemblyEnd(Fu);
CHKERR VecScale(Fu, -1);
CHKERR VecDestroy(&Du);
CHKERR VecDestroy(&Fu);
const Problem *prb_ptr;
CHKERR m_field.get_problem("ELASTIC_PROB", &prb_ptr);
boost::shared_ptr<Problem::SubProblemData> sub_data =
prb_ptr->getSubData();
CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
sub_nested_matrices(0, 0) = Kuu;
IS isUpsilon;
->isCreateFromProblemFieldToOtherProblemField(
"ELASTIC_PROB", "U", ROW, "SUB_CONTROL_PROB", "UPSILON", ROW,
PETSC_NULL, &isUpsilon);
sub_nested_is_rows[1] = isUpsilon;
sub_nested_is_cols[1] = isUpsilon;
sub_nested_matrices(1, 1) = Kuu;
PetscObjectReference((PetscObject)Kuu);
PetscObjectReference((PetscObject)isUpsilon);
// Matrix View
if (debug) {
cerr << "Kuu" << endl;
MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
CHKERR DMDestroy(&dm_sub_sub_elastic);
}
{
DM dm_sub_disp_penalty;
// Craete dm_control instance
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
// set dm_sub_disp_penalty data structure which created mofem
// datastructures
CHKERR DMMoFEMCreateSubDM(dm_sub_disp_penalty, dm_sub_volume_control,
"S_PROB");
CHKERR DMMoFEMSetSquareProblem(dm_sub_disp_penalty, PETSC_FALSE);
CHKERR DMMoFEMAddSubFieldRow(dm_sub_disp_penalty, "UPSILON");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_disp_penalty, "U");
// add elements to dm_sub_disp_penalty
CHKERR DMMoFEMAddElement(dm_sub_disp_penalty, "DISPLACEMENTS_PENALTY");
CHKERR DMSetUp(dm_sub_disp_penalty);
Mat S;
CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
CHKERR MatZeroEntries(S);
CellEngineering::FaceElement face_element(m_field);
face_element.getOpPtrVector().push_back(new OpCellS(S, eps_u));
CHKERR DMoFEMLoopFiniteElements(dm_sub_disp_penalty,
"DISPLACEMENTS_PENALTY", &face_element);
CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
// // Matrix View
if (debug) {
cerr << "S" << endl;
MatView(S, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
const Problem *problem_ptr;
CHKERR m_field.get_problem("S_PROB", &problem_ptr);
boost::shared_ptr<Problem::SubProblemData> sub_data =
problem_ptr->getSubData();
// CHKERR sub_data->getRowIs(&sub_nested_is_rows[1]);
// CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
sub_nested_matrices(1, 0) = S;
CHKERR DMDestroy(&dm_sub_disp_penalty);
}
// Calculate penalty matrix
{
DM dm_sub_force_penalty;
// Craete dm_control instance
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
CHKERR DMSetType(dm_sub_force_penalty, dm_name);
// set dm_sub_force_penalty data structure which created mofem
// datastructures
CHKERR DMMoFEMCreateSubDM(dm_sub_force_penalty, dm_control, "D_PROB");
CHKERR DMMoFEMSetSquareProblem(dm_sub_force_penalty, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dm_sub_force_penalty, "RHO");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_force_penalty, "RHO");
// add elements to dm_sub_force_penalty
CHKERR DMMoFEMAddElement(dm_sub_force_penalty, "D");
CHKERR DMSetUp(dm_sub_force_penalty);
Mat D;
CHKERR DMCreateMatrix(dm_sub_force_penalty, &D);
CHKERR MatZeroEntries(D);
{
CellEngineering::FaceElement face_d_matrix(m_field);
face_d_matrix.getOpPtrVector().push_back(
new OpCalculateInvJacForFace(inv_jac));
if (is_curl) {
face_d_matrix.getOpPtrVector().push_back(
new OpSetInvJacHcurlFace(inv_jac));
face_d_matrix.getOpPtrVector().push_back(
new OpCellCurlD(D, eps_rho / eps_u, eps_l * eps_u));
} else {
face_d_matrix.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac));
face_d_matrix.getOpPtrVector().push_back(
new OpCellPotentialD(D, eps_rho / eps_u));
}
CHKERR DMoFEMLoopFiniteElements(dm_sub_force_penalty, "D",
&face_d_matrix);
}
CHKERR MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);
const Problem *problem_ptr;
CHKERR m_field.get_problem("D_PROB", &problem_ptr);
// Zero rows, force field is given by gradients of potential field, so one
// of the values has to be fixed like for rigid body motion.
if (is_curl == PETSC_FALSE) {
int nb_dofs_to_fix = 0;
int index_to_fix = 0;
if (!vertex_to_fix.empty()) {
boost::shared_ptr<NumeredDofEntity> dof_ptr;
m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
dof_ptr);
if (dof_ptr) {
if (dof_ptr->getPart() == m_field.get_comm_rank()) {
nb_dofs_to_fix = 1;
index_to_fix = dof_ptr->getPetscGlobalDofIdx();
cerr << *dof_ptr << endl;
}
}
}
CHKERR MatZeroRowsColumns(D, nb_dofs_to_fix, &index_to_fix,
eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
} else {
std::vector<int> dofs_to_fix;
for (auto p_eit = edges_to_fix.pair_begin();
p_eit != edges_to_fix.pair_end(); ++p_eit) {
auto bit_number = m_field.get_field_bit_number("RHO");
auto row_dofs = problem_ptr->numeredRowDofsPtr;
auto lo = row_dofs->lower_bound(
FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
auto hi =
bit_number, p_eit->second));
for (; lo != hi; ++lo)
dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
}
CHKERR MatZeroRowsColumns(D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
}
// Matrix View
if (debug) {
cerr << "D" << endl;
MatView(D, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
boost::shared_ptr<Problem::SubProblemData> sub_data =
problem_ptr->getSubData();
CHKERR sub_data->getRowIs(&nested_is_rows[1]);
CHKERR sub_data->getColIs(&nested_is_cols[1]);
nested_matrices(1, 1) = D;
CHKERR DMDestroy(&dm_sub_force_penalty);
}
{
DM dm_sub_force;
// Craete dm_control instance
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
CHKERR DMSetType(dm_sub_force, dm_name);
// set dm_sub_force data structure which created mofem data structures
CHKERR DMMoFEMCreateSubDM(dm_sub_force, dm_control, "FORCES_PROB");
CHKERR DMMoFEMSetSquareProblem(dm_sub_force, PETSC_FALSE);
CHKERR DMMoFEMAddSubFieldRow(dm_sub_force, "RHO");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "U");
CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "UPSILON");
// add elements to dm_sub_force
CHKERR DMMoFEMAddElement(dm_sub_force, "B");
CHKERR DMSetUp(dm_sub_force);
Mat UB, UPSILONB;
CHKERR DMCreateMatrix(dm_sub_force, &UB);
CHKERR MatZeroEntries(UB);
// note be will be transposed in place latter on
CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
CHKERR MatZeroEntries(UPSILONB);
{
CellEngineering::FaceElement face_b_matrices(m_field);
face_b_matrices.getOpPtrVector().push_back(
new OpCalculateInvJacForFace(inv_jac));
if (is_curl) {
face_b_matrices.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac));
face_b_matrices.getOpPtrVector().push_back(
new OpSetInvJacHcurlFace(inv_jac));
face_b_matrices.getOpPtrVector().push_back(new OpCellCurlB(UB, "U"));
face_b_matrices.getOpPtrVector().push_back(
new OpCellCurlB(UPSILONB, "UPSILON"));
} else {
face_b_matrices.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac));
face_b_matrices.getOpPtrVector().push_back(
new OpCellPotentialB(UB, "U"));
face_b_matrices.getOpPtrVector().push_back(
new OpCellPotentialB(UPSILONB, "UPSILON"));
}
CHKERR DMoFEMLoopFiniteElements(dm_sub_force, "B", &face_b_matrices);
}
CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
const Problem *problem_ptr;
CHKERR m_field.get_problem("FORCES_PROB", &problem_ptr);
// Zero rows, force field is given by gradients of potential field, so one
// of the values has to be fixed like for rigid body motion.
if (is_curl == PETSC_FALSE) {
int nb_dofs_to_fix = 0;
int index_to_fix = 0;
if (!vertex_to_fix.empty()) {
boost::shared_ptr<NumeredDofEntity> dof_ptr;
m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
dof_ptr);
if (dof_ptr) {
if (dof_ptr->getPart() == m_field.get_comm_rank()) {
nb_dofs_to_fix = 1;
index_to_fix = dof_ptr->getPetscGlobalDofIdx();
cerr << *dof_ptr << endl;
}
}
}
CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
PETSC_NULL);
CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
PETSC_NULL, PETSC_NULL);
} else {
std::vector<int> dofs_to_fix;
for (auto p_eit = edges_to_fix.pair_begin();
p_eit != edges_to_fix.pair_end(); ++p_eit) {
auto bit_number = m_field.get_field_bit_number("RHO");
auto row_dofs = problem_ptr->numeredRowDofsPtr;
auto lo = row_dofs->lower_bound(
FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
auto hi =
bit_number, p_eit->second));
for (; lo != hi; ++lo)
dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
}
CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
PETSC_NULL, PETSC_NULL);
CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
0, PETSC_NULL, PETSC_NULL);
}
Mat UBT;
CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
CHKERR MatDestroy(&UB);
// Matrix View
if (debug) {
cerr << "UBT" << endl;
MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
boost::shared_ptr<Problem::SubProblemData> sub_data =
problem_ptr->getSubData();
// CHKERR sub_data->getColIs(&nested_is_rows[0]);
// CHKERR sub_data->getRowIs(&nested_is_cols[1]);
nested_matrices(0, 1) = UBT;
if (debug) {
cerr << "UPSILONB" << endl;
MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
// CHKERR sub_data->getRowIs(&nested_is_rows[1]);
// CHKERR sub_data->getColIs(&nested_is_cols[0]);
nested_matrices(1, 0) = UPSILONB;
CHKERR DMDestroy(&dm_sub_force);
}
Mat SubA;
CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
&sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
&SubA);
nested_matrices(0, 0) = SubA;
CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
if (debug) {
cerr << "Nested SubA" << endl;
MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
}
Mat A;
CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
&nested_is_cols[0], &nested_matrices(0, 0), &A);
CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
if (debug) {
cerr << "Nested A" << endl;
MatView(A, PETSC_VIEWER_STDOUT_WORLD);
}
Vec D, F;
CHKERR DMCreateGlobalVector(dm_control, &D);
CHKERR DMCreateGlobalVector(dm_control, &F);
// Asseble the right hand side vector
{
CellEngineering::FaceElement face_element(m_field);
face_element.getOpPtrVector().push_back(new OpGetDispX(common_data));
face_element.getOpPtrVector().push_back(new OpGetDispY(common_data));
face_element.getOpPtrVector().push_back(
new OpCell_g(F, eps_u, common_data));
CHKERR DMoFEMLoopFiniteElements(dm_control, "DISPLACEMENTS_PENALTY",
&face_element);
CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(F);
CHKERR VecAssemblyEnd(F);
}
KSP solver;
// KSP solver;
{
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetDM(solver, dm_control);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetOperators(solver, A, A);
CHKERR KSPSetDMActive(solver, PETSC_FALSE);
CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
PC pc;
CHKERR KSPGetPC(solver, &pc);
CHKERR PCSetType(pc, PCFIELDSPLIT);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
CHKERR PCSetOperators(pc, A, A);
CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
CHKERR PCSetUp(pc);
KSP *sub_ksp;
PetscInt n;
CHKERR PCFieldSplitGetSubKSP(pc, &n, &sub_ksp);
{
PC sub_pc_0;
CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
// CHKERR
// PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER);
// CHKERR PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR);
CHKERR PCSetUp(sub_pc_0);
}
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Only works with pre-conditioner PCFIELDSPLIT");
}
CHKERR KSPSetUp(solver);
}
// Solve system of equations
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm_control, D, INSERT_VALUES,
SCATTER_REVERSE);
if (debug) {
CHKERR VecView(D, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
// Clean sub matrices and sub indices
for (int i = 0; i != 2; i++) {
if (sub_nested_is_rows[i]) {
CHKERR ISDestroy(&sub_nested_is_rows[i]);
}
if (sub_nested_is_cols[i]) {
CHKERR ISDestroy(&sub_nested_is_cols[i]);
}
for (int j = 0; j != 2; j++) {
if (sub_nested_matrices(i, j)) {
CHKERR MatDestroy(&sub_nested_matrices(i, j));
}
}
}
for (int i = 0; i != 2; i++) {
if (nested_is_rows[i]) {
CHKERR ISDestroy(&nested_is_rows[i]);
}
if (nested_is_cols[i]) {
CHKERR ISDestroy(&nested_is_cols[i]);
}
for (int j = 0; j != 2; j++) {
if (nested_matrices(i, j)) {
CHKERR MatDestroy(&nested_matrices(i, j));
}
}
}
CHKERR MatDestroy(&SubA);
CHKERR MatDestroy(&A);
CHKERR VecDestroy(&D);
CHKERR VecDestroy(&F);
CHKERR DMDestroy(&dm_sub_volume_control);
PostProcVolumeOnRefinedMesh post_proc(m_field);
{
CHKERR post_proc.generateReferenceElementMesh();
CHKERR post_proc.addFieldValuesPostProc("U");
CHKERR post_proc.addFieldValuesGradientPostProc("U");
// add postpocessing for sresses
post_proc.getOpPtrVector().push_back(new PostProcHookStress(
m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "U",
post_proc.commonData, &elastic.setOfBlocks));
CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC", &post_proc);
CHKERR post_proc.writeFile("out.h5m");
elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
elastic.getLoopFeEnergy().eNergy = 0;
CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC",
&elastic.getLoopFeEnergy());
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
elastic.getLoopFeEnergy().eNergy);
}
{
PostProcFaceOnRefinedMesh post_proc_face(m_field);
CHKERR post_proc_face.generateReferenceElementMesh();
post_proc_face.getOpPtrVector().push_back(
new OpCalculateInvJacForFace(inv_jac));
if (is_curl) {
post_proc_face.getOpPtrVector().push_back(
new OpSetInvJacHcurlFace(inv_jac));
post_proc_face.getOpPtrVector().push_back(
new OpVirtualCurlRho("RHO", common_data));
} else {
post_proc_face.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac));
post_proc_face.getOpPtrVector().push_back(
new OpVirtualPotentialRho("RHO", common_data));
}
post_proc_face.getOpPtrVector().push_back(
new PostProcTraction(m_field, post_proc_face.postProcMesh,
post_proc_face.mapGaussPts, common_data));
CHKERR DMoFEMLoopFiniteElements(dm_control, "D", &post_proc_face);
CHKERR post_proc_face.writeFile("out_tractions.h5m");
}
CHKERR DMDestroy(&dm_control);
}
return 0;
}
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
SIDESET
@ SIDESET
Definition: definitions.h:160
debug
static int debug
Definition: cell_forces.cpp:430
CellEngineering::PostProcTraction
Shave results on mesh tags for post-processing.
Definition: CellForces.hpp:600
CellEngineering
Definition: cell_forces.cpp:412
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition: UserDataOperators.hpp:2978
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
MoFEM::UnknownInterface::query_interface
virtual MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
CellEngineering::FaceElement
Definition: CellForces.hpp:48
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::FieldEntity::getLoLocalEntityBitNumber
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition: FieldEntsMultiIndices.hpp:247
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
CellEngineering::OpCellCurlB
Calculate and assemble Z matrix.
Definition: CellForces.hpp:392
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
scale
double scale
Definition: plastic.cpp:119
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3145
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
help
static char help[]
Definition: cell_forces.cpp:409
CellEngineering::OpVirtualCurlRho
Post-process tractions.
Definition: CellForces.hpp:528
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:304
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3172
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::MeshsetsManager::getEntitiesByDimension
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimension
Definition: MeshsetsManager.cpp:669
main
int main(int argc, char *argv[])
Definition: cell_forces.cpp:432
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
Hooke.hpp
Implementation of linear elastic material.
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
CellEngineering::OpCellCurlD
Calculate and assemble Z matrix.
Definition: CellForces.hpp:315
CellEngineering::OpCellPotentialD
Calculate and assemble D matrix.
Definition: CellForces.hpp:73
MoFEM::Problem::getSubData
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
Definition: ProblemsMultiIndices.hpp:119
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
Definition: LoopMethods.hpp:145
MoFEM::DeprecatedCoreInterface::synchronise_entities
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:41
CellEngineering::FatPrism
Definition: CellForces.hpp:21
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
CellEngineering::OpCellS
Calculate and assemble S matrix.
Definition: CellForces.hpp:201
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
CellEngineering::OpGetDispX
Definition: CellForces.hpp:57
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:114
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:325
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
convert.n
n
Definition: convert.py:82
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:73
MoFEM::OpSetInvJacH1ForFace
Definition: UserDataOperators.hpp:2918
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
MoFEM::DeprecatedCoreInterface::seed_ref_level
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Definition: DeprecatedCoreInterface.cpp:24
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:290
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
std
Definition: enable_if.hpp:5
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
CellEngineering::OpCellPotentialB
Calculate and assemble B matrix.
Definition: CellForces.hpp:140
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
CellEngineering::CommonData
Definition: CellForces.hpp:31
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
CellEngineering::OpGetDispY
Definition: CellForces.hpp:65
MoFEM::DeprecatedCoreInterface::synchronise_field_entities
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:47
MoFEM::FieldEntity::getHiLocalEntityBitNumber
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition: FieldEntsMultiIndices.hpp:258
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
BlockOptionData
Definition: elasticity.cpp:62
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::Problem::getDofByNameEntAndEntDofIdx
MoFEMErrorCode getDofByNameEntAndEntDofIdx(const int field_bit_number, const EntityHandle ent, const int ent_dof_idx, const RowColData row_or_col, boost::shared_ptr< NumeredDofEntity > &dof_ptr) const
get DOFs from problem
Definition: ProblemsMultiIndices.cpp:132
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:360
CellEngineering::OpCell_g
Calculate and assemble g vector.
Definition: CellForces.hpp:267
CellEngineering::OpVirtualPotentialRho
Post-process tractions.
Definition: CellForces.hpp:460
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
CellForces.hpp