v0.14.0
arc_length_nonlinear_elasticity.cpp
Go to the documentation of this file.
1 /** \file arc_length_nonlinear_elasticity.cpp
2  * \ingroup nonlinear_elastic_elem
3  * \brief nonlinear elasticity (arc-length control)
4  *
5  * Solves nonlinear elastic problem. Using arc length control.
6  */
7 
8 
9 
10 static char help[] = "\
11  -my_file mesh file name\n\
12  -my_sr reduction of step size\n\
13  -my_ms maximal number of steps\n\n";
14 
15 #include <BasicFiniteElements.hpp>
16 using namespace MoFEM;
17 
18 #include <ElasticMaterials.hpp>
19 #include <NeoHookean.hpp>
20 
21 #include <SurfacePressureComplexForLazy.hpp>
22 
23 int main(int argc, char *argv[]) {
24 
25  const string default_options = "-ksp_type fgmres \n"
26  "-pc_type lu \n"
27  "-pc_factor_mat_solver_type mumps \n"
28  "-mat_mumps_icntl_20 0 \n"
29  "-ksp_atol 1e-10 \n"
30  "-ksp_rtol 1e-10 \n"
31  "-snes_monitor \n"
32  "-snes_type newtonls \n"
33  "-snes_linesearch_type basic \n"
34  "-snes_max_it 100 \n"
35  "-snes_atol 1e-7 \n"
36  "-snes_rtol 1e-7 \n"
37  "-ts_monitor \n"
38  "-ts_type alpha \n";
39 
40  string param_file = "param_file.petsc";
41  if (!static_cast<bool>(ifstream(param_file))) {
42  std::ofstream file(param_file.c_str(), std::ios::ate);
43  if (file.is_open()) {
44  file << default_options;
45  file.close();
46  }
47  }
48 
49  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
50 
51 
52  try {
53 
54  moab::Core mb_instance;
55  moab::Interface &moab = mb_instance;
56 
57  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
58  auto moab_comm_wrap =
59  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
60  if (pcomm == NULL)
61  pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
62 
63  PetscBool flg = PETSC_TRUE;
64  char mesh_file_name[255];
65  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
66  mesh_file_name, 255, &flg);
67  if (flg != PETSC_TRUE) {
68  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
69  }
70 
71  PetscInt order;
72  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
73  &flg);
74  if (flg != PETSC_TRUE) {
75  order = 2;
76  }
77 
78  // use this if your mesh is partitioned and you run code on parts,
79  // you can solve very big problems
80  PetscBool is_partitioned = PETSC_FALSE;
81  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_partitioned",
82  &is_partitioned, &flg);
83 
84  if (is_partitioned == PETSC_TRUE) {
85  // Read mesh to MOAB
86  const char *option;
87  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
88  "PARALLEL_PARTITION;";
89  CHKERR moab.load_file(mesh_file_name, 0, option);
90  } else {
91  const char *option;
92  option = "";
93  CHKERR moab.load_file(mesh_file_name, 0, option);
94  }
95 
96  // data stored on mesh for restart
97  Tag th_step_size, th_step;
98  double def_step_size = 1;
99  CHKERR moab.tag_get_handle("_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
100  MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
101  if (rval == MB_ALREADY_ALLOCATED)
102  CHKERR MB_SUCCESS;
103 
104  int def_step = 1;
105  CHKERR moab.tag_get_handle("_STEP", 1, MB_TYPE_INTEGER, th_step,
106  MB_TAG_CREAT | MB_TAG_MESH, &def_step);
107  if (rval == MB_ALREADY_ALLOCATED)
108  CHKERR MB_SUCCESS;
109 
110  const void *tag_data_step_size[1];
111  EntityHandle root = moab.get_root_set();
112  CHKERR moab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
113  double &step_size = *(double *)tag_data_step_size[0];
114  const void *tag_data_step[1];
115  CHKERR moab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
116  int &step = *(int *)tag_data_step[0];
117  // end of data stored for restart
118  CHKERR PetscPrintf(PETSC_COMM_WORLD,
119  "Start step %D and step_size = %6.4e\n", step,
120  step_size);
121 
122  MoFEM::Core core(moab);
123  MoFEM::Interface &m_field = core;
124 
125  // ref meshset ref level 0
126  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
127  0, 3, BitRefLevel().set(0));
128  std::vector<BitRefLevel> bit_levels;
129  bit_levels.push_back(BitRefLevel().set(0));
130  BitRefLevel problem_bit_level;
131 
132  if (step == 1) {
133 
134  problem_bit_level = bit_levels.back();
135 
136  // Fields
137  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE,
138  3);
139  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1,
140  AINSWORTH_LEGENDRE_BASE, 3);
141 
142  CHKERR m_field.add_field("LAMBDA", NOFIELD, NOBASE, 1);
143 
144  // Field for ArcLength
145  CHKERR m_field.add_field("X0_SPATIAL_POSITION", H1,
146  AINSWORTH_LEGENDRE_BASE, 3);
147 
148  // FE
149  CHKERR m_field.add_finite_element("ELASTIC");
150  CHKERR m_field.add_finite_element("ARC_LENGTH");
151 
152  // Add spring boundary condition applied on surfaces.
153  // This is only declaration not implementation.
154  CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
155  "MESH_NODE_POSITIONS");
156 
157  // Define rows/cols and element data
158  CHKERR m_field.modify_finite_element_add_field_row("ELASTIC",
159  "SPATIAL_POSITION");
160  CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
161  CHKERR m_field.modify_finite_element_add_field_col("ELASTIC",
162  "SPATIAL_POSITION");
163  CHKERR m_field.modify_finite_element_add_field_col(
164  "ELASTIC", "LAMBDA"); // this is for parmetis
165  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
166  "SPATIAL_POSITION");
167  CHKERR m_field.modify_finite_element_add_field_data(
168  "ELASTIC", "MESH_NODE_POSITIONS");
169  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "LAMBDA");
170 
171  // Define rows/cols and element data
172  CHKERR m_field.modify_finite_element_add_field_row("ARC_LENGTH",
173  "LAMBDA");
174  CHKERR m_field.modify_finite_element_add_field_col("ARC_LENGTH",
175  "LAMBDA");
176  // elem data
177  CHKERR m_field.modify_finite_element_add_field_data("ARC_LENGTH",
178  "LAMBDA");
179 
180  // define problems
181  CHKERR m_field.add_problem("ELASTIC_MECHANICS");
182 
183  // set finite elements for problems
184  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
185  "ELASTIC");
186  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
187  "ARC_LENGTH");
188 
189  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
190  "SPRING");
191 
192  // set refinement level for problem
193  CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
194  problem_bit_level);
195 
196  // add entities (by tets) to the field
197  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
198  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "X0_SPATIAL_POSITION");
199  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
200 
201  // Setting up LAMBDA field and ARC_LENGTH interface
202  {
203  // Add dummy no-field vertex
204  EntityHandle no_field_vertex;
205  {
206  const double coords[] = {0, 0, 0};
207  CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
208  Range range_no_field_vertex;
209  range_no_field_vertex.insert(no_field_vertex);
210  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
211  range_no_field_vertex, BitRefLevel().set());
212  EntityHandle lambda_meshset = m_field.get_field_meshset("LAMBDA");
213  CHKERR m_field.get_moab().add_entities(lambda_meshset,
214  range_no_field_vertex);
215  }
216  // this entity will carry data for this finite element
217  EntityHandle meshset_fe_arc_length;
218  {
219  CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
220  CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
221  CHKERR m_field.getInterface<BitRefManager>()->setBitLevelToMeshset(
222  meshset_fe_arc_length, BitRefLevel().set());
223  }
224  // finally add created meshset to the ARC_LENGTH finite element
225  CHKERR m_field.add_ents_to_finite_element_by_MESHSET(
226  meshset_fe_arc_length, "ARC_LENGTH", false);
227  }
228 
229  // set app. order
230  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
231  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
232  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
233  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
234 
235  CHKERR m_field.set_field_order(0, MBTET, "X0_SPATIAL_POSITION", order);
236  CHKERR m_field.set_field_order(0, MBTRI, "X0_SPATIAL_POSITION", order);
237  CHKERR m_field.set_field_order(0, MBEDGE, "X0_SPATIAL_POSITION", order);
238  CHKERR m_field.set_field_order(0, MBVERTEX, "X0_SPATIAL_POSITION", 1);
239 
240  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
241  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
242  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
243  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
244 
245  // add Neumman finite elements to add static boundary conditions
246  CHKERR m_field.add_finite_element("NEUMANN_FE");
247  CHKERR m_field.modify_finite_element_add_field_row("NEUMANN_FE",
248  "SPATIAL_POSITION");
249  CHKERR m_field.modify_finite_element_add_field_col("NEUMANN_FE",
250  "SPATIAL_POSITION");
251  CHKERR m_field.modify_finite_element_add_field_data("NEUMANN_FE",
252  "SPATIAL_POSITION");
253  CHKERR m_field.modify_finite_element_add_field_data(
254  "NEUMANN_FE", "MESH_NODE_POSITIONS");
255  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
256  "NEUMANN_FE");
257  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
258  NODESET | FORCESET, it)) {
259  Range tris;
260  CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, true);
261  CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI,
262  "NEUMANN_FE");
263  }
264  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
265  m_field, SIDESET | PRESSURESET, it)) {
266  Range tris;
267  CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, true);
268  CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI,
269  "NEUMANN_FE");
270  }
271  // add nodal force element
272  CHKERR MetaNodalForces::addElement(m_field, "SPATIAL_POSITION");
273  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
274  "FORCE_FE");
275  }
276 
277  // Implementation of spring element
278  // Create new instances of face elements for springs
279  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
280  new FaceElementForcesAndSourcesCore(m_field));
281  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
282  new FaceElementForcesAndSourcesCore(m_field));
283 
284  CHKERR MetaSpringBC::setSpringOperators(
285  m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
286  "MESH_NODE_POSITIONS");
287 
288  PetscBool linear;
289  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-is_linear", &linear,
290  &linear);
291 
292  NonlinearElasticElement elastic(m_field, 2);
293  ElasticMaterials elastic_materials(m_field);
294  CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
295  CHKERR elastic.addElement("ELASTIC", "SPATIAL_POSITION");
296  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeRhs(), true, false,
297  false, false);
298  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeLhs(), true, false,
299  false, false);
300  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeEnergy(), true,
301  false, false, false);
302  CHKERR elastic.setOperators("SPATIAL_POSITION");
303 
304  // post_processing
305  PostProcVolumeOnRefinedMesh post_proc(m_field);
306  CHKERR post_proc.generateReferenceElementMesh();
307  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, true, false, false,
308  false);
309  CHKERR post_proc.addFieldValuesPostProc("SPATIAL_POSITION");
310  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
311  CHKERR post_proc.addFieldValuesGradientPostProc("SPATIAL_POSITION");
312  std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
313  elastic.setOfBlocks.begin();
314  for (; sit != elastic.setOfBlocks.end(); sit++) {
315  post_proc.getOpPtrVector().push_back(new PostProcStress(
316  post_proc.postProcMesh, post_proc.mapGaussPts, "SPATIAL_POSITION",
317  sit->second, post_proc.commonData));
318  }
319 
320  // build field
321  CHKERR m_field.build_fields();
322  if (step == 1) {
323  // 10 node tets
324  Projection10NodeCoordsOnField ent_method_material(m_field,
325  "MESH_NODE_POSITIONS");
326  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material, 0);
327  CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBVERTEX,
328  "SPATIAL_POSITION");
329  CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBEDGE,
330  "SPATIAL_POSITION");
331  CHKERR m_field.getInterface<FieldBlas>()->fieldAxpy(
332  1., "MESH_NODE_POSITIONS", "SPATIAL_POSITION");
333  CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBTRI,
334  "SPATIAL_POSITION");
335  CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBTET,
336  "SPATIAL_POSITION");
337  }
338 
339  // build finite elements
340  CHKERR m_field.build_finite_elements();
341 
342  // build adjacencies
343  CHKERR m_field.build_adjacencies(problem_bit_level);
344 
345  ProblemsManager *prb_mng_ptr;
346  CHKERR m_field.getInterface(prb_mng_ptr);
347  // build database
348  if (is_partitioned) {
349  SETERRQ(PETSC_COMM_SELF, 1,
350  "Not implemented, problem with arc-length force multiplayer");
351  } else {
352  CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
353  CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
354  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS");
355  }
356  CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
357 
358  // print bcs
359  MeshsetsManager *mmanager_ptr;
360  CHKERR m_field.getInterface(mmanager_ptr);
361  CHKERR mmanager_ptr->printDisplacementSet();
362  CHKERR mmanager_ptr->printForceSet();
363  // print block sets with materials
364  CHKERR mmanager_ptr->printMaterialsSet();
365 
366  // create matrices
367  Vec F;
368  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
369  "ELASTIC_MECHANICS", COL, &F);
370  Vec D;
371  CHKERR VecDuplicate(F, &D);
372  Mat Aij;
373  CHKERR m_field.getInterface<MatrixManager>()
374  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
375  &Aij);
376 
377  boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
378  new ArcLengthCtx(m_field, "ELASTIC_MECHANICS"));
379 
380  PetscInt M, N;
381  CHKERR MatGetSize(Aij, &M, &N);
382  PetscInt m, n;
383  CHKERR MatGetLocalSize(Aij, &m, &n);
384  boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
385  new ArcLengthMatShell(Aij, arc_ctx, "ELASTIC_MECHANICS"));
386 
387  Mat ShellAij;
388  CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, mat_ctx.get(),
389  &ShellAij);
390  CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
391  (void (*)(void))ArcLengthMatMultShellOp);
392 
393  ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
394  ///< do not very if element of given name exist when do loop over elements
395  snes_ctx.bH = MF_ZERO;
396 
397  Range node_set;
398  for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field, "LoadPath", cit)) {
399  EntityHandle meshset = cit->getMeshset();
400  Range nodes;
401  CHKERR moab.get_entities_by_type(meshset, MBVERTEX, nodes, true);
402  MOAB_THROW(rval);
403  node_set.merge(nodes);
404  }
405  PetscPrintf(PETSC_COMM_WORLD, "Nb. nodes in load path: %u\n",
406  node_set.size());
407 
408  SphericalArcLengthControl arc_method(arc_ctx);
409 
410  double scaled_reference_load = 1;
411  double *scale_lhs = &(arc_ctx->getFieldData());
412  double *scale_rhs = &(scaled_reference_load);
413  NeumannForcesSurfaceComplexForLazy neumann_forces(
414  m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
415  NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE &fe_neumann =
416  neumann_forces.getLoopSpatialFe();
417  if (linear) {
418  fe_neumann.typeOfForces = NeumannForcesSurfaceComplexForLazy::
419  MyTriangleSpatialFE::NONCONSERVATIVE;
420  }
421  fe_neumann.uSeF = true;
422  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
423  it)) {
424  CHKERR fe_neumann.addForce(it->getMeshsetId());
425  }
426  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
427  m_field, SIDESET | PRESSURESET, it)) {
428  CHKERR fe_neumann.addPressure(it->getMeshsetId());
429  }
430 
431  boost::shared_ptr<FEMethod> my_dirichlet_bc =
432  boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
433  m_field, "SPATIAL_POSITION", Aij, D, F));
434  CHKERR m_field.get_problem("ELASTIC_MECHANICS",
435  &(my_dirichlet_bc->problemPtr));
436  CHKERR dynamic_cast<DirichletSpatialPositionsBc *>(my_dirichlet_bc.get())
437  ->iNitialize();
438 
439  struct AssembleRhsVectors : public FEMethod {
440 
441  boost::shared_ptr<ArcLengthCtx> arcPtr;
442  Range &nodeSet;
443 
444  AssembleRhsVectors(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
445  Range &node_set)
446  : arcPtr(arc_ptr), nodeSet(node_set) {}
447 
448  MoFEMErrorCode preProcess() {
449  MoFEMFunctionBeginHot;
450 
451  // PetscAttachDebugger();
452  switch (snes_ctx) {
453  case CTX_SNESSETFUNCTION: {
454  CHKERR VecZeroEntries(snes_f);
455  CHKERR VecGhostUpdateBegin(snes_f, INSERT_VALUES, SCATTER_FORWARD);
456  CHKERR VecGhostUpdateEnd(snes_f, INSERT_VALUES, SCATTER_FORWARD);
457  CHKERR VecZeroEntries(arcPtr->F_lambda);
458  CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, INSERT_VALUES,
459  SCATTER_FORWARD);
460  CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, INSERT_VALUES,
461  SCATTER_FORWARD);
462  } break;
463  default:
464  SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
465  }
466 
467  MoFEMFunctionReturnHot(0);
468  }
469 
470  MoFEMErrorCode postProcess() {
471  MoFEMFunctionBeginHot;
472  switch (snes_ctx) {
473  case CTX_SNESSETFUNCTION: {
474  // snes_f
475  CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
476  CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
477  CHKERR VecAssemblyBegin(snes_f);
478  CHKERR VecAssemblyEnd(snes_f);
479  } break;
480  default:
481  SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
482  }
483  MoFEMFunctionReturnHot(0);
484  }
485 
486  MoFEMErrorCode potsProcessLoadPath() {
487  MoFEMFunctionBeginHot;
488  boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
489  problemPtr->getNumeredRowDofsPtr();
490  Range::iterator nit = nodeSet.begin();
491  for (; nit != nodeSet.end(); nit++) {
492  NumeredDofEntityByEnt::iterator dit, hi_dit;
493  dit = numered_dofs_rows->get<Ent_mi_tag>().lower_bound(*nit);
494  hi_dit = numered_dofs_rows->get<Ent_mi_tag>().upper_bound(*nit);
495  for (; dit != hi_dit; dit++) {
496  PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e -> ", "LAMBDA", 0,
497  arcPtr->getFieldData());
498  PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e\n",
499  dit->get()->getName().c_str(),
500  dit->get()->getDofCoeffIdx(),
501  dit->get()->getFieldData());
502  }
503  }
504  MoFEMFunctionReturnHot(0);
505  }
506  };
507 
508  struct AddLambdaVectorToFInternal : public FEMethod {
509 
510  boost::shared_ptr<ArcLengthCtx> arcPtr;
511  boost::shared_ptr<DirichletSpatialPositionsBc> bC;
512 
513  AddLambdaVectorToFInternal(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
514  boost::shared_ptr<FEMethod> &bc)
515  : arcPtr(arc_ptr),
516  bC(boost::shared_ptr<DirichletSpatialPositionsBc>(
517  bc, dynamic_cast<DirichletSpatialPositionsBc *>(bc.get()))) {}
518 
519  MoFEMErrorCode preProcess() {
520  MoFEMFunctionBeginHot;
521  MoFEMFunctionReturnHot(0);
522  }
523  MoFEMErrorCode operator()() {
524  MoFEMFunctionBeginHot;
525  MoFEMFunctionReturnHot(0);
526  }
527  MoFEMErrorCode postProcess() {
528  MoFEMFunctionBeginHot;
529  switch (snes_ctx) {
530  case CTX_SNESSETFUNCTION: {
531  // F_lambda
532  CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES,
533  SCATTER_REVERSE);
534  CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES,
535  SCATTER_REVERSE);
536  CHKERR VecAssemblyBegin(arcPtr->F_lambda);
537  CHKERR VecAssemblyEnd(arcPtr->F_lambda);
538  for (std::vector<int>::iterator vit = bC->dofsIndices.begin();
539  vit != bC->dofsIndices.end(); vit++) {
540  CHKERR VecSetValue(arcPtr->F_lambda, *vit, 0, INSERT_VALUES);
541  }
542  CHKERR VecAssemblyBegin(arcPtr->F_lambda);
543  CHKERR VecAssemblyEnd(arcPtr->F_lambda);
544  CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
545  PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n",
546  arcPtr->F_lambda2);
547  // add F_lambda
548  CHKERR VecAXPY(snes_f, arcPtr->getFieldData(), arcPtr->F_lambda);
549  PetscPrintf(PETSC_COMM_WORLD, "\tlambda = %6.4e\n",
550  arcPtr->getFieldData());
551  double fnorm;
552  CHKERR VecNorm(snes_f, NORM_2, &fnorm);
553  PetscPrintf(PETSC_COMM_WORLD, "\tfnorm = %6.4e\n", fnorm);
554  } break;
555  default:
556  SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
557  }
558  MoFEMFunctionReturnHot(0);
559  }
560  };
561 
562  AssembleRhsVectors pre_post_method(arc_ctx, node_set);
563  AddLambdaVectorToFInternal assemble_F_lambda(arc_ctx, my_dirichlet_bc);
564 
565  SNES snes;
566  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
567  CHKERR SNESSetApplicationContext(snes, &snes_ctx);
568  CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
569  CHKERR SNESSetJacobian(snes, ShellAij, Aij, SnesMat, &snes_ctx);
570  CHKERR SNESSetFromOptions(snes);
571 
572  PetscReal my_tol;
573  CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_tol", &my_tol,
574  &flg);
575  if (flg == PETSC_TRUE) {
576  PetscReal atol, rtol, stol;
577  PetscInt maxit, maxf;
578  CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
579  atol = my_tol;
580  rtol = atol * 1e2;
581  CHKERR SNESSetTolerances(snes, atol, rtol, stol, maxit, maxf);
582  }
583 
584  KSP ksp;
585  CHKERR SNESGetKSP(snes, &ksp);
586  PC pc;
587  CHKERR KSPGetPC(ksp, &pc);
588  boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
589  new PCArcLengthCtx(ShellAij, Aij, arc_ctx));
590  CHKERR PCSetType(pc, PCSHELL);
591  CHKERR PCShellSetContext(pc, pc_ctx.get());
592  CHKERR PCShellSetApply(pc, PCApplyArcLength);
593  CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
594 
595  if (flg == PETSC_TRUE) {
596  PetscReal rtol, atol, dtol;
597  PetscInt maxits;
598  CHKERR KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
599  atol = my_tol * 1e-2;
600  rtol = atol * 1e-2;
601  CHKERR KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
602  }
603 
604  SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
605  snes_ctx.getComputeRhs();
606  snes_ctx.getPreProcComputeRhs().push_back(my_dirichlet_bc);
607  snes_ctx.getPreProcComputeRhs().push_back(&pre_post_method);
608  loops_to_do_Rhs.push_back(
609  SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
610 
611  loops_to_do_Rhs.push_back(
612  SnesCtx::PairNameFEMethodPtr("SPRING", fe_spring_rhs_ptr.get()));
613 
614  // surface forces and pressures
615  loops_to_do_Rhs.push_back(
616  SnesCtx::PairNameFEMethodPtr("NEUMANN_FE", &fe_neumann));
617 
618  // edge forces
619  boost::ptr_map<std::string, EdgeForce> edge_forces;
620  string fe_name_str = "FORCE_FE";
621  edge_forces.insert(fe_name_str, new EdgeForce(m_field));
622  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
623  it)) {
624  CHKERR edge_forces.at(fe_name_str)
625  .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
626  }
627  for (boost::ptr_map<std::string, EdgeForce>::iterator eit =
628  edge_forces.begin();
629  eit != edge_forces.end(); eit++) {
630  loops_to_do_Rhs.push_back(
631  SnesCtx::PairNameFEMethodPtr(eit->first, &eit->second->getLoopFe()));
632  }
633 
634  // nodal forces
635  boost::ptr_map<std::string, NodalForce> nodal_forces;
636  // string fe_name_str ="FORCE_FE";
637  nodal_forces.insert(fe_name_str, new NodalForce(m_field));
638  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
639  it)) {
640  CHKERR nodal_forces.at(fe_name_str)
641  .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
642  }
643  for (boost::ptr_map<std::string, NodalForce>::iterator fit =
644  nodal_forces.begin();
645  fit != nodal_forces.end(); fit++) {
646  loops_to_do_Rhs.push_back(
647  SnesCtx::PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
648  }
649 
650  // arc length
651  loops_to_do_Rhs.push_back(
652  SnesCtx::PairNameFEMethodPtr("NONE", &assemble_F_lambda));
653  loops_to_do_Rhs.push_back(
654  SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", &arc_method));
655  snes_ctx.getPostProcComputeRhs().push_back(&pre_post_method);
656  snes_ctx.getPostProcComputeRhs().push_back(my_dirichlet_bc);
657 
658  SnesCtx::FEMethodsSequence &loops_to_do_Mat =
659  snes_ctx.getSetOperators();
660  snes_ctx.getPreProcSetOperators().push_back(my_dirichlet_bc);
661  loops_to_do_Mat.push_back(
662  SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
663 
664  loops_to_do_Mat.push_back(
665  SnesCtx::PairNameFEMethodPtr("SPRING", fe_spring_lhs_ptr.get()));
666 
667  loops_to_do_Mat.push_back(
668  SnesCtx::PairNameFEMethodPtr("NEUMANN_FE", &fe_neumann));
669  loops_to_do_Mat.push_back(
670  SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", &arc_method));
671  snes_ctx.getPostProcSetOperators().push_back(my_dirichlet_bc);
672 
673  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
674  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
675  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
676  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
677 
678  PetscScalar step_size_reduction;
679  CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_sr",
680  &step_size_reduction, &flg);
681  if (flg != PETSC_TRUE) {
682  step_size_reduction = 1.;
683  }
684 
685  PetscInt max_steps;
686  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_ms", &max_steps,
687  &flg);
688  if (flg != PETSC_TRUE) {
689  max_steps = 5;
690  }
691 
692  int its_d;
693  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_its_d", &its_d,
694  &flg);
695  if (flg != PETSC_TRUE) {
696  its_d = 4;
697  }
698  PetscScalar max_reduction = 10, min_reduction = 0.1;
699  CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_max_step_reduction",
700  &max_reduction, &flg);
701  CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_min_step_reduction",
702  &min_reduction, &flg);
703 
704  double gamma = 0.5, reduction = 1;
705  // step = 1;
706  if (step == 1) {
707  step_size = step_size_reduction;
708  } else {
709  reduction = step_size_reduction;
710  step++;
711  }
712  double step_size0 = step_size;
713 
714  if (step > 1) {
715  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
716  "ELASTIC_MECHANICS", "SPATIAL_POSITION", "X0_SPATIAL_POSITION", COL,
717  arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
718  double x0_nrm;
719  CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
720  CHKERR PetscPrintf(PETSC_COMM_WORLD,
721  "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
722  arc_ctx->dLambda);
723  CHKERR arc_ctx->setAlphaBeta(1, 0);
724  } else {
725  CHKERR arc_ctx->setS(step_size);
726  CHKERR arc_ctx->setAlphaBeta(0, 1);
727  }
728 
729  CHKERR SnesRhs(snes, D, F, &snes_ctx);
730 
731  Vec D0, x00;
732  CHKERR VecDuplicate(D, &D0);
733  CHKERR VecDuplicate(arc_ctx->x0, &x00);
734  bool converged_state = false;
735 
736  for (int jj = 0; step < max_steps; step++, jj++) {
737 
738  CHKERR VecCopy(D, D0);
739  CHKERR VecCopy(arc_ctx->x0, x00);
740 
741  if (step == 1) {
742 
743  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load Step %D step_size = %6.4e\n",
744  step, step_size);
745  CHKERR arc_ctx->setS(step_size);
746  CHKERR arc_ctx->setAlphaBeta(0, 1);
747  CHKERR VecCopy(D, arc_ctx->x0);
748  double dlambda;
749  CHKERR arc_method.calculateInitDlambda(&dlambda);
750  CHKERR arc_method.setDlambdaToX(D, dlambda);
751 
752  } else if (step == 2) {
753 
754  CHKERR arc_ctx->setAlphaBeta(1, 0);
755  CHKERR arc_method.calculateDxAndDlambda(D);
756  step_size = std::sqrt(arc_method.calculateLambdaInt());
757  step_size0 = step_size;
758  CHKERR arc_ctx->setS(step_size);
759  double dlambda = arc_ctx->dLambda;
760  double dx_nrm;
761  CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
762  CHKERR PetscPrintf(PETSC_COMM_WORLD,
763  "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
764  "dx_nrm = %6.4e dx2 = %6.4e\n",
765  step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
766  CHKERR VecCopy(D, arc_ctx->x0);
767  CHKERR VecAXPY(D, 1., arc_ctx->dx);
768  CHKERR arc_method.setDlambdaToX(D, dlambda);
769 
770  } else {
771 
772  if (jj == 0) {
773  step_size0 = step_size;
774  }
775 
776  CHKERR arc_method.calculateDxAndDlambda(D);
777  step_size *= reduction;
778  if (step_size > max_reduction * step_size0) {
779  step_size = max_reduction * step_size0;
780  } else if (step_size < min_reduction * step_size0) {
781  step_size = min_reduction * step_size0;
782  }
783  CHKERR arc_ctx->setS(step_size);
784  double dlambda = reduction * arc_ctx->dLambda;
785  double dx_nrm;
786  CHKERR VecScale(arc_ctx->dx, reduction);
787  CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
788  CHKERR PetscPrintf(PETSC_COMM_WORLD,
789  "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
790  "dx_nrm = %6.4e dx2 = %6.4e\n",
791  step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
792  CHKERR VecCopy(D, arc_ctx->x0);
793  CHKERR VecAXPY(D, 1., arc_ctx->dx);
794  CHKERR arc_method.setDlambdaToX(D, dlambda);
795  }
796 
797  CHKERR SNESSolve(snes, PETSC_NULL, D);
798  int its;
799  CHKERR SNESGetIterationNumber(snes, &its);
800  CHKERR PetscPrintf(PETSC_COMM_WORLD, "number of Newton iterations = %D\n",
801  its);
802 
803  SNESConvergedReason reason;
804  CHKERR SNESGetConvergedReason(snes, &reason);
805  if (reason < 0) {
806 
807  CHKERR VecCopy(D0, D);
808  CHKERR VecCopy(x00, arc_ctx->x0);
809 
810  double x0_nrm;
811  CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
812  CHKERR PetscPrintf(PETSC_COMM_WORLD,
813  "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
814  arc_ctx->dLambda);
815  CHKERR arc_ctx->setAlphaBeta(1, 0);
816 
817  reduction = 0.1;
818  converged_state = false;
819 
820  continue;
821 
822  } else {
823 
824  if (step > 1 && converged_state) {
825 
826  reduction = pow((double)its_d / (double)(its + 1), gamma);
827  if (step_size >= max_reduction * step_size0 && reduction > 1) {
828  reduction = 1;
829  } else if (step_size <= min_reduction * step_size0 && reduction < 1) {
830  reduction = 1;
831  }
832  CHKERR PetscPrintf(PETSC_COMM_WORLD, "reduction step_size = %6.4e\n",
833  reduction);
834  }
835 
836  // Save data on mesh
837  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
838  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
839  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
840  "ELASTIC_MECHANICS", "SPATIAL_POSITION", "X0_SPATIAL_POSITION", COL,
841  arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
842  converged_state = true;
843  }
844 
845  if (step % 1 == 0) {
846  // Save restart file
847  // #ifdef MOAB_HDF5_PARALLEL
848  // std::ostringstream sss;
849  // sss << "restart_" << step << ".h5m";
850  // CHKERR
851  // moab.write_file(sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART");
852  // #else
853  // #warning "No parallel HDF5, no writing restart file"
854  // #endif
855  // Save data on mesh
856  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
857  post_proc);
858  std::ostringstream o1;
859  o1 << "out_" << step << ".h5m";
860  CHKERR post_proc.writeFile(o1.str().c_str());
861  }
862 
863  CHKERR pre_post_method.potsProcessLoadPath();
864  }
865 
866  CHKERR VecDestroy(&D0);
867  CHKERR VecDestroy(&x00);
868 
869  // detroy matrices
870  CHKERR VecDestroy(&F);
871  CHKERR VecDestroy(&D);
872  CHKERR MatDestroy(&Aij);
873  CHKERR MatDestroy(&ShellAij);
874  CHKERR SNESDestroy(&snes);
875  }
876  CATCH_ERRORS;
877 
878  MoFEM::Core::Finalize();
879 
880  return 0;
881 }
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
MoFEM::Ent_mi_tag
Definition: TagMultiIndices.hpp:21
SIDESET
@ SIDESET
Definition: definitions.h:160
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.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
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
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
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
main
int main(int argc, char *argv[])
Definition: arc_length_nonlinear_elasticity.cpp:23
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:205
EdgeForce
Force on edges and lines.
Definition: EdgeForce.hpp:13
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:235
SphericalArcLengthControl::calculateLambdaInt
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
Definition: ArcLengthTools.cpp:685
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NeumannForcesSurfaceComplexForLazy
NonLinear surface pressure element (obsolete implementation)
Definition: SurfacePressureComplexForLazy.hpp:16
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
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
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
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::uSeF
bool uSeF
Definition: SurfacePressureComplexForLazy.hpp:47
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
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addForce
MoFEMErrorCode addForce(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:516
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
ElasticMaterials
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Definition: ElasticMaterials.hpp:24
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
SphericalArcLengthControl
Implementation of spherical arc-length method.
Definition: ArcLengthTools.hpp:374
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
help
static char help[]
Definition: arc_length_nonlinear_elasticity.cpp:10
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
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
SurfacePressureComplexForLazy.hpp
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
PostProcStress
Definition: PostProcStresses.hpp:17
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
SphericalArcLengthControl::calculateDxAndDlambda
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition: ArcLengthTools.cpp:742
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
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition: ArcLengthTools.cpp:184
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.
PostProcTemplateVolumeOnRefinedMesh::commonData
CommonData commonData
Definition: PostProcOnRefMesh.hpp:287
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
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
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
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
NonlinearElasticElement::addElement
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition: NonLinearElasticElement.cpp:1120
SphericalArcLengthControl::calculateInitDlambda
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
Definition: ArcLengthTools.cpp:769
ElasticMaterials.hpp
Elastic materials.
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE
Definition: SurfacePressureComplexForLazy.hpp:39
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::typeOfForces
FORCES typeOfForces
Definition: SurfacePressureComplexForLazy.hpp:45
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
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
NeumannForcesSurfaceComplexForLazy::getLoopSpatialFe
MyTriangleSpatialFE & getLoopSpatialFe()
Definition: SurfacePressureComplexForLazy.hpp:163
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
MetaSpringBC::setSpringOperators
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
Definition: SpringElement.cpp:1178
SphericalArcLengthControl::setDlambdaToX
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
Definition: ArcLengthTools.cpp:783
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
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
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::NONCONSERVATIVE
@ NONCONSERVATIVE
Definition: SurfacePressureComplexForLazy.hpp:43
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp: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
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
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
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addPressure
MoFEMErrorCode addPressure(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:531
_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
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
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
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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.
NeoHookean.hpp
Implementation of Neo-Hookean elastic material.
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
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
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
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.
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
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
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
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