v0.14.0
Functions | Variables
mortar_contact.cpp File Reference
#include <Mortar.hpp>
#include <MortarCtestFunctions.hpp>
#include <ElasticMaterials.hpp>
#include <NonlinearElasticElementInterface.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
mortar_contact.cpp.

Definition at line 33 of file mortar_contact.cpp.

33  {
34 
35  const string default_options = "-ksp_type fgmres \n"
36  "-pc_type lu \n"
37  "-pc_factor_mat_solver_type mumps \n"
38  "-mat_mumps_icntl_13 1 \n"
39  "-mat_mumps_icntl_14 800 \n"
40  "-mat_mumps_icntl_20 0 \n"
41  "-mat_mumps_icntl_24 1 \n"
42  "-snes_type newtonls \n"
43  "-snes_linesearch_type basic \n"
44  "-snes_divergence_tolerance 0 \n"
45  "-snes_max_it 50 \n"
46  "-snes_atol 1e-8 \n"
47  "-snes_rtol 1e-10 \n"
48  "-snes_monitor \n"
49  "-ksp_monitor \n"
50  "-snes_converged_reason \n"
51  "-order 1 \n"
52  "-order_lambda 1 \n"
53  "-order_contact 2 \n"
54  "-ho_levels_num 1 \n"
55  "-step_num 1 \n"
56  "-cn_value 1. \n"
57  "-alm_flag 0 \n";
58 
59  string param_file = "param_file.petsc";
60  if (!static_cast<bool>(ifstream(param_file))) {
61  std::ofstream file(param_file.c_str(), std::ios::ate);
62  if (file.is_open()) {
63  file << default_options;
64  file.close();
65  }
66  }
67 
68  // Initialize MoFEM
69  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
70 
71  // Create mesh database
72  moab::Core mb_instance; // create database
73  moab::Interface &moab = mb_instance; // create interface to database
74 
75  try {
76  PetscBool flg_file;
77 
78  char mesh_file_name[255];
79  PetscInt order = 1;
80  PetscInt nb_sub_steps = 1;
81  PetscBool is_partitioned = PETSC_FALSE;
82  PetscInt test_num = 0;
83  PetscBool wave_surf_flag = PETSC_FALSE;
84  PetscInt wave_dim = 2;
85  PetscBool is_displacement_field = PETSC_FALSE;
86  PetscInt wave_surf_block_id = 1;
87  PetscReal wave_length = 1.0;
88  PetscReal wave_ampl = 0.01;
89  PetscReal mesh_height = 1.0;
90 
91  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
92 
93  CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
94  mesh_file_name, 255, &flg_file);
95 
96  CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
97  "", 1, &order, PETSC_NULL);
98 
99  CHKERR PetscOptionsInt("-step_num", "number of steps", "", nb_sub_steps,
100  &nb_sub_steps, PETSC_NULL);
101 
102  CHKERR PetscOptionsBool("-is_partitioned",
103  "set if mesh is partitioned (this result that each "
104  "process keeps only part of the mesh",
105  "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
106 
107  CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
108  PETSC_NULL);
109 
110  CHKERR PetscOptionsBool("-is_displacement_field", "use displacement field",
111  "", PETSC_FALSE, &is_displacement_field,
112  PETSC_NULL);
113  CHKERR PetscOptionsBool("-wave_surf",
114  "if set true, make one of the surfaces wavy", "",
115  PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
116  CHKERR PetscOptionsInt("-wave_surf_block_id",
117  "make wavy surface of the block with this id", "",
118  wave_surf_block_id, &wave_surf_block_id, PETSC_NULL);
119  CHKERR PetscOptionsInt("-wave_dim", "dimension (2 or 3)", "", wave_dim,
120  &wave_dim, PETSC_NULL);
121  CHKERR PetscOptionsReal("-wave_length", "profile wavelength", "",
122  wave_length, &wave_length, PETSC_NULL);
123  CHKERR PetscOptionsReal("-wave_ampl", "profile amplitude", "", wave_ampl,
124  &wave_ampl, PETSC_NULL);
125  CHKERR PetscOptionsReal("-mesh_height", "vertical dimension of the mesh ",
126  "", mesh_height, &mesh_height, PETSC_NULL);
127 
128  ierr = PetscOptionsEnd();
129  CHKERRQ(ierr);
130 
131  // Check if mesh file was provided
132  if (flg_file != PETSC_TRUE) {
133  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
134  }
135  // FIXME: really???
136  if (is_partitioned == PETSC_TRUE)
137  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
138  "Partitioned mesh is not supported");
139 
140  // CHKERR moab.load_file(mesh_file_name, 0, option);
141 
142  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
143  if (pcomm == NULL)
144  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
145 
146  // Create MoFEM database and link it to MoAB
147  MoFEM::Core core(moab);
148  MoFEM::Interface &m_field = core;
149 
150  auto simple = m_field.getInterface<Simple>();
151  simple->loadFile("", mesh_file_name);
152 
153  MeshsetsManager *mmanager_ptr;
154  CHKERR m_field.getInterface(mmanager_ptr);
155  CHKERR mmanager_ptr->printDisplacementSet();
156  CHKERR mmanager_ptr->printForceSet();
157  // print block sets with materials
158  CHKERR mmanager_ptr->printMaterialsSet();
159 
160  string primary_field_name;
161 
162  if (!is_displacement_field)
163  primary_field_name = "SPATIAL_POSITION";
164  else
165  primary_field_name = "U";
166  // create interfaces (could be based on cmake file ?)
167 
168  MortarContactInterface m_contact(m_field, primary_field_name,
169  "MESH_NODE_POSITIONS",
170  is_displacement_field);
171  NonlinearElasticElementInterface m_elastic(m_field, primary_field_name,
172  "MESH_NODE_POSITIONS",
173  is_displacement_field);
174 
175  BasicBoundaryConditionsInterface m_boundary(
176  m_field, primary_field_name, "MESH_NODE_POSITIONS", "CONTACT_PROB",
177  "ELASTIC", is_displacement_field, true,
178  &MortarContactProblem::LoadScale::lAmbda, is_partitioned);
179 
180  CHKERR m_boundary.getCommandLineParameters();
181  CHKERR m_contact.getCommandLineParameters();
182  CHKERR m_elastic.getCommandLineParameters();
183 
184  CHKERR m_boundary.addElementFields();
185  CHKERR m_contact.addElementFields();
186  CHKERR m_elastic.addElementFields();
187 
188  // build field
189  CHKERR m_field.build_fields();
190 
191  CHKERR m_boundary.createElements();
192  CHKERR m_contact.createElements();
193  CHKERR m_elastic.createElements();
194 
195  // Projection on "x" field
196  if (!is_displacement_field) {
197  Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
198  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
199  }
200  // Projection on "X" field
201  {
202  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
203  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
204  }
205 
206  // build finite elemnts
207  CHKERR m_field.build_finite_elements();
208 
209  auto bit_ref = m_contact.getBitRefLevel();
210  // build adjacencies
211  CHKERR m_field.build_adjacencies(bit_ref);
212 
213  // define problems
214  CHKERR m_field.add_problem("CONTACT_PROB");
215 
216  // set refinement level for problem
217  CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB", bit_ref);
218 
219  DMType dm_name = "DMMOFEM";
220  CHKERR DMRegister_MoFEM(dm_name);
221 
222  SmartPetscObj<DM> dm;
223  dm = createSmartDM(m_field.get_comm(), dm_name);
224 
225  // create dm instance
226  CHKERR DMSetType(dm, dm_name);
227 
228  // set dm datastruture which created mofem datastructures
229  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_ref);
230  CHKERR DMSetFromOptions(dm);
231  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
232  // add elements to dm
233 
234  CHKERR m_boundary.setOperators();
235  CHKERR m_contact.setOperators();
236  CHKERR m_elastic.setOperators();
237 
238  CHKERR m_boundary.addElementsToDM(dm);
239  CHKERR m_contact.addElementsToDM(dm);
240  CHKERR m_elastic.addElementsToDM(dm);
241 
242  CHKERR DMSetUp(dm);
243 
244  CHKERR m_boundary.setupSolverFunctionSNES();
245  CHKERR m_contact.setupSolverFunctionSNES();
246  CHKERR m_elastic.setupSolverFunctionSNES();
247 
248  CHKERR m_boundary.setupSolverJacobianSNES();
249  CHKERR m_contact.setupSolverJacobianSNES();
250  CHKERR m_elastic.setupSolverJacobianSNES();
251 
252 // #ifndef NDEBUG
253 // cout << "<<< checkMPIAIJWithArraysMatrixFillIn <<< " << endl;
254 // CHKERR m_field.getInterface<MatrixManager>()
255 // ->checkMPIAIJWithArraysMatrixFillIn<PetscGlobalIdx_mi_tag>(
256 // "CONTACT_PROB", -2, -2, 0);
257 // #endif // NDEBUG
258 
259  auto D = smartCreateDMVector(dm);
260  if (test_num) {
261  char testing_options[] = "-ksp_type fgmres "
262  "-pc_type lu "
263  "-pc_factor_mat_solver_type mumps "
264  "-snes_type newtonls "
265  "-snes_linesearch_type basic "
266  "-snes_max_it 20 "
267  "-snes_atol 1e-8 "
268  "-snes_rtol 1e-8 ";
269  CHKERR PetscOptionsInsertString(NULL, testing_options);
270  }
271 
272  auto snes = MoFEM::createSNES(m_field.get_comm());
273  CHKERR SNESSetDM(snes, dm);
274  CHKERR SNESSetFromOptions(snes);
275 
276  for (int ss = 0; ss != nb_sub_steps; ++ss) {
277 
278  MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
279  MOFEM_LOG_C("WORLD", Sev::inform, "Load lambda %6.4e",
280  MortarContactProblem::LoadScale::lAmbda);
281 
282  CHKERR SNESSolve(snes, PETSC_NULL, D);
283 
284  // save on mesh
285  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
286  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
287  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
288 
289  CHKERR m_elastic.postProcessElement(ss);
290 
291  auto common_data_post = m_contact.commonDataPostProc;
292  std::ofstream ofs((std::string("test_simple_contact") + ".txt").c_str());
293 
294  m_contact.getPostProcSimpleContactPipeline().push_back(
295  new MortarContactProblem::OpMakeTestTextFile(
296  m_field, primary_field_name, common_data_post, ofs));
297 
298  CHKERR m_contact.postProcessElement(ss);
299  CHKERR m_boundary.postProcessElement(ss);
300 
301  double energy = m_elastic.elasticElementPtr->getLoopFeEnergy().eNergy;
302 
303  ofs << "Elastic energy: " << energy << endl;
304  ofs.flush();
305  ofs.close();
306 
307  if (test_num) {
308 
309  auto &nb_gauss_pts = m_contact.nbGaussPts;
310  auto &contact_area = m_contact.contactArea;
311 
312  int expected_nb_gauss_pts;
313  double expected_energy, expected_contact_area;
314 
315  MortarCtestFunctions::check_tests(ss, test_num, expected_energy,
316  expected_contact_area,
317  expected_nb_gauss_pts);
318 
319  constexpr double eps = 1e-6;
320  if (m_field.get_comm_rank() == 0) {
321  if ((int)nb_gauss_pts[0] != expected_nb_gauss_pts) {
322  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
323  "Wrong number of active gauss pts %d != %d",
324  expected_nb_gauss_pts, (int)nb_gauss_pts[0]);
325  }
326  if (std::abs(contact_area[0] - expected_contact_area) > eps) {
327  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
328  "Wrong active contact area %6.4e != %6.4e",
329  expected_contact_area, contact_area[0]);
330  }
331  }
332 
333  if (std::abs(energy - expected_energy) > eps)
334  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
335  "Wrong energy %6.4e != %6.4e", expected_energy, energy);
336  }
337  }
338  }
339  CATCH_ERRORS;
340 
341  // finish work cleaning memory, getting statistics, etc
342  MoFEM::Core::Finalize();
343 
344  return 0;
345 }

