Testing implementation of simple contact element (for contact between surfaces with matching meshes) by verifying its tangent matrix
static char help[] =
"\n";
int main(
int argc,
char *argv[]) {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
try {
PetscBool flg_file;
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;
PetscBool eigen_pos_flag = PETSC_FALSE;
PetscBool use_reference_coordinates = PETSC_TRUE;
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",
CHKERR PetscOptionsInt(
"-my_order",
"approximation order for spatial positions", "", 1,
"-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 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 PetscOptionsGetBool(PETSC_NULL,
"",
"-test_ale", &test_ale,
PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL,
"",
"-my_use_reference_coordinates",
&use_reference_coordinates, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
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);
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);
}
};
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 (eigen_pos_flag) {
}
PetscRandom rctx;
PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
auto set_coord = [&](
VectorAdaptor &&field_data,
double *x,
double *y,
double *z) {
double value;
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;
PetscRandomGetValueReal(rctx, &value);
field_data[0] = value *
scale;
};
"SPATIAL_POSITION");
"LAGMULT");
if (eigen_pos_flag) {
set_coord, "SPATIAL_POSITION");
}
if (test_ale == PETSC_TRUE) {
set_coord, "MESH_NODE_POSITIONS");
} else {
{
}
}
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, eigen_pos_flag, "EIGEN_POSITIONS",
use_reference_coordinates);
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, eigen_pos_flag, "EIGEN_POSITIONS",
use_reference_coordinates);
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, eigen_pos_flag, "EIGEN_POSITIONS",
use_reference_coordinates);
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, eigen_pos_flag, "EIGEN_POSITIONS",
use_reference_coordinates);
return fe_lhs_simple_contact;
};
auto get_contact_material_rhs = [&](auto contact_problem, auto make_element,
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", eigen_pos_flag,
"EIGEN_POSITIONS");
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;
};
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");
if (test_ale == PETSC_TRUE) {
if (!eigen_pos_flag)
contact_problem->addContactElementALE(
"ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
"LAGMULT", contact_prisms);
else
contact_problem->addContactElementALE(
"ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
"LAGMULT", contact_prisms, eigen_pos_flag, "EIGEN_POSITIONS");
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);
"SPATIAL_POSITION");
"SPATIAL_POSITION");
"MESH_NODE_POSITIONS");
"MESH_NODE_POSITIONS");
"SPATIAL_POSITION");
"MATERIAL", "MESH_NODE_POSITIONS");
"MATERIAL");
}
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, PETSC_FALSE);
CHKERR DMMoFEMAddElement(dm,
"CONTACT_ELEM");
if (test_ale == PETSC_TRUE) {
CHKERR DMMoFEMAddElement(dm,
"ALE_CONTACT_ELEM");
CHKERR DMMoFEMAddElement(dm,
"MATERIAL");
}
auto D = smartCreateDMVector(dm);
auto F = smartVectorDuplicate(
D);
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 VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(
A, MAT_SPD, PETSC_TRUE);
auto fdA = smartMatDuplicate(
A, MAT_COPY_VALUES);
if (convect_pts == PETSC_TRUE) {
CHKERR DMMoFEMSNESSetFunction(
dm, "CONTACT_ELEM",
get_contact_rhs(contact_problem, make_convective_master_element),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetJacobian(
dm, "CONTACT_ELEM",
get_master_contact_lhs(contact_problem,
make_convective_master_element),
NULL, NULL);
CHKERR DMMoFEMSNESSetFunction(
dm, "CONTACT_ELEM",
get_master_traction_rhs(contact_problem,
make_convective_slave_element),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetJacobian(
dm, "CONTACT_ELEM",
get_master_traction_lhs(contact_problem,
make_convective_slave_element),
NULL, 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);
CHKERR DMMoFEMSNESSetJacobian(dm,
"CONTACT_ELEM",
get_master_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);
}
if (test_ale == PETSC_TRUE) {
CHKERR moab.get_connectivity(all_tets, nodes,
false);
CHKERR DMMoFEMSNESSetFunction(
dm, "ALE_CONTACT_ELEM",
get_contact_material_rhs(contact_problem, make_contact_element,
nodes),
PETSC_NULL, PETSC_NULL);
CHKERR DMMoFEMSNESSetJacobian(
dm, "ALE_CONTACT_ELEM",
get_simple_contact_ale_lhs(contact_problem, make_contact_element),
NULL, NULL);
CHKERR DMMoFEMSNESSetJacobian(
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);
}
SNESConvergedReason snes_reason;
{
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 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;
"Difference between hand-calculated tangent matrix and finite "
"difference matrix is too big");
}
}
}
return 0;
}
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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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
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.
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_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
VectorShallowArrayAdaptor< double > VectorAdaptor
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)
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.
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.
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.