v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
arc_length_nonlinear_elasticity.cpp File Reference

nonlinear elasticity (arc-length control) More...

#include <BasicFiniteElements.hpp>
#include <ElasticMaterials.hpp>
#include <NeoHookean.hpp>
#include <SurfacePressureComplexForLazy.hpp>

Go to the source code of this file.

Functions

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

Variables

static char help []
 

Detailed Description

nonlinear elasticity (arc-length control)

Solves nonlinear elastic problem. Using arc length control.

Definition in file arc_length_nonlinear_elasticity.cpp.

Function Documentation

◆ main()

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

< do not very if element of given name exist when do loop over elements

Definition at line 23 of file arc_length_nonlinear_elasticity.cpp.

23 {
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,
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,
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
159 "SPATIAL_POSITION");
160 CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
162 "SPATIAL_POSITION");
164 "ELASTIC", "LAMBDA"); // this is for parmetis
166 "SPATIAL_POSITION");
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
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");
254 "NEUMANN_FE", "MESH_NODE_POSITIONS");
255 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
256 "NEUMANN_FE");
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 }
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(
281 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
283
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
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;
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);
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);
414 m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
416 neumann_forces.getLoopSpatialFe();
417 if (linear) {
420 }
421 fe_neumann.uSeF = true;
423 it)) {
424 CHKERR fe_neumann.addForce(it->getMeshsetId());
425 }
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
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
468 }
469
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 }
484 }
485
486 MoFEMErrorCode potsProcessLoadPath() {
488 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
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 }
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
522 }
526 }
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 }
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));
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));
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 }
877
879
880 return 0;
881}
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
MoFEMErrorCode PCSetupArcLength(PC pc)
static Index< 'M', 3 > M
const std::string default_options
std::string param_file
static char help[]
@ COL
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ NOBASE
Definition: definitions.h:59
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ PRESSURESET
Definition: definitions.h:152
@ FORCESET
Definition: definitions.h:151
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ F
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode add_ents_to_finite_element_by_MESHSET(const EntityHandle meshset, const std::string &name, const bool recursive=false)=0
add MESHSET element to finite element database given by name
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
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
double D
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const int N
Definition: speed_test.cpp:3
Store variables for ArcLength analysis.
shell matrix for arc-length method
PetscErrorCode preProcess()
function is run at the beginning of loop
PetscErrorCode postProcess()
function is run at the end of loop
Set Dirichlet boundary conditions on spatial displacements.
Force on edges and lines.
Definition: EdgeForce.hpp:13
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
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
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.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
virtual MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
const Problem * problemPtr
raw pointer to problem
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
structure for User Loop Methods on finite elements
Basic algebra on fields.
Definition: FieldBlas.hpp:21
Matrix manager is used to build and partition problems.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:27
Vec & snes_f
residual
SNESContext snes_ctx
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
NonLinear surface pressure element (obsolete implementation)
Force applied to nodes.
Definition: NodalForce.hpp:13
structure grouping operators and data used for calculation of nonlinear elastic element
structure for Arc Length pre-conditioner
Post processing.
Implementation of spherical arc-length method.

Variable Documentation

◆ help

char help[]
static
Initial value:
= "\
-my_file mesh file name\n\
-my_sr reduction of step size\n\
-my_ms maximal number of steps\n\n"

Definition at line 10 of file arc_length_nonlinear_elasticity.cpp.