elasticity.cpp File Reference

#include <BasicFiniteElements.hpp>
#include <boost/program_options.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes
struct	BlockOptionData
	stores some optional data More...
struct	VolRule
	Set integration rule. More...
struct	PrismFE

Typedefs
using	BlockData = NonlinearElasticElement::BlockData
using	MassBlockData = ConvectiveMassElement::BlockData
using	VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator

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

Variables
static char	help []

Typedef Documentation

BlockData

using BlockData = NonlinearElasticElement::BlockData

Examples

elasticity.cpp, testing_jacobian_of_hook_element.cpp, and testing_jacobian_of_hook_scaled_with_density_element.cpp.

Definition at line 89 of file elasticity.cpp.

MassBlockData

using MassBlockData = ConvectiveMassElement::BlockData

Definition at line 90 of file elasticity.cpp.

VolUserDataOperator

using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator

Examples

testing_jacobian_of_hook_scaled_with_density_element.cpp.

Definition at line 91 of file elasticity.cpp.

Function Documentation

main()

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

Examples

elasticity.cpp.

Definition at line 109 of file elasticity.cpp.

                                 {

  const string default_options = "-ksp_type gmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_package mumps \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);

  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;

    // Select base
    enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
    const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier"};
    PetscInt choice_base_value = LOBATTO;

    // Read options from command line
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
    CHKERR(ierr);
    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_NULL);

    CHKERR PetscOptionsEList("-base", "approximation base", "", list_bases,
                             LASBASETOP, list_bases[choice_base_value],
                             &choice_base_value, PETSC_NULL);

    CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "", test_nb,
                           &test_nb, PETSC_NULL);

    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_NULL);

    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_NULL);

    ierr = PetscOptionsEnd();
    CHKERRG(ierr);

    // 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_prisms = 0;

    CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets, true);
    CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms, true);

    mesh_has_tets = nb_tets > 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;
    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, AINSWORTH_LEGENDRE_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, MBTRI, "DISPLACEMENT", order);
    // CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
    CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);

    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;
              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);
              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++) {
              CHKERR PetscPrintf(PETSC_COMM_WORLD,
                                 "** WARNING Unrecognized option %s\n",
                                 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;
            CHKERR PetscPrintf(PETSC_COMM_WORLD, "Block %d\n", id);
            CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung modulus %3.4g\n",
                               bd.E);
            CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson ratio %3.4g\n",
                               bd.PoissonRatio);
          }

          MoFEMFunctionReturn(0);
        };

    // Add elastic element

    boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
        boost::make_shared<std::map<int, 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);

    boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
        new VolumeElementForcesAndSourcesCore(m_field));
    boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
        new VolumeElementForcesAndSourcesCore(m_field));
    fe_lhs_ptr->getRuleHook = VolRule();
    fe_rhs_ptr->getRuleHook = VolRule();

    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);
    if (mesh_has_tets) {
      CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
                                        "DISPLACEMENT", "MESH_NODE_POSITIONS",
                                        false, true);
    }
    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);
    }

    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
        new VolumeElementForcesAndSourcesCore(m_field));

    boost::shared_ptr<HookeElement::DataAtIntegrationPts> data_at_pts(
        new HookeElement::DataAtIntegrationPts());

    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<EdgeElementForcesAndSourcesCore> fe_simple_rod_lhs_ptr(
        new EdgeElementForcesAndSourcesCore(m_field));
    boost::shared_ptr<EdgeElementForcesAndSourcesCore> fe_simple_rod_rhs_ptr(
        new EdgeElementForcesAndSourcesCore(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) {
          PetscPrintf(PETSC_COMM_WORLD, "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
    DM dm;
    CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
    CHKERR DMSetType(dm, "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 DMSetUp(dm);

    // Create matrices & vectors. Note that native PETSc DM interface is used,
    // but under the PETSc interface MoFEM implementation is running.
    Vec F, D, D0;
    CHKERR DMCreateGlobalVector(dm, &F);
    CHKERR VecDuplicate(F, &D);
    CHKERR VecDuplicate(F, &D0);
    Mat Aij;
    CHKERR DMCreateMatrix(dm, &Aij);
    CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);

    // Initialise mass matrix
    Mat Mij;
    if (is_calculating_frequency == PETSC_TRUE) {
      CHKERR MatDuplicate(Aij, MAT_DO_NOT_COPY_VALUES, &Mij);
      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);

    // This controls how kinematic constrains are set, by blockset or nodeset.
    // Cubit sets kinetic boundary conditions by blockset.
    bool flag_cubit_disp = false;
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, NODESET | DISPLACEMENTSET, it)) {
      flag_cubit_disp = true;
    }

    // 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 stets constrains
    // to MatrixDouble and the right hand side vector.
    boost::shared_ptr<FEMethod> dirichlet_bc_ptr;

    // if normally defined boundary conditions are not found, try to use
    // DISPLACEMENT blockset. To implementations available here, depending how
    // BC is defined on mesh file.
    if (!flag_cubit_disp) {
      dirichlet_bc_ptr =
          boost::shared_ptr<FEMethod>(new DirichletSetFieldFromBlockWithFlags(
              m_field, "DISPLACEMENT", "DISPLACEMENT", Aij, D0, F));
    } else {
      dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(
          new 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;
    PetscPrintf(PETSC_COMM_WORLD, "Assemble external force vector  ...");
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
    PetscPrintf(PETSC_COMM_WORLD, " done\n");
    // Assemble matrix
    fe_lhs_ptr->snes_B = Aij;
    prism_fe_lhs_ptr->snes_B = Aij;
    PetscPrintf(PETSC_COMM_WORLD, "Calculate stiffness matrix  ...");
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
    PetscPrintf(PETSC_COMM_WORLD, " done\n");

    // 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;
      PetscPrintf(PETSC_COMM_WORLD, "Calculate mass matrix  ...");
      CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_mass_ptr);
      PetscPrintf(PETSC_COMM_WORLD, " done\n");
    }

    // 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");
    {
      boost::ptr_map<std::string, EdgeForce>::iterator fit =
          edge_forces.begin();
      for (; fit != edge_forces.end(); fit++) {
        CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
                                        &fit->second->getLoopFe());
      }
    }
    // 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 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
    KSP solver;
    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
    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) {
        PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
        CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
        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);

    // Set up post-processor. It is some generic implementation of finite
    // element.
    PostProcVolumeOnRefinedMesh post_proc(m_field);
    // Add operators to the elements, starting with some generic operators
    CHKERR post_proc.generateReferenceElementMesh();
    CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
    CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
    CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
    // Add problem specific operator on element to post-process stresses
    post_proc.getOpPtrVector().push_back(new PostProcHookStress(
        m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "DISPLACEMENT",
        post_proc.commonData, block_sets_ptr.get()));

    PostProcEdgeOnRefinedMesh post_proc_edge(m_field);
    CHKERR post_proc_edge.generateReferenceElementMesh();
    CHKERR post_proc_edge.addFieldValuesPostProc("DISPLACEMENT");

    PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
    CHKERR prism_post_proc.generateReferenceElementMesh();

    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));

    CHKERR prism_post_proc.addFieldValuesPostProc("DISPLACEMENT");
    CHKERR prism_post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
    CHKERR prism_post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
    prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
        m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
        "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));

    // 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)) {
          PetscPrintf(PETSC_COMM_WORLD, "Process step %d\n",
                      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);

          PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
          PetscReal nrm_F_thermal;
          CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
          PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
                      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());
          // Post-process results
          CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);

          std::ostringstream o1;
          o1 << "out_" << sit->step_number << ".h5m";
          CHKERR post_proc.writeFile(o1.str().c_str());
        }
      } 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);

        PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
        PetscReal nrm_F_thermal;
        CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);

        PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
                    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);

        // Save data on mesh
        CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
        CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
        // Save results to file
        PetscPrintf(PETSC_COMM_WORLD, "Write output file ...");
        CHKERR post_proc.writeFile("out.h5m");
        PetscPrintf(PETSC_COMM_WORLD, " done\n");
      }

      // 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
      CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
      CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
      // CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", &post_proc);
      CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
      // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
      PetscPrintf(PETSC_COMM_WORLD, "Write output file ...");
      if (mesh_has_tets) {
        CHKERR post_proc.writeFile("out.h5m");
      }
      if (mesh_has_prisms) {
        CHKERR prism_post_proc.writeFile("prism_out.h5m");
      }
      if (!edges_in_simple_rod.empty())
        CHKERR post_proc_edge.writeFile("out_edge.h5m");
      PetscPrintf(PETSC_COMM_WORLD, " done\n");
    }

    if (is_calculating_frequency == PETSC_TRUE) {
      // Calculate model mass, m = u^T * M * u
      Vec u1;
      VecDuplicate(D, &u1);
      CHKERR MatMult(Mij, D, u1);
      double model_mass;
      CHKERR VecDot(u1, D, &model_mass);
      // PetscPrintf(PETSC_COMM_WORLD, "Model mass  %6.4e\n", model_mass);

      // Calculate model stiffness, k = v^T * K * v where v = u / norm(u)
      double u_norm;
      CHKERR VecNorm(D, NORM_2, &u_norm);

      // CHKERR VecScale(D,1./u_norm);     // v obtained
      CHKERR VecScale(D, 1. / 1.); // v obtained
      // PetscPrintf(PETSC_COMM_WORLD, "u_norm  %6.4e\n", u_norm);

      Vec v1;
      VecDuplicate(D, &v1);
      CHKERR MatMult(Aij, D, v1);

      double model_stiffness;
      CHKERR VecDot(v1, D, &model_stiffness);
      // PetscPrintf(PETSC_COMM_WORLD, "Model stiffness  %6.4e\n",
      //             model_stiffness);

      double frequency;
      double pi = 3.14159265359;
      frequency = sqrt(model_stiffness / model_mass) / (2 * pi);
      PetscPrintf(PETSC_COMM_WORLD, "Frequency  %6.4e\n", frequency);
    }

    // Calculate elastic energy
    auto calculate_strain_energy = [dm, &block_sets_ptr, test_nb]() {
      MoFEMFunctionBegin;

      Vec v_energy;
      CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
                                           "MESH_NODE_POSITIONS", false, true,
                                           &v_energy);

      // Print elastic energy
      const double *eng_ptr;
      CHKERR VecGetArrayRead(v_energy, &eng_ptr);
      PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n", *eng_ptr);

      switch (test_nb) {
      case 1:
        if (fabs(*eng_ptr - 17.129) > 1e-3)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
        break;
      case 2:
        if (fabs(*eng_ptr - 5.6475e-03) > 1e-3)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "atom test diverged!");
        break;

      default:
        break;
      }

      CHKERR VecRestoreArrayRead(v_energy, &eng_ptr);
      CHKERR VecDestroy(&v_energy);
      MoFEMFunctionReturn(0);
    };
    CHKERR calculate_strain_energy();

    // Destroy matrices, vecors, solver and DM
    CHKERR VecDestroy(&F);
    CHKERR VecDestroy(&D);
    CHKERR VecDestroy(&D0);
    CHKERR MatDestroy(&Aij);
    CHKERR KSPDestroy(&solver);
    CHKERR DMDestroy(&dm);

    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"

Examples

elasticity.cpp.

Definition at line 72 of file elasticity.cpp.