v0.16.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 AddHOOps<3, 3, 3>::add(elastic.getLoopFeRhs().getOpPtrVector(), {H1},
297 "MESH_NODE_POSITIONS");
298 CHKERR AddHOOps<3, 3, 3>::add(elastic.getLoopFeLhs().getOpPtrVector(), {H1},
299 "MESH_NODE_POSITIONS");
300 CHKERR AddHOOps<3, 3, 3>::add(
301 elastic.getLoopFeEnergy().getOpPtrVector(), {H1},
302 "MESH_NODE_POSITIONS");
303 CHKERR elastic.setOperators("SPATIAL_POSITION");
304
305 // post_processing
306 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
307 m_field);
308 CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {H1},
309 "MESH_NODE_POSITIONS");
310 auto spatial_pos_ptr = boost::make_shared<MatrixDouble>();
311 auto mesh_pos_ptr = boost::make_shared<MatrixDouble>();
312 auto spatial_pos_grad_ptr = boost::make_shared<MatrixDouble>();
313 post_proc.getOpPtrVector().push_back(
314 new OpCalculateVectorFieldValues<3>("SPATIAL_POSITION",
315 spatial_pos_ptr));
316 post_proc.getOpPtrVector().push_back(
317 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
318 mesh_pos_ptr));
319 post_proc.getOpPtrVector().push_back(
320 new OpCalculateVectorFieldGradient<3, 3>("SPATIAL_POSITION",
321 spatial_pos_grad_ptr));
322 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
323 elastic.setOfBlocks.begin();
324 for (; sit != elastic.setOfBlocks.end(); sit++) {
325 post_proc.getOpPtrVector().push_back(new PostProcStress(
326 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
327 post_proc.getPostProcElements(), "SPATIAL_POSITION", sit->second,
328 spatial_pos_ptr, mesh_pos_ptr, spatial_pos_grad_ptr));
329 }
330 using OpPPMap = OpPostProcMapInMoab<3, 3>;
331 post_proc.getOpPtrVector().push_back(new OpPPMap(
332 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), {},
333 {{"SPATIAL_POSITION", spatial_pos_ptr},
334 {"MESH_NODE_POSITIONS", mesh_pos_ptr}},
335 {{"SPATIAL_POSITION_GRAD", spatial_pos_grad_ptr}}, {}));
336
337 // build field
338 CHKERR m_field.build_fields();
339 if (step == 1) {
340 // 10 node tets
341 Projection10NodeCoordsOnField ent_method_material(m_field,
342 "MESH_NODE_POSITIONS");
343 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material, 0);
344 CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBVERTEX,
345 "SPATIAL_POSITION");
346 CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBEDGE,
347 "SPATIAL_POSITION");
348 CHKERR m_field.getInterface<FieldBlas>()->fieldAxpy(
349 1., "MESH_NODE_POSITIONS", "SPATIAL_POSITION");
350 CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBTRI,
351 "SPATIAL_POSITION");
352 CHKERR m_field.getInterface<FieldBlas>()->setField(0, MBTET,
353 "SPATIAL_POSITION");
354 }
355
356 // build finite elements
357 CHKERR m_field.build_finite_elements();
358
359 // build adjacencies
360 CHKERR m_field.build_adjacencies(problem_bit_level);
361
362 ProblemsManager *prb_mng_ptr;
363 CHKERR m_field.getInterface(prb_mng_ptr);
364 // build database
365 if (is_partitioned) {
366 SETERRQ(PETSC_COMM_SELF, 1,
367 "Not implemented, problem with arc-length force multiplayer");
368 } else {
369 CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
370 CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
371 CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS");
372 }
373 CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
374
375 // print bcs
376 MeshsetsManager *mmanager_ptr;
377 CHKERR m_field.getInterface(mmanager_ptr);
378 CHKERR mmanager_ptr->printDisplacementSet();
379 CHKERR mmanager_ptr->printForceSet();
380 // print block sets with materials
381 CHKERR mmanager_ptr->printMaterialsSet();
382
383 // create matrices
384 Vec F;
385 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
386 "ELASTIC_MECHANICS", COL, &F);
387 Vec D;
388 CHKERR VecDuplicate(F, &D);
389 Mat Aij;
390 CHKERR m_field.getInterface<MatrixManager>()
391 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
392 &Aij);
393
394 boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
395 new ArcLengthCtx(m_field, "ELASTIC_MECHANICS"));
396
397 PetscInt M, N;
398 CHKERR MatGetSize(Aij, &M, &N);
399 PetscInt m, n;
400 CHKERR MatGetLocalSize(Aij, &m, &n);
401 boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
402 new ArcLengthMatShell(Aij, arc_ctx, "ELASTIC_MECHANICS"));
403
404 Mat ShellAij;
405 CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, mat_ctx.get(),
406 &ShellAij);
407 CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
408 (void (*)(void))ArcLengthMatMultShellOp);
409
410 ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
411
412 Range node_set;
413 for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field, "LoadPath", cit)) {
414 EntityHandle meshset = cit->getMeshset();
415 Range nodes;
416 CHKERR moab.get_entities_by_type(meshset, MBVERTEX, nodes, true);
417 MOAB_THROW(rval);
418 node_set.merge(nodes);
419 }
420 PetscPrintf(PETSC_COMM_WORLD, "Nb. nodes in load path: %u\n",
421 node_set.size());
422
423 SphericalArcLengthControl arc_method(arc_ctx);
424
425 double scaled_reference_load = 1;
426 double *scale_lhs = &(arc_ctx->getFieldData());
427 double *scale_rhs = &(scaled_reference_load);
428 NeumannForcesSurfaceComplexForLazy neumann_forces(
429 m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
430 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE &fe_neumann =
431 neumann_forces.getLoopSpatialFe();
432 if (linear) {
433 fe_neumann.typeOfForces = NeumannForcesSurfaceComplexForLazy::
434 MyTriangleSpatialFE::NONCONSERVATIVE;
435 }
436 fe_neumann.uSeF = true;
437 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
438 it)) {
439 CHKERR fe_neumann.addForce(it->getMeshsetId());
440 }
441 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
442 m_field, SIDESET | PRESSURESET, it)) {
443 CHKERR fe_neumann.addPressure(it->getMeshsetId());
444 }
445
446 boost::shared_ptr<FEMethod> my_dirichlet_bc =
447 boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
448 m_field, "SPATIAL_POSITION", Aij, D, F));
449 CHKERR m_field.get_problem("ELASTIC_MECHANICS",
450 &(my_dirichlet_bc->problemPtr));
451 CHKERR dynamic_cast<DirichletSpatialPositionsBc *>(my_dirichlet_bc.get())
452 ->iNitialize();
453
454 struct AssembleRhsVectors : public FEMethod {
455
456 boost::shared_ptr<ArcLengthCtx> arcPtr;
457 Range &nodeSet;
458
459 AssembleRhsVectors(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
460 Range &node_set)
461 : arcPtr(arc_ptr), nodeSet(node_set) {}
462
463 MoFEMErrorCode preProcess() {
464 MoFEMFunctionBeginHot;
465
466 // PetscAttachDebugger();
467 switch (snes_ctx) {
468 case CTX_SNESSETFUNCTION: {
469 CHKERR VecZeroEntries(snes_f);
470 CHKERR VecGhostUpdateBegin(snes_f, INSERT_VALUES, SCATTER_FORWARD);
471 CHKERR VecGhostUpdateEnd(snes_f, INSERT_VALUES, SCATTER_FORWARD);
472 CHKERR VecZeroEntries(arcPtr->F_lambda);
473 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, INSERT_VALUES,
474 SCATTER_FORWARD);
475 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, INSERT_VALUES,
476 SCATTER_FORWARD);
477 } break;
478 default:
479 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
480 }
481
482 MoFEMFunctionReturnHot(0);
483 }
484
485 MoFEMErrorCode postProcess() {
486 MoFEMFunctionBeginHot;
487 switch (snes_ctx) {
488 case CTX_SNESSETFUNCTION: {
489 // snes_f
490 CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
491 CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
492 CHKERR VecAssemblyBegin(snes_f);
493 CHKERR VecAssemblyEnd(snes_f);
494 } break;
495 default:
496 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
497 }
498 MoFEMFunctionReturnHot(0);
499 }
500
501 MoFEMErrorCode potsProcessLoadPath() {
502 MoFEMFunctionBeginHot;
503 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
504 problemPtr->getNumeredRowDofsPtr();
505 Range::iterator nit = nodeSet.begin();
506 for (; nit != nodeSet.end(); nit++) {
507 NumeredDofEntityByEnt::iterator dit, hi_dit;
508 dit = numered_dofs_rows->get<Ent_mi_tag>().lower_bound(*nit);
509 hi_dit = numered_dofs_rows->get<Ent_mi_tag>().upper_bound(*nit);
510 for (; dit != hi_dit; dit++) {
511 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e -> ", "LAMBDA", 0,
512 arcPtr->getFieldData());
513 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e\n",
514 dit->get()->getName().c_str(),
515 dit->get()->getDofCoeffIdx(),
516 dit->get()->getFieldData());
517 }
518 }
519 MoFEMFunctionReturnHot(0);
520 }
521 };
522
523 struct AddLambdaVectorToFInternal : public FEMethod {
524
525 boost::shared_ptr<ArcLengthCtx> arcPtr;
526 boost::shared_ptr<DirichletSpatialPositionsBc> bC;
527
528 AddLambdaVectorToFInternal(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
529 boost::shared_ptr<FEMethod> &bc)
530 : arcPtr(arc_ptr),
531 bC(boost::shared_ptr<DirichletSpatialPositionsBc>(
532 bc, dynamic_cast<DirichletSpatialPositionsBc *>(bc.get()))) {}
533
534 MoFEMErrorCode preProcess() {
535 MoFEMFunctionBeginHot;
536 MoFEMFunctionReturnHot(0);
537 }
538 MoFEMErrorCode operator()() {
539 MoFEMFunctionBeginHot;
540 MoFEMFunctionReturnHot(0);
541 }
542 MoFEMErrorCode postProcess() {
543 MoFEMFunctionBeginHot;
544 switch (snes_ctx) {
545 case CTX_SNESSETFUNCTION: {
546 // F_lambda
547 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES,
548 SCATTER_REVERSE);
549 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES,
550 SCATTER_REVERSE);
551 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
552 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
553 for (std::vector<int>::iterator vit = bC->dofsIndices.begin();
554 vit != bC->dofsIndices.end(); vit++) {
555 CHKERR VecSetValue(arcPtr->F_lambda, *vit, 0, INSERT_VALUES);
556 }
557 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
558 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
559 CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
560 PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n",
561 arcPtr->F_lambda2);
562 // add F_lambda
563 CHKERR VecAssemblyBegin(snes_f);
564 CHKERR VecAssemblyEnd(snes_f);
565 CHKERR VecAXPY(snes_f, arcPtr->getFieldData(), arcPtr->F_lambda);
566 PetscPrintf(PETSC_COMM_WORLD, "\tlambda = %6.4e\n",
567 arcPtr->getFieldData());
568 double fnorm;
569 CHKERR VecNorm(snes_f, NORM_2, &fnorm);
570 PetscPrintf(PETSC_COMM_WORLD, "\tfnorm = %6.4e\n", fnorm);
571 } break;
572 default:
573 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
574 }
575 MoFEMFunctionReturnHot(0);
576 }
577 };
578
579 AssembleRhsVectors pre_post_method(arc_ctx, node_set);
580 AddLambdaVectorToFInternal assemble_F_lambda(arc_ctx, my_dirichlet_bc);
581
582 SNES snes;
583 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
584 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
585 CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
586 CHKERR SNESSetJacobian(snes, ShellAij, Aij, SnesMat, &snes_ctx);
587 CHKERR SNESSetFromOptions(snes);
588 ///< do not very if element of given name exist when do loop over elements
589 CHKERR SnesMoFEMSetBehavior(snes, MF_ZERO);
590
591 PetscReal my_tol;
592 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_tol", &my_tol,
593 &flg);
594 if (flg == PETSC_TRUE) {
595 PetscReal atol, rtol, stol;
596 PetscInt maxit, maxf;
597 CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
598 atol = my_tol;
599 rtol = atol * 1e2;
600 CHKERR SNESSetTolerances(snes, atol, rtol, stol, maxit, maxf);
601 }
602
603 KSP ksp;
604 CHKERR SNESGetKSP(snes, &ksp);
605 PC pc;
606 CHKERR KSPGetPC(ksp, &pc);
607 boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
608 new PCArcLengthCtx(ShellAij, Aij, arc_ctx));
609 CHKERR PCSetType(pc, PCSHELL);
610 CHKERR PCShellSetContext(pc, pc_ctx.get());
611 CHKERR PCShellSetApply(pc, PCApplyArcLength);
612 CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
613
614 if (flg == PETSC_TRUE) {
615 PetscReal rtol, atol, dtol;
616 PetscInt maxits;
617 CHKERR KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
618 atol = my_tol * 1e-2;
619 rtol = atol * 1e-2;
620 CHKERR KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
621 }
622
623 SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
624 snes_ctx.getComputeRhs();
625 snes_ctx.getPreProcComputeRhs().push_back(my_dirichlet_bc);
626 snes_ctx.getPreProcComputeRhs().push_back(&pre_post_method);
627 loops_to_do_Rhs.push_back(
628 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
629
630 loops_to_do_Rhs.push_back(
631 SnesCtx::PairNameFEMethodPtr("SPRING", fe_spring_rhs_ptr.get()));
632
633 // surface forces and pressures
634 loops_to_do_Rhs.push_back(
635 SnesCtx::PairNameFEMethodPtr("NEUMANN_FE", &fe_neumann));
636
637 // edge forces
638 boost::ptr_map<std::string, EdgeForce> edge_forces;
639 string fe_name_str = "FORCE_FE";
640 edge_forces.insert(fe_name_str, new EdgeForce(m_field));
641 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
642 it)) {
643 CHKERR edge_forces.at(fe_name_str)
644 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
645 }
646 for (boost::ptr_map<std::string, EdgeForce>::iterator eit =
647 edge_forces.begin();
648 eit != edge_forces.end(); eit++) {
649 loops_to_do_Rhs.push_back(
650 SnesCtx::PairNameFEMethodPtr(eit->first, &eit->second->getLoopFe()));
651 }
652
653 // nodal forces
654 boost::ptr_map<std::string, NodalForce> nodal_forces;
655 // string fe_name_str ="FORCE_FE";
656 nodal_forces.insert(fe_name_str, new NodalForce(m_field));
657 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
658 it)) {
659 CHKERR nodal_forces.at(fe_name_str)
660 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
661 }
662 for (boost::ptr_map<std::string, NodalForce>::iterator fit =
663 nodal_forces.begin();
664 fit != nodal_forces.end(); fit++) {
665 loops_to_do_Rhs.push_back(
666 SnesCtx::PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
667 }
668
669 // arc length
670 loops_to_do_Rhs.push_back(
671 SnesCtx::PairNameFEMethodPtr("NONE", &assemble_F_lambda));
672 loops_to_do_Rhs.push_back(
673 SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", &arc_method));
674 snes_ctx.getPostProcComputeRhs().push_back(&pre_post_method);
675 snes_ctx.getPostProcComputeRhs().push_back(my_dirichlet_bc);
676
677 SnesCtx::FEMethodsSequence &loops_to_do_Mat =
678 snes_ctx.getSetOperators();
679 snes_ctx.getPreProcSetOperators().push_back(my_dirichlet_bc);
680 loops_to_do_Mat.push_back(
681 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
682
683 loops_to_do_Mat.push_back(
684 SnesCtx::PairNameFEMethodPtr("SPRING", fe_spring_lhs_ptr.get()));
685
686 loops_to_do_Mat.push_back(
687 SnesCtx::PairNameFEMethodPtr("NEUMANN_FE", &fe_neumann));
688 loops_to_do_Mat.push_back(
689 SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", &arc_method));
690 snes_ctx.getPostProcSetOperators().push_back(my_dirichlet_bc);
691
692 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
693 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
694 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
695 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
696
697 PetscScalar step_size_reduction;
698 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_sr",
699 &step_size_reduction, &flg);
700 if (flg != PETSC_TRUE) {
701 step_size_reduction = 1.;
702 }
703
704 PetscInt max_steps;
705 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_ms", &max_steps,
706 &flg);
707 if (flg != PETSC_TRUE) {
708 max_steps = 5;
709 }
710
711 int its_d;
712 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_its_d", &its_d,
713 &flg);
714 if (flg != PETSC_TRUE) {
715 its_d = 4;
716 }
717 PetscScalar max_reduction = 10, min_reduction = 0.1;
718 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_max_step_reduction",
719 &max_reduction, &flg);
720 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_min_step_reduction",
721 &min_reduction, &flg);
722
723 double gamma = 0.5, reduction = 1;
724 // step = 1;
725 if (step == 1) {
726 step_size = step_size_reduction;
727 } else {
728 reduction = step_size_reduction;
729 step++;
730 }
731 double step_size0 = step_size;
732
733 if (step > 1) {
734 CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
735 "ELASTIC_MECHANICS", "SPATIAL_POSITION", "X0_SPATIAL_POSITION", COL,
736 arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
737 double x0_nrm;
738 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
739 CHKERR PetscPrintf(PETSC_COMM_WORLD,
740 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
741 arc_ctx->dLambda);
742 CHKERR arc_ctx->setAlphaBeta(1, 0);
743 } else {
744 CHKERR arc_ctx->setS(step_size);
745 CHKERR arc_ctx->setAlphaBeta(0, 1);
746 }
747
748 CHKERR SnesRhs(snes, D, F, &snes_ctx);
749
750 Vec D0, x00;
751 CHKERR VecDuplicate(D, &D0);
752 CHKERR VecDuplicate(arc_ctx->x0, &x00);
753 bool converged_state = false;
754
755 for (int jj = 0; step < max_steps; step++, jj++) {
756
757 CHKERR VecCopy(D, D0);
758 CHKERR VecCopy(arc_ctx->x0, x00);
759
760 if (step == 1) {
761
762 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load Step %D step_size = %6.4e\n",
763 step, step_size);
764 CHKERR arc_ctx->setS(step_size);
765 CHKERR arc_ctx->setAlphaBeta(0, 1);
766 CHKERR VecCopy(D, arc_ctx->x0);
767 double dlambda;
768 CHKERR arc_method.calculateInitDlambda(&dlambda);
769 CHKERR arc_method.setDlambdaToX(D, dlambda);
770
771 } else if (step == 2) {
772
773 CHKERR arc_ctx->setAlphaBeta(1, 0);
774 CHKERR arc_method.calculateDxAndDlambda(D);
775 step_size = std::sqrt(arc_method.calculateLambdaInt());
776 step_size0 = step_size;
777 CHKERR arc_ctx->setS(step_size);
778 double dlambda = arc_ctx->dLambda;
779 double dx_nrm;
780 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
781 CHKERR PetscPrintf(PETSC_COMM_WORLD,
782 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
783 "dx_nrm = %6.4e dx2 = %6.4e\n",
784 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
785 CHKERR VecCopy(D, arc_ctx->x0);
786 CHKERR VecAXPY(D, 1., arc_ctx->dx);
787 CHKERR arc_method.setDlambdaToX(D, dlambda);
788
789 } else {
790
791 if (jj == 0) {
792 step_size0 = step_size;
793 }
794
795 CHKERR arc_method.calculateDxAndDlambda(D);
796 step_size *= reduction;
797 if (step_size > max_reduction * step_size0) {
798 step_size = max_reduction * step_size0;
799 } else if (step_size < min_reduction * step_size0) {
800 step_size = min_reduction * step_size0;
801 }
802 CHKERR arc_ctx->setS(step_size);
803 double dlambda = reduction * arc_ctx->dLambda;
804 double dx_nrm;
805 CHKERR VecScale(arc_ctx->dx, reduction);
806 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
807 CHKERR PetscPrintf(PETSC_COMM_WORLD,
808 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
809 "dx_nrm = %6.4e dx2 = %6.4e\n",
810 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
811 CHKERR VecCopy(D, arc_ctx->x0);
812 CHKERR VecAXPY(D, 1., arc_ctx->dx);
813 CHKERR arc_method.setDlambdaToX(D, dlambda);
814 }
815
816 CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
817 int its;
818 CHKERR SNESGetIterationNumber(snes, &its);
819 CHKERR PetscPrintf(PETSC_COMM_WORLD, "number of Newton iterations = %D\n",
820 its);
821
822 SNESConvergedReason reason;
823 CHKERR SNESGetConvergedReason(snes, &reason);
824 if (reason < 0) {
825
826 CHKERR VecCopy(D0, D);
827 CHKERR VecCopy(x00, arc_ctx->x0);
828
829 double x0_nrm;
830 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
831 CHKERR PetscPrintf(PETSC_COMM_WORLD,
832 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
833 arc_ctx->dLambda);
834 CHKERR arc_ctx->setAlphaBeta(1, 0);
835
836 reduction = 0.1;
837 converged_state = false;
838
839 continue;
840
841 } else {
842
843 if (step > 1 && converged_state) {
844
845 reduction = pow((double)its_d / (double)(its + 1), gamma);
846 if (step_size >= max_reduction * step_size0 && reduction > 1) {
847 reduction = 1;
848 } else if (step_size <= min_reduction * step_size0 && reduction < 1) {
849 reduction = 1;
850 }
851 CHKERR PetscPrintf(PETSC_COMM_WORLD, "reduction step_size = %6.4e\n",
852 reduction);
853 }
854
855 // Save data on mesh
856 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
857 "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
858 CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
859 "ELASTIC_MECHANICS", "SPATIAL_POSITION", "X0_SPATIAL_POSITION", COL,
860 arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
861 converged_state = true;
862 }
863
864 if (step % 1 == 0) {
865 // Save restart file
866 // #ifdef MOAB_HDF5_PARALLEL
867 // std::ostringstream sss;
868 // sss << "restart_" << step << ".h5m";
869 // CHKERR
870 // moab.write_file(sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART");
871 // #else
872 // #warning "No parallel HDF5, no writing restart file"
873 // #endif
874 // Save data on mesh
875 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
876 post_proc);
877 std::ostringstream o1;
878 o1 << "out_" << step << ".h5m";
879 CHKERR post_proc.writeFile(o1.str().c_str());
880 }
881
882 CHKERR pre_post_method.potsProcessLoadPath();
883 }
884
885 CHKERR VecDestroy(&D0);
886 CHKERR VecDestroy(&x00);
887
888 // detroy matrices
889 CHKERR VecDestroy(&F);
890 CHKERR VecDestroy(&D);
891 CHKERR MatDestroy(&Aij);
892 CHKERR MatDestroy(&ShellAij);
893 CHKERR SNESDestroy(&snes);
894 }
895 CATCH_ERRORS;
896
897 MoFEM::Core::Finalize();
898
899 return 0;
900}
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 .
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:38
F
@ F
Definition free_surface.cpp:630
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::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::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::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
Print meshsets with force boundary conditions.
Definition MeshsetsManager.cpp:425
_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:117
_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:72
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
Print meshsets with material properties.
Definition MeshsetsManager.cpp:446
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
Print meshsets with displacement boundary conditions.
Definition MeshsetsManager.cpp:411
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:45
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:39
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:491
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:227
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:586
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:40
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:63
ArcLengthMatShell
shell matrix for arc-length method
Definition ArcLengthTools.hpp:203
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition DirichletBC.hpp:210
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::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition HODataOperators.hpp:642
MoFEM::BasicMethod::postProcess
virtual MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Definition LoopMethods.cpp:162
MoFEM::BasicMethod::operator()
virtual MoFEMErrorCode operator()()
Main operator function executed for each loop iteration.
Definition LoopMethods.cpp:173
MoFEM::BasicMethod::preProcess
virtual MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
Definition LoopMethods.cpp:151
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::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:83
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:68
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
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:657
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:121
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:208
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1711
MoFEM::OpCalculateVectorFieldValues
Specialization for MatrixDouble vector field values calculation.
Definition UserDataOperators.hpp:589
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:745
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition ProblemsManager.hpp:133
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::SnesMethod::CTX_SNESSETFUNCTION
@ CTX_SNESSETFUNCTION
Setting up nonlinear function evaluation.
Definition LoopMethods.hpp:213
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:1118
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:1151
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition ArcLengthTools.hpp:233
SphericalArcLengthControl
Implementation of spherical arc-length method.
Definition ArcLengthTools.hpp:372
SphericalArcLengthControl::calculateDxAndDlambda
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition ArcLengthTools.cpp:738
SphericalArcLengthControl::calculateLambdaInt
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
Definition ArcLengthTools.cpp:680
SphericalArcLengthControl::calculateInitDlambda
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
Definition ArcLengthTools.cpp:766
SphericalArcLengthControl::setDlambdaToX
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
Definition ArcLengthTools.cpp:780
PCApplyArcLength
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Definition ArcLengthTools.cpp:229
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition ArcLengthTools.cpp:174
PCSetupArcLength
MoFEMErrorCode PCSetupArcLength(PC pc)
Definition ArcLengthTools.cpp:316
Generated by Doxygen 1.9.8 and hosted at University of Glasgow