v0.13.2
Loading...
Searching...
No Matches
Functions | Variables
arc_length_nonlinear_elasticity.cpp File Reference
User modules » NonLinear Elastic Element

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,
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.get_loops_to_do_Rhs();
606 snes_ctx.get_preProcess_to_do_Rhs().push_back(my_dirichlet_bc);
607 snes_ctx.get_preProcess_to_do_Rhs().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.get_postProcess_to_do_Rhs().push_back(&pre_post_method);
656 snes_ctx.get_postProcess_to_do_Rhs().push_back(my_dirichlet_bc);
657
658 SnesCtx::FEMethodsSequence &loops_to_do_Mat =
659 snes_ctx.get_loops_to_do_Mat();
660 snes_ctx.get_preProcess_to_do_Mat().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.get_postProcess_to_do_Mat().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}
PCApplyArcLength
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Definition: ArcLengthTools.cpp:243
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition: ArcLengthTools.cpp:184
PCSetupArcLength
MoFEMErrorCode PCSetupArcLength(PC pc)
Definition: ArcLengthTools.cpp:326
M
static Index< 'M', 3 > M
Definition: BasicFeTools.hpp:96
default_options
const std::string default_options
Definition: DefaultParams.hpp:1
param_file
std::string param_file
Definition: DefaultParams.hpp:113
help
static char help[]
Definition: arc_length_nonlinear_elasticity.cpp:10
COL
@ COL
Definition: definitions.h:123
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
NOBASE
@ NOBASE
Definition: definitions.h:59
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
H1
@ H1
continuous field
Definition: definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
PRESSURESET
@ PRESSURESET
Definition: definitions.h:152
FORCESET
@ FORCESET
Definition: definitions.h:151
NODESET
@ NODESET
Definition: definitions.h:146
SIDESET
@ SIDESET
Definition: definitions.h:147
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
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::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
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::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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::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
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
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::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.
_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
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2341
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::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1685
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:2169
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::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::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
D
double D
Definition: initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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::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:136
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:730
N
const int N
Definition: speed_test.cpp:3
ArcLengthCtx
Store variables for ArcLength analysis.
Definition: ArcLengthTools.hpp:65
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:205
AssembleRhsVectors::preProcess
PetscErrorCode preProcess()
function is run at the beginning of loop
Definition: solid_shell_nonlinear.cpp:102
AssembleRhsVectors::postProcess
PetscErrorCode postProcess()
function is run at the end of loop
Definition: solid_shell_nonlinear.cpp:118
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211
EdgeForce
Force on edges and lines.
Definition: EdgeForce.hpp:13
ElasticMaterials
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Definition: ElasticMaterials.hpp:24
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
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
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::BasicMethod::postProcess
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: LoopMethods.cpp:149
MoFEM::BasicMethod::operator()
virtual MoFEMErrorCode operator()()
function is run for every finite element
Definition: LoopMethods.cpp:160
MoFEM::BasicMethod::preProcess
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: LoopMethods.cpp:138
MoFEM::BasicMethod::problemPtr
const Problem * problemPtr
raw pointer to problem
Definition: LoopMethods.hpp:250
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Ent_mi_tag
Definition: TagMultiIndices.hpp:21
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:216
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:237
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:202
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:68
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::SnesCtx::FEMethodsSequence
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:27
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
MoFEM::SnesMethod::snes_ctx
SNESContext snes_ctx
Definition: LoopMethods.hpp:118
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
NeumannForcesSurfaceComplexForLazy
NonLinear surface pressure element (obsolete implementation)
Definition: SurfacePressureComplexForLazy.hpp:16
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:13
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:235
PostProcStress
Definition: PostProcStresses.hpp:18
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:957
SphericalArcLengthControl
Implementation of spherical arc-length method.
Definition: ArcLengthTools.hpp:374

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.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow