v0.14.0
simple_contact.cpp
Go to the documentation of this file.
1 /** \file simple_contact.cpp
2  * \example simple_contact.cpp
3  *
4  * Implementation of simple contact between surfaces with matching meshes
5  *
6  **/
7 
8 /* MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
22 
23 #include <Mortar.hpp>
24 #include <Hooke.hpp>
25 
26 using namespace MoFEM;
27 using namespace std;
28 
29 
30 static char help[] = "\n";
31 double SimpleContactProblem::LoadScale::lAmbda = 1;
32 
33 int main(int argc, char *argv[]) {
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  "-my_order 1 \n"
52  "-my_order_lambda 1 \n"
53  "-my_order_contact 2 \n"
54  "-my_ho_levels_num 1 \n"
55  "-my_step_num 1 \n"
56  "-my_cn_value 1. \n"
57  "-my_r_value 1. \n"
58  "-my_alm_flag 0 \n";
59 
60  string param_file = "param_file.petsc";
61  if (!static_cast<bool>(ifstream(param_file))) {
62  std::ofstream file(param_file.c_str(), std::ios::ate);
63  if (file.is_open()) {
64  file << default_options;
65  file.close();
66  }
67  }
68 
69  enum contact_tests {
70  EIGHT_CUBE = 1,
71  FOUR_SEASONS = 2,
72  T_INTERFACE = 3,
73  PUNCH_TOP_AND_MID = 4,
74  PUNCH_TOP_ONLY = 5,
75  PLANE_AXI = 6,
76  ARC_THREE_SURF = 7,
77  SMILING_FACE = 8,
78  SMILING_FACE_CONVECT = 9,
79  WAVE_2D = 10,
80  WAVE_2D_ALM = 11,
81  LAST_TEST
82  };
83 
84  // Initialize MoFEM
85  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
86 
87  // Create mesh database
88  moab::Core mb_instance; // create database
89  moab::Interface &moab = mb_instance; // create interface to database
90 
91  try {
92  PetscBool flg_file;
93 
94  char mesh_file_name[255];
95  PetscInt order = 1;
96  PetscInt order_contact = 1;
97  PetscInt nb_ho_levels = 0;
98  PetscInt order_lambda = 1;
99  PetscReal r_value = 1.;
100  PetscReal cn_value = -1;
101  PetscInt nb_sub_steps = 1;
102  PetscBool is_partitioned = PETSC_FALSE;
103  PetscBool is_newton_cotes = PETSC_FALSE;
104  PetscInt test_num = 0;
105  PetscBool convect_pts = PETSC_FALSE;
106  PetscBool print_contact_state = PETSC_FALSE;
107  PetscBool alm_flag = PETSC_FALSE;
108  PetscBool wave_surf_flag = PETSC_FALSE;
109  PetscInt wave_dim = 2;
110  PetscInt wave_surf_block_id = 1;
111  PetscReal wave_length = 1.0;
112  PetscReal wave_ampl = 0.01;
113  PetscReal mesh_height = 1.0;
114 
115  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
116 
117  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
118  mesh_file_name, 255, &flg_file);
119 
120  CHKERR PetscOptionsInt("-my_order",
121  "approximation order of spatial positions", "", 1,
122  &order, PETSC_NULL);
123  CHKERR PetscOptionsInt(
124  "-my_order_contact",
125  "approximation order of spatial positions in contact interface", "", 1,
126  &order_contact, PETSC_NULL);
127  CHKERR PetscOptionsInt("-my_ho_levels_num", "number of higher order levels",
128  "", 0, &nb_ho_levels, PETSC_NULL);
129  CHKERR PetscOptionsInt("-my_order_lambda",
130  "approximation order of Lagrange multipliers", "", 1,
131  &order_lambda, PETSC_NULL);
132 
133  CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", nb_sub_steps,
134  &nb_sub_steps, PETSC_NULL);
135 
136  CHKERR PetscOptionsBool("-my_is_partitioned",
137  "set if mesh is partitioned (this result that each "
138  "process keeps only part of the mes",
139  "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
140 
141  CHKERR PetscOptionsReal("-my_cn_value", "default regularisation cn value",
142  "", 1., &cn_value, PETSC_NULL);
143 
144  CHKERR PetscOptionsBool("-my_is_newton_cotes",
145  "set if Newton-Cotes quadrature rules are used", "",
146  PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
147  CHKERR PetscOptionsInt("-my_test_num", "test number", "", 0, &test_num,
148  PETSC_NULL);
149  CHKERR PetscOptionsBool("-my_convect", "set to convect integration pts", "",
150  PETSC_FALSE, &convect_pts, PETSC_NULL);
151  CHKERR PetscOptionsBool("-my_print_contact_state",
152  "output number of active gp at every iteration", "",
153  PETSC_FALSE, &print_contact_state, PETSC_NULL);
154  CHKERR PetscOptionsBool("-my_alm_flag",
155  "if set use ALM, if not use C-function", "",
156  PETSC_FALSE, &alm_flag, PETSC_NULL);
157 
158  CHKERR PetscOptionsBool("-my_wave_surf",
159  "if set true, make one of the surfaces wavy", "",
160  PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
161  CHKERR PetscOptionsInt("-my_wave_surf_block_id",
162  "make wavy surface of the block with this id", "",
163  wave_surf_block_id, &wave_surf_block_id, PETSC_NULL);
164  CHKERR PetscOptionsInt("-my_wave_dim", "dimension (2 or 3)", "", wave_dim,
165  &wave_dim, PETSC_NULL);
166  CHKERR PetscOptionsReal("-my_wave_length", "profile wavelength", "",
167  wave_length, &wave_length, PETSC_NULL);
168  CHKERR PetscOptionsReal("-my_wave_ampl", "profile amplitude", "", wave_ampl,
169  &wave_ampl, PETSC_NULL);
170  CHKERR PetscOptionsReal("-my_mesh_height",
171  "vertical dimension of the mesh ", "", mesh_height,
172  &mesh_height, PETSC_NULL);
173 
174  ierr = PetscOptionsEnd();
175  CHKERRQ(ierr);
176 
177  // Check if mesh file was provided
178  if (flg_file != PETSC_TRUE) {
179  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
180  }
181 
182  if (is_partitioned == PETSC_TRUE) {
183  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
184  "Partitioned mesh is not supported");
185  }
186 
187  const char *option;
188  option = "";
189  CHKERR moab.load_file(mesh_file_name, 0, option);
190 
191  // Create MoFEM database and link it to MoAB
192  MoFEM::Core core(moab);
193  MoFEM::Interface &m_field = core;
194 
195  MeshsetsManager *mmanager_ptr;
196  CHKERR m_field.getInterface(mmanager_ptr);
197  CHKERR mmanager_ptr->printDisplacementSet();
198  CHKERR mmanager_ptr->printForceSet();
199  // print block sets with materials
200  CHKERR mmanager_ptr->printMaterialsSet();
201 
202  auto add_prism_interface = [&](Range &contact_prisms, Range &master_tris,
203  Range &slave_tris,
204  std::vector<BitRefLevel> &bit_levels) {
205  MoFEMFunctionBegin;
206  PrismInterface *interface;
207  CHKERR m_field.getInterface(interface);
208 
209  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
210  if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
211  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
212  cit->getName().c_str(), cit->getMeshsetId());
213  EntityHandle cubit_meshset = cit->getMeshset();
214 
215  // get tet entities from back bit_level
216  EntityHandle ref_level_meshset;
217  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
218  CHKERR m_field.getInterface<BitRefManager>()
219  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
220  BitRefLevel().set(), MBTET,
221  ref_level_meshset);
222  CHKERR m_field.getInterface<BitRefManager>()
223  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
224  BitRefLevel().set(), MBPRISM,
225  ref_level_meshset);
226 
227  // get faces and tets to split
228  CHKERR interface->getSides(cubit_meshset, bit_levels.back(), true, 0);
229  // set new bit level
230  bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
231  // split faces and tets
232  CHKERR interface->splitSides(ref_level_meshset, bit_levels.back(),
233  cubit_meshset, true, true, 0);
234  // clean meshsets
235  CHKERR moab.delete_entities(&ref_level_meshset, 1);
236 
237  // update cubit meshsets
238  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
239  EntityHandle cubit_meshset = ciit->meshset;
240  CHKERR m_field.getInterface<BitRefManager>()
241  ->updateMeshsetByEntitiesChildren(
242  cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
243  true);
244  CHKERR m_field.getInterface<BitRefManager>()
245  ->updateMeshsetByEntitiesChildren(cubit_meshset,
246  bit_levels.back(),
247  cubit_meshset, MBEDGE, true);
248  CHKERR m_field.getInterface<BitRefManager>()
249  ->updateMeshsetByEntitiesChildren(cubit_meshset,
250  bit_levels.back(),
251  cubit_meshset, MBTRI, true);
252  CHKERR m_field.getInterface<BitRefManager>()
253  ->updateMeshsetByEntitiesChildren(cubit_meshset,
254  bit_levels.back(),
255  cubit_meshset, MBTET, true);
256  }
257 
258  CHKERR m_field.getInterface<BitRefManager>()->shiftRightBitRef(1);
259  bit_levels.pop_back();
260  }
261  }
262 
263  EntityHandle meshset_prisms;
264  CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
265  CHKERR m_field.getInterface<BitRefManager>()
266  ->getEntitiesByTypeAndRefLevel(bit_levels.back(), BitRefLevel().set(),
267  MBPRISM, meshset_prisms);
268  CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
269  CHKERR moab.delete_entities(&meshset_prisms, 1);
270 
271  EntityHandle tri;
272  for (Range::iterator pit = contact_prisms.begin();
273  pit != contact_prisms.end(); pit++) {
274  CHKERR moab.side_element(*pit, 2, 3, tri);
275  master_tris.insert(tri);
276  CHKERR moab.side_element(*pit, 2, 4, tri);
277  slave_tris.insert(tri);
278  }
279 
280  MoFEMFunctionReturn(0);
281  };
282 
283  auto make_wavy_surface = [&](int block_id, int dim, double lambda,
284  double delta, double height) {
285  MoFEMFunctionBegin;
286  Range all_tets, all_nodes;
287  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
288  if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
289  const int id = bit->getMeshsetId();
290  Range tets;
291  if (id == block_id) {
292  CHKERR m_field.get_moab().get_entities_by_dimension(
293  bit->getMeshset(), 3, tets, true);
294  all_tets.merge(tets);
295  }
296  }
297  }
298  CHKERR m_field.get_moab().get_connectivity(all_tets, all_nodes);
299  double coords[3];
300  for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
301  nit++) {
302  CHKERR moab.get_coords(&*nit, 1, coords);
303  double x = coords[0];
304  double y = coords[1];
305  double z = coords[2];
306  double coef = (height + z) / height;
307  switch (dim) {
308  case 2:
309  coords[2] -= coef * delta * (1. - cos(2. * M_PI * x / lambda));
310  break;
311  case 3:
312  coords[2] -=
313  coef * delta *
314  (1. - cos(2. * M_PI * x / lambda) * cos(2. * M_PI * y / lambda));
315  break;
316  default:
317  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
318  "Wrong dimension = %d", dim);
319  }
320 
321  CHKERR moab.set_coords(&*nit, 1, coords);
322  }
323  MoFEMFunctionReturn(0);
324  };
325 
326  auto set_contact_order = [&](Range &contact_prisms, int order_contact,
327  int nb_ho_levels) {
328  MoFEMFunctionBegin;
329  Range contact_tris, contact_edges;
330  CHKERR moab.get_adjacencies(contact_prisms, 2, false, contact_tris,
331  moab::Interface::UNION);
332  contact_tris = contact_tris.subset_by_type(MBTRI);
333  CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
334  moab::Interface::UNION);
335  Range ho_ents;
336  ho_ents.merge(contact_tris);
337  ho_ents.merge(contact_edges);
338  for (int ll = 0; ll < nb_ho_levels; ll++) {
339  Range ents, verts, tets;
340  CHKERR moab.get_connectivity(ho_ents, verts, true);
341  CHKERR moab.get_adjacencies(verts, 3, false, tets,
342  moab::Interface::UNION);
343  tets = tets.subset_by_type(MBTET);
344  for (auto d : {1, 2}) {
345  CHKERR moab.get_adjacencies(tets, d, false, ents,
346  moab::Interface::UNION);
347  }
348  ho_ents = unite(ho_ents, ents);
349  ho_ents = unite(ho_ents, tets);
350  }
351 
352  CHKERR m_field.set_field_order(ho_ents, "SPATIAL_POSITION",
353  order_contact);
354 
355  MoFEMFunctionReturn(0);
356  };
357 
358  Range contact_prisms, master_tris, slave_tris;
359  std::vector<BitRefLevel> bit_levels;
360 
361  bit_levels.push_back(BitRefLevel().set(0));
362  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
363  0, 3, bit_levels.back());
364 
365  CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
366  bit_levels);
367 
368  if (wave_surf_flag) {
369  CHKERR make_wavy_surface(wave_surf_block_id, wave_dim, wave_length,
370  wave_ampl, mesh_height);
371  }
372 
373  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
374  MB_TAG_SPARSE, MF_ZERO);
375 
376  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
377  3, MB_TAG_SPARSE, MF_ZERO);
378 
379  CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
380  MB_TAG_SPARSE, MF_ZERO);
381 
382  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
383  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
384  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
385  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
386  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
387 
388  // Declare problem add entities (by tets) to the field
389  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
390  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
391  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
392  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
393  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
394 
395  CHKERR m_field.add_ents_to_field_by_type(slave_tris, MBTRI, "LAGMULT");
396  CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
397  CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
398  CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
399 
400  if (order_contact > order) {
401  CHKERR set_contact_order(contact_prisms, order_contact, nb_ho_levels);
402  }
403 
404  // build field
405  CHKERR m_field.build_fields();
406 
407  // Projection on "x" field
408  {
409  Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
410  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
411  }
412  // Projection on "X" field
413  {
414  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
415  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
416  }
417 
418  // Add elastic element
419  boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(new Hooke<adouble>());
420  boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>());
421  NonlinearElasticElement elastic(m_field, 2);
422  CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
423  CHKERR elastic.addElement("ELASTIC", "SPATIAL_POSITION");
424 
425  CHKERR elastic.setOperators("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
426  false, false);
427 
428  auto make_contact_element = [&]() {
429  return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
430  m_field);
431  };
432 
433  auto make_convective_master_element = [&]() {
434  return boost::make_shared<
435  SimpleContactProblem::ConvectMasterContactElement>(
436  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
437  };
438 
439  auto make_convective_slave_element = [&]() {
440  return boost::make_shared<
441  SimpleContactProblem::ConvectSlaveContactElement>(
442  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
443  };
444 
445  auto make_contact_common_data = [&]() {
446  return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
447  m_field);
448  };
449 
450  auto get_contact_rhs = [&](auto contact_problem, auto make_element,
451  bool is_alm = false) {
452  auto fe_rhs_simple_contact = make_element();
453  auto common_data_simple_contact = make_contact_common_data();
454  if (print_contact_state) {
455  fe_rhs_simple_contact->contactStateVec =
456  common_data_simple_contact->gaussPtsStateVec;
457  }
458  contact_problem->setContactOperatorsRhs(
459  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
460  "LAGMULT", is_alm);
461  return fe_rhs_simple_contact;
462  };
463 
464  auto get_master_traction_rhs = [&](auto contact_problem, auto make_element,
465  bool is_alm = false) {
466  auto fe_rhs_simple_contact = make_element();
467  auto common_data_simple_contact = make_contact_common_data();
468  contact_problem->setMasterForceOperatorsRhs(
469  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
470  "LAGMULT", is_alm);
471  return fe_rhs_simple_contact;
472  };
473 
474  auto get_master_traction_lhs = [&](auto contact_problem, auto make_element,
475  bool is_alm = false) {
476  auto fe_lhs_simple_contact = make_element();
477  auto common_data_simple_contact = make_contact_common_data();
478  contact_problem->setMasterForceOperatorsLhs(
479  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
480  "LAGMULT", is_alm);
481  return fe_lhs_simple_contact;
482  };
483 
484  auto get_contact_lhs = [&](auto contact_problem, auto make_element,
485  bool is_alm = false) {
486  auto fe_lhs_simple_contact = make_element();
487  auto common_data_simple_contact = make_contact_common_data();
488  contact_problem->setContactOperatorsLhs(
489  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
490  "LAGMULT", is_alm);
491  return fe_lhs_simple_contact;
492  };
493 
494  auto cn_value_ptr = boost::make_shared<double>(cn_value);
495  auto contact_problem = boost::make_shared<SimpleContactProblem>(
496  m_field, cn_value_ptr, is_newton_cotes);
497 
498  // add fields to the global matrix by adding the element
499  contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
500  "LAGMULT", contact_prisms);
501 
502  contact_problem->addPostProcContactElement(
503  "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
504  "MESH_NODE_POSITIONS", slave_tris);
505 
506  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "SPATIAL_POSITION");
507 
508  // Add spring boundary condition applied on surfaces.
509  CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
510  "MESH_NODE_POSITIONS");
511 
512  // build finite elemnts
513  CHKERR m_field.build_finite_elements();
514 
515  // build adjacencies
516  CHKERR m_field.build_adjacencies(bit_levels.back());
517 
518  // define problems
519  CHKERR m_field.add_problem("CONTACT_PROB");
520 
521  // set refinement level for problem
522  CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
523  bit_levels.back());
524 
525  DMType dm_name = "DMMOFEM";
526  CHKERR DMRegister_MoFEM(dm_name);
527 
528  SmartPetscObj<DM> dm;
529  dm = createSmartDM(m_field.get_comm(), dm_name);
530 
531  // create dm instance
532  CHKERR DMSetType(dm, dm_name);
533 
534  // set dm datastruture which created mofem datastructures
535  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_levels.back());
536  CHKERR DMSetFromOptions(dm);
537  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
538  // add elements to dm
539  CHKERR DMMoFEMAddElement(dm, "CONTACT_ELEM");
540  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
541  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
542  CHKERR DMMoFEMAddElement(dm, "SPRING");
543  CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
544 
545  CHKERR DMSetUp(dm);
546 
547  // Vector of DOFs and the RHS
548  auto D = smartCreateDMVector(dm);
549  auto F = smartVectorDuplicate(D);
550 
551  // Stiffness matrix
552  auto Aij = smartCreateDMMatrix(dm);
553 
554  CHKERR VecZeroEntries(D);
555  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
556  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
557  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
558 
559  CHKERR VecZeroEntries(F);
560  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
561  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
562 
563  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
564  CHKERR MatZeroEntries(Aij);
565 
566  // Dirichlet BC
567  boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
568  boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
569  m_field, "SPATIAL_POSITION", Aij, D, F));
570 
571  dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
572  dirichlet_bc_ptr->snes_x = D;
573 
574  // Assemble pressure and traction forces
575  boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
576  CHKERR MetaNeumannForces::setMomentumFluxOperators(
577  m_field, neumann_forces, NULL, "SPATIAL_POSITION");
578 
579  boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
580  neumann_forces.begin();
581  for (; mit != neumann_forces.end(); mit++) {
582  mit->second->methodsOp.push_back(new SimpleContactProblem::LoadScale());
583  CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
584  &mit->second->getLoopFe(), NULL, NULL);
585  }
586 
587  // Implementation of spring element
588  // Create new instances of face elements for springs
589  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
590  new FaceElementForcesAndSourcesCore(m_field));
591  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
592  new FaceElementForcesAndSourcesCore(m_field));
593 
594  CHKERR MetaSpringBC::setSpringOperators(
595  m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
596  "MESH_NODE_POSITIONS");
597 
598  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
599  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
600  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
601  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
602  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL,
603  dirichlet_bc_ptr.get(), NULL);
604  if (convect_pts == PETSC_TRUE) {
605  CHKERR DMMoFEMSNESSetFunction(
606  dm, "CONTACT_ELEM",
607  get_contact_rhs(contact_problem, make_convective_master_element),
608  PETSC_NULL, PETSC_NULL);
609  CHKERR DMMoFEMSNESSetFunction(
610  dm, "CONTACT_ELEM",
611  get_master_traction_rhs(contact_problem,
612  make_convective_slave_element),
613  PETSC_NULL, PETSC_NULL);
614  } else {
615  CHKERR DMMoFEMSNESSetFunction(
616  dm, "CONTACT_ELEM",
617  get_contact_rhs(contact_problem, make_contact_element, alm_flag),
618  PETSC_NULL, PETSC_NULL);
619  CHKERR DMMoFEMSNESSetFunction(
620  dm, "CONTACT_ELEM",
621  get_master_traction_rhs(contact_problem, make_contact_element,
622  alm_flag),
623  PETSC_NULL, PETSC_NULL);
624  }
625 
626  CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", &elastic.getLoopFeRhs(),
627  PETSC_NULL, PETSC_NULL);
628  CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
629  PETSC_NULL);
630  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, PETSC_NULL, PETSC_NULL,
631  dirichlet_bc_ptr.get());
632 
633  boost::shared_ptr<FEMethod> fe_null;
634  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
635  fe_null);
636  if (convect_pts == PETSC_TRUE) {
637  CHKERR DMMoFEMSNESSetJacobian(
638  dm, "CONTACT_ELEM",
639  get_contact_lhs(contact_problem, make_convective_master_element),
640  PETSC_NULL, PETSC_NULL);
641  CHKERR DMMoFEMSNESSetJacobian(
642  dm, "CONTACT_ELEM",
643  get_master_traction_lhs(contact_problem,
644  make_convective_slave_element),
645  PETSC_NULL, PETSC_NULL);
646  } else {
647  CHKERR DMMoFEMSNESSetJacobian(
648  dm, "CONTACT_ELEM",
649  get_contact_lhs(contact_problem, make_contact_element, alm_flag),
650  PETSC_NULL, PETSC_NULL);
651  CHKERR DMMoFEMSNESSetJacobian(
652  dm, "CONTACT_ELEM",
653  get_master_traction_lhs(contact_problem, make_contact_element,
654  alm_flag),
655  PETSC_NULL, PETSC_NULL);
656  }
657  CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", &elastic.getLoopFeLhs(),
658  PETSC_NULL, PETSC_NULL);
659  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
660  PETSC_NULL);
661  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, fe_null,
662  dirichlet_bc_ptr);
663 
664  if (test_num) {
665  char testing_options[] = "-ksp_type fgmres "
666  "-pc_type lu "
667  "-pc_factor_mat_solver_type mumps "
668  "-snes_type newtonls "
669  "-snes_linesearch_type basic "
670  "-snes_max_it 10 "
671  "-snes_atol 1e-8 "
672  "-snes_rtol 1e-8 ";
673  CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
674  }
675 
676  auto snes = MoFEM::createSNES(m_field.get_comm());
677  CHKERR SNESSetDM(snes, dm);
678  SNESConvergedReason snes_reason;
679  SnesCtx *snes_ctx;
680  // create snes nonlinear solver
681  {
682  CHKERR SNESSetDM(snes, dm);
683  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
684  CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
685  CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
686  CHKERR SNESSetFromOptions(snes);
687  }
688 
689  PostProcVolumeOnRefinedMesh post_proc(m_field);
690  // Add operators to the elements, starting with some generic
691  CHKERR post_proc.generateReferenceElementMesh();
692  CHKERR post_proc.addFieldValuesPostProc("SPATIAL_POSITION");
693  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
694  CHKERR post_proc.addFieldValuesGradientPostProc("SPATIAL_POSITION");
695 
696  std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
697  elastic.setOfBlocks.begin();
698  for (; sit != elastic.setOfBlocks.end(); sit++) {
699  post_proc.getOpPtrVector().push_back(new PostProcStress(
700  post_proc.postProcMesh, post_proc.mapGaussPts, "SPATIAL_POSITION",
701  sit->second, post_proc.commonData));
702  }
703 
704  for (int ss = 0; ss != nb_sub_steps; ++ss) {
705  SimpleContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
706  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
707  SimpleContactProblem::LoadScale::lAmbda);
708 
709  CHKERR SNESSolve(snes, PETSC_NULL, D);
710 
711  CHKERR SNESGetConvergedReason(snes, &snes_reason);
712 
713  int its;
714  CHKERR SNESGetIterationNumber(snes, &its);
715  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Number of Newton iterations = %D\n",
716  its);
717 
718  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
719  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
720  }
721 
722  // save on mesh
723  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
724  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
725  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
726 
727  PetscPrintf(PETSC_COMM_WORLD, "Loop post proc\n");
728  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
729 
730  elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
731  elastic.getLoopFeEnergy().eNergy = 0;
732  PetscPrintf(PETSC_COMM_WORLD, "Loop energy\n");
733  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeEnergy());
734  // Print elastic energy
735  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %9.9f\n",
736  elastic.getLoopFeEnergy().eNergy);
737 
738  {
739  string out_file_name;
740  std::ostringstream stm;
741  stm << "out"
742  << ".h5m";
743  out_file_name = stm.str();
744  CHKERR
745  PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n", out_file_name.c_str());
746  CHKERR post_proc.postProcMesh.write_file(out_file_name.c_str(), "MOAB",
747  "PARALLEL=WRITE_PART");
748  }
749 
750  // moab_instance
751  moab::Core mb_post; // create database
752  moab::Interface &moab_proc = mb_post; // create interface to database
753 
754  auto common_data_simple_contact = make_contact_common_data();
755 
756  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
757  fe_post_proc_simple_contact;
758  if (convect_pts == PETSC_TRUE) {
759  fe_post_proc_simple_contact = make_convective_master_element();
760  } else {
761  fe_post_proc_simple_contact = make_contact_element();
762  }
763 
764  contact_problem->setContactOperatorsForPostProc(
765  fe_post_proc_simple_contact, common_data_simple_contact, m_field,
766  "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag);
767 
768  mb_post.delete_mesh();
769 
770  CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
771  CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
772 
773  CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_ELEM",
774  fe_post_proc_simple_contact);
775 
776  std::array<double, 2> nb_gauss_pts;
777  std::array<double, 2> contact_area;
778 
779  auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
780  MoFEMFunctionBegin;
781  CHKERR VecAssemblyBegin(vec);
782  CHKERR VecAssemblyEnd(vec);
783  const double *array;
784  CHKERR VecGetArrayRead(vec, &array);
785  if (m_field.get_comm_rank() == 0) {
786  for (int i : {0, 1})
787  data[i] = array[i];
788  }
789  CHKERR VecRestoreArrayRead(vec, &array);
790  MoFEMFunctionReturn(0);
791  };
792 
793  CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
794  nb_gauss_pts);
795  CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
796  contact_area);
797 
798  if (m_field.get_comm_rank() == 0) {
799  PetscPrintf(PETSC_COMM_SELF, "Active gauss pts: %d out of %d\n",
800  (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
801 
802  PetscPrintf(PETSC_COMM_SELF, "Active contact area: %9.9f out of %9.9f\n",
803  contact_area[0], contact_area[1]);
804  }
805 
806  if (test_num) {
807  double expected_energy, expected_contact_area;
808  int expected_nb_gauss_pts;
809  constexpr double eps = 1e-6;
810  switch (test_num) {
811  case EIGHT_CUBE:
812  expected_energy = 3.0e-04;
813  expected_contact_area = 3.0;
814  expected_nb_gauss_pts = 576;
815  break;
816  case FOUR_SEASONS:
817  expected_energy = 1.2e-01;
818  expected_contact_area = 106.799036701;
819  expected_nb_gauss_pts = 672;
820  break;
821  case T_INTERFACE:
822  expected_energy = 3.0e-04;
823  expected_contact_area = 1.75;
824  expected_nb_gauss_pts = 336;
825  break;
826  case PUNCH_TOP_AND_MID:
827  expected_energy = 3.125e-04;
828  expected_contact_area = 0.25;
829  expected_nb_gauss_pts = 84;
830  break;
831  case PUNCH_TOP_ONLY:
832  expected_energy = 0.000096432;
833  expected_contact_area = 0.25;
834  expected_nb_gauss_pts = 336;
835  break;
836  case PLANE_AXI:
837  expected_energy = 0.000438889;
838  expected_contact_area = 0.784409608;
839  expected_nb_gauss_pts = 300;
840  break;
841  case ARC_THREE_SURF:
842  expected_energy = 0.002573411;
843  expected_contact_area = 2.831455633;
844  expected_nb_gauss_pts = 228;
845  break;
846  case SMILING_FACE:
847  expected_energy = 0.000733553;
848  expected_contact_area = 3.0;
849  expected_nb_gauss_pts = 144;
850  break;
851  case SMILING_FACE_CONVECT:
852  expected_energy = 0.000733621;
853  expected_contact_area = 3.0;
854  expected_nb_gauss_pts = 144;
855  break;
856  case WAVE_2D:
857  expected_energy = 0.008537863;
858  expected_contact_area = 0.125;
859  expected_nb_gauss_pts = 384;
860  break;
861  case WAVE_2D_ALM:
862  expected_energy = 0.008538894;
863  expected_contact_area = 0.125;
864  expected_nb_gauss_pts = 384;
865  break;
866  default:
867  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
868  "Unknown test number %d", test_num);
869  }
870  if (std::abs(elastic.getLoopFeEnergy().eNergy - expected_energy) > eps) {
871  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
872  "Wrong energy %6.4e != %6.4e", expected_energy,
873  elastic.getLoopFeEnergy().eNergy);
874  }
875  if (m_field.get_comm_rank() == 0) {
876  if ((int)nb_gauss_pts[0] != expected_nb_gauss_pts) {
877  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
878  "Wrong number of active gauss pts %d != %d",
879  expected_nb_gauss_pts, (int)nb_gauss_pts[0]);
880  }
881  if (std::abs(contact_area[0] - expected_contact_area) > eps) {
882  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
883  "Wrong active contact area %6.4e != %6.4e",
884  expected_contact_area, contact_area[0]);
885  }
886  }
887  }
888 
889  {
890  string out_file_name;
891  std::ostringstream strm;
892  strm << "out_contact_integ_pts"
893  << ".h5m";
894  out_file_name = strm.str();
895  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
896  out_file_name.c_str());
897  CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
898  "PARALLEL=WRITE_PART");
899  }
900 
901  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
902  new PostProcFaceOnRefinedMesh(m_field));
903 
904  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *post_proc_contact_ptr, false,
905  false);
906  CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
907  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("LAGMULT");
908  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("SPATIAL_POSITION");
909  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
910 
911  CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
912  post_proc_contact_ptr);
913 
914  {
915  string out_file_name;
916  std::ostringstream stm;
917  stm << "out_contact"
918  << ".h5m";
919  out_file_name = stm.str();
920  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
921  out_file_name.c_str());
922  CHKERR post_proc_contact_ptr->postProcMesh.write_file(
923  out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
924  }
925  }
926  CATCH_ERRORS;
927 
928  // finish work cleaning memory, getting statistics, etc
929  MoFEM::Core::Finalize();
930 
931  return 0;
932 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:699
MortarCtestFunctions::LAST_TEST
@ LAST_TEST
Definition: MortarCtestFunctions.hpp:37
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
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
H1
@ H1
continuous field
Definition: definitions.h:85
EntityHandle
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
SimpleContactProblem::ConvectSlaveContactElement
Element used to integrate on master surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:206
MoFEM::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition: PrismInterface.hpp:23
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::createSNES
auto createSNES(MPI_Comm comm)
Definition: PetscSmartObj.hpp:255
Mortar.hpp
MoFEM::createSmartDM
DEPRECATED auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Definition: PetscSmartObj.hpp:149
MortarCtestFunctions::PUNCH_TOP_AND_MID
@ PUNCH_TOP_AND_MID
Definition: MortarCtestFunctions.hpp:27
MortarCtestFunctions::contact_tests
contact_tests
Definition: MortarCtestFunctions.hpp:24
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
SimpleContactProblem::LoadScale::lAmbda
static double lAmbda
Definition: SimpleContact.hpp:37
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
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
out_file_name
char out_file_name[255]
Definition: initial_diffusion.cpp:53
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MortarCtestFunctions::WAVE_2D_ALM
@ WAVE_2D_ALM
Definition: MortarCtestFunctions.hpp:33
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
help
static char help[]
Definition: simple_contact.cpp:30
MortarCtestFunctions::SMILING_FACE
@ SMILING_FACE
Definition: MortarCtestFunctions.hpp:30
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
main
int main(int argc, char *argv[])
Definition: simple_contact.cpp:33
PostProcStress
Definition: PostProcStresses.hpp:17
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
Hooke.hpp
Implementation of linear elastic material.
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
PostProcTemplateVolumeOnRefinedMesh::commonData
CommonData commonData
Definition: PostProcOnRefMesh.hpp:287
MortarCtestFunctions::T_INTERFACE
@ T_INTERFACE
Definition: MortarCtestFunctions.hpp:26
MoFEM::PrismInterface::splitSides
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...
Definition: PrismInterface.cpp:519
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
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::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1094
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::smartCreateDMMatrix
DEPRECATED auto smartCreateDMMatrix(DM dm)
Definition: DMMoFEM.hpp:1092
MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition: PetscSmartObj.hpp:230
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::DMMoFEMSNESSetJacobian
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: DMMoFEM.cpp:759
MoFEM::smartCreateDMVector
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1107
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:195
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
NonlinearElasticElement::addElement
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition: NonLinearElasticElement.cpp:1120
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
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
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
NonlinearElasticElement::setOperators
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Definition: NonLinearElasticElement.cpp:1153
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MortarCtestFunctions::WAVE_2D
@ WAVE_2D
Definition: MortarCtestFunctions.hpp:32
SimpleContactProblem::LoadScale
Definition: SimpleContact.hpp:35
MoFEM::SnesRhs
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:27
Range
MetaSpringBC::setSpringOperators
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
Definition: SpringElement.cpp:1178
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:290
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
std
Definition: enable_if.hpp:5
MoFEM::PrismInterface::getSides
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...
Definition: PrismInterface.cpp:56
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MortarCtestFunctions::EIGHT_CUBE
@ EIGHT_CUBE
Definition: MortarCtestFunctions.hpp:25
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
SimpleContactProblem::ConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:181
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
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::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MortarCtestFunctions::SMILING_FACE_CONVECT
@ SMILING_FACE_CONVECT
Definition: MortarCtestFunctions.hpp:31
MoFEM::CoreInterface::set_field_order
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.
MoFEM::SmartPetscObj< DM >
MoFEM::SnesMat
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:139
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MortarCtestFunctions::PUNCH_TOP_ONLY
@ PUNCH_TOP_ONLY
Definition: MortarCtestFunctions.hpp:28
MoFEM::CoreInterface::add_field
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.
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
NonlinearElasticElement::setBlocks
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
Definition: NonLinearElasticElement.cpp:1086
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MoFEM::DMMoFEMSNESSetFunction
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: DMMoFEM.cpp:718
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:153
Generated by Doxygen 1.8.17 and hosted at University of Glasgow