v0.14.0
mortar_contact_thermal.cpp
Go to the documentation of this file.
1 /** \file mortar_contact_thermal.cpp
2  * \example mortar_contact_thermal.cpp
3  *
4  * Implementation of mortar contact between surfaces with non-matching meshes
5  * taking into account internal stress resulting from the thermal expansion
6  *
7  **/
8 
9 /* MoFEM is free software: you can redistribute it and/or modify it under
10  * the terms of the GNU Lesser General Public License as published by the
11  * Free Software Foundation, either version 3 of the License, or (at your
12  * option) any later version.
13  *
14  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  * License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #include <Mortar.hpp>
24 #include <HookeInternalStressElement.hpp>
25 
26 static char help[] = "\n";
27 double MortarContactProblem::LoadScale::lAmbda = 1;
28 using BlockData = NonlinearElasticElement::BlockData;
29 using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide;
30 
31 struct VolRule {
32  int operator()(int, int, int order) const { return 2 * (order - 1); }
33 };
34 
35 int main(int argc, char *argv[]) {
36 
37  const string default_options = "-ksp_type fgmres \n"
38  "-pc_type lu \n"
39  "-pc_factor_mat_solver_type mumps \n"
40  "-mat_mumps_icntl_13 1 \n"
41  "-mat_mumps_icntl_14 800 \n"
42  "-mat_mumps_icntl_20 0 \n"
43  "-mat_mumps_icntl_24 1 \n"
44  "-snes_type newtonls \n"
45  "-snes_linesearch_type basic \n"
46  "-snes_divergence_tolerance 0 \n"
47  "-snes_max_it 50 \n"
48  "-snes_atol 1e-8 \n"
49  "-snes_rtol 1e-10 \n"
50  "-snes_monitor \n"
51  "-ksp_monitor \n"
52  "-snes_converged_reason \n"
53  "-order 1 \n"
54  "-order_lambda 1 \n"
55  "-order_contact 2 \n"
56  "-ho_levels_num 1 \n"
57  "-step_num 1 \n"
58  "-cn_value 1. \n"
59  "-r_value 1. \n"
60  "-alm_flag 0 \n"
61  "-eigen_pos_flag 0 \n";
62 
63  string param_file = "param_file.petsc";
64  if (!static_cast<bool>(ifstream(param_file))) {
65  std::ofstream file(param_file.c_str(), std::ios::ate);
66  if (file.is_open()) {
67  file << default_options;
68  file.close();
69  }
70  }
71 
72  // Initialize MoFEM
73  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
74 
75  // Create mesh database
76  moab::Core mb_instance; // create database
77  moab::Interface &moab = mb_instance; // create interface to database
78 
79  try {
80  PetscBool flg_file;
81  PetscBool flg_file_out;
82 
83  char mesh_file_name[255];
84  char output_mesh_name[255];
85  PetscInt order = 1;
86  PetscInt order_contact = 1;
87  PetscInt nb_ho_levels = 0;
88  PetscInt order_lambda = 1;
89  PetscReal r_value = 1.;
90  PetscReal cn_value = -1;
91  PetscInt nb_sub_steps = 1;
92  PetscBool is_partitioned = PETSC_FALSE;
93  PetscBool is_newton_cotes = PETSC_FALSE;
94  PetscInt test_num = 0;
95  PetscBool convect_pts = PETSC_FALSE;
96  PetscBool print_contact_state = PETSC_FALSE;
97  PetscBool alm_flag = PETSC_FALSE;
98  PetscBool eigen_pos_flag = PETSC_FALSE;
99 
100  PetscReal thermal_expansion_coef = 1e-5;
101  PetscReal init_temp = 250.0;
102  PetscReal scale_factor = 1.0;
103  PetscBool analytical_input = PETSC_TRUE;
104  char stress_tag_name[255];
105  PetscBool flg_tag_name;
106  PetscBool save_mean_stress = PETSC_TRUE;
107  PetscBool ignore_contact = PETSC_FALSE;
108  PetscBool ignore_pressure = PETSC_FALSE;
109 
110  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
111 
112  CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
113  mesh_file_name, 255, &flg_file);
114  CHKERR PetscOptionsString("-output_mesh_file", "output mesh file name", "",
115  "mesh.h5m", output_mesh_name, 255, &flg_file_out);
116 
117  CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
118  "", 1, &order, PETSC_NULL);
119  CHKERR PetscOptionsInt(
120  "-order_contact",
121  "approximation order of spatial positions in contact interface", "", 1,
122  &order_contact, PETSC_NULL);
123  CHKERR PetscOptionsInt("-ho_levels_num", "number of higher order levels",
124  "", 0, &nb_ho_levels, PETSC_NULL);
125  CHKERR PetscOptionsInt("-order_lambda",
126  "approximation order of Lagrange multipliers", "", 1,
127  &order_lambda, PETSC_NULL);
128 
129  CHKERR PetscOptionsInt("-step_num", "number of steps", "", nb_sub_steps,
130  &nb_sub_steps, PETSC_NULL);
131 
132  CHKERR PetscOptionsBool("-is_partitioned",
133  "set if mesh is partitioned (this result that each "
134  "process keeps only part of the mes",
135  "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
136 
137  CHKERR PetscOptionsReal("-cn_value", "default regularisation cn value", "",
138  1., &cn_value, PETSC_NULL);
139 
140  CHKERR PetscOptionsBool("-is_newton_cotes",
141  "set if Newton-Cotes quadrature rules are used", "",
142  PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
143  CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
144  PETSC_NULL);
145  CHKERR PetscOptionsBool("-convect", "set to convect integration pts", "",
146  PETSC_FALSE, &convect_pts, PETSC_NULL);
147  CHKERR PetscOptionsBool("-print_contact_state",
148  "output number of active gp at every iteration", "",
149  PETSC_FALSE, &print_contact_state, PETSC_NULL);
150  CHKERR PetscOptionsBool("-alm_flag",
151  "if set use ALM, if not use C-function", "",
152  PETSC_FALSE, &alm_flag, PETSC_NULL);
153  CHKERR PetscOptionsBool("-eigen_pos_flag",
154  "if set use eigen spatial positions are taken into "
155  "account for predeformed configuration",
156  "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
157 
158  CHKERR PetscOptionsReal("-scale_factor", "scale factor", "", scale_factor,
159  &scale_factor, PETSC_NULL);
160 
161  CHKERR PetscOptionsBool("-ignore_contact", "if set true, ignore contact",
162  "", PETSC_FALSE, &ignore_contact, PETSC_NULL);
163  CHKERR PetscOptionsBool("-ignore_pressure", "if set true, ignore pressure",
164  "", PETSC_FALSE, &ignore_pressure, PETSC_NULL);
165  CHKERR PetscOptionsBool("-analytical_input",
166  "if set true, use analytical strain", "",
167  PETSC_TRUE, &analytical_input, PETSC_NULL);
168  CHKERR PetscOptionsBool("-save_mean_stress",
169  "if set true, save mean stress", "", PETSC_TRUE,
170  &save_mean_stress, PETSC_NULL);
171  CHKERR PetscOptionsString("-stress_tag_name", "stress tag name file name",
172  "", "INTERNAL_STRESS", stress_tag_name, 255,
173  &flg_tag_name);
174  CHKERR PetscOptionsReal(
175  "-thermal_expansion_coef", "thermal expansion coef ", "",
176  thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
177  CHKERR PetscOptionsReal("-init_temp", "init_temp ", "", init_temp,
178  &init_temp, PETSC_NULL);
179 
180  ierr = PetscOptionsEnd();
181  CHKERRQ(ierr);
182 
183  // Check if mesh file was provided
184  if (flg_file != PETSC_TRUE) {
185  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
186  }
187 
188  if (is_partitioned) {
189  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
190  "Partitioned mesh is not supported");
191  }
192 
193  const char *option;
194  option = "";
195  CHKERR moab.load_file(mesh_file_name, 0, option);
196  Version file_ver;
197  CHKERR MoFEM::Interface::getFileVersion(moab, file_ver);
198  MOFEM_LOG("WORLD", Sev::inform) << "File version " << file_ver.strVersion();
199 
200  // Create MoFEM database and link it to MoAB
201  MoFEM::Core core(moab);
202  MoFEM::Interface &m_field = core;
203 
204  std::vector<BitRefLevel> bit_levels;
205  bit_levels.push_back(BitRefLevel().set(0));
206  auto bit_ref_manager = m_field.getInterface<BitRefManager>();
207  CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
208 
209  auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
210  MoFEMFunctionBegin;
211  auto prism_interface = m_field.getInterface<PrismInterface>();
212 
213  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
214  if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
215  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
216  cit->getName().c_str(), cit->getMeshsetId());
217  EntityHandle cubit_meshset = cit->getMeshset();
218 
219  // get tet entities from back bit_level
220  EntityHandle ref_level_meshset;
221  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
222  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
223  bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
224  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
225  bit_levels.back(), BitRefLevel().set(), MBPRISM,
226  ref_level_meshset);
227 
228  // get faces and tets to split
229  CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
230  true, 0);
231  // set new bit level
232  bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
233  // split faces and tets
234  CHKERR prism_interface->splitSides(ref_level_meshset,
235  bit_levels.back(), cubit_meshset,
236  true, true, 0);
237  // clean meshsets
238  CHKERR moab.delete_entities(&ref_level_meshset, 1);
239 
240  // update cubit meshsets
241  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
242  EntityHandle cubit_meshset = ciit->meshset;
243  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
244  cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
245  true);
246  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
247  cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE, true);
248  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
249  cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
250  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
251  cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
252  }
253 
254  CHKERR bit_ref_manager->shiftRightBitRef(1);
255  bit_levels.pop_back();
256  }
257  }
258 
259  MoFEMFunctionReturn(0);
260  };
261 
262  auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
263  Range &simple_contact_prisms,
264  Range &simple_master_tris,
265  Range &simple_slave_tris) {
266  MoFEMFunctionBegin;
267 
268  EntityHandle meshset_prisms;
269  CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
270  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
271  bit_levels.back(), BitRefLevel().set(), MBPRISM, meshset_prisms);
272  CHKERR moab.get_entities_by_handle(meshset_prisms, simple_contact_prisms);
273  CHKERR moab.delete_entities(&meshset_prisms, 1);
274 
275  EntityHandle tri;
276  for (Range::iterator pit = simple_contact_prisms.begin();
277  pit != simple_contact_prisms.end(); pit++) {
278  CHKERR moab.side_element(*pit, 2, 3, tri);
279  simple_master_tris.insert(tri);
280  CHKERR moab.side_element(*pit, 2, 4, tri);
281  simple_slave_tris.insert(tri);
282  }
283 
284  MoFEMFunctionReturn(0);
285  };
286 
287  Range simple_contact_prisms, simple_master_tris, simple_slave_tris;
288 
289  if (!ignore_contact) {
290  CHKERR add_prism_interface(bit_levels);
291  CHKERR find_contact_prisms(bit_levels, simple_contact_prisms,
292  simple_master_tris, simple_slave_tris);
293  }
294 
295  Range mortar_contact_prisms, mortar_master_tris, mortar_slave_tris;
296 
297  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
298  if (it->getName().compare(0, 13, "MORTAR_MASTER") == 0) {
299  CHKERR m_field.get_moab().get_entities_by_type(
300  it->meshset, MBTRI, mortar_master_tris, true);
301  }
302  }
303 
304  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
305  if (it->getName().compare(0, 12, "MORTAR_SLAVE") == 0) {
306  CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI,
307  mortar_slave_tris, true);
308  }
309  }
310 
311  ContactSearchKdTree contact_search_kd_tree(m_field);
312 
313  boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
314  contact_commondata_multi_index;
315 
316  contact_commondata_multi_index =
317  boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
318  new ContactSearchKdTree::ContactCommonData_multiIndex());
319 
320  CHKERR contact_search_kd_tree.buildTree(mortar_master_tris);
321 
322  CHKERR contact_search_kd_tree.contactSearchAlgorithm(
323  mortar_master_tris, mortar_slave_tris, contact_commondata_multi_index,
324  mortar_contact_prisms);
325 
326  EntityHandle meshset_slave_master_prisms;
327  CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
328 
329  CHKERR
330  moab.add_entities(meshset_slave_master_prisms, mortar_contact_prisms);
331 
332  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
333  meshset_slave_master_prisms, 3, bit_levels.back());
334 
335  Range tris_from_prism;
336  // find actual masters and slave used
337  CHKERR m_field.get_moab().get_adjacencies(mortar_contact_prisms, 2, true,
338  tris_from_prism,
339  moab::Interface::UNION);
340 
341  tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
342  mortar_slave_tris = intersect(tris_from_prism, mortar_slave_tris);
343 
344  Range all_contact_prisms, all_slave_tris;
345  all_contact_prisms = unite(simple_contact_prisms, mortar_contact_prisms);
346  all_slave_tris = unite(simple_slave_tris, mortar_slave_tris);
347 
348  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
349  MB_TAG_SPARSE, MF_ZERO);
350 
351  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
352  3, MB_TAG_SPARSE, MF_ZERO);
353 
354  if (!eigen_pos_flag) {
355  CHKERR m_field.add_field("EIGEN_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
356  3, MB_TAG_SPARSE, MF_ZERO);
357  }
358 
359  CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
360  MB_TAG_SPARSE, MF_ZERO);
361 
362  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
363  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
364  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
365  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
366  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
367 
368  // Declare problem add entities (by tets) to the field
369  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
370  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
371  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
372  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
373  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
374 
375  if (!eigen_pos_flag) {
376  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
377  CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
378  CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
379  CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
380  CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_POSITIONS", 1);
381  }
382 
383  CHKERR m_field.add_ents_to_field_by_type(all_slave_tris, MBTRI, "LAGMULT");
384 
385  CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
386  CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
387  CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
388 
389  auto set_contact_order = [&](Range &contact_prisms, int order_contact,
390  int nb_ho_levels) {
391  MoFEMFunctionBegin;
392  Range contact_tris, contact_edges;
393  CHKERR moab.get_adjacencies(contact_prisms, 2, false, contact_tris,
394  moab::Interface::UNION);
395  contact_tris = contact_tris.subset_by_type(MBTRI);
396  CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
397  moab::Interface::UNION);
398  Range ho_ents;
399  ho_ents.merge(contact_tris);
400  ho_ents.merge(contact_edges);
401  for (int ll = 0; ll < nb_ho_levels; ll++) {
402  Range ents, verts, tets;
403  CHKERR moab.get_connectivity(ho_ents, verts, true);
404  CHKERR moab.get_adjacencies(verts, 3, false, tets,
405  moab::Interface::UNION);
406  tets = tets.subset_by_type(MBTET);
407  for (auto d : {1, 2}) {
408  CHKERR moab.get_adjacencies(tets, d, false, ents,
409  moab::Interface::UNION);
410  }
411  ho_ents = unite(ho_ents, ents);
412  ho_ents = unite(ho_ents, tets);
413  }
414 
415  CHKERR m_field.set_field_order(ho_ents, "SPATIAL_POSITION",
416  order_contact);
417 
418  MoFEMFunctionReturn(0);
419  };
420 
421  if (!ignore_contact && order_contact > order) {
422  CHKERR set_contact_order(all_contact_prisms, order_contact, nb_ho_levels);
423  }
424 
425  // build field
426  CHKERR m_field.build_fields();
427 
428  // Projection on "x" field
429  {
430  Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
431  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
432  }
433  // Projection on "X" field
434  {
435  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
436  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
437  }
438 
439  Range slave_verts;
440  CHKERR moab.get_connectivity(all_slave_tris, slave_verts, true);
441  CHKERR m_field.getInterface<FieldBlas>()->setField(0.0, MBVERTEX,
442  slave_verts, "LAGMULT");
443 
444  // Add elastic element
445  boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
446  boost::make_shared<std::map<int, BlockData>>();
447  CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
448 
449  boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
450  new VolumeElementForcesAndSourcesCore(m_field));
451  boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
452  new VolumeElementForcesAndSourcesCore(m_field));
453  fe_elastic_lhs_ptr->getRuleHook = VolRule();
454  fe_elastic_rhs_ptr->getRuleHook = VolRule();
455 
456  CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
457  "SPATIAL_POSITION",
458  "MESH_NODE_POSITIONS", false);
459 
460  auto data_hooke_element_at_pts =
461  boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
462  CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
463  block_sets_ptr, "SPATIAL_POSITION",
464  "MESH_NODE_POSITIONS", false, false,
465  MBTET, data_hooke_element_at_pts);
466  auto thermal_strain =
467  [&thermal_expansion_coef, &init_temp, &test_num](
468  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords) {
469  FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
470  FTensor::Index<'i', 3> i;
471  FTensor::Index<'j', 3> j;
472 
473  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
474  double temp;
475 
476  t_thermal_strain(i, j) = 0.0;
477 
478  switch (test_num) {
479  case 0:
480  // Put here analytical formula which may depend on coordinates
481  temp = init_temp - 1.0;
482  t_thermal_strain(i, j) =
483  thermal_expansion_coef * (temp - init_temp) * t_kd(i, j);
484  break;
485  case 1:
486  case 2:
487  t_thermal_strain(i, j) = -thermal_expansion_coef * t_kd(i, j);
488  break;
489  case 3:
490  t_thermal_strain(2, 2) = thermal_expansion_coef;
491  break;
492  case 4:
493  t_thermal_strain(i, j) = thermal_expansion_coef * t_kd(i, j);
494  break;
495  default:
496  break;
497  }
498 
499  return t_thermal_strain;
500  };
501 
502  if (analytical_input) {
503  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
504  new HookeElement::OpAnalyticalInternalStrain_dx<0>(
505  "SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
506  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
507  new HookeInternalStressElement::OpGetAnalyticalInternalStress<0>(
508  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
509  thermal_strain));
510  } else {
511  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
512  new HookeInternalStressElement::OpGetInternalStress(
513  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
514  moab, stress_tag_name));
515  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
516  new HookeInternalStressElement::OpInternalStrain_dx(
517  "SPATIAL_POSITION", data_hooke_element_at_pts));
518  }
519 
520  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
521  new OpCalculateVectorFieldGradient<3, 3>(
522  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
523  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
524  new OpCalculateVectorFieldGradient<3, 3>(
525  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
526  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
527  new HookeInternalStressElement::OpSaveStress(
528  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
529  *block_sets_ptr.get(), moab, scale_factor, save_mean_stress, false,
530  false));
531 
532  Range all_tets;
533  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
534  if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
535  Range tets;
536  CHKERR moab.get_entities_by_dimension(bit->getMeshset(), 3, tets, true);
537  all_tets.merge(tets);
538  }
539  }
540  Skinner skinner(&moab);
541  Range skin_tris;
542  CHKERR skinner.find_skin(0, all_tets, false, skin_tris);
543 
544  CHKERR m_field.add_finite_element("SKIN", MF_ZERO);
545  CHKERR m_field.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
546  CHKERR m_field.modify_finite_element_add_field_row("SKIN",
547  "SPATIAL_POSITION");
548  CHKERR m_field.modify_finite_element_add_field_col("SKIN",
549  "SPATIAL_POSITION");
550  CHKERR m_field.modify_finite_element_add_field_data("SKIN",
551  "SPATIAL_POSITION");
552  CHKERR m_field.modify_finite_element_add_field_data("SKIN",
553  "MESH_NODE_POSITIONS");
554  CHKERR m_field.modify_finite_element_add_field_data("SKIN",
555  "EIGEN_POSITIONS");
556 
557  auto make_mortar_contact_element = [&]() {
558  return boost::make_shared<MortarContactProblem::MortarContactElement>(
559  m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
560  "MESH_NODE_POSITIONS");
561  };
562 
563  auto make_convective_mortar_master_element = [&]() {
564  return boost::make_shared<
565  MortarContactProblem::MortarConvectMasterContactElement>(
566  m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
567  "MESH_NODE_POSITIONS");
568  };
569 
570  auto make_convective_mortar_slave_element = [&]() {
571  return boost::make_shared<
572  MortarContactProblem::MortarConvectSlaveContactElement>(
573  m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
574  "MESH_NODE_POSITIONS");
575  };
576 
577  auto make_mortar_contact_common_data = [&]() {
578  return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
579  m_field);
580  };
581 
582  auto get_mortar_contact_rhs = [&](auto contact_problem, auto make_element,
583  bool is_alm = false) {
584  auto fe_rhs_mortar_contact = make_element();
585  auto common_data_mortar_contact = make_mortar_contact_common_data();
586  if (print_contact_state) {
587  fe_rhs_mortar_contact->contactStateVec =
588  common_data_mortar_contact->gaussPtsStateVec;
589  }
590  contact_problem->setContactOperatorsRhs(
591  fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
592  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
593  return fe_rhs_mortar_contact;
594  };
595 
596  auto get_mortar_master_traction_rhs = [&](auto contact_problem,
597  auto make_element,
598  bool is_alm = false) {
599  auto fe_rhs_mortar_contact = make_element();
600  auto common_data_mortar_contact = make_mortar_contact_common_data();
601  contact_problem->setMasterForceOperatorsRhs(
602  fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
603  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
604  return fe_rhs_mortar_contact;
605  };
606 
607  auto get_mortar_master_traction_lhs = [&](auto contact_problem,
608  auto make_element,
609  bool is_alm = false) {
610  auto fe_lhs_mortar_contact = make_element();
611  auto common_data_mortar_contact = make_mortar_contact_common_data();
612  contact_problem->setMasterForceOperatorsLhs(
613  fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
614  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
615  return fe_lhs_mortar_contact;
616  };
617 
618  auto get_mortar_master_help_traction_lhs =
619  [&](auto contact_problem, auto make_element, bool is_alm = false) {
620  auto fe_lhs_mortar_contact = make_element();
621  auto common_data_mortar_contact = make_mortar_contact_common_data();
622  contact_problem->setMasterForceOperatorsLhs(
623  fe_lhs_mortar_contact, common_data_mortar_contact,
624  "SPATIAL_POSITION", "LAGMULT", is_alm);
625  return fe_lhs_mortar_contact;
626  };
627 
628  auto get_mortar_contact_lhs = [&](auto contact_problem, auto make_element,
629  bool is_alm = false) {
630  auto fe_lhs_mortar_contact = make_element();
631  auto common_data_mortar_contact = make_mortar_contact_common_data();
632  contact_problem->setContactOperatorsLhs(
633  fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
634  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
635  return fe_lhs_mortar_contact;
636  };
637 
638  auto get_mortar_contact_help_lhs =
639  [&](auto contact_problem, auto make_element, bool is_alm = false) {
640  auto fe_lhs_mortar_contact = make_element();
641  auto common_data_mortar_contact = make_mortar_contact_common_data();
642  contact_problem->setContactOperatorsLhs(
643  fe_lhs_mortar_contact, common_data_mortar_contact,
644  "SPATIAL_POSITION", "LAGMULT", is_alm);
645  return fe_lhs_mortar_contact;
646  };
647 
648  auto cn_value_ptr = boost::make_shared<double>(cn_value);
649  auto mortar_contact_problem = boost::make_shared<MortarContactProblem>(
650  m_field, contact_commondata_multi_index, cn_value_ptr, is_newton_cotes);
651 
652  // add fields to the global matrix by adding the element
653 
654  if (!eigen_pos_flag)
655  mortar_contact_problem->addMortarContactElement(
656  "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
657  mortar_contact_prisms);
658  else
659  mortar_contact_problem->addMortarContactElement(
660  "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
661  mortar_contact_prisms, "MESH_NODE_POSITIONS", eigen_pos_flag,
662  "EIGEN_POSITIONS");
663 
664  auto make_simple_contact_element = [&]() {
665  return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
666  m_field);
667  };
668 
669  auto make_convective_simple_master_element = [&]() {
670  return boost::make_shared<
671  SimpleContactProblem::ConvectMasterContactElement>(
672  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
673  };
674 
675  auto make_convective_simple_slave_element = [&]() {
676  return boost::make_shared<
677  SimpleContactProblem::ConvectSlaveContactElement>(
678  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
679  };
680 
681  auto make_simple_contact_common_data = [&]() {
682  return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
683  m_field);
684  };
685 
686  auto get_simple_contact_rhs = [&](auto contact_problem, auto make_element,
687  bool is_alm = false) {
688  auto fe_rhs_simple_contact = make_element();
689  auto common_data_simple_contact = make_simple_contact_common_data();
690  if (print_contact_state) {
691  fe_rhs_simple_contact->contactStateVec =
692  common_data_simple_contact->gaussPtsStateVec;
693  }
694  contact_problem->setContactOperatorsRhs(
695  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
696  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
697  return fe_rhs_simple_contact;
698  };
699 
700  auto get_simple_master_traction_rhs = [&](auto contact_problem,
701  auto make_element,
702  bool is_alm = false) {
703  auto fe_rhs_simple_contact = make_element();
704  auto common_data_simple_contact = make_simple_contact_common_data();
705  contact_problem->setMasterForceOperatorsRhs(
706  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
707  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
708  return fe_rhs_simple_contact;
709  };
710 
711  auto get_simple_master_traction_lhs =
712  [&](auto contact_problem, auto make_element, bool is_alm = false) {
713  auto fe_lhs_simple_contact = make_element();
714  auto common_data_simple_contact = make_simple_contact_common_data();
715  contact_problem->setMasterForceOperatorsLhs(
716  fe_lhs_simple_contact, common_data_simple_contact,
717  "SPATIAL_POSITION", "LAGMULT", is_alm);
718  return fe_lhs_simple_contact;
719  };
720 
721  auto get_contact_lhs = [&](auto contact_problem, auto make_element,
722  bool is_alm = false) {
723  auto fe_lhs_simple_contact = make_element();
724  auto common_data_simple_contact = make_simple_contact_common_data();
725  contact_problem->setContactOperatorsLhs(
726  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
727  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
728  return fe_lhs_simple_contact;
729  };
730 
731  auto simple_contact_problem = boost::make_shared<SimpleContactProblem>(
732  m_field, cn_value_ptr, is_newton_cotes);
733  simple_contact_problem->addContactElement("SIMPLE_CONTACT_ELEM",
734  "SPATIAL_POSITION", "LAGMULT",
735  simple_contact_prisms);
736 
737  simple_contact_problem->addPostProcContactElement(
738  "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
739  "MESH_NODE_POSITIONS", all_slave_tris);
740 
741  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "SPATIAL_POSITION");
742 
743  // Add spring boundary condition applied on surfaces.
744  CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
745  "MESH_NODE_POSITIONS");
746 
747  // build finite elemnts
748  CHKERR m_field.build_finite_elements();
749 
750  // build adjacencies
751  CHKERR m_field.build_adjacencies(bit_levels.back());
752 
753  // define problems
754  CHKERR m_field.add_problem("CONTACT_PROB", MF_ZERO);
755 
756  // set refinement level for problem
757  CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
758  bit_levels.back());
759 
760  DMType dm_name = "DMMOFEM";
761  CHKERR DMRegister_MoFEM(dm_name);
762 
763  SmartPetscObj<DM> dm;
764  dm = createSmartDM(m_field.get_comm(), dm_name);
765 
766  // create dm instance
767  CHKERR DMSetType(dm, dm_name);
768 
769  // set dm datastruture which created mofem datastructures
770  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_levels.back());
771  CHKERR DMSetFromOptions(dm);
772  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
773  // add elements to dm
774  CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT_ELEM");
775  CHKERR DMMoFEMAddElement(dm, "SIMPLE_CONTACT_ELEM");
776  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
777  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
778  CHKERR DMMoFEMAddElement(dm, "SPRING");
779  CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
780  CHKERR DMMoFEMAddElement(dm, "SKIN");
781 
782  CHKERR DMSetUp(dm);
783 
784  // Vector of DOFs and the RHS
785  auto D = smartCreateDMVector(dm);
786  auto F = smartVectorDuplicate(D);
787 
788  // Stiffness matrix
789  auto Aij = smartCreateDMMatrix(dm);
790 
791  CHKERR VecZeroEntries(D);
792  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
793  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
794  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
795 
796  CHKERR VecZeroEntries(F);
797  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
798  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
799 
800  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
801  CHKERR MatZeroEntries(Aij);
802 
803  fe_elastic_rhs_ptr->snes_f = F;
804  fe_elastic_lhs_ptr->snes_B = Aij;
805 
806  // Dirichlet BC
807  boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
808  boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
809  m_field, "SPATIAL_POSITION", Aij, D, F));
810 
811  dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
812  dirichlet_bc_ptr->snes_x = D;
813 
814  // Assemble pressure and traction forces
815  boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
816  CHKERR MetaNeumannForces::setMomentumFluxOperators(
817  m_field, neumann_forces, NULL, "SPATIAL_POSITION");
818 
819  boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
820  neumann_forces.begin();
821  for (; mit != neumann_forces.end(); mit++) {
822  mit->second->methodsOp.push_back(new MortarContactProblem::LoadScale());
823  CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
824  &mit->second->getLoopFe(), NULL, NULL);
825  }
826 
827  // Implementation of spring element
828  // Create new instances of face elements for springs
829  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
830  new FaceElementForcesAndSourcesCore(m_field));
831  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
832  new FaceElementForcesAndSourcesCore(m_field));
833 
834  CHKERR MetaSpringBC::setSpringOperators(
835  m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
836  "MESH_NODE_POSITIONS");
837 
838  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
839  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
840  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
841  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
842  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL,
843  dirichlet_bc_ptr.get(), NULL);
844  if (convect_pts == PETSC_TRUE) {
845  CHKERR DMMoFEMSNESSetFunction(
846  dm, "MORTAR_CONTACT_ELEM",
847  get_mortar_contact_rhs(mortar_contact_problem,
848  make_convective_mortar_master_element),
849  PETSC_NULL, PETSC_NULL);
850  CHKERR DMMoFEMSNESSetFunction(
851  dm, "MORTAR_CONTACT_ELEM",
852  get_mortar_master_traction_rhs(mortar_contact_problem,
853  make_convective_mortar_slave_element),
854  PETSC_NULL, PETSC_NULL);
855  } else {
856  CHKERR DMMoFEMSNESSetFunction(
857  dm, "MORTAR_CONTACT_ELEM",
858  get_mortar_contact_rhs(mortar_contact_problem,
859  make_mortar_contact_element, alm_flag),
860  PETSC_NULL, PETSC_NULL);
861  CHKERR DMMoFEMSNESSetFunction(
862  dm, "MORTAR_CONTACT_ELEM",
863  get_mortar_master_traction_rhs(mortar_contact_problem,
864  make_mortar_contact_element, alm_flag),
865  PETSC_NULL, PETSC_NULL);
866  }
867 
868  if (convect_pts == PETSC_TRUE) {
869  CHKERR DMMoFEMSNESSetFunction(
870  dm, "SIMPLE_CONTACT_ELEM",
871  get_simple_contact_rhs(simple_contact_problem,
872  make_convective_simple_master_element),
873  PETSC_NULL, PETSC_NULL);
874  CHKERR DMMoFEMSNESSetFunction(
875  dm, "SIMPLE_CONTACT_ELEM",
876  get_simple_master_traction_rhs(simple_contact_problem,
877  make_convective_simple_slave_element),
878  PETSC_NULL, PETSC_NULL);
879  } else {
880  CHKERR DMMoFEMSNESSetFunction(
881  dm, "SIMPLE_CONTACT_ELEM",
882  get_simple_contact_rhs(simple_contact_problem,
883  make_simple_contact_element, alm_flag),
884  PETSC_NULL, PETSC_NULL);
885  CHKERR DMMoFEMSNESSetFunction(
886  dm, "SIMPLE_CONTACT_ELEM",
887  get_simple_master_traction_rhs(simple_contact_problem,
888  make_simple_contact_element, alm_flag),
889  PETSC_NULL, PETSC_NULL);
890  }
891 
892  CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", fe_elastic_rhs_ptr, PETSC_NULL,
893  PETSC_NULL);
894 
895  CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
896  PETSC_NULL);
897  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, PETSC_NULL, PETSC_NULL,
898  dirichlet_bc_ptr.get());
899 
900  boost::shared_ptr<FEMethod> fe_null;
901  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
902  fe_null);
903  if (convect_pts == PETSC_TRUE) {
904  CHKERR DMMoFEMSNESSetJacobian(
905  dm, "MORTAR_CONTACT_ELEM",
906  get_mortar_contact_help_lhs(mortar_contact_problem,
907  make_convective_mortar_master_element),
908  PETSC_NULL, PETSC_NULL);
909  CHKERR DMMoFEMSNESSetJacobian(
910  dm, "MORTAR_CONTACT_ELEM",
911  get_mortar_master_help_traction_lhs(
912  mortar_contact_problem, make_convective_mortar_slave_element),
913  PETSC_NULL, PETSC_NULL);
914  } else {
915  CHKERR DMMoFEMSNESSetJacobian(
916  dm, "MORTAR_CONTACT_ELEM",
917  get_mortar_contact_lhs(mortar_contact_problem,
918  make_mortar_contact_element, alm_flag),
919  PETSC_NULL, PETSC_NULL);
920  CHKERR DMMoFEMSNESSetJacobian(
921  dm, "MORTAR_CONTACT_ELEM",
922  get_mortar_master_traction_lhs(mortar_contact_problem,
923  make_mortar_contact_element, alm_flag),
924  PETSC_NULL, PETSC_NULL);
925  }
926 
927  if (convect_pts == PETSC_TRUE) {
928  CHKERR DMMoFEMSNESSetJacobian(
929  dm, "SIMPLE_CONTACT_ELEM",
930  get_contact_lhs(simple_contact_problem,
931  make_convective_simple_master_element),
932  PETSC_NULL, PETSC_NULL);
933  CHKERR DMMoFEMSNESSetJacobian(
934  dm, "SIMPLE_CONTACT_ELEM",
935  get_simple_master_traction_lhs(simple_contact_problem,
936  make_convective_simple_slave_element),
937  PETSC_NULL, PETSC_NULL);
938  } else {
939  CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ELEM",
940  get_contact_lhs(simple_contact_problem,
941  make_simple_contact_element,
942  alm_flag),
943  PETSC_NULL, PETSC_NULL);
944  CHKERR DMMoFEMSNESSetJacobian(
945  dm, "SIMPLE_CONTACT_ELEM",
946  get_simple_master_traction_lhs(simple_contact_problem,
947  make_simple_contact_element, alm_flag),
948  PETSC_NULL, PETSC_NULL);
949  }
950  CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", fe_elastic_lhs_ptr, PETSC_NULL,
951  PETSC_NULL);
952  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
953  PETSC_NULL);
954  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, fe_null,
955  dirichlet_bc_ptr);
956 
957  if (test_num) {
958  char testing_options[] = "-ksp_type fgmres "
959  "-pc_type lu "
960  "-pc_factor_mat_solver_type mumps "
961  "-snes_type newtonls "
962  "-snes_linesearch_type basic "
963  "-snes_max_it 10 "
964  "-snes_atol 1e-8 "
965  "-snes_rtol 1e-8 ";
966  CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
967  }
968 
969  auto snes = MoFEM::createSNES(m_field.get_comm());
970  CHKERR SNESSetDM(snes, dm);
971  SnesCtx *snes_ctx;
972  // create snes nonlinear solver
973  {
974  CHKERR SNESSetDM(snes, dm);
975  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
976  CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
977  CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
978  CHKERR SNESSetFromOptions(snes);
979  }
980 
981  /// Post proc on the skin
982  PostProcFaceOnRefinedMesh post_proc_skin(m_field);
983  CHKERR post_proc_skin.generateReferenceElementMesh();
984  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false, false);
985  CHKERR post_proc_skin.addFieldValuesPostProc("SPATIAL_POSITION");
986  CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
987  CHKERR post_proc_skin.addFieldValuesPostProc("EIGEN_POSITIONS");
988 
989  struct OpGetFieldGradientValuesOnSkin
990  : public FaceElementForcesAndSourcesCore::UserDataOperator {
991 
992  const std::string feVolName;
993  boost::shared_ptr<VolSideFe> sideOpFe;
994 
995  OpGetFieldGradientValuesOnSkin(const std::string field_name,
996  const std::string vol_fe_name,
997  boost::shared_ptr<VolSideFe> side_fe)
998  : FaceElementForcesAndSourcesCore::UserDataOperator(
999  field_name, UserDataOperator::OPROW),
1000  feVolName(vol_fe_name), sideOpFe(side_fe) {}
1001 
1002  MoFEMErrorCode doWork(int side, EntityType type,
1003  EntitiesFieldData::EntData &data) {
1004  MoFEMFunctionBegin;
1005  if (type != MBVERTEX)
1006  MoFEMFunctionReturnHot(0);
1007  CHKERR loopSideVolumes(feVolName, *sideOpFe);
1008  MoFEMFunctionReturn(0);
1009  }
1010  };
1011 
1012  boost::shared_ptr<VolSideFe> vol_side_fe_ptr =
1013  boost::make_shared<VolSideFe>(m_field);
1014  vol_side_fe_ptr->getOpPtrVector().push_back(
1015  new OpCalculateVectorFieldGradient<3, 3>(
1016  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
1017  vol_side_fe_ptr->getOpPtrVector().push_back(
1018  new OpCalculateVectorFieldGradient<3, 3>(
1019  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
1020 
1021  post_proc_skin.getOpPtrVector().push_back(
1022  new OpGetFieldGradientValuesOnSkin("SPATIAL_POSITION", "ELASTIC",
1023  vol_side_fe_ptr));
1024  post_proc_skin.getOpPtrVector().push_back(
1025  new OpCalculateVectorFieldValues<3>(
1026  "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
1027  post_proc_skin.getOpPtrVector().push_back(
1028  new OpCalculateVectorFieldValues<3>(
1029  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
1030 
1031  post_proc_skin.getOpPtrVector().push_back(
1032  new HookeInternalStressElement::OpPostProcHookeElement<
1033  FaceElementForcesAndSourcesCore>(
1034  "SPATIAL_POSITION", data_hooke_element_at_pts,
1035  *block_sets_ptr.get(), post_proc_skin.postProcMesh,
1036  post_proc_skin.mapGaussPts, false, false));
1037 
1038  for (int ss = 0; ss != nb_sub_steps; ++ss) {
1039  if (!ignore_pressure) {
1040  MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
1041  } else {
1042  MortarContactProblem::LoadScale::lAmbda = 0;
1043  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Ignoring pressure...\n");
1044  }
1045  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
1046  MortarContactProblem::LoadScale::lAmbda);
1047 
1048  CHKERR SNESSolve(snes, PETSC_NULL, D);
1049 
1050  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1051  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1052  }
1053 
1054  // save on mesh
1055  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1056 
1057  SmartPetscObj<Vec> v_energy;
1058  CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "SPATIAL_POSITION",
1059  "MESH_NODE_POSITIONS", false, false,
1060  v_energy);
1061  double energy;
1062  CHKERR VecSum(v_energy, &energy);
1063  // Print elastic energy
1064  MOFEM_LOG_C("WORLD", Sev::inform, "Elastic energy: %8.8e", energy);
1065 
1066  // moab_instance
1067  moab::Core mb_post; // create database
1068  moab::Interface &moab_proc = mb_post; // create interface to database
1069 
1070  auto common_data_mortar_contact = make_mortar_contact_common_data();
1071 
1072  boost::shared_ptr<MortarContactProblem::MortarContactElement>
1073  fe_post_proc_mortar_contact;
1074  if (convect_pts == PETSC_TRUE) {
1075  fe_post_proc_mortar_contact = make_convective_mortar_master_element();
1076  } else {
1077  fe_post_proc_mortar_contact = make_mortar_contact_element();
1078  }
1079 
1080  std::array<double, 2> nb_gauss_pts;
1081  std::array<double, 2> contact_area;
1082 
1083  auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
1084  MoFEMFunctionBegin;
1085  CHKERR VecAssemblyBegin(vec);
1086  CHKERR VecAssemblyEnd(vec);
1087  const double *array;
1088  CHKERR VecGetArrayRead(vec, &array);
1089  if (m_field.get_comm_rank() == 0) {
1090  for (int i : {0, 1})
1091  data[i] = array[i];
1092  }
1093  CHKERR VecRestoreArrayRead(vec, &array);
1094  MoFEMFunctionReturn(0);
1095  };
1096 
1097  mb_post.delete_mesh();
1098 
1099  if (!mortar_contact_prisms.empty()) {
1100  mortar_contact_problem->setContactOperatorsForPostProc(
1101  fe_post_proc_mortar_contact, common_data_mortar_contact, m_field,
1102  "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1103  "EIGEN_POSITIONS");
1104 
1105  CHKERR VecZeroEntries(common_data_mortar_contact->gaussPtsStateVec);
1106  CHKERR VecZeroEntries(common_data_mortar_contact->contactAreaVec);
1107 
1108  CHKERR DMoFEMLoopFiniteElements(dm, "MORTAR_CONTACT_ELEM",
1109  fe_post_proc_mortar_contact);
1110 
1111  CHKERR get_contact_data(common_data_mortar_contact->gaussPtsStateVec,
1112  nb_gauss_pts);
1113  CHKERR get_contact_data(common_data_mortar_contact->contactAreaVec,
1114  contact_area);
1115 
1116  if (m_field.get_comm_rank() == 0) {
1117  PetscPrintf(PETSC_COMM_SELF,
1118  "Active gauss pts (mortar): %d out of %d\n",
1119  (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
1120 
1121  PetscPrintf(
1122  PETSC_COMM_SELF,
1123  "Active contact area (mortar): %8.8f out of %8.8f (%8.8f% %)\n",
1124  contact_area[0], contact_area[1],
1125  contact_area[0] / contact_area[1] * 100.);
1126  }
1127  }
1128 
1129  auto common_data_simple_contact = make_simple_contact_common_data();
1130 
1131  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1132  fe_post_proc_simple_contact;
1133  if (convect_pts == PETSC_TRUE) {
1134  fe_post_proc_simple_contact = make_convective_simple_master_element();
1135  } else {
1136  fe_post_proc_simple_contact = make_simple_contact_element();
1137  }
1138 
1139  if (!simple_contact_prisms.empty()) {
1140  simple_contact_problem->setContactOperatorsForPostProc(
1141  fe_post_proc_simple_contact, common_data_simple_contact, m_field,
1142  "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1143  "EIGEN_POSITIONS");
1144 
1145  CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
1146  CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
1147 
1148  CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_CONTACT_ELEM",
1149  fe_post_proc_simple_contact);
1150 
1151  CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
1152  nb_gauss_pts);
1153  CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
1154  contact_area);
1155 
1156  if (m_field.get_comm_rank() == 0) {
1157  MOFEM_LOG_C("SELF", Sev::inform,
1158  "Active gauss pts (matching): %d out of %d",
1159  (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
1160  MOFEM_LOG_C(
1161  "SELF", Sev::inform,
1162  "Active contact area (matching): %8.8f out of %8.8f (%8.8f% %)",
1163  contact_area[0], contact_area[1],
1164  contact_area[0] / contact_area[1] * 100.);
1165  }
1166  }
1167 
1168  if (!ignore_contact) {
1169  string out_file_name;
1170  std::ostringstream strm;
1171  strm << "out_contact_integ_pts"
1172  << ".h5m";
1173  out_file_name = strm.str();
1174  MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s\n",
1175  out_file_name.c_str());
1176  CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
1177  "PARALLEL=WRITE_PART");
1178  }
1179 
1180  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
1181  new PostProcFaceOnRefinedMesh(m_field));
1182 
1183  CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
1184  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *post_proc_contact_ptr, false,
1185  false);
1186  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("LAGMULT");
1187  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("SPATIAL_POSITION");
1188  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
1189 
1190  if (!all_slave_tris.empty()) {
1191  CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
1192  post_proc_contact_ptr);
1193  string out_file_name;
1194  std::ostringstream stm;
1195  stm << "out_contact_surface"
1196  << ".h5m";
1197  out_file_name = stm.str();
1198  MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
1199  out_file_name.c_str());
1200  CHKERR post_proc_contact_ptr->postProcMesh.write_file(
1201  out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
1202  }
1203 
1204  CHKERR m_field.getInterface<FieldBlas>()->fieldCopy(1., "SPATIAL_POSITION",
1205  "EIGEN_POSITIONS");
1206 
1207  {
1208  PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
1209  CHKERR DMoFEMLoopFiniteElements(dm, "SKIN", &post_proc_skin);
1210  ostringstream stm;
1211  string out_file_name;
1212  stm << "out_skin.h5m";
1213  out_file_name = stm.str();
1214  MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
1215  out_file_name.c_str());
1216  CHKERR post_proc_skin.writeFile(stm.str());
1217  }
1218 
1219  const int n_parts = m_field.get_comm_size();
1220  if (m_field.get_comm_rank() == 0) {
1221  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
1222  dm, "ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
1223 
1224  Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
1225 
1226  CHKERR moab.get_adjacencies(all_contact_prisms, 2, true, faces,
1227  moab::Interface::UNION);
1228  tris = faces.subset_by_type(MBTRI);
1229  quads = faces.subset_by_type(MBQUAD);
1230  CHKERR moab.get_adjacencies(tris, 1, true, tris_edges,
1231  moab::Interface::UNION);
1232  CHKERR moab.get_adjacencies(quads, 1, true, quads_edges,
1233  moab::Interface::UNION);
1234 
1235  ents_to_delete.merge(all_contact_prisms);
1236  ents_to_delete.merge(quads);
1237  ents_to_delete.merge(subtract(quads_edges, tris_edges));
1238 
1239  CHKERR moab.delete_entities(ents_to_delete);
1240 
1241  MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s", output_mesh_name);
1242  CHKERR moab.write_file(output_mesh_name, "MOAB");
1243 
1244  auto get_tag_handle = [&](auto name, auto size) {
1245  Tag th;
1246  std::vector<double> def_vals(size, 0.0);
1247  CHKERR moab.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1248  MB_TAG_CREAT | MB_TAG_SPARSE,
1249  def_vals.data());
1250  return th;
1251  };
1252 
1253  if (test_num) {
1254  Range tets;
1255  CHKERR moab.get_entities_by_dimension(0, 3, tets);
1256  EntityHandle tet = tets.front();
1257  std::array<double, 9> internal_stress, actual_stress;
1258  std::array<double, 9> internal_stress_ref, actual_stress_ref;
1259  std::array<double, 2> nb_gauss_pts_ref, contact_area_ref;
1260  switch (test_num) {
1261  case 3:
1262  actual_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1263  if (strcmp(stress_tag_name, "INTERNAL_STRESS") == 0)
1264  internal_stress_ref = {0., 0., -200., 0., 0., 0., 0., 0., 0.};
1265  else
1266  internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1267  break;
1268  case 4:
1269  nb_gauss_pts_ref = {216, 492};
1270  contact_area_ref = {0.125, 0.25};
1271  break;
1272  default:
1273  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Test number %d not found",
1274  test_num);
1275  }
1276 
1277  auto th_internal_stress = get_tag_handle("MED_INTERNAL_STRESS", 9);
1278  auto th_actual_stress = get_tag_handle("MED_ACTUAL_STRESS", 9);
1279  CHKERR moab.tag_get_data(th_internal_stress, &tet, 1,
1280  internal_stress.data());
1281  CHKERR moab.tag_get_data(th_actual_stress, &tet, 1,
1282  actual_stress.data());
1283  if (test_num == 4) {
1284  const double eps = 1e-3;
1285  if (std::abs(nb_gauss_pts_ref[0] - nb_gauss_pts[0]) > eps) {
1286  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1287  "Wrong number of active contact gauss pts: should be %d "
1288  "but is %d",
1289  (int)nb_gauss_pts_ref[0], (int)nb_gauss_pts[0]);
1290  }
1291  if (std::abs(nb_gauss_pts_ref[1] - nb_gauss_pts[1]) > eps) {
1292  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1293  "Wrong total number of contact gauss pts: should be %d "
1294  "but is %d",
1295  (int)nb_gauss_pts_ref[1], (int)nb_gauss_pts[1]);
1296  }
1297  if (std::abs(contact_area_ref[0] - contact_area[0]) > eps) {
1298  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1299  "Wrong active contact area: should be %g "
1300  "but is %g",
1301  contact_area_ref[0], contact_area[0]);
1302  }
1303  if (std::abs(contact_area_ref[1] - contact_area[1]) > eps) {
1304  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1305  "Wrong potential contact area: should be %g "
1306  "but is %g",
1307  contact_area_ref[1], contact_area[1]);
1308  }
1309  } else {
1310  if (save_mean_stress) {
1311  const double eps = 1e-10;
1312  for (int i = 0; i < 9; i++) {
1313  if (std::abs(internal_stress[i] - internal_stress_ref[i]) > eps) {
1314  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1315  "Wrong component %d of internal stress: should be %g "
1316  "but is %g",
1317  i, internal_stress_ref[i], internal_stress[i]);
1318  }
1319  if (std::abs(actual_stress[i] - actual_stress_ref[i]) > eps) {
1320  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1321  "Wrong component %d of actual stress: should be %g "
1322  "but is %g",
1323  i, actual_stress_ref[i], actual_stress[i]);
1324  }
1325  }
1326  }
1327  }
1328  }
1329  }
1330  }
1331  CATCH_ERRORS;
1332 
1333  // finish work cleaning memory, getting statistics, etc
1334  MoFEM::Core::Finalize();
1335 
1336  return 0;
1337 }
HookeInternalStressElement::OpGetInternalStress
Definition: HookeInternalStressElement.hpp:54
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:699
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MortarContactProblem::MortarConvectSlaveContactElement
Element used to integrate on master surfaces. It convects integration points on slaves,...
Definition: MortarContactProblem.hpp:311
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
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
ContactSearchKdTree::ContactCommonData_multiIndex
multi_index_container< boost::shared_ptr< ContactCommonData >, indexed_by< ordered_unique< tag< Prism_tag >, member< ContactCommonData, EntityHandle, &ContactCommonData::pRism > > > > ContactCommonData_multiIndex
Definition: ContactSearchKdTree.hpp:51
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
HookeInternalStressElement::OpSaveStress
Definition: HookeInternalStressElement.hpp:117
BlockData
NonlinearElasticElement::BlockData BlockData
Definition: mortar_contact_thermal.cpp:28
SimpleContactProblem::ConvectSlaveContactElement
Element used to integrate on master surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:206
MoFEM::Version::strVersion
std::string strVersion()
Definition: UnknownInterface.hpp:24
help
static char help[]
Definition: mortar_contact_thermal.cpp:26
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
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
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
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
ContactSearchKdTree
Definition: ContactSearchKdTree.hpp:18
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::UnknownInterface::getFileVersion
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.
Definition: UnknownInterface.cpp:16
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
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
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
MoFEM::DataOperator::doWork
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: DataOperators.hpp:40
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
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.
MortarContactProblem::LoadScale
Definition: MortarContactProblem.hpp:82
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
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
VolSideFe
VolumeElementForcesAndSourcesCoreOnSide VolSideFe
Definition: photon_diffusion.cpp:31
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMoFEMLoopFiniteElementsUpAndLowRank
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:567
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
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
VolRule::operator()
int operator()(int, int, int order) const
Definition: mortar_contact_thermal.cpp:32
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method)
Definition: FaceElementForcesAndSourcesCore.cpp:413
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
main
int main(int argc, char *argv[])
Definition: mortar_contact_thermal.cpp:35
MoFEM::Version
Definition: UnknownInterface.hpp:12
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
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcOnRefMesh.hpp:253
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
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
temp
void temp(int x, int y=10)
Definition: simple.cpp:4
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
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
init_temp
double init_temp
Definition: thermo_elastic.cpp:125
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::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
HookeInternalStressElement::OpPostProcHookeElement
Definition: HookeInternalStressElement.hpp:142
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
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
ContactSearchKdTree::contactSearchAlgorithm
PetscErrorCode contactSearchAlgorithm(Range &range_surf_master, Range &range_surf_slave, boost::shared_ptr< ContactCommonData_multiIndex > contact_commondata_multi_index, Range &range_slave_master_prisms)
Definition: ContactSearchKdTree.hpp:114
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
_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
HookeInternalStressElement::OpGetAnalyticalInternalStress
Definition: HookeInternalStressElement.hpp:88
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
ContactSearchKdTree::buildTree
MoFEMErrorCode buildTree(Range &range_surf_master)
Definition: ContactSearchKdTree.hpp:28
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
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
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
SimpleContactProblem::ConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:181
MoFEM::BlockData
Definition: MeshsetsManager.cpp:755
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
addElasticElement
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
Definition: HookeElement.cpp:533
setOperators
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
HookeInternalStressElement::OpInternalStrain_dx
Definition: HookeInternalStressElement.hpp:74
MortarContactProblem::MortarConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: MortarContactProblem.hpp:282
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
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
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
MoFEM::VolumeElementForcesAndSourcesCoreOnSide
Base volume element used to integrate on skeleton.
Definition: VolumeElementForcesAndSourcesCoreOnSide.hpp:22
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
HookeInternalStressElement.hpp
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
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
PostProcFaceOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: PostProcOnRefMesh.cpp:539
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
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
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