v0.13.0
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/>.
*/
#include <Mortar.hpp>
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;
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",
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 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();
CHKERRQ(ierr);
// 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);
CHKERR m_field.getInterface<BitRefManager>()
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(), MBTET,
ref_level_meshset);
CHKERR m_field.getInterface<BitRefManager>()
->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);
CHKERR m_field.getInterface<BitRefManager>()
->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);
if (eigen_pos_flag) {
CHKERR m_field.add_field("EIGEN_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
3, MB_TAG_SPARSE, MF_ZERO);
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_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 (eigen_pos_flag) {
CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(
set_coord, "SPATIAL_POSITION");
}
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, 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,
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", 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;
};
// add fields to the global matrix by adding the element
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");
Range all_tets;
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");
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
SmartPetscObj<DM> dm;
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;
}
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:676
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
auto smartCreateDMMatrix
Get smart matrix from DM.
Definition: DMMoFEM.hpp:944
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:717
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1032
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 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 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_col(const std::string &fe_name, const std::string &name_row)=0
set field col 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
double tol
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:126
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto createSNES
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:148
auto createSmartDM
Creates smart DM object.
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:39
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
double A
double scale
Definition: plastic.cpp:103
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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 MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Element used to integrate on master surfaces. It convects integration points on slaves,...
int main(int argc, char *argv[])