v0.8.17
Classes | Namespaces | Macros | Functions | Variables
arc_length_interface.cpp File Reference

Example of arc-length with interface element. More...

#include <BasicFiniteElements.hpp>
#include <CohesiveInterfaceElement.hpp>
#include <Hooke.hpp>
#include <InterfaceGapArcLengthControl.hpp>

Go to the source code of this file.

Classes

struct  CohesiveElement::ArcLengthElement
 
struct  CohesiveElement::AssembleRhsVectors
 

Namespaces

 CohesiveElement
 

Macros

#define DATAFILENAME   "load_disp.txt"
 

Functions

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

Variables

static char help [] = "...\n\n"
 

Detailed Description

Example of arc-length with interface element.

Todo:
Mak it work with multi-grid and distributed mesh
Todo:
Make more clever step adaptation

Definition in file arc_length_interface.cpp.

Macro Definition Documentation

◆ DATAFILENAME

#define DATAFILENAME   "load_disp.txt"

Definition at line 35 of file arc_length_interface.cpp.

Function Documentation

◆ main()

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

Definition at line 156 of file arc_length_interface.cpp.

156  {
157 
158  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
159 
160  try {
161 
162  moab::Core mb_instance;
163  moab::Interface &moab = mb_instance;
164  int rank;
165  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
166 
167  // Reade parameters from line command
168  PetscBool flg = PETSC_TRUE;
169  char mesh_file_name[255];
170  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
171  mesh_file_name, 255, &flg);
172  if (flg != PETSC_TRUE) {
173  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
174  "*** ERROR -my_file (MESH FILE NEEDED)");
175  }
176 
177  PetscScalar step_size_reduction;
178  CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_sr",
179  &step_size_reduction, &flg);
180  if (flg != PETSC_TRUE) {
181  step_size_reduction = 1.;
182  }
183 
184  PetscInt max_steps;
185  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_ms", &max_steps,
186  &flg);
187  if (flg != PETSC_TRUE) {
188  max_steps = 5;
189  }
190 
191  int its_d;
192  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-my_its_d", &its_d, &flg);
193  if (flg != PETSC_TRUE) {
194  its_d = 6;
195  }
196 
197  PetscInt order;
198  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
199  &flg);
200  if (flg != PETSC_TRUE) {
201  order = 2;
202  }
203 
204  // Check if new start or restart. If new start, delete previous
205  // load_disp.txt
206  if (std::string(mesh_file_name).find("restart") == std::string::npos) {
207  remove(DATAFILENAME);
208  }
209 
210  // Read mesh to MOAB
211  const char *option;
212  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
213  CHKERR moab.load_file(mesh_file_name, 0, option);
214  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
215  if (pcomm == NULL)
216  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
217 
218  // Data stored on mesh for restart
219  Tag th_step_size, th_step;
220  double def_step_size = 1;
221  rval = moab.tag_get_handle("_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
222  MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
223  if (rval == MB_ALREADY_ALLOCATED)
224  rval = MB_SUCCESS;
225  CHKERRG(rval);
226  int def_step = 1;
227  rval = moab.tag_get_handle("_STEP", 1, MB_TYPE_INTEGER, th_step,
228  MB_TAG_CREAT | MB_TAG_MESH, &def_step);
229  if (rval == MB_ALREADY_ALLOCATED)
230  rval = MB_SUCCESS;
231  CHKERRG(rval);
232  const void *tag_data_step_size[1];
233  EntityHandle root = 0;
234  CHKERR moab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
235  double &step_size = *(double *)tag_data_step_size[0];
236  const void *tag_data_step[1];
237  CHKERR moab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
238  int &step = *(int *)tag_data_step[0];
239  // end of data stored for restart
240  CHKERR PetscPrintf(PETSC_COMM_WORLD,
241  "Start step %D and step_size = %6.4e\n", step,
242  step_size);
243 
244  // Create MoFEM 2database
245  MoFEM::Core core(moab);
246  MoFEM::Interface &m_field = core;
247 
248  PrismInterface *interface_ptr;
249  CHKERR m_field.getInterface(interface_ptr);
250 
251  Tag th_my_ref_level;
252  BitRefLevel def_bit_level = 0;
253  CHKERR m_field.get_moab().tag_get_handle(
254  "_MY_REFINEMENT_LEVEL", sizeof(BitRefLevel), MB_TYPE_OPAQUE,
255  th_my_ref_level, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_BYTES,
256  &def_bit_level);
257  const EntityHandle root_meshset = m_field.get_moab().get_root_set();
258  BitRefLevel *ptr_bit_level0;
259  CHKERR m_field.get_moab().tag_get_by_ptr(th_my_ref_level, &root_meshset, 1,
260  (const void **)&ptr_bit_level0);
261  BitRefLevel &bit_level0 = *ptr_bit_level0;
262  BitRefLevel problem_bit_level = bit_level0;
263 
264  if (step == 1) {
265 
266  // ref meshset ref level 0
267  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
268  0, 3, BitRefLevel().set(0));
269 
270  std::vector<BitRefLevel> bit_levels;
271  bit_levels.push_back(BitRefLevel().set(0));
272 
273  int ll = 1;
274 
276  m_field, SIDESET | INTERFACESET, cit)) {
277  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert Interface %d\n",
278  cit->getMeshsetId());
279  EntityHandle cubit_meshset = cit->getMeshset();
280  {
281  // get tet entities form back bit_level
282  EntityHandle ref_level_meshset = 0;
283  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
285  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
286  BitRefLevel().set(), MBTET,
287  ref_level_meshset);
289  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
290  BitRefLevel().set(), MBPRISM,
291  ref_level_meshset);
292  Range ref_level_tets;
293  CHKERR moab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
294  true);
295  // get faces and test to split
296  CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
297  0);
298  // set new bit level
299  bit_levels.push_back(BitRefLevel().set(ll++));
300  // split faces and
301  CHKERR interface_ptr->splitSides(ref_level_meshset, bit_levels.back(),
302  cubit_meshset, true, true, 0);
303  // clean meshsets
304  CHKERR moab.delete_entities(&ref_level_meshset, 1);
305  }
306  // Update cubit meshsets
307  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
308  EntityHandle cubit_meshset = ciit->meshset;
310  ->updateMeshsetByEntitiesChildren(cubit_meshset,
311  bit_levels.back(),
312  cubit_meshset, MBVERTEX, true);
314  ->updateMeshsetByEntitiesChildren(cubit_meshset,
315  bit_levels.back(),
316  cubit_meshset, MBEDGE, true);
318  ->updateMeshsetByEntitiesChildren(
319  cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
321  ->updateMeshsetByEntitiesChildren(
322  cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
323  }
324  }
325 
326  bit_level0 = bit_levels.back();
327  problem_bit_level = bit_level0;
328 
329  /***/
330  // Define problem
331 
332  // Fields
333  CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
334  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1,
336 
337  CHKERR m_field.add_field("LAMBDA", NOFIELD, NOBASE, 1);
338 
339  // Field for ArcLength
340  CHKERR
341  m_field.add_field("X0_DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
342 
343  // FE
344  CHKERR m_field.add_finite_element("ELASTIC");
345 
346  // Define rows/cols and element data
348  "DISPLACEMENT");
350  "DISPLACEMENT");
352  "DISPLACEMENT");
354  "ELASTIC", "MESH_NODE_POSITIONS");
355  CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
356  CHKERR m_field.modify_finite_element_add_field_col("ELASTIC", "LAMBDA");
357  // this is for paremtis
358  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "LAMBDA");
359 
360  // FE Interface
361  CHKERR m_field.add_finite_element("INTERFACE");
362  CHKERR m_field.modify_finite_element_add_field_row("INTERFACE",
363  "DISPLACEMENT");
364  CHKERR m_field.modify_finite_element_add_field_col("INTERFACE",
365  "DISPLACEMENT");
366  CHKERR m_field.modify_finite_element_add_field_data("INTERFACE",
367  "DISPLACEMENT");
369  "INTERFACE", "MESH_NODE_POSITIONS");
370 
371  // FE ArcLength
372  CHKERR m_field.add_finite_element("ARC_LENGTH");
373 
374  // Define rows/cols and element data
375  CHKERR m_field.modify_finite_element_add_field_row("ARC_LENGTH",
376  "LAMBDA");
377  CHKERR m_field.modify_finite_element_add_field_col("ARC_LENGTH",
378  "LAMBDA");
379 
380  // elem data
381  CHKERR
382  m_field.modify_finite_element_add_field_data("ARC_LENGTH", "LAMBDA");
383 
384  // define problems
385  CHKERR m_field.add_problem("ELASTIC_MECHANICS");
386 
387  // set finite elements for problem
388  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
389  "ELASTIC");
390  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
391  "INTERFACE");
392  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
393  "ARC_LENGTH");
394  // set refinement level for problem
395  CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
396  problem_bit_level);
397 
398  /***/
399  // Declare problem
400 
401  // add entities (by tets) to the field
402  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
403  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
404 
405  // add finite elements entities
407  problem_bit_level, BitRefLevel().set(), "ELASTIC", MBTET);
409  problem_bit_level, BitRefLevel().set(), "INTERFACE", MBPRISM);
410 
411  // Setting up LAMBDA field and ARC_LENGTH interface
412  {
413  // Add dummy no-field vertex
414  EntityHandle no_field_vertex;
415  {
416  const double coords[] = {0, 0, 0};
417  CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
418  CHKERRG(rval);
419  Range range_no_field_vertex;
420  range_no_field_vertex.insert(no_field_vertex);
421  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
422  range_no_field_vertex, BitRefLevel().set());
423 
424  EntityHandle lambda_meshset = m_field.get_field_meshset("LAMBDA");
425  CHKERR m_field.get_moab().add_entities(lambda_meshset,
426  range_no_field_vertex);
427  CHKERRG(rval);
428  }
429  // this entity will carray data for this finite element
430  EntityHandle meshset_fe_arc_length;
431  {
432  CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
433  CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
434  CHKERR m_field.getInterface<BitRefManager>()->setBitLevelToMeshset(
435  meshset_fe_arc_length, BitRefLevel().set());
436  }
437  // finally add created meshset to the ARC_LENGTH finite element
439  meshset_fe_arc_length, "ARC_LENGTH", false);
440  }
441 
442  // set app. order
443  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
444  // (Mark Ainsworth & Joe Coyle)
445  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
446  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
447  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
448  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
449 
450  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
451  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
452  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
453  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
454 
455  // Elements with boundary conditions
456  CHKERR MetaNeummanForces::addNeumannBCElements(m_field, "DISPLACEMENT");
457  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
458 
459  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
460  "FORCE_FE");
461  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
462  "PRESSURE_FE");
463  CHKERR m_field.modify_finite_element_add_field_row("FORCE_FE", "LAMBDA");
464  CHKERR m_field.modify_finite_element_add_field_col("FORCE_FE", "LAMBDA");
465  CHKERR m_field.modify_finite_element_add_field_data("FORCE_FE", "LAMBDA");
466 
467  CHKERR m_field.modify_finite_element_add_field_row("PRESSURE_FE",
468  "LAMBDA");
469  CHKERR m_field.modify_finite_element_add_field_col("PRESSURE_FE",
470  "LAMBDA");
471  CHKERR m_field.modify_finite_element_add_field_data("PRESSURE_FE",
472  "LAMBDA");
473  CHKERR m_field.add_finite_element("BODY_FORCE");
474  CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
475  "DISPLACEMENT");
476  CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
477  "DISPLACEMENT");
478  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
479  "DISPLACEMENT");
480  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
481  "BODY_FORCE");
482 
484  m_field, BLOCKSET | BODYFORCESSET, it)) {
485  Range tets;
486  CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTET, tets,
487  true);
488  CHKERR m_field.add_ents_to_finite_element_by_type(tets, MBTET,
489  "BODY_FORCE");
490  }
491  }
492 
493  /****/
494  // build database
495 
496  // build field
497  CHKERR m_field.build_fields();
498  Projection10NodeCoordsOnField ent_method_material(m_field,
499  "MESH_NODE_POSITIONS");
500  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
501 
502  // build finite elemnts
503  CHKERR m_field.build_finite_elements();
504 
505  // build adjacencies
506  CHKERR m_field.build_adjacencies(problem_bit_level);
507 
508  /****/
509  ProblemsManager *prb_mng_ptr;
510  CHKERR m_field.getInterface(prb_mng_ptr);
511  // build problem
512  CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
513  // partition
514  CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
515  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS", false, 0,
516  pcomm->size());
517  // what are ghost nodes, see Petsc Manual
518  CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
519 
520  // print bcs
521  MeshsetsManager *mmanager_ptr;
522  CHKERR m_field.getInterface(mmanager_ptr);
523  CHKERR mmanager_ptr->printDisplacementSet();
524  CHKERR mmanager_ptr->printForceSet();
525  // print block sets with materials
526  CHKERR mmanager_ptr->printMaterialsSet();
527 
528  // create matrices
529  Vec F, F_body_force, D;
530  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
531  "ELASTIC_MECHANICS", COL, &F);
532  CHKERR VecDuplicate(F, &D);
533  CHKERR VecDuplicate(F, &F_body_force);
534  Mat Aij;
535  CHKERR m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS", &Aij);
536 
537  // Assemble F and Aij
538  double young_modulus = 1;
539  double poisson_ratio = 0.0;
540 
541  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
542  interface_materials;
543 
544  // FIXME this in fact allow only for one type of interface,
545  // problem is Young Modulus in interface mayoung_modulusterial
547  cout << std::endl << *it << std::endl;
548 
549  // Get block name
550  string name = it->getName();
551 
552  if (name.compare(0, 11, "MAT_ELASTIC") == 0) {
553  Mat_Elastic mydata;
554  CHKERR it->getAttributeDataStructure(mydata);
555  cout << mydata;
556  young_modulus = mydata.data.Young;
557  poisson_ratio = mydata.data.Poisson;
558  } else if (name.compare(0, 10, "MAT_INTERF") == 0) {
559  Mat_Interf mydata;
560  CHKERR it->getAttributeDataStructure(mydata);
561  cout << mydata;
562 
563  interface_materials.push_back(
565  interface_materials.back().h = 1;
566  interface_materials.back().youngModulus = mydata.data.alpha;
567  interface_materials.back().beta = mydata.data.beta;
568  interface_materials.back().ft = mydata.data.ft;
569  interface_materials.back().Gf = mydata.data.Gf;
570 
571  EntityHandle meshset = it->getMeshset();
572  Range tris;
573  rval = moab.get_entities_by_type(meshset, MBTRI, tris, true);
574  CHKERRG(rval);
575  Range ents3d;
576  rval = moab.get_adjacencies(tris, 3, false, ents3d,
577  moab::Interface::UNION);
578  CHKERRG(rval);
579  interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
580  }
581  }
582 
583  { // FIXME
584  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator
585  pit = interface_materials.begin();
586  for (; pit != interface_materials.end(); pit++) {
587  pit->youngModulus = young_modulus;
588  }
589  }
590 
591  boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
592  new ArcLengthCtx(m_field, "ELASTIC_MECHANICS"));
593  boost::scoped_ptr<ArcLengthElement> my_arc_method_ptr(
594  new ArcLengthElement(m_field, arc_ctx));
595  ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
596  AssembleRhsVectors pre_post_proc_fe(m_field, F_body_force, arc_ctx);
597 
598  DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT", Aij, D, F);
599  CHKERR m_field.get_problem("ELASTIC_MECHANICS",
600  &my_dirichlet_bc.problemPtr);
601  CHKERR my_dirichlet_bc.iNitalize();
602  boost::shared_ptr<Hooke<adouble> > hooke_adouble_ptr(new Hooke<adouble>);
603  boost::shared_ptr<Hooke<double> > hooke_double_ptr(new Hooke<double>);
604  NonlinearElasticElement elastic(m_field, 2);
605  {
606  int id = 0;
607  elastic.setOfBlocks[id].iD = id;
608  elastic.setOfBlocks[id].E = young_modulus;
609  elastic.setOfBlocks[id].PoissonRatio = poisson_ratio;
610  elastic.setOfBlocks[id].materialDoublePtr = hooke_double_ptr;
611  elastic.setOfBlocks[id].materialAdoublePtr = hooke_adouble_ptr;
612  rval = m_field.get_moab().get_entities_by_type(
613  m_field.get_finite_element_meshset("ELASTIC"), MBTET,
614  elastic.setOfBlocks[id].tEts, true);
615  CHKERRG(rval);
616  CHKERR elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false,
617  true);
618  }
619  CohesiveInterfaceElement cohesive_elements(m_field);
620  CHKERR cohesive_elements.addOps("DISPLACEMENT", interface_materials);
621 
622  PetscInt M, N;
623  CHKERR MatGetSize(Aij, &M, &N);
624  PetscInt m, n;
625  MatGetLocalSize(Aij, &m, &n);
626  boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
627  new ArcLengthMatShell(Aij, arc_ctx, "ELASTIC_MECHANICS"));
628  Mat ShellAij;
629  CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, (void *)mat_ctx.get(),
630  &ShellAij);
631  CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
632  (void (*)(void))ArcLengthMatMultShellOp);
633 
634  // body forces
635  BodyForceConstantField body_forces_methods(m_field);
637  m_field, BLOCKSET | BODYFORCESSET, it)) {
638  CHKERR body_forces_methods.addBlock("DISPLACEMENT", F_body_force,
639  it->getMeshsetId());
640  }
641  CHKERR VecZeroEntries(F_body_force);
642  CHKERR VecGhostUpdateBegin(F_body_force, INSERT_VALUES, SCATTER_FORWARD);
643  CHKERR VecGhostUpdateEnd(F_body_force, INSERT_VALUES, SCATTER_FORWARD);
644  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "BODY_FORCE",
645  body_forces_methods.getLoopFe());
646  CHKERR VecGhostUpdateBegin(F_body_force, ADD_VALUES, SCATTER_REVERSE);
647  CHKERR VecGhostUpdateEnd(F_body_force, ADD_VALUES, SCATTER_REVERSE);
648  CHKERR VecAssemblyBegin(F_body_force);
649  CHKERR VecAssemblyEnd(F_body_force);
650 
651  // surface forces
652  boost::ptr_map<std::string, NeummanForcesSurface> neumann_forces;
653  string fe_name_str = "FORCE_FE";
654  neumann_forces.insert(fe_name_str, new NeummanForcesSurface(m_field));
656  it)) {
657  CHKERR neumann_forces.at(fe_name_str)
658  .addForce("DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
659  }
660  fe_name_str = "PRESSURE_FE";
661  neumann_forces.insert(fe_name_str, new NeummanForcesSurface(m_field));
663  m_field, SIDESET | PRESSURESET, it)) {
664  CHKERR neumann_forces.at(fe_name_str)
665  .addPressure("DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
666  }
667  // add nodal forces
668  boost::ptr_map<std::string, NodalForce> nodal_forces;
669  fe_name_str = "FORCE_FE";
670  nodal_forces.insert(fe_name_str, new NodalForce(m_field));
672  it)) {
673  CHKERR nodal_forces.at(fe_name_str)
674  .addForce("DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
675  }
676 
677  SNES snes;
678  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
679  CHKERR SNESSetApplicationContext(snes, &snes_ctx);
680  CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
681  CHKERR SNESSetJacobian(snes, ShellAij, Aij, SnesMat, &snes_ctx);
682  CHKERR SNESSetFromOptions(snes);
683 
684  KSP ksp;
685  CHKERR SNESGetKSP(snes, &ksp);
686  PC pc;
687  CHKERR KSPGetPC(ksp, &pc);
688  boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
689  new PCArcLengthCtx(ShellAij, Aij, arc_ctx));
690  CHKERR PCSetType(pc, PCSHELL);
691  CHKERR PCShellSetContext(pc, pc_ctx.get());
692  CHKERR PCShellSetApply(pc, PCApplyArcLength);
693  CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
694 
695  // Rhs
696  SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
697  snes_ctx.get_loops_to_do_Rhs();
698  snes_ctx.get_preProcess_to_do_Rhs().push_back(&my_dirichlet_bc);
699  snes_ctx.get_preProcess_to_do_Rhs().push_back(&pre_post_proc_fe);
700  loops_to_do_Rhs.push_back(SnesCtx::PairNameFEMethodPtr(
701  "INTERFACE", &cohesive_elements.getFeRhs()));
702  loops_to_do_Rhs.push_back(
703  SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
704  loops_to_do_Rhs.push_back(
705  SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", my_arc_method_ptr.get()));
706  snes_ctx.get_postProcess_to_do_Rhs().push_back(&pre_post_proc_fe);
707  snes_ctx.get_postProcess_to_do_Rhs().push_back(&my_dirichlet_bc);
708 
709  // Mat
710  SnesCtx::FEMethodsSequence &loops_to_do_Mat =
711  snes_ctx.get_loops_to_do_Mat();
712  snes_ctx.get_preProcess_to_do_Mat().push_back(&my_dirichlet_bc);
713  loops_to_do_Mat.push_back(SnesCtx::PairNameFEMethodPtr(
714  "INTERFACE", &cohesive_elements.getFeLhs()));
715  loops_to_do_Mat.push_back(
716  SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
717  loops_to_do_Mat.push_back(
718  SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", my_arc_method_ptr.get()));
719  snes_ctx.get_postProcess_to_do_Mat().push_back(&my_dirichlet_bc);
720 
721  double gamma = 0.5, reduction = 1;
722  // step = 1;
723  if (step == 1) {
724  step_size = step_size_reduction;
725  } else {
726  reduction = step_size_reduction;
727  step++;
728  }
729 
730  boost::ptr_map<std::string, NeummanForcesSurface>::iterator mit =
731  neumann_forces.begin();
732  CHKERR VecZeroEntries(arc_ctx->F_lambda);
733  CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, INSERT_VALUES,
734  SCATTER_FORWARD);
735  CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
736  for (; mit != neumann_forces.end(); mit++) {
737  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", mit->first,
738  mit->second->getLoopFe());
739  }
740  CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
741  CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
742  CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
743  CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
744  for (std::vector<int>::iterator vit = my_dirichlet_bc.dofsIndices.begin();
745  vit != my_dirichlet_bc.dofsIndices.end(); vit++) {
746  CHKERR VecSetValue(arc_ctx->F_lambda, *vit, 0, INSERT_VALUES);
747  }
748  CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
749  CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
750  // F_lambda2
751  CHKERR VecDot(arc_ctx->F_lambda, arc_ctx->F_lambda, &arc_ctx->F_lambda2);
752  PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n", arc_ctx->F_lambda2);
753 
754  if (step > 1) {
755  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
756  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
757  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
758  "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL,
759  arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
760  double x0_nrm;
761  CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
762  CHKERR PetscPrintf(PETSC_COMM_WORLD,
763  "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
764  arc_ctx->dLambda);
765  CHKERR arc_ctx->setAlphaBeta(1, 0);
766  } else {
767  CHKERR arc_ctx->setS(0);
768  CHKERR arc_ctx->setAlphaBeta(0, 1);
769  }
770  CHKERR SnesRhs(snes, D, F, &snes_ctx);
771 
772  PostProcVolumeOnRefinedMesh post_proc(m_field);
773  CHKERR post_proc.generateReferenceElementMesh();
774  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
775  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
776  // add postpocessing for sresses
777  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
778  m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "DISPLACEMENT",
779  post_proc.commonData));
780 
781  bool converged_state = false;
782  for (; step < max_steps; step++) {
783 
784  if (step == 1) {
785  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load Step %D step_size = %6.4e\n",
786  step, step_size);
787  CHKERR arc_ctx->setS(step_size);
788  CHKERR arc_ctx->setAlphaBeta(0, 1);
789  CHKERR VecCopy(D, arc_ctx->x0);
790  double dlambda;
791  CHKERR my_arc_method_ptr->calculate_init_dlambda(&dlambda);
792  CHKERR my_arc_method_ptr->set_dlambda_to_x(D, dlambda);
793  } else if (step == 2) {
794  CHKERR arc_ctx->setAlphaBeta(1, 0);
795  CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(D);
796  CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
797  CHKERR arc_ctx->setS(step_size);
798  double dlambda = arc_ctx->dLambda;
799  double dx_nrm;
800  CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
801  CHKERR PetscPrintf(PETSC_COMM_WORLD,
802  "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
803  "dx_nrm = %6.4e dx2 = %6.4e\n",
804  step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
805  CHKERR VecCopy(D, arc_ctx->x0);
806  CHKERR VecAXPY(D, 1., arc_ctx->dx);
807  CHKERR my_arc_method_ptr->set_dlambda_to_x(D, dlambda);
808  } else {
809  CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(D);
810  CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
811  // step_size0_1/step_size0 = step_stize1/step_size
812  // step_size0_1 = step_size0*(step_stize1/step_size)
813  step_size *= reduction;
814  CHKERR arc_ctx->setS(step_size);
815  double dlambda = reduction * arc_ctx->dLambda;
816  CHKERR VecScale(arc_ctx->dx, reduction);
817  double dx_nrm;
818  CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
819  CHKERR PetscPrintf(PETSC_COMM_WORLD,
820  "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
821  "dx_nrm = %6.4e dx2 = %6.4e\n",
822  step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
823  CHKERR VecCopy(D, arc_ctx->x0);
824  CHKERR VecAXPY(D, 1., arc_ctx->dx);
825  CHKERR my_arc_method_ptr->set_dlambda_to_x(D, dlambda);
826  }
827 
828  CHKERR SNESSolve(snes, PETSC_NULL, D);
829 
830  // Distribute displacements on all processors
831  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
832  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
833  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "INTERFACE",
834  cohesive_elements.getFeHistory(), 0,
835  m_field.get_comm_size());
836  // Remove nodes of damaged prisms
837  CHKERR my_arc_method_ptr->remove_damaged_prisms_nodes();
838 
839  int its;
840  CHKERR SNESGetIterationNumber(snes, &its);
841  CHKERR PetscPrintf(PETSC_COMM_WORLD, "number of Newton iterations = %D\n",
842  its);
843 
844  SNESConvergedReason reason;
845  CHKERR SNESGetConvergedReason(snes, &reason);
846 
847  if (reason < 0) {
848  CHKERR arc_ctx->setAlphaBeta(1, 0);
849  reduction = 0.1;
850  converged_state = false;
851  continue;
852  } else {
853  if (step > 1 && converged_state) {
854  reduction = pow((double)its_d / (double)(its + 1), gamma);
855  CHKERR PetscPrintf(PETSC_COMM_WORLD, "reduction step_size = %6.4e\n",
856  reduction);
857  }
858 
859  // Save data on mesh
860  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
861  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
862  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
863  "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL,
864  arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
865  converged_state = true;
866  }
867  //
868  if (reason > 0) {
869  FILE *datafile;
870  PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "a+", &datafile);
871  PetscFPrintf(PETSC_COMM_WORLD, datafile, "%d %d ", reason, its);
872  fclose(datafile);
873  CHKERR my_arc_method_ptr->postProcessLoadPath();
874  }
875 
876  if (step % 1 == 0) {
877 
878  // if(pcomm->rank()==0) {
879  // std::ostringstream sss;
880  // sss << "restart_" << step << ".h5m";
881  // rval = moab.write_file(sss.str().c_str()); CHKERRG(rval);
882  // }
883 
884  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
885  post_proc);
886  std::ostringstream ss;
887  ss << "out_values_" << step << ".h5m";
888  CHKERR post_proc.writeFile(ss.str().c_str());
889  }
890  }
891 
892  // Save data on mesh
893  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
894  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
895 
896  // detroy matrices
897  CHKERR VecDestroy(&F);
898  CHKERR VecDestroy(&D);
899  CHKERR VecDestroy(&F_body_force);
900  CHKERR MatDestroy(&Aij);
901  CHKERR SNESDestroy(&snes);
902  CHKERR MatDestroy(&ShellAij);
903  }
904  CATCH_ERRORS;
905 
907 
908  return 0;
909 }
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 modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
#define DATAFILENAME
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
Problem manager is used to build and partition problems .
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode buildProblem(const std::string &name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionGhostDofs(const std::string &name, int verb=VERBOSE)
determine ghost nodes
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.
virtual moab::Interface & get_moab()=0
Hook equation.
Definition: Hooke.hpp:31
MoFEMErrorCode partitionFiniteElements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elementsFunction which partition finite elements based on dofs partitioning. In addition it sets information about local row and cols dofs at given element on partition.
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
MoFEMErrorCode partitionProblem(const std::string &name, int verb=VERBOSE)
partition problem dofs (collective)
scalar or vector of scalars describe (no true field)
Definition: definitions.h:167
virtual int get_comm_size() const =0
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:128
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 ...
Elastic material data structure.
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 tetrahedral in children sets and add prism elements ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:534
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.
Interface for managing meshsets containing materials and boundary conditions.
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.
Core (interface) class.
Definition: Core.hpp:50
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
char mesh_file_name[255]
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Projection of edge entities with one mid-node on hierarchical basis.
shell matrix for arc-length method
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
Operator post-procesing stresses for Hook isotropic material.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
virtual EntityHandle get_field_meshset(const std::string &name) const =0
get field meshset
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
Body forces elements.
Definition: BodyForce.hpp:25
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:43
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:142
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Linear interface data structure.
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:104
Make interface on given faces and create flat prism in that space.
Post processing.
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:50
block name is "BODY_FORCES"
Definition: definitions.h:222
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 modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Vector manager is used to create vectors .
Definition: VecManager.hpp:35
structure for Arc Length pre-conditioner
Force applied to nodes.
Definition: NodalForce.hpp:25
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
create two children meshsets in the meshset containing tetrahedral on two sides of faces ...
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#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...
#define CHKERR
Inline error check.
Definition: definitions.h:586
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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 MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)=0
create Mat (MPIAIJ) for problem (collective)
static char help[]
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:282
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEMErrorCode PCSetupArcLength(PC pc)
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
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
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:430
continuous field
Definition: definitions.h:168
Constitutive (physical) equation for interface.
Finite element and operators to apply force/pressures applied to surfaces.
const int N
Definition: speed_test.cpp:3
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:55
Store variables for ArcLength analysis.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 33 of file arc_length_interface.cpp.