v0.14.0
Loading...
Searching...
No Matches
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 mesh",
"", 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");
// 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;
auto simple = m_field.getInterface<Simple>();
simple->loadFile("", mesh_file_name);
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);
m_field, primary_field_name, "MESH_NODE_POSITIONS", "CONTACT_PROB",
"ELASTIC", is_displacement_field, true,
CHKERR m_boundary.addElementFields();
CHKERR m_contact.addElementFields();
CHKERR m_elastic.addElementFields();
// build field
CHKERR m_field.build_fields();
CHKERR m_boundary.createElements();
CHKERR m_contact.createElements();
CHKERR m_elastic.createElements();
// 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);
}
// 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";
CHKERR DMRegister_MoFEM(dm_name);
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 m_boundary.setOperators();
CHKERR m_contact.setOperators();
CHKERR m_elastic.setOperators();
CHKERR m_boundary.addElementsToDM(dm);
CHKERR m_contact.addElementsToDM(dm);
CHKERR m_elastic.addElementsToDM(dm);
CHKERR DMSetUp(dm);
// #ifndef NDEBUG
// cout << "<<< checkMPIAIJWithArraysMatrixFillIn <<< " << endl;
// CHKERR m_field.getInterface<MatrixManager>()
// ->checkMPIAIJWithArraysMatrixFillIn<PetscGlobalIdx_mi_tag>(
// "CONTACT_PROB", -2, -2, 0);
// #endif // NDEBUG
auto D = smartCreateDMVector(dm);
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);
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",
CHKERR SNESSolve(snes, PETSC_NULL, D);
// 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);
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);
CHKERR m_boundary.postProcessElement(ss);
double energy = m_elastic.elasticElementPtr->getLoopFeEnergy().eNergy;
ofs << "Elastic energy: " << energy << endl;
ofs.flush();
ofs.close();
if (test_num) {
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-6;
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;
}
const std::string default_options
std::string param_file
Elastic materials.
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
static char help[]
int main()
static PetscErrorCode ierr
static const double eps
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define CHKERR
Inline error check.
constexpr int order
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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
double D
char mesh_file_name[255]
auto createSNES(MPI_Comm comm)
MoFEMErrorCode check_tests(int ss, int test_num, double &expected_energy, double &expected_contact_area, int &expected_nb_gauss_pts)
Set of functions declaring elements and setting operators for basic boundary conditions interface.
MoFEMErrorCode addElementsToDM(SmartPetscObj< DM > dm) override
MoFEMErrorCode postProcessElement(int step) override
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:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode getCommandLineParameters() override
MoFEMErrorCode createElements() override
MoFEMErrorCode setupSolverJacobianSNES() override
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::ptr_deque< MoFEM::ForcesAndSourcesCore::UserDataOperator > & getPostProcSimpleContactPipeline()
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)