Variable Documentation

◆ help

char help[] = "\n"
static
Examples
mortar_contact.cpp.

Definition at line 30 of file mortar_contact.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MortarCtestFunctions::check_tests
MoFEMErrorCode check_tests(int ss, int test_num, double &expected_energy, double &expected_contact_area, int &expected_nb_gauss_pts)
Definition: MortarCtestFunctions.hpp:229
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
BasicBoundaryConditionsInterface
Set of functions declaring elements and setting operators for basic boundary conditions interface.
Definition: BasicBoundaryConditionsInterface.hpp:20
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::createSNES
auto createSNES(MPI_Comm comm)
Definition: PetscSmartObj.hpp:255
MoFEM::createSmartDM
DEPRECATED auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Definition: PetscSmartObj.hpp:149
SimpleContactProblem::OpMakeTestTextFile
Definition: SimpleContact.hpp:2207
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MortarContactInterface
Definition: MortarContactInterface.hpp:22
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:304
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
help
static char help[]
Definition: mortar_contact.cpp:30
MoFEM::smartCreateDMVector
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1107
MoFEM::DMMoFEMCreateMoFEM
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: DMMoFEM.cpp:114
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:325
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:290
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MortarContactProblem::LoadScale::lAmbda
static double lAmbda
Definition: MortarContactProblem.hpp:84
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::SmartPetscObj< DM >
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
NonlinearElasticElementInterface
Set of functions declaring elements and setting operators for generic element interface.
Definition: NonlinearElasticElementInterface.hpp:15
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
Generated by Doxygen 1.8.17 and hosted at University of Glasgow