Identify cell tractions from measured displacements. More...

#include <MoFEM.hpp>
#include <CellForces.hpp>
#include <HookeOps.hpp>
#include <boost/program_options.hpp>

Classes
struct	CellEngineering::BlockOptionData

struct	CellEngineering::ElasticVolume

struct	CellEngineering::OpElasticEnergy

Namespaces
namespace	CellEngineering

Typedefs
using	DomainEleOp = VolumeElementForcesAndSourcesCore::UserDataOperator

using	CellEngineering::ElasticBlockMap = std::map< int, Range >

Functions
int	main (int argc, char *argv[])

Variables
constexpr int	SPACE_DIM = 3

static char	help []

static int	debug = 0

Detailed Description

Identify cell tractions from measured displacements.

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 \(x\)- \(y\) plane. An extension to curved surfaces would require intrinsic surface differential operators and a well-defined model for any normal traction.

How to cite us?

Inverse traction formulation

Let \(V\) be the elastic body, \(S_u\) the measurement surface, and \(S_\rho\) the surface on which the unknown traction acts. The state displacement is denoted by \(\mathbf{u}\), the measured displacement by \(\mathbf{u}_d\), the Lagrange multiplier (adjoint field) by \(\Upsilon\), and the traction by \(\mathbf{t}\). The small-strain constitutive equation is

\[ \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). \]

For a symmetric elasticity tensor \(\mathcal{A}\), define

\[ a(\mathbf{u},\mathbf{v}) = \left(\varepsilon(\mathbf{v}), \mathcal{A}:\varepsilon(\mathbf{u})\right)_V. \]

The inverse problem is

\[ \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} \]

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

\[ a(\mathbf{u},\mathbf{v})-(\mathbf{t},\mathbf{v})_{S_\rho}=0. \]

It follows from \(-(\Upsilon,\operatorname{div}\sigma)_V\); boundary terms on portions other than \(S_\rho\) vanish because of the corresponding essential or natural boundary conditions. The Lagrangian is therefore

\[ \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}. \]

Thus, traction is coupled to equilibrium in the weak sense.

For the quadratic regularisation

\[ R(\mathbf{t})= \frac{\epsilon_\rho^\mathrm{abs}}{2} (\mathbf{t},\mathbf{t})_{S_\rho}, \]

stationarity gives

\[ \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} \]

Finite element discretisation gives

\[ \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], \]

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:

\[ \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]. \]

Only the measured \(x\)- and \(y\)-components enter the displacement misfit. Consequently,

\[ 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}. \]

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:

\[ \mathbf{t}=\nabla_s\Phi , \]

where \(\nabla_s\) denotes the surface gradient. The regularisation and traction coupling become

\[ R(\Phi)=\frac{\epsilon_\rho^\mathrm{abs}}{2} (\nabla_s\Phi,\nabla_s\Phi)_{S_\rho},\qquad -(\nabla_s\Phi,\Upsilon)_{S_\rho}. \]

Variation with respect to the potential gives

\[ \epsilon_\rho^\mathrm{abs} (\nabla_s\Phi,\nabla_s\delta\Phi)_{S_\rho} -(\Upsilon,\nabla_s\delta\Phi)_{S_\rho}=0. \]

The third algebraic unknown in this formulation is \(\Phi\), not \(\mathbf{t}\). The matrix \(D\) is the surface stiffness matrix for the potential, and one value of \(\Phi\) must be fixed to remove its arbitrary additive constant.

Approximation

The state and adjoint fields satisfy \(\mathbf{u},\Upsilon\in[H^1(V)]^3\), with the appropriate essential boundary conditions. The potential satisfies \(\Phi\in H^1(S_\rho)\), and the recovered tangential traction \(\nabla_s\Phi\) belongs to \([L^2(S_\rho)]^2\).

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

\[ 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}. \]

The traction optimality equation is

\[ \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. \]

The command-line parameters are scaled in the implementation as

\[ \epsilon_\rho^\mathrm{abs} =\epsilon_\rho^\mathrm{ref}\epsilon_u^{-1},\qquad \epsilon_l^\mathrm{abs} =\epsilon_l^\mathrm{ref}\epsilon_u. \]

For positive \(\epsilon_l^\mathrm{ref}\), decreasing it strengthens the curl-free penalty. The exact value \(\epsilon_l^\mathrm{ref}=0\) disables the curl term in the current implementation to avoid division by zero; it does not select the limiting exact-potential problem. Use -my_curl 0 to select the exact force-potential formulation.

Approximation

The state and adjoint fields satisfy \(\mathbf{u},\Upsilon\in[H^1(V)]^3\). The traction belongs to the surface space \(\mathbf{t}\in H(\operatorname{curl}_s;S_\rho)\).

Field-split preconditioner

The implementation uses the block structure of the matrix through PCFIELDSPLIT and stores the matrix in the nested MATNEST format.

The nested matrix containing the two upper-left blocks is

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

The complete matrix is then written as

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

where

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

and

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

The system associated with \(X\) is solved using PCFIELDSPLIT and a multiplicative relaxation scheme. The submatrix \(K\) is factorised only once. The system associated with \(A\), composed of \(X\), \(D\), \(V\), and \(H\), is solved using a Schur complement. MoFEM subproblems are used to construct the required blocks.

Running code

Mapping measured displacements:

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

If only one layer is specified, a second thin polymer layer can be created automatically using prism elements:

./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

Configuration file:

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

Calculating tractions:

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

The resulting traction field is shown below:

Example: results

Installation with Docker

First, install Docker according to the instructions at: https://docs.docker.com/installation/#installation
Install the cell user module: docker pull 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.

Definition in file cell_forces.cpp.

Typedef Documentation

◆ DomainEleOp

using DomainEleOp = VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 377 of file cell_forces.cpp.

Function Documentation

◆ main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 453 of file cell_forces.cpp.

                                 {
 
  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;
}

Variable Documentation

◆ debug

int debug = 0

static

Definition at line 451 of file cell_forces.cpp.

◆ help

char help[]

static

Initial value:

= "-my_block_config set block data\n"

"\n"

Definition at line 381 of file cell_forces.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM = 3

constexpr

Definition at line 376 of file cell_forces.cpp.

Classes

Namespaces

Typedefs

Functions

Variables

Detailed Description

Force-potential implementation

Inverse traction formulation

Exact curl-free force-potential formulation

Approximation

Weak curl-free formulation

Approximation

Field-split preconditioner

Running code

Installation with Docker

Typedef Documentation

◆ DomainEleOp

Function Documentation

◆ main()

Variable Documentation

◆ debug

◆ help

◆ SPACE_DIM