v0.5.86
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{div}^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

cell_engineering_forces_example.gif
Example: results

Installation with docker

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.
/**
* \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{div}^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 {
struct BlockOptionData {
int oRder;
double yOung;
double pOisson;
oRder(-1),
yOung(-1),
pOisson(-2) {
}
};
}
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[]) {
PetscInitialize(&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;
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Elastic Config","none"); CHKERRQ(ierr);
ierr = PetscOptionsString(
"-my_file",
"mesh file name","",
"mesh.h5m",mesh_file_name,255,&flg_file
ierr = PetscOptionsInt(
"-my_order",
"default approximation order","",
order,&order,PETSC_NULL
ierr = PetscOptionsInt(
"-my_order_force",
"default approximation order for force approximation","",
order_force,&order_force,&flg_order_force
ierr = PetscOptionsString("-my_block_config",
"elastic configure file name","",
"block_conf.in",block_config_file,255,&flg_block_config
ierr = PetscOptionsReal(
"-my_eps_rho",
"foce optimality parameter ","",
eps_rho,&eps_rho,&flg_eps_rho
ierr = PetscOptionsReal(
"-my_eps_u",
"displacement optimality parameter ","",
eps_u,&eps_u,&flg_eps_u
ierr = PetscOptionsReal(
"-my_eps_l",
"displacement optimality parameter ","",
eps_l,&eps_l,&flg_eps_l
ierr = 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 = "";
rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
//Create MoFEM (Joseph) database
MoFEM::Core core(moab);
MoFEM::Interface& m_field = core;
MeshsetsManager *mmanager_ptr;
ierr = m_field.query_interface(mmanager_ptr); CHKERRQ(ierr);
//print bcs
ierr = mmanager_ptr->printDisplacementSet(); CHKERRQ(ierr);
ierr = mmanager_ptr->printForceSet(); CHKERRQ(ierr);
//print block sets with materials
ierr = mmanager_ptr->printMaterialsSet(); CHKERRQ(ierr);
// stl::bitset see for more details
BitRefLevel bit_level0;
bit_level0.set(0);
{
Range ents3d;
rval = moab.get_entities_by_dimension(0,3,ents3d); CHKERRQ_MOAB(rval);
ierr = m_field.seed_ref_level(ents3d,bit_level0,false); CHKERRQ(ierr);
}
//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;
rval = moab.get_entities_by_handle(it->getMeshset(),block_ents,true); CHKERRQ_MOAB(rval);
// block_ents = block_ents.subset_by_type(MBTET);
Range nodes;
rval = moab.get_connectivity(block_ents,nodes,true); CHKERRQ(ierr);
Range ents_to_set_order,ents3d;
rval = moab.get_adjacencies(nodes,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
rval = moab.get_adjacencies(ents3d,2,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
rval = moab.get_adjacencies(ents3d,1,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
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);
}
ierr = m_field.synchronise_entities(set_order_ents,0); CHKERRQ(ierr);
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++) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"** WARNING Unrecognized option %s\n",vit->c_str()); CHKERRQ(ierr);
}
} 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
ierr = m_field.add_field("U",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
ierr = m_field.add_field("UPSILON",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
if(is_curl) {
} else {
}
//add entitities (by tets) to the field
ierr = m_field.add_ents_to_field_by_type(0,MBTET,"U"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_field_by_type(0,MBTET,"UPSILON"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_field_by_type(0,MBPRISM,"U"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_field_by_type(0,MBPRISM,"UPSILON"); CHKERRQ(ierr);
ierr = m_field.synchronise_field_entities("UPSILON"); CHKERRQ(ierr);
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)) {
ierr = mmanager_ptr->getEntitiesByDimension(202,SIDESET,2,ents_1st_layer,true); CHKERRQ(ierr);
ierr = mmanager_ptr->getEntitiesByDimension(202,SIDESET,0,vertex_to_fix,false); CHKERRQ(ierr);
ierr = mmanager_ptr->getEntitiesByDimension(202,SIDESET,1,edges_to_fix,false); CHKERRQ(ierr);
if(vertex_to_fix.size()!=1 && !vertex_to_fix.empty()) {
SETERRQ1(
PETSC_COMM_WORLD,
"Should be one vertex only, but is %d",
vertex_to_fix.size()
);
}
}
ierr = m_field.synchronise_entities(ents_1st_layer,0); CHKERRQ(ierr);
ierr = m_field.add_ents_to_field_by_type(ents_1st_layer.subset_by_type(MBTRI),MBTRI,"RHO"); CHKERRQ(ierr);
Range ents_2nd_layer;
// If problem is partitioned meshset culd not exist on this part
if(mmanager_ptr->checkMeshset(101,SIDESET)) {
ierr = mmanager_ptr->getEntitiesByDimension(101,SIDESET,2,ents_2nd_layer,true); CHKERRQ(ierr);
}
ierr = m_field.synchronise_entities(ents_2nd_layer,0); CHKERRQ(ierr);
for(int oo = 2;oo!=setOrderToEnts.size();oo++) {
if(setOrderToEnts[oo].size()>0) {
ierr = m_field.synchronise_entities(setOrderToEnts[oo],0); CHKERRQ(ierr);
ierr = m_field.set_field_order(setOrderToEnts[oo],"U",oo); CHKERRQ(ierr);
ierr = m_field.set_field_order(setOrderToEnts[oo],"UPSILON",oo); CHKERRQ(ierr);
}
}
const int through_thickness_order = 2;
{
Range ents3d;
rval = moab.get_entities_by_dimension(0,3,ents3d); CHKERRQ_MOAB(rval);
Range ents;
rval = moab.get_adjacencies(ents3d,2,false,ents,moab::Interface::UNION); CHKERRQ_MOAB(rval);
rval = moab.get_adjacencies(ents3d,1,false,ents,moab::Interface::UNION); CHKERRQ_MOAB(rval);
Range prisms;
rval = moab.get_entities_by_type(0,MBPRISM,prisms); CHKERRQ_MOAB(rval);
{
Range quads;
rval = moab.get_adjacencies(prisms,2,false,quads,moab::Interface::UNION); CHKERRQ_MOAB(rval);
Range prism_tris;
prism_tris = quads.subset_by_type(MBTRI);
quads = subtract(quads,prism_tris);
Range quads_edges;
rval = moab.get_adjacencies(quads,1,false,quads_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
Range prism_tris_edges;
rval = moab.get_adjacencies(prism_tris,1,false,prism_tris_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
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);
ierr = m_field.synchronise_entities(ents,0); CHKERRQ(ierr);
ierr = m_field.synchronise_entities(prisms,0); CHKERRQ(ierr);
ierr = m_field.set_field_order(ents,"U",order); CHKERRQ(ierr);
ierr = m_field.set_field_order(ents,"UPSILON",order); CHKERRQ(ierr);
// approx. order through thickness to 2
ierr = m_field.set_field_order(prisms,"U",through_thickness_order); CHKERRQ(ierr);
ierr = m_field.set_field_order(prisms,"UPSILON",through_thickness_order); CHKERRQ(ierr);
}
ierr = m_field.set_field_order(0,MBVERTEX,"U",1); CHKERRQ(ierr);
ierr = m_field.set_field_order(0,MBVERTEX,"UPSILON",1); CHKERRQ(ierr);
if(is_curl) {
ierr = m_field.set_field_order(0,MBTRI,"RHO",order_force); CHKERRQ(ierr);
ierr = m_field.set_field_order(0,MBEDGE,"RHO",order_force); CHKERRQ(ierr);
} else {
ierr = m_field.set_field_order(0,MBTRI,"RHO",order_force); CHKERRQ(ierr);
ierr = m_field.set_field_order(0,MBEDGE,"RHO",order_force); CHKERRQ(ierr);
ierr = m_field.set_field_order(0,MBVERTEX,"RHO",1); CHKERRQ(ierr);
}
// 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);
ierr = elastic.setBlocks(hooke_double_ptr,hooke_adouble_ptr); CHKERRQ(ierr);
ierr = elastic.addElement("ELASTIC","U"); CHKERRQ(ierr);
ierr = elastic.setOperators("U","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
// Add prisms
CellEngineering::FatPrism fat_prism_rhs(m_field);
CellEngineering::FatPrism fat_prism_lhs(m_field);
{
ierr = m_field.add_finite_element("ELASTIC_PRISM",MF_ZERO); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("ELASTIC_PRISM","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("ELASTIC_PRISM","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("ELASTIC_PRISM","U"); CHKERRQ(ierr);
Range ents_2nd_layer;
ierr = mmanager_ptr->getEntitiesByDimension(2,BLOCKSET,3,elastic.setOfBlocks[2].tEts); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(elastic.setOfBlocks[2].tEts,MBPRISM,"ELASTIC_PRISM"); CHKERRQ(ierr);
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
fat_prism_rhs.getOpPtrVector().push_back
"U",elastic.setOfBlocks[2],elastic.commonData,2,false,false,true)
);
fat_prism_rhs.getOpPtrVector().push_back
// Left hand side operators
fat_prism_lhs.getOpPtrVector().push_back
fat_prism_lhs.getOpPtrVector().push_back
fat_prism_lhs.getOpPtrVector().push_back
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;
ierr = PetscPrintf(
PETSC_COMM_WORLD,
"Block %d set Young modulus %3.4g\n",
id,elastic.setOfBlocks[id].E
); CHKERRQ(ierr);
}
if(block_data[id].pOisson>=-1) {
elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
ierr = PetscPrintf(
PETSC_COMM_WORLD,"Block %d set Poisson ratio %3.4g\n",id,elastic.setOfBlocks[id].PoissonRatio
); CHKERRQ(ierr);
}
}
//build field
ierr = m_field.build_fields(); CHKERRQ(ierr);
// Control elements
ierr = m_field.add_finite_element("KUPSUPS"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("KUPSUPS","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("KUPSUPS","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("KUPSUPS","UPSILON"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(0,MBTET,"KUPSUPS"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(0,MBPRISM,"KUPSUPS"); CHKERRQ(ierr);
ierr = m_field.add_finite_element("DISPLACEMENTS_PENALTY"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("DISPLACEMENTS_PENALTY","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("DISPLACEMENTS_PENALTY","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY","DISP_X"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY","DISP_Y"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(ents_2nd_layer,MBTRI,"DISPLACEMENTS_PENALTY"); CHKERRQ(ierr);
// Add element to calculate residual on 1st layer
ierr = m_field.add_finite_element("BT"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("BT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("BT","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("BT","RHO"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("BT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("BT","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("BT","RHO"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(ents_1st_layer,MBTRI,"BT"); CHKERRQ(ierr);
ierr = m_field.add_finite_element("B"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("B","RHO"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("B","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("B","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("B","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("B","UPSILON"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("B","RHO"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(ents_1st_layer,MBTRI,"B"); CHKERRQ(ierr);
// Add element to calculate residual on 1st layer
ierr = m_field.add_finite_element("D"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("D","RHO"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("D","RHO"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("D","RHO"); CHKERRQ(ierr);
ierr = m_field.add_ents_to_finite_element_by_type(ents_1st_layer,MBTRI,"D"); CHKERRQ(ierr);
// FIXME This is dirty fix, If there are entities in the field on partition, but
// there is no element owned by given partition that DOFs will be not added to the problem.
// This solve that issue. Mesh partitioning should be done better. This is general problem
// when 3d elements are mixed with 1d/2d elements on the surface where surface elements
// uses field which is not present in 3d elements.
// {
// Range force_field_ents;
// ierr = m_field.get_field_entities_by_handle("RHO",force_field_ents); CHKERRQ(ierr);
// Range vforce_ents;
// ierr = m_field.get_finite_element_entities_by_handle("B",vforce_ents); CHKERRQ(ierr);
// Range fe_ents_adj;
// for(int d = 2;d>=0;d--) {
// rval = moab.get_adjacencies(vforce_ents,d,false,fe_ents_adj,moab::Interface::UNION); CHKERRQ_MOAB(rval);
// }
// fe_ents_adj.merge(vforce_ents);
// force_field_ents = subtract(force_field_ents,fe_ents_adj);
// ierr = m_field.synchronise_entities(fe_ents_adj,0); CHKERRQ(ierr);
// ierr = m_field.add_ents_to_finite_element_by_type(fe_ents_adj,MBVERTEX,"B"); CHKERRQ(ierr);
// ierr = m_field.add_ents_to_finite_element_by_type(fe_ents_adj,MBEDGE,"B"); CHKERRQ(ierr);
// ierr = m_field.add_ents_to_finite_element_by_type(fe_ents_adj,MBTRI,"B"); CHKERRQ(ierr);
// }
//build finite elemnts
ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
//build adjacencies
ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
// Register MOFEM DM
DMType dm_name = "MOFEM";
ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
DM dm_control;
{
// Craete dm_control instance
ierr = DMCreate(PETSC_COMM_WORLD,&dm_control);CHKERRQ(ierr);
ierr = DMSetType(dm_control,dm_name);CHKERRQ(ierr);
//set dm_control data structure which created mofem datastructures
ierr = DMMoFEMCreateMoFEM(dm_control,&m_field,"CONTROL_PROB",bit_level0); CHKERRQ(ierr);
ierr = DMSetFromOptions(dm_control); CHKERRQ(ierr);
ierr = DMMoFEMSetSquareProblem(dm_control,PETSC_TRUE); CHKERRQ(ierr);
ierr = DMMoFEMSetIsPartitioned(dm_control,PETSC_TRUE); CHKERRQ(ierr);
//add elements to dm_control
ierr = DMMoFEMAddElement(dm_control,"ELASTIC"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_control,"ELASTIC_PRISM"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_control,"KUPSUPS"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_control,"DISPLACEMENTS_PENALTY"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_control,"B"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_control,"BT"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_control,"D"); CHKERRQ(ierr);
ierr = DMSetUp(dm_control); CHKERRQ(ierr);
}
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;
{
ierr = DMCreate(PETSC_COMM_WORLD,&dm_sub_volume_control);CHKERRQ(ierr);
ierr = DMSetType(dm_sub_volume_control,dm_name);CHKERRQ(ierr);
//set dm_sub_volume_control data structure which created mofem data structures
ierr = DMMoFEMCreateSubDM(dm_sub_volume_control,dm_control,"SUB_CONTROL_PROB"); CHKERRQ(ierr);
ierr = DMMoFEMSetSquareProblem(dm_sub_volume_control,PETSC_TRUE); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldRow(dm_sub_volume_control,"U"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldRow(dm_sub_volume_control,"UPSILON"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_volume_control,"U"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_volume_control,"UPSILON"); CHKERRQ(ierr);
// add elements to dm_sub_volume_control
ierr = DMSetUp(dm_sub_volume_control); CHKERRQ(ierr);
const Problem *prb_ptr;
ierr = m_field.get_problem("SUB_CONTROL_PROB",&prb_ptr); CHKERRQ(ierr);
boost::shared_ptr<Problem::SubProblemData> sub_data = prb_ptr->getSubData();
ierr = sub_data->getRowIs(&nested_is_rows[0]); CHKERRQ(ierr);
ierr = sub_data->getColIs(&nested_is_cols[0]); CHKERRQ(ierr);
// That will be filled at the end
nested_matrices(0,0) = PETSC_NULL;
}
{
DM dm_sub_sub_elastic;
ierr = DMCreate(PETSC_COMM_WORLD,&dm_sub_sub_elastic);CHKERRQ(ierr);
ierr = DMSetType(dm_sub_sub_elastic,dm_name);CHKERRQ(ierr);
//set dm_sub_sub_elastic data structure which created mofem data structures
ierr = DMMoFEMCreateSubDM(dm_sub_sub_elastic,dm_sub_volume_control,"ELASTIC_PROB"); CHKERRQ(ierr);
ierr = DMMoFEMSetSquareProblem(dm_sub_sub_elastic,PETSC_TRUE); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_sub_sub_elastic,"ELASTIC"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm_sub_sub_elastic,"ELASTIC_PRISM"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldRow(dm_sub_sub_elastic,"U"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_sub_elastic,"U"); CHKERRQ(ierr);
// add elements to dm_sub_sub_elastic
ierr = DMSetUp(dm_sub_sub_elastic); CHKERRQ(ierr);
Mat Kuu;
Vec Du,Fu;
ierr = DMCreateMatrix(dm_sub_sub_elastic,&Kuu); CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm_sub_sub_elastic,&Du); CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm_sub_sub_elastic,&Fu); CHKERRQ(ierr);
ierr = MatZeroEntries(Kuu); CHKERRQ(ierr);
ierr = VecZeroEntries(Du); CHKERRQ(ierr);
ierr = VecZeroEntries(Fu); CHKERRQ(ierr);
ierr = DMoFEMMeshToLocalVector(dm_sub_sub_elastic,Du,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
// 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
ierr = DMoFEMPreProcessFiniteElements(dm_sub_sub_elastic,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
// internal force vector (to take into account Dirichlet boundary conditions)
elastic.getLoopFeRhs().snes_f = Fu;
fat_prism_rhs.snes_f = Fu;
ierr = DMoFEMLoopFiniteElements(dm_sub_sub_elastic,"ELASTIC",&elastic.getLoopFeRhs()); CHKERRQ(ierr);
ierr = DMoFEMLoopFiniteElements(dm_sub_sub_elastic,"ELASTIC_PRISM",&fat_prism_rhs); CHKERRQ(ierr);
//elastic element matrix
elastic.getLoopFeLhs().snes_B = Kuu;
fat_prism_lhs.snes_B = Kuu;
ierr = DMoFEMLoopFiniteElements(dm_sub_sub_elastic,"ELASTIC",&elastic.getLoopFeLhs()); CHKERRQ(ierr);
ierr = DMoFEMLoopFiniteElements(dm_sub_sub_elastic,"ELASTIC_PRISM",&fat_prism_lhs); CHKERRQ(ierr);
//postproc
ierr = DMoFEMPostProcessFiniteElements(dm_sub_sub_elastic,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
ierr = VecGhostUpdateBegin(Fu,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(Fu,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = VecAssemblyBegin(Fu); CHKERRQ(ierr);
ierr = VecAssemblyEnd(Fu); CHKERRQ(ierr);
ierr = VecScale(Fu,-1); CHKERRQ(ierr);
ierr = VecDestroy(&Du); CHKERRQ(ierr);
ierr = VecDestroy(&Fu); CHKERRQ(ierr);
const Problem *prb_ptr;
ierr = m_field.get_problem("ELASTIC_PROB",&prb_ptr); CHKERRQ(ierr);
boost::shared_ptr<Problem::SubProblemData> sub_data = prb_ptr->getSubData();
ierr = sub_data->getRowIs(&sub_nested_is_rows[0]); CHKERRQ(ierr);
ierr = sub_data->getColIs(&sub_nested_is_cols[0]); CHKERRQ(ierr);
sub_nested_matrices(0,0) = Kuu;
IS isUpsilon;
ierr = m_field.query_interface<ISManager>()->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;
}
ierr = DMDestroy(&dm_sub_sub_elastic); CHKERRQ(ierr);
}
{
DM dm_sub_disp_penalty;
// Craete dm_control instance
ierr = DMCreate(PETSC_COMM_WORLD,&dm_sub_disp_penalty);CHKERRQ(ierr);
ierr = DMSetType(dm_sub_disp_penalty,dm_name);CHKERRQ(ierr);
//set dm_sub_disp_penalty data structure which created mofem datastructures
ierr = DMMoFEMCreateSubDM(dm_sub_disp_penalty,dm_sub_volume_control,"S_PROB"); CHKERRQ(ierr);
ierr = DMMoFEMSetSquareProblem(dm_sub_disp_penalty,PETSC_FALSE); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldRow(dm_sub_disp_penalty,"UPSILON"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_disp_penalty,"U"); CHKERRQ(ierr);
//add elements to dm_sub_disp_penalty
ierr = DMMoFEMAddElement(dm_sub_disp_penalty,"DISPLACEMENTS_PENALTY"); CHKERRQ(ierr);
ierr = DMSetUp(dm_sub_disp_penalty); CHKERRQ(ierr);
Mat S;
ierr = DMCreateMatrix(dm_sub_disp_penalty,&S); CHKERRQ(ierr);
ierr = MatZeroEntries(S); CHKERRQ(ierr);
CellEngineering::FaceElement face_element(m_field);
face_element.getOpPtrVector().push_back(new OpCellS(S,eps_u));
ierr = DMoFEMLoopFiniteElements(dm_sub_disp_penalty,"DISPLACEMENTS_PENALTY",&face_element); CHKERRQ(ierr);
ierr = MatAssemblyBegin(S,MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(S,MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
// // Matrix View
if(debug) {
cerr << "S" << endl;
MatView(S,PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
const Problem *problem_ptr;
ierr = m_field.get_problem("S_PROB",&problem_ptr); CHKERRQ(ierr);
boost::shared_ptr<Problem::SubProblemData> sub_data = problem_ptr->getSubData();
// ierr = sub_data->getRowIs(&sub_nested_is_rows[1]); CHKERRQ(ierr);
// ierr = sub_data->getColIs(&sub_nested_is_cols[0]); CHKERRQ(ierr);
sub_nested_matrices(1,0) = S;
ierr = DMDestroy(&dm_sub_disp_penalty); CHKERRQ(ierr);
}
// Calculate penalty matrix
{
DM dm_sub_force_penalty;
// Craete dm_control instance
ierr = DMCreate(PETSC_COMM_WORLD,&dm_sub_force_penalty);CHKERRQ(ierr);
ierr = DMSetType(dm_sub_force_penalty,dm_name);CHKERRQ(ierr);
//set dm_sub_force_penalty data structure which created mofem datastructures
ierr = DMMoFEMCreateSubDM(dm_sub_force_penalty,dm_control,"D_PROB"); CHKERRQ(ierr);
ierr = DMMoFEMSetSquareProblem(dm_sub_force_penalty,PETSC_TRUE); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldRow(dm_sub_force_penalty,"RHO"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_force_penalty,"RHO"); CHKERRQ(ierr);
//add elements to dm_sub_force_penalty
ierr = DMMoFEMAddElement(dm_sub_force_penalty,"D"); CHKERRQ(ierr);
ierr = DMSetUp(dm_sub_force_penalty); CHKERRQ(ierr);
Mat D;
ierr = DMCreateMatrix(dm_sub_force_penalty,&D); CHKERRQ(ierr);
ierr = MatZeroEntries(D); CHKERRQ(ierr);
{
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));
}
ierr = DMoFEMLoopFiniteElements(dm_sub_force_penalty,"D",&face_d_matrix); CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
const Problem *problem_ptr;
ierr = m_field.get_problem("D_PROB",&problem_ptr); CHKERRQ(ierr);
// 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;
ierr = problem_ptr->getDofByNameEntAndEntDofIdx(
"RHO",vertex_to_fix[0],0,ROW,dof_ptr
); CHKERRQ(ierr);
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;
}
}
}
ierr = MatZeroRowsColumns(
D,nb_dofs_to_fix,&index_to_fix,eps_rho/eps_u,PETSC_NULL,PETSC_NULL
); CHKERRQ(ierr);
} else {
std::vector<int> dofs_to_fix;
for(Range::iterator eit = edges_to_fix.begin();eit!=edges_to_fix.end();eit++) {
for(_IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(problem_ptr,"RHO",*eit,m_field.get_comm_rank(),dof_it)) {
dofs_to_fix.push_back(dof_it->get()->getPetscGlobalDofIdx());
}
}
ierr = MatZeroRowsColumns(
D,dofs_to_fix.size(),&*dofs_to_fix.begin(),eps_rho/eps_u,PETSC_NULL,PETSC_NULL
); CHKERRQ(ierr);
}
// 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();
ierr = sub_data->getRowIs(&nested_is_rows[1]); CHKERRQ(ierr);
ierr = sub_data->getColIs(&nested_is_cols[1]); CHKERRQ(ierr);
nested_matrices(1,1) = D;
ierr = DMDestroy(&dm_sub_force_penalty); CHKERRQ(ierr);
}
{
DM dm_sub_force;
// Craete dm_control instance
ierr = DMCreate(PETSC_COMM_WORLD,&dm_sub_force);CHKERRQ(ierr);
ierr = DMSetType(dm_sub_force,dm_name);CHKERRQ(ierr);
//set dm_sub_force data structure which created mofem data structures
ierr = DMMoFEMCreateSubDM(dm_sub_force,dm_control,"FORCES_PROB"); CHKERRQ(ierr);
ierr = DMMoFEMSetSquareProblem(dm_sub_force,PETSC_FALSE); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldRow(dm_sub_force,"RHO"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_force,"U"); CHKERRQ(ierr);
ierr = DMMoFEMAddSubFieldCol(dm_sub_force,"UPSILON"); CHKERRQ(ierr);
//add elements to dm_sub_force
ierr = DMMoFEMAddElement(dm_sub_force,"B"); CHKERRQ(ierr);
ierr = DMSetUp(dm_sub_force); CHKERRQ(ierr);
Mat UB,UPSILONB;
ierr = DMCreateMatrix(dm_sub_force,&UB); CHKERRQ(ierr);
ierr = MatZeroEntries(UB); CHKERRQ(ierr);
// note be will be transposed in place latter on
ierr = DMCreateMatrix(dm_sub_force,&UPSILONB); CHKERRQ(ierr);
ierr = MatZeroEntries(UPSILONB); CHKERRQ(ierr);
{
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"));
}
ierr = DMoFEMLoopFiniteElements(dm_sub_force,"B",&face_b_matrices); CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(UB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyBegin(UPSILONB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(UB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(UPSILONB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
const Problem *problem_ptr;
ierr = m_field.get_problem("FORCES_PROB",&problem_ptr); CHKERRQ(ierr);
// 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;
ierr = problem_ptr->getDofByNameEntAndEntDofIdx(
"RHO",vertex_to_fix[0],0,ROW,dof_ptr
); CHKERRQ(ierr);
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;
}
}
}
ierr = MatZeroRows(UB,nb_dofs_to_fix,&index_to_fix,0,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
ierr = MatZeroRows(UPSILONB,nb_dofs_to_fix,&index_to_fix,0,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
} else {
std::vector<int> dofs_to_fix;
for(Range::iterator eit = edges_to_fix.begin();eit!=edges_to_fix.end();eit++) {
for(_IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(problem_ptr,"RHO",*eit,m_field.get_comm_rank(),dof_it)) {
dofs_to_fix.push_back(dof_it->get()->getPetscGlobalDofIdx());
}
}
ierr = MatZeroRows(UB,dofs_to_fix.size(),&*dofs_to_fix.begin(),0,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
ierr = MatZeroRows(UPSILONB,dofs_to_fix.size(),&*dofs_to_fix.begin(),0,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
}
Mat UBT;
ierr = MatTranspose(UB,MAT_INITIAL_MATRIX,&UBT); CHKERRQ(ierr);
ierr = MatDestroy(&UB); CHKERRQ(ierr);
// 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();
// ierr = sub_data->getColIs(&nested_is_rows[0]); CHKERRQ(ierr);
// ierr = sub_data->getRowIs(&nested_is_cols[1]); CHKERRQ(ierr);
nested_matrices(0,1) = UBT;
if(debug) {
cerr << "UPSILONB" << endl;
MatView(UPSILONB,PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
// ierr = sub_data->getRowIs(&nested_is_rows[1]); CHKERRQ(ierr);
// ierr = sub_data->getColIs(&nested_is_cols[0]); CHKERRQ(ierr);
nested_matrices(1,0) = UPSILONB;
ierr = DMDestroy(&dm_sub_force); CHKERRQ(ierr);
}
Mat SubA;
ierr = MatCreateNest(
PETSC_COMM_WORLD,
2,&sub_nested_is_rows[0],
2,&sub_nested_is_cols[0],
&sub_nested_matrices(0,0),
&SubA
); CHKERRQ(ierr);
nested_matrices(0,0) = SubA;
ierr = MatAssemblyBegin(SubA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(SubA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
if(debug) {
cerr << "Nested SubA" << endl;
MatView(SubA,PETSC_VIEWER_STDOUT_WORLD);
}
Mat A;
ierr = MatCreateNest(
PETSC_COMM_WORLD,
2,&nested_is_rows[0],
2,&nested_is_cols[0],
&nested_matrices(0,0),
&A
); CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
if(debug) {
cerr << "Nested A" << endl;
MatView(A,PETSC_VIEWER_STDOUT_WORLD);
}
Vec D,F;
ierr = DMCreateGlobalVector(dm_control,&D); CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm_control,&F); CHKERRQ(ierr);
// 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));
dm_control,"DISPLACEMENTS_PENALTY",&face_element
); CHKERRQ(ierr);
ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
}
KSP solver;
// KSP solver;
{
ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
ierr = KSPSetDM(solver,dm_control); CHKERRQ(ierr);
ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
ierr = KSPSetDMActive(solver,PETSC_FALSE); CHKERRQ(ierr);
ierr = KSPSetInitialGuessKnoll(solver,PETSC_FALSE); CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(solver,PETSC_FALSE); CHKERRQ(ierr);
PC pc;
ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
ierr = PCSetType(pc,PCFIELDSPLIT); CHKERRQ(ierr);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);
if(is_pcfs) {
ierr = PCSetOperators(pc,A,A); CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc,NULL,nested_is_rows[0]); CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc,NULL,nested_is_rows[1]); CHKERRQ(ierr);
ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR); CHKERRQ(ierr);
ierr = PCSetUp(pc); CHKERRQ(ierr);
KSP *sub_ksp;
PetscInt n;
ierr = PCFieldSplitGetSubKSP(pc,&n,&sub_ksp); CHKERRQ(ierr);
{
PC sub_pc_0;
ierr = KSPGetPC(sub_ksp[0],&sub_pc_0); CHKERRQ(ierr);
ierr = PCSetOperators(sub_pc_0,SubA,SubA); CHKERRQ(ierr);
ierr = PCSetType(sub_pc_0,PCFIELDSPLIT); CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(sub_pc_0,NULL,sub_nested_is_rows[0]); CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(sub_pc_0,NULL,sub_nested_is_rows[1]); CHKERRQ(ierr);
ierr = PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_MULTIPLICATIVE); CHKERRQ(ierr);
// ierr = PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER); CHKERRQ(ierr);
// ierr = PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR); CHKERRQ(ierr);
ierr = PCSetUp(sub_pc_0); CHKERRQ(ierr);
}
} else {
SETERRQ(
PETSC_COMM_WORLD,
"Only works with pre-conditioner PCFIELDSPLIT"
);
}
ierr = KSPSetUp(solver); CHKERRQ(ierr);
}
// Solve system of equations
ierr = KSPSolve(solver,F,D); CHKERRQ(ierr);
ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = DMoFEMMeshToLocalVector(dm_control,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
if(debug) {
ierr = VecView(D,PETSC_VIEWER_DRAW_WORLD); CHKERRQ(ierr);
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]) {
ierr = ISDestroy(&sub_nested_is_rows[i]); CHKERRQ(ierr);
}
if(sub_nested_is_cols[i]) {
ierr = ISDestroy(&sub_nested_is_cols[i]); CHKERRQ(ierr);
}
for(int j = 0;j!=2;j++) {
if(sub_nested_matrices(i,j)) {
ierr = MatDestroy(&sub_nested_matrices(i,j)); CHKERRQ(ierr);
}
}
}
for(int i = 0;i!=2;i++) {
if(nested_is_rows[i]) {
ierr = ISDestroy(&nested_is_rows[i]); CHKERRQ(ierr);
}
if(nested_is_cols[i]) {
ierr = ISDestroy(&nested_is_cols[i]); CHKERRQ(ierr);
}
for(int j = 0;j!=2;j++) {
if(nested_matrices(i,j)) {
ierr = MatDestroy(&nested_matrices(i,j)); CHKERRQ(ierr);
}
}
}
ierr = MatDestroy(&SubA); CHKERRQ(ierr);
ierr = MatDestroy(&A); CHKERRQ(ierr);
ierr = VecDestroy(&D); CHKERRQ(ierr);
ierr = VecDestroy(&F); CHKERRQ(ierr);
ierr = DMDestroy(&dm_sub_volume_control); CHKERRQ(ierr);
PostProcVolumeOnRefinedMesh post_proc(m_field);
{
ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
ierr = post_proc.addFieldValuesPostProc("U"); CHKERRQ(ierr);
ierr = post_proc.addFieldValuesGradientPostProc("U"); CHKERRQ(ierr);
//add postpocessing for sresses
post_proc.getOpPtrVector().push_back(
m_field,
post_proc.postProcMesh,
post_proc.mapGaussPts,
"U",
post_proc.commonData,
&elastic.setOfBlocks
)
);
ierr = DMoFEMLoopFiniteElements(dm_control,"ELASTIC",&post_proc); CHKERRQ(ierr);
ierr = post_proc.writeFile("out.h5m"); CHKERRQ(ierr);
elastic.getLoopFeEnergy().eNergy = 0;
ierr = DMoFEMLoopFiniteElements(dm_control,"ELASTIC",&elastic.getLoopFeEnergy()); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"Elastic energy %6.4e\n",elastic.getLoopFeEnergy().eNergy);
}
{
PostProcFaceOnRefinedMesh post_proc_face(m_field);
ierr = post_proc_face.generateReferenceElementMesh(); CHKERRQ(ierr);
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(
m_field,
post_proc_face.postProcMesh,
post_proc_face.mapGaussPts,
common_data
)
);
ierr = DMoFEMLoopFiniteElements(dm_control,"D",&post_proc_face); CHKERRQ(ierr);
ierr = post_proc_face.writeFile("out_tractions.h5m"); CHKERRQ(ierr);
}
ierr = DMDestroy(&dm_control); CHKERRQ(ierr);
} catch (MoFEMException const &e) {
SETERRQ(PETSC_COMM_SELF,e.errorCode,e.errorMessage);
}
PetscFinalize();
return 0;
}