v0.14.0
test_jacobian_of_simple_contact_element.cpp
Go to the documentation of this file.
1 /** \file test_jacobian_of_simple_contact_element.cpp
2  * \example test_jacobian_of_simple_contact_element.cpp
3  *
4  * Testing implementation of simple contact element (for contact between
5  * surfaces with matching meshes) by verifying its tangent matrix
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 
25 using namespace std;
26 using namespace MoFEM;
27 
28 static char help[] = "\n";
29 int main(int argc, char *argv[]) {
30 
31  // Initialize MoFEM
32  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
33 
34  // Create mesh database
35  moab::Core mb_instance; // create database
36  moab::Interface &moab = mb_instance; // create interface to database
37 
38  try {
39  PetscBool flg_file;
40  char mesh_file_name[255];
41  PetscInt order = 1;
42  PetscInt order_lambda = 1;
43  PetscReal r_value = 1.;
44  PetscReal cn_value = 1.;
45  PetscBool is_newton_cotes = PETSC_FALSE;
46  PetscBool test_jacobian = PETSC_FALSE;
47  PetscBool convect_pts = PETSC_FALSE;
48  PetscBool test_ale = PETSC_FALSE;
49  PetscBool alm_flag = PETSC_FALSE;
50  PetscBool eigen_pos_flag = PETSC_FALSE;
51  PetscBool use_reference_coordinates = PETSC_TRUE;
52 
53  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
54 
55  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_jacobian", &test_jacobian,
56  PETSC_NULL);
57 
58  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
59  mesh_file_name, 255, &flg_file);
60 
61  CHKERR PetscOptionsInt("-my_order",
62  "approximation order for spatial positions", "", 1,
63  &order, PETSC_NULL);
64  CHKERR PetscOptionsInt(
65  "-my_order_lambda",
66  "default approximation order of Lagrange multipliers", "", 1,
67  &order_lambda, PETSC_NULL);
68 
69  CHKERR PetscOptionsReal("-my_cn_value", "default regularisation cn value",
70  "", 1., &cn_value, PETSC_NULL);
71 
72  CHKERR PetscOptionsBool("-my_is_newton_cotes",
73  "set if Newton-Cotes integration rules are used",
74  "", PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
75  CHKERR PetscOptionsBool("-my_convect", "set to convect integration pts", "",
76  PETSC_FALSE, &convect_pts, PETSC_NULL);
77  CHKERR PetscOptionsBool("-my_alm_flag", "set to convect integration pts",
78  "", PETSC_FALSE, &alm_flag, PETSC_NULL);
79 
80  CHKERR PetscOptionsBool("-my_eigen_pos_flag",
81  "if set use eigen spatial positions are taken into "
82  "account for predeformed configuration",
83  "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
84 
85  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_ale", &test_ale,
86  PETSC_NULL);
87 
88  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-my_use_reference_coordinates",
89  &use_reference_coordinates, PETSC_NULL);
90 
91  ierr = PetscOptionsEnd();
92  CHKERRQ(ierr);
93 
94  // Check if mesh file was provided
95  if (flg_file != PETSC_TRUE) {
96  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
97  }
98 
99  // Read mesh to MOAB
100  const char *option;
101  option = "";
102  CHKERR moab.load_file(mesh_file_name, 0, option);
103 
104  // Create MoFEM database and link it to MoAB
105  MoFEM::Core core(moab);
106  MoFEM::Interface &m_field = core;
107 
108  auto add_prism_interface = [&](Range &contact_prisms, Range &master_tris,
109  Range &slave_tris,
110  std::vector<BitRefLevel> &bit_levels) {
111  MoFEMFunctionBegin;
112  PrismInterface *interface;
113  CHKERR m_field.getInterface(interface);
114 
115  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
116  if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
117  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
118  cit->getName().c_str(), cit->getMeshsetId());
119  EntityHandle cubit_meshset = cit->getMeshset();
120 
121  // get tet entities from back bit_level
122  EntityHandle ref_level_meshset;
123  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
124  CHKERR m_field.getInterface<BitRefManager>()
125  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
126  BitRefLevel().set(), MBTET,
127  ref_level_meshset);
128  CHKERR m_field.getInterface<BitRefManager>()
129  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
130  BitRefLevel().set(), MBPRISM,
131  ref_level_meshset);
132 
133  // get faces and tets to split
134  CHKERR interface->getSides(cubit_meshset, bit_levels.back(), true, 0);
135  // set new bit level
136  bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
137  // split faces and tets
138  CHKERR interface->splitSides(ref_level_meshset, bit_levels.back(),
139  cubit_meshset, true, true, 0);
140  // clean meshsets
141  CHKERR moab.delete_entities(&ref_level_meshset, 1);
142 
143  CHKERR m_field.getInterface<BitRefManager>()->shiftRightBitRef(1);
144  bit_levels.pop_back();
145  }
146  }
147 
148  EntityHandle meshset_prisms;
149  CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
150  CHKERR m_field.getInterface<BitRefManager>()
151  ->getEntitiesByTypeAndRefLevel(bit_levels.back(), BitRefLevel().set(),
152  MBPRISM, meshset_prisms);
153  CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
154  CHKERR moab.delete_entities(&meshset_prisms, 1);
155 
156  EntityHandle tri;
157  for (Range::iterator pit = contact_prisms.begin();
158  pit != contact_prisms.end(); pit++) {
159  CHKERR moab.side_element(*pit, 2, 3, tri);
160  master_tris.insert(tri);
161  CHKERR moab.side_element(*pit, 2, 4, tri);
162  slave_tris.insert(tri);
163  }
164 
165  MoFEMFunctionReturn(0);
166  };
167 
168  Range contact_prisms, master_tris, slave_tris;
169  std::vector<BitRefLevel> bit_levels;
170 
171  bit_levels.push_back(BitRefLevel().set(0));
172  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
173  0, 3, bit_levels.back());
174 
175  CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
176  bit_levels);
177 
178  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
179  MB_TAG_SPARSE, MF_ZERO);
180 
181  // Declare problem
182  // add entities (by tets) to the field
183  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
184  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
185  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
186  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
187  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
188 
189  CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
190  MB_TAG_SPARSE, MF_ZERO);
191 
192  CHKERR m_field.add_ents_to_field_by_type(slave_tris, MBTRI, "LAGMULT");
193  CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
194  CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
195  CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
196 
197  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
198  3, MB_TAG_SPARSE, MF_ZERO);
199 
200  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
201  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
202  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
203  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
204  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
205 
206  if (eigen_pos_flag) {
207  CHKERR m_field.add_field("EIGEN_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
208  3, MB_TAG_SPARSE, MF_ZERO);
209  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
210  CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
211  CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
212  CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
213  CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_POSITIONS", 1);
214  }
215 
216  // build field
217  CHKERR m_field.build_fields();
218 
219  PetscRandom rctx;
220  PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
221 
222  auto set_coord = [&](VectorAdaptor &&field_data, double *x, double *y,
223  double *z) {
224  MoFEMFunctionBegin;
225  double value;
226  double scale = 0.5;
227  PetscRandomGetValue(rctx, &value);
228  field_data[0] = (*x) + (value - 0.5) * scale;
229  PetscRandomGetValue(rctx, &value);
230  field_data[1] = (*y) + (value - 0.5) * scale;
231  PetscRandomGetValue(rctx, &value);
232  field_data[2] = (*z) + (value - 0.5) * scale;
233  MoFEMFunctionReturn(0);
234  };
235 
236  auto set_pressure = [&](VectorAdaptor &&field_data, double *x, double *y,
237  double *z) {
238  MoFEMFunctionBegin;
239  double value;
240  double scale = 1.0;
241  PetscRandomGetValueReal(rctx, &value);
242  field_data[0] = value * scale;
243  MoFEMFunctionReturn(0);
244  };
245 
246  CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(set_coord,
247  "SPATIAL_POSITION");
248  CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(set_pressure,
249  "LAGMULT");
250 
251  if (eigen_pos_flag) {
252  CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(
253  set_coord, "SPATIAL_POSITION");
254  }
255 
256  if (test_ale == PETSC_TRUE) {
257  CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(
258  set_coord, "MESH_NODE_POSITIONS");
259  } else {
260  // MESH_NODE_POSITIONS
261  {
262  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
263  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
264  }
265  }
266 
267  PetscRandomDestroy(&rctx);
268 
269  auto cn_value_ptr = boost::make_shared<double>(cn_value);
270  auto contact_problem = boost::make_shared<SimpleContactProblem>(
271  m_field, cn_value_ptr, is_newton_cotes);
272 
273  auto make_contact_element = [&]() {
274  return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
275  m_field);
276  };
277 
278  auto make_convective_master_element = [&]() {
279  return boost::make_shared<
280  SimpleContactProblem::ConvectMasterContactElement>(
281  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
282  };
283 
284  auto make_convective_slave_element = [&]() {
285  return boost::make_shared<
286  SimpleContactProblem::ConvectSlaveContactElement>(
287  m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
288  };
289 
290  auto make_contact_common_data = [&]() {
291  return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
292  m_field);
293  };
294 
295  auto get_contact_rhs = [&](auto contact_problem, auto make_element,
296  bool is_alm = false) {
297  auto fe_rhs_simple_contact = make_element();
298  auto common_data_simple_contact = make_contact_common_data();
299  contact_problem->setContactOperatorsRhs(
300  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
301  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS",
302  use_reference_coordinates);
303  return fe_rhs_simple_contact;
304  };
305 
306  auto get_master_contact_lhs = [&](auto contact_problem, auto make_element,
307  bool is_alm = false) {
308  auto fe_lhs_simple_contact = make_element();
309  auto common_data_simple_contact = make_contact_common_data();
310  contact_problem->setContactOperatorsLhs(
311  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
312  "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS",
313  use_reference_coordinates);
314  return fe_lhs_simple_contact;
315  };
316 
317  auto get_master_traction_rhs = [&](auto contact_problem, auto make_element,
318  bool alm_flag = false) {
319  auto fe_rhs_simple_contact = make_element();
320  auto common_data_simple_contact = make_contact_common_data();
321  contact_problem->setMasterForceOperatorsRhs(
322  fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
323  "LAGMULT", alm_flag, eigen_pos_flag, "EIGEN_POSITIONS",
324  use_reference_coordinates);
325  return fe_rhs_simple_contact;
326  };
327 
328  auto get_master_traction_lhs = [&](auto contact_problem, auto make_element,
329  bool alm_flag = false) {
330  auto fe_lhs_simple_contact = make_element();
331  auto common_data_simple_contact = make_contact_common_data();
332  contact_problem->setMasterForceOperatorsLhs(
333  fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
334  "LAGMULT", alm_flag, eigen_pos_flag, "EIGEN_POSITIONS",
335  use_reference_coordinates);
336  return fe_lhs_simple_contact;
337  };
338 
339  auto get_contact_material_rhs = [&](auto contact_problem, auto make_element,
340  Range &ale_nodes) {
341  auto fe_rhs_simple_contact_ale_material = make_element();
342  auto common_data_simple_contact = make_contact_common_data();
343  common_data_simple_contact->forcesOnlyOnEntitiesRow.clear();
344  common_data_simple_contact->forcesOnlyOnEntitiesRow = ale_nodes;
345  contact_problem->setContactOperatorsRhsALEMaterial(
346  fe_rhs_simple_contact_ale_material, common_data_simple_contact,
347  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", "MATERIAL");
348  return fe_rhs_simple_contact_ale_material;
349  };
350 
351  auto get_simple_contact_ale_lhs = [&](auto contact_problem,
352  auto make_element) {
353  auto fe_lhs_simple_contact_ale = make_element();
354  auto common_data_simple_contact = make_contact_common_data();
355  contact_problem->setContactOperatorsLhsALE(
356  fe_lhs_simple_contact_ale, common_data_simple_contact,
357  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", eigen_pos_flag,
358  "EIGEN_POSITIONS");
359  return fe_lhs_simple_contact_ale;
360  };
361 
362  auto get_simple_contact_ale_material_lhs =
363  [&](auto contact_problem, auto make_element, Range &ale_nodes) {
364  auto fe_lhs_simple_contact_material_ale = make_element();
365  auto common_data_simple_contact = make_contact_common_data();
366  common_data_simple_contact->forcesOnlyOnEntitiesRow.clear();
367  common_data_simple_contact->forcesOnlyOnEntitiesRow = ale_nodes;
368  contact_problem->setContactOperatorsLhsALEMaterial(
369  fe_lhs_simple_contact_material_ale, common_data_simple_contact,
370  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", "MATERIAL");
371  return fe_lhs_simple_contact_material_ale;
372  };
373 
374  // add fields to the global matrix by adding the element
375  if (!eigen_pos_flag)
376  contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
377  "LAGMULT", contact_prisms);
378  else
379  contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
380  "LAGMULT", contact_prisms,
381  eigen_pos_flag, "EIGEN_POSITIONS");
382 
383  Range all_tets;
384  if (test_ale == PETSC_TRUE) {
385  if (!eigen_pos_flag)
386  contact_problem->addContactElementALE(
387  "ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
388  "LAGMULT", contact_prisms);
389  else
390  contact_problem->addContactElementALE(
391  "ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
392  "LAGMULT", contact_prisms, eigen_pos_flag, "EIGEN_POSITIONS");
393 
394  Range faces;
395  CHKERR moab.get_adjacencies(contact_prisms, 2, false, faces,
396  moab::Interface::UNION);
397  Range tris = faces.subset_by_type(MBTRI);
398 
399  CHKERR moab.get_adjacencies(tris, 3, false, all_tets,
400  moab::Interface::UNION);
401 
402  // Add finite elements
403  CHKERR m_field.add_finite_element("MATERIAL", MF_ZERO);
404  CHKERR m_field.modify_finite_element_add_field_row("MATERIAL",
405  "SPATIAL_POSITION");
406  CHKERR m_field.modify_finite_element_add_field_col("MATERIAL",
407  "SPATIAL_POSITION");
408  CHKERR m_field.modify_finite_element_add_field_row("MATERIAL",
409  "MESH_NODE_POSITIONS");
410  CHKERR m_field.modify_finite_element_add_field_col("MATERIAL",
411  "MESH_NODE_POSITIONS");
412  CHKERR m_field.modify_finite_element_add_field_data("MATERIAL",
413  "SPATIAL_POSITION");
414  CHKERR m_field.modify_finite_element_add_field_data(
415  "MATERIAL", "MESH_NODE_POSITIONS");
416  CHKERR m_field.add_ents_to_finite_element_by_type(all_tets, MBTET,
417  "MATERIAL");
418  CHKERR m_field.build_finite_elements("MATERIAL", &all_tets);
419  }
420 
421  // build finite elemnts
422  CHKERR m_field.build_finite_elements();
423 
424  // build adjacencies
425  CHKERR m_field.build_adjacencies(bit_levels.back());
426 
427  // define problems
428  CHKERR m_field.add_problem("CONTACT_PROB");
429 
430  // set refinement level for problem
431  CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
432  bit_levels.back());
433 
434  DMType dm_name = "DMMOFEM";
435  CHKERR DMRegister_MoFEM(dm_name);
436 
437  // create dm instance
438  SmartPetscObj<DM> dm;
439  dm = createSmartDM(m_field.get_comm(), dm_name);
440  CHKERR DMSetType(dm, dm_name);
441 
442  // set dm datastruture which created mofem datastructures
443  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_levels.back());
444  CHKERR DMSetFromOptions(dm);
445  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_FALSE);
446  // add elements to dm
447  CHKERR DMMoFEMAddElement(dm, "CONTACT_ELEM");
448 
449  if (test_ale == PETSC_TRUE) {
450  CHKERR DMMoFEMAddElement(dm, "ALE_CONTACT_ELEM");
451  CHKERR DMMoFEMAddElement(dm, "MATERIAL");
452  }
453 
454  CHKERR DMSetUp(dm);
455 
456  // Vector of DOFs and the RHS
457  auto D = smartCreateDMVector(dm);
458  auto F = smartVectorDuplicate(D);
459 
460  // Stiffness matrix
461  auto A = smartCreateDMMatrix(dm);
462 
463  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
464  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
465  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
466 
467  CHKERR VecZeroEntries(F);
468  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
469  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
470 
471  CHKERR MatSetOption(A, MAT_SPD, PETSC_TRUE);
472  CHKERR MatZeroEntries(A);
473 
474  auto fdA = smartMatDuplicate(A, MAT_COPY_VALUES);
475 
476  if (convect_pts == PETSC_TRUE) {
477  CHKERR DMMoFEMSNESSetFunction(
478  dm, "CONTACT_ELEM",
479  get_contact_rhs(contact_problem, make_convective_master_element),
480  PETSC_NULL, PETSC_NULL);
481  CHKERR DMMoFEMSNESSetJacobian(
482  dm, "CONTACT_ELEM",
483  get_master_contact_lhs(contact_problem,
484  make_convective_master_element),
485  NULL, NULL);
486  CHKERR DMMoFEMSNESSetFunction(
487  dm, "CONTACT_ELEM",
488  get_master_traction_rhs(contact_problem,
489  make_convective_slave_element),
490  PETSC_NULL, PETSC_NULL);
491  CHKERR DMMoFEMSNESSetJacobian(
492  dm, "CONTACT_ELEM",
493  get_master_traction_lhs(contact_problem,
494  make_convective_slave_element),
495  NULL, NULL);
496  } else {
497  CHKERR DMMoFEMSNESSetFunction(
498  dm, "CONTACT_ELEM",
499  get_contact_rhs(contact_problem, make_contact_element, alm_flag),
500  PETSC_NULL, PETSC_NULL);
501  CHKERR DMMoFEMSNESSetFunction(
502  dm, "CONTACT_ELEM",
503  get_master_traction_rhs(contact_problem, make_contact_element,
504  alm_flag),
505  PETSC_NULL, PETSC_NULL);
506  CHKERR DMMoFEMSNESSetJacobian(dm, "CONTACT_ELEM",
507  get_master_contact_lhs(contact_problem,
508  make_contact_element,
509  alm_flag),
510  PETSC_NULL, PETSC_NULL);
511  CHKERR DMMoFEMSNESSetJacobian(
512  dm, "CONTACT_ELEM",
513  get_master_traction_lhs(contact_problem, make_contact_element,
514  alm_flag),
515  PETSC_NULL, PETSC_NULL);
516  }
517 
518  if (test_ale == PETSC_TRUE) {
519  Range nodes;
520  CHKERR moab.get_connectivity(all_tets, nodes, false);
521 
522  CHKERR DMMoFEMSNESSetFunction(
523  dm, "ALE_CONTACT_ELEM",
524  get_contact_material_rhs(contact_problem, make_contact_element,
525  nodes),
526  PETSC_NULL, PETSC_NULL);
527 
528  CHKERR DMMoFEMSNESSetJacobian(
529  dm, "ALE_CONTACT_ELEM",
530  get_simple_contact_ale_lhs(contact_problem, make_contact_element),
531  NULL, NULL);
532 
533  CHKERR DMMoFEMSNESSetJacobian(
534  dm, "ALE_CONTACT_ELEM",
535  get_simple_contact_ale_material_lhs(contact_problem,
536  make_contact_element, nodes),
537  NULL, NULL);
538  }
539 
540  if (test_jacobian == PETSC_TRUE) {
541  char testing_options[] =
542  "-snes_test_jacobian -snes_test_jacobian_display "
543  "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 "
544  "-snes_max_it "
545  "1 ";
546  CHKERR PetscOptionsInsertString(NULL, testing_options);
547  } else {
548  char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
549  "-snes_rtol 0 "
550  "-snes_max_it 1 ";
551  CHKERR PetscOptionsInsertString(NULL, testing_options);
552  }
553 
554  auto snes = MoFEM::createSNES(m_field.get_comm());
555  SNESConvergedReason snes_reason;
556  SnesCtx *snes_ctx;
557 
558  // create snes nonlinear solver
559  {
560  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
561  CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
562  CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
563  CHKERR SNESSetFromOptions(snes);
564  }
565 
566  CHKERR SNESSolve(snes, PETSC_NULL, D);
567 
568  if (test_jacobian == PETSC_FALSE) {
569  double nrm_A0;
570  CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
571 
572  char testing_options_fd[] = "-snes_fd";
573  CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
574 
575  CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
576  CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
577  CHKERR SNESSetFromOptions(snes);
578 
579  CHKERR SNESSolve(snes, NULL, D);
580  CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
581 
582  double nrm_A;
583  CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
584  PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
585  nrm_A / nrm_A0);
586  nrm_A /= nrm_A0;
587 
588  constexpr double tol = 1e-6;
589  if (nrm_A > tol) {
590  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
591  "Difference between hand-calculated tangent matrix and finite "
592  "difference matrix is too big");
593  }
594  }
595  }
596  CATCH_ERRORS;
597 
598  // finish work cleaning memory, getting statistics, etc
599  MoFEM::Core::Finalize();
600 
601  return 0;
602 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
EntityHandle
SimpleContactProblem::ConvectSlaveContactElement
Element used to integrate on master surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:206
MoFEM::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition: PrismInterface.hpp:23
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
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:251
Mortar.hpp
MoFEM::createSmartDM
DEPRECATED auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Definition: PetscSmartObj.hpp:145
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
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:527
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::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:501
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::smartMatDuplicate
DEPRECATED SmartPetscObj< Mat > smartMatDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:240
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
main
int main(int argc, char *argv[])
Definition: test_jacobian_of_simple_contact_element.cpp:29
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
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
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
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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.
MoFEM::PrismInterface::splitSides
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
Definition: PrismInterface.cpp:519
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1098
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEM::smartCreateDMMatrix
DEPRECATED auto smartCreateDMMatrix(DM dm)
Definition: DMMoFEM.hpp:1011
MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition: PetscSmartObj.hpp:226
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:763
MoFEM::smartCreateDMVector
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1026
help
static char help[]
Definition: test_jacobian_of_simple_contact_element.cpp:28
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
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:118
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
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
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:98
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_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:1094
std
Definition: enable_if.hpp:5
MoFEM::PrismInterface::getSides
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
Definition: PrismInterface.cpp:56
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
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
SimpleContactProblem::ConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:181
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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_filed)=0
set finite element field data
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
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.
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
tol
double tol
Definition: mesh_smoothing.cpp:27
F
@ F
Definition: free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
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:722
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow