v0.13.1
arc_length_interface.cpp
Go to the documentation of this file.
1/** \file arc_length_interface.cpp
2 * \brief Example of arc-length with interface element
3
4 \todo Mak it work with multi-grid and distributed mesh
5
6 \todo Make more clever step adaptation
7
8*/
9
10
11
12#include <BasicFiniteElements.hpp>
13using namespace MoFEM;
14
16#include <Hooke.hpp>
18
19using namespace boost::numeric;
20
21static char help[] = "\
22 -my_file mesh file name\n\
23 -my_sr reduction of step size\n\
24 -my_its_d desired number of steps\n\
25 -my_ms maximal number of steps\n\n";
26
27#define DATAFILENAME "load_disp.txt"
28
29namespace CohesiveElement {
30
35 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
36 : ArcLengthIntElemFEMethod(m_field.get_moab(), arc_ptr), mField(m_field) {
37
38 for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(mField, "LoadPath", cit)) {
39 EntityHandle meshset = cit->getMeshset();
40 Range nodes;
41 rval = mOab.get_entities_by_type(meshset, MBVERTEX, nodes, true);
43 postProcNodes.merge(nodes);
44 }
45
46 PetscPrintf(PETSC_COMM_WORLD, "Nb. PostProcNodes %lu\n",
47 postProcNodes.size());
48 };
49
52 FILE *datafile;
53 PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "a+", &datafile);
54 const auto bit_number_lambda = mField.get_field_bit_number("LAMBDA");
55
56 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
57 problemPtr->getNumeredRowDofsPtr();
58 auto lit = numered_dofs_rows->lower_bound(
59 FieldEntity::getLoBitNumberUId(bit_number_lambda));
60 auto hi_lit = numered_dofs_rows->upper_bound(
61 FieldEntity::getHiBitNumberUId(bit_number_lambda));
62
63 if (lit == numered_dofs_rows->end()) {
64 fclose(datafile);
66 } else if (std::distance(lit, hi_lit) != 1) {
67 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
68 "Only one DOF is expected");
69 }
70
71 Range::iterator nit = postProcNodes.begin();
72 for (; nit != postProcNodes.end(); nit++) {
73 NumeredDofEntityByEnt::iterator dit, hi_dit;
74 dit = numered_dofs_rows->get<Ent_mi_tag>().lower_bound(*nit);
75 hi_dit = numered_dofs_rows->get<Ent_mi_tag>().upper_bound(*nit);
76 double coords[3];
77 CHKERR mOab.get_coords(&*nit, 1, coords);
78 for (; dit != hi_dit; dit++) {
79 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e -> ",
80 lit->get()->getName().c_str(), lit->get()->getDofCoeffIdx(),
81 lit->get()->getFieldData());
82 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e ",
83 dit->get()->getName().c_str(), dit->get()->getDofCoeffIdx(),
84 dit->get()->getFieldData());
85 PetscPrintf(PETSC_COMM_WORLD, "-> %3.4f %3.4f %3.4f\n", coords[0],
86 coords[1], coords[2]);
87 PetscFPrintf(PETSC_COMM_WORLD, datafile, "%6.4e %6.4e ",
88 dit->get()->getFieldData(), lit->get()->getFieldData());
89 }
90 }
91 PetscFPrintf(PETSC_COMM_WORLD, datafile, "\n");
92 fclose(datafile);
94 }
95};
96
98
101 boost::shared_ptr<ArcLengthCtx> arcPtr;
102
104 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
105 : mField(m_field), bodyForce(body_force), arcPtr(arc_ptr) {}
106
109
110 switch (snes_ctx) {
111 case CTX_SNESNONE: {
112 } break;
113 case CTX_SNESSETFUNCTION: {
114 CHKERR VecZeroEntries(snes_f);
115 CHKERR VecGhostUpdateBegin(snes_f, INSERT_VALUES, SCATTER_FORWARD);
116 CHKERR VecGhostUpdateEnd(snes_f, INSERT_VALUES, SCATTER_FORWARD);
117 } break;
118 default:
119 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
120 }
121
123 }
124
127 switch (snes_ctx) {
128 case CTX_SNESNONE: {
129 } break;
130 case CTX_SNESSETFUNCTION: {
131 CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
132 CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
133 CHKERR VecAssemblyBegin(snes_f);
134 CHKERR VecAssemblyEnd(snes_f);
135 // add F_lambda
136 CHKERR VecAXPY(snes_f, arcPtr->getFieldData(), arcPtr->F_lambda);
137 CHKERR VecAXPY(snes_f, -1., bodyForce);
138 PetscPrintf(PETSC_COMM_WORLD, "\tlambda = %6.4e\n",
139 arcPtr->getFieldData());
140 // snes_f norm
141 double fnorm;
142 CHKERR VecNorm(snes_f, NORM_2, &fnorm);
143 PetscPrintf(PETSC_COMM_WORLD, "\tfnorm = %6.4e\n", fnorm);
144 } break;
145 default:
146 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
147 }
148
150 }
151};
152
153} // namespace CohesiveElement
154
155using namespace CohesiveElement;
156
157int main(int argc, char *argv[]) {
158
159 const string default_options = "-ksp_type fgmres \n"
160 "-pc_type lu \n"
161 "-pc_factor_mat_solver_type mumps\n"
162 "-mat_mumps_icntl_20 0\n"
163 "-ksp_monitor \n"
164 "-ksp_atol 1e-10 \n"
165 "-ksp_rtol 1e-10 \n"
166 "-snes_monitor \n"
167 "-snes_type newtonls \n"
168 "-snes_linesearch_type l2 \n"
169 "-snes_linesearch_monitor \n"
170 "-snes_max_it 16 \n"
171 "-snes_atol 1e-8 \n"
172 "-snes_rtol 1e-8 \n"
173 "-snes_converged_reason \n";
174
175 string param_file = "param_file.petsc";
176 if (!static_cast<bool>(ifstream(param_file))) {
177 std::ofstream file(param_file.c_str(), std::ios::ate);
178 if (file.is_open()) {
179 file << default_options;
180 file.close();
181 }
182 }
183 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
184
185 try {
186
187 moab::Core mb_instance;
188 moab::Interface &moab = mb_instance;
189 int rank;
190 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
191
192 // Reade parameters from line command
193 PetscBool flg = PETSC_TRUE;
194 char mesh_file_name[255];
195 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
196 mesh_file_name, 255, &flg);
197 if (flg != PETSC_TRUE) {
198 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
199 "*** ERROR -my_file (MESH FILE NEEDED)");
200 }
201
202 PetscScalar step_size_reduction;
203 CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_sr",
204 &step_size_reduction, &flg);
205 if (flg != PETSC_TRUE) {
206 step_size_reduction = 1.;
207 }
208
209 PetscInt max_steps;
210 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_ms", &max_steps,
211 &flg);
212 if (flg != PETSC_TRUE) {
213 max_steps = 5;
214 }
215
216 int its_d;
217 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-my_its_d", &its_d, &flg);
218 if (flg != PETSC_TRUE) {
219 its_d = 6;
220 }
221
222 PetscInt order;
223 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
224 &flg);
225 if (flg != PETSC_TRUE) {
226 order = 2;
227 }
228
229 // Check if new start or restart. If new start, delete previous
230 // load_disp.txt
231 if (std::string(mesh_file_name).find("restart") == std::string::npos) {
232 remove(DATAFILENAME);
233 }
234
235 // Read mesh to MOAB
236 const char *option;
237 option = "";
238 CHKERR moab.load_file(mesh_file_name, 0, option);
239
240 // Data stored on mesh for restart
241 Tag th_step_size, th_step;
242 double def_step_size = 1;
243 rval = moab.tag_get_handle("_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
244 MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
245 if (rval == MB_ALREADY_ALLOCATED)
246 rval = MB_SUCCESS;
247 CHKERRG(rval);
248 int def_step = 1;
249 rval = moab.tag_get_handle("_STEP", 1, MB_TYPE_INTEGER, th_step,
250 MB_TAG_CREAT | MB_TAG_MESH, &def_step);
251 if (rval == MB_ALREADY_ALLOCATED)
252 rval = MB_SUCCESS;
253 CHKERRG(rval);
254 const void *tag_data_step_size[1];
255 EntityHandle root = 0;
256 CHKERR moab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
257 double &step_size = *(double *)tag_data_step_size[0];
258 const void *tag_data_step[1];
259 CHKERR moab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
260 int &step = *(int *)tag_data_step[0];
261 // end of data stored for restart
262 CHKERR PetscPrintf(PETSC_COMM_WORLD,
263 "Start step %D and step_size = %6.4e\n", step,
264 step_size);
265
266 // Create MoFEM 2database
267 MoFEM::Core core(moab);
268 MoFEM::Interface &m_field = core;
269
270 PrismInterface *interface_ptr;
271 CHKERR m_field.getInterface(interface_ptr);
272
273 Tag th_my_ref_level;
274 BitRefLevel def_bit_level = 0;
275 CHKERR m_field.get_moab().tag_get_handle(
276 "_MY_REFINEMENT_LEVEL", sizeof(BitRefLevel), MB_TYPE_OPAQUE,
277 th_my_ref_level, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_BYTES,
278 &def_bit_level);
279 const EntityHandle root_meshset = m_field.get_moab().get_root_set();
280 BitRefLevel *ptr_bit_level0;
281 CHKERR m_field.get_moab().tag_get_by_ptr(th_my_ref_level, &root_meshset, 1,
282 (const void **)&ptr_bit_level0);
283 BitRefLevel &bit_level0 = *ptr_bit_level0;
284 BitRefLevel problem_bit_level = bit_level0;
285
286 if (step == 1) {
287
288 // ref meshset ref level 0
289 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
290 0, 3, BitRefLevel().set(0));
291
292 std::vector<BitRefLevel> bit_levels;
293 bit_levels.push_back(BitRefLevel().set(0));
294
295 int ll = 1;
296
298 m_field, SIDESET | INTERFACESET, cit)) {
299 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert Interface %d\n",
300 cit->getMeshsetId());
301 EntityHandle cubit_meshset = cit->getMeshset();
302 {
303 // get tet entities form back bit_level
304 EntityHandle ref_level_meshset = 0;
305 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
307 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
308 BitRefLevel().set(), MBTET,
309 ref_level_meshset);
311 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
312 BitRefLevel().set(), MBPRISM,
313 ref_level_meshset);
314 Range ref_level_tets;
315 CHKERR moab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
316 true);
317 // get faces and test to split
318 CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
319 0);
320 // set new bit level
321 bit_levels.push_back(BitRefLevel().set(ll++));
322 // split faces and
323 CHKERR interface_ptr->splitSides(ref_level_meshset, bit_levels.back(),
324 cubit_meshset, true, true, 0);
325 // clean meshsets
326 CHKERR moab.delete_entities(&ref_level_meshset, 1);
327 }
328 // Update cubit meshsets
329 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
330 EntityHandle cubit_meshset = ciit->meshset;
332 ->updateMeshsetByEntitiesChildren(cubit_meshset,
333 bit_levels.back(),
334 cubit_meshset, MBMAXTYPE, true);
335 }
336 }
337
338 bit_level0 = bit_levels.back();
339 problem_bit_level = bit_level0;
340
341 /***/
342 // Define problem
343
344 // Fields
345 CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
346 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1,
348
349 CHKERR m_field.add_field("LAMBDA", NOFIELD, NOBASE, 1);
350
351 // Field for ArcLength
352 CHKERR
353 m_field.add_field("X0_DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
354
355 // FE
356 CHKERR m_field.add_finite_element("ELASTIC");
357
358 // Define rows/cols and element data
360 "DISPLACEMENT");
362 "DISPLACEMENT");
364 "DISPLACEMENT");
366 "ELASTIC", "MESH_NODE_POSITIONS");
367 CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
368 CHKERR m_field.modify_finite_element_add_field_col("ELASTIC", "LAMBDA");
369 // this is for paremtis
370 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "LAMBDA");
371
372 // FE Interface
373 CHKERR m_field.add_finite_element("INTERFACE");
375 "DISPLACEMENT");
377 "DISPLACEMENT");
379 "DISPLACEMENT");
381 "INTERFACE", "MESH_NODE_POSITIONS");
382
383 // FE ArcLength
384 CHKERR m_field.add_finite_element("ARC_LENGTH");
385
386 // Define rows/cols and element data
387 CHKERR m_field.modify_finite_element_add_field_row("ARC_LENGTH",
388 "LAMBDA");
389 CHKERR m_field.modify_finite_element_add_field_col("ARC_LENGTH",
390 "LAMBDA");
391
392 // elem data
393 CHKERR
394 m_field.modify_finite_element_add_field_data("ARC_LENGTH", "LAMBDA");
395
396 // define problems
397 CHKERR m_field.add_problem("ELASTIC_MECHANICS");
398
399 // set finite elements for problem
400 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
401 "ELASTIC");
402 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
403 "INTERFACE");
404 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
405 "ARC_LENGTH");
406 // set refinement level for problem
407 CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
408 problem_bit_level);
409
410 /***/
411 // Declare problem
412
413 // add entities (by tets) to the field
414 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
415 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "X0_DISPLACEMENT");
416 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
417
418 // add finite elements entities
420 problem_bit_level, BitRefLevel().set(), "ELASTIC", MBTET);
422 problem_bit_level, BitRefLevel().set(), "INTERFACE", MBPRISM);
423
424 // Setting up LAMBDA field and ARC_LENGTH interface
425 {
426 // Add dummy no-field vertex
427 EntityHandle no_field_vertex;
428 {
429 const double coords[] = {0, 0, 0};
430 CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
431 Range range_no_field_vertex;
432 range_no_field_vertex.insert(no_field_vertex);
433 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
434 range_no_field_vertex, BitRefLevel().set());
435
436 EntityHandle lambda_meshset = m_field.get_field_meshset("LAMBDA");
437 CHKERR m_field.get_moab().add_entities(lambda_meshset,
438 range_no_field_vertex);
439 }
440 // this entity will carray data for this finite element
441 EntityHandle meshset_fe_arc_length;
442 {
443 CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
444 CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
445 CHKERR m_field.getInterface<BitRefManager>()->setBitLevelToMeshset(
446 meshset_fe_arc_length, BitRefLevel().set());
447 }
448 // finally add created meshset to the ARC_LENGTH finite element
450 meshset_fe_arc_length, "ARC_LENGTH", false);
451 }
452
453 // set app. order
454 // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
455 // (Mark Ainsworth & Joe Coyle)
456 CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
457 CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
458 CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
459 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
460
461 CHKERR m_field.set_field_order(0, MBTET, "X0_DISPLACEMENT", order);
462 CHKERR m_field.set_field_order(0, MBTRI, "X0_DISPLACEMENT", order);
463 CHKERR m_field.set_field_order(0, MBEDGE, "X0_DISPLACEMENT", order);
464 CHKERR m_field.set_field_order(0, MBVERTEX, "X0_DISPLACEMENT", 1);
465
466 CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
467 CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
468 CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
469 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
470
471 // Elements with boundary conditions
472 CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
473 CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
474
475 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
476 "FORCE_FE");
477 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
478 "PRESSURE_FE");
479 CHKERR m_field.modify_finite_element_add_field_row("FORCE_FE", "LAMBDA");
480 CHKERR m_field.modify_finite_element_add_field_col("FORCE_FE", "LAMBDA");
481 CHKERR m_field.modify_finite_element_add_field_data("FORCE_FE", "LAMBDA");
482
483 CHKERR m_field.modify_finite_element_add_field_row("PRESSURE_FE",
484 "LAMBDA");
485 CHKERR m_field.modify_finite_element_add_field_col("PRESSURE_FE",
486 "LAMBDA");
487 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE_FE",
488 "LAMBDA");
489 CHKERR m_field.add_finite_element("BODY_FORCE");
490 CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
491 "DISPLACEMENT");
492 CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
493 "DISPLACEMENT");
494 CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
495 "DISPLACEMENT");
496 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
497 "BODY_FORCE");
498
500 m_field, BLOCKSET | BODYFORCESSET, it)) {
501 Range tets;
502 CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTET, tets,
503 true);
504 CHKERR m_field.add_ents_to_finite_element_by_type(tets, MBTET,
505 "BODY_FORCE");
506 }
507 }
508
509 /****/
510 // build database
511
512 // build field
513 CHKERR m_field.build_fields();
514 Projection10NodeCoordsOnField ent_method_material(m_field,
515 "MESH_NODE_POSITIONS");
516 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
517
518 // build finite elemnts
520
521 // build adjacencies
522 CHKERR m_field.build_adjacencies(problem_bit_level);
523
524 /****/
525 ProblemsManager *prb_mng_ptr;
526 CHKERR m_field.getInterface(prb_mng_ptr);
527 // build problem
528 CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
529 // partition
530 CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
531 CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS", false, 0,
532 m_field.get_comm_size());
533 // what are ghost nodes, see Petsc Manual
534 CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
535
536 // print bcs
537 MeshsetsManager *mmanager_ptr;
538 CHKERR m_field.getInterface(mmanager_ptr);
539 CHKERR mmanager_ptr->printDisplacementSet();
540 CHKERR mmanager_ptr->printForceSet();
541 // print block sets with materials
542 CHKERR mmanager_ptr->printMaterialsSet();
543
544 // create matrices
545 Vec F, F_body_force, D;
546 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
547 "ELASTIC_MECHANICS", COL, &F);
548 CHKERR VecDuplicate(F, &D);
549 CHKERR VecDuplicate(F, &F_body_force);
550 Mat Aij;
552 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
553 &Aij);
554
555 // Assemble F and Aij
556 double young_modulus = 1;
557 double poisson_ratio = 0.0;
558
559 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
560 interface_materials;
561
562 // FIXME this in fact allow only for one type of interface,
563 // problem is Young Modulus in interface mayoung_modulusterial
565 cout << std::endl << *it << std::endl;
566
567 // Get block name
568 string name = it->getName();
569
570 if (name.compare(0, 11, "MAT_ELASTIC") == 0) {
571 Mat_Elastic mydata;
572 CHKERR it->getAttributeDataStructure(mydata);
573 cout << mydata;
574 young_modulus = mydata.data.Young;
575 poisson_ratio = mydata.data.Poisson;
576 } else if (name.compare(0, 10, "MAT_INTERF") == 0) {
577 Mat_Interf mydata;
578 CHKERR it->getAttributeDataStructure(mydata);
579 cout << mydata;
580
581 interface_materials.push_back(
583 interface_materials.back().h = 1;
584 interface_materials.back().youngModulus = mydata.data.alpha;
585 interface_materials.back().beta = mydata.data.beta;
586 interface_materials.back().ft = mydata.data.ft;
587 interface_materials.back().Gf = mydata.data.Gf;
588
589 EntityHandle meshset = it->getMeshset();
590 Range tris;
591 rval = moab.get_entities_by_type(meshset, MBTRI, tris, true);
592 CHKERRG(rval);
593 Range ents3d;
594 rval = moab.get_adjacencies(tris, 3, false, ents3d,
595 moab::Interface::UNION);
596 CHKERRG(rval);
597 interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
598 }
599 }
600
601 { // FIXME
602 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator
603 pit = interface_materials.begin();
604 for (; pit != interface_materials.end(); pit++) {
605 pit->youngModulus = young_modulus;
606 }
607 }
608
609 boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
610 new ArcLengthCtx(m_field, "ELASTIC_MECHANICS"));
611 boost::scoped_ptr<ArcLengthElement> my_arc_method_ptr(
612 new ArcLengthElement(m_field, arc_ctx));
613 ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
614 AssembleRhsVectors pre_post_proc_fe(m_field, F_body_force, arc_ctx);
615
616 DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT", Aij, D, F);
617 CHKERR m_field.get_problem("ELASTIC_MECHANICS",
618 &my_dirichlet_bc.problemPtr);
619 CHKERR my_dirichlet_bc.iNitialize();
620 boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(new Hooke<adouble>);
621 boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>);
622 NonlinearElasticElement elastic(m_field, 2);
623 {
624 int id = 0;
625 elastic.setOfBlocks[id].iD = id;
626 elastic.setOfBlocks[id].E = young_modulus;
627 elastic.setOfBlocks[id].PoissonRatio = poisson_ratio;
628 elastic.setOfBlocks[id].materialDoublePtr = hooke_double_ptr;
629 elastic.setOfBlocks[id].materialAdoublePtr = hooke_adouble_ptr;
630 CHKERR m_field.get_moab().get_entities_by_type(
631 m_field.get_finite_element_meshset("ELASTIC"), MBTET,
632 elastic.setOfBlocks[id].tEts, true);
633 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeRhs(), true,
634 false, false, false);
635 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeRhs(), true,
636 false, false, false);
637 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeEnergy(), true,
638 false, false, false);
639 CHKERR elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false,
640 true);
641 }
642 CohesiveInterfaceElement cohesive_elements(m_field);
643 CHKERR cohesive_elements.addOps("DISPLACEMENT", interface_materials);
644
645 PetscInt M, N;
646 CHKERR MatGetSize(Aij, &M, &N);
647 PetscInt m, n;
648 MatGetLocalSize(Aij, &m, &n);
649 boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
650 new ArcLengthMatShell(Aij, arc_ctx, "ELASTIC_MECHANICS"));
651 Mat ShellAij;
652 CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, (void *)mat_ctx.get(),
653 &ShellAij);
654 CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
655 (void (*)(void))ArcLengthMatMultShellOp);
656
657 // body forces
658 BodyForceConstantField body_forces_methods(m_field);
660 m_field, BLOCKSET | BODYFORCESSET, it)) {
661 CHKERR body_forces_methods.addBlock("DISPLACEMENT", F_body_force,
662 it->getMeshsetId());
663 }
664 CHKERR VecZeroEntries(F_body_force);
665 CHKERR VecGhostUpdateBegin(F_body_force, INSERT_VALUES, SCATTER_FORWARD);
666 CHKERR VecGhostUpdateEnd(F_body_force, INSERT_VALUES, SCATTER_FORWARD);
667 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "BODY_FORCE",
668 body_forces_methods.getLoopFe());
669 CHKERR VecGhostUpdateBegin(F_body_force, ADD_VALUES, SCATTER_REVERSE);
670 CHKERR VecGhostUpdateEnd(F_body_force, ADD_VALUES, SCATTER_REVERSE);
671 CHKERR VecAssemblyBegin(F_body_force);
672 CHKERR VecAssemblyEnd(F_body_force);
673
674 // surface forces
675 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
676 string fe_name_str = "FORCE_FE";
677 neumann_forces.insert(fe_name_str, new NeumannForcesSurface(m_field));
678 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS",
679 neumann_forces.at(fe_name_str).getLoopFe(), false,
680 false);
682 it)) {
683 CHKERR neumann_forces.at(fe_name_str)
684 .addForce("DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
685 }
686 fe_name_str = "PRESSURE_FE";
687 neumann_forces.insert(fe_name_str, new NeumannForcesSurface(m_field));
688 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS",
689 neumann_forces.at(fe_name_str).getLoopFe(), false,
690 false);
692 m_field, SIDESET | PRESSURESET, it)) {
693 CHKERR neumann_forces.at(fe_name_str)
694 .addPressure("DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
695 }
696 // add nodal forces
697 boost::ptr_map<std::string, NodalForce> nodal_forces;
698 fe_name_str = "FORCE_FE";
699 nodal_forces.insert(fe_name_str, new NodalForce(m_field));
701 it)) {
702 CHKERR nodal_forces.at(fe_name_str)
703 .addForce("DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
704 }
705
706 SNES snes;
707 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
708 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
709 CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
710 CHKERR SNESSetJacobian(snes, ShellAij, Aij, SnesMat, &snes_ctx);
711 CHKERR SNESSetFromOptions(snes);
712
713 KSP ksp;
714 CHKERR SNESGetKSP(snes, &ksp);
715 PC pc;
716 CHKERR KSPGetPC(ksp, &pc);
717 boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
718 new PCArcLengthCtx(ShellAij, Aij, arc_ctx));
719 CHKERR PCSetType(pc, PCSHELL);
720 CHKERR PCShellSetContext(pc, pc_ctx.get());
721 CHKERR PCShellSetApply(pc, PCApplyArcLength);
722 CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
723
724 // Rhs
725 SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
726 snes_ctx.get_loops_to_do_Rhs();
727 snes_ctx.get_preProcess_to_do_Rhs().push_back(&my_dirichlet_bc);
728 snes_ctx.get_preProcess_to_do_Rhs().push_back(&pre_post_proc_fe);
729 loops_to_do_Rhs.push_back(SnesCtx::PairNameFEMethodPtr(
730 "INTERFACE", &cohesive_elements.getFeRhs()));
731 loops_to_do_Rhs.push_back(
732 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
733 loops_to_do_Rhs.push_back(
734 SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", my_arc_method_ptr.get()));
735 snes_ctx.get_postProcess_to_do_Rhs().push_back(&pre_post_proc_fe);
736 snes_ctx.get_postProcess_to_do_Rhs().push_back(&my_dirichlet_bc);
737
738 // Mat
739 SnesCtx::FEMethodsSequence &loops_to_do_Mat =
740 snes_ctx.get_loops_to_do_Mat();
741 snes_ctx.get_preProcess_to_do_Mat().push_back(&my_dirichlet_bc);
742 loops_to_do_Mat.push_back(SnesCtx::PairNameFEMethodPtr(
743 "INTERFACE", &cohesive_elements.getFeLhs()));
744 loops_to_do_Mat.push_back(
745 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
746 loops_to_do_Mat.push_back(
747 SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", my_arc_method_ptr.get()));
748 snes_ctx.get_postProcess_to_do_Mat().push_back(&my_dirichlet_bc);
749
750 double gamma = 0.5, reduction = 1;
751 // step = 1;
752 if (step == 1) {
753 step_size = step_size_reduction;
754 } else {
755 reduction = step_size_reduction;
756 step++;
757 }
758
759 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
760 neumann_forces.begin();
761 CHKERR VecZeroEntries(arc_ctx->F_lambda);
762 CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, INSERT_VALUES,
763 SCATTER_FORWARD);
764 CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
765 for (; mit != neumann_forces.end(); mit++) {
766 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", mit->first,
767 mit->second->getLoopFe());
768 }
769 CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
770 CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
771 CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
772 CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
773 for (std::vector<int>::iterator vit = my_dirichlet_bc.dofsIndices.begin();
774 vit != my_dirichlet_bc.dofsIndices.end(); vit++) {
775 CHKERR VecSetValue(arc_ctx->F_lambda, *vit, 0, INSERT_VALUES);
776 }
777 CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
778 CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
779 // F_lambda2
780 CHKERR VecDot(arc_ctx->F_lambda, arc_ctx->F_lambda, &arc_ctx->F_lambda2);
781 PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n", arc_ctx->F_lambda2);
782
783 if (step > 1) {
784 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
785 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
786 CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
787 "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL,
788 arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
789 double x0_nrm;
790 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
791 CHKERR PetscPrintf(PETSC_COMM_WORLD,
792 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
793 arc_ctx->dLambda);
794 CHKERR arc_ctx->setAlphaBeta(1, 0);
795 } else {
796 CHKERR arc_ctx->setS(0);
797 CHKERR arc_ctx->setAlphaBeta(0, 1);
798 }
799 CHKERR SnesRhs(snes, D, F, &snes_ctx);
800
801 PostProcVolumeOnRefinedMesh post_proc(m_field);
803 CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
804 CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
805 // add postpocessing for sresses
806 post_proc.getOpPtrVector().push_back(new PostProcHookStress(
807 m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "DISPLACEMENT",
808 post_proc.commonData));
809
810 bool converged_state = false;
811 for (; step < max_steps; step++) {
812
813 if (step == 1) {
814 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load Step %D step_size = %6.4e\n",
815 step, step_size);
816 CHKERR arc_ctx->setS(step_size);
817 CHKERR arc_ctx->setAlphaBeta(0, 1);
818 CHKERR VecCopy(D, arc_ctx->x0);
819 double dlambda;
820 CHKERR my_arc_method_ptr->calculate_init_dlambda(&dlambda);
821 CHKERR my_arc_method_ptr->set_dlambda_to_x(D, dlambda);
822 } else if (step == 2) {
823 CHKERR arc_ctx->setAlphaBeta(1, 0);
824 CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(D);
825 CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
826 CHKERR arc_ctx->setS(step_size);
827 double dlambda = arc_ctx->dLambda;
828 double dx_nrm;
829 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
830 CHKERR PetscPrintf(PETSC_COMM_WORLD,
831 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
832 "dx_nrm = %6.4e dx2 = %6.4e\n",
833 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
834 CHKERR VecCopy(D, arc_ctx->x0);
835 CHKERR VecAXPY(D, 1., arc_ctx->dx);
836 CHKERR my_arc_method_ptr->set_dlambda_to_x(D, dlambda);
837 } else {
838 CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(D);
839 CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
840 // step_size0_1/step_size0 = step_stize1/step_size
841 // step_size0_1 = step_size0*(step_stize1/step_size)
842 step_size *= reduction;
843 CHKERR arc_ctx->setS(step_size);
844 double dlambda = reduction * arc_ctx->dLambda;
845 CHKERR VecScale(arc_ctx->dx, reduction);
846 double dx_nrm;
847 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
848 CHKERR PetscPrintf(PETSC_COMM_WORLD,
849 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
850 "dx_nrm = %6.4e dx2 = %6.4e\n",
851 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
852 CHKERR VecCopy(D, arc_ctx->x0);
853 CHKERR VecAXPY(D, 1., arc_ctx->dx);
854 CHKERR my_arc_method_ptr->set_dlambda_to_x(D, dlambda);
855 }
856
857 CHKERR SNESSolve(snes, PETSC_NULL, D);
858
859 // Distribute displacements on all processors
860 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
861 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
862 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "INTERFACE",
863 cohesive_elements.getFeHistory(), 0,
864 m_field.get_comm_size());
865 // Remove nodes of damaged prisms
866 CHKERR my_arc_method_ptr->remove_damaged_prisms_nodes();
867
868 int its;
869 CHKERR SNESGetIterationNumber(snes, &its);
870 CHKERR PetscPrintf(PETSC_COMM_WORLD, "number of Newton iterations = %D\n",
871 its);
872
873 SNESConvergedReason reason;
874 CHKERR SNESGetConvergedReason(snes, &reason);
875
876 if (reason < 0) {
877 CHKERR arc_ctx->setAlphaBeta(1, 0);
878 reduction = 0.1;
879 converged_state = false;
880 continue;
881 } else {
882 if (step > 1 && converged_state) {
883 reduction = pow((double)its_d / (double)(its + 1), gamma);
884 CHKERR PetscPrintf(PETSC_COMM_WORLD, "reduction step_size = %6.4e\n",
885 reduction);
886 }
887
888 // Save data on mesh
889 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
890 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
891 CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
892 "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL,
893 arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
894 converged_state = true;
895 }
896 //
897 if (reason > 0) {
898 FILE *datafile;
899 PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "a+", &datafile);
900 PetscFPrintf(PETSC_COMM_WORLD, datafile, "%d %d ", reason, its);
901 fclose(datafile);
902 CHKERR my_arc_method_ptr->postProcessLoadPath();
903 }
904
905 if (step % 1 == 0) {
906
907 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
908 post_proc);
909 std::ostringstream ss;
910 ss << "out_values_" << step << ".h5m";
911 CHKERR post_proc.writeFile(ss.str().c_str());
912 }
913 }
914
915 // Save data on mesh
916 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
917 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
918
919 // detroy matrices
920 CHKERR VecDestroy(&F);
921 CHKERR VecDestroy(&D);
922 CHKERR VecDestroy(&F_body_force);
923 CHKERR MatDestroy(&Aij);
924 CHKERR SNESDestroy(&snes);
925 CHKERR MatDestroy(&ShellAij);
926 }
928
930
931 return 0;
932}
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
MoFEMErrorCode PCSetupArcLength(PC pc)
static Index< 'M', 3 > M
Implementation of linear interface element.
const std::string default_options
std::string param_file
Implementation of linear elastic material.
Implementation of arc-lebgth control for cohesive elements.
int main(int argc, char *argv[])
static char help[]
#define DATAFILENAME
@ COL
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ NOBASE
Definition: definitions.h:59
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ 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 ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ PRESSURESET
Definition: definitions.h:152
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
@ FORCESET
Definition: definitions.h:151
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
@ BLOCKSET
Definition: definitions.h:148
@ INTERFACESET
Definition: definitions.h:157
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
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 add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_bit_ref(const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name, EntityType type, int verb=DEFAULT_VERBOSITY)=0
add TET entities from given refinement level to finite element database given by name
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
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 add_ents_to_finite_element_by_MESHSET(const EntityHandle meshset, const std::string &name, const bool recursive=false)=0
add MESHSET element to finite element database given by name
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_col(const std::string &fe_name, const std::string &name_row)=0
set field col 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.
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
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.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#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.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
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
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:136
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
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
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
double young_modulus
Definition: plastic.cpp:87
const int N
Definition: speed_test.cpp:3
Store variables for ArcLength analysis.
shell matrix for arc-length method
Body forces elements.
Definition: BodyForce.hpp:12
MyVolumeFE & getLoopFe()
Definition: BodyForce.hpp:23
MoFEMErrorCode addBlock(const std::string field_name, Vec F, int ms_id)
Definition: BodyForce.hpp:86
ArcLengthElement(MoFEM::Interface &m_field, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
MoFEMErrorCode postProcess()
function is run at the end of loop
AssembleRhsVectors(MoFEM::Interface &m_field, Vec &body_force, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
boost::shared_ptr< ArcLengthCtx > arcPtr
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode addOps(const std::string field_name, boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation > &interfaces)
Driver function settting all operators needed for interface element.
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:71
virtual MoFEMErrorCode iNitialize()
Hook equation.
Definition: Hooke.hpp:21
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.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
const Problem * problemPtr
raw pointer to problem
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
structure for User Loop Methods on finite elements
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Elastic material data structure.
Linear interface data structure.
Matrix manager is used to build and partition problems.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
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...
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:29
Vec & snes_f
residual
SNESContext snes_ctx
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
Finite element and operators to apply force/pressures applied to surfaces.
Force applied to nodes.
Definition: NodalForce.hpp:13
structure grouping operators and data used for calculation of nonlinear elastic element
MyVolumeFE & getLoopFeLhs()
get lhs volume element
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE & getLoopFeRhs()
get rhs volume element
MyVolumeFE & getLoopFeEnergy()
get energy fe
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
structure for Arc Length pre-conditioner
Operator post-procesing stresses for Hook isotropic material.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
CommonData commonData
Post processing.