v0.15.0
Loading...
Searching...
No Matches
arc_length_nonlinear_elasticity.cpp
Go to the documentation of this file.
1/** \file arc_length_nonlinear_elasticity.cpp
2 * \ingroup nonlinear_elastic_elem
3 * \brief nonlinear elasticity (arc-length control)
4 *
5 * Solves nonlinear elastic problem. Using arc length control.
6 */
7
8
9
10static char help[] = "\
11 -my_file mesh file name\n\
12 -my_sr reduction of step size\n\
13 -my_ms maximal number of steps\n\n";
14
15#include <BasicFiniteElements.hpp>
16using namespace MoFEM;
17
18#include <ElasticMaterials.hpp>
19#include <NeoHookean.hpp>
20
21#include <SurfacePressureComplexForLazy.hpp>
22
23int main(int argc, char *argv[]) {
24
25 const string default_options = "-ksp_type fgmres \n"
26 "-pc_type lu \n"
27 "-pc_factor_mat_solver_type mumps \n"
28 "-mat_mumps_icntl_20 0 \n"
29 "-ksp_atol 1e-10 \n"
30 "-ksp_rtol 1e-10 \n"
31 "-snes_monitor \n"
32 "-snes_type newtonls \n"
33 "-snes_linesearch_type basic \n"
34 "-snes_max_it 100 \n"
35 "-snes_atol 1e-7 \n"
36 "-snes_rtol 1e-7 \n"
37 "-ts_monitor \n"
38 "-ts_type alpha \n";
39
40 string param_file = "param_file.petsc";
41 if (!static_cast<bool>(ifstream(param_file))) {
42 std::ofstream file(param_file.c_str(), std::ios::ate);
43 if (file.is_open()) {
44 file << default_options;
45 file.close();
46 }
47 }
48
49 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
50
51
52 try {
53
54 moab::Core mb_instance;
55 moab::Interface &moab = mb_instance;
56
57 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
58 auto moab_comm_wrap =
59 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
60 if (pcomm == NULL)
61 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
62
63 PetscBool flg = PETSC_TRUE;
64 char mesh_file_name[255];
65 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-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_NULLPTR, PETSC_NULLPTR, "-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_NULLPTR, PETSC_NULLPTR, "-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_NULLPTR, PETSC_NULLPTR, "-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
395 Range node_set;
396 for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field, "LoadPath", cit)) {
397 EntityHandle meshset = cit->getMeshset();
398 Range nodes;
399 CHKERR moab.get_entities_by_type(meshset, MBVERTEX, nodes, true);
400 MOAB_THROW(rval);
401 node_set.merge(nodes);
402 }
403 PetscPrintf(PETSC_COMM_WORLD, "Nb. nodes in load path: %u\n",
404 node_set.size());
405
406 SphericalArcLengthControl arc_method(arc_ctx);
407
408 double scaled_reference_load = 1;
409 double *scale_lhs = &(arc_ctx->getFieldData());
410 double *scale_rhs = &(scaled_reference_load);
411 NeumannForcesSurfaceComplexForLazy neumann_forces(
412 m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
413 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE &fe_neumann =
414 neumann_forces.getLoopSpatialFe();
415 if (linear) {
416 fe_neumann.typeOfForces = NeumannForcesSurfaceComplexForLazy::
417 MyTriangleSpatialFE::NONCONSERVATIVE;
418 }
419 fe_neumann.uSeF = true;
420 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
421 it)) {
422 CHKERR fe_neumann.addForce(it->getMeshsetId());
423 }
424 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
425 m_field, SIDESET | PRESSURESET, it)) {
426 CHKERR fe_neumann.addPressure(it->getMeshsetId());
427 }
428
429 boost::shared_ptr<FEMethod> my_dirichlet_bc =
430 boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
431 m_field, "SPATIAL_POSITION", Aij, D, F));
432 CHKERR m_field.get_problem("ELASTIC_MECHANICS",
433 &(my_dirichlet_bc->problemPtr));
434 CHKERR dynamic_cast<DirichletSpatialPositionsBc *>(my_dirichlet_bc.get())
435 ->iNitialize();
436
437 struct AssembleRhsVectors : public FEMethod {
438
439 boost::shared_ptr<ArcLengthCtx> arcPtr;
440 Range &nodeSet;
441
442 AssembleRhsVectors(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
443 Range &node_set)
444 : arcPtr(arc_ptr), nodeSet(node_set) {}
445
446 MoFEMErrorCode preProcess() {
447 MoFEMFunctionBeginHot;
448
449 // PetscAttachDebugger();
450 switch (snes_ctx) {
451 case CTX_SNESSETFUNCTION: {
452 CHKERR VecZeroEntries(snes_f);
453 CHKERR VecGhostUpdateBegin(snes_f, INSERT_VALUES, SCATTER_FORWARD);
454 CHKERR VecGhostUpdateEnd(snes_f, INSERT_VALUES, SCATTER_FORWARD);
455 CHKERR VecZeroEntries(arcPtr->F_lambda);
456 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, INSERT_VALUES,
457 SCATTER_FORWARD);
458 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, INSERT_VALUES,
459 SCATTER_FORWARD);
460 } break;
461 default:
462 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
463 }
464
465 MoFEMFunctionReturnHot(0);
466 }
467
468 MoFEMErrorCode postProcess() {
469 MoFEMFunctionBeginHot;
470 switch (snes_ctx) {
471 case CTX_SNESSETFUNCTION: {
472 // snes_f
473 CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
474 CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
475 CHKERR VecAssemblyBegin(snes_f);
476 CHKERR VecAssemblyEnd(snes_f);
477 } break;
478 default:
479 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
480 }
481 MoFEMFunctionReturnHot(0);
482 }
483
484 MoFEMErrorCode potsProcessLoadPath() {
485 MoFEMFunctionBeginHot;
486 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
487 problemPtr->getNumeredRowDofsPtr();
488 Range::iterator nit = nodeSet.begin();
489 for (; nit != nodeSet.end(); nit++) {
490 NumeredDofEntityByEnt::iterator dit, hi_dit;
491 dit = numered_dofs_rows->get<Ent_mi_tag>().lower_bound(*nit);
492 hi_dit = numered_dofs_rows->get<Ent_mi_tag>().upper_bound(*nit);
493 for (; dit != hi_dit; dit++) {
494 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e -> ", "LAMBDA", 0,
495 arcPtr->getFieldData());
496 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e\n",
497 dit->get()->getName().c_str(),
498 dit->get()->getDofCoeffIdx(),
499 dit->get()->getFieldData());
500 }
501 }
502 MoFEMFunctionReturnHot(0);
503 }
504 };
505
506 struct AddLambdaVectorToFInternal : public FEMethod {
507
508 boost::shared_ptr<ArcLengthCtx> arcPtr;
509 boost::shared_ptr<DirichletSpatialPositionsBc> bC;
510
511 AddLambdaVectorToFInternal(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
512 boost::shared_ptr<FEMethod> &bc)
513 : arcPtr(arc_ptr),
514 bC(boost::shared_ptr<DirichletSpatialPositionsBc>(
515 bc, dynamic_cast<DirichletSpatialPositionsBc *>(bc.get()))) {}
516
517 MoFEMErrorCode preProcess() {
518 MoFEMFunctionBeginHot;
519 MoFEMFunctionReturnHot(0);
520 }
521 MoFEMErrorCode operator()() {
522 MoFEMFunctionBeginHot;
523 MoFEMFunctionReturnHot(0);
524 }
525 MoFEMErrorCode postProcess() {
526 MoFEMFunctionBeginHot;
527 switch (snes_ctx) {
528 case CTX_SNESSETFUNCTION: {
529 // F_lambda
530 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES,
531 SCATTER_REVERSE);
532 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES,
533 SCATTER_REVERSE);
534 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
535 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
536 for (std::vector<int>::iterator vit = bC->dofsIndices.begin();
537 vit != bC->dofsIndices.end(); vit++) {
538 CHKERR VecSetValue(arcPtr->F_lambda, *vit, 0, INSERT_VALUES);
539 }
540 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
541 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
542 CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
543 PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n",
544 arcPtr->F_lambda2);
545 // add F_lambda
546 CHKERR VecAssemblyBegin(snes_f);
547 CHKERR VecAssemblyEnd(snes_f);
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 ///< do not very if element of given name exist when do loop over elements
572 CHKERR SnesMoFEMSetBehavior(snes, MF_ZERO);
573
574 PetscReal my_tol;
575 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_tol", &my_tol,
576 &flg);
577 if (flg == PETSC_TRUE) {
578 PetscReal atol, rtol, stol;
579 PetscInt maxit, maxf;
580 CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
581 atol = my_tol;
582 rtol = atol * 1e2;
583 CHKERR SNESSetTolerances(snes, atol, rtol, stol, maxit, maxf);
584 }
585
586 KSP ksp;
587 CHKERR SNESGetKSP(snes, &ksp);
588 PC pc;
589 CHKERR KSPGetPC(ksp, &pc);
590 boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
591 new PCArcLengthCtx(ShellAij, Aij, arc_ctx));
592 CHKERR PCSetType(pc, PCSHELL);
593 CHKERR PCShellSetContext(pc, pc_ctx.get());
594 CHKERR PCShellSetApply(pc, PCApplyArcLength);
595 CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
596
597 if (flg == PETSC_TRUE) {
598 PetscReal rtol, atol, dtol;
599 PetscInt maxits;
600 CHKERR KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
601 atol = my_tol * 1e-2;
602 rtol = atol * 1e-2;
603 CHKERR KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
604 }
605
606 SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
607 snes_ctx.getComputeRhs();
608 snes_ctx.getPreProcComputeRhs().push_back(my_dirichlet_bc);
609 snes_ctx.getPreProcComputeRhs().push_back(&pre_post_method);
610 loops_to_do_Rhs.push_back(
611 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
612
613 loops_to_do_Rhs.push_back(
614 SnesCtx::PairNameFEMethodPtr("SPRING", fe_spring_rhs_ptr.get()));
615
616 // surface forces and pressures
617 loops_to_do_Rhs.push_back(
618 SnesCtx::PairNameFEMethodPtr("NEUMANN_FE", &fe_neumann));
619
620 // edge forces
621 boost::ptr_map<std::string, EdgeForce> edge_forces;
622 string fe_name_str = "FORCE_FE";
623 edge_forces.insert(fe_name_str, new EdgeForce(m_field));
624 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
625 it)) {
626 CHKERR edge_forces.at(fe_name_str)
627 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
628 }
629 for (boost::ptr_map<std::string, EdgeForce>::iterator eit =
630 edge_forces.begin();
631 eit != edge_forces.end(); eit++) {
632 loops_to_do_Rhs.push_back(
633 SnesCtx::PairNameFEMethodPtr(eit->first, &eit->second->getLoopFe()));
634 }
635
636 // nodal forces
637 boost::ptr_map<std::string, NodalForce> nodal_forces;
638 // string fe_name_str ="FORCE_FE";
639 nodal_forces.insert(fe_name_str, new NodalForce(m_field));
640 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
641 it)) {
642 CHKERR nodal_forces.at(fe_name_str)
643 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
644 }
645 for (boost::ptr_map<std::string, NodalForce>::iterator fit =
646 nodal_forces.begin();
647 fit != nodal_forces.end(); fit++) {
648 loops_to_do_Rhs.push_back(
649 SnesCtx::PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
650 }
651
652 // arc length
653 loops_to_do_Rhs.push_back(
654 SnesCtx::PairNameFEMethodPtr("NONE", &assemble_F_lambda));
655 loops_to_do_Rhs.push_back(
656 SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", &arc_method));
657 snes_ctx.getPostProcComputeRhs().push_back(&pre_post_method);
658 snes_ctx.getPostProcComputeRhs().push_back(my_dirichlet_bc);
659
660 SnesCtx::FEMethodsSequence &loops_to_do_Mat =
661 snes_ctx.getSetOperators();
662 snes_ctx.getPreProcSetOperators().push_back(my_dirichlet_bc);
663 loops_to_do_Mat.push_back(
664 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
665
666 loops_to_do_Mat.push_back(
667 SnesCtx::PairNameFEMethodPtr("SPRING", fe_spring_lhs_ptr.get()));
668
669 loops_to_do_Mat.push_back(
670 SnesCtx::PairNameFEMethodPtr("NEUMANN_FE", &fe_neumann));
671 loops_to_do_Mat.push_back(
672 SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", &arc_method));
673 snes_ctx.getPostProcSetOperators().push_back(my_dirichlet_bc);
674
675 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
676 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
677 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
678 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
679
680 PetscScalar step_size_reduction;
681 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_sr",
682 &step_size_reduction, &flg);
683 if (flg != PETSC_TRUE) {
684 step_size_reduction = 1.;
685 }
686
687 PetscInt max_steps;
688 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_ms", &max_steps,
689 &flg);
690 if (flg != PETSC_TRUE) {
691 max_steps = 5;
692 }
693
694 int its_d;
695 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_its_d", &its_d,
696 &flg);
697 if (flg != PETSC_TRUE) {
698 its_d = 4;
699 }
700 PetscScalar max_reduction = 10, min_reduction = 0.1;
701 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_max_step_reduction",
702 &max_reduction, &flg);
703 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_min_step_reduction",
704 &min_reduction, &flg);
705
706 double gamma = 0.5, reduction = 1;
707 // step = 1;
708 if (step == 1) {
709 step_size = step_size_reduction;
710 } else {
711 reduction = step_size_reduction;
712 step++;
713 }
714 double step_size0 = step_size;
715
716 if (step > 1) {
717 CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
718 "ELASTIC_MECHANICS", "SPATIAL_POSITION", "X0_SPATIAL_POSITION", COL,
719 arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
720 double x0_nrm;
721 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
722 CHKERR PetscPrintf(PETSC_COMM_WORLD,
723 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
724 arc_ctx->dLambda);
725 CHKERR arc_ctx->setAlphaBeta(1, 0);
726 } else {
727 CHKERR arc_ctx->setS(step_size);
728 CHKERR arc_ctx->setAlphaBeta(0, 1);
729 }
730
731 CHKERR SnesRhs(snes, D, F, &snes_ctx);
732
733 Vec D0, x00;
734 CHKERR VecDuplicate(D, &D0);
735 CHKERR VecDuplicate(arc_ctx->x0, &x00);
736 bool converged_state = false;
737
738 for (int jj = 0; step < max_steps; step++, jj++) {
739
740 CHKERR VecCopy(D, D0);
741 CHKERR VecCopy(arc_ctx->x0, x00);
742
743 if (step == 1) {
744
745 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load Step %D step_size = %6.4e\n",
746 step, step_size);
747 CHKERR arc_ctx->setS(step_size);
748 CHKERR arc_ctx->setAlphaBeta(0, 1);
749 CHKERR VecCopy(D, arc_ctx->x0);
750 double dlambda;
751 CHKERR arc_method.calculateInitDlambda(&dlambda);
752 CHKERR arc_method.setDlambdaToX(D, dlambda);
753
754 } else if (step == 2) {
755
756 CHKERR arc_ctx->setAlphaBeta(1, 0);
757 CHKERR arc_method.calculateDxAndDlambda(D);
758 step_size = std::sqrt(arc_method.calculateLambdaInt());
759 step_size0 = step_size;
760 CHKERR arc_ctx->setS(step_size);
761 double dlambda = arc_ctx->dLambda;
762 double dx_nrm;
763 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
764 CHKERR PetscPrintf(PETSC_COMM_WORLD,
765 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
766 "dx_nrm = %6.4e dx2 = %6.4e\n",
767 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
768 CHKERR VecCopy(D, arc_ctx->x0);
769 CHKERR VecAXPY(D, 1., arc_ctx->dx);
770 CHKERR arc_method.setDlambdaToX(D, dlambda);
771
772 } else {
773
774 if (jj == 0) {
775 step_size0 = step_size;
776 }
777
778 CHKERR arc_method.calculateDxAndDlambda(D);
779 step_size *= reduction;
780 if (step_size > max_reduction * step_size0) {
781 step_size = max_reduction * step_size0;
782 } else if (step_size < min_reduction * step_size0) {
783 step_size = min_reduction * step_size0;
784 }
785 CHKERR arc_ctx->setS(step_size);
786 double dlambda = reduction * arc_ctx->dLambda;
787 double dx_nrm;
788 CHKERR VecScale(arc_ctx->dx, reduction);
789 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
790 CHKERR PetscPrintf(PETSC_COMM_WORLD,
791 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
792 "dx_nrm = %6.4e dx2 = %6.4e\n",
793 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
794 CHKERR VecCopy(D, arc_ctx->x0);
795 CHKERR VecAXPY(D, 1., arc_ctx->dx);
796 CHKERR arc_method.setDlambdaToX(D, dlambda);
797 }
798
799 CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
800 int its;
801 CHKERR SNESGetIterationNumber(snes, &its);
802 CHKERR PetscPrintf(PETSC_COMM_WORLD, "number of Newton iterations = %D\n",
803 its);
804
805 SNESConvergedReason reason;
806 CHKERR SNESGetConvergedReason(snes, &reason);
807 if (reason < 0) {
808
809 CHKERR VecCopy(D0, D);
810 CHKERR VecCopy(x00, arc_ctx->x0);
811
812 double x0_nrm;
813 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
814 CHKERR PetscPrintf(PETSC_COMM_WORLD,
815 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
816 arc_ctx->dLambda);
817 CHKERR arc_ctx->setAlphaBeta(1, 0);
818
819 reduction = 0.1;
820 converged_state = false;
821
822 continue;
823
824 } else {
825
826 if (step > 1 && converged_state) {
827
828 reduction = pow((double)its_d / (double)(its + 1), gamma);
829 if (step_size >= max_reduction * step_size0 && reduction > 1) {
830 reduction = 1;
831 } else if (step_size <= min_reduction * step_size0 && reduction < 1) {
832 reduction = 1;
833 }
834 CHKERR PetscPrintf(PETSC_COMM_WORLD, "reduction step_size = %6.4e\n",
835 reduction);
836 }
837
838 // Save data on mesh
839 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
840 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
841 CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
842 "ELASTIC_MECHANICS", "SPATIAL_POSITION", "X0_SPATIAL_POSITION", COL,
843 arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
844 converged_state = true;
845 }
846
847 if (step % 1 == 0) {
848 // Save restart file
849 // #ifdef MOAB_HDF5_PARALLEL
850 // std::ostringstream sss;
851 // sss << "restart_" << step << ".h5m";
852 // CHKERR
853 // moab.write_file(sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART");
854 // #else
855 // #warning "No parallel HDF5, no writing restart file"
856 // #endif
857 // Save data on mesh
858 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
859 post_proc);
860 std::ostringstream o1;
861 o1 << "out_" << step << ".h5m";
862 CHKERR post_proc.writeFile(o1.str().c_str());
863 }
864
865 CHKERR pre_post_method.potsProcessLoadPath();
866 }
867
868 CHKERR VecDestroy(&D0);
869 CHKERR VecDestroy(&x00);
870
871 // detroy matrices
872 CHKERR VecDestroy(&F);
873 CHKERR VecDestroy(&D);
874 CHKERR MatDestroy(&Aij);
875 CHKERR MatDestroy(&ShellAij);
876 CHKERR SNESDestroy(&snes);
877 }
878 CATCH_ERRORS;
879
880 MoFEM::Core::Finalize();
881
882 return 0;
883}
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
ElasticMaterials.hpp
Elastic materials.
NeoHookean.hpp
Implementation of Neo-Hookean elastic material.
main
int main()
Definition adol-c_atom.cpp:46
help
static char help[]
Definition arc_length_nonlinear_elasticity.cpp:10
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 nme:nme847.
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:557
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
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: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:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
F
@ F
Definition free_surface.cpp:396
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_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::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_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::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::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.
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition PostProcOnRefMesh.hpp:153
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
Definition PostProcOnRefMesh.hpp:195
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
Definition PostProcOnRefMesh.hpp:253
MoFEM::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:2344
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition ProblemsManager.cpp:86
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition ProblemsManager.cpp:1679
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:2163
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
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
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
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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:483
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:223
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::SnesMoFEMSetBehavior
MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
Definition SnesCtx.cpp:578
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:678
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
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
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:170
MoFEM::BasicMethod::operator()
virtual MoFEMErrorCode operator()()
function is run for every finite element
Definition LoopMethods.cpp:181
MoFEM::BasicMethod::preProcess
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition LoopMethods.cpp:159
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:118
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:404
MoFEM::FieldBlas
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition ForcesAndSourcesCore.hpp:83
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:315
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition MeshsetsManager.cpp:336
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition MeshsetsManager.cpp:301
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:18
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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
NonlinearElasticElement::getLoopFeLhs
MyVolumeFE & getLoopFeLhs()
get lhs volume element
Definition NonLinearElasticElement.hpp:68
NonlinearElasticElement::addElement
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition NonLinearElasticElement.cpp:1119
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition NonLinearElasticElement.hpp:100
NonlinearElasticElement::getLoopFeRhs
MyVolumeFE & getLoopFeRhs()
get rhs volume element
Definition NonLinearElasticElement.hpp:66
NonlinearElasticElement::getLoopFeEnergy
MyVolumeFE & getLoopFeEnergy()
get energy fe
Definition NonLinearElasticElement.hpp:71
NonlinearElasticElement::setOperators
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Definition NonLinearElasticElement.cpp:1152
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition ArcLengthTools.hpp:235
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition PostProcOnRefMesh.hpp:125
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition PostProcOnRefMesh.hpp:301
SphericalArcLengthControl
Implementation of spherical arc-length method.
Definition ArcLengthTools.hpp:374
SphericalArcLengthControl::calculateDxAndDlambda
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition ArcLengthTools.cpp:742
SphericalArcLengthControl::calculateLambdaInt
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
Definition ArcLengthTools.cpp:685
SphericalArcLengthControl::calculateInitDlambda
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
Definition ArcLengthTools.cpp:769
SphericalArcLengthControl::setDlambdaToX
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
Definition ArcLengthTools.cpp:783
Generated by Doxygen 1.12.0 and hosted at University of Glasgow