v0.14.0
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>
13 using namespace MoFEM;
14 
15 #include <CohesiveInterfaceElement.hpp>
16 #include <Hooke.hpp>
17 #include <InterfaceGapArcLengthControl.hpp>
18 
19 using namespace boost::numeric;
20 
21 static 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 
29 namespace CohesiveElement {
30 
31 struct ArcLengthElement : public ArcLengthIntElemFEMethod {
32  MoFEM::Interface &mField;
33  Range postProcNodes;
34  ArcLengthElement(MoFEM::Interface &m_field,
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);
42  MOAB_THROW(rval);
43  postProcNodes.merge(nodes);
44  }
45 
46  PetscPrintf(PETSC_COMM_WORLD, "Nb. PostProcNodes %lu\n",
47  postProcNodes.size());
48  };
49 
50  MoFEMErrorCode postProcessLoadPath() {
51  MoFEMFunctionBegin;
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);
65  MoFEMFunctionReturnHot(0);
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);
93  MoFEMFunctionReturn(0);
94  }
95 };
96 
97 struct AssembleRhsVectors : public FEMethod {
98 
99  MoFEM::Interface &mField;
100  Vec &bodyForce;
101  boost::shared_ptr<ArcLengthCtx> arcPtr;
102 
103  AssembleRhsVectors(MoFEM::Interface &m_field, Vec &body_force,
104  boost::shared_ptr<ArcLengthCtx> &arc_ptr)
105  : mField(m_field), bodyForce(body_force), arcPtr(arc_ptr) {}
106 
107  MoFEMErrorCode preProcess() {
108  MoFEMFunctionBegin;
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 
122  MoFEMFunctionReturn(0);
123  }
124 
125  MoFEMErrorCode postProcess() {
126  MoFEMFunctionBegin;
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 
149  MoFEMFunctionReturn(0);
150  }
151 };
152 
153 } // namespace CohesiveElement
154 
155 using namespace CohesiveElement;
156 
157 int 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 
297  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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);
306  CHKERR m_field.getInterface<BitRefManager>()
307  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
308  BitRefLevel().set(), MBTET,
309  ref_level_meshset);
310  CHKERR m_field.getInterface<BitRefManager>()
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;
331  CHKERR m_field.getInterface<BitRefManager>()
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,
347  AINSWORTH_LEGENDRE_BASE, 3);
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
359  CHKERR m_field.modify_finite_element_add_field_row("ELASTIC",
360  "DISPLACEMENT");
361  CHKERR m_field.modify_finite_element_add_field_col("ELASTIC",
362  "DISPLACEMENT");
363  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
364  "DISPLACEMENT");
365  CHKERR m_field.modify_finite_element_add_field_data(
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");
374  CHKERR m_field.modify_finite_element_add_field_row("INTERFACE",
375  "DISPLACEMENT");
376  CHKERR m_field.modify_finite_element_add_field_col("INTERFACE",
377  "DISPLACEMENT");
378  CHKERR m_field.modify_finite_element_add_field_data("INTERFACE",
379  "DISPLACEMENT");
380  CHKERR m_field.modify_finite_element_add_field_data(
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
419  CHKERR m_field.add_ents_to_finite_element_by_bit_ref(
420  problem_bit_level, BitRefLevel().set(), "ELASTIC", MBTET);
421  CHKERR m_field.add_ents_to_finite_element_by_bit_ref(
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
449  CHKERR m_field.add_ents_to_finite_element_by_MESHSET(
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 
499  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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
519  CHKERR m_field.build_finite_elements();
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;
551  CHKERR m_field.getInterface<MatrixManager>()
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
564  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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(
582  new CohesiveInterfaceElement::PhysicalEquation(m_field));
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);
659  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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);
681  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
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);
691  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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));
700  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
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.getComputeRhs();
727  snes_ctx.getPreProcComputeRhs().push_back(&my_dirichlet_bc);
728  snes_ctx.getPreProcComputeRhs().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.getPostProcComputeRhs().push_back(&pre_post_proc_fe);
736  snes_ctx.getPostProcComputeRhs().push_back(&my_dirichlet_bc);
737 
738  // Mat
739  SnesCtx::FEMethodsSequence &loops_to_do_Mat =
740  snes_ctx.getSetOperators();
741  snes_ctx.getPreProcSetOperators().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.getPostProcSetOperators().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);
802  CHKERR post_proc.generateReferenceElementMesh();
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  }
927  CATCH_ERRORS;
928 
929  MoFEM::Core::Finalize();
930 
931  return 0;
932 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:699
MoFEM::Ent_mi_tag
Definition: TagMultiIndices.hpp:21
SIDESET
@ SIDESET
Definition: definitions.h:160
main
int main(int argc, char *argv[])
Definition: arc_length_interface.cpp:157
MoFEM::CoreInterface::loop_finite_elements
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.
MoFEM::CoreInterface::add_ents_to_finite_element_by_MESHSET
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
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
CohesiveInterfaceElement.hpp
Implementation of linear interface element.
MoFEM::CoreInterface::get_finite_element_meshset
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
BodyForceConstantField::addBlock
MoFEMErrorCode addBlock(const std::string field_name, Vec F, int ms_id)
Definition: BodyForce.hpp:94
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:205
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:235
CohesiveElement::ArcLengthIntElemFEMethod
Definition: InterfaceGapArcLengthControl.hpp:12
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation
Constitutive (physical) equation for interface.
Definition: CohesiveInterfaceElement.hpp:51
EntityHandle
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
CohesiveElement::ArcLengthElement::postProcessLoadPath
MoFEMErrorCode postProcessLoadPath()
Definition: arc_length_interface.cpp:50
MetaNodalForces::addElement
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
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:674
MoFEM::Mat_Interf::data
_data_ data
Definition: MaterialBlocks.hpp:441
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
MoFEM::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition: PrismInterface.hpp:23
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
ArcLengthCtx
Store variables for ArcLength analysis.
Definition: ArcLengthTools.hpp:65
NOBASE
@ NOBASE
Definition: definitions.h:59
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
help
static char help[]
Definition: arc_length_interface.cpp:21
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:13
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::SnesCtx::FEMethodsSequence
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:27
MoFEM::CoreInterface::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
BodyForceConstantField::getLoopFe
MyVolumeFE & getLoopFe()
Definition: BodyForce.hpp:23
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
CohesiveElement::AssembleRhsVectors::bodyForce
Vec & bodyForce
Definition: arc_length_interface.cpp:100
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
order
constexpr int order
Definition: dg_projection.cpp:18
PCSetupArcLength
MoFEMErrorCode PCSetupArcLength(PC pc)
Definition: ArcLengthTools.cpp:326
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::Mat_Interf
Linear interface data structure.
Definition: MaterialBlocks.hpp:430
NODESET
@ NODESET
Definition: definitions.h:159
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CohesiveElement::AssembleRhsVectors::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: arc_length_interface.cpp:107
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:304
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
BodyForceConstantField
Body forces elements.
Definition: BodyForce.hpp:12
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
CohesiveElement::ArcLengthElement
Definition: arc_length_interface.cpp:31
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition: ArcLengthTools.cpp:184
CohesiveElement::CohesiveInterfaceElement::getFeLhs
MyPrism & getFeLhs()
Definition: CohesiveInterfaceElement.hpp:39
CohesiveElement::ArcLengthElement::ArcLengthElement
ArcLengthElement(MoFEM::Interface &m_field, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
Definition: arc_length_interface.cpp:34
FORCESET
@ FORCESET
Definition: definitions.h:164
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::add_ents_to_finite_element_by_bit_ref
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
CohesiveElement::CohesiveInterfaceElement::addOps
MoFEMErrorCode addOps(const std::string field_name, boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation > &interfaces)
Driver function settting all operators needed for interface element.
Definition: CohesiveInterfaceElement.hpp:568
Hooke.hpp
Implementation of linear elastic material.
BODYFORCESSET
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:175
InterfaceGapArcLengthControl.hpp
Implementation of arc-lebgth control for cohesive elements.
PostProcTemplateVolumeOnRefinedMesh::commonData
CommonData commonData
Definition: PostProcOnRefMesh.hpp:287
CohesiveElement::CohesiveInterfaceElement
Cohesive element implementation.
Definition: CohesiveInterfaceElement.hpp:14
MoFEM::PrismInterface::splitSides
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
Definition: PrismInterface.cpp:519
DATAFILENAME
#define DATAFILENAME
Definition: arc_length_interface.cpp:27
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
COL
@ COL
Definition: definitions.h:136
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcOnRefMesh.hpp:253
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
INTERFACESET
@ INTERFACESET
Definition: definitions.h:170
CohesiveElement::AssembleRhsVectors::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: arc_length_interface.cpp:101
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:195
CohesiveElement::CohesiveInterfaceElement::getFeRhs
MyPrism & getFeRhs()
Definition: CohesiveInterfaceElement.hpp:38
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::ProblemsManager::partitionFiniteElements
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
Definition: ProblemsManager.cpp:2167
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:325
CohesiveElement::AssembleRhsVectors::AssembleRhsVectors
AssembleRhsVectors(MoFEM::Interface &m_field, Vec &body_force, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
Definition: arc_length_interface.cpp:103
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
NonlinearElasticElement::setOperators
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.
Definition: NonLinearElasticElement.cpp:1153
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
MoFEM::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:290
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
CohesiveElement::ArcLengthElement::mField
MoFEM::Interface & mField
Definition: arc_length_interface.cpp:32
MoFEM::PrismInterface::getSides
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
Definition: PrismInterface.cpp:56
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
CohesiveElement::CohesiveInterfaceElement::getFeHistory
MyPrism & getFeHistory()
Definition: CohesiveInterfaceElement.hpp:40
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
CohesiveElement::AssembleRhsVectors::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: arc_length_interface.cpp:125
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
CohesiveElement::ArcLengthElement::postProcNodes
Range postProcNodes
Definition: arc_length_interface.cpp:33
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::BasicMethod::problemPtr
const Problem * problemPtr
raw pointer to problem
Definition: LoopMethods.hpp:250
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
Definition: MeshsetsManager.hpp:94
CohesiveElement::AssembleRhsVectors::mField
MoFEM::Interface & mField
Definition: arc_length_interface.cpp:99
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
CohesiveElement
Definition: arc_length_interface.cpp:29
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreInterface::modify_problem_add_finite_element
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
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
PCApplyArcLength
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Definition: ArcLengthTools.cpp:243
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
NeumannForcesSurface
Finite element and operators to apply force/pressures applied to surfaces.
Definition: SurfacePressure.hpp:14
MoFEM::SnesMat
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:139
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
CohesiveElement::AssembleRhsVectors
Definition: arc_length_interface.cpp:97
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
DirichletDisplacementBc::iNitialize
virtual MoFEMErrorCode iNitialize()
Definition: DirichletBC.cpp:198
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
DirichletDisplacementBc::dofsIndices
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:71
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1683
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:153
Generated by Doxygen 1.8.17 and hosted at University of Glasgow