simple_elasticity.cpp File Reference

#include <MoFEM.hpp>

Go to the source code of this file.

Classes
struct	OpK
struct	OpPressure
struct	ApplyDirichletBc
struct	VolRule
	Set integration rule to volume elements. More...
struct	FaceRule
	Set integration rule to boundary elements. More...

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

Variables
static char	help []

Function Documentation

main()

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

Examples

simple_elasticity.cpp.

Definition at line 390 of file simple_elasticity.cpp.

                                 {
 
  const string default_options = "-ksp_type fgmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_type mumps \n"
                                 "-mat_mumps_icntl_20 0 \n"
                                 "-ksp_monitor\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);
 
  // Create mesh database
  moab::Core mb_instance;              // create database
  moab::Interface &moab = mb_instance; // create interface to database
 
  try {
    // Create MoFEM database and link it to MoAB
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    CHKERR DMRegister_MoFEM("DMMOFEM");
 
    // Get command line options
    int order = 3;                    // default approximation order
    PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
                             "none");
 
    // Set approximation order
    CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
                           PETSC_NULL);
 
    // Set testing (used by CTest)
    CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
                            &flg_test, PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
 
    CHKERR simple_interface->getOptions();
    CHKERR simple_interface->loadFile();
 
    Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
 
    enum MyBcTypes {
      FIX_BRICK_FACES = 1,
      FIX_NODES = 2,
      BRICK_PRESSURE_FACES = 3,
      FIX_NODES_Y_DIR = 4
    };
 
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
      EntityHandle meshset = bit->getMeshset();
      int id = bit->getMeshsetId();
 
      if (id == FIX_BRICK_FACES) { // brick-faces
 
        CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
                                                            fix_faces, true);
 
        Range adj_ents;
        CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
                                                  moab::Interface::UNION);
 
        CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
                                                  moab::Interface::UNION);
        fix_faces.merge(adj_ents);
      } else if (id == FIX_NODES) { // node(s)
 
        CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
                                                            fix_nodes, true);
 
      } else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
        CHKERR m_field.get_moab().get_entities_by_dimension(
            meshset, 2, pressure_faces, true);
 
      } else if (id ==
                 FIX_NODES_Y_DIR) { // restrained second node in y direction
        CHKERR m_field.get_moab().get_entities_by_dimension(
            meshset, 0, fix_second_node, true);
 
      } else {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");
      }
    }
 
    CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
                                            3);
    CHKERR simple_interface->setFieldOrder("U", order);
 
    CHKERR simple_interface->defineFiniteElements();
 
    // Add pressure element
    CHKERR m_field.add_finite_element("PRESSURE");
    CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
    CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
    CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
 
    CHKERR simple_interface->defineProblem();
 
    DM dm;
    CHKERR simple_interface->getDM(&dm);
 
    CHKERR DMMoFEMAddElement(dm, "PRESSURE");
    CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
 
    CHKERR simple_interface->buildFields();
    CHKERR simple_interface->buildFiniteElements();
 
    CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
                                                     "PRESSURE");
    CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
 
    CHKERR simple_interface->buildProblem();
 
    // Create elements instances
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
        new VolumeElementForcesAndSourcesCore(m_field));
    boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
        new FaceElementForcesAndSourcesCore(m_field));
 
    // Set integration rule to elements instances
    elastic_fe->getRuleHook = VolRule();
    pressure_fe->getRuleHook = FaceRule();
 
    // Add operators to element instances
    // push operators to elastic_fe
    elastic_fe->getOpPtrVector().push_back(new OpK());
    // push operators to pressure_fe
    pressure_fe->getOpPtrVector().push_back(new OpPressure());
 
    boost::shared_ptr<FEMethod> fix_dofs_fe(
        new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));
 
    boost::shared_ptr<FEMethod> null_fe;
 
    // Set operators for KSP solver
    CHKERR DMMoFEMKSPSetComputeOperators(
        dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);
 
    CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_fe, null_fe,
                                   null_fe);
 
    // initialise matrix A used as the global stiffness matrix
    Mat A;
 
    // initialise left hand side vector x and right hand side vector f
    Vec x, f;
 
    // allocate memory handled by MoFEM discrete manager for matrix A
    CHKERR DMCreateMatrix(dm, &A);
 
    // allocate memory handled by MoFEM discrete manager for vector x
    CHKERR DMCreateGlobalVector(dm, &x);
 
    // allocate memory handled by MoFEM discrete manager for vector f of the
    // same size as x
    CHKERR VecDuplicate(x, &f);
 
    // precondition matrix A according to fix_dofs_fe  and elastic_fe finite
    // elements
    elastic_fe->ksp_B = A;
    fix_dofs_fe->ksp_B = A;
 
    // precondition the right hand side vector f according to fix_dofs_fe  and
    // elastic_fe finite elements
    fix_dofs_fe->ksp_f = f;
    pressure_fe->ksp_f = f;
 
    CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                    elastic_fe);
 
    CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_fe);
 
    // This is done because only post processor is implemented in the
    // ApplyDirichletBc struct
    CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());
 
    // make available a KSP solver
    KSP solver;
 
    // make the solver available for parallel computing by determining its MPI
    // communicator
    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
 
    // making available running all options available for KSP solver in running
    // command
    CHKERR KSPSetFromOptions(solver);
 
    // set A matrix with preconditioner
    CHKERR KSPSetOperators(solver, A, A);
 
    // set up the solver data strucure for the iterative solver
    CHKERR KSPSetUp(solver);
 
    // solve the system of linear equations
    CHKERR KSPSolve(solver, f, x);
 
    // destroy solver no needed any more
    CHKERR KSPDestroy(&solver);
 
    // make vector x available for parallel computations for visualization
    // context
    VecView(x, PETSC_VIEWER_STDOUT_WORLD);
 
    // save solution in vector x on mesh
    CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
 
    // Set up post-processor. It is some generic implementation of finite
    // element
    auto get_post_proc_ele = [&]() {
      auto jac_ptr = boost::make_shared<MatrixDouble>();
      auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
      auto det_ptr = boost::make_shared<VectorDouble>();
 
      auto post_proc_ele = boost::make_shared<
          PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>>(m_field);
      // Add operators to the elements, starting with some generic
      constexpr int SPACE_DIM = 3;
      post_proc_ele->getOpPtrVector().push_back(
          new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
      post_proc_ele->getOpPtrVector().push_back(
          new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
      post_proc_ele->getOpPtrVector().push_back(
          new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      auto grad_ptr = boost::make_shared<MatrixDouble>();
      post_proc_ele->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      post_proc_ele->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
                                                                   grad_ptr));
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      post_proc_ele->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
              {},
 
              {{"U", u_ptr}},
 
              {{"GRAD", grad_ptr}},
 
              {})
 
      );
 
      return post_proc_ele;
    };
 
    if (auto post_proc = get_post_proc_ele()) {
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                      post_proc);
      // write output
      CHKERR post_proc->writeFile("out.h5m");
    }
 
    {
      if (flg_test == PETSC_TRUE) {
        const double x_vec_norm_const = 0.4;
        // Check norm_1  value
        double norm_check;
        // Takes maximal element of the vector, that should be maximal
        // displacement at the end of the bar
        CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
        if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
            1.e-10) {
          SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                   "test failed (nrm 2 %6.4e)", norm_check);
        }
      }
    }
 
    // free memory handled by mofem discrete manager for A, x and f
    CHKERR MatDestroy(&A);
    CHKERR VecDestroy(&x);
    CHKERR VecDestroy(&f);
 
    // free memory allocated for mofem discrete manager
    CHKERR DMDestroy(&dm);
 
    // This is a good reference for the future
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

help

char help[]

static

Initial value:

= "-order approximation order\n"
                     "\n"

Examples

simple_elasticity.cpp.

Definition at line 13 of file simple_elasticity.cpp.