v0.14.0
simple_contact_thermal.cpp

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

/** \file simple_contact_thermal.cpp
* \example simple_contact_thermal.cpp
*
* Implementation of simple contact between surfaces with 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>
using namespace MoFEM;
using namespace std;
static char help[] = "\n";
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"
"-my_order 1 \n"
"-my_order_lambda 1 \n"
"-my_order_contact 2 \n"
"-my_ho_levels_num 1 \n"
"-my_step_num 1 \n"
"-my_cn_value 1. \n"
"-my_r_value 1. \n"
"-my_alm_flag 0 \n"
"-my_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;
PetscBool deform_flat_flag = PETSC_FALSE;
PetscReal flat_shift = 1.0;
PetscReal mesh_height = 1.0;
PetscBool wave_surf_flag = PETSC_FALSE;
PetscInt wave_dim = 2;
PetscReal wave_length = 1.0;
PetscReal wave_ampl = 0.01;
PetscBool delete_prisms = 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 PetscOptionsString("-my_output_mesh_file", "output mesh file name",
"", "mesh.h5m", output_mesh_name, 255,
&flg_file_out);
CHKERR PetscOptionsInt("-my_order",
"approximation order of spatial positions", "", 1,
&order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-my_order_contact",
"approximation order of spatial positions in contact interface", "", 1,
&order_contact, PETSC_NULL);
CHKERR PetscOptionsInt("-my_ho_levels_num", "number of higher order levels",
"", 0, &nb_ho_levels, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_lambda",
"approximation order of Lagrange multipliers", "", 1,
&order_lambda, PETSC_NULL);
CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", nb_sub_steps,
&nb_sub_steps, 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 PetscOptionsInt("-my_test_num", "test number", "", 0, &test_num,
PETSC_NULL);
CHKERR PetscOptionsBool("-my_convect", "set to convect integration pts", "",
PETSC_FALSE, &convect_pts, PETSC_NULL);
CHKERR PetscOptionsBool("-my_print_contact_state",
"output number of active gp at every iteration", "",
PETSC_FALSE, &print_contact_state, PETSC_NULL);
CHKERR PetscOptionsBool("-my_alm_flag",
"if set use ALM, if not use C-function", "",
PETSC_FALSE, &alm_flag, PETSC_NULL);
CHKERR PetscOptionsBool("-my_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("-my_scale_factor", "scale factor", "",
scale_factor, &scale_factor, PETSC_NULL);
CHKERR PetscOptionsBool("-my_ignore_contact", "if set true, ignore contact",
"", PETSC_FALSE, &ignore_contact, PETSC_NULL);
CHKERR PetscOptionsBool("-my_ignore_pressure",
"if set true, ignore pressure", "", PETSC_FALSE,
&ignore_pressure, PETSC_NULL);
CHKERR PetscOptionsBool("-my_analytical_input",
"if set true, use analytical strain", "",
PETSC_TRUE, &analytical_input, PETSC_NULL);
CHKERR PetscOptionsBool("-my_save_mean_stress",
"if set true, save mean stress", "", PETSC_TRUE,
&save_mean_stress, PETSC_NULL);
CHKERR PetscOptionsString(
"-my_stress_tag_name", "stress tag name file name", "",
"INTERNAL_STRESS", stress_tag_name, 255, &flg_tag_name);
CHKERR PetscOptionsReal(
"-my_thermal_expansion_coef", "thermal expansion coef ", "",
thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
CHKERR PetscOptionsReal("-my_init_temp", "init_temp ", "", init_temp,
&init_temp, PETSC_NULL);
CHKERR PetscOptionsReal("-my_mesh_height",
"vertical dimension of the mesh ", "", mesh_height,
&mesh_height, PETSC_NULL);
CHKERR PetscOptionsBool("-my_deform_flat", "if set true, deform flat", "",
PETSC_FALSE, &deform_flat_flag, PETSC_NULL);
CHKERR PetscOptionsReal("-my_flat_shift", "flat shift ", "", flat_shift,
&flat_shift, PETSC_NULL);
CHKERR PetscOptionsBool("-my_wave_surf",
"if set true, make one of the surfaces wavy", "",
PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
CHKERR PetscOptionsInt("-my_wave_dim", "dimension (2 or 3)", "", wave_dim,
&wave_dim, PETSC_NULL);
CHKERR PetscOptionsReal("-my_wave_length", "profile wavelength", "",
wave_length, &wave_length, PETSC_NULL);
CHKERR PetscOptionsReal("-my_wave_ampl", "profile amplitude", "", wave_ampl,
&wave_ampl, PETSC_NULL);
CHKERR PetscOptionsBool("-my_delete_prisms", "if set true, delete prisms",
"", PETSC_FALSE, &delete_prisms, PETSC_NULL);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
// 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) {
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);
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
Version file_ver;
MOFEM_LOG("WORLD", Sev::inform) << "File version " << file_ver.strVersion();
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 &contact_prisms, Range &master_tris,
Range &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, 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;
if (!ignore_contact) {
if (analytical_input) {
CHKERR add_prism_interface(bit_levels);
}
CHKERR find_contact_prisms(bit_levels, contact_prisms, master_tris,
slave_tris);
}
auto deform_flat_surface = [&](int block_id, double shift, double height) {
Range all_tets, all_nodes;
if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
const int id = bit->getMeshsetId();
Range tets;
if (id == block_id) {
CHKERR m_field.get_moab().get_entities_by_dimension(
bit->getMeshset(), 3, tets, true);
all_tets.merge(tets);
}
}
}
CHKERR m_field.get_moab().get_connectivity(all_tets, all_nodes);
double coords[3];
for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
nit++) {
CHKERR moab.get_coords(&*nit, 1, coords);
double x = coords[0];
double y = coords[1];
double z = coords[2];
coords[2] -= shift;
CHKERR moab.set_coords(&*nit, 1, coords);
}
};
auto make_wavy_surface = [&](int block_id, int dim, double lambda,
double delta, double height) {
Range all_tets, all_nodes;
if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
const int id = bit->getMeshsetId();
Range tets;
if (id == block_id) {
CHKERR m_field.get_moab().get_entities_by_dimension(
bit->getMeshset(), 3, tets, true);
all_tets.merge(tets);
}
}
}
CHKERR m_field.get_moab().get_connectivity(all_tets, all_nodes);
double coords[3];
for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
nit++) {
CHKERR moab.get_coords(&*nit, 1, coords);
double x = coords[0];
double y = coords[1];
double z = coords[2];
double coef = (height + z) / height;
switch (dim) {
case 2:
coords[2] -= coef * delta * (1. - cos(2. * M_PI * x / lambda));
break;
case 3:
coords[2] -=
coef * delta *
(1. - cos(2. * M_PI * x / lambda) * cos(2. * M_PI * y / lambda));
break;
default:
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong dimension = %d", dim);
}
CHKERR moab.set_coords(&*nit, 1, coords);
}
};
if (deform_flat_flag && analytical_input) {
CHKERR deform_flat_surface(1, flat_shift, mesh_height);
CHKERR deform_flat_surface(2, -flat_shift, mesh_height);
}
if (wave_surf_flag && analytical_input) {
CHKERR make_wavy_surface(1, wave_dim, wave_length, wave_ampl,
mesh_height);
// CHKERR make_wavy_surface(2, wave_dim, wave_length, -wave_ampl,
// mesh_height);
}
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);
// 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(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);
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(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(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(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(
boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
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_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,
bool is_alm = false) {
auto fe_rhs_simple_contact = make_element();
auto common_data_simple_contact = make_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", true);
return fe_rhs_simple_contact;
};
auto get_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_contact_common_data();
contact_problem->setMasterForceOperatorsRhs(
fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
return fe_rhs_simple_contact;
};
auto get_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_contact_common_data();
contact_problem->setMasterForceOperatorsLhs(
fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
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_contact_common_data();
contact_problem->setContactOperatorsLhs(
fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
return fe_lhs_simple_contact;
};
auto cn_value_ptr = boost::make_shared<double>(cn_value);
auto contact_problem = boost::make_shared<SimpleContactProblem>(
m_field, cn_value_ptr, is_newton_cotes);
// add fields to the global matrix by adding the element
if (!eigen_pos_flag)
contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
"LAGMULT", contact_prisms);
else
contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
"LAGMULT", contact_prisms,
eigen_pos_flag, "EIGEN_POSITIONS");
contact_problem->addPostProcContactElement(
"CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
"MESH_NODE_POSITIONS", 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";
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");
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 SimpleContactProblem::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, "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, alm_flag),
PETSC_NULL, PETSC_NULL);
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem, make_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, "CONTACT_ELEM",
get_contact_lhs(contact_problem, make_convective_master_element),
PETSC_NULL, PETSC_NULL);
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem,
make_convective_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
dm, "CONTACT_ELEM",
get_contact_lhs(contact_problem, make_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem, make_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
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,
if (type != MBVERTEX)
CHKERR loopSideVolumes(feVolName, *sideOpFe);
}
};
boost::shared_ptr<VolSideFe> my_vol_side_fe_ptr =
boost::make_shared<VolSideFe>(m_field);
my_vol_side_fe_ptr->getOpPtrVector().push_back(
"SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
my_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",
my_vol_side_fe_ptr));
post_proc_skin.getOpPtrVector().push_back(
"SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
post_proc_skin.getOpPtrVector().push_back(
"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) {
SimpleContactProblem::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);
CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "SPATIAL_POSITION",
"MESH_NODE_POSITIONS", false, false,
v_energy);
const double *eng_ptr;
CHKERR VecGetArrayRead(v_energy, &eng_ptr);
// Print elastic energy
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy: %8.8e\n", *eng_ptr);
{
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();
PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n", out_file_name.c_str());
CHKERR post_proc_skin.writeFile(stm.str());
}
// 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();
}
std::array<double, 2> nb_gauss_pts;
std::array<double, 2> contact_area;
if (!ignore_contact) {
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", true);
mb_post.delete_mesh();
CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_ELEM",
fe_post_proc_simple_contact);
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);
};
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) {
PetscPrintf(PETSC_COMM_SELF, "Active gauss pts: %d out of %d\n",
(int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
PetscPrintf(PETSC_COMM_SELF,
"Active contact area: %8.8f out of %8.8f (%8.8f% %)\n",
contact_area[0], contact_area[1],
contact_area[0] / contact_area[1] * 100.);
}
string out_file_name;
std::ostringstream strm;
strm << "out_contact_integ_pts"
<< ".h5m";
out_file_name = strm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "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 (!ignore_contact) {
CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
post_proc_contact_ptr);
string out_file_name;
std::ostringstream stm;
stm << "out_contact"
<< ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write 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");
}
CHKERR m_field.getInterface<FieldBlas>()->fieldCopy(1., "SPATIAL_POSITION",
"EIGEN_POSITIONS");
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);
if (delete_prisms) {
Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
CHKERR moab.get_adjacencies(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(contact_prisms);
ents_to_delete.merge(quads);
ents_to_delete.merge(subtract(quads_edges, tris_edges));
CHKERR moab.delete_entities(ents_to_delete);
}
if (flg_file_out) {
PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n", 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 1:
internal_stress_ref = {5., 5., 5., 0., 0., 0., 0., 0., 0.};
actual_stress_ref = {0., 0., 1., 0., 0., 0., 0., 0., 0.};
break;
case 2:
internal_stress_ref = {5., 5., 5., 0., 0., 0., 0., 0., 0.};
actual_stress_ref = {0., 5. / 3., 5. / 3., 0., 0., 0., 0., 0., 0.};
break;
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 = {96, 192};
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());
const double eps = 1e-10;
if (test_num == 4) {
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) {
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;
}
HookeInternalStressElement::OpGetInternalStress
Definition: HookeInternalStressElement.hpp:54
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:699
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MetaNeumannForces::addNeumannBCElements
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.
Definition: SurfacePressure.cpp:1974
MoFEM::CoreInterface::loop_dofs
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.
main
int main(int argc, char *argv[])
Definition: simple_contact_thermal.cpp:39
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
HookeInternalStressElement::OpSaveStress
Definition: HookeInternalStressElement.hpp:117
SimpleContactProblem::ConvectSlaveContactElement
Element used to integrate on master surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:206
MoFEM::Version::strVersion
std::string strVersion()
Definition: UnknownInterface.hpp:24
MoFEM::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition: PrismInterface.hpp:23
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::createSNES
auto createSNES(MPI_Comm comm)
Definition: PetscSmartObj.hpp:255
Mortar.hpp
MoFEM::createSmartDM
DEPRECATED auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Definition: PetscSmartObj.hpp:149
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
SimpleContactProblem::LoadScale::lAmbda
static double lAmbda
Definition: SimpleContact.hpp:37
MoFEM::UnknownInterface::getFileVersion
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.
Definition: UnknownInterface.cpp:16
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
out_file_name
char out_file_name[255]
Definition: initial_diffusion.cpp:53
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
BlockData
NonlinearElasticElement::BlockData BlockData
Definition: simple_contact_thermal.cpp:32
VolSideFe
VolumeElementForcesAndSourcesCoreOnSide VolSideFe
Definition: photon_diffusion.cpp:31
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMoFEMLoopFiniteElementsUpAndLowRank
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: DMMoFEM.cpp:567
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MetaNeumannForces::setMomentumFluxOperators
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.
Definition: SurfacePressure.cpp:2069
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MoFEM::Version
Definition: UnknownInterface.hpp:12
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1094
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::smartCreateDMMatrix
DEPRECATED auto smartCreateDMMatrix(DM dm)
Definition: DMMoFEM.hpp:1092
MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition: PetscSmartObj.hpp:230
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
temp
void temp(int x, int y=10)
Definition: simple.cpp:4
MoFEM::DMMoFEMSNESSetJacobian
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: DMMoFEM.cpp:759
MoFEM::smartCreateDMVector
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1107
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
init_temp
double init_temp
Definition: thermo_elastic.cpp:125
MoFEM::DMMoFEMCreateMoFEM
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: DMMoFEM.cpp:114
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
HookeInternalStressElement::OpPostProcHookeElement
Definition: HookeInternalStressElement.hpp:142
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
help
static char help[]
Definition: simple_contact_thermal.cpp:30
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
SimpleContactProblem::LoadScale
Definition: SimpleContact.hpp:35
MoFEM::SnesRhs
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:27
Range
MetaSpringBC::setSpringOperators
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.
Definition: SpringElement.cpp:1178
MoFEM::CoreTmp< 0 >::Initialize
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
HookeInternalStressElement::OpGetAnalyticalInternalStress
Definition: HookeInternalStressElement.hpp:88
std
Definition: enable_if.hpp:5
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
SimpleContactProblem::ConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:181
MoFEM::BlockData
Definition: MeshsetsManager.cpp:755
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
addElasticElement
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
Definition: HookeElement.cpp:533
setOperators
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
HookeInternalStressElement::OpInternalStrain_dx
Definition: HookeInternalStressElement.hpp:74
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEM::VolumeElementForcesAndSourcesCoreOnSide
Base volume element used to integrate on skeleton.
Definition: VolumeElementForcesAndSourcesCoreOnSide.hpp:22
MoFEM::CoreInterface::set_field_order
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.
MoFEM::SmartPetscObj< DM >
MoFEM::SnesMat
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:139
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
HookeInternalStressElement.hpp
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::CoreInterface::add_field
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.
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
F
@ F
Definition: free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MoFEM::DMMoFEMSNESSetFunction
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: DMMoFEM.cpp:718
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211