static char help[] =
"\n";
int main(
int argc,
char *argv[]) {
"-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";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
enum contact_tests {
EIGHT_CUBE = 1,
FOUR_SEASONS = 2,
T_INTERFACE = 3,
PUNCH_TOP_AND_MID = 4,
PUNCH_TOP_ONLY = 5,
PLANE_AXI = 6,
ARC_THREE_SURF = 7,
SMILING_FACE = 8,
SMILING_FACE_CONVECT = 9,
WAVE_2D = 10,
WAVE_2D_ALM = 11,
LAST_TEST
};
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
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;
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);
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);
}
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 DMRegister_MoFEM(dm_name);
dm = createSmartDM(m_field.
get_comm(), dm_name);
CHKERR DMSetType(dm, dm_name);
CHKERR DMMoFEMCreateMoFEM(dm, &m_field,
"CONTACT_PROB", bit_levels.back());
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
CHKERR DMMoFEMAddElement(dm,
"CONTACT_ELEM");
CHKERR DMMoFEMAddElement(dm,
"ELASTIC");
CHKERR DMMoFEMAddElement(dm,
"PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm,
"SPRING");
CHKERR DMMoFEMAddElement(dm,
"CONTACT_POST_PROC");
auto D = smartCreateDMVector(dm);
auto F = smartVectorDuplicate(
D);
auto Aij = smartCreateDMMatrix(dm);
CHKERR DMoFEMMeshToLocalVector(dm,
D, INSERT_VALUES, SCATTER_FORWARD);
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_ctx = SnesMethod::CTX_SNESNONE;
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++) {
CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
&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 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) {
CHKERR DMMoFEMSNESSetFunction(
dm, "CONTACT_ELEM",
get_contact_rhs(contact_problem, make_convective_master_element),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetFunction(
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem,
make_convective_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
CHKERR DMMoFEMSNESSetFunction(
dm, "CONTACT_ELEM",
get_contact_rhs(contact_problem, make_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetFunction(
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem, make_contact_element,
alm_flag),
PETSC_NULL, PETSC_NULL);
}
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetFunction(dm,
"SPRING", fe_spring_rhs_ptr, PETSC_NULL,
PETSC_NULL);
dirichlet_bc_ptr.get());
boost::shared_ptr<FEMethod> fe_null;
fe_null);
if (convect_pts == PETSC_TRUE) {
CHKERR DMMoFEMSNESSetJacobian(
dm, "CONTACT_ELEM",
get_contact_lhs(contact_problem, make_convective_master_element),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetJacobian(
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem,
make_convective_slave_element),
PETSC_NULL, PETSC_NULL);
} else {
CHKERR DMMoFEMSNESSetJacobian(
dm, "CONTACT_ELEM",
get_contact_lhs(contact_problem, make_contact_element, alm_flag),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetJacobian(
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem, make_contact_element,
alm_flag),
PETSC_NULL, PETSC_NULL);
}
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);
}
SNESConvergedReason snes_reason;
{
CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
CHKERR SNESSetFunction(snes,
F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
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);
CHKERR DMoFEMMeshToGlobalVector(dm,
D, INSERT_VALUES, SCATTER_REVERSE);
PetscPrintf(PETSC_COMM_WORLD, "Loop post proc\n");
CHKERR DMoFEMLoopFiniteElements(dm,
"ELASTIC", &post_proc);
PetscPrintf(PETSC_COMM_WORLD, "Loop energy\n");
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %9.9f\n",
{
std::ostringstream stm;
stm << "out"
<< ".h5m";
PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
out_file_name.c_str());
"PARALLEL=WRITE_PART");
}
moab::Core mb_post;
moab::Interface &moab_proc = mb_post;
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);
CHKERR DMoFEMLoopFiniteElements(dm,
"CONTACT_ELEM",
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) {
case EIGHT_CUBE:
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;
case T_INTERFACE:
expected_energy = 3.0e-04;
expected_contact_area = 1.75;
expected_nb_gauss_pts = 336;
break;
case PUNCH_TOP_AND_MID:
expected_energy = 3.125e-04;
expected_contact_area = 0.25;
expected_nb_gauss_pts = 84;
break;
case PUNCH_TOP_ONLY:
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;
case SMILING_FACE:
expected_energy = 0.000733553;
expected_contact_area = 3.0;
expected_nb_gauss_pts = 144;
break;
case SMILING_FACE_CONVECT:
expected_energy = 0.000733621;
expected_contact_area = 3.0;
expected_nb_gauss_pts = 144;
break;
case WAVE_2D:
expected_energy = 0.008537863;
expected_contact_area = 0.125;
expected_nb_gauss_pts = 384;
break;
case WAVE_2D_ALM:
expected_energy = 0.008538894;
expected_contact_area = 0.125;
expected_nb_gauss_pts = 384;
break;
default:
"Unknown test number %d", test_num);
}
"Wrong energy %6.4e != %6.4e", expected_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(
CHKERR addHOOpsFace3D(
"MESH_NODE_POSITIONS", *post_proc_contact_ptr,
false,
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");
CHKERR DMoFEMLoopFiniteElements(dm,
"CONTACT_POST_PROC",
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;
}
const std::string default_options
Implementation of linear elastic material.
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr double lambda
surface tension
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
FTensor::Index< 'i', SPACE_DIM > i
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto createSNES(MPI_Comm comm)
static constexpr double delta
Set Dirichlet boundary conditions on spatial displacements.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
structure grouping operators and data used for calculation of nonlinear elastic element
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE & getLoopFeRhs()
get rhs volume element
MyVolumeFE & getLoopFeEnergy()
get energy fe
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.