v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
test_jacobian_of_simple_contact_element.cpp File Reference
#include <Mortar.hpp>

Go to the source code of this file.

Functions

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

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

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

Definition at line 29 of file test_jacobian_of_simple_contact_element.cpp.

29 {
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) {
112 PrismInterface *interface;
113 CHKERR m_field.getInterface(interface);
114
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);
125 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
126 BitRefLevel().set(), MBTET,
127 ref_level_meshset);
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);
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
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) {
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;
234 };
235
236 auto set_pressure = [&](VectorAdaptor &&field_data, double *x, double *y,
237 double *z) {
239 double value;
240 double scale = 1.0;
241 PetscRandomGetValueReal(rctx, &value);
242 field_data[0] = value * scale;
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<
281 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
282 };
283
284 auto make_convective_slave_element = [&]() {
285 return boost::make_shared<
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);
405 "SPATIAL_POSITION");
407 "SPATIAL_POSITION");
409 "MESH_NODE_POSITIONS");
411 "MESH_NODE_POSITIONS");
413 "SPATIAL_POSITION");
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
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
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) {
478 dm, "CONTACT_ELEM",
479 get_contact_rhs(contact_problem, make_convective_master_element),
480 PETSC_NULL, PETSC_NULL);
482 dm, "CONTACT_ELEM",
483 get_master_contact_lhs(contact_problem,
484 make_convective_master_element),
485 NULL, NULL);
487 dm, "CONTACT_ELEM",
488 get_master_traction_rhs(contact_problem,
489 make_convective_slave_element),
490 PETSC_NULL, PETSC_NULL);
492 dm, "CONTACT_ELEM",
493 get_master_traction_lhs(contact_problem,
494 make_convective_slave_element),
495 NULL, NULL);
496 } else {
498 dm, "CONTACT_ELEM",
499 get_contact_rhs(contact_problem, make_contact_element, alm_flag),
500 PETSC_NULL, PETSC_NULL);
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);
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
523 dm, "ALE_CONTACT_ELEM",
524 get_contact_material_rhs(contact_problem, make_contact_element,
525 nodes),
526 PETSC_NULL, PETSC_NULL);
527
529 dm, "ALE_CONTACT_ELEM",
530 get_simple_contact_ale_lhs(contact_problem, make_contact_element),
531 NULL, NULL);
532
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 }
597
598 // finish work cleaning memory, getting statistics, etc
600
601 return 0;
602}
#define CATCH_ERRORS
Catch errors.
@ MF_ZERO
Definition definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
@ F
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
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:704
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
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:509
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:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1080
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
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.
#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.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
double D
double tol
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
DEPRECATED auto smartCreateDMMatrix(DM dm)
Definition DMMoFEM.hpp:996
auto createSNES(MPI_Comm comm)
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
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
DEPRECATED auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
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
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
DEPRECATED auto smartCreateDMVector(DM dm)
Definition DMMoFEM.hpp:1011
DEPRECATED SmartPetscObj< Mat > smartMatDuplicate(Mat mat, MatDuplicateOption op)
double scale
Definition plastic.cpp:170
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
Basic algebra on fields.
Definition FieldBlas.hpp:21
Create interface from given surface and insert flat prisms in-between.
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...
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...
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:13
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Element used to integrate on master surfaces. It convects integration points on slaves,...

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 28 of file test_jacobian_of_simple_contact_element.cpp.