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

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,
152 AINSWORTH_LEGENDRE_BASE, 3);
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,
158 AINSWORTH_LEGENDRE_BASE, 3);
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
170 CHKERR m_field.modify_finite_element_add_field_row("ELASTIC",
171 "SPATIAL_POSITION");
172 CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
173 CHKERR m_field.modify_finite_element_add_field_col("ELASTIC",
174 "SPATIAL_POSITION");
175 CHKERR m_field.modify_finite_element_add_field_col(
176 "ELASTIC", "LAMBDA"); // this is for parmetis
177 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
178 "SPATIAL_POSITION");
179 CHKERR m_field.modify_finite_element_add_field_data(
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
237 CHKERR m_field.add_ents_to_finite_element_by_MESHSET(
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");
265 CHKERR m_field.modify_finite_element_add_field_data(
266 "NEUMANN_FE", "MESH_NODE_POSITIONS");
267 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
268 "NEUMANN_FE");
269 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
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 }
276 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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(
292 new FaceElementForcesAndSourcesCore(m_field));
293 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
294 new FaceElementForcesAndSourcesCore(m_field));
295
296 CHKERR MetaSpringBC::setSpringOperators(
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
352 CHKERR m_field.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;
385 CHKERR m_field.getInterface<MatrixManager>()
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);
414 MOAB_THROW(rval);
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);
425 NeumannForcesSurfaceComplexForLazy neumann_forces(
426 m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
427 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE &fe_neumann =
428 neumann_forces.getLoopSpatialFe();
429 if (linear) {
430 fe_neumann.typeOfForces = NeumannForcesSurfaceComplexForLazy::
431 MyTriangleSpatialFE::NONCONSERVATIVE;
432 }
433 fe_neumann.uSeF = true;
434 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
435 it)) {
436 CHKERR fe_neumann.addForce(it->getMeshsetId());
437 }
438 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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
460 MoFEMErrorCode preProcess() {
461 MoFEMFunctionBeginHot;
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
479 MoFEMFunctionReturnHot(0);
480 }
481
482 MoFEMErrorCode postProcess() {
483 MoFEMFunctionBeginHot;
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 }
495 MoFEMFunctionReturnHot(0);
496 }
497
498 MoFEMErrorCode potsProcessLoadPath() {
499 MoFEMFunctionBeginHot;
500 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
501 problemPtr->getNumeredRowDofsPtr();
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 }
516 MoFEMFunctionReturnHot(0);
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
531 MoFEMErrorCode preProcess() {
532 MoFEMFunctionBeginHot;
533 MoFEMFunctionReturnHot(0);
534 }
535 MoFEMErrorCode operator()() {
536 MoFEMFunctionBeginHot;
537 MoFEMFunctionReturnHot(0);
538 }
539 MoFEMErrorCode postProcess() {
540 MoFEMFunctionBeginHot;
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 }
570 MoFEMFunctionReturnHot(0);
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));
634 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
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));
650 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
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 }
888 CATCH_ERRORS;
889
890 MoFEM::Core::Finalize();
891
892 return 0;
893}
PCApplyArcLength
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Definition: ArcLengthTools.cpp:255
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition: ArcLengthTools.cpp:196
PCSetupArcLength
MoFEMErrorCode PCSetupArcLength(PC pc)
Definition: ArcLengthTools.cpp:338
M
static Index< 'M', 3 > M
Definition: BasicFeTools.hpp:96
n
static Index< 'n', 3 > n
Definition: BasicFeTools.hpp:88
default_options
const std::string default_options
Definition: DefaultParams.hpp:1
param_file
std::string param_file
Definition: DefaultParams.hpp:94
help
static char help[]
Definition: arc_length_nonlinear_elasticity.cpp:22
COL
@ COL
Definition: definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
NOBASE
@ NOBASE
Definition: definitions.h:72
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
H1
@ H1
continuous field
Definition: definitions.h:98
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
FORCESET
@ FORCESET
Definition: definitions.h:164
NODESET
@ NODESET
Definition: definitions.h:159
SIDESET
@ SIDESET
Definition: definitions.h:160
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
order
int order
Definition: fluid_structure_eigenproblem.cpp:68
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
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::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::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::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_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::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::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:104
_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:59
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2343
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:101
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1687
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:2171
MoFEM::CoreInterface::modify_problem_add_finite_element
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
MoFEM::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
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:35
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:78
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
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:148
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:153
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:163
MoFEM::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:39
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:193
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1096
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:183
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:484
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:32
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:95
N
const int N
Definition: speed_test.cpp:3
ArcLengthCtx
Store variables for ArcLength analysis.
Definition: ArcLengthTools.hpp:77
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:217
AssembleRhsVectors::preProcess
PetscErrorCode preProcess()
function is run at the beginning of loop
Definition: solid_shell_nonlinear.cpp:102
AssembleRhsVectors::AssembleRhsVectors
AssembleRhsVectors(ArcLengthCtx *arc_ptr)
Definition: solid_shell_nonlinear.cpp:101
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:217
EdgeForce
Force on edges and lines.
Definition: EdgeForce.hpp:25
ElasticMaterials
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Definition: ElasticMaterials.hpp:36
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:104
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:1221
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:1170
MoFEM::BasicMethod::postProcess
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: LoopMethods.cpp:161
MoFEM::BasicMethod::operator()
virtual MoFEMErrorCode operator()()
function is run for every finite element
Definition: LoopMethods.cpp:172
MoFEM::BasicMethod::preProcess
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: LoopMethods.cpp:150
MoFEM::BasicMethod::problemPtr
const Problem * problemPtr
raw pointer to problem
Definition: LoopMethods.hpp:259
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:31
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:92
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:85
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MoFEM::Ent_mi_tag
Definition: TagMultiIndices.hpp:34
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:378
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:31
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:114
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:233
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:254
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:219
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:31
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:36
MoFEM::SnesCtx::FEMethodsSequence
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:41
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:132
MoFEM::SnesMethod::snes_ctx
SNESContext snes_ctx
Definition: LoopMethods.hpp:128
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:103
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:33
NeumannForcesSurfaceComplexForLazy
NonLinear surface pressure element (obsolete implementation)
Definition: SurfacePressureComplexForLazy.hpp:27
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:25
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:40
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:247
PostProcStress
Definition: PostProcStresses.hpp:31
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:968
SphericalArcLengthControl
Implementation of spherical arc-length method.
Definition: ArcLengthTools.hpp:386

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.

Generated on Tue Jul 5 2022 22:52:24 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow