v0.9.2
test_jacobian_of_simple_contact_element.cpp

Testing implementation of simple contact element (for contact between surfaces with matching meshes) by verifying its tangent matrix

/** \file test_jacobian_of_simple_contact_element.cpp
* \example test_jacobian_of_simple_contact_element.cpp
*
* Testing implementation of simple contact element (for contact between
* surfaces with matching meshes) by verifying its tangent matrix
*
**/
/* 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/>.
*/
using namespace std;
using namespace MoFEM;
static char help[] = "\n";
int main(int argc, char *argv[]) {
// Initialize MoFEM
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
// Create mesh database
moab::Core mb_instance; // create database
moab::Interface &moab = mb_instance; // create interface to database
try {
PetscBool flg_file;
char mesh_file_name[255];
PetscInt order = 1;
PetscInt order_lambda = 1;
PetscReal r_value = 1.;
PetscReal cn_value = 1.;
PetscBool is_newton_cotes = PETSC_FALSE;
PetscBool test_jacobian = PETSC_FALSE;
PetscBool convect_pts = PETSC_FALSE;
PetscBool test_ale = PETSC_FALSE;
PetscBool alm_flag = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_jacobian", &test_jacobian,
PETSC_NULL);
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-my_order",
"approximation order for spatial positions", "", 1,
&order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-my_order_lambda",
"default approximation order of Lagrange multipliers", "", 1,
&order_lambda, 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 integration rules are used",
"", PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
CHKERR PetscOptionsBool("-my_convect", "set to convect integration pts", "",
PETSC_FALSE, &convect_pts, PETSC_NULL);
CHKERR PetscOptionsBool("-my_alm_flag", "set to convect integration pts",
"", PETSC_FALSE, &alm_flag, PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_ale", &test_ale,
PETSC_NULL);
ierr = PetscOptionsEnd();
// Check if mesh file was provided
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
// Read mesh to MOAB
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;
auto add_prism_interface = [&](Range &contact_prisms, Range &master_tris,
Range &slave_tris,
std::vector<BitRefLevel> &bit_levels) {
PrismInterface *interface;
CHKERR m_field.getInterface(interface);
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);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(), MBTET,
ref_level_meshset);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(), MBPRISM,
ref_level_meshset);
// get faces and tets to split
CHKERR 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 interface->splitSides(ref_level_meshset, bit_levels.back(),
cubit_meshset, true, true, 0);
// clean meshsets
CHKERR moab.delete_entities(&ref_level_meshset, 1);
CHKERR m_field.getInterface<BitRefManager>()->shiftRightBitRef(1);
bit_levels.pop_back();
}
}
EntityHandle meshset_prisms;
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);
}
};
Range contact_prisms, master_tris, slave_tris;
std::vector<BitRefLevel> bit_levels;
bit_levels.push_back(BitRefLevel().set(0));
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 3, bit_levels.back());
CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
bit_levels);
CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
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_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
MB_TAG_SPARSE, MF_ZERO);
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);
CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
3, MB_TAG_SPARSE, MF_ZERO);
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
// build field
CHKERR m_field.build_fields();
PetscRandom rctx;
PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
auto set_coord = [&](VectorAdaptor &&field_data, double *x, double *y,
double *z) {
double value;
double scale = 0.5;
PetscRandomGetValue(rctx, &value);
field_data[0] = (*x) + (value - 0.5) * scale;
PetscRandomGetValue(rctx, &value);
field_data[1] = (*y) + (value - 0.5) * scale;
PetscRandomGetValue(rctx, &value);
field_data[2] = (*z) + (value - 0.5) * scale;
};
auto set_pressure = [&](VectorAdaptor &&field_data, double *x, double *y,
double *z) {
double value;
double scale = 1.0;
PetscRandomGetValueReal(rctx, &value);
field_data[0] = value * scale;
};
CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(set_coord,
"SPATIAL_POSITION");
CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(set_pressure,
"LAGMULT");
if (test_ale == PETSC_TRUE) {
CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(
set_coord, "MESH_NODE_POSITIONS");
} else {
// MESH_NODE_POSITIONS
{
Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
}
}
PetscRandomDestroy(&rctx);
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);
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();
contact_problem->setContactOperatorsRhs(
fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
"LAGMULT", is_alm);
return fe_rhs_simple_contact;
};
auto get_master_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 get_master_traction_rhs = [&](auto contact_problem, auto make_element,
bool alm_flag = 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", alm_flag);
return fe_rhs_simple_contact;
};
auto get_master_traction_lhs = [&](auto contact_problem, auto make_element,
bool alm_flag = 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", alm_flag);
return fe_lhs_simple_contact;
};
auto get_contact_material_rhs = [&](auto contact_problem, auto make_element,
Range &ale_nodes) {
auto fe_rhs_simple_contact_ale_material = make_element();
auto common_data_simple_contact = make_contact_common_data();
common_data_simple_contact->forcesOnlyOnEntitiesRow.clear();
common_data_simple_contact->forcesOnlyOnEntitiesRow = ale_nodes;
contact_problem->setContactOperatorsRhsALEMaterial(
fe_rhs_simple_contact_ale_material, common_data_simple_contact,
"SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", "MATERIAL");
return fe_rhs_simple_contact_ale_material;
};
auto get_simple_contact_ale_lhs = [&](auto contact_problem,
auto make_element) {
auto fe_lhs_simple_contact_ale = make_element();
auto common_data_simple_contact = make_contact_common_data();
contact_problem->setContactOperatorsLhsALE(
fe_lhs_simple_contact_ale, common_data_simple_contact,
"SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT");
return fe_lhs_simple_contact_ale;
};
auto get_simple_contact_ale_material_lhs =
[&](auto contact_problem, auto make_element, Range &ale_nodes) {
auto fe_lhs_simple_contact_material_ale = make_element();
auto common_data_simple_contact = make_contact_common_data();
common_data_simple_contact->forcesOnlyOnEntitiesRow.clear();
common_data_simple_contact->forcesOnlyOnEntitiesRow = ale_nodes;
contact_problem->setContactOperatorsLhsALEMaterial(
fe_lhs_simple_contact_material_ale, common_data_simple_contact,
"SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", "MATERIAL");
return fe_lhs_simple_contact_material_ale;
};
// add fields to the global matrix by adding the element
contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
"LAGMULT", contact_prisms);
Range all_tets;
if (test_ale == PETSC_TRUE) {
contact_problem->addContactElementALE(
"ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
"LAGMULT", contact_prisms);
Range faces;
CHKERR moab.get_adjacencies(contact_prisms, 2, false, faces,
moab::Interface::UNION);
Range tris = faces.subset_by_type(MBTRI);
CHKERR moab.get_adjacencies(tris, 3, false, all_tets,
moab::Interface::UNION);
// Add finite elements
CHKERR m_field.add_finite_element("MATERIAL", MF_ZERO);
"SPATIAL_POSITION");
"SPATIAL_POSITION");
"MESH_NODE_POSITIONS");
"MESH_NODE_POSITIONS");
"SPATIAL_POSITION");
"MATERIAL", "MESH_NODE_POSITIONS");
CHKERR m_field.add_ents_to_finite_element_by_type(all_tets, MBTET,
"MATERIAL");
CHKERR m_field.build_finite_elements("MATERIAL", &all_tets);
}
// build finite elemnts
// build adjacencies
CHKERR m_field.build_adjacencies(bit_levels.back());
// define problems
CHKERR m_field.add_problem("CONTACT_PROB");
// set refinement level for problem
CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
bit_levels.back());
DMType dm_name = "DMMOFEM";
// create dm instance
dm = createSmartDM(m_field.get_comm(), dm_name);
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, PETSC_FALSE);
// add elements to dm
CHKERR DMMoFEMAddElement(dm, "CONTACT_ELEM");
if (test_ale == PETSC_TRUE) {
CHKERR DMMoFEMAddElement(dm, "ALE_CONTACT_ELEM");
CHKERR DMMoFEMAddElement(dm, "MATERIAL");
}
CHKERR DMSetUp(dm);
// Vector of DOFs and the RHS
auto D = smartCreateDMVector(dm);
// Stiffness matrix
auto A = 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 VecZeroEntries(F);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(A, MAT_SPD, PETSC_TRUE);
CHKERR MatZeroEntries(A);
auto fdA = smartMatDuplicate(A, MAT_COPY_VALUES);
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_contact_lhs(contact_problem,
make_convective_master_element),
NULL, NULL);
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem,
make_convective_slave_element),
PETSC_NULL, PETSC_NULL);
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem,
make_convective_slave_element),
NULL, 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 DMMoFEMSNESSetJacobian(dm, "CONTACT_ELEM",
get_master_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);
}
if (test_ale == PETSC_TRUE) {
Range nodes;
CHKERR moab.get_connectivity(all_tets, nodes, false);
dm, "ALE_CONTACT_ELEM",
get_contact_material_rhs(contact_problem, make_contact_element,
nodes),
PETSC_NULL, PETSC_NULL);
dm, "ALE_CONTACT_ELEM",
get_simple_contact_ale_lhs(contact_problem, make_contact_element),
NULL, NULL);
dm, "ALE_CONTACT_ELEM",
get_simple_contact_ale_material_lhs(contact_problem,
make_contact_element, nodes),
NULL, NULL);
}
if (test_jacobian == PETSC_TRUE) {
char testing_options[] =
"-snes_test_jacobian -snes_test_jacobian_display "
"-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 "
"-snes_max_it "
"1 ";
CHKERR PetscOptionsInsertString(NULL, testing_options);
} else {
char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
"-snes_rtol 0 "
"-snes_max_it 1 ";
CHKERR PetscOptionsInsertString(NULL, testing_options);
}
auto snes = MoFEM::createSNES(m_field.get_comm());
SNESConvergedReason snes_reason;
SnesCtx *snes_ctx;
// create snes nonlinear solver
{
CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
}
CHKERR SNESSolve(snes, PETSC_NULL, D);
if (test_jacobian == PETSC_FALSE) {
double nrm_A0;
CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
char testing_options_fd[] = "-snes_fd";
CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
CHKERR SNESSolve(snes, NULL, D);
CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
double nrm_A;
CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
nrm_A / nrm_A0);
nrm_A /= nrm_A0;
constexpr double tol = 1e-6;
if (nrm_A > tol) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Difference between hand-calculated tangent matrix and finite "
"difference matrix is too big");
}
}
}
// finish work cleaning memory, getting statistics, etc
return 0;
}