v0.13.1
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[] 
)

Definition at line 35 of file arc_length_nonlinear_elasticity.cpp.

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