v0.13.0
mortar_contact.cpp

Implementation of mortar contact between surfaces with non-matching meshes

/** \file mortar_contact.cpp
* \example mortar_contact.cpp
*
* Implementation of mortar contact between surfaces with non-matching meshes
*
**/
/* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <Mortar.hpp>
// #include <boost/program_options.hpp>
// namespace po = boost::program_options;
static char help[] = "\n";
int main(int argc, char *argv[]) {
const string default_options = "-ksp_type fgmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_13 1 \n"
"-mat_mumps_icntl_14 800 \n"
"-mat_mumps_icntl_20 0 \n"
"-mat_mumps_icntl_24 1 \n"
"-snes_type newtonls \n"
"-snes_linesearch_type basic \n"
"-snes_divergence_tolerance 0 \n"
"-snes_max_it 50 \n"
"-snes_atol 1e-8 \n"
"-snes_rtol 1e-10 \n"
"-snes_monitor \n"
"-ksp_monitor \n"
"-snes_converged_reason \n"
"-order 1 \n"
"-order_lambda 1 \n"
"-order_contact 2 \n"
"-ho_levels_num 1 \n"
"-step_num 1 \n"
"-cn_value 1. \n"
"-alm_flag 0 \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();
}
}
// Initialize MoFEM
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 {
PetscBool flg_file;
char mesh_file_name[255];
PetscInt order = 1;
PetscInt nb_sub_steps = 1;
PetscBool is_partitioned = PETSC_FALSE;
PetscInt test_num = 0;
PetscBool wave_surf_flag = PETSC_FALSE;
PetscInt wave_dim = 2;
PetscBool is_displacement_field = PETSC_FALSE;
PetscInt wave_surf_block_id = 1;
PetscReal wave_length = 1.0;
PetscReal wave_ampl = 0.01;
PetscReal mesh_height = 1.0;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
"", 1, &order, PETSC_NULL);
CHKERR PetscOptionsInt("-step_num", "number of steps", "", nb_sub_steps,
&nb_sub_steps, PETSC_NULL);
CHKERR PetscOptionsBool("-is_partitioned",
"set if mesh is partitioned (this result that each "
"process keeps only part of the mes",
"", PETSC_FALSE, &is_partitioned, PETSC_NULL);
CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
PETSC_NULL);
CHKERR PetscOptionsBool("-is_displacement_field",
"use displacement field", "",
PETSC_FALSE, &is_displacement_field, PETSC_NULL);
CHKERR PetscOptionsBool("-wave_surf",
"if set true, make one of the surfaces wavy", "",
PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
CHKERR PetscOptionsInt("-wave_surf_block_id",
"make wavy surface of the block with this id", "",
wave_surf_block_id, &wave_surf_block_id, PETSC_NULL);
CHKERR PetscOptionsInt("-wave_dim", "dimension (2 or 3)", "", wave_dim,
&wave_dim, PETSC_NULL);
CHKERR PetscOptionsReal("-wave_length", "profile wavelength", "",
wave_length, &wave_length, PETSC_NULL);
CHKERR PetscOptionsReal("-wave_ampl", "profile amplitude", "", wave_ampl,
&wave_ampl, PETSC_NULL);
CHKERR PetscOptionsReal("-mesh_height", "vertical dimension of the mesh ",
"", mesh_height, &mesh_height, PETSC_NULL);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
// Check if mesh file was provided
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
}
// FIXME: really???
if (is_partitioned == PETSC_TRUE)
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Partitioned mesh is not supported");
const char *option;
option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
MeshsetsManager *mmanager_ptr;
CHKERR m_field.getInterface(mmanager_ptr);
CHKERR mmanager_ptr->printDisplacementSet();
CHKERR mmanager_ptr->printForceSet();
// print block sets with materials
CHKERR mmanager_ptr->printMaterialsSet();
string primary_field_name;
if(!is_displacement_field)
primary_field_name = "SPATIAL_POSITION";
else
primary_field_name = "U";
// create interfaces (could be based on cmake file ?)
MortarContactInterface m_contact(m_field, primary_field_name,
"MESH_NODE_POSITIONS",
is_displacement_field);
NonlinearElasticElementInterface m_elastic(m_field, primary_field_name,
"MESH_NODE_POSITIONS",
is_displacement_field);
CHKERR m_contact.addElementFields();
CHKERR m_elastic.addElementFields();
// build field
CHKERR m_field.build_fields();
CHKERR m_contact.createElements();
CHKERR m_elastic.createElements();
CHKERR m_contact.setOperators();
CHKERR m_elastic.setOperators();
// Projection on "x" field
if (!is_displacement_field) {
Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
}
// Projection on "X" field
{
Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
}
CHKERR MetaNeumannForces::addNeumannBCElements(m_field, primary_field_name);
// Add spring boundary condition applied on surfaces.
CHKERR MetaSpringBC::addSpringElements(m_field, primary_field_name,
"MESH_NODE_POSITIONS");
// build finite elemnts
auto bit_ref = m_contact.getBitRefLevel();
// build adjacencies
CHKERR m_field.build_adjacencies(bit_ref);
// define problems
CHKERR m_field.add_problem("CONTACT_PROB");
// set refinement level for problem
CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB", bit_ref);
DMType dm_name = "DMMOFEM";
SmartPetscObj<DM> dm;
dm = createSmartDM(m_field.get_comm(), dm_name);
// create dm instance
CHKERR DMSetType(dm, dm_name);
// set dm datastruture which created mofem datastructures
CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_ref);
CHKERR DMSetFromOptions(dm);
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
// add elements to dm
CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm, "SPRING");
CHKERR m_contact.addElementsToDM(dm);
CHKERR m_elastic.addElementsToDM(dm);
CHKERR DMSetUp(dm);
// Vector of DOFs and the RHS
auto D = smartCreateDMVector(dm);
// Stiffness matrix
auto Aij = smartCreateDMMatrix(dm);
CHKERR VecZeroEntries(D);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecZeroEntries(F);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
CHKERR MatZeroEntries(Aij);
// Dirichlet BC
boost::shared_ptr<FEMethod> dirichlet_bc_ptr;
if (!is_displacement_field)
dirichlet_bc_ptr =
boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
m_field, primary_field_name, Aij, D, F));
else
dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(
new DirichletDisplacementBc(m_field, primary_field_name, Aij, D, F));
dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
if (!is_displacement_field) {
dirichlet_bc_ptr->snes_x = D;
} else {
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);
}
// Assemble pressure and traction forces
boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
m_field, neumann_forces, NULL, primary_field_name);
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumann_forces.begin();
for (; mit != neumann_forces.end(); mit++) {
mit->second->methodsOp.push_back(new MortarContactProblem::LoadScale());
CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
&mit->second->getLoopFe(), NULL, NULL);
}
// Implementation of spring element
// Create new instances of face elements for springs
auto fe_spring_lhs_ptr =
boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
auto fe_spring_rhs_ptr =
boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, primary_field_name,
"MESH_NODE_POSITIONS");
CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
dirichlet_bc_ptr.get(), NULL);
CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
PETSC_NULL);
// setup contact methods
dirichlet_bc_ptr.get());
boost::shared_ptr<FEMethod> fe_null;
CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
fe_null);
CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, NULL, NULL);
// setup contact methods
dirichlet_bc_ptr);
if (test_num) {
char testing_options[] = "-ksp_type fgmres "
"-pc_type lu "
"-pc_factor_mat_solver_type mumps "
"-snes_type newtonls "
"-snes_linesearch_type basic "
"-snes_max_it 20 "
"-snes_atol 1e-8 "
"-snes_rtol 1e-8 ";
CHKERR PetscOptionsInsertString(NULL, testing_options);
}
auto snes = MoFEM::createSNES(m_field.get_comm());
CHKERR SNESSetDM(snes, dm);
SNESConvergedReason snes_reason;
SnesCtx *snes_ctx;
// create snes nonlinear solver
{
CHKERR SNESSetDM(snes, dm);
CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
}
for (int ss = 0; ss != nb_sub_steps; ++ss) {
MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
MOFEM_LOG_C("WORLD", Sev::inform, "Load lambda %6.4e\n",
CHKERR SNESSolve(snes, PETSC_NULL, D);
CHKERR SNESGetConvergedReason(snes, &snes_reason);
int its;
CHKERR SNESGetIterationNumber(snes, &its);
MOFEM_LOG_C("WORLD", Sev::inform, "number of Newton iterations = %D\n\n",
its);
// save on mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
MOFEM_LOG("WORLD", Sev::inform) << "Loop post proc\n";
CHKERR m_elastic.postProcessElement(ss);
auto common_data_post = m_contact.commonDataPostProc;
std::ofstream ofs((std::string("test_simple_contact") + ".txt").c_str());
m_contact.getPostProcSimpleContactPipeline().push_back(
m_field, primary_field_name, common_data_post, ofs));
CHKERR m_contact.postProcessElement(ss);
double energy = m_elastic.elasticElementPtr->getLoopFeEnergy().eNergy;
ofs << "Elastic energy: " << energy << endl;
ofs.flush();
ofs.close();
auto &nb_gauss_pts = m_contact.nbGaussPts;
auto &contact_area = m_contact.contactArea;
int expected_nb_gauss_pts;
double expected_energy, expected_contact_area;
MortarCtestFunctions::check_tests(ss, test_num, expected_energy,
expected_contact_area,
expected_nb_gauss_pts);
constexpr double eps = 1e-8;
if (m_field.get_comm_rank() == 0) {
if ((int)nb_gauss_pts[0] != expected_nb_gauss_pts) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong number of active gauss pts %d != %d",
expected_nb_gauss_pts, (int)nb_gauss_pts[0]);
}
if (std::abs(contact_area[0] - expected_contact_area) > eps) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong active contact area %6.4e != %6.4e",
expected_contact_area, contact_area[0]);
}
}
if (std::abs(energy - expected_energy) > eps)
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong energy %6.4e != %6.4e", expected_energy, energy);
}
}
// finish work cleaning memory, getting statistics, etc
return 0;
}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
const std::string default_options
std::string param_file
Elastic materials.
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
static const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:676
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
auto smartCreateDMMatrix
Get smart matrix from DM.
Definition: DMMoFEM.hpp:944
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:717
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1032
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:493
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:504
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
char mesh_file_name[255]
int main(int argc, char *argv[])
static char help[]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
auto createSNES
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:148
auto createSmartDM
Creates smart DM object.
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:39
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode check_tests(int ss, int test_num, double &expected_energy, double &expected_contact_area, int &expected_nb_gauss_pts)
const double D
diffusivity
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
Set Dirichlet boundary conditions on spatial displacements.
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode getCommandLineParameters() override
MoFEMErrorCode createElements() override
MoFEMErrorCode setupSolverJacobianSNES() override
boost::ptr_vector< MoFEM::ForcesAndSourcesCore::UserDataOperator > & getPostProcSimpleContactPipeline()
std::array< double, 2 > nbGaussPts
MoFEMErrorCode postProcessElement(int step) override
MoFEMErrorCode setOperators() override
MoFEMErrorCode addElementsToDM(SmartPetscObj< DM > dm) override
MoFEMErrorCode setupSolverFunctionSNES() override
BitRefLevel getBitRefLevel() override
boost::shared_ptr< MortarContactProblem::CommonDataMortarContact > commonDataPostProc
std::array< double, 2 > contactArea
MoFEMErrorCode addElementFields() override
Set of functions declaring elements and setting operators for generic element interface.
boost::shared_ptr< NonlinearElasticElement > elasticElementPtr
MoFEMErrorCode addElementsToDM(SmartPetscObj< DM > dm)