Implementation of simple contact between surfaces with matching meshes taking into account internal stress resulting from the thermal expansion
static char help[] =
"\n";
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();
}
}
try {
PetscBool flg_file;
PetscBool flg_file_out;
char output_mesh_name[255];
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 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",
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,
"-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);
"-my_stress_tag_name", "stress tag name file name", "",
"INTERNAL_STRESS", stress_tag_name, 255, &flg_tag_name);
"-my_thermal_expansion_coef", "thermal expansion coef ", "",
thermal_expansion_coef, &thermal_expansion_coef, 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();
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
if (is_partitioned == PETSC_TRUE) {
"Partitioned mesh is not supported");
}
const char *option;
option = "";
std::vector<BitRefLevel> bit_levels;
CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
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());
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(
ref_level_meshset);
CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
true, 0);
bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
CHKERR prism_interface->splitSides(ref_level_meshset,
bit_levels.back(), cubit_meshset,
true, true, 0);
CHKERR moab.delete_entities(&ref_level_meshset, 1);
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,
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();
if (id == block_id) {
bit->getMeshset(), 3, tets,
true);
all_tets.merge(tets);
}
}
}
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();
if (id == block_id) {
bit->getMeshset(), 3, tets,
true);
all_tets.merge(tets);
}
}
}
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] -=
(1. - cos(2. * M_PI * x /
lambda) * cos(2. * M_PI * y /
lambda));
break;
default:
"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);
}
if (!eigen_pos_flag) {
}
if (!eigen_pos_flag) {
}
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);
ho_ents.merge(contact_tris);
ho_ents.merge(contact_edges);
for (int ll = 0; ll < nb_ho_levels; ll++) {
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);
CHKERR moab.get_adjacencies(tets,
d,
false, ents,
moab::Interface::UNION);
}
ho_ents = unite(ho_ents, ents);
ho_ents = unite(ho_ents, tets);
}
order_contact);
};
if (!ignore_contact && order_contact >
order) {
CHKERR set_contact_order(contact_prisms, order_contact, nb_ho_levels);
}
{
}
{
}
CHKERR moab.get_connectivity(slave_tris, slave_verts,
true);
slave_verts, "LAGMULT");
boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
boost::make_shared<std::map<int, BlockData>>();
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();
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS", false);
auto data_hooke_element_at_pts =
boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
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](
t_thermal_strain(
i,
j) = 0.0;
switch (test_num) {
case 0:
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));
if (
bit->getName().compare(0, 11,
"MAT_ELASTIC") == 0) {
CHKERR moab.get_entities_by_dimension(
bit->getMeshset(), 3, tets,
true);
all_tets.merge(tets);
}
}
Skinner skinner(&moab);
CHKERR skinner.find_skin(0, all_tets,
false, skin_tris);
"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);
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);
"MESH_NODE_POSITIONS");
bit_levels.back());
DMType dm_name = "DMMOFEM";
CHKERR DMSetType(dm, dm_name);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
fe_elastic_rhs_ptr->snes_f =
F;
fe_elastic_lhs_ptr->snes_B = Aij;
boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
m_field,
"SPATIAL_POSITION", Aij,
D,
F));
dirichlet_bc_ptr->snes_x =
D;
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->getLoopFe(), NULL, NULL);
}
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 VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
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);
}
PETSC_NULL);
PETSC_NULL);
dirichlet_bc_ptr.get());
boost::shared_ptr<FEMethod> fe_null;
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);
}
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);
}
{
CHKERR SNESSetFromOptions(snes);
}
CHKERR post_proc_skin.generateReferenceElementMesh();
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) {}
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) {
} 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);
}
"MESH_NODE_POSITIONS", false, false,
v_energy);
const double *eng_ptr;
CHKERR VecGetArrayRead(v_energy, &eng_ptr);
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy: %8.8e\n", *eng_ptr);
{
PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
ostringstream stm;
stm << "out_skin.h5m";
PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
out_file_name.c_str());
CHKERR post_proc_skin.writeFile(stm.str());
}
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);
fe_post_proc_simple_contact);
auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
const double *array;
CHKERR VecGetArrayRead(vec, &array);
}
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);
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.);
}
std::ostringstream strm;
strm << "out_contact_integ_pts"
<< ".h5m";
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
"PARALLEL=WRITE_PART");
}
boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
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) {
post_proc_contact_ptr);
std::ostringstream stm;
stm << "out_contact"
<< ".h5m";
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
CHKERR post_proc_contact_ptr->postProcMesh.write_file(
}
"EIGEN_POSITIONS");
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) {
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());
};
if (test_num) {
CHKERR moab.get_entities_by_dimension(0, 3, tets);
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:
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) {
"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) {
"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) {
"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) {
"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) {
"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) {
"Wrong component %d of actual stress: should be %g "
"but is %g",
i, actual_stress_ref[
i], actual_stress[
i]);
}
}
}
}
}
}
}
return 0;
}