cell_forces.cpp

Complete example for identifying cell forces from measured displacements.

See Implementation with forces potential for the mathematical formulation, implementation details, and command-line usage.

/**
  \file cell_forces.cpp
  \brief Calculate cell forces
 
  \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.
 
Noting 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. Utilizing 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 generalization of the local model to a weakly enforced force
curl-free cell force field is considered. This generalization 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.
 
*/
 
/**
  \example cell_forces.cpp
  Complete example for identifying cell forces from measured displacements.
 
  See \ref cell_potential for the mathematical formulation, implementation
  details, and command-line usage.
*/
 
/* Copyright (c) 2025 Authors under the MIT License.
 * See LICENSE.md for details.
 * SPDX-License-Identifier: MIT
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <CellForces.hpp>
 
constexpr int SPACE_DIM = 3;
using DomainEleOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
 
#include <HookeOps.hpp>
 
static char help[] = "-my_block_config set block data\n"
                     "\n";
 
namespace CellEngineering {
 
struct BlockOptionData {
  int oRder;
  double yOung;
  double pOisson;
  BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2) {}
};
 
struct ElasticVolume : public VolumeElementForcesAndSourcesCore {
  using VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore;
  int getRule(int order) override { return 2 * order + 1; }
};
 
using ElasticBlockMap = std::map<int, Range>;
 
struct OpElasticEnergy
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  OpElasticEnergy(boost::shared_ptr<MatrixDouble> strain_ptr,
                  boost::shared_ptr<MatrixDouble> stress_ptr, double &energy)
      : VolumeElementForcesAndSourcesCore::UserDataOperator(
            NOSPACE, UserDataOperator::OPSPACE),
        strainPtr(strain_ptr), stressPtr(stress_ptr), energy(energy) {}
 
  MoFEMErrorCode doWork(int, EntityType,
                        EntitiesFieldData::EntData &) override {
    MoFEMFunctionBegin;
 
    const size_t nb_gauss_pts = getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto get_strain =
        MatrixSizeHelper<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>,
                         DL>::get(*strainPtr, nb_gauss_pts);
    auto get_stress =
        MatrixSizeHelper<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>,
                         DL>::get(*stressPtr, nb_gauss_pts);
    auto t_strain = get_strain();
    auto t_stress = get_stress();
    FTensor::Index<'i', SPACE_DIM> i;
    FTensor::Index<'j', SPACE_DIM> j;
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      const double alpha = getMeasure() * getGaussPts()(SPACE_DIM, gg);
      energy += 0.5 * alpha * t_strain(i, j) * t_stress(i, j);
      ++t_strain;
      ++t_stress;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> strainPtr;
  boost::shared_ptr<MatrixDouble> stressPtr;
  double &energy;
};
 
} // 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;
    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_NULLPTR);
    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, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    // 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 == nullptr)
      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;
 
    auto mmanager_ptr = m_field.getInterface<MeshsetsManager>();
    auto comm_interface_ptr = m_field.getInterface<CommInterface>();
    // 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.getInterface<BitRefManager>()->setBitRefLevel(
          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));
        for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
          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);
        for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
          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 comm_interface_ptr->synchroniseEntities(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,
                             MF_ZERO);
    CHKERR m_field.add_field("UPSILON", H1, AINSWORTH_LEGENDRE_BASE, 3,
                             MB_TAG_SPARSE, MF_ZERO);
    if (is_curl) {
      CHKERR m_field.add_field("RHO", HCURL, AINSWORTH_LEGENDRE_BASE, 1);
    } else {
      CHKERR m_field.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1);
    }
 
    // 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 comm_interface_ptr->synchroniseFieldEntities("U");
    CHKERR comm_interface_ptr->synchroniseFieldEntities("UPSILON");
    CHKERR comm_interface_ptr->synchroniseFieldEntities("RHO");
 
    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()) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Should be one vertex only, but is %zu", vertex_to_fix.size());
      }
    }
    CHKERR comm_interface_ptr->synchroniseEntities(ents_1st_layer, 0);
    CHKERR m_field.add_ents_to_field_by_type(
        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 comm_interface_ptr->synchroniseEntities(ents_2nd_layer, 0);
 
    for (int oo = 2; oo != setOrderToEnts.size(); oo++) {
      if (setOrderToEnts[oo].size() > 0) {
        CHKERR comm_interface_ptr->synchroniseEntities(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 comm_interface_ptr->synchroniseEntities(ents, 0);
      CHKERR comm_interface_ptr->synchroniseEntities(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);
    }
 
    // The vec-0 Hooke operators read material data directly from
    // MAT_ELASTIC* meshsets.
    ElasticBlockMap elastic_blocks;
    int default_block_id = -1;
    Tag block_id_tag;
    CHKERR moab.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, block_id_tag,
                               MB_TAG_CREAT | MB_TAG_SPARSE, &default_block_id);
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, BLOCKSET | MAT_ELASTICSET, it)) {
      Mat_Elastic material;
      CHKERR it->getAttributeDataStructure(material);
      const int block_id = it->getMeshsetId();
      CHKERR moab.get_entities_by_handle(it->getMeshset(),
                                         elastic_blocks[block_id], true);
      const auto block_elements =
          elastic_blocks[block_id].subset_by_type(MBTET);
      CHKERR moab.tag_clear_data(block_id_tag, block_elements, &block_id);
 
      const auto options_it = block_data.find(block_id);
      if (options_it != block_data.end() && options_it->second.yOung > 0) {
        material.data.Young = options_it->second.yOung;
        CHKERR PetscPrintf(PETSC_COMM_WORLD,
                           "Block %d set Young modulus %3.4g\n", block_id,
                           material.data.Young);
      }
      if (options_it != block_data.end() && options_it->second.pOisson >= -1) {
        material.data.Poisson = options_it->second.pOisson;
        CHKERR PetscPrintf(PETSC_COMM_WORLD,
                           "Block %d set Poisson ratio %3.4g\n", block_id,
                           material.data.Poisson);
      }
      CHKERR mmanager_ptr->setAttributesByDataStructure(BLOCKSET, block_id,
                                                        material);
    }
 
    CHKERR m_field.add_finite_element("ELASTIC", MF_ZERO);
    CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "U");
    CHKERR m_field.modify_finite_element_add_field_col("ELASTIC", "U");
    CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "U");
    for (const auto &[id, entities] : elastic_blocks)
      CHKERR m_field.add_ents_to_finite_element_by_type(entities, MBTET,
                                                        "ELASTIC");
 
    ElasticVolume elastic_rhs(m_field);
    ElasticVolume elastic_lhs(m_field);
    ElasticVolume elastic_energy_fe(m_field);
    double elastic_energy = 0;
    auto elastic_rhs_common =
        HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
            m_field, elastic_rhs.getOpPtrVector(), "U", "MAT_ELASTIC",
            Sev::verbose);
    CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        m_field, elastic_rhs.getOpPtrVector(), "U", elastic_rhs_common,
        Sev::verbose, true);
    auto elastic_lhs_common =
        HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
            m_field, elastic_lhs.getOpPtrVector(), "U", "MAT_ELASTIC",
            Sev::verbose);
    CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        m_field, elastic_lhs.getOpPtrVector(), "U", elastic_lhs_common,
        Sev::verbose);
    auto elastic_energy_common =
        HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
            m_field, elastic_energy_fe.getOpPtrVector(), "U", "MAT_ELASTIC",
            Sev::verbose);
    elastic_energy_fe.getOpPtrVector().push_back(new OpElasticEnergy(
        elastic_energy_common->getMatStrain(),
        elastic_energy_common->getMatCauchyStress(), elastic_energy));
 
    // 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");
      auto block_it = elastic_blocks.find(2);
      if (block_it == elastic_blocks.end()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Elastic material block 2 is required for prism elements");
      }
      CHKERR mmanager_ptr->getEntitiesByDimension(2, BLOCKSET, 3,
                                                  block_it->second);
      CHKERR m_field.add_ents_to_finite_element_by_type(
          block_it->second, MBPRISM, "ELASTIC_PRISM");
      // right hand side operators
      auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
      fat_prism_rhs.getOpPtrVector().push_back(
          new MoFEM::OpCalculateInvJacForFatPrism(inv_jac_ptr));
      fat_prism_rhs.getOpPtrVector().push_back(
          new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac_ptr));
      using FatPrismOp = FatPrismElementForcesAndSourcesCore::UserDataOperator;
      auto prism_rhs_common =
          HookeOps::commonDataFactory<SPACE_DIM, GAUSS, FatPrismOp>(
              m_field, fat_prism_rhs.getOpPtrVector(), "U", "MAT_ELASTIC",
              Sev::verbose);
      CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, FatPrismOp>(
          m_field, fat_prism_rhs.getOpPtrVector(), "U", prism_rhs_common,
          Sev::verbose, true);
      // Left hand side operators
      fat_prism_lhs.getOpPtrVector().push_back(
          new MoFEM::OpCalculateInvJacForFatPrism(inv_jac_ptr));
      fat_prism_lhs.getOpPtrVector().push_back(
          new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac_ptr));
      auto prism_lhs_common =
          HookeOps::commonDataFactory<SPACE_DIM, GAUSS, FatPrismOp>(
              m_field, fat_prism_lhs.getOpPtrVector(), "U", "MAT_ELASTIC",
              Sev::verbose);
      CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, FatPrismOp>(
          m_field, fat_prism_lhs.getOpPtrVector(), "U", prism_lhs_common,
          Sev::verbose);
    }
 
    // 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", "U");
    CHKERR m_field.modify_finite_element_add_field_row("BT", "UPSILON");
    CHKERR m_field.modify_finite_element_add_field_col("BT", "RHO");
    CHKERR m_field.modify_finite_element_add_field_data("BT", "U");
    CHKERR m_field.modify_finite_element_add_field_data("BT", "UPSILON");
    CHKERR m_field.modify_finite_element_add_field_data("BT", "RHO");
    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.modify_finite_element_add_field_row("B", "RHO");
    CHKERR m_field.modify_finite_element_add_field_col("B", "U");
    CHKERR m_field.modify_finite_element_add_field_col("B", "UPSILON");
    CHKERR m_field.modify_finite_element_add_field_data("B", "U");
    CHKERR m_field.modify_finite_element_add_field_data("B", "UPSILON");
    CHKERR m_field.modify_finite_element_add_field_data("B", "RHO");
    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.modify_finite_element_add_field_row("D", "RHO");
    CHKERR m_field.modify_finite_element_add_field_col("D", "RHO");
    CHKERR m_field.modify_finite_element_add_field_data("D", "RHO");
    CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
                                                      "D");
 
    // build finite elemnts
    CHKERR m_field.build_finite_elements();
    // build adjacencies
    CHKERR m_field.build_adjacencies(bit_level0);
 
    // Register MOFEM DM
    DMType dm_name = "MOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    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);
    }
 
    CellEngineering::CommonData common_data;
 
    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_NULLPTR;
      nested_is_cols[i] = PETSC_NULLPTR;
      for (int j = 0; j != 2; j++) {
        nested_matrices(i, j) = PETSC_NULLPTR;
      }
    }
 
    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_NULLPTR;
      sub_nested_is_cols[i] = PETSC_NULLPTR;
      for (int j = 0; j != 2; j++) {
        sub_nested_matrices(i, j) = PETSC_NULLPTR;
      }
    }
 
    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_NULLPTR;
    }
 
    {
      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);
      CHKERR m_field.getInterface<BcManager>()
          ->pushMarkDOFsOnEntities<DisplacementCubitBcData>("ELASTIC_PROB",
                                                            "U");
 
      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);
 
      // Apply displacement constraints with the current core BC machinery.
      auto dirichlet_bc_ptr = boost::make_shared<FEMethod>();
      dirichlet_bc_ptr->vecAssembleSwitch =
          boost::movelib::make_unique<bool>(false);
      dirichlet_bc_ptr->matAssembleSwitch =
          boost::movelib::make_unique<bool>(false);
      dirichlet_bc_ptr->preProcessHook =
          EssentialPreProc<DisplacementCubitBcData>(
              m_field, dirichlet_bc_ptr,
              std::vector<boost::shared_ptr<ScalingMethod>>{}, false);
      dirichlet_bc_ptr->postProcessHook = [&]() {
        MoFEMFunctionBegin;
        CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecAssemblyBegin(Fu);
        CHKERR VecAssemblyEnd(Fu);
        CHKERR MatAssemblyBegin(Kuu, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(Kuu, MAT_FINAL_ASSEMBLY);
        CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
            m_field, dirichlet_bc_ptr, 0., SmartPetscObj<Vec>(Fu, true))();
        CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
            m_field, dirichlet_bc_ptr, 1., SmartPetscObj<Mat>(Kuu, true))();
        MoFEMFunctionReturn(0);
      };
      dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
      dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
      // preproc
      CHKERR DMoFEMPreProcessFiniteElements(dm_sub_sub_elastic,
                                            dirichlet_bc_ptr.get());
      CHKERR DMoFEMMeshToLocalVector(dm_sub_sub_elastic, Du, INSERT_VALUES,
                                     SCATTER_REVERSE);
      // internal force vector (to take into account  Dirichlet boundary
      // conditions)
      elastic_rhs.snes_f = Fu;
      fat_prism_rhs.snes_f = Fu;
      CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
                                      &elastic_rhs);
      CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
                                      &fat_prism_rhs);
      // elastic element matrix
      elastic_lhs.snes_B = Kuu;
      fat_prism_lhs.snes_B = Kuu;
      CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
                                      &elastic_lhs);
      CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
                                      &fat_prism_lhs);
      // postproc
      CHKERR DMoFEMPostProcessFiniteElements(dm_sub_sub_elastic,
                                             dirichlet_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;
      CHKERR m_field.getInterface<ISManager>()
          ->isCreateFromProblemFieldToOtherProblemField(
              "ELASTIC_PROB", "U", ROW, "SUB_CONTROL_PROB", "UPSILON", ROW,
              PETSC_NULLPTR, &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);
      CHKERR AddHOOps<2, 2, 2>::add(face_element.getOpPtrVector(), {H1});
      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);
 
        if (is_curl) {
          CHKERR AddHOOps<2, 2, 2>::add(face_d_matrix.getOpPtrVector(),
                                        {HCURL});
          face_d_matrix.getOpPtrVector().push_back(
              new OpCellCurlD(D, eps_rho / eps_u, eps_l * eps_u));
        } else {
          CHKERR AddHOOps<2, 2, 2>::add(face_d_matrix.getOpPtrVector(), {H1});
          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;
          CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
              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_NULLPTR,
                                  PETSC_NULLPTR);
      } 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 =
              row_dofs->upper_bound(FieldEntity::getHiLocalEntityBitNumber(
                  bit_number, p_eit->second));
          for (; lo != hi; ++lo)
            if ((*lo)->getPart() == m_field.get_comm_rank())
              dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
        }
        CHKERR MatZeroRowsColumns(D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
                                  eps_rho / eps_u, PETSC_NULLPTR,
                                  PETSC_NULLPTR);
      }
 
      // 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);
        if (is_curl) {
          CHKERR AddHOOps<2, 2, 2>::add(face_b_matrices.getOpPtrVector(),
                                        {H1, HCURL});
          face_b_matrices.getOpPtrVector().push_back(new OpCellCurlB(UB, "U"));
          face_b_matrices.getOpPtrVector().push_back(
              new OpCellCurlB(UPSILONB, "UPSILON"));
        } else {
          CHKERR AddHOOps<2, 2, 2>::add(face_b_matrices.getOpPtrVector(),
                                        {H1});
          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;
          CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
              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_NULLPTR,
                           PETSC_NULLPTR);
        CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
                           PETSC_NULLPTR, PETSC_NULLPTR);
      } 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 =
              row_dofs->upper_bound(FieldEntity::getHiLocalEntityBitNumber(
                  bit_number, p_eit->second));
          for (; lo != hi; ++lo)
            if ((*lo)->getPart() == m_field.get_comm_rank())
              dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
        }
        CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
                           PETSC_NULLPTR, PETSC_NULLPTR);
        CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
                           0, PETSC_NULLPTR, PETSC_NULLPTR);
      }
 
      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);
      CHKERR AddHOOps<2, 2, 2>::add(face_element.getOpPtrVector(), {H1});
      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);
 
    using PostProcVolume =
        PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
    using OpPPVolume = OpPostProcMapInMoab<3, 3>;
 
    PostProcVolume post_proc(m_field, "my");
    {
      CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {H1});
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("U", u_ptr));
      auto hooke_common =
          HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
              m_field, post_proc.getOpPtrVector(), "U", "MAT_ELASTIC",
              Sev::verbose);
      post_proc.getOpPtrVector().push_back(new OpPPVolume(
          post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
          OpPPVolume::DataMapVec{}, OpPPVolume::DataMapMat{{"U", u_ptr}},
          OpPPVolume::DataMapMat{{"U_GRAD", hooke_common->matGradPtr}},
          OpPPVolume::DataMapMat{
              {"STRESS", hooke_common->getMatCauchyStress()}}));
      CHKERR post_proc.setTagsToTransfer({block_id_tag});
 
      CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC", &post_proc);
      CHKERR post_proc.writeFile("out.h5m");
      elastic_energy = 0;
      CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC",
                                      &elastic_energy_fe);
      double global_elastic_energy = 0;
      CHKERR MPI_Allreduce(&elastic_energy, &global_elastic_energy, 1,
                           MPI_DOUBLE, MPI_SUM, m_field.get_comm());
      PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
                  global_elastic_energy);
    }
 
    {
      using PostProcFace =
          PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>;
      PostProcFace post_proc_face(m_field, "my");
      if (is_curl) {
        CHKERR AddHOOps<2, 2, 2>::add(post_proc_face.getOpPtrVector(),
                                      {HCURL});
        post_proc_face.getOpPtrVector().push_back(
            new OpVirtualCurlRho("RHO", common_data));
      } else {
        CHKERR AddHOOps<2, 2, 2>::add(post_proc_face.getOpPtrVector(), {H1});
        post_proc_face.getOpPtrVector().push_back(
            new OpVirtualPotentialRho("RHO", common_data));
      }
      post_proc_face.getOpPtrVector().push_back(
          new PostProcTraction(m_field, post_proc_face.getPostProcMesh(),
                               post_proc_face.getMapGaussPts(), common_data));
      CHKERR DMoFEMLoopFiniteElements(dm_control, "D", &post_proc_face);
      CHKERR post_proc_face.writeFile("out_tractions.h5m");
    }
 
    CHKERR DMDestroy(&dm_control);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}