elasticity.cpp File Reference

#include <BasicFiniteElements.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes
struct	BlockOptionData
struct	VolRule
	Set integration rule to volume elements. More...
struct	PrismFE

Typedefs
using	BlockData = NonlinearElasticElement::BlockData
using	MassBlockData = ConvectiveMassElement::BlockData
using	VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
using	EdgeEle = MoFEM::EdgeElementForcesAndSourcesCore

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

Variables
static char	help []

Typedef Documentation

BlockData

using BlockData = NonlinearElasticElement::BlockData

Definition at line 74 of file elasticity.cpp.

EdgeEle

using EdgeEle = MoFEM::EdgeElementForcesAndSourcesCore

Definition at line 94 of file elasticity.cpp.

MassBlockData

using MassBlockData = ConvectiveMassElement::BlockData

Definition at line 75 of file elasticity.cpp.

VolUserDataOperator

using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 76 of file elasticity.cpp.

Function Documentation

main()

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

PETSC_VIEWER_ASCII_DENSE, PETSC_VIEWER_ASCII_MATLAB

Definition at line 96 of file elasticity.cpp.

                                 {
 
  const string default_options = "-ksp_type gmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_type mumps \n"
                                 "-mat_mumps_icntl_20 0 \n"
                                 "-ksp_monitor \n"
                                 "-snes_type newtonls \n"
                                 "-snes_linesearch_type basic \n"
                                 "-snes_atol 1e-8 \n"
                                 "-snes_rtol 1e-8 \n"
                                 "-snes_monitor \n"
                                 "-ts_monitor \n"
                                 "-ts_type beuler \n";
 
  string param_file = "param_file.petsc";
  if (!static_cast<bool>(ifstream(param_file))) {
    std::ofstream file(param_file.c_str(), std::ios::ate);
    if (file.is_open()) {
      file << default_options;
      file.close();
    }
  }
 
  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "ELASTIC"));
  LogManager::setLog("ELASTIC");
  MOFEM_LOG_TAG("ELASTIC", "elasticity")
 
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "ELASTIC_SYNC"));
  LogManager::setLog("ELASTIC_SYNC");
  MOFEM_LOG_TAG("ELASTIC_SYNC", "elastic_sync");
 
  try {
 
    PetscBool flg_block_config, flg_file;
    char mesh_file_name[255];
    char block_config_file[255];
    PetscInt test_nb = 0;
    PetscInt order = 2;
    PetscBool is_partitioned = PETSC_FALSE;
    PetscBool is_calculating_frequency = PETSC_FALSE;
    PetscBool is_post_proc_volume = PETSC_TRUE;
 
    // Select base
    enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
    const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
                                "jacobi"};
    PetscInt choice_base_value = LOBATTO;
 
    // Read options from command line
    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 PetscOptionsEList("-base", "approximation base", "", list_bases,
                             LASBASETOP, list_bases[choice_base_value],
                             &choice_base_value, PETSC_NULLPTR);
 
    CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "", test_nb,
                           &test_nb, PETSC_NULLPTR);
 
    CHKERR PetscOptionsBool("-my_is_partitioned",
                            "set if mesh is partitioned (this result that each "
                            "process keeps only one part of the mesh)",
                            "", is_partitioned, &is_partitioned, PETSC_NULLPTR);
 
    CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
                              "", "block_conf.in", block_config_file, 255,
                              &flg_block_config);
 
    CHKERR PetscOptionsBool(
        "-my_is_calculating_frequency", "set if frequency will be calculated",
        "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULLPTR);
 
    CHKERR PetscOptionsBool("-my_is_post_proc_volume",
                            "if true post proc volume", "", is_post_proc_volume,
                            &is_post_proc_volume, PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    // Throw error if file with mesh is not provided
    if (flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
    }
 
    // Create mesh database
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    // Create moab communicator
    // Create separate MOAB communicator, it is duplicate of PETSc communicator.
    // NOTE That this should eliminate potential communication problems between
    // MOAB and PETSC functions.
    MPI_Comm moab_comm_world;
    MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, moab_comm_world);
 
    // Read whole mesh or part of it if partitioned
    if (is_partitioned == PETSC_TRUE) {
      // This is a case of distributed mesh and algebra. In this case each
      // processor keeps only one part of the problem.
      const char *option;
      option = "PARALLEL=READ_PART;"
               "PARALLEL_RESOLVE_SHARED_ENTS;"
               "PARTITION=PARALLEL_PARTITION;";
      CHKERR moab.load_file(mesh_file_name, 0, option);
    } else {
      // If that case we have distributed algebra, i.e. assembly of vectors and
      // matrices is in parallel, but whole mesh is stored on all processors.
      // Solver and matrix scales well, however problem set-up of problem is
      // not fully parallel.
      const char *option;
      option = "";
      CHKERR moab.load_file(mesh_file_name, 0, option);
    }
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    // Print boundary conditions and material parameters
    MeshsetsManager *meshsets_mng_ptr;
    CHKERR m_field.getInterface(meshsets_mng_ptr);
    CHKERR meshsets_mng_ptr->printDisplacementSet();
    CHKERR meshsets_mng_ptr->printForceSet();
    CHKERR meshsets_mng_ptr->printMaterialsSet();
 
    bool mesh_has_tets = false;
    bool mesh_has_prisms = false;
    int nb_tets = 0;
    int nb_hexs = 0;
    int nb_prisms = 0;
 
    CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets, true);
    CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs, true);
    CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms, true);
 
    mesh_has_tets = (nb_tets + nb_hexs) > 0;
    mesh_has_prisms = nb_prisms > 0;
 
    // Set bit refinement level to all entities (we do not refine mesh in
    // this example so all entities have the same bit refinement level)
    BitRefLevel bit_level0;
    bit_level0.set(0);
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, bit_level0);
 
    // CHECK IF EDGE BLOCKSET EXIST AND IF IT IS ADD ALL ENTITIES FROM IT
    // CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
    // MESHSET_OF_EDGE_BLOCKSET, 1, bit_level0);
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
      if (bit->getName().compare(0, 3, "ROD") == 0) {
        CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
            0, 1, bit_level0);
      }
    }
 
    // Declare approximation fields
    FieldApproximationBase base = NOBASE;
    switch (choice_base_value) {
    case LEGENDRE:
      base = AINSWORTH_LEGENDRE_BASE;
      break;
    case LOBATTO:
      base = AINSWORTH_LOBATTO_BASE;
      break;
    case BERNSTEIN_BEZIER:
      base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
      break;
    case JACOBI:
      base = DEMKOWICZ_JACOBI_BASE;
      break;
    default:
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
    };
    CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_DENSE,
                             MF_ZERO);
 
    // We can use higher oder geometry to define body
    CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, base, 3, MB_TAG_DENSE,
                             MF_ZERO);
 
    // Declare problem
 
    // Add entities (by tets) to the field (all entities in the mesh, root_set
    // = 0 )
    CHKERR m_field.add_ents_to_field_by_dim(0, 3, "DISPLACEMENT");
    CHKERR m_field.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
 
    // Get all edges in the mesh
    Range all_edges;
    CHKERR m_field.get_moab().get_entities_by_type(0, MBEDGE, all_edges, true);
 
    // Get edges associated with simple rod
    Range edges_in_simple_rod;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
      if (bit->getName().compare(0, 3, "ROD") == 0) {
        Range edges;
        CHKERR m_field.get_moab().get_entities_by_type(bit->getMeshset(),
                                                       MBEDGE, edges, true);
        edges_in_simple_rod.merge(edges);
      }
    }
 
    CHKERR m_field.add_ents_to_field_by_type(edges_in_simple_rod, MBEDGE,
                                             "DISPLACEMENT");
 
    // Set order of edge in rod to be 1
    CHKERR m_field.set_field_order(edges_in_simple_rod, "DISPLACEMENT", 1);
 
    // Get remaining edges (not associated with simple rod) to set order
    Range edges_to_set_order;
    edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
 
    // Set approximation order.
    // See Hierarchical Finite Element Bases on Unstructured Tetrahedral
    // Meshes.
    CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBHEX, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
 
    if (base == AINSWORTH_BERNSTEIN_BEZIER_BASE)
      CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
    else
      CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
 
    // Set order of approximation of geometry.
    // Apply 2nd order only on skin (or in whole body)
    auto setting_second_order_geometry = [&m_field]() {
      MoFEMFunctionBegin;
      // Setting geometry order everywhere
      Range tets, edges;
      CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
      CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
                                                moab::Interface::UNION);
 
      // Setting 2nd geometry order only on skin
      // Range tets, faces, edges;
      // Skinner skin(&m_field.get_moab());
      // CHKERR skin.find_skin(0,tets,false,faces);
      // CHKERR m_field.get_moab().get_adjacencies(
      //   faces,1,false,edges,moab::Interface::UNION
      // );
      // CHKERR
      // m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
 
      CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
      CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
 
      MoFEMFunctionReturn(0);
    };
    CHKERR setting_second_order_geometry();
 
    // Configure blocks by parsing config file. It allows setting
    // approximation order for each block independently.
    std::map<int, BlockOptionData> block_data;
    auto setting_blocks_data_and_order_from_config_file =
        [&m_field, &moab, &block_data, flg_block_config, block_config_file,
         order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
          MoFEMFunctionBegin;
          if (flg_block_config) {
            ifstream ini_file(block_config_file);
            po::variables_map vm;
            po::options_description config_file_options;
            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));
 
              std::ostringstream str_init_temp;
              str_init_temp << "block_" << it->getMeshsetId()
                            << ".initial_temperature";
              config_file_options.add_options()(
                  str_init_temp.str().c_str(),
                  po::value<double>(&block_data[it->getMeshsetId()].initTemp)
                      ->default_value(0));
            }
            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;
              MOFEM_LOG_C("ELASTIC", Sev::inform, "Set block %d order to %d",
                          it->getMeshsetId(),
                          block_data[it->getMeshsetId()].oRder);
              Range block_ents;
              CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
                                                 true);
              Range ents_to_set_order;
              CHKERR moab.get_adjacencies(block_ents, 3, false,
                                          ents_to_set_order,
                                          moab::Interface::UNION);
              ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
              CHKERR moab.get_adjacencies(block_ents, 2, false,
                                          ents_to_set_order,
                                          moab::Interface::UNION);
              CHKERR moab.get_adjacencies(block_ents, 1, false,
                                          ents_to_set_order,
                                          moab::Interface::UNION);
              CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
                  ents_to_set_order);
 
              CHKERR m_field.set_field_order(
                  ents_to_set_order, "DISPLACEMENT",
                  block_data[it->getMeshsetId()].oRder);
            }
            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++) {
              MOFEM_LOG_C("ELASTIC", Sev::warning, "Unrecognized option %s",
                          vit->c_str());
            }
          }
 
          // Update material parameters. Set material parameters block by
          // block.
          for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
                   m_field, BLOCKSET | MAT_ELASTICSET, it)) {
            const int id = it->getMeshsetId();
            auto &bd = (*block_sets_ptr)[id];
            if (block_data[id].yOung > 0)
              bd.E = block_data[id].yOung;
            if (block_data[id].pOisson >= -1)
              bd.PoissonRatio = block_data[id].pOisson;
            MOFEM_LOG_C("ELASTIC", Sev::inform, "Block %d", id);
            MOFEM_LOG_C("ELASTIC", Sev::inform, "\tYoung modulus %3.4g", bd.E);
            MOFEM_LOG_C("ELASTIC", Sev::inform, "\tPoisson ratio %3.4g",
                        bd.PoissonRatio);
          }
 
          MoFEMFunctionReturn(0);
        };
 
    // Add elastic element
 
    boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
        boost::make_shared<std::map<int, HookeElement::BlockData>>();
    CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
    CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
 
    boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
        boost::make_shared<std::map<int, MassBlockData>>();
    CHKERR ConvectiveMassElement::setBlocks(m_field, mass_block_sets_ptr);
 
    auto fe_lhs_ptr =
        boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
    auto fe_rhs_ptr =
        boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
    fe_lhs_ptr->getRuleHook = VolRule();
    fe_rhs_ptr->getRuleHook = VolRule();
 
    CHKERR AddHOOps<3, 3, 3>::add(fe_lhs_ptr->getOpPtrVector(), {H1},
                                  "MESH_NODE_POSITIONS");
    CHKERR AddHOOps<3, 3, 3>::add(fe_rhs_ptr->getOpPtrVector(), {H1},
                                  "MESH_NODE_POSITIONS");
 
    boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
        new PrismFE(m_field));
    boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
        new PrismFE(m_field));
 
    CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
                                           "DISPLACEMENT",
                                           "MESH_NODE_POSITIONS", false);
 
    auto add_skin_element_for_post_processing = [&]() {
      MoFEMFunctionBegin;
      Range elastic_element_ents;
      CHKERR m_field.get_finite_element_entities_by_dimension(
          "ELASTIC", 3, elastic_element_ents);
      Skinner skin(&m_field.get_moab());
      Range skin_faces; // skin faces from 3d ents
      CHKERR skin.find_skin(0, elastic_element_ents, false, skin_faces);
      Range proc_skin;
      if (is_partitioned) {
        CHKERR pcomm->filter_pstatus(skin_faces,
                                     PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                     PSTATUS_NOT, -1, &proc_skin);
      } else {
        proc_skin = skin_faces;
      }
      CHKERR m_field.add_finite_element("POST_PROC_SKIN");
      CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
                                                          "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
                                                         "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
                                                          "DISPLACEMENT");
      if (m_field.check_field("TEMP")) {
        // CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
        //                                                    "TEMP");
        // CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
        //                                                    "TEMP");
        CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
                                                            "TEMP");
      }                                                    
      CHKERR m_field.modify_finite_element_add_field_data(
          "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
      CHKERR m_field.add_ents_to_finite_element_by_dim(proc_skin, 2,
                                                       "POST_PROC_SKIN");
      MoFEMFunctionReturn(0);
    };
    CHKERR add_skin_element_for_post_processing();
 
    auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
    if (mesh_has_tets) {
      CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
                                        "DISPLACEMENT", "MESH_NODE_POSITIONS",
                                        false, true, MBTET, data_at_pts);
    }
    if (mesh_has_prisms) {
      CHKERR HookeElement::setOperators(
          prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
          "MESH_NODE_POSITIONS", false, true, MBPRISM, data_at_pts);
    }
 
    if (test_nb == 4) {
 
      auto thermal_strain =
          [](FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords) {
            FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
            constexpr double alpha = 1;
            FTensor::Index<'i', 3> i;
            FTensor::Index<'k', 3> j;
            constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
            t_thermal_strain(i, j) = alpha * t_coords(2) * t_kd(i, j);
            return t_thermal_strain;
          };
 
      fe_rhs_ptr->getOpPtrVector().push_back(
          new HookeElement::OpAnalyticalInternalStrain_dx<0>(
              "DISPLACEMENT", data_at_pts, thermal_strain));
    }
 
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
        new VolumeElementForcesAndSourcesCore(m_field));
 
    for (auto &sit : *block_sets_ptr) {
      for (auto &mit : *mass_block_sets_ptr) {
        fe_mass_ptr->getOpPtrVector().push_back(
            new HookeElement::OpCalculateMassMatrix("DISPLACEMENT",
                                                    "DISPLACEMENT", sit.second,
                                                    mit.second, data_at_pts));
      }
    }
 
    // Add spring boundary condition applied on surfaces.
    // This is only declaration not implementation.
    CHKERR MetaSpringBC::addSpringElements(m_field, "DISPLACEMENT",
                                           "MESH_NODE_POSITIONS");
 
    // Implementation of spring element
    // Create new instances of face elements for springs
    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
        new FaceElementForcesAndSourcesCore(m_field));
    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
        new FaceElementForcesAndSourcesCore(m_field));
 
    CHKERR MetaSpringBC::setSpringOperators(m_field, fe_spring_lhs_ptr,
                                            fe_spring_rhs_ptr, "DISPLACEMENT",
                                            "MESH_NODE_POSITIONS");
 
    // Add Simple Rod elements
    // This is only declaration not implementation.
    CHKERR MetaSimpleRodElement::addSimpleRodElements(m_field, "DISPLACEMENT",
                                                      "MESH_NODE_POSITIONS");
 
    // CHKERR m_field.add_ents_to_finite_element_by_type(edges_in_simple_rod,
    //                                                   MBEDGE, "SIMPLE_ROD");
 
    // Implementation of Simple Rod element
    // Create new instances of edge elements for Simple Rod
    boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(new EdgeEle(m_field));
    boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(new EdgeEle(m_field));
 
    
    CHKERR MetaSimpleRodElement::setSimpleRodOperators(
        m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
        "MESH_NODE_POSITIONS");
 
    // Add body force element. This is only declaration of element. not its
    // implementation.
    CHKERR m_field.add_finite_element("BODY_FORCE");
    CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
                                                        "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
                                                        "MESH_NODE_POSITIONS");
 
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, BLOCKSET | BODYFORCESSET, it)) {
      Range tets;
      CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
                                                          true);
      CHKERR m_field.add_ents_to_finite_element_by_dim(tets, 3, "BODY_FORCE");
    }
    CHKERR m_field.build_finite_elements("BODY_FORCE");
 
    // Add Neumann forces, i.e. pressure or traction forces applied on body
    // surface. This is only declaration not implementation.
    CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
    CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
    CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
 
    // Add fluid pressure finite elements. This is special pressure on the
    // surface from fluid, i.e. pressure which linearly change with the depth.
    FluidPressure fluid_pressure_fe(m_field);
    // This function only declare element. Element is implemented by operators
    // in class FluidPressure.
    fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
 
    // Add elements for thermo elasticity if temperature field is defined.
    ThermalStressElement thermal_stress_elem(m_field);
    // Check if TEMP field exist, and then add element.
    if (!m_field.check_field("TEMP")) {
      bool add_temp_field = false;
      for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
        if (block_data[it->getMeshsetId()].initTemp != 0) {
          add_temp_field = true;
          break;
        }
      }
      if (add_temp_field) {
        CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
                                 MB_TAG_SPARSE, MF_ZERO);
 
        CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
        CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
      }
    }
    if (m_field.check_field("TEMP")) {
      CHKERR thermal_stress_elem.addThermalStressElement(
          "ELASTIC", "DISPLACEMENT", "TEMP");
    }
 
    // All is declared, at this point build fields first,
    CHKERR m_field.build_fields();
    // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
    // Projection10NodeCoordsOnField
    // do the trick.
    Projection10NodeCoordsOnField ent_method_material(m_field,
                                                      "MESH_NODE_POSITIONS");
    CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
    if (m_field.check_field("TEMP")) {
      for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
        if (block_data[it->getMeshsetId()].initTemp != 0) {
          MOFEM_LOG_C("ELASTIC", Sev::inform,
                      "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
                      block_data[it->getMeshsetId()].initTemp);
          Range block_ents;
          CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
          Range vertices;
          CHKERR moab.get_connectivity(block_ents, vertices, true);
          CHKERR m_field.getInterface<FieldBlas>()->setField(
              block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
              "TEMP");
        }
      }
    }
 
    // Build database for elements. Actual implementation of element is not need
    // here, only elements has to be declared.
    CHKERR m_field.build_finite_elements();
    // Build adjacencies between elements and field entities
    CHKERR m_field.build_adjacencies(bit_level0);
 
    // Register MOFEM DM implementation in PETSc
    CHKERR DMRegister_MGViaApproxOrders("MOFEM");
 
    // Create DM manager
    auto dm = createDM(PETSC_COMM_WORLD, "MOFEM");
    CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
    CHKERR DMSetFromOptions(dm);
    CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
    // Add elements to DM manager
    CHKERR DMMoFEMAddElement(dm, "ELASTIC");
    CHKERR DMMoFEMAddElement(dm, "SPRING");
    CHKERR DMMoFEMAddElement(dm, "SIMPLE_ROD");
    CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
    CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
    CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
    CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
    CHKERR DMMoFEMAddElement(dm, "POST_PROC_SKIN");
    CHKERR DMSetUp(dm);
 
    // Create matrices & vectors. Note that native PETSc DM interface is used,
    // but under the PETSc interface MoFEM implementation is running.
    SmartPetscObj<Vec> F;
    CHKERR DMCreateGlobalVector_MoFEM(dm, F);
    auto D = vectorDuplicate(F);
    auto D0 = vectorDuplicate(F);
    SmartPetscObj<Mat> Aij;
    CHKERR DMCreateMatrix_MoFEM(dm, Aij);
    CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
 
    // Initialise mass matrix
    SmartPetscObj<Mat> Mij;
    if (is_calculating_frequency == PETSC_TRUE) {
      Mij = matDuplicate(Aij, MAT_DO_NOT_COPY_VALUES);
      CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
      // MatView(Mij, PETSC_VIEWER_STDOUT_SELF);
    }
 
    // Assign global matrix/vector contributed by springs
    fe_spring_lhs_ptr->ksp_B = Aij;
    fe_spring_rhs_ptr->ksp_f = F;
 
    // Assign global matrix/vector contributed by Simple Rod
    fe_simple_rod_lhs_ptr->ksp_B = Aij;
    fe_simple_rod_rhs_ptr->ksp_f = F;
 
    // Zero vectors and matrices
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecZeroEntries(D);
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR MatZeroEntries(Aij);
 
    // Below particular implementations of finite elements are used to assemble
    // problem matrixes and vectors.  Implementation of element does not change
    // how element is declared.
 
    // Assemble Aij and F. Define Dirichlet bc element, which sets constrains
    // to MatrixDouble and the right hand side vector.
 
    // if normally defined boundary conditions are not found,
    // DirichletDisplacementBc will try to use DISPLACEMENT blockset. Two
    // implementations are available, depending how BC is defined on mesh file.
    auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
        m_field, "DISPLACEMENT", Aij, D0, F);
 
    // That sets Dirichlet bc objects that problem is linear, i.e. no newton
    // (SNES) solver is run for this problem.
    dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
    dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
 
    // D0 vector will store initial displacements
    CHKERR VecZeroEntries(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
    // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
    // implementation
    // of particular dirichlet_bc.
    CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
    // Set values from D0 on the field (on the mesh)
 
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
 
    // Calculate residual forces as result of applied kinematic constrains. Run
    // implementation
    // of particular finite element implementation. Look how
    // NonlinearElasticElement is implemented,
    // in that case. We run NonlinearElasticElement with hook material.
    // Calculate right hand side vector
    fe_rhs_ptr->snes_f = F;
    prism_fe_rhs_ptr->snes_f = F;
    MOFEM_LOG("ELASTIC", Sev::inform) << "Assemble external force vector  ...";
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
    MOFEM_LOG("ELASTIC", Sev::inform) << "done";
    // Assemble matrix
    fe_lhs_ptr->snes_B = Aij;
    prism_fe_lhs_ptr->snes_B = Aij;
    MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate stiffness matrix  ...";
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
    MOFEM_LOG("ELASTIC", Sev::inform) << "done";
 
    // Assemble springs
    CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_lhs_ptr);
    CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_rhs_ptr);
 
    // Assemble Simple Rod
    CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_lhs_ptr);
    CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_rhs_ptr);
 
    if (is_calculating_frequency == PETSC_TRUE) {
      // Assemble mass matrix
      fe_mass_ptr->snes_B = Mij;
      MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate mass matrix  ...";
      CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_mass_ptr);
      MOFEM_LOG("ELASTIC", Sev::inform) << "done";
    }
 
    // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
 
    // Assemble pressure and traction forces. Run particular implemented for do
    // this, see
    // MetaNeumannForces how this is implemented.
    boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
    CHKERR MetaNeumannForces::setMomentumFluxOperators(m_field, neumann_forces,
                                                       F, "DISPLACEMENT");
 
    {
      boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
          neumann_forces.begin();
      for (; mit != neumann_forces.end(); mit++) {
        CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
                                        &mit->second->getLoopFe());
      }
    }
    // Assemble forces applied to nodes, see implementation in NodalForce
    boost::ptr_map<std::string, NodalForce> nodal_forces;
    CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
                                         "DISPLACEMENT");
 
    {
      boost::ptr_map<std::string, NodalForce>::iterator fit =
          nodal_forces.begin();
      for (; fit != nodal_forces.end(); fit++) {
        CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
                                        &fit->second->getLoopFe());
      }
    }
    // Assemble edge forces
    boost::ptr_map<std::string, EdgeForce> edge_forces;
    CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
                                        "DISPLACEMENT");
    {
      auto fit = edge_forces.begin();
      for (; fit != edge_forces.end(); fit++) {
        auto &fe = fit->second->getLoopFe();
        CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(), &fe);
      }
    }
    // Assemble body forces, implemented in BodyForceConstantField
    BodyForceConstantField body_forces_methods(m_field);
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, BLOCKSET | BODYFORCESSET, it)) {
      CHKERR body_forces_methods.addBlock("DISPLACEMENT", F,
                                          it->getMeshsetId());
    }
    CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
                                    &body_forces_methods.getLoopFe());
    // Assemble fluid pressure forces
    CHKERR AddHOOps<2, 3, 3>::add(
        fluid_pressure_fe.getLoopFe().getOpPtrVector(), {},
        "MESH_NODE_POSITIONS");
    CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
        "DISPLACEMENT", F, false, true);
 
    CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
                                    &fluid_pressure_fe.getLoopFe());
    // Apply kinematic constrain to right hand side vector and matrix
    CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
 
    // Matrix View
    PetscViewerPushFormat(
        PETSC_VIEWER_STDOUT_SELF,
        PETSC_VIEWER_ASCII_MATLAB); /// PETSC_VIEWER_ASCII_DENSE,
                                    /// PETSC_VIEWER_ASCII_MATLAB
    // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
    // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
    // std::string wait;
    // std::cin >> wait;
 
    if (is_calculating_frequency == PETSC_TRUE) {
      CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
    }
 
    // Set matrix positive defined and symmetric for Cholesky and icc
    // pre-conditioner
 
    CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
    CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
    CHKERR VecScale(F, -1);
 
    // Create solver
    auto solver = createKSP(PETSC_COMM_WORLD);
    CHKERR KSPSetDM(solver, dm);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetOperators(solver, Aij, Aij);
 
    // Setup multi-grid pre-conditioner if set from command line
    {
      // from PETSc example ex42.c
      PetscBool same = PETSC_FALSE;
      PC pc;
      CHKERR KSPGetPC(solver, &pc);
      PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
      if (same) {
        CHKERR PCMGSetUpViaApproxOrders(
            pc, createPCMGSetUpViaApproxOrdersCtx(dm, Aij, true));
        CHKERR PCSetFromOptions(pc);
      } else {
        // Operators are already set, do not use DM for doing that
        CHKERR KSPSetDMActive(solver, PETSC_FALSE);
      }
    }
    CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
    CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
    // Set up solver
    CHKERR KSPSetUp(solver);
 
    using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
    auto set_post_proc_skin = [&](auto &post_proc_skin) {
      MoFEMFunctionBegin;
      CHKERR AddHOOps<2, 3, 3>::add(post_proc_skin.getOpPtrVector(), {},
                                    "MESH_NODE_POSITIONS");
      auto disp_ptr = boost::make_shared<MatrixDouble>();
      auto mesh_pos_ptr = boost::make_shared<MatrixDouble>();
      post_proc_skin.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("DISPLACEMENT", disp_ptr));
      post_proc_skin.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
                                               mesh_pos_ptr));
 
      auto side_loop =
          new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
              m_field, "ELASTIC", 3);
      CHKERR AddHOOps<3, 3, 3>::add(side_loop->getOpPtrVector(), {H1},
                                    "MESH_NODE_POSITIONS");
      side_loop->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENT",
                                                   data_at_pts->hMat));
      side_loop->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
                                                   data_at_pts->HMat));
      post_proc_skin.getOpPtrVector().push_back(side_loop);
 
      OpPPMap::DataMapVec scalar_map;
      OpPPMap::DataMapMat vector_map{
          {"DISPLACEMENT", disp_ptr},
          {"MESH_NODE_POSITIONS", mesh_pos_ptr}};
      OpPPMap::DataMapMat tensor_map{
          {"DISPLACEMENT_GRAD", data_at_pts->hMat}};
      if (m_field.check_field("TEMP")) {
        auto temp_ptr = boost::make_shared<VectorDouble>();
        auto temp_grad_ptr = boost::make_shared<MatrixDouble>();
        post_proc_skin.getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("TEMP", temp_ptr));
        post_proc_skin.getOpPtrVector().push_back(
            new OpCalculateScalarFieldGradient<3>("TEMP", temp_grad_ptr));
        scalar_map["TEMP"] = temp_ptr;
        vector_map["TEMP_GRAD"] = temp_grad_ptr;
      }
      post_proc_skin.getOpPtrVector().push_back(
          new HookeElement::OpPostProcHookeElement<
              FaceElementForcesAndSourcesCore>(
              "DISPLACEMENT", data_at_pts, *block_sets_ptr,
              post_proc_skin.getPostProcMesh(),
              post_proc_skin.getMapGaussPts(), true, true));
      post_proc_skin.getOpPtrVector().push_back(new OpPPMap(
          post_proc_skin.getPostProcMesh(), post_proc_skin.getMapGaussPts(),
          scalar_map, vector_map, tensor_map, {}));
      MoFEMFunctionReturn(0);
    };
 
    auto set_post_proc_tets = [&](auto &post_proc) {
      MoFEMFunctionBegin;
      // Add operators to the elements, starting with some generic operators
      CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {H1},
                                    "MESH_NODE_POSITIONS");
      auto disp_ptr = boost::make_shared<MatrixDouble>();
      auto mesh_pos_ptr = boost::make_shared<MatrixDouble>();
      auto disp_grad_ptr = boost::make_shared<MatrixDouble>();
      post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("DISPLACEMENT", disp_ptr));
      post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
                                               mesh_pos_ptr));
      post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENT",
                                                   disp_grad_ptr));
      OpPPMap::DataMapVec scalar_map;
      OpPPMap::DataMapMat vector_map{
          {"DISPLACEMENT", disp_ptr},
          {"MESH_NODE_POSITIONS", mesh_pos_ptr}};
      OpPPMap::DataMapMat tensor_map{
          {"DISPLACEMENT_GRAD", disp_grad_ptr}};
      if (m_field.check_field("TEMP")) {
        auto temp_ptr = boost::make_shared<VectorDouble>();
        auto temp_grad_ptr = boost::make_shared<MatrixDouble>();
        post_proc.getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("TEMP", temp_ptr));
        post_proc.getOpPtrVector().push_back(
            new OpCalculateScalarFieldGradient<3>("TEMP", temp_grad_ptr));
        scalar_map["TEMP"] = temp_ptr;
        vector_map["TEMP_GRAD"] = temp_grad_ptr;
      }
      // Add problem specific operator on element to post-process stresses
      post_proc.getOpPtrVector().push_back(new PostProcHookStress(
          m_field, post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
          post_proc.getPostProcElements(), "DISPLACEMENT", disp_grad_ptr,
          block_sets_ptr.get()));
      post_proc.getOpPtrVector().push_back(new OpPPMap(
          post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), scalar_map,
          vector_map, tensor_map, {}));
      MoFEMFunctionReturn(0);
    };
 
    auto set_post_proc_edge = [&](auto &post_proc_edge) {
      MoFEMFunctionBegin;
      auto disp_ptr = boost::make_shared<MatrixDouble>();
      post_proc_edge.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("DISPLACEMENT", disp_ptr));
      post_proc_edge.getOpPtrVector().push_back(new OpPPMap(
          post_proc_edge.getPostProcMesh(), post_proc_edge.getMapGaussPts(), {},
          {{"DISPLACEMENT", disp_ptr}}, {}, {}));
      MoFEMFunctionReturn(0);
    };
 
    auto set_post_proc_prisms = [&](auto &prism_post_proc) {
      MoFEMFunctionBegin;
      boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
      prism_post_proc.getOpPtrVector().push_back(
          new OpCalculateInvJacForFatPrism(inv_jac_ptr));
      prism_post_proc.getOpPtrVector().push_back(
          new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
      auto disp_ptr = boost::make_shared<MatrixDouble>();
      auto mesh_pos_ptr = boost::make_shared<MatrixDouble>();
      auto disp_grad_ptr = boost::make_shared<MatrixDouble>();
      prism_post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("DISPLACEMENT", disp_ptr));
      prism_post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
                                               mesh_pos_ptr));
      prism_post_proc.getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENT",
                                                   disp_grad_ptr));
      prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
          m_field, prism_post_proc.getPostProcMesh(),
          prism_post_proc.getMapGaussPts(),
          prism_post_proc.getPostProcElements(), "DISPLACEMENT",
          disp_grad_ptr, block_sets_ptr.get()));
      prism_post_proc.getOpPtrVector().push_back(new OpPPMap(
          prism_post_proc.getPostProcMesh(),
          prism_post_proc.getMapGaussPts(), {},
          {{"DISPLACEMENT", disp_ptr},
           {"MESH_NODE_POSITIONS", mesh_pos_ptr}},
          {{"DISPLACEMENT_GRAD", disp_grad_ptr}}, {}));
      MoFEMFunctionReturn(0);
    };
 
    PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore> post_proc_skin(
        m_field);
    BasicPostProcFatPrism prism_post_proc(m_field);
    PostProcBrokenMeshInMoab<EdgeElementForcesAndSourcesCore> post_proc_edge(
        m_field);
    PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
        m_field);
 
    auto is_prism = [](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getEntType() == MBPRISM;
    };
    prism_post_proc.exeTestHook = is_prism;
    post_proc.exeTestHook = [is_prism](FEMethod *fe_ptr) {
      return !is_prism(fe_ptr);
    };
 
    CHKERR set_post_proc_skin(post_proc_skin);
    CHKERR set_post_proc_tets(post_proc);
    CHKERR set_post_proc_prisms(prism_post_proc);
    CHKERR set_post_proc_edge(post_proc_edge);
 
    PetscBool field_eval_flag = PETSC_FALSE;
    std::array<double, 3> field_eval_coords;
    boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> field_eval_data;
    PetscInt coords_dim = 3;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                    field_eval_coords.data(), &coords_dim,
                                    &field_eval_flag);
 
    auto scalar_field_ptr = boost::make_shared<VectorDouble>();
    auto vector_field_ptr = boost::make_shared<MatrixDouble>();
    auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
 
    if (field_eval_flag) {
      field_eval_data = m_field.getInterface<FieldEvaluatorInterface>()
                            ->getData<VolumeElementForcesAndSourcesCore>();
      CHKERR m_field.getInterface<FieldEvaluatorInterface>()->buildTree<3>(
          field_eval_data, "ELASTIC");
 
      field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
      auto no_rule = [](int, int, int) { return -1; };
 
      auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
      field_eval_fe_ptr->getRuleHook = no_rule;
 
      if (m_field.check_field("TEMP")) {
        field_eval_fe_ptr->getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("TEMP", scalar_field_ptr));
      }
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("DISPLACEMENT", vector_field_ptr));
       field_eval_fe_ptr->getOpPtrVector().push_back(
         new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENT", tensor_field_ptr));
    }
 
    // Temperature field is defined on the mesh
    if (m_field.check_field("TEMP")) {
 
      // Create thermal vector
      Vec F_thermal;
      CHKERR VecDuplicate(F, &F_thermal);
 
      // Set up implementation for calculation of thermal stress vector. Look
      // how thermal stresses and vector is assembled in ThermalStressElement.
      CHKERR thermal_stress_elem.setThermalStressRhsOperators(
          "DISPLACEMENT", "TEMP", F_thermal);
 
      SeriesRecorder *recorder_ptr;
      CHKERR m_field.getInterface(recorder_ptr);
 
      // Read time series and do thermo-elastic analysis, this is when time
      // dependent
      // temperature problem was run before on the mesh. It means that before
      // non-stationary
      // problem was solved for temperature and filed "TEMP" is stored for
      // subsequent time
      // steps in the recorder.
      if (recorder_ptr->check_series("THEMP_SERIES")) {
        // This is time dependent case, so loop of data series stored by tape
        // recorder.
        // Loop over time steps
        for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
                                                sit)) {
          MOFEM_LOG_C("ELASTIC", Sev::inform, "Process step %d",
                      sit->get_step_number());
          // Load field data for this time step
          CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
                                                sit->get_step_number());
 
          CHKERR VecZeroEntries(F_thermal);
          CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
          CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
 
          // Calculate the right-hand side vector as result of thermal stresses.
          // It MetaNodalForces
          // that on "ELASTIC" element data structure the element implementation
          // in thermal_stress_elem
          // is executed.
          CHKERR DMoFEMLoopFiniteElements(
              dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
 
          // Assemble vector
          CHKERR VecAssemblyBegin(F_thermal);
          CHKERR VecAssemblyEnd(F_thermal);
          // Accumulate ghost dofs
          CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
          CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
 
          // Calculate norm of vector and print values
          PetscReal nrm_F;
          CHKERR VecNorm(F, NORM_2, &nrm_F);
 
          MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
          PetscReal nrm_F_thermal;
          CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
          MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
                      nrm_F_thermal);
 
          CHKERR VecScale(F_thermal, -1);
          // check this !!!
          CHKERR VecAXPY(F_thermal, 1, F);
 
          // Set dirichlet boundary to thermal stresses vector
          dirichlet_bc_ptr->snes_x = D;
          dirichlet_bc_ptr->snes_f = F_thermal;
          CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
 
          // Solve problem
          CHKERR KSPSolve(solver, F_thermal, D);
          // Add boundary conditions vector
          CHKERR VecAXPY(D, 1., D0);
          CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
          CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
          // Save data on the mesh
          CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
          // Save data on mesh
          CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
 
          if (field_eval_flag) {
            CHKERR m_field.getInterface<FieldEvaluatorInterface>()
                ->evalFEAtThePoint<3>(
                    field_eval_coords.data(), 1e-12, "ELASTIC_PROB", "ELASTIC",
                    field_eval_data, m_field.get_comm_rank(),
                    m_field.get_comm_rank(), nullptr, MF_EXIST, QUIET);
            if (scalar_field_ptr->size()) {
              auto t_temp = getFTensor0FromVec(*scalar_field_ptr);
              MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
                  << "Eval point TEMP: " << t_temp;
            }
            if (vector_field_ptr->size1()) {
              FTensor::Index<'i', 3> i;
              auto t_disp = getFTensor1FromMat<3>(*vector_field_ptr);
              MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
                  << "Eval point DISPLACEMENT magnitude: "
                  << sqrt(t_disp(i) * t_disp(i));
            }
            if (tensor_field_ptr->size1()) {
              FTensor::Index<'i', 3> i;
              auto t_disp_grad = getFTensor2FromMat<3, 3>(*tensor_field_ptr);
              MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
                  << "Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(i, i);
            }
 
            MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
          }
 
          // Post-process results
          if (is_post_proc_volume == PETSC_TRUE) {
            MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
            CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
            std::ostringstream o1;
            o1 << "out_" << sit->step_number << ".h5m";
            if (!test_nb)
              CHKERR post_proc.writeFile(o1.str().c_str());
            MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
          }
 
          MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
          CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN",
                                          &post_proc_skin);
          std::ostringstream o1_skin;
          o1_skin << "out_skin_" << sit->step_number << ".h5m";
          if (!test_nb)
            CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
          MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
        }
      } else {
 
        // This is a case when stationary problem for temperature was solved.
        CHKERR VecZeroEntries(F_thermal);
        CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
 
        // Calculate the right-hand side vector with thermal stresses
        CHKERR DMoFEMLoopFiniteElements(
            dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
 
        // Assemble vector
        CHKERR VecAssemblyBegin(F_thermal);
        CHKERR VecAssemblyEnd(F_thermal);
 
        // Accumulate ghost dofs
        CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
 
        // Calculate norm
        PetscReal nrm_F;
        CHKERR VecNorm(F, NORM_2, &nrm_F);
 
        MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
        PetscReal nrm_F_thermal;
        CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
 
        MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
                    nrm_F_thermal);
 
        // Add thermal stress vector and other forces vector
        CHKERR VecScale(F_thermal, -1);
        CHKERR VecAXPY(F_thermal, 1, F);
 
        // Apply kinetic boundary conditions
        dirichlet_bc_ptr->snes_x = D;
        dirichlet_bc_ptr->snes_f = F_thermal;
        CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
 
        // Solve problem
        CHKERR KSPSolve(solver, F_thermal, D);
        CHKERR VecAXPY(D, 1., D0);
 
        // Update ghost values for solution vector
        CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
        // Save data on mesh
        if (is_post_proc_volume == PETSC_TRUE) {
          MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
          CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
          // Save results to file
          if (!test_nb)
            CHKERR post_proc.writeFile("out.h5m");
          MOFEM_LOG("ELASTIC", Sev::inform) << "done";
        }
 
        MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
        CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
        if (!test_nb)
          CHKERR post_proc_skin.writeFile("out_skin.h5m");
        MOFEM_LOG("ELASTIC", Sev::inform) << "done";
      }
 
      // Destroy vector, no needed any more
      CHKERR VecDestroy(&F_thermal);
    } else {
      // Elastic analysis no temperature field
      // VecView(F, PETSC_VIEWER_STDOUT_WORLD);
      // Solve for vector D
      CHKERR KSPSolve(solver, F, D);
 
      // VecView(D, PETSC_VIEWER_STDOUT_WORLD);
      // cerr << F;
 
      // Add kinetic boundary conditions
      CHKERR VecAXPY(D, 1., D0);
      // Update ghost values
      CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
      // Save data from vector on mesh
      CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
      // Post-process results
      MOFEM_LOG("ELASTIC", Sev::inform) << "Post-process start ...";
      if (is_post_proc_volume == PETSC_TRUE) {
        CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
      }
      CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
      CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
      CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
      MOFEM_LOG("ELASTIC", Sev::inform) << "done";
 
      // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
      MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
      if (mesh_has_tets) {
        if (is_post_proc_volume == PETSC_TRUE) {
          if (!test_nb)
            CHKERR post_proc.writeFile("out.h5m");
        }
        if (!test_nb)
          CHKERR post_proc_skin.writeFile("out_skin.h5m");
      }
      if (mesh_has_prisms) {
        if (!test_nb)
          CHKERR prism_post_proc.writeFile("prism_out.h5m");
      }
      if (!edges_in_simple_rod.empty())
        if (!test_nb)
          CHKERR post_proc_edge.writeFile("out_edge.h5m");
      MOFEM_LOG("ELASTIC", Sev::inform) << "done";
    }
 
    if (is_calculating_frequency == PETSC_TRUE) {
      // Calculate mode mass, m = u^T * M * u
      Vec u1;
      VecDuplicate(D, &u1);
      CHKERR MatMult(Mij, D, u1);
      double mode_mass;
      CHKERR VecDot(u1, D, &mode_mass);
      MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode mass  %6.4e\n", mode_mass);
 
      Vec v1;
      VecDuplicate(D, &v1);
      CHKERR MatMult(Aij, D, v1);
 
      double mode_stiffness;
      CHKERR VecDot(v1, D, &mode_stiffness);
      MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode stiffness  %6.4e\n",
                  mode_stiffness);
 
      double frequency;
      double pi = 3.14159265359;
      frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
      MOFEM_LOG_C("ELASTIC", Sev::inform, "Frequency  %6.4e", frequency);
    }
 
    // Calculate elastic energy
    auto calculate_strain_energy = [&]() {
      MoFEMFunctionBegin;
 
      SmartPetscObj<Vec> v_energy;
      CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
                                           "MESH_NODE_POSITIONS", false, true,
                                           v_energy);
 
      // Print elastic energy
      double energy;
      CHKERR VecSum(v_energy, &energy);
      MOFEM_LOG_C("ELASTIC", Sev::inform, "Elastic energy %6.4e", energy);
 
      switch (test_nb) {
      case 1:
        if (fabs(energy - 17.129) > 1e-3)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
        break;
      case 2:
        if (fabs(energy - 5.6475e-03) > 1e-4)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
        break;
      case 3:
        if (fabs(energy - 7.4679e-03) > 1e-4)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
        break;
      case 4:
        if (fabs(energy - 2.4992e+00) > 1e-3)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
        break;
      // FIXME: Here are missing regersion tests
      case 8: {
        double min;
        CHKERR VecMin(D, PETSC_NULLPTR, &min);
        constexpr double expected_val = 0.10001;
        if (fabs(min + expected_val) > 1e-10)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged! %3.4e != %3.4e", min, expected_val);
      } break;
      case 9: {
        if (fabs(energy - 4.7416e-04) > 1e-8)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
      }
      default:
        break;
      }
 
      MoFEMFunctionReturn(0);
    };
    CHKERR calculate_strain_energy();
 
    MPI_Comm_free(&moab_comm_world);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

help

char help[]

static

Initial value:

= "-my_block_config set block data\n"
                     "-my_order approximation order\n"
                     "-my_is_partitioned set if mesh is partitioned\n"
                     "\n"

Definition at line 57 of file elasticity.cpp.