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"
"-order 1 \n"
"-order_lambda 1 \n"
"-order_contact 2 \n"
"-ho_levels_num 1 \n"
"-step_num 1 \n"
"-cn_value 1. \n"
"-alm_flag 0 \n";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
try {
PetscBool flg_file;
PetscInt nb_sub_steps = 1;
PetscBool is_partitioned = PETSC_FALSE;
PetscInt test_num = 0;
PetscBool wave_surf_flag = PETSC_FALSE;
PetscInt wave_dim = 2;
PetscBool is_displacement_field = PETSC_FALSE;
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(
"-file_name",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsInt(
"-order",
"approximation order of spatial positions",
"", 1, &
order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-step_num",
"number of steps",
"", nb_sub_steps,
&nb_sub_steps, PETSC_NULL);
CHKERR PetscOptionsBool(
"-is_partitioned",
"set if mesh is partitioned (this result that each "
"process keeps only part of the mesh",
"", PETSC_FALSE, &is_partitioned, PETSC_NULL);
CHKERR PetscOptionsInt(
"-test_num",
"test number",
"", 0, &test_num,
PETSC_NULL);
CHKERR PetscOptionsBool(
"-is_displacement_field",
"use displacement field",
"", PETSC_FALSE, &is_displacement_field,
PETSC_NULL);
CHKERR PetscOptionsBool(
"-wave_surf",
"if set true, make one of the surfaces wavy", "",
PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
CHKERR PetscOptionsInt(
"-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(
"-wave_dim",
"dimension (2 or 3)",
"", wave_dim,
&wave_dim, PETSC_NULL);
CHKERR PetscOptionsReal(
"-wave_length",
"profile wavelength",
"",
wave_length, &wave_length, PETSC_NULL);
CHKERR PetscOptionsReal(
"-wave_ampl",
"profile amplitude",
"", wave_ampl,
&wave_ampl, PETSC_NULL);
CHKERR PetscOptionsReal(
"-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 -file_name (MESH FILE NEEDED)");
}
if (is_partitioned == PETSC_TRUE)
"Partitioned mesh is not supported");
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
MeshsetsManager *mmanager_ptr;
CHKERR mmanager_ptr->printDisplacementSet();
CHKERR mmanager_ptr->printForceSet();
CHKERR mmanager_ptr->printMaterialsSet();
string primary_field_name;
if (!is_displacement_field)
primary_field_name = "SPATIAL_POSITION";
else
primary_field_name = "U";
"MESH_NODE_POSITIONS",
is_displacement_field);
"MESH_NODE_POSITIONS",
is_displacement_field);
m_field, primary_field_name, "MESH_NODE_POSITIONS", "CONTACT_PROB",
"ELASTIC", is_displacement_field, true,
if (!is_displacement_field) {
Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
}
{
Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
}
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
SmartPetscObj<DM> dm;
dm = createSmartDM(m_field.
get_comm(), dm_name);
CHKERR DMSetType(dm, dm_name);
CHKERR DMMoFEMCreateMoFEM(dm, &m_field,
"CONTACT_PROB", bit_ref);
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
auto D = smartCreateDMVector(dm);
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 20 "
"-snes_atol 1e-8 "
"-snes_rtol 1e-8 ";
CHKERR PetscOptionsInsertString(NULL, testing_options);
}
CHKERR SNESSetFromOptions(snes);
for (int ss = 0; ss != nb_sub_steps; ++ss) {
CHKERR SNESSolve(snes, PETSC_NULL,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToGlobalVector(dm,
D, INSERT_VALUES, SCATTER_REVERSE);
std::ofstream ofs((std::string("test_simple_contact") + ".txt").c_str());
m_field, primary_field_name, common_data_post, ofs));
ofs << "Elastic energy: " << energy << endl;
ofs.flush();
ofs.close();
if (test_num) {
int expected_nb_gauss_pts;
double expected_energy, expected_contact_area;
expected_contact_area,
expected_nb_gauss_pts);
constexpr double eps = 1e-6;
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]);
}
}
if (std::abs(energy - expected_energy) >
eps)
"Wrong energy %6.4e != %6.4e", expected_energy, energy);
}
}
}
return 0;
}
const std::string default_options
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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
auto createSNES(MPI_Comm comm)
MoFEMErrorCode check_tests(int ss, int test_num, double &expected_energy, double &expected_contact_area, int &expected_nb_gauss_pts)
Set of functions declaring elements and setting operators for basic boundary conditions interface.
MoFEMErrorCode setupSolverFunctionSNES() override
MoFEMErrorCode getCommandLineParameters() override
MoFEMErrorCode addElementFields() override
MoFEMErrorCode setOperators() override
MoFEMErrorCode createElements() override
MoFEMErrorCode addElementsToDM(SmartPetscObj< DM > dm) override
MoFEMErrorCode postProcessElement(int step) override
MoFEMErrorCode setupSolverJacobianSNES() override
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Set of functions declaring elements and setting operators for generic element interface.
MoFEMErrorCode createElements()
MoFEMErrorCode setupSolverFunctionSNES()
MoFEMErrorCode setOperators()
MoFEMErrorCode setupSolverJacobianSNES()
MoFEMErrorCode addElementFields()
boost::shared_ptr< NonlinearElasticElement > elasticElementPtr
MoFEMErrorCode addElementsToDM(SmartPetscObj< DM > dm)
MoFEMErrorCode getCommandLineParameters()
MoFEMErrorCode postProcessElement(int step)