v0.9.1
simple_contact.cpp

Implementation of simple contact between surfaces with matching meshes

/** \file simple_contact.cpp
* \example simple_contact.cpp
*
* Implementation of simple contact between surfaces with 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 <Hooke.hpp>
using namespace std;
using namespace MoFEM;
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_package mumps \n"
"-snes_type newtonls \n"
"-snes_linesearch_type basic \n"
"-snes_max_it 10 \n"
"-snes_atol 1e-8 \n"
"-snes_rtol 1e-8 \n"
"-my_order 1 \n"
"-my_order_lambda 1 \n"
"-my_cn_value 1. \n"
"-my_is_newton_cotes 0 \n"
"-my_is_test 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 order_lambda = 1;
PetscReal r_value = 1.;
PetscReal cn_value = -1;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool is_newton_cotes = PETSC_FALSE;
PetscBool is_test = PETSC_FALSE;
PetscBool convect_pts = PETSC_FALSE;
CHKERR 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", "", 1,
&order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-my_order_lambda",
"default approximation order of Lagrange multipliers", "", 1,
&order_lambda, PETSC_NULL);
CHKERR PetscOptionsBool("-my_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 PetscOptionsReal("-my_cn_value", "default regularisation cn value",
"", 1., &cn_value, PETSC_NULL);
CHKERR PetscOptionsBool("-my_is_newton_cotes",
"set if Newton-Cotes quadrature rules are used", "",
PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
CHKERR PetscOptionsBool("-my_is_test", "set if run as test", "",
PETSC_FALSE, &is_test, PETSC_NULL);
CHKERR PetscOptionsBool("-my_convect", "set to convect integration pts", "",
PETSC_FALSE, &convect_pts, PETSC_NULL);
ierr = PetscOptionsEnd();
// Check if mesh file was provided
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
if (is_partitioned == PETSC_TRUE) {
// Read mesh to MOAB
const char *option;
option = "PARALLEL=BCAST_DELETE;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR moab.load_file(mesh_file_name, 0, option);
} else {
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();
auto add_prism_interface = [&](Range &contact_prisms, Range &master_tris,
Range &slave_tris,
std::vector<BitRefLevel> &bit_levels) {
PrismInterface *interface;
CHKERR m_field.getInterface(interface);
if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
cit->getName().c_str(), cit->getMeshsetId());
EntityHandle cubit_meshset = cit->getMeshset();
// get tet entities from back bit_level
EntityHandle ref_level_meshset;
CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
ref_level_meshset);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(), MBTET,
ref_level_meshset);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(), MBPRISM,
ref_level_meshset);
// get faces and tets to split
CHKERR interface->getSides(cubit_meshset, bit_levels.back(), true, 0);
// set new bit level
bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
// split faces and tets
CHKERR interface->splitSides(ref_level_meshset, bit_levels.back(),
cubit_meshset, true, true, 0);
// clean meshsets
CHKERR moab.delete_entities(&ref_level_meshset, 1);
// update cubit meshsets
for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
EntityHandle cubit_meshset = ciit->meshset;
->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
true);
->updateMeshsetByEntitiesChildren(cubit_meshset,
bit_levels.back(),
cubit_meshset, MBEDGE, true);
->updateMeshsetByEntitiesChildren(cubit_meshset,
bit_levels.back(),
cubit_meshset, MBTRI, true);
->updateMeshsetByEntitiesChildren(cubit_meshset,
bit_levels.back(),
cubit_meshset, MBTET, true);
}
CHKERR m_field.getInterface<BitRefManager>()->shiftRightBitRef(1);
bit_levels.pop_back();
}
}
EntityHandle meshset_prisms;
CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
->getEntitiesByTypeAndRefLevel(bit_levels.back(), BitRefLevel().set(),
MBPRISM, meshset_prisms);
CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
CHKERR moab.delete_entities(&meshset_prisms, 1);
for (Range::iterator pit = contact_prisms.begin();
pit != contact_prisms.end(); pit++) {
CHKERR moab.side_element(*pit, 2, 3, tri);
master_tris.insert(tri);
CHKERR moab.side_element(*pit, 2, 4, tri);
slave_tris.insert(tri);
}
};
Range contact_prisms, master_tris, slave_tris;
std::vector<BitRefLevel> bit_levels;
bit_levels.push_back(BitRefLevel().set(0));
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 3, bit_levels.back());
CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
bit_levels);
CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
MB_TAG_SPARSE, MF_ZERO);
CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
3, MB_TAG_SPARSE, MF_ZERO);
CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
MB_TAG_SPARSE, MF_ZERO);
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
// Declare problem add entities (by tets) to the field
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
CHKERR m_field.add_ents_to_field_by_type(slave_tris, MBTRI, "LAGMULT");
CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
// build field
CHKERR m_field.build_fields();
// Projection on "x" 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);
}
// Add elastic element
boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(new Hooke<adouble>());
boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>());
NonlinearElasticElement elastic(m_field, 2);
CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
CHKERR elastic.addElement("ELASTIC", "SPATIAL_POSITION");
CHKERR elastic.setOperators("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
false, false);
if (true) {
CHKERR m_field.add_finite_element("CONTACT_VTK");
"SPATIAL_POSITION");
"SPATIAL_POSITION");
"SPATIAL_POSITION");
"LAGMULT");
"LAGMULT");
"LAGMULT");
CHKERR m_field.add_ents_to_finite_element_by_type(slave_tris, MBTRI,
"CONTACT_VTK");
}
auto make_contact_element = [&]() {
return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
m_field);
};
auto make_convective_master_element = [&]() {
return boost::make_shared<
m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
};
auto make_convective_slave_element = [&]() {
return boost::make_shared<
m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
};
auto make_contact_common_data = [&]() {
return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
m_field);
};
auto get_contact_rhs = [&](auto contact_problem, auto make_element) {
auto fe_rhs_simple_contact = make_element();
auto common_data_simple_contact = make_contact_common_data();
contact_problem->setContactOperatorsRhs(fe_rhs_simple_contact,
common_data_simple_contact,
"SPATIAL_POSITION", "LAGMULT");
return fe_rhs_simple_contact;
};
auto get_master_traction_rhs = [&](auto contact_problem,
auto make_element) {
auto fe_rhs_simple_contact = make_element();
auto common_data_simple_contact = make_contact_common_data();
contact_problem->setMasterForceOperatorsRhs(
fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT");
return fe_rhs_simple_contact;
};
auto get_master_traction_lhs = [&](auto contact_problem,
auto make_element) {
auto fe_lhs_simple_contact = make_element();
auto common_data_simple_contact = make_contact_common_data();
contact_problem->setMasterForceOperatorsLhs(
fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT");
return fe_lhs_simple_contact;
};
auto get_master_contact_lhs = [&](auto contact_problem, auto make_element) {
auto fe_lhs_simple_contact = make_element();
auto common_data_simple_contact = make_contact_common_data();
contact_problem->setContactOperatorsLhs(fe_lhs_simple_contact,
common_data_simple_contact,
"SPATIAL_POSITION", "LAGMULT");
return fe_lhs_simple_contact;
};
auto contact_problem = boost::make_shared<SimpleContactProblem>(
m_field, cn_value, is_newton_cotes);
// add fields to the global matrix by adding the element
contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
"LAGMULT", contact_prisms);
CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "SPATIAL_POSITION");
// Add spring boundary condition applied on surfaces.
CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
// build finite elemnts
// build adjacencies
CHKERR m_field.build_adjacencies(bit_levels.back());
// 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_levels.back());
DMType dm_name = "DMMOFEM";
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_levels.back());
CHKERR DMSetFromOptions(dm);
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
// add elements to dm
CHKERR DMMoFEMAddElement(dm, "CONTACT_ELEM");
CHKERR DMMoFEMAddElement(dm, "ELASTIC");
CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm, "SPRING");
if (true) {
CHKERR DMMoFEMAddElement(dm, "CONTACT_VTK");
}
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 =
boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
m_field, "SPATIAL_POSITION", Aij, D, F));
dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
dirichlet_bc_ptr->snes_x = D;
// Assemble pressure and traction forces
boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
m_field, neumann_forces, NULL, "SPATIAL_POSITION");
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumann_forces.begin();
for (; mit != neumann_forces.end(); mit++) {
CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
&mit->second->getLoopFe(), NULL, NULL);
}
// Implementation of spring element
// Create new instances of face elements for springs
boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
"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);
if (convect_pts == PETSC_TRUE) {
dm, "CONTACT_ELEM",
get_contact_rhs(contact_problem, make_convective_master_element),
PETSC_NULL, PETSC_NULL);
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem,
make_convective_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
dm, "CONTACT_ELEM",
get_contact_rhs(contact_problem, make_contact_element), PETSC_NULL,
PETSC_NULL);
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem, make_contact_element),
PETSC_NULL, PETSC_NULL);
}
CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", &elastic.getLoopFeRhs(),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
PETSC_NULL);
dirichlet_bc_ptr.get());
boost::shared_ptr<FEMethod> fe_null;
CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
fe_null);
if (convect_pts == PETSC_TRUE) {
dm, "CONTACT_ELEM",
get_master_contact_lhs(contact_problem,
make_convective_master_element),
NULL, NULL);
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem,
make_convective_slave_element),
NULL, NULL);
} else {
dm, "CONTACT_ELEM",
get_master_contact_lhs(contact_problem, make_contact_element), NULL,
NULL);
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem, make_contact_element), NULL,
NULL);
}
CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", &elastic.getLoopFeLhs(), NULL,
NULL);
CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, NULL, NULL);
dirichlet_bc_ptr);
if (is_test == PETSC_TRUE) {
char testing_options[] = "-ksp_type fgmres "
"-pc_type lu "
"-pc_factor_mat_solver_package mumps "
"-snes_type newtonls "
"-snes_linesearch_type basic "
"-snes_max_it 10 "
"-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);
}
PostProcVolumeOnRefinedMesh post_proc(m_field);
// Add operators to the elements, starting with some generic
CHKERR post_proc.addFieldValuesPostProc("SPATIAL_POSITION");
CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
CHKERR post_proc.addFieldValuesGradientPostProc("SPATIAL_POSITION");
std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
elastic.setOfBlocks.begin();
for (; sit != elastic.setOfBlocks.end(); sit++) {
post_proc.getOpPtrVector().push_back(new PostProcStress(
post_proc.postProcMesh, post_proc.mapGaussPts, "SPATIAL_POSITION",
sit->second, post_proc.commonData));
}
CHKERR SNESSolve(snes, PETSC_NULL, D);
CHKERR SNESGetConvergedReason(snes, &snes_reason);
int its;
CHKERR SNESGetIterationNumber(snes, &its);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "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);
PetscPrintf(PETSC_COMM_WORLD, "Loop post proc\n");
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
elastic.getLoopFeEnergy().eNergy = 0;
PetscPrintf(PETSC_COMM_WORLD, "Loop energy\n");
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeEnergy());
// Print elastic energy
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
elastic.getLoopFeEnergy().eNergy);
string out_file_name;
std::ostringstream stm;
stm << "out"
<< ".h5m";
out_file_name = stm.str();
PetscPrintf(PETSC_COMM_WORLD, "out file %s\n", out_file_name.c_str());
CHKERR post_proc.postProcMesh.write_file(out_file_name.c_str(), "MOAB",
"PARALLEL=WRITE_PART");
// moab_instance
moab::Core mb_post; // create database
moab::Interface &moab_proc = mb_post; // create interface to database
auto common_data_simple_contact = make_contact_common_data();
boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
fe_post_proc_simple_contact;
if (convect_pts == PETSC_TRUE) {
fe_post_proc_simple_contact = make_convective_master_element();
} else {
fe_post_proc_simple_contact = make_contact_element();
}
contact_problem->setContactOperatorsForPostProc(
fe_post_proc_simple_contact, common_data_simple_contact, m_field,
"SPATIAL_POSITION", "LAGMULT", mb_post);
mb_post.delete_mesh();
if (is_test == PETSC_TRUE) {
std::ofstream ofs((std ::string("test_simple_contact") + ".txt").c_str());
fe_post_proc_simple_contact->getOpPtrVector().push_back(
m_field, "SPATIAL_POSITION", common_data_simple_contact, ofs));
CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_ELEM",
fe_post_proc_simple_contact);
ofs << "Elastic energy: " << elastic.getLoopFeEnergy().eNergy << endl;
ofs.flush();
ofs.close();
} else {
CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_ELEM",
fe_post_proc_simple_contact);
}
std::ostringstream ostrm;
ostrm << "out_contact"
<< ".h5m";
out_file_name = ostrm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
out_file_name.c_str());
CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
"PARALLEL=WRITE_PART");
if (true) {
boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
auto common_post_proc_data_simple_contact = make_contact_common_data();
CHKERR contact_problem->setPostProcContactOperators(
post_proc_contact_ptr, "SPATIAL_POSITION", "LAGMULT",
common_post_proc_data_simple_contact);
CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_VTK", post_proc_contact_ptr);
std::ostringstream stm;
std::string file_name_for_lagrange = "out_lagrange_for_vtk_";
stm << "file_name_for_lagrange"
<< ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
out_file_name.c_str());
CHKERR post_proc_contact_ptr->postProcMesh.write_file(
out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
}
EntityHandle out_meshset_slave_tris;
EntityHandle out_meshset_master_tris;
CHKERR moab.create_meshset(MESHSET_SET, out_meshset_slave_tris);
CHKERR moab.create_meshset(MESHSET_SET, out_meshset_master_tris);
CHKERR moab.add_entities(out_meshset_slave_tris, slave_tris);
CHKERR moab.add_entities(out_meshset_master_tris, master_tris);
CHKERR moab.write_file("out_slave_tris.vtk", "VTK", "",
&out_meshset_slave_tris, 1);
CHKERR moab.write_file("out_master_tris.vtk", "VTK", "",
&out_meshset_master_tris, 1);
CHKERR moab.delete_entities(&out_meshset_slave_tris, 1);
CHKERR moab.delete_entities(&out_meshset_master_tris, 1);
}
// finish work cleaning memory, getting statistics, etc
return 0;
}