v0.13.1
Loading...
Searching...
No Matches
mortar_contact_thermal.cpp

Implementation of mortar contact between surfaces with non-matching meshes taking into account internal stress resulting from the thermal expansion

/** \file mortar_contact_thermal.cpp
* \example mortar_contact_thermal.cpp
*
* Implementation of mortar contact between surfaces with non-matching meshes
* taking into account internal stress resulting from the thermal expansion
*
**/
/* 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>
static char help[] = "\n";
using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide;
struct VolRule {
int operator()(int, int, int order) const { return 2 * (order - 1); }
};
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"
"-r_value 1. \n"
"-alm_flag 0 \n"
"-eigen_pos_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;
PetscBool flg_file_out;
char mesh_file_name[255];
char output_mesh_name[255];
PetscInt order = 1;
PetscInt order_contact = 1;
PetscInt nb_ho_levels = 0;
PetscInt order_lambda = 1;
PetscReal r_value = 1.;
PetscReal cn_value = -1;
PetscInt nb_sub_steps = 1;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool is_newton_cotes = PETSC_FALSE;
PetscInt test_num = 0;
PetscBool convect_pts = PETSC_FALSE;
PetscBool print_contact_state = PETSC_FALSE;
PetscBool alm_flag = PETSC_FALSE;
PetscBool eigen_pos_flag = PETSC_FALSE;
PetscReal thermal_expansion_coef = 1e-5;
PetscReal init_temp = 250.0;
PetscReal scale_factor = 1.0;
PetscBool analytical_input = PETSC_TRUE;
char stress_tag_name[255];
PetscBool flg_tag_name;
PetscBool save_mean_stress = PETSC_TRUE;
PetscBool ignore_contact = PETSC_FALSE;
PetscBool ignore_pressure = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsString("-output_mesh_file", "output mesh file name", "",
"mesh.h5m", output_mesh_name, 255, &flg_file_out);
CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
"", 1, &order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-order_contact",
"approximation order of spatial positions in contact interface", "", 1,
&order_contact, PETSC_NULL);
CHKERR PetscOptionsInt("-ho_levels_num", "number of higher order levels",
"", 0, &nb_ho_levels, PETSC_NULL);
CHKERR PetscOptionsInt("-order_lambda",
"approximation order of Lagrange multipliers", "", 1,
&order_lambda, 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 PetscOptionsReal("-cn_value", "default regularisation cn value", "",
1., &cn_value, PETSC_NULL);
CHKERR PetscOptionsBool("-is_newton_cotes",
"set if Newton-Cotes quadrature rules are used", "",
PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
PETSC_NULL);
CHKERR PetscOptionsBool("-convect", "set to convect integration pts", "",
PETSC_FALSE, &convect_pts, PETSC_NULL);
CHKERR PetscOptionsBool("-print_contact_state",
"output number of active gp at every iteration", "",
PETSC_FALSE, &print_contact_state, PETSC_NULL);
CHKERR PetscOptionsBool("-alm_flag",
"if set use ALM, if not use C-function", "",
PETSC_FALSE, &alm_flag, PETSC_NULL);
CHKERR PetscOptionsBool("-eigen_pos_flag",
"if set use eigen spatial positions are taken into "
"account for predeformed configuration",
"", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
CHKERR PetscOptionsReal("-scale_factor", "scale factor", "", scale_factor,
&scale_factor, PETSC_NULL);
CHKERR PetscOptionsBool("-ignore_contact", "if set true, ignore contact",
"", PETSC_FALSE, &ignore_contact, PETSC_NULL);
CHKERR PetscOptionsBool("-ignore_pressure", "if set true, ignore pressure",
"", PETSC_FALSE, &ignore_pressure, PETSC_NULL);
CHKERR PetscOptionsBool("-analytical_input",
"if set true, use analytical strain", "",
PETSC_TRUE, &analytical_input, PETSC_NULL);
CHKERR PetscOptionsBool("-save_mean_stress",
"if set true, save mean stress", "", PETSC_TRUE,
&save_mean_stress, PETSC_NULL);
CHKERR PetscOptionsString("-stress_tag_name", "stress tag name file name",
"", "INTERNAL_STRESS", stress_tag_name, 255,
&flg_tag_name);
CHKERR PetscOptionsReal(
"-thermal_expansion_coef", "thermal expansion coef ", "",
thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
CHKERR PetscOptionsReal("-init_temp", "init_temp ", "", init_temp,
&init_temp, 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)");
}
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);
Version file_ver;
MOFEM_LOG("WORLD", Sev::inform) << "File version " << file_ver.strVersion();
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
std::vector<BitRefLevel> bit_levels;
bit_levels.push_back(BitRefLevel().set(0));
auto bit_ref_manager = m_field.getInterface<BitRefManager>();
CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
auto prism_interface = m_field.getInterface<PrismInterface>();
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, ref_level_meshset);
CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
bit_levels.back(), BitRefLevel().set(), MBPRISM,
ref_level_meshset);
// get faces and tets to split
CHKERR prism_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 prism_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;
CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
true);
CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE, true);
CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
}
CHKERR bit_ref_manager->shiftRightBitRef(1);
bit_levels.pop_back();
}
}
};
auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
Range &simple_contact_prisms,
Range &simple_master_tris,
Range &simple_slave_tris) {
EntityHandle meshset_prisms;
CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
bit_levels.back(), BitRefLevel().set(), MBPRISM, meshset_prisms);
CHKERR moab.get_entities_by_handle(meshset_prisms, simple_contact_prisms);
CHKERR moab.delete_entities(&meshset_prisms, 1);
for (Range::iterator pit = simple_contact_prisms.begin();
pit != simple_contact_prisms.end(); pit++) {
CHKERR moab.side_element(*pit, 2, 3, tri);
simple_master_tris.insert(tri);
CHKERR moab.side_element(*pit, 2, 4, tri);
simple_slave_tris.insert(tri);
}
};
Range simple_contact_prisms, simple_master_tris, simple_slave_tris;
if (!ignore_contact) {
CHKERR add_prism_interface(bit_levels);
CHKERR find_contact_prisms(bit_levels, simple_contact_prisms,
simple_master_tris, simple_slave_tris);
}
Range mortar_contact_prisms, mortar_master_tris, mortar_slave_tris;
if (it->getName().compare(0, 13, "MORTAR_MASTER") == 0) {
CHKERR m_field.get_moab().get_entities_by_type(
it->meshset, MBTRI, mortar_master_tris, true);
}
}
if (it->getName().compare(0, 12, "MORTAR_SLAVE") == 0) {
CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI,
mortar_slave_tris, true);
}
}
ContactSearchKdTree contact_search_kd_tree(m_field);
boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
contact_commondata_multi_index;
contact_commondata_multi_index =
boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
CHKERR contact_search_kd_tree.buildTree(mortar_master_tris);
CHKERR contact_search_kd_tree.contactSearchAlgorithm(
mortar_master_tris, mortar_slave_tris, contact_commondata_multi_index,
mortar_contact_prisms);
EntityHandle meshset_slave_master_prisms;
CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
moab.add_entities(meshset_slave_master_prisms, mortar_contact_prisms);
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
meshset_slave_master_prisms, 3, bit_levels.back());
Range tris_from_prism;
// find actual masters and slave used
CHKERR m_field.get_moab().get_adjacencies(mortar_contact_prisms, 2, true,
tris_from_prism,
moab::Interface::UNION);
tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
mortar_slave_tris = intersect(tris_from_prism, mortar_slave_tris);
Range all_contact_prisms, all_slave_tris;
all_contact_prisms = unite(simple_contact_prisms, mortar_contact_prisms);
all_slave_tris = unite(simple_slave_tris, mortar_slave_tris);
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);
if (!eigen_pos_flag) {
CHKERR m_field.add_field("EIGEN_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);
if (!eigen_pos_flag) {
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_POSITIONS", 1);
}
CHKERR m_field.add_ents_to_field_by_type(all_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);
auto set_contact_order = [&](Range &contact_prisms, int order_contact,
int nb_ho_levels) {
Range contact_tris, contact_edges;
CHKERR moab.get_adjacencies(contact_prisms, 2, false, contact_tris,
moab::Interface::UNION);
contact_tris = contact_tris.subset_by_type(MBTRI);
CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
moab::Interface::UNION);
Range ho_ents;
ho_ents.merge(contact_tris);
ho_ents.merge(contact_edges);
for (int ll = 0; ll < nb_ho_levels; ll++) {
Range ents, verts, tets;
CHKERR moab.get_connectivity(ho_ents, verts, true);
CHKERR moab.get_adjacencies(verts, 3, false, tets,
moab::Interface::UNION);
tets = tets.subset_by_type(MBTET);
for (auto d : {1, 2}) {
CHKERR moab.get_adjacencies(tets, d, false, ents,
moab::Interface::UNION);
}
ho_ents = unite(ho_ents, ents);
ho_ents = unite(ho_ents, tets);
}
CHKERR m_field.set_field_order(ho_ents, "SPATIAL_POSITION",
order_contact);
};
if (!ignore_contact && order_contact > order) {
CHKERR set_contact_order(all_contact_prisms, order_contact, nb_ho_levels);
}
// 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);
}
Range slave_verts;
CHKERR moab.get_connectivity(all_slave_tris, slave_verts, true);
CHKERR m_field.getInterface<FieldBlas>()->setField(0.0, MBVERTEX,
slave_verts, "LAGMULT");
// 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);
boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
new VolumeElementForcesAndSourcesCore(m_field));
boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
new VolumeElementForcesAndSourcesCore(m_field));
fe_elastic_lhs_ptr->getRuleHook = VolRule();
fe_elastic_rhs_ptr->getRuleHook = VolRule();
CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS", false);
auto data_hooke_element_at_pts =
boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
block_sets_ptr, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS", false, false,
MBTET, data_hooke_element_at_pts);
auto thermal_strain =
[&thermal_expansion_coef, &init_temp, &test_num](
double temp;
t_thermal_strain(i, j) = 0.0;
switch (test_num) {
case 0:
// Put here analytical formula which may depend on coordinates
temp = init_temp - 1.0;
t_thermal_strain(i, j) =
thermal_expansion_coef * (temp - init_temp) * t_kd(i, j);
break;
case 1:
case 2:
t_thermal_strain(i, j) = -thermal_expansion_coef * t_kd(i, j);
break;
case 3:
t_thermal_strain(2, 2) = thermal_expansion_coef;
break;
case 4:
t_thermal_strain(i, j) = thermal_expansion_coef * t_kd(i, j);
break;
default:
break;
}
return t_thermal_strain;
};
if (analytical_input) {
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAnalyticalInternalStrain_dx<0>(
"SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
thermal_strain));
} else {
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
moab, stress_tag_name));
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", data_hooke_element_at_pts));
}
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
"MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
fe_elastic_rhs_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
*block_sets_ptr.get(), moab, scale_factor, save_mean_stress, false,
false));
Range all_tets;
if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
Range tets;
CHKERR moab.get_entities_by_dimension(bit->getMeshset(), 3, tets, true);
all_tets.merge(tets);
}
}
Skinner skinner(&moab);
Range skin_tris;
CHKERR skinner.find_skin(0, all_tets, false, skin_tris);
CHKERR m_field.add_finite_element("SKIN", MF_ZERO);
CHKERR m_field.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
"SPATIAL_POSITION");
"SPATIAL_POSITION");
"SPATIAL_POSITION");
"MESH_NODE_POSITIONS");
"EIGEN_POSITIONS");
auto make_mortar_contact_element = [&]() {
return boost::make_shared<MortarContactProblem::MortarContactElement>(
m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
};
auto make_convective_mortar_master_element = [&]() {
return boost::make_shared<
m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
};
auto make_convective_mortar_slave_element = [&]() {
return boost::make_shared<
m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
};
auto make_mortar_contact_common_data = [&]() {
return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
m_field);
};
auto get_mortar_contact_rhs = [&](auto contact_problem, auto make_element,
bool is_alm = false) {
auto fe_rhs_mortar_contact = make_element();
auto common_data_mortar_contact = make_mortar_contact_common_data();
if (print_contact_state) {
fe_rhs_mortar_contact->contactStateVec =
common_data_mortar_contact->gaussPtsStateVec;
}
contact_problem->setContactOperatorsRhs(
fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_rhs_mortar_contact;
};
auto get_mortar_master_traction_rhs = [&](auto contact_problem,
auto make_element,
bool is_alm = false) {
auto fe_rhs_mortar_contact = make_element();
auto common_data_mortar_contact = make_mortar_contact_common_data();
contact_problem->setMasterForceOperatorsRhs(
fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_rhs_mortar_contact;
};
auto get_mortar_master_traction_lhs = [&](auto contact_problem,
auto make_element,
bool is_alm = false) {
auto fe_lhs_mortar_contact = make_element();
auto common_data_mortar_contact = make_mortar_contact_common_data();
contact_problem->setMasterForceOperatorsLhs(
fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_lhs_mortar_contact;
};
auto get_mortar_master_help_traction_lhs =
[&](auto contact_problem, auto make_element, bool is_alm = false) {
auto fe_lhs_mortar_contact = make_element();
auto common_data_mortar_contact = make_mortar_contact_common_data();
contact_problem->setMasterForceOperatorsLhs(
fe_lhs_mortar_contact, common_data_mortar_contact,
"SPATIAL_POSITION", "LAGMULT", is_alm);
return fe_lhs_mortar_contact;
};
auto get_mortar_contact_lhs = [&](auto contact_problem, auto make_element,
bool is_alm = false) {
auto fe_lhs_mortar_contact = make_element();
auto common_data_mortar_contact = make_mortar_contact_common_data();
contact_problem->setContactOperatorsLhs(
fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_lhs_mortar_contact;
};
auto get_mortar_contact_help_lhs =
[&](auto contact_problem, auto make_element, bool is_alm = false) {
auto fe_lhs_mortar_contact = make_element();
auto common_data_mortar_contact = make_mortar_contact_common_data();
contact_problem->setContactOperatorsLhs(
fe_lhs_mortar_contact, common_data_mortar_contact,
"SPATIAL_POSITION", "LAGMULT", is_alm);
return fe_lhs_mortar_contact;
};
auto cn_value_ptr = boost::make_shared<double>(cn_value);
auto mortar_contact_problem = boost::make_shared<MortarContactProblem>(
m_field, contact_commondata_multi_index, cn_value_ptr, is_newton_cotes);
// add fields to the global matrix by adding the element
if (!eigen_pos_flag)
mortar_contact_problem->addMortarContactElement(
"MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
mortar_contact_prisms);
else
mortar_contact_problem->addMortarContactElement(
"MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
mortar_contact_prisms, "MESH_NODE_POSITIONS", eigen_pos_flag,
"EIGEN_POSITIONS");
auto make_simple_contact_element = [&]() {
return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
m_field);
};
auto make_convective_simple_master_element = [&]() {
return boost::make_shared<
m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
};
auto make_convective_simple_slave_element = [&]() {
return boost::make_shared<
m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
};
auto make_simple_contact_common_data = [&]() {
return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
m_field);
};
auto get_simple_contact_rhs = [&](auto contact_problem, auto make_element,
bool is_alm = false) {
auto fe_rhs_simple_contact = make_element();
auto common_data_simple_contact = make_simple_contact_common_data();
if (print_contact_state) {
fe_rhs_simple_contact->contactStateVec =
common_data_simple_contact->gaussPtsStateVec;
}
contact_problem->setContactOperatorsRhs(
fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_rhs_simple_contact;
};
auto get_simple_master_traction_rhs = [&](auto contact_problem,
auto make_element,
bool is_alm = false) {
auto fe_rhs_simple_contact = make_element();
auto common_data_simple_contact = make_simple_contact_common_data();
contact_problem->setMasterForceOperatorsRhs(
fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_rhs_simple_contact;
};
auto get_simple_master_traction_lhs =
[&](auto contact_problem, auto make_element, bool is_alm = false) {
auto fe_lhs_simple_contact = make_element();
auto common_data_simple_contact = make_simple_contact_common_data();
contact_problem->setMasterForceOperatorsLhs(
fe_lhs_simple_contact, common_data_simple_contact,
"SPATIAL_POSITION", "LAGMULT", is_alm);
return fe_lhs_simple_contact;
};
auto get_contact_lhs = [&](auto contact_problem, auto make_element,
bool is_alm = false) {
auto fe_lhs_simple_contact = make_element();
auto common_data_simple_contact = make_simple_contact_common_data();
contact_problem->setContactOperatorsLhs(
fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
return fe_lhs_simple_contact;
};
auto simple_contact_problem = boost::make_shared<SimpleContactProblem>(
m_field, cn_value_ptr, is_newton_cotes);
simple_contact_problem->addContactElement("SIMPLE_CONTACT_ELEM",
"SPATIAL_POSITION", "LAGMULT",
simple_contact_prisms);
simple_contact_problem->addPostProcContactElement(
"CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
"MESH_NODE_POSITIONS", all_slave_tris);
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", MF_ZERO);
// set refinement level for problem
CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
bit_levels.back());
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_levels.back());
CHKERR DMSetFromOptions(dm);
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
// add elements to dm
CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT_ELEM");
CHKERR DMMoFEMAddElement(dm, "SIMPLE_CONTACT_ELEM");
CHKERR DMMoFEMAddElement(dm, "ELASTIC");
CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm, "SPRING");
CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
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);
fe_elastic_rhs_ptr->snes_f = F;
fe_elastic_lhs_ptr->snes_B = 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++) {
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
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, "MORTAR_CONTACT_ELEM",
get_mortar_contact_rhs(mortar_contact_problem,
make_convective_mortar_master_element),
PETSC_NULL, PETSC_NULL);
dm, "MORTAR_CONTACT_ELEM",
get_mortar_master_traction_rhs(mortar_contact_problem,
make_convective_mortar_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
dm, "MORTAR_CONTACT_ELEM",
get_mortar_contact_rhs(mortar_contact_problem,
make_mortar_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
dm, "MORTAR_CONTACT_ELEM",
get_mortar_master_traction_rhs(mortar_contact_problem,
make_mortar_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
}
if (convect_pts == PETSC_TRUE) {
dm, "SIMPLE_CONTACT_ELEM",
get_simple_contact_rhs(simple_contact_problem,
make_convective_simple_master_element),
PETSC_NULL, PETSC_NULL);
dm, "SIMPLE_CONTACT_ELEM",
get_simple_master_traction_rhs(simple_contact_problem,
make_convective_simple_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
dm, "SIMPLE_CONTACT_ELEM",
get_simple_contact_rhs(simple_contact_problem,
make_simple_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
dm, "SIMPLE_CONTACT_ELEM",
get_simple_master_traction_rhs(simple_contact_problem,
make_simple_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
}
CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", fe_elastic_rhs_ptr, PETSC_NULL,
PETSC_NULL);
CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
PETSC_NULL);
CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, 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, "MORTAR_CONTACT_ELEM",
get_mortar_contact_help_lhs(mortar_contact_problem,
make_convective_mortar_master_element),
PETSC_NULL, PETSC_NULL);
dm, "MORTAR_CONTACT_ELEM",
get_mortar_master_help_traction_lhs(
mortar_contact_problem, make_convective_mortar_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
dm, "MORTAR_CONTACT_ELEM",
get_mortar_contact_lhs(mortar_contact_problem,
make_mortar_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
dm, "MORTAR_CONTACT_ELEM",
get_mortar_master_traction_lhs(mortar_contact_problem,
make_mortar_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
}
if (convect_pts == PETSC_TRUE) {
dm, "SIMPLE_CONTACT_ELEM",
get_contact_lhs(simple_contact_problem,
make_convective_simple_master_element),
PETSC_NULL, PETSC_NULL);
dm, "SIMPLE_CONTACT_ELEM",
get_simple_master_traction_lhs(simple_contact_problem,
make_convective_simple_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ELEM",
get_contact_lhs(simple_contact_problem,
make_simple_contact_element,
alm_flag),
PETSC_NULL, PETSC_NULL);
dm, "SIMPLE_CONTACT_ELEM",
get_simple_master_traction_lhs(simple_contact_problem,
make_simple_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
}
CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", fe_elastic_lhs_ptr, PETSC_NULL,
PETSC_NULL);
CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
PETSC_NULL);
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 10 "
"-snes_atol 1e-8 "
"-snes_rtol 1e-8 ";
CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
}
auto snes = MoFEM::createSNES(m_field.get_comm());
CHKERR SNESSetDM(snes, dm);
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);
}
/// Post proc on the skin
PostProcFaceOnRefinedMesh post_proc_skin(m_field);
CHKERR post_proc_skin.generateReferenceElementMesh();
CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false, false);
CHKERR post_proc_skin.addFieldValuesPostProc("SPATIAL_POSITION");
CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
CHKERR post_proc_skin.addFieldValuesPostProc("EIGEN_POSITIONS");
struct OpGetFieldGradientValuesOnSkin
: public FaceElementForcesAndSourcesCore::UserDataOperator {
const std::string feVolName;
boost::shared_ptr<VolSideFe> sideOpFe;
OpGetFieldGradientValuesOnSkin(const std::string field_name,
const std::string vol_fe_name,
boost::shared_ptr<VolSideFe> side_fe)
feVolName(vol_fe_name), sideOpFe(side_fe) {}
MoFEMErrorCode doWork(int side, EntityType type,
DataForcesAndSourcesCore::EntData &data) {
if (type != MBVERTEX)
CHKERR loopSideVolumes(feVolName, *sideOpFe);
}
};
boost::shared_ptr<VolSideFe> vol_side_fe_ptr =
boost::make_shared<VolSideFe>(m_field);
vol_side_fe_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
vol_side_fe_ptr->getOpPtrVector().push_back(
"MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
post_proc_skin.getOpPtrVector().push_back(
new OpGetFieldGradientValuesOnSkin("SPATIAL_POSITION", "ELASTIC",
vol_side_fe_ptr));
post_proc_skin.getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
"SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
post_proc_skin.getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
"MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
post_proc_skin.getOpPtrVector().push_back(
"SPATIAL_POSITION", data_hooke_element_at_pts,
*block_sets_ptr.get(), post_proc_skin.postProcMesh,
post_proc_skin.mapGaussPts, false, false));
for (int ss = 0; ss != nb_sub_steps; ++ss) {
if (!ignore_pressure) {
MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
} else {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Ignoring pressure...\n");
}
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
CHKERR SNESSolve(snes, PETSC_NULL, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
}
// save on mesh
CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
SmartPetscObj<Vec> v_energy;
CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS", false, false,
v_energy);
double energy;
CHKERR VecSum(v_energy, &energy);
// Print elastic energy
MOFEM_LOG_C("WORLD", Sev::inform, "Elastic energy: %8.8e", energy);
// moab_instance
moab::Core mb_post; // create database
moab::Interface &moab_proc = mb_post; // create interface to database
auto common_data_mortar_contact = make_mortar_contact_common_data();
boost::shared_ptr<MortarContactProblem::MortarContactElement>
fe_post_proc_mortar_contact;
if (convect_pts == PETSC_TRUE) {
fe_post_proc_mortar_contact = make_convective_mortar_master_element();
} else {
fe_post_proc_mortar_contact = make_mortar_contact_element();
}
std::array<double, 2> nb_gauss_pts;
std::array<double, 2> contact_area;
auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
CHKERR VecAssemblyBegin(vec);
CHKERR VecAssemblyEnd(vec);
const double *array;
CHKERR VecGetArrayRead(vec, &array);
if (m_field.get_comm_rank() == 0) {
for (int i : {0, 1})
data[i] = array[i];
}
CHKERR VecRestoreArrayRead(vec, &array);
};
mb_post.delete_mesh();
if (!mortar_contact_prisms.empty()) {
mortar_contact_problem->setContactOperatorsForPostProc(
fe_post_proc_mortar_contact, common_data_mortar_contact, m_field,
"SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
"EIGEN_POSITIONS");
CHKERR VecZeroEntries(common_data_mortar_contact->gaussPtsStateVec);
CHKERR VecZeroEntries(common_data_mortar_contact->contactAreaVec);
CHKERR DMoFEMLoopFiniteElements(dm, "MORTAR_CONTACT_ELEM",
fe_post_proc_mortar_contact);
CHKERR get_contact_data(common_data_mortar_contact->gaussPtsStateVec,
nb_gauss_pts);
CHKERR get_contact_data(common_data_mortar_contact->contactAreaVec,
contact_area);
if (m_field.get_comm_rank() == 0) {
PetscPrintf(PETSC_COMM_SELF,
"Active gauss pts (mortar): %d out of %d\n",
(int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
PetscPrintf(
PETSC_COMM_SELF,
"Active contact area (mortar): %8.8f out of %8.8f (%8.8f% %)\n",
contact_area[0], contact_area[1],
contact_area[0] / contact_area[1] * 100.);
}
}
auto common_data_simple_contact = make_simple_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_simple_master_element();
} else {
fe_post_proc_simple_contact = make_simple_contact_element();
}
if (!simple_contact_prisms.empty()) {
simple_contact_problem->setContactOperatorsForPostProc(
fe_post_proc_simple_contact, common_data_simple_contact, m_field,
"SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
"EIGEN_POSITIONS");
CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_CONTACT_ELEM",
fe_post_proc_simple_contact);
CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
nb_gauss_pts);
CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
contact_area);
if (m_field.get_comm_rank() == 0) {
MOFEM_LOG_C("SELF", Sev::inform,
"Active gauss pts (matching): %d out of %d",
(int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
"SELF", Sev::inform,
"Active contact area (matching): %8.8f out of %8.8f (%8.8f% %)",
contact_area[0], contact_area[1],
contact_area[0] / contact_area[1] * 100.);
}
}
if (!ignore_contact) {
string out_file_name;
std::ostringstream strm;
strm << "out_contact_integ_pts"
<< ".h5m";
out_file_name = strm.str();
MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s\n",
out_file_name.c_str());
CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
"PARALLEL=WRITE_PART");
}
boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *post_proc_contact_ptr, false,
false);
CHKERR post_proc_contact_ptr->addFieldValuesPostProc("LAGMULT");
CHKERR post_proc_contact_ptr->addFieldValuesPostProc("SPATIAL_POSITION");
CHKERR post_proc_contact_ptr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
if (!all_slave_tris.empty()) {
CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
post_proc_contact_ptr);
string out_file_name;
std::ostringstream stm;
stm << "out_contact_surface"
<< ".h5m";
out_file_name = stm.str();
MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
out_file_name.c_str());
CHKERR post_proc_contact_ptr->postProcMesh.write_file(
out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
}
CHKERR m_field.getInterface<FieldBlas>()->fieldCopy(1., "SPATIAL_POSITION",
"EIGEN_POSITIONS");
{
PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
CHKERR DMoFEMLoopFiniteElements(dm, "SKIN", &post_proc_skin);
ostringstream stm;
string out_file_name;
stm << "out_skin.h5m";
out_file_name = stm.str();
MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
out_file_name.c_str());
CHKERR post_proc_skin.writeFile(stm.str());
}
const int n_parts = m_field.get_comm_size();
if (m_field.get_comm_rank() == 0) {
dm, "ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
CHKERR moab.get_adjacencies(all_contact_prisms, 2, true, faces,
moab::Interface::UNION);
tris = faces.subset_by_type(MBTRI);
quads = faces.subset_by_type(MBQUAD);
CHKERR moab.get_adjacencies(tris, 1, true, tris_edges,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(quads, 1, true, quads_edges,
moab::Interface::UNION);
ents_to_delete.merge(all_contact_prisms);
ents_to_delete.merge(quads);
ents_to_delete.merge(subtract(quads_edges, tris_edges));
CHKERR moab.delete_entities(ents_to_delete);
MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s", output_mesh_name);
CHKERR moab.write_file(output_mesh_name, "MOAB");
auto get_tag_handle = [&](auto name, auto size) {
Tag th;
std::vector<double> def_vals(size, 0.0);
CHKERR moab.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE,
def_vals.data());
return th;
};
if (test_num) {
Range tets;
CHKERR moab.get_entities_by_dimension(0, 3, tets);
EntityHandle tet = tets.front();
std::array<double, 9> internal_stress, actual_stress;
std::array<double, 9> internal_stress_ref, actual_stress_ref;
std::array<double, 2> nb_gauss_pts_ref, contact_area_ref;
switch (test_num) {
case 3:
actual_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
if (strcmp(stress_tag_name, "INTERNAL_STRESS") == 0)
internal_stress_ref = {0., 0., -200., 0., 0., 0., 0., 0., 0.};
else
internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
break;
case 4:
nb_gauss_pts_ref = {216, 492};
contact_area_ref = {0.125, 0.25};
break;
default:
SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Test number %d not found",
test_num);
}
auto th_internal_stress = get_tag_handle("MED_INTERNAL_STRESS", 9);
auto th_actual_stress = get_tag_handle("MED_ACTUAL_STRESS", 9);
CHKERR moab.tag_get_data(th_internal_stress, &tet, 1,
internal_stress.data());
CHKERR moab.tag_get_data(th_actual_stress, &tet, 1,
actual_stress.data());
if (test_num == 4) {
const double eps = 1e-3;
if (std::abs(nb_gauss_pts_ref[0] - nb_gauss_pts[0]) > eps) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong number of active contact gauss pts: should be %d "
"but is %d",
(int)nb_gauss_pts_ref[0], (int)nb_gauss_pts[0]);
}
if (std::abs(nb_gauss_pts_ref[1] - nb_gauss_pts[1]) > eps) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong total number of contact gauss pts: should be %d "
"but is %d",
(int)nb_gauss_pts_ref[1], (int)nb_gauss_pts[1]);
}
if (std::abs(contact_area_ref[0] - contact_area[0]) > eps) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong active contact area: should be %g "
"but is %g",
contact_area_ref[0], contact_area[0]);
}
if (std::abs(contact_area_ref[1] - contact_area[1]) > eps) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong potential contact area: should be %g "
"but is %g",
contact_area_ref[1], contact_area[1]);
}
} else {
if (save_mean_stress) {
const double eps = 1e-10;
for (int i = 0; i < 9; i++) {
if (std::abs(internal_stress[i] - internal_stress_ref[i]) > eps) {
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong component %d of internal stress: should be %g "
"but is %g",
i, internal_stress_ref[i], internal_stress[i]);
}
if (std::abs(actual_stress[i] - actual_stress_ref[i]) > eps) {
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong component %d of actual stress: should be %g "
"but is %g",
i, actual_stress_ref[i], actual_stress[i]);
}
}
}
}
}
}
}
// finish work cleaning memory, getting statistics, etc
return 0;
}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
const std::string default_options
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static PetscErrorCode ierr
static const double eps
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:953
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1070
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:665
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:118
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:450
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:470
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:706
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1041
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:482
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:514
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:493
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
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.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
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
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
char out_file_name[255]
double D
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
VolumeElementForcesAndSourcesCoreOnSide VolSideFe
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
auto createSNES(MPI_Comm comm)
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
constexpr auto field_name
multi_index_container< boost::shared_ptr< ContactCommonData >, indexed_by< ordered_unique< tag< Prism_tag >, member< ContactCommonData, EntityHandle, &ContactCommonData::pRism > > > > ContactCommonData_multiIndex
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 int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Element used to integrate on master surfaces. It convects integration points on slaves,...
data for calculation heat conductivity and heat capacity elements
Postprocess on face.
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Element used to integrate on master surfaces. It convects integration points on slaves,...
Set integration rule.
int operator()(int, int, int) const