cell_forces.cpp

Complete example for identifying cell forces from measured displacements.

See Force-potential implementation for the mathematical formulation, implementation details, and command-line usage.

/**
  \file cell_forces.cpp
  \brief Identify cell tractions from measured displacements
 
  \tableofcontents
 
  \section cell_potential Force-potential implementation
 
This example identifies surface tractions that produce measured displacements
on an interface within an elastic body. In the experiment, a transparent gel
layer covers the body, and the motion of marker beads at the interface is
measured. Cells on the upper surface deform the gel. The objective is to infer
the tractions exerted by the cells from the measured displacements.
 
The current implementation assumes that the measurement and traction surfaces
are planar and parallel to the \f$x\f$-\f$y\f$ plane. An extension to curved
surfaces would require intrinsic surface differential operators and a
well-defined model for any normal traction.
 
\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_formulation Inverse traction formulation
 
Let \f$V\f$ be the elastic body, \f$S_u\f$ the measurement surface, and
\f$S_\rho\f$ the surface on which the unknown traction acts. The state
displacement is denoted by \f$\mathbf{u}\f$, the measured displacement by
\f$\mathbf{u}_d\f$, the Lagrange multiplier (adjoint field) by
\f$\Upsilon\f$, and the traction by \f$\mathbf{t}\f$. The small-strain
constitutive equation is
\f[
\sigma(\mathbf{u}) =
\mathcal{A}:\varepsilon(\mathbf{u}), \qquad
\varepsilon(\mathbf{u}) =
\frac{1}{2}\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm T}\right).
\f]
For a symmetric elasticity tensor \f$\mathcal{A}\f$, define
\f[
a(\mathbf{u},\mathbf{v}) =
\left(\varepsilon(\mathbf{v}),
\mathcal{A}:\varepsilon(\mathbf{u})\right)_V.
\f]
 
The inverse problem is
\f[
\begin{aligned}
\underset{\mathbf{u},\mathbf{t}}{\operatorname{minimise}}\quad
J(\mathbf{u},\mathbf{t}) &=
\frac{1}{2}
\left(\mathbf{u}-\mathbf{u}_d,
S(\mathbf{u}-\mathbf{u}_d)\right)_{S_u}
+R(\mathbf{t}),\\
-\operatorname{div}\sigma(\mathbf{u})&=0
\quad\text{in }V,\\
\sigma(\mathbf{u})\mathbf{n}&=\mathbf{t}
\quad\text{on }S_\rho .
\end{aligned}
\f]
The remaining boundary conditions must make the elasticity problem unique. In
particular, sufficient displacement constraints are required to eliminate
rigid-body modes.
 
After integration by parts, the weak equilibrium constraint is
\f[
a(\mathbf{u},\mathbf{v})-(\mathbf{t},\mathbf{v})_{S_\rho}=0.
\f]
It follows from
\f$-(\Upsilon,\operatorname{div}\sigma)_V\f$; boundary
terms on portions other than \f$S_\rho\f$ vanish because of the corresponding
essential or natural boundary conditions. The Lagrangian is therefore
\f[
\mathcal{L}(\mathbf{u},\mathbf{t},\Upsilon)=
\frac{1}{2}
\left(\mathbf{u}-\mathbf{u}_d,
S(\mathbf{u}-\mathbf{u}_d)\right)_{S_u}
+R(\mathbf{t})
+a(\mathbf{u},\Upsilon)
-(\mathbf{t},\Upsilon)_{S_\rho}.
\f]
Thus, traction is coupled to equilibrium in the weak sense.
 
For the quadratic regularisation
\f[
R(\mathbf{t})=
\frac{\epsilon_\rho^\mathrm{abs}}{2}
(\mathbf{t},\mathbf{t})_{S_\rho},
\f]
stationarity gives
\f[
\begin{aligned}
\left(S(\mathbf{u}-\mathbf{u}_d),\delta\mathbf{u}\right)_{S_u}
+a(\delta\mathbf{u},\Upsilon)&=0,\\
a(\mathbf{u},\delta\Upsilon)
-(\mathbf{t},\delta\Upsilon)_{S_\rho}&=0,\\
\epsilon_\rho^\mathrm{abs}
(\mathbf{t},\delta\mathbf{t})_{S_\rho}
-(\Upsilon,\delta\mathbf{t})_{S_\rho}&=0.
\end{aligned}
\f]
 
Finite element discretisation gives
\f[
\left[
\begin{array}{ccc}
S & K & 0 \\
K & 0 & -B^\mathrm{T} \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
\mathbf{u} \\ \Upsilon \\ \mathbf{t}
\end{array}
\right]
=
\left[
\begin{array}{c}
S\mathbf{u}_d \\ 0 \\ 0
\end{array}
\right],
\f]
where the two elastic blocks are identical because the linear-elastic bilinear
form is symmetric. Swapping the first two block rows produces the system used
by the implementation:
\f[
\left[
\begin{array}{ccc}
K & 0 & -B^\mathrm{T} \\
S & K & 0 \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
\mathbf{u} \\ \Upsilon \\ \mathbf{t}
\end{array}
\right]
=
\left[
\begin{array}{c}
0 \\ S\mathbf{u}_d \\ 0
\end{array}
\right].
\f]
 
Only the measured \f$x\f$- and \f$y\f$-components enter the displacement
misfit. Consequently,
\f[
S=\epsilon_u^{-1}P_{xy},\qquad
P_{xy}=\operatorname{diag}(1,1,0),\qquad
\epsilon_\rho^\mathrm{abs}
=\epsilon_\rho^\mathrm{ref}\epsilon_u^{-1}.
\f]
 
\subsection cell_local Exact curl-free force-potential formulation
 
The local model assumes that the tangential traction is curl-free. This is a
modelling assumption motivated by an idealised network of short, prestressed
fibres; it is not a consequence of equilibrium alone. On a simply connected
planar surface, the traction can then be represented by a scalar potential:
\f[
\mathbf{t}=\nabla_s\Phi ,
\f]
where \f$\nabla_s\f$ denotes the surface gradient. The regularisation and
traction coupling become
\f[
R(\Phi)=\frac{\epsilon_\rho^\mathrm{abs}}{2}
(\nabla_s\Phi,\nabla_s\Phi)_{S_\rho},\qquad
-(\nabla_s\Phi,\Upsilon)_{S_\rho}.
\f]
Variation with respect to the potential gives
\f[
\epsilon_\rho^\mathrm{abs}
(\nabla_s\Phi,\nabla_s\delta\Phi)_{S_\rho}
-(\Upsilon,\nabla_s\delta\Phi)_{S_\rho}=0.
\f]
The third algebraic unknown in this formulation is \f$\Phi\f$, not
\f$\mathbf{t}\f$. The matrix \f$D\f$ is the surface stiffness matrix for the
potential, and one value of \f$\Phi\f$ must be fixed to remove its arbitrary
additive constant.
 
  \subsubsection cell_local_approx Approximation
 
The state and adjoint fields satisfy
\f$\mathbf{u},\Upsilon\in[H^1(V)]^3\f$, with the appropriate
essential boundary conditions. The potential satisfies
\f$\Phi\in H^1(S_\rho)\f$, and the recovered tangential traction
\f$\nabla_s\Phi\f$ belongs to \f$[L^2(S_\rho)]^2\f$.
 
  \subsection cell_nonlocal Weak curl-free formulation
 
The nonlocal model represents the traction directly and penalises, rather than
eliminates, its surface curl. It is intended to account for the finite size and
complex prestressed structure of a cell. The regularisation is
\f[
R(\mathbf{t})=
\frac{\epsilon_\rho^\mathrm{abs}}{2}
(\mathbf{t},\mathbf{t})_{S_\rho}
+\frac{1}{2\epsilon_l^\mathrm{abs}}
(\operatorname{curl}_s\mathbf{t},
\operatorname{curl}_s\mathbf{t})_{S_\rho}.
\f]
The traction optimality equation is
\f[
\epsilon_\rho^\mathrm{abs}
(\mathbf{t},\delta\mathbf{t})_{S_\rho}
+(\epsilon_l^\mathrm{abs})^{-1}
(\operatorname{curl}_s\mathbf{t},
\operatorname{curl}_s\delta\mathbf{t})_{S_\rho}
-(\Upsilon,\delta\mathbf{t})_{S_\rho}=0.
\f]
The command-line parameters are scaled in the implementation as
\f[
\epsilon_\rho^\mathrm{abs}
=\epsilon_\rho^\mathrm{ref}\epsilon_u^{-1},\qquad
\epsilon_l^\mathrm{abs}
=\epsilon_l^\mathrm{ref}\epsilon_u.
\f]
For positive \f$\epsilon_l^\mathrm{ref}\f$, decreasing it strengthens the
curl-free penalty. The exact value \f$\epsilon_l^\mathrm{ref}=0\f$ disables the
curl term in the current implementation to avoid division by zero; it does not
select the limiting exact-potential problem. Use \c -my_curl \c 0 to select the
exact force-potential formulation.
 
  \subsubsection cell_nonlocal_approx Approximation
 
The state and adjoint fields satisfy
\f$\mathbf{u},\Upsilon\in[H^1(V)]^3\f$. The traction belongs to the
surface space
\f$\mathbf{t}\in H(\operatorname{curl}_s;S_\rho)\f$.
 
\section cell_solution Field-split preconditioner
 
The implementation uses the block structure of the matrix through
\e PCFIELDSPLIT and stores the matrix in the nested \e MATNEST format.
 
The nested matrix containing the two upper-left blocks is
\f[
X = \left[
\begin{array}{cc}
K & 0  \\
S & K
\end{array}
\right]
\f]
The complete matrix is then written as
\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
\end{array}
\right]
\f]
 
The system associated with \f$X\f$ is solved using \e PCFIELDSPLIT and a
multiplicative relaxation scheme. The submatrix \f$K\f$ is factorised only
once. The system associated with \f$A\f$, composed of
\f$X\f$, \f$D\f$, \f$V\f$, and \f$H\f$, is solved using a Schur complement.
MoFEM subproblems are used to construct the required blocks.
 
\section cell_running_code Running code
 
Mapping measured 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, a second thin polymer layer can be created
automatically using prism elements:
\code
./map_disp_prism \
-my_thickness 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
 
Configuration file:
\include users_modules/cell_engineering/examples/block_config.in
 
Calculating tractions:
\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
 
The resulting traction field is shown below:
\image html cell_engineering_forces_example.gif "Example: results" width=600px
 
\section cell_install Installation with Docker
 
- First, install Docker according to the instructions at:
<https://docs.docker.com/installation/#installation>
- Install the cell user module: \e docker \e pull \e likask/cell_engineering
 
\todo Improve documentation
 
\todo Generalise the implementation to curved surfaces using intrinsic surface
differential operators and an appropriate model for normal traction.
 
\todo Replace the elastic material with the GEL model developed in another
module to account for drying and other rheological effects.
 
\todo Use an SNES solver when introducing a nonlinear material model.
 
\todo Add equations for cell mechanics either directly to the minimised
functional or as additional constraints.
 
*/
 
/**
  \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 traction approximation", "",
        order_force, &order_force, &flg_order_force);
    CHKERR PetscOptionsString("-my_block_config",
                              "elastic configuration file name", "",
                              "block_conf.in", block_config_file, 255,
                              &flg_block_config);
    CHKERR PetscOptionsReal("-my_eps_rho", "traction regularisation parameter",
                            "", eps_rho, &eps_rho, &flg_eps_rho);
    CHKERR PetscOptionsReal(
        "-my_eps_u", "displacement-misfit regularisation parameter", "", eps_u,
        &eps_u, &flg_eps_u);
    CHKERR PetscOptionsReal("-my_eps_l", "curl regularisation parameter", "",
                            eps_l, &eps_l, &flg_eps_l);
    CHKERR PetscOptionsBool(
        "-my_curl", "use H(curl) space to approximate tractions", "", is_curl,
        &is_curl, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    // Read command-line parameters
    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 the approximation 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 entities (by tetrahedra) 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;
    // In a partitioned problem, the meshset might not exist on this rank
    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;
    // In a partitioned problem, the meshset might not exist on this rank
    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 elements
    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;
    {
      // Create the dm_control instance
      CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
      CHKERR DMSetType(dm_control, dm_name);
      // Set the dm_control data structure that created the MoFEM data structures
      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;
 
      // Create the 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
      // data structures
      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;
 
      // Create the 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
      // data structures
      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;
 
      // Create the 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);
      // This matrix will be transposed later
      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);
 
    // Assemble 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,
                "This solver requires the PCFIELDSPLIT preconditioner");
      }
      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;
}