v0.14.0
simple_contact_thermal.cpp
Go to the documentation of this file.
1 /** \file simple_contact_thermal.cpp
2  * \example simple_contact_thermal.cpp
3  *
4  * Implementation of simple contact between surfaces with 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 using namespace MoFEM;
27 using namespace std;
28 
29 
30 static char help[] = "\n";
31 double SimpleContactProblem::LoadScale::lAmbda = 1;
32 using BlockData = NonlinearElasticElement::BlockData;
33 using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide;
34 
35 struct VolRule {
36  int operator()(int, int, int order) const { return 2 * (order - 1); }
37 };
38 
39 int main(int argc, char *argv[]) {
40 
41  const string default_options = "-ksp_type fgmres \n"
42  "-pc_type lu \n"
43  "-pc_factor_mat_solver_type mumps \n"
44  "-mat_mumps_icntl_13 1 \n"
45  "-mat_mumps_icntl_14 800 \n"
46  "-mat_mumps_icntl_20 0 \n"
47  "-mat_mumps_icntl_24 1 \n"
48  "-snes_type newtonls \n"
49  "-snes_linesearch_type basic \n"
50  "-snes_divergence_tolerance 0 \n"
51  "-snes_max_it 50 \n"
52  "-snes_atol 1e-8 \n"
53  "-snes_rtol 1e-10 \n"
54  "-snes_monitor \n"
55  "-ksp_monitor \n"
56  "-snes_converged_reason \n"
57  "-my_order 1 \n"
58  "-my_order_lambda 1 \n"
59  "-my_order_contact 2 \n"
60  "-my_ho_levels_num 1 \n"
61  "-my_step_num 1 \n"
62  "-my_cn_value 1. \n"
63  "-my_r_value 1. \n"
64  "-my_alm_flag 0 \n"
65  "-my_eigen_pos_flag 0 \n";
66 
67  string param_file = "param_file.petsc";
68  if (!static_cast<bool>(ifstream(param_file))) {
69  std::ofstream file(param_file.c_str(), std::ios::ate);
70  if (file.is_open()) {
71  file << default_options;
72  file.close();
73  }
74  }
75 
76  // Initialize MoFEM
77  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
78 
79  // Create mesh database
80  moab::Core mb_instance; // create database
81  moab::Interface &moab = mb_instance; // create interface to database
82 
83  try {
84  PetscBool flg_file;
85  PetscBool flg_file_out;
86 
87  char mesh_file_name[255];
88  char output_mesh_name[255];
89  PetscInt order = 1;
90  PetscInt order_contact = 1;
91  PetscInt nb_ho_levels = 0;
92  PetscInt order_lambda = 1;
93  PetscReal r_value = 1.;
94  PetscReal cn_value = -1;
95  PetscInt nb_sub_steps = 1;
96  PetscBool is_partitioned = PETSC_FALSE;
97  PetscBool is_newton_cotes = PETSC_FALSE;
98  PetscInt test_num = 0;
99  PetscBool convect_pts = PETSC_FALSE;
100  PetscBool print_contact_state = PETSC_FALSE;
101  PetscBool alm_flag = PETSC_FALSE;
102  PetscBool eigen_pos_flag = PETSC_FALSE;
103 
104  PetscReal thermal_expansion_coef = 1e-5;
105  PetscReal init_temp = 250.0;
106  PetscReal scale_factor = 1.0;
107  PetscBool analytical_input = PETSC_TRUE;
108  char stress_tag_name[255];
109  PetscBool flg_tag_name;
110  PetscBool save_mean_stress = PETSC_TRUE;
111  PetscBool ignore_contact = PETSC_FALSE;
112  PetscBool ignore_pressure = PETSC_FALSE;
113 
114  PetscBool deform_flat_flag = PETSC_FALSE;
115  PetscReal flat_shift = 1.0;
116  PetscReal mesh_height = 1.0;
117 
118  PetscBool wave_surf_flag = PETSC_FALSE;
119  PetscInt wave_dim = 2;
120  PetscReal wave_length = 1.0;
121  PetscReal wave_ampl = 0.01;
122 
123  PetscBool delete_prisms = PETSC_FALSE;
124 
125  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
126 
127  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
128  mesh_file_name, 255, &flg_file);
129  CHKERR PetscOptionsString("-my_output_mesh_file", "output mesh file name",
130  "", "mesh.h5m", output_mesh_name, 255,
131  &flg_file_out);
132 
133  CHKERR PetscOptionsInt("-my_order",
134  "approximation order of spatial positions", "", 1,
135  &order, PETSC_NULL);
136  CHKERR PetscOptionsInt(
137  "-my_order_contact",
138  "approximation order of spatial positions in contact interface", "", 1,
139  &order_contact, PETSC_NULL);
140  CHKERR PetscOptionsInt("-my_ho_levels_num", "number of higher order levels",
141  "", 0, &nb_ho_levels, PETSC_NULL);
142  CHKERR PetscOptionsInt("-my_order_lambda",
143  "approximation order of Lagrange multipliers", "", 1,
144  &order_lambda, PETSC_NULL);
145 
146  CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", nb_sub_steps,
147  &nb_sub_steps, PETSC_NULL);
148 
149  CHKERR PetscOptionsBool("-my_is_partitioned",
150  "set if mesh is partitioned (this result that each "
151  "process keeps only part of the mes",
152  "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
153 
154  CHKERR PetscOptionsReal("-my_cn_value", "default regularisation cn value",
155  "", 1., &cn_value, PETSC_NULL);
156 
157  CHKERR PetscOptionsBool("-my_is_newton_cotes",
158  "set if Newton-Cotes quadrature rules are used", "",
159  PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
160  CHKERR PetscOptionsInt("-my_test_num", "test number", "", 0, &test_num,
161  PETSC_NULL);
162  CHKERR PetscOptionsBool("-my_convect", "set to convect integration pts", "",
163  PETSC_FALSE, &convect_pts, PETSC_NULL);
164  CHKERR PetscOptionsBool("-my_print_contact_state",
165  "output number of active gp at every iteration", "",
166  PETSC_FALSE, &print_contact_state, PETSC_NULL);
167  CHKERR PetscOptionsBool("-my_alm_flag",
168  "if set use ALM, if not use C-function", "",
169  PETSC_FALSE, &alm_flag, PETSC_NULL);
170  CHKERR PetscOptionsBool("-my_eigen_pos_flag",
171  "if set use eigen spatial positions are taken into "
172  "account for predeformed configuration",
173  "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
174 
175  CHKERR PetscOptionsReal("-my_scale_factor", "scale factor", "",
176  scale_factor, &scale_factor, PETSC_NULL);
177 
178  CHKERR PetscOptionsBool("-my_ignore_contact", "if set true, ignore contact",
179  "", PETSC_FALSE, &ignore_contact, PETSC_NULL);
180  CHKERR PetscOptionsBool("-my_ignore_pressure",
181  "if set true, ignore pressure", "", PETSC_FALSE,
182  &ignore_pressure, PETSC_NULL);
183  CHKERR PetscOptionsBool("-my_analytical_input",
184  "if set true, use analytical strain", "",
185  PETSC_TRUE, &analytical_input, PETSC_NULL);
186  CHKERR PetscOptionsBool("-my_save_mean_stress",
187  "if set true, save mean stress", "", PETSC_TRUE,
188  &save_mean_stress, PETSC_NULL);
189  CHKERR PetscOptionsString(
190  "-my_stress_tag_name", "stress tag name file name", "",
191  "INTERNAL_STRESS", stress_tag_name, 255, &flg_tag_name);
192  CHKERR PetscOptionsReal(
193  "-my_thermal_expansion_coef", "thermal expansion coef ", "",
194  thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
195  CHKERR PetscOptionsReal("-my_init_temp", "init_temp ", "", init_temp,
196  &init_temp, PETSC_NULL);
197 
198  CHKERR PetscOptionsReal("-my_mesh_height",
199  "vertical dimension of the mesh ", "", mesh_height,
200  &mesh_height, PETSC_NULL);
201  CHKERR PetscOptionsBool("-my_deform_flat", "if set true, deform flat", "",
202  PETSC_FALSE, &deform_flat_flag, PETSC_NULL);
203  CHKERR PetscOptionsReal("-my_flat_shift", "flat shift ", "", flat_shift,
204  &flat_shift, PETSC_NULL);
205 
206  CHKERR PetscOptionsBool("-my_wave_surf",
207  "if set true, make one of the surfaces wavy", "",
208  PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
209  CHKERR PetscOptionsInt("-my_wave_dim", "dimension (2 or 3)", "", wave_dim,
210  &wave_dim, PETSC_NULL);
211  CHKERR PetscOptionsReal("-my_wave_length", "profile wavelength", "",
212  wave_length, &wave_length, PETSC_NULL);
213  CHKERR PetscOptionsReal("-my_wave_ampl", "profile amplitude", "", wave_ampl,
214  &wave_ampl, PETSC_NULL);
215 
216  CHKERR PetscOptionsBool("-my_delete_prisms", "if set true, delete prisms",
217  "", PETSC_FALSE, &delete_prisms, PETSC_NULL);
218 
219  ierr = PetscOptionsEnd();
220  CHKERRQ(ierr);
221 
222  // Check if mesh file was provided
223  if (flg_file != PETSC_TRUE) {
224  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
225  }
226 
227  if (is_partitioned == PETSC_TRUE) {
228  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
229  "Partitioned mesh is not supported");
230  }
231 
232  const char *option;
233  option = "";
234  CHKERR moab.load_file(mesh_file_name, 0, option);
235 
236  // Create MoFEM database and link it to MoAB
237  MoFEM::Core core(moab);
238  MoFEM::Interface &m_field = core;
239 
240  Version file_ver;
241  CHKERR MoFEM::Interface::getFileVersion(moab, file_ver);
242  MOFEM_LOG("WORLD", Sev::inform) << "File version " << file_ver.strVersion();
243 
244  std::vector<BitRefLevel> bit_levels;
245  bit_levels.push_back(BitRefLevel().set(0));
246  auto bit_ref_manager = m_field.getInterface<BitRefManager>();
247  CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
248 
249  auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
250  MoFEMFunctionBegin;
251  auto prism_interface = m_field.getInterface<PrismInterface>();
252 
253  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
254  if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
255  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
256  cit->getName().c_str(), cit->getMeshsetId());
257  EntityHandle cubit_meshset = cit->getMeshset();
258 
259  // get tet entities from back bit_level
260  EntityHandle ref_level_meshset;
261  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
262  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
263  bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
264  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
265  bit_levels.back(), BitRefLevel().set(), MBPRISM,
266  ref_level_meshset);
267 
268  // get faces and tets to split
269  CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
270  true, 0);
271  // set new bit level
272  bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
273  // split faces and tets
274  CHKERR prism_interface->splitSides(ref_level_meshset,
275  bit_levels.back(), cubit_meshset,
276  true, true, 0);
277  // clean meshsets
278  CHKERR moab.delete_entities(&ref_level_meshset, 1);
279 
280  // update cubit meshsets
281  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
282  EntityHandle cubit_meshset = ciit->meshset;
283  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
284  cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
285  true);
286  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
287  cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE, true);
288  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
289  cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
290  CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
291  cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
292  }
293 
294  CHKERR bit_ref_manager->shiftRightBitRef(1);
295  bit_levels.pop_back();
296  }
297  }
298 
299  MoFEMFunctionReturn(0);
300  };
301 
302  auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
303  Range &contact_prisms, Range &master_tris,
304  Range &slave_tris) {
305  MoFEMFunctionBegin;
306 
307  EntityHandle meshset_prisms;
308  CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
309  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
310  bit_levels.back(), BitRefLevel().set(), MBPRISM, meshset_prisms);
311  CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
312  CHKERR moab.delete_entities(&meshset_prisms, 1);
313 
314  EntityHandle tri;
315  for (Range::iterator pit = contact_prisms.begin();
316  pit != contact_prisms.end(); pit++) {
317  CHKERR moab.side_element(*pit, 2, 3, tri);
318  master_tris.insert(tri);
319  CHKERR moab.side_element(*pit, 2, 4, tri);
320  slave_tris.insert(tri);
321  }
322 
323  MoFEMFunctionReturn(0);
324  };
325 
326  Range contact_prisms, master_tris, slave_tris;
327 
328  if (!ignore_contact) {
329  if (analytical_input) {
330  CHKERR add_prism_interface(bit_levels);
331  }
332  CHKERR find_contact_prisms(bit_levels, contact_prisms, master_tris,
333  slave_tris);
334  }
335 
336  auto deform_flat_surface = [&](int block_id, double shift, double height) {
337  MoFEMFunctionBegin;
338  Range all_tets, all_nodes;
339  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
340  if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
341  const int id = bit->getMeshsetId();
342  Range tets;
343  if (id == block_id) {
344  CHKERR m_field.get_moab().get_entities_by_dimension(
345  bit->getMeshset(), 3, tets, true);
346  all_tets.merge(tets);
347  }
348  }
349  }
350  CHKERR m_field.get_moab().get_connectivity(all_tets, all_nodes);
351  double coords[3];
352  for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
353  nit++) {
354  CHKERR moab.get_coords(&*nit, 1, coords);
355  double x = coords[0];
356  double y = coords[1];
357  double z = coords[2];
358  coords[2] -= shift;
359  CHKERR moab.set_coords(&*nit, 1, coords);
360  }
361  MoFEMFunctionReturn(0);
362  };
363 
364  auto make_wavy_surface = [&](int block_id, int dim, double lambda,
365  double delta, double height) {
366  MoFEMFunctionBegin;
367  Range all_tets, all_nodes;
368  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
369  if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
370  const int id = bit->getMeshsetId();
371  Range tets;
372  if (id == block_id) {
373  CHKERR m_field.get_moab().get_entities_by_dimension(
374  bit->getMeshset(), 3, tets, true);
375  all_tets.merge(tets);
376  }
377  }
378  }
379  CHKERR m_field.get_moab().get_connectivity(all_tets, all_nodes);
380  double coords[3];
381  for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
382  nit++) {
383  CHKERR moab.get_coords(&*nit, 1, coords);
384  double x = coords[0];
385  double y = coords[1];
386  double z = coords[2];
387  double coef = (height + z) / height;
388  switch (dim) {
389  case 2:
390  coords[2] -= coef * delta * (1. - cos(2. * M_PI * x / lambda));
391  break;
392  case 3:
393  coords[2] -=
394  coef * delta *
395  (1. - cos(2. * M_PI * x / lambda) * cos(2. * M_PI * y / lambda));
396  break;
397  default:
398  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
399  "Wrong dimension = %d", dim);
400  }
401 
402  CHKERR moab.set_coords(&*nit, 1, coords);
403  }
404  MoFEMFunctionReturn(0);
405  };
406 
407  if (deform_flat_flag && analytical_input) {
408  CHKERR deform_flat_surface(1, flat_shift, mesh_height);
409  CHKERR deform_flat_surface(2, -flat_shift, mesh_height);
410  }
411 
412  if (wave_surf_flag && analytical_input) {
413  CHKERR make_wavy_surface(1, wave_dim, wave_length, wave_ampl,
414  mesh_height);
415  // CHKERR make_wavy_surface(2, wave_dim, wave_length, -wave_ampl,
416  // mesh_height);
417  }
418 
419  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
420  MB_TAG_SPARSE, MF_ZERO);
421 
422  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
423  3, MB_TAG_SPARSE, MF_ZERO);
424  if (!eigen_pos_flag) {
425  CHKERR m_field.add_field("EIGEN_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
426  3, MB_TAG_SPARSE, MF_ZERO);
427  }
428 
429  CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
430  MB_TAG_SPARSE, MF_ZERO);
431 
432  // Declare problem add entities (by tets) to the field
433  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
434  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
435  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
436  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
437  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
438 
439  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
440  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
441  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
442  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
443  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
444 
445  if (!eigen_pos_flag) {
446  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
447  CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
448  CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
449  CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
450  CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_POSITIONS", 1);
451  }
452 
453  CHKERR m_field.add_ents_to_field_by_type(slave_tris, MBTRI, "LAGMULT");
454  CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
455  CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
456  CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
457 
458  auto set_contact_order = [&](Range &contact_prisms, int order_contact,
459  int nb_ho_levels) {
460  MoFEMFunctionBegin;
461  Range contact_tris, contact_edges;
462  CHKERR moab.get_adjacencies(contact_prisms, 2, false, contact_tris,
463  moab::Interface::UNION);
464  contact_tris = contact_tris.subset_by_type(MBTRI);
465  CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
466  moab::Interface::UNION);
467  Range ho_ents;
468  ho_ents.merge(contact_tris);
469  ho_ents.merge(contact_edges);
470  for (int ll = 0; ll < nb_ho_levels; ll++) {
471  Range ents, verts, tets;
472  CHKERR moab.get_connectivity(ho_ents, verts, true);
473  CHKERR moab.get_adjacencies(verts, 3, false, tets,
474  moab::Interface::UNION);
475  tets = tets.subset_by_type(MBTET);
476  for (auto d : {1, 2}) {
477  CHKERR moab.get_adjacencies(tets, d, false, ents,
478  moab::Interface::UNION);
479  }
480  ho_ents = unite(ho_ents, ents);
481  ho_ents = unite(ho_ents, tets);
482  }
483 
484  CHKERR m_field.set_field_order(ho_ents, "SPATIAL_POSITION",
485  order_contact);
486 
487  MoFEMFunctionReturn(0);
488  };
489 
490  if (!ignore_contact && order_contact > order) {
491  CHKERR set_contact_order(contact_prisms, order_contact, nb_ho_levels);
492  }
493 
494  // build field
495  CHKERR m_field.build_fields();
496 
497  // Projection on "x" field
498  {
499  Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
500  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
501  }
502  // Projection on "X" field
503  {
504  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
505  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
506  }
507 
508  Range slave_verts;
509  CHKERR moab.get_connectivity(slave_tris, slave_verts, true);
510  CHKERR m_field.getInterface<FieldBlas>()->setField(0.0, MBVERTEX,
511  slave_verts, "LAGMULT");
512 
513  // Add elastic element
514  boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
515  boost::make_shared<std::map<int, BlockData>>();
516  CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
517 
518  boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
519  new VolumeElementForcesAndSourcesCore(m_field));
520  boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
521  new VolumeElementForcesAndSourcesCore(m_field));
522  fe_elastic_lhs_ptr->getRuleHook = VolRule();
523  fe_elastic_rhs_ptr->getRuleHook = VolRule();
524 
525  CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
526  "SPATIAL_POSITION",
527  "MESH_NODE_POSITIONS", false);
528 
529  auto data_hooke_element_at_pts =
530  boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
531  CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
532  block_sets_ptr, "SPATIAL_POSITION",
533  "MESH_NODE_POSITIONS", false, false,
534  MBTET, data_hooke_element_at_pts);
535  auto thermal_strain =
536  [&thermal_expansion_coef, &init_temp, &test_num](
537  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords) {
538  FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
539  FTensor::Index<'i', 3> i;
540  FTensor::Index<'j', 3> j;
541 
542  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
543  double temp;
544  t_thermal_strain(i, j) = 0.0;
545 
546  switch (test_num) {
547  case 0:
548  // Put here analytical formula which may depend on coordinates
549  temp = init_temp - 1.0;
550  t_thermal_strain(i, j) =
551  thermal_expansion_coef * (temp - init_temp) * t_kd(i, j);
552  break;
553  case 1:
554  case 2:
555  t_thermal_strain(i, j) = -thermal_expansion_coef * t_kd(i, j);
556  break;
557  case 3:
558  t_thermal_strain(2, 2) = thermal_expansion_coef;
559  break;
560  case 4:
561  t_thermal_strain(i, j) = thermal_expansion_coef * t_kd(i, j);
562  break;
563  default:
564  break;
565  }
566  return t_thermal_strain;
567  };
568 
569  if (analytical_input) {
570  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
571  new HookeElement::OpAnalyticalInternalStrain_dx<0>(
572  "SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
573  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
574  new HookeInternalStressElement::OpGetAnalyticalInternalStress<0>(
575  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
576  thermal_strain));
577  } else {
578  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
579  new HookeInternalStressElement::OpGetInternalStress(
580  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
581  moab, stress_tag_name));
582  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
583  new HookeInternalStressElement::OpInternalStrain_dx(
584  "SPATIAL_POSITION", data_hooke_element_at_pts));
585  }
586 
587  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
588  new OpCalculateVectorFieldGradient<3, 3>(
589  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
590  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
591  new OpCalculateVectorFieldGradient<3, 3>(
592  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
593  fe_elastic_rhs_ptr->getOpPtrVector().push_back(
594  new HookeInternalStressElement::OpSaveStress(
595  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
596  *block_sets_ptr.get(), moab, scale_factor, save_mean_stress, false,
597  false));
598 
599  Range all_tets;
600  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
601  if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
602  Range tets;
603  CHKERR moab.get_entities_by_dimension(bit->getMeshset(), 3, tets, true);
604  all_tets.merge(tets);
605  }
606  }
607  Skinner skinner(&moab);
608  Range skin_tris;
609  CHKERR skinner.find_skin(0, all_tets, false, skin_tris);
610 
611  CHKERR m_field.add_finite_element("SKIN", MF_ZERO);
612  CHKERR m_field.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
613  CHKERR m_field.modify_finite_element_add_field_row("SKIN",
614  "SPATIAL_POSITION");
615  CHKERR m_field.modify_finite_element_add_field_col("SKIN",
616  "SPATIAL_POSITION");
617  CHKERR m_field.modify_finite_element_add_field_data("SKIN",
618  "SPATIAL_POSITION");
619  CHKERR m_field.modify_finite_element_add_field_data("SKIN",
620  "MESH_NODE_POSITIONS");
621  CHKERR m_field.modify_finite_element_add_field_data("SKIN",
622  "EIGEN_POSITIONS");
623 
624  auto make_contact_element = [&]() {
625  return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
626  m_field);
627  };
628 
629  auto make_convective_master_element = [&]() {
630  return boost::make_shared<
631  SimpleContactProblem::ConvectMasterContactElement>(
632  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
633  };
634 
635  auto make_convective_slave_element = [&]() {
636  return boost::make_shared<
637  SimpleContactProblem::ConvectSlaveContactElement>(
638  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
639  };
640 
641  auto make_contact_common_data = [&]() {
642  return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
643  m_field);
644  };
645 
646  auto get_contact_rhs = [&](auto contact_problem, auto make_element,
647  bool is_alm = false) {
648  auto fe_rhs_simple_contact = make_element();
649  auto common_data_simple_contact = make_contact_common_data();
650  if (print_contact_state) {
651  fe_rhs_simple_contact->contactStateVec =
652  common_data_simple_contact->gaussPtsStateVec;
653  }
654  contact_problem->setContactOperatorsRhs(
655  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
656  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
657  return fe_rhs_simple_contact;
658  };
659 
660  auto get_master_traction_rhs = [&](auto contact_problem, auto make_element,
661  bool is_alm = false) {
662  auto fe_rhs_simple_contact = make_element();
663  auto common_data_simple_contact = make_contact_common_data();
664  contact_problem->setMasterForceOperatorsRhs(
665  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
666  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
667  return fe_rhs_simple_contact;
668  };
669 
670  auto get_master_traction_lhs = [&](auto contact_problem, auto make_element,
671  bool is_alm = false) {
672  auto fe_lhs_simple_contact = make_element();
673  auto common_data_simple_contact = make_contact_common_data();
674  contact_problem->setMasterForceOperatorsLhs(
675  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
676  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
677  return fe_lhs_simple_contact;
678  };
679 
680  auto get_contact_lhs = [&](auto contact_problem, auto make_element,
681  bool is_alm = false) {
682  auto fe_lhs_simple_contact = make_element();
683  auto common_data_simple_contact = make_contact_common_data();
684  contact_problem->setContactOperatorsLhs(
685  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
686  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS", true);
687  return fe_lhs_simple_contact;
688  };
689 
690  auto cn_value_ptr = boost::make_shared<double>(cn_value);
691  auto contact_problem = boost::make_shared<SimpleContactProblem>(
692  m_field, cn_value_ptr, is_newton_cotes);
693 
694  // add fields to the global matrix by adding the element
695  if (!eigen_pos_flag)
696  contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
697  "LAGMULT", contact_prisms);
698  else
699  contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
700  "LAGMULT", contact_prisms,
701  eigen_pos_flag, "EIGEN_POSITIONS");
702 
703  contact_problem->addPostProcContactElement(
704  "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
705  "MESH_NODE_POSITIONS", slave_tris);
706 
707  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "SPATIAL_POSITION");
708 
709  // Add spring boundary condition applied on surfaces.
710  CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
711  "MESH_NODE_POSITIONS");
712 
713  // build finite elemnts
714  CHKERR m_field.build_finite_elements();
715 
716  // build adjacencies
717  CHKERR m_field.build_adjacencies(bit_levels.back());
718 
719  // define problems
720  CHKERR m_field.add_problem("CONTACT_PROB", MF_ZERO);
721 
722  // set refinement level for problem
723  CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
724  bit_levels.back());
725 
726  DMType dm_name = "DMMOFEM";
727  CHKERR DMRegister_MoFEM(dm_name);
728 
729  SmartPetscObj<DM> dm;
730  dm = createSmartDM(m_field.get_comm(), dm_name);
731 
732  // create dm instance
733  CHKERR DMSetType(dm, dm_name);
734 
735  // set dm datastruture which created mofem datastructures
736  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_levels.back());
737  CHKERR DMSetFromOptions(dm);
738  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
739  // add elements to dm
740  CHKERR DMMoFEMAddElement(dm, "CONTACT_ELEM");
741  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
742  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
743  CHKERR DMMoFEMAddElement(dm, "SPRING");
744  CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
745  CHKERR DMMoFEMAddElement(dm, "SKIN");
746 
747  CHKERR DMSetUp(dm);
748 
749  // Vector of DOFs and the RHS
750  auto D = smartCreateDMVector(dm);
751  auto F = smartVectorDuplicate(D);
752 
753  // Stiffness matrix
754  auto Aij = smartCreateDMMatrix(dm);
755 
756  CHKERR VecZeroEntries(D);
757  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
758  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
759  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
760 
761  CHKERR VecZeroEntries(F);
762  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
763  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
764 
765  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
766  CHKERR MatZeroEntries(Aij);
767 
768  fe_elastic_rhs_ptr->snes_f = F;
769  fe_elastic_lhs_ptr->snes_B = Aij;
770 
771  // Dirichlet BC
772  boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
773  boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
774  m_field, "SPATIAL_POSITION", Aij, D, F));
775 
776  dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
777  dirichlet_bc_ptr->snes_x = D;
778 
779  // Assemble pressure and traction forces
780  boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
781  CHKERR MetaNeumannForces::setMomentumFluxOperators(
782  m_field, neumann_forces, NULL, "SPATIAL_POSITION");
783 
784  boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
785  neumann_forces.begin();
786  for (; mit != neumann_forces.end(); mit++) {
787  mit->second->methodsOp.push_back(new SimpleContactProblem::LoadScale());
788  CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
789  &mit->second->getLoopFe(), NULL, NULL);
790  }
791 
792  // Implementation of spring element
793  // Create new instances of face elements for springs
794  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
795  new FaceElementForcesAndSourcesCore(m_field));
796  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
797  new FaceElementForcesAndSourcesCore(m_field));
798 
799  CHKERR MetaSpringBC::setSpringOperators(
800  m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
801  "MESH_NODE_POSITIONS");
802 
803  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
804  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
805  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
806  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
807  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL,
808  dirichlet_bc_ptr.get(), NULL);
809  if (convect_pts == PETSC_TRUE) {
810  CHKERR DMMoFEMSNESSetFunction(
811  dm, "CONTACT_ELEM",
812  get_contact_rhs(contact_problem, make_convective_master_element),
813  PETSC_NULL, PETSC_NULL);
814  CHKERR DMMoFEMSNESSetFunction(
815  dm, "CONTACT_ELEM",
816  get_master_traction_rhs(contact_problem,
817  make_convective_slave_element),
818  PETSC_NULL, PETSC_NULL);
819  } else {
820  CHKERR DMMoFEMSNESSetFunction(
821  dm, "CONTACT_ELEM",
822  get_contact_rhs(contact_problem, make_contact_element, alm_flag),
823  PETSC_NULL, PETSC_NULL);
824  CHKERR DMMoFEMSNESSetFunction(
825  dm, "CONTACT_ELEM",
826  get_master_traction_rhs(contact_problem, make_contact_element,
827  alm_flag),
828  PETSC_NULL, PETSC_NULL);
829  }
830 
831  CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", fe_elastic_rhs_ptr, PETSC_NULL,
832  PETSC_NULL);
833 
834  CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
835  PETSC_NULL);
836  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, PETSC_NULL, PETSC_NULL,
837  dirichlet_bc_ptr.get());
838 
839  boost::shared_ptr<FEMethod> fe_null;
840  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
841  fe_null);
842  if (convect_pts == PETSC_TRUE) {
843  CHKERR DMMoFEMSNESSetJacobian(
844  dm, "CONTACT_ELEM",
845  get_contact_lhs(contact_problem, make_convective_master_element),
846  PETSC_NULL, PETSC_NULL);
847  CHKERR DMMoFEMSNESSetJacobian(
848  dm, "CONTACT_ELEM",
849  get_master_traction_lhs(contact_problem,
850  make_convective_slave_element),
851  PETSC_NULL, PETSC_NULL);
852  } else {
853  CHKERR DMMoFEMSNESSetJacobian(
854  dm, "CONTACT_ELEM",
855  get_contact_lhs(contact_problem, make_contact_element, alm_flag),
856  PETSC_NULL, PETSC_NULL);
857  CHKERR DMMoFEMSNESSetJacobian(
858  dm, "CONTACT_ELEM",
859  get_master_traction_lhs(contact_problem, make_contact_element,
860  alm_flag),
861  PETSC_NULL, PETSC_NULL);
862  }
863  CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", fe_elastic_lhs_ptr, PETSC_NULL,
864  PETSC_NULL);
865  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
866  PETSC_NULL);
867  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, fe_null,
868  dirichlet_bc_ptr);
869 
870  if (test_num) {
871  char testing_options[] = "-ksp_type fgmres "
872  "-pc_type lu "
873  "-pc_factor_mat_solver_type mumps "
874  "-snes_type newtonls "
875  "-snes_linesearch_type basic "
876  "-snes_max_it 10 "
877  "-snes_atol 1e-8 "
878  "-snes_rtol 1e-8 ";
879  CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
880  }
881 
882  auto snes = MoFEM::createSNES(m_field.get_comm());
883  CHKERR SNESSetDM(snes, dm);
884  SnesCtx *snes_ctx;
885  // create snes nonlinear solver
886  {
887  CHKERR SNESSetDM(snes, dm);
888  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
889  CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
890  CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
891  CHKERR SNESSetFromOptions(snes);
892  }
893 
894  /// Post proc on the skin
895  PostProcFaceOnRefinedMesh post_proc_skin(m_field);
896  CHKERR post_proc_skin.generateReferenceElementMesh();
897  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false, false);
898  CHKERR post_proc_skin.addFieldValuesPostProc("SPATIAL_POSITION");
899  CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
900  CHKERR post_proc_skin.addFieldValuesPostProc("EIGEN_POSITIONS");
901 
902  struct OpGetFieldGradientValuesOnSkin
903  : public FaceElementForcesAndSourcesCore::UserDataOperator {
904 
905  const std::string feVolName;
906  boost::shared_ptr<VolSideFe> sideOpFe;
907 
908  OpGetFieldGradientValuesOnSkin(const std::string field_name,
909  const std::string vol_fe_name,
910  boost::shared_ptr<VolSideFe> side_fe)
911  : FaceElementForcesAndSourcesCore::UserDataOperator(
912  field_name, UserDataOperator::OPROW),
913  feVolName(vol_fe_name), sideOpFe(side_fe) {}
914 
915  MoFEMErrorCode doWork(int side, EntityType type,
916  EntitiesFieldData::EntData &data) {
917  MoFEMFunctionBegin;
918  if (type != MBVERTEX)
919  MoFEMFunctionReturnHot(0);
920  CHKERR loopSideVolumes(feVolName, *sideOpFe);
921  MoFEMFunctionReturn(0);
922  }
923  };
924 
925  boost::shared_ptr<VolSideFe> my_vol_side_fe_ptr =
926  boost::make_shared<VolSideFe>(m_field);
927  my_vol_side_fe_ptr->getOpPtrVector().push_back(
928  new OpCalculateVectorFieldGradient<3, 3>(
929  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
930  my_vol_side_fe_ptr->getOpPtrVector().push_back(
931  new OpCalculateVectorFieldGradient<3, 3>(
932  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
933 
934  post_proc_skin.getOpPtrVector().push_back(
935  new OpGetFieldGradientValuesOnSkin("SPATIAL_POSITION", "ELASTIC",
936  my_vol_side_fe_ptr));
937  post_proc_skin.getOpPtrVector().push_back(
938  new OpCalculateVectorFieldValues<3>(
939  "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
940  post_proc_skin.getOpPtrVector().push_back(
941  new OpCalculateVectorFieldValues<3>(
942  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
943 
944  post_proc_skin.getOpPtrVector().push_back(
945  new HookeInternalStressElement::OpPostProcHookeElement<
946  FaceElementForcesAndSourcesCore>(
947  "SPATIAL_POSITION", data_hooke_element_at_pts,
948  *block_sets_ptr.get(), post_proc_skin.postProcMesh,
949  post_proc_skin.mapGaussPts, false, false));
950 
951  for (int ss = 0; ss != nb_sub_steps; ++ss) {
952  if (!ignore_pressure) {
953  SimpleContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
954  } else {
955  SimpleContactProblem::LoadScale::lAmbda = 0;
956  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Ignoring pressure...\n");
957  }
958  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
959  SimpleContactProblem::LoadScale::lAmbda);
960 
961  CHKERR SNESSolve(snes, PETSC_NULL, D);
962 
963  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
964  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
965  }
966 
967  // save on mesh
968  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
969 
970  SmartPetscObj<Vec> v_energy;
971  CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "SPATIAL_POSITION",
972  "MESH_NODE_POSITIONS", false, false,
973  v_energy);
974  const double *eng_ptr;
975  CHKERR VecGetArrayRead(v_energy, &eng_ptr);
976  // Print elastic energy
977  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy: %8.8e\n", *eng_ptr);
978 
979  {
980  PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
981  CHKERR DMoFEMLoopFiniteElements(dm, "SKIN", &post_proc_skin);
982  ostringstream stm;
983  string out_file_name;
984  stm << "out_skin.h5m";
985  out_file_name = stm.str();
986  PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n", out_file_name.c_str());
987  CHKERR post_proc_skin.writeFile(stm.str());
988  }
989 
990  // moab_instance
991  moab::Core mb_post; // create database
992  moab::Interface &moab_proc = mb_post; // create interface to database
993 
994  auto common_data_simple_contact = make_contact_common_data();
995 
996  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
997  fe_post_proc_simple_contact;
998  if (convect_pts == PETSC_TRUE) {
999  fe_post_proc_simple_contact = make_convective_master_element();
1000  } else {
1001  fe_post_proc_simple_contact = make_contact_element();
1002  }
1003 
1004  std::array<double, 2> nb_gauss_pts;
1005  std::array<double, 2> contact_area;
1006 
1007  if (!ignore_contact) {
1008  contact_problem->setContactOperatorsForPostProc(
1009  fe_post_proc_simple_contact, common_data_simple_contact, m_field,
1010  "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1011  "EIGEN_POSITIONS", true);
1012 
1013  mb_post.delete_mesh();
1014 
1015  CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
1016  CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
1017 
1018  CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_ELEM",
1019  fe_post_proc_simple_contact);
1020 
1021  auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
1022  MoFEMFunctionBegin;
1023  CHKERR VecAssemblyBegin(vec);
1024  CHKERR VecAssemblyEnd(vec);
1025  const double *array;
1026  CHKERR VecGetArrayRead(vec, &array);
1027  if (m_field.get_comm_rank() == 0) {
1028  for (int i : {0, 1})
1029  data[i] = array[i];
1030  }
1031  CHKERR VecRestoreArrayRead(vec, &array);
1032  MoFEMFunctionReturn(0);
1033  };
1034 
1035  CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
1036  nb_gauss_pts);
1037  CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
1038  contact_area);
1039 
1040  if (m_field.get_comm_rank() == 0) {
1041  PetscPrintf(PETSC_COMM_SELF, "Active gauss pts: %d out of %d\n",
1042  (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
1043 
1044  PetscPrintf(PETSC_COMM_SELF,
1045  "Active contact area: %8.8f out of %8.8f (%8.8f% %)\n",
1046  contact_area[0], contact_area[1],
1047  contact_area[0] / contact_area[1] * 100.);
1048  }
1049 
1050  string out_file_name;
1051  std::ostringstream strm;
1052  strm << "out_contact_integ_pts"
1053  << ".h5m";
1054  out_file_name = strm.str();
1055  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
1056  out_file_name.c_str());
1057  CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
1058  "PARALLEL=WRITE_PART");
1059  }
1060 
1061  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
1062  new PostProcFaceOnRefinedMesh(m_field));
1063 
1064  CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
1065  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *post_proc_contact_ptr, false,
1066  false);
1067  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("LAGMULT");
1068  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("SPATIAL_POSITION");
1069  CHKERR post_proc_contact_ptr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
1070 
1071  if (!ignore_contact) {
1072  CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
1073  post_proc_contact_ptr);
1074  string out_file_name;
1075  std::ostringstream stm;
1076  stm << "out_contact"
1077  << ".h5m";
1078  out_file_name = stm.str();
1079  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
1080  out_file_name.c_str());
1081  CHKERR post_proc_contact_ptr->postProcMesh.write_file(
1082  out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
1083  }
1084 
1085  CHKERR m_field.getInterface<FieldBlas>()->fieldCopy(1., "SPATIAL_POSITION",
1086  "EIGEN_POSITIONS");
1087 
1088  const int n_parts = m_field.get_comm_size();
1089  if (m_field.get_comm_rank() == 0) {
1090  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
1091  dm, "ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
1092 
1093  if (delete_prisms) {
1094  Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
1095 
1096  CHKERR moab.get_adjacencies(contact_prisms, 2, true, faces,
1097  moab::Interface::UNION);
1098  tris = faces.subset_by_type(MBTRI);
1099  quads = faces.subset_by_type(MBQUAD);
1100  CHKERR moab.get_adjacencies(tris, 1, true, tris_edges,
1101  moab::Interface::UNION);
1102  CHKERR moab.get_adjacencies(quads, 1, true, quads_edges,
1103  moab::Interface::UNION);
1104 
1105  ents_to_delete.merge(contact_prisms);
1106  ents_to_delete.merge(quads);
1107  ents_to_delete.merge(subtract(quads_edges, tris_edges));
1108 
1109  CHKERR moab.delete_entities(ents_to_delete);
1110  }
1111  if (flg_file_out) {
1112  PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n", output_mesh_name);
1113  CHKERR moab.write_file(output_mesh_name, "MOAB");
1114  }
1115 
1116  auto get_tag_handle = [&](auto name, auto size) {
1117  Tag th;
1118  std::vector<double> def_vals(size, 0.0);
1119  CHKERR moab.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1120  MB_TAG_CREAT | MB_TAG_SPARSE,
1121  def_vals.data());
1122  return th;
1123  };
1124 
1125  if (test_num) {
1126  Range tets;
1127  CHKERR moab.get_entities_by_dimension(0, 3, tets);
1128  EntityHandle tet = tets.front();
1129  std::array<double, 9> internal_stress, actual_stress;
1130  std::array<double, 9> internal_stress_ref, actual_stress_ref;
1131  std::array<double, 2> nb_gauss_pts_ref, contact_area_ref;
1132  switch (test_num) {
1133  case 1:
1134  internal_stress_ref = {5., 5., 5., 0., 0., 0., 0., 0., 0.};
1135  actual_stress_ref = {0., 0., 1., 0., 0., 0., 0., 0., 0.};
1136  break;
1137  case 2:
1138  internal_stress_ref = {5., 5., 5., 0., 0., 0., 0., 0., 0.};
1139  actual_stress_ref = {0., 5. / 3., 5. / 3., 0., 0., 0., 0., 0., 0.};
1140  break;
1141  case 3:
1142  actual_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1143  if (strcmp(stress_tag_name, "INTERNAL_STRESS") == 0)
1144  internal_stress_ref = {0., 0., -200., 0., 0., 0., 0., 0., 0.};
1145  else
1146  internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1147  break;
1148  case 4:
1149  nb_gauss_pts_ref = {96, 192};
1150  contact_area_ref = {0.125, 0.25};
1151  break;
1152  default:
1153  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Test number %d not found",
1154  test_num);
1155  }
1156 
1157  auto th_internal_stress = get_tag_handle("MED_INTERNAL_STRESS", 9);
1158  auto th_actual_stress = get_tag_handle("MED_ACTUAL_STRESS", 9);
1159  CHKERR moab.tag_get_data(th_internal_stress, &tet, 1,
1160  internal_stress.data());
1161  CHKERR moab.tag_get_data(th_actual_stress, &tet, 1,
1162  actual_stress.data());
1163  const double eps = 1e-10;
1164  if (test_num == 4) {
1165  if (std::abs(nb_gauss_pts_ref[0] - nb_gauss_pts[0]) > eps) {
1166  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1167  "Wrong number of active contact gauss pts: should be %d "
1168  "but is %d",
1169  (int)nb_gauss_pts_ref[0], (int)nb_gauss_pts[0]);
1170  }
1171  if (std::abs(nb_gauss_pts_ref[1] - nb_gauss_pts[1]) > eps) {
1172  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1173  "Wrong total number of contact gauss pts: should be %d "
1174  "but is %d",
1175  (int)nb_gauss_pts_ref[1], (int)nb_gauss_pts[1]);
1176  }
1177  if (std::abs(contact_area_ref[0] - contact_area[0]) > eps) {
1178  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1179  "Wrong active contact area: should be %g "
1180  "but is %g",
1181  contact_area_ref[0], contact_area[0]);
1182  }
1183  if (std::abs(contact_area_ref[1] - contact_area[1]) > eps) {
1184  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1185  "Wrong potential contact area: should be %g "
1186  "but is %g",
1187  contact_area_ref[1], contact_area[1]);
1188  }
1189  } else {
1190  if (save_mean_stress) {
1191  for (int i = 0; i < 9; i++) {
1192  if (std::abs(internal_stress[i] - internal_stress_ref[i]) > eps) {
1193  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1194  "Wrong component %d of internal stress: should be %g "
1195  "but is %g",
1196  i, internal_stress_ref[i], internal_stress[i]);
1197  }
1198  if (std::abs(actual_stress[i] - actual_stress_ref[i]) > eps) {
1199  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1200  "Wrong component %d of actual stress: should be %g "
1201  "but is %g",
1202  i, actual_stress_ref[i], actual_stress[i]);
1203  }
1204  }
1205  }
1206  }
1207  }
1208  }
1209  }
1210  CATCH_ERRORS;
1211 
1212  // finish work cleaning memory, getting statistics, etc
1213  MoFEM::Core::Finalize();
1214 
1215  return 0;
1216 }
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:703
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
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.
main
int main(int argc, char *argv[])
Definition: simple_contact_thermal.cpp:39
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
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
HookeInternalStressElement::OpSaveStress
Definition: HookeInternalStressElement.hpp:117
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
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:468
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
SimpleContactProblem::LoadScale::lAmbda
static double lAmbda
Definition: SimpleContact.hpp:37
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
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.
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
BlockData
NonlinearElasticElement::BlockData BlockData
Definition: simple_contact_thermal.cpp:32
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::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
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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: simple_contact_thermal.cpp:36
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::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:110
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
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:128
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:139
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
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
help
static char help[]
Definition: simple_contact_thermal.cpp:30
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:1561
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
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
std
Definition: enable_if.hpp:5
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
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
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
SimpleContactProblem::ConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:181
MoFEM::BlockData
Definition: MeshsetsManager.cpp:766
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
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:142
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:396
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