static char help[] =
"\n";
int main(
int argc,
char *argv[]) {
const string default_options = "-ksp_type fgmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_13 1 \n"
"-mat_mumps_icntl_14 800 \n"
"-mat_mumps_icntl_20 0 \n"
"-mat_mumps_icntl_24 1 \n"
"-snes_type newtonls \n"
"-snes_linesearch_type basic \n"
"-snes_divergence_tolerance 0 \n"
"-snes_max_it 50 \n"
"-snes_atol 1e-8 \n"
"-snes_rtol 1e-10 \n"
"-snes_monitor \n"
"-ksp_monitor \n"
"-snes_converged_reason \n"
"-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";
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();
}
}
FOUR_SEASONS = 2,
PLANE_AXI = 6,
ARC_THREE_SURF = 7,
};
try {
PetscBool flg_file;
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 wave_surf_flag = PETSC_FALSE;
PetscInt wave_dim = 2;
PetscInt wave_surf_block_id = 1;
PetscReal wave_length = 1.0;
PetscReal wave_ampl = 0.01;
PetscReal mesh_height = 1.0;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
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_wave_surf",
"if set true, make one of the surfaces wavy", "",
PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
CHKERR PetscOptionsInt(
"-my_wave_surf_block_id",
"make wavy surface of the block with this id", "",
wave_surf_block_id, &wave_surf_block_id, PETSC_NULL);
CHKERR PetscOptionsInt(
"-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 PetscOptionsReal(
"-my_mesh_height",
"vertical dimension of the mesh ", "", mesh_height,
&mesh_height, 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 = "";
auto add_prism_interface = [&](
Range &contact_prisms,
Range &master_tris,
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);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
ref_level_meshset);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
ref_level_meshset);
CHKERR interface->
getSides(cubit_meshset, bit_levels.back(),
true, 0);
bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
cubit_meshset, true, true, 0);
CHKERR moab.delete_entities(&ref_level_meshset, 1);
->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
true);
->updateMeshsetByEntitiesChildren(cubit_meshset,
bit_levels.back(),
cubit_meshset, MBEDGE, true);
->updateMeshsetByEntitiesChildren(cubit_meshset,
bit_levels.back(),
cubit_meshset, MBTRI, true);
->updateMeshsetByEntitiesChildren(cubit_meshset,
bit_levels.back(),
cubit_meshset, MBTET, true);
}
bit_levels.pop_back();
}
}
CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(),
MBPRISM, meshset_prisms);
CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
CHKERR moab.delete_entities(&meshset_prisms, 1);
for (Range::iterator pit = contact_prisms.begin();
pit != contact_prisms.end(); pit++) {
CHKERR moab.side_element(*pit, 2, 3, tri);
master_tris.insert(tri);
CHKERR moab.side_element(*pit, 2, 4, tri);
slave_tris.insert(tri);
}
};
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);
}
};
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);
};
Range contact_prisms, master_tris, slave_tris;
std::vector<BitRefLevel> bit_levels;
0, 3, bit_levels.back());
CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
bit_levels);
if (wave_surf_flag) {
CHKERR make_wavy_surface(wave_surf_block_id, wave_dim, wave_length,
wave_ampl, mesh_height);
}
if (order_contact >
order) {
CHKERR set_contact_order(contact_prisms, order_contact, nb_ho_levels);
}
{
}
{
}
boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(
new Hooke<adouble>());
boost::shared_ptr<Hooke<double>> hooke_double_ptr(
new Hooke<double>());
false, false);
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);
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);
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);
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);
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);
contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
"LAGMULT", contact_prisms);
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);
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);
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);
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);
}
SNESConvergedReason snes_reason;
{
CHKERR SNESSetFromOptions(snes);
}
std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
}
for (int ss = 0; ss != nb_sub_steps; ++ss) {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load scale: %6.4e\n",
CHKERR SNESSolve(snes, PETSC_NULL,
D);
CHKERR SNESGetConvergedReason(snes, &snes_reason);
int its;
CHKERR SNESGetIterationNumber(snes, &its);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Number of Newton iterations = %D\n",
its);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
PetscPrintf(PETSC_COMM_WORLD, "Loop post proc\n");
elastic.getLoopFeEnergy().eNergy = 0;
PetscPrintf(PETSC_COMM_WORLD, "Loop energy\n");
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %9.9f\n",
elastic.getLoopFeEnergy().eNergy);
{
std::ostringstream stm;
stm << "out"
<< ".h5m";
PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
out_file_name.c_str());
"PARALLEL=WRITE_PART");
}
auto common_data_simple_contact = make_contact_common_data();
boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
fe_post_proc_simple_contact;
if (convect_pts == PETSC_TRUE) {
fe_post_proc_simple_contact = make_convective_master_element();
} else {
fe_post_proc_simple_contact = make_contact_element();
}
contact_problem->setContactOperatorsForPostProc(
fe_post_proc_simple_contact, common_data_simple_contact, m_field,
"SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag);
mb_post.delete_mesh();
CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
fe_post_proc_simple_contact);
std::array<double, 2> nb_gauss_pts;
std::array<double, 2> contact_area;
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: %9.9f out of %9.9f\n",
contact_area[0], contact_area[1]);
}
if (test_num) {
double expected_energy, expected_contact_area;
int expected_nb_gauss_pts;
constexpr
double eps = 1e-6;
switch (test_num) {
expected_energy = 3.0e-04;
expected_contact_area = 3.0;
expected_nb_gauss_pts = 576;
break;
case FOUR_SEASONS:
expected_energy = 1.2e-01;
expected_contact_area = 106.799036701;
expected_nb_gauss_pts = 672;
break;
expected_energy = 3.0e-04;
expected_contact_area = 1.75;
expected_nb_gauss_pts = 336;
break;
expected_energy = 3.125e-04;
expected_contact_area = 0.25;
expected_nb_gauss_pts = 84;
break;
expected_energy = 0.000096432;
expected_contact_area = 0.25;
expected_nb_gauss_pts = 336;
break;
case PLANE_AXI:
expected_energy = 0.000438889;
expected_contact_area = 0.784409608;
expected_nb_gauss_pts = 300;
break;
case ARC_THREE_SURF:
expected_energy = 0.002573411;
expected_contact_area = 2.831455633;
expected_nb_gauss_pts = 228;
break;
expected_energy = 0.000733553;
expected_contact_area = 3.0;
expected_nb_gauss_pts = 144;
break;
expected_energy = 0.000733621;
expected_contact_area = 3.0;
expected_nb_gauss_pts = 144;
break;
expected_energy = 0.008537863;
expected_contact_area = 0.125;
expected_nb_gauss_pts = 384;
break;
expected_energy = 0.008538894;
expected_contact_area = 0.125;
expected_nb_gauss_pts = 384;
break;
default:
"Unknown test number %d", test_num);
}
if (std::abs(elastic.getLoopFeEnergy().eNergy - expected_energy) >
eps) {
"Wrong energy %6.4e != %6.4e", expected_energy,
elastic.getLoopFeEnergy().eNergy);
}
if ((int)nb_gauss_pts[0] != expected_nb_gauss_pts) {
"Wrong number of active gauss pts %d != %d",
expected_nb_gauss_pts, (int)nb_gauss_pts[0]);
}
if (std::abs(contact_area[0] - expected_contact_area) >
eps) {
"Wrong active contact area %6.4e != %6.4e",
expected_contact_area, contact_area[0]);
}
}
}
{
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(
false);
CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"LAGMULT");
CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"SPATIAL_POSITION");
CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
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(
}
}
return 0;
}