v0.14.0
Searching...
No Matches
elasticity.cpp
Go to the documentation of this file.
1/** \file elasticity.cpp
2 * \ingroup nonlinear_elastic_elem
3 * \example elasticity.cpp
4
5 The example shows how to solve the linear elastic problem. An example can read
6 file with temperature field, then thermal stresses are included.
7
8 What example can do:
9 - take into account temperature field, i.e. calculate thermal stresses and
10deformation
11 - stationary and time depend field is considered
12 - take into account gravitational body forces
13 - take in account fluid pressure
14 - can work with higher order geometry definition
15 - works on distributed meshes
16 - multi-grid solver where each grid level is approximation order level
17 - each mesh block can have different material parameters and approximation
18order
19
20See example how code can be used \cite jordi:2017,
21 \image html SquelaDamExampleByJordi.png "Example what you can do with this
22code. Analysis of the arch dam of Susqueda, located in Catalonia (Spain)"
23width=800px
24
25 This is an example of application code; it does not show how elements are
26implemented. Example presents how to:
28 - set-up problem
29 - run finite elements on the problem
30 - assemble matrices and vectors
31 - solve the linear system of equations
32 - save results
33
34
35 If you like to see how to implement finite elements, material, are other parts
36of the code, look here;
37 - Hooke material, see \ref Hooke
38 - Thermal-stress assembly, see \ref ThermalElement
39 - Body forces element, see \ref BodyForceConstantField
40 - Fluid pressure element, see \ref FluidPressure
41 - The general implementation of an element for arbitrary Lagrangian-Eulerian
42 formulation for a nonlinear elastic problem is here \ref
43 NonlinearElasticElement. Here we limit ourselves to Hooke equation and fix
44 mesh, so the problem becomes linear. Not that elastic element is implemented
45 with automatic differentiation.
46
47*/
48
49
50
52using namespace MoFEM;
53
54#include <Hooke.hpp>
55using namespace boost::numeric;
56
57static char help[] = "-my_block_config set block data\n"
58 "-my_order approximation order\n"
59 "-my_is_partitioned set if mesh is partitioned\n"
60 "\n";
61
63 int oRder;
64 int iD;
65 double yOung;
66 double pOisson;
67 double initTemp;
68
70
71 BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2), initTemp(0) {}
72};
73
77
78/// Set integration rule
79struct VolRule {
80 int operator()(int, int, int order) const { return 2 * order; }
81};
82
84
86 FatPrismElementForcesAndSourcesCore;
89};
90
91int PrismFE::getRuleTrianglesOnly(int order) { return 2 * order; };
92int PrismFE::getRuleThroughThickness(int order) { return 2 * order; };
93
95
96int main(int argc, char *argv[]) {
97
98 const string default_options = "-ksp_type gmres \n"
99 "-pc_type lu \n"
100 "-pc_factor_mat_solver_type mumps \n"
101 "-mat_mumps_icntl_20 0 \n"
102 "-ksp_monitor \n"
103 "-snes_type newtonls \n"
104 "-snes_linesearch_type basic \n"
105 "-snes_atol 1e-8 \n"
106 "-snes_rtol 1e-8 \n"
107 "-snes_monitor \n"
108 "-ts_monitor \n"
109 "-ts_type beuler \n";
110
111 string param_file = "param_file.petsc";
112 if (!static_cast<bool>(ifstream(param_file))) {
113 std::ofstream file(param_file.c_str(), std::ios::ate);
114 if (file.is_open()) {
115 file << default_options;
116 file.close();
117 }
118 }
119
120 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
121 auto core_log = logging::core::get();
124 LogManager::setLog("ELASTIC");
125 MOFEM_LOG_TAG("ELASTIC", "elasticity")
126
129 LogManager::setLog("ELASTIC_SYNC");
130 MOFEM_LOG_TAG("ELASTIC_SYNC", "elastic_sync");
131
132 try {
133
134 PetscBool flg_block_config, flg_file;
135 char mesh_file_name[255];
136 char block_config_file[255];
137 PetscInt test_nb = 0;
138 PetscInt order = 2;
139 PetscBool is_partitioned = PETSC_FALSE;
140 PetscBool is_calculating_frequency = PETSC_FALSE;
141 PetscBool is_post_proc_volume = PETSC_TRUE;
142
143 // Select base
144 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
145 const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
146 "jacobi"};
147 PetscInt choice_base_value = LOBATTO;
148
149 // Read options from command line
150 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
151 CHKERR(ierr);
152 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
153 mesh_file_name, 255, &flg_file);
154
155 CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
156 order, &order, PETSC_NULL);
157
158 CHKERR PetscOptionsEList("-base", "approximation base", "", list_bases,
159 LASBASETOP, list_bases[choice_base_value],
160 &choice_base_value, PETSC_NULL);
161
162 CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "", test_nb,
163 &test_nb, PETSC_NULL);
164
165 CHKERR PetscOptionsBool("-my_is_partitioned",
166 "set if mesh is partitioned (this result that each "
167 "process keeps only one part of the mesh)",
168 "", is_partitioned, &is_partitioned, PETSC_NULL);
169
170 CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
171 "", "block_conf.in", block_config_file, 255,
172 &flg_block_config);
173
174 CHKERR PetscOptionsBool(
175 "-my_is_calculating_frequency", "set if frequency will be calculated",
176 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
177
178 CHKERR PetscOptionsBool("-my_is_post_proc_volume",
179 "if true post proc volume", "", is_post_proc_volume,
180 &is_post_proc_volume, PETSC_NULL);
181
182 ierr = PetscOptionsEnd();
183 CHKERRG(ierr);
184
185 // Throw error if file with mesh is not provided
186 if (flg_file != PETSC_TRUE) {
187 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
188 }
189
190 // Create mesh database
191 moab::Core mb_instance;
192 moab::Interface &moab = mb_instance;
193
194 // Create moab communicator
195 // Create separate MOAB communicator, it is duplicate of PETSc communicator.
196 // NOTE That this should eliminate potential communication problems between
197 // MOAB and PETSC functions.
198 MPI_Comm moab_comm_world;
199 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
200 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
201 if (pcomm == NULL)
202 pcomm = new ParallelComm(&moab, moab_comm_world);
203
204 // Read whole mesh or part of it if partitioned
205 if (is_partitioned == PETSC_TRUE) {
206 // This is a case of distributed mesh and algebra. In this case each
207 // processor keeps only one part of the problem.
208 const char *option;
210 "PARALLEL_RESOLVE_SHARED_ENTS;"
211 "PARTITION=PARALLEL_PARTITION;";
213 } else {
214 // If that case we have distributed algebra, i.e. assembly of vectors and
215 // matrices is in parallel, but whole mesh is stored on all processors.
216 // Solver and matrix scales well, however problem set-up of problem is
217 // not fully parallel.
218 const char *option;
219 option = "";
221 }
222
223 // Create MoFEM database and link it to MoAB
224 MoFEM::Core core(moab);
225 MoFEM::Interface &m_field = core;
226
227 // Print boundary conditions and material parameters
228 MeshsetsManager *meshsets_mng_ptr;
229 CHKERR m_field.getInterface(meshsets_mng_ptr);
230 CHKERR meshsets_mng_ptr->printDisplacementSet();
231 CHKERR meshsets_mng_ptr->printForceSet();
232 CHKERR meshsets_mng_ptr->printMaterialsSet();
233
234 bool mesh_has_tets = false;
235 bool mesh_has_prisms = false;
236 int nb_tets = 0;
237 int nb_hexs = 0;
238 int nb_prisms = 0;
239
240 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets, true);
241 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs, true);
242 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms, true);
243
244 mesh_has_tets = (nb_tets + nb_hexs) > 0;
245 mesh_has_prisms = nb_prisms > 0;
246
247 // Set bit refinement level to all entities (we do not refine mesh in
248 // this example so all entities have the same bit refinement level)
249 BitRefLevel bit_level0;
250 bit_level0.set(0);
251 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
252 0, 3, bit_level0);
253
254 // CHECK IF EDGE BLOCKSET EXIST AND IF IT IS ADD ALL ENTITIES FROM IT
255 // CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
256 // MESHSET_OF_EDGE_BLOCKSET, 1, bit_level0);
257
259 if (bit->getName().compare(0, 3, "ROD") == 0) {
260 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
261 0, 1, bit_level0);
262 }
263 }
264
265 // Declare approximation fields
267 switch (choice_base_value) {
268 case LEGENDRE:
270 break;
271 case LOBATTO:
273 break;
274 case BERNSTEIN_BEZIER:
276 break;
277 case JACOBI:
279 break;
280 default:
281 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
282 };
283 CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_DENSE,
284 MF_ZERO);
285
286 // We can use higher oder geometry to define body
287 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, base, 3, MB_TAG_DENSE,
288 MF_ZERO);
289
290 // Declare problem
291
292 // Add entities (by tets) to the field (all entities in the mesh, root_set
293 // = 0 )
296
297 // Get all edges in the mesh
298 Range all_edges;
299 CHKERR m_field.get_moab().get_entities_by_type(0, MBEDGE, all_edges, true);
300
301 // Get edges associated with simple rod
302 Range edges_in_simple_rod;
304 if (bit->getName().compare(0, 3, "ROD") == 0) {
305 Range edges;
306 CHKERR m_field.get_moab().get_entities_by_type(bit->getMeshset(),
307 MBEDGE, edges, true);
308 edges_in_simple_rod.merge(edges);
309 }
310 }
311
313 "DISPLACEMENT");
314
315 // Set order of edge in rod to be 1
316 CHKERR m_field.set_field_order(edges_in_simple_rod, "DISPLACEMENT", 1);
317
318 // Get remaining edges (not associated with simple rod) to set order
319 Range edges_to_set_order;
320 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
321
322 // Set approximation order.
323 // See Hierarchical Finite Element Bases on Unstructured Tetrahedral
324 // Meshes.
325 CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
326 CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
327 CHKERR m_field.set_field_order(0, MBHEX, "DISPLACEMENT", order);
328 CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
329 CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
330 CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
331 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
332
334 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
335 else
336 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
337
338 // Set order of approximation of geometry.
339 // Apply 2nd order only on skin (or in whole body)
340 auto setting_second_order_geometry = [&m_field]() {
342 // Setting geometry order everywhere
343 Range tets, edges;
344 CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
345 CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
346 moab::Interface::UNION);
347
348 // Setting 2nd geometry order only on skin
349 // Range tets, faces, edges;
350 // Skinner skin(&m_field.get_moab());
351 // CHKERR skin.find_skin(0,tets,false,faces);
353 // faces,1,false,edges,moab::Interface::UNION
354 // );
355 // CHKERR
356 // m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
357
358 CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
359 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
360
362 };
363 CHKERR setting_second_order_geometry();
364
365 // Configure blocks by parsing config file. It allows setting
366 // approximation order for each block independently.
367 std::map<int, BlockOptionData> block_data;
368 auto setting_blocks_data_and_order_from_config_file =
369 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
370 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
372 if (flg_block_config) {
373 ifstream ini_file(block_config_file);
374 po::variables_map vm;
375 po::options_description config_file_options;
377 it)) {
378 std::ostringstream str_order;
379 str_order << "block_" << it->getMeshsetId()
380 << ".displacement_order";
382 str_order.str().c_str(),
383 po::value<int>(&block_data[it->getMeshsetId()].oRder)
384 ->default_value(order));
385
386 std::ostringstream str_cond;
387 str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
389 str_cond.str().c_str(),
390 po::value<double>(&block_data[it->getMeshsetId()].yOung)
391 ->default_value(-1));
392
393 std::ostringstream str_capa;
394 str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
396 str_capa.str().c_str(),
397 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
398 ->default_value(-2));
399
400 std::ostringstream str_init_temp;
401 str_init_temp << "block_" << it->getMeshsetId()
402 << ".initial_temperature";
404 str_init_temp.str().c_str(),
405 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
406 ->default_value(0));
407 }
408 po::parsed_options parsed =
409 parse_config_file(ini_file, config_file_options, true);
410 store(parsed, vm);
411 po::notify(vm);
413 it)) {
414 if (block_data[it->getMeshsetId()].oRder == -1)
415 continue;
416 if (block_data[it->getMeshsetId()].oRder == order)
417 continue;
418 MOFEM_LOG_C("ELASTIC", Sev::inform, "Set block %d order to %d",
419 it->getMeshsetId(),
420 block_data[it->getMeshsetId()].oRder);
421 Range block_ents;
422 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
423 true);
424 Range ents_to_set_order;
426 ents_to_set_order,
427 moab::Interface::UNION);
428 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
430 ents_to_set_order,
431 moab::Interface::UNION);
433 ents_to_set_order,
434 moab::Interface::UNION);
435 CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
436 ents_to_set_order);
437
438 CHKERR m_field.set_field_order(
439 ents_to_set_order, "DISPLACEMENT",
440 block_data[it->getMeshsetId()].oRder);
441 }
444 collect_unrecognized(parsed.options, po::include_positional);
445 for (std::vector<std::string>::iterator vit =
447 vit != additional_parameters.end(); vit++) {
448 MOFEM_LOG_C("ELASTIC", Sev::warning, "Unrecognized option %s",
449 vit->c_str());
450 }
451 }
452
453 // Update material parameters. Set material parameters block by
454 // block.
456 m_field, BLOCKSET | MAT_ELASTICSET, it)) {
457 const int id = it->getMeshsetId();
458 auto &bd = (*block_sets_ptr)[id];
459 if (block_data[id].yOung > 0)
460 bd.E = block_data[id].yOung;
461 if (block_data[id].pOisson >= -1)
462 bd.PoissonRatio = block_data[id].pOisson;
463 MOFEM_LOG_C("ELASTIC", Sev::inform, "Block %d", id);
464 MOFEM_LOG_C("ELASTIC", Sev::inform, "\tYoung modulus %3.4g", bd.E);
465 MOFEM_LOG_C("ELASTIC", Sev::inform, "\tPoisson ratio %3.4g",
466 bd.PoissonRatio);
467 }
468
470 };
471
473
474 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
475 boost::make_shared<std::map<int, HookeElement::BlockData>>();
476 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
477 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
478
479 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
480 boost::make_shared<std::map<int, MassBlockData>>();
481 CHKERR ConvectiveMassElement::setBlocks(m_field, mass_block_sets_ptr);
482
483 auto fe_lhs_ptr =
484 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
485 auto fe_rhs_ptr =
486 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
487 fe_lhs_ptr->getRuleHook = VolRule();
488 fe_rhs_ptr->getRuleHook = VolRule();
489
490 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_lhs_ptr, true, false, false,
491 false);
492 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_rhs_ptr, true, false, false,
493 false);
494
495 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
496 new PrismFE(m_field));
497 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
498 new PrismFE(m_field));
499
501 "DISPLACEMENT",
502 "MESH_NODE_POSITIONS", false);
503
504 auto add_skin_element_for_post_processing = [&]() {
506 Range elastic_element_ents;
508 "ELASTIC", 3, elastic_element_ents);
509 Skinner skin(&m_field.get_moab());
510 Range skin_faces; // skin faces from 3d ents
511 CHKERR skin.find_skin(0, elastic_element_ents, false, skin_faces);
512 Range proc_skin;
513 if (is_partitioned) {
514 CHKERR pcomm->filter_pstatus(skin_faces,
515 PSTATUS_SHARED | PSTATUS_MULTISHARED,
516 PSTATUS_NOT, -1, &proc_skin);
517 } else {
518 proc_skin = skin_faces;
519 }
522 "DISPLACEMENT");
524 "DISPLACEMENT");
526 "DISPLACEMENT");
527 if (m_field.check_field("TEMP")) {
529 // "TEMP");
531 // "TEMP");
533 "TEMP");
534 }
536 "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
538 "POST_PROC_SKIN");
540 };
542
543 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
544 if (mesh_has_tets) {
545 CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
546 "DISPLACEMENT", "MESH_NODE_POSITIONS",
547 false, true, MBTET, data_at_pts);
548 }
549 if (mesh_has_prisms) {
550 CHKERR HookeElement::setOperators(
551 prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
552 "MESH_NODE_POSITIONS", false, true, MBPRISM, data_at_pts);
553 }
554
555 if (test_nb == 4) {
556
557 auto thermal_strain =
560 constexpr double alpha = 1;
564 t_thermal_strain(i, j) = alpha * t_coords(2) * t_kd(i, j);
565 return t_thermal_strain;
566 };
567
568 fe_rhs_ptr->getOpPtrVector().push_back(
569 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
570 "DISPLACEMENT", data_at_pts, thermal_strain));
571 }
572
573 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
575
576 for (auto &sit : *block_sets_ptr) {
577 for (auto &mit : *mass_block_sets_ptr) {
578 fe_mass_ptr->getOpPtrVector().push_back(
579 new HookeElement::OpCalculateMassMatrix("DISPLACEMENT",
580 "DISPLACEMENT", sit.second,
581 mit.second, data_at_pts));
582 }
583 }
584
585 // Add spring boundary condition applied on surfaces.
586 // This is only declaration not implementation.
588 "MESH_NODE_POSITIONS");
589
590 // Implementation of spring element
591 // Create new instances of face elements for springs
592 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
594 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
596
597 CHKERR MetaSpringBC::setSpringOperators(m_field, fe_spring_lhs_ptr,
598 fe_spring_rhs_ptr, "DISPLACEMENT",
599 "MESH_NODE_POSITIONS");
600
601 // Add Simple Rod elements
602 // This is only declaration not implementation.
604 "MESH_NODE_POSITIONS");
605
607 // MBEDGE, "SIMPLE_ROD");
608
609 // Implementation of Simple Rod element
610 // Create new instances of edge elements for Simple Rod
611 boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(new EdgeEle(m_field));
612 boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(new EdgeEle(m_field));
613
614
616 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
617 "MESH_NODE_POSITIONS");
618
619 // Add body force element. This is only declaration of element. not its
620 // implementation.
623 "DISPLACEMENT");
625 "DISPLACEMENT");
627 "DISPLACEMENT");
629 "MESH_NODE_POSITIONS");
630
632 m_field, BLOCKSET | BODYFORCESSET, it)) {
633 Range tets;
634 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
635 true);
637 }
638 CHKERR m_field.build_finite_elements("BODY_FORCE");
639
640 // Add Neumann forces, i.e. pressure or traction forces applied on body
641 // surface. This is only declaration not implementation.
645
646 // Add fluid pressure finite elements. This is special pressure on the
647 // surface from fluid, i.e. pressure which linearly change with the depth.
648 FluidPressure fluid_pressure_fe(m_field);
649 // This function only declare element. Element is implemented by operators
650 // in class FluidPressure.
652
653 // Add elements for thermo elasticity if temperature field is defined.
654 ThermalStressElement thermal_stress_elem(m_field);
655 // Check if TEMP field exist, and then add element.
656 if (!m_field.check_field("TEMP")) {
659 if (block_data[it->getMeshsetId()].initTemp != 0) {
661 break;
662 }
663 }
665 CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
666 MB_TAG_SPARSE, MF_ZERO);
667
669 CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
670 }
671 }
672 if (m_field.check_field("TEMP")) {
674 "ELASTIC", "DISPLACEMENT", "TEMP");
675 }
676
677 // All is declared, at this point build fields first,
678 CHKERR m_field.build_fields();
679 // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
680 // Projection10NodeCoordsOnField
681 // do the trick.
682 Projection10NodeCoordsOnField ent_method_material(m_field,
683 "MESH_NODE_POSITIONS");
684 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
685 if (m_field.check_field("TEMP")) {
687 if (block_data[it->getMeshsetId()].initTemp != 0) {
688 MOFEM_LOG_C("ELASTIC", Sev::inform,
689 "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
690 block_data[it->getMeshsetId()].initTemp);
691 Range block_ents;
692 CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
693 Range vertices;
694 CHKERR moab.get_connectivity(block_ents, vertices, true);
695 CHKERR m_field.getInterface<FieldBlas>()->setField(
696 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
697 "TEMP");
698 }
699 }
700 }
701
702 // Build database for elements. Actual implementation of element is not need
703 // here, only elements has to be declared.
705 // Build adjacencies between elements and field entities
707
708 // Register MOFEM DM implementation in PETSc
710
711 // Create DM manager
712 auto dm = createDM(PETSC_COMM_WORLD, "MOFEM");
713 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
714 CHKERR DMSetFromOptions(dm);
715 CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
716 // Add elements to DM manager
725 CHKERR DMSetUp(dm);
726
727 // Create matrices & vectors. Note that native PETSc DM interface is used,
728 // but under the PETSc interface MoFEM implementation is running.
731 auto D = vectorDuplicate(F);
732 auto D0 = vectorDuplicate(F);
735 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
736
737 // Initialise mass matrix
739 if (is_calculating_frequency == PETSC_TRUE) {
740 Mij = matDuplicate(Aij, MAT_DO_NOT_COPY_VALUES);
741 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
742 // MatView(Mij, PETSC_VIEWER_STDOUT_SELF);
743 }
744
745 // Assign global matrix/vector contributed by springs
746 fe_spring_lhs_ptr->ksp_B = Aij;
747 fe_spring_rhs_ptr->ksp_f = F;
748
749 // Assign global matrix/vector contributed by Simple Rod
750 fe_simple_rod_lhs_ptr->ksp_B = Aij;
751 fe_simple_rod_rhs_ptr->ksp_f = F;
752
753 // Zero vectors and matrices
754 CHKERR VecZeroEntries(F);
755 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
756 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
757 CHKERR VecZeroEntries(D);
758 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
759 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
760 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
761 CHKERR MatZeroEntries(Aij);
762
763 // Below particular implementations of finite elements are used to assemble
764 // problem matrixes and vectors. Implementation of element does not change
765 // how element is declared.
766
767 // Assemble Aij and F. Define Dirichlet bc element, which sets constrains
768 // to MatrixDouble and the right hand side vector.
769
771 // DirichletDisplacementBc will try to use DISPLACEMENT blockset. Two
772 // implementations are available, depending how BC is defined on mesh file.
773 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
774 m_field, "DISPLACEMENT", Aij, D0, F);
775
776 // That sets Dirichlet bc objects that problem is linear, i.e. no newton
777 // (SNES) solver is run for this problem.
778 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
779 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
780
781 // D0 vector will store initial displacements
782 CHKERR VecZeroEntries(D0);
783 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
784 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
785 CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
786 // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
787 // implementation
788 // of particular dirichlet_bc.
789 CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
790 // Set values from D0 on the field (on the mesh)
791
792 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
793 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
794 CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
795
796 // Calculate residual forces as result of applied kinematic constrains. Run
797 // implementation
798 // of particular finite element implementation. Look how
799 // NonlinearElasticElement is implemented,
800 // in that case. We run NonlinearElasticElement with hook material.
801 // Calculate right hand side vector
802 fe_rhs_ptr->snes_f = F;
803 prism_fe_rhs_ptr->snes_f = F;
804 MOFEM_LOG("ELASTIC", Sev::inform) << "Assemble external force vector ...";
805 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
806 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
807 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
808 // Assemble matrix
809 fe_lhs_ptr->snes_B = Aij;
810 prism_fe_lhs_ptr->snes_B = Aij;
811 MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate stiffness matrix ...";
812 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
813 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
814 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
815
816 // Assemble springs
817 CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_lhs_ptr);
818 CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_rhs_ptr);
819
820 // Assemble Simple Rod
821 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_lhs_ptr);
822 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_rhs_ptr);
823
824 if (is_calculating_frequency == PETSC_TRUE) {
825 // Assemble mass matrix
826 fe_mass_ptr->snes_B = Mij;
827 MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate mass matrix ...";
828 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_mass_ptr);
829 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
830 }
831
832 // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
833
834 // Assemble pressure and traction forces. Run particular implemented for do
835 // this, see
836 // MetaNeumannForces how this is implemented.
837 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
839 F, "DISPLACEMENT");
840
841 {
842 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
843 neumann_forces.begin();
844 for (; mit != neumann_forces.end(); mit++) {
845 CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
846 &mit->second->getLoopFe());
847 }
848 }
849 // Assemble forces applied to nodes, see implementation in NodalForce
850 boost::ptr_map<std::string, NodalForce> nodal_forces;
851 CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
852 "DISPLACEMENT");
853
854 {
855 boost::ptr_map<std::string, NodalForce>::iterator fit =
856 nodal_forces.begin();
857 for (; fit != nodal_forces.end(); fit++) {
858 CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
859 &fit->second->getLoopFe());
860 }
861 }
862 // Assemble edge forces
863 boost::ptr_map<std::string, EdgeForce> edge_forces;
864 CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
865 "DISPLACEMENT");
866 {
867 auto fit = edge_forces.begin();
868 for (; fit != edge_forces.end(); fit++) {
869 auto &fe = fit->second->getLoopFe();
870 CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(), &fe);
871 }
872 }
873 // Assemble body forces, implemented in BodyForceConstantField
874 BodyForceConstantField body_forces_methods(m_field);
876 m_field, BLOCKSET | BODYFORCESSET, it)) {
878 it->getMeshsetId());
879 }
880 CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
881 &body_forces_methods.getLoopFe());
882 // Assemble fluid pressure forces
884 false, false);
886 "DISPLACEMENT", F, false, true);
887
888 CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
889 &fluid_pressure_fe.getLoopFe());
890 // Apply kinematic constrain to right hand side vector and matrix
891 CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
892
893 // Matrix View
894 PetscViewerPushFormat(
895 PETSC_VIEWER_STDOUT_SELF,
896 PETSC_VIEWER_ASCII_MATLAB); /// PETSC_VIEWER_ASCII_DENSE,
897 /// PETSC_VIEWER_ASCII_MATLAB
898 // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
899 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
900 // std::string wait;
901 // std::cin >> wait;
902
903 if (is_calculating_frequency == PETSC_TRUE) {
904 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
905 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
906 }
907
908 // Set matrix positive defined and symmetric for Cholesky and icc
909 // pre-conditioner
910
911 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
914 CHKERR VecAssemblyBegin(F);
915 CHKERR VecAssemblyEnd(F);
916 CHKERR VecScale(F, -1);
917
918 // Create solver
919 auto solver = createKSP(PETSC_COMM_WORLD);
920 CHKERR KSPSetDM(solver, dm);
921 CHKERR KSPSetFromOptions(solver);
922 CHKERR KSPSetOperators(solver, Aij, Aij);
923
924 // Setup multi-grid pre-conditioner if set from command line
925 {
926 // from PETSc example ex42.c
927 PetscBool same = PETSC_FALSE;
928 PC pc;
929 CHKERR KSPGetPC(solver, &pc);
930 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
931 if (same) {
932 PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
933 CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
934 CHKERR PCSetFromOptions(pc);
935 } else {
936 // Operators are already set, do not use DM for doing that
937 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
938 }
939 }
940 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
941 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
942 // Set up solver
943 CHKERR KSPSetUp(solver);
944
945 auto set_post_proc_skin = [&](auto &post_proc_skin) {
948 false);
949 CHKERR post_proc_skin.generateReferenceElementMesh();
953 "DISPLACEMENT", "ELASTIC", data_at_pts->hMat, true);
955 "MESH_NODE_POSITIONS", "ELASTIC", data_at_pts->HMat, false);
956 if (m_field.check_field("TEMP")) {
959 }
960 post_proc_skin.getOpPtrVector().push_back(
961 new HookeElement::OpPostProcHookeElement<
963 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
964 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
965 true));
967 };
968
969 auto set_post_proc_tets = [&](auto &post_proc) {
971 // Add operators to the elements, starting with some generic operators
972 CHKERR post_proc.generateReferenceElementMesh();
973 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, true, false, false,
974 false);
978 if (m_field.check_field("TEMP")) {
981 }
982 // Add problem specific operator on element to post-process stresses
983 post_proc.getOpPtrVector().push_back(new PostProcHookStress(
984 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
985 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
987 };
988
989 auto set_post_proc_edge = [&](auto &post_proc_edge) {
991 CHKERR post_proc_edge.generateReferenceElementMesh();
994 };
995
996 auto set_post_proc_prisms = [&](auto &prism_post_proc) {
998 CHKERR prism_post_proc.generateReferenceElementMesh();
999 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
1000 prism_post_proc.getOpPtrVector().push_back(
1001 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
1002 prism_post_proc.getOpPtrVector().push_back(
1003 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
1007 prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
1008 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
1009 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
1011 };
1012
1013 PostProcFaceOnRefinedMesh post_proc_skin(m_field);
1014 PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
1015 PostProcEdgeOnRefinedMesh post_proc_edge(m_field);
1016 PostProcVolumeOnRefinedMesh post_proc(m_field);
1017
1018 CHKERR set_post_proc_skin(post_proc_skin);
1019 CHKERR set_post_proc_tets(post_proc);
1020 CHKERR set_post_proc_prisms(prism_post_proc);
1021 CHKERR set_post_proc_edge(post_proc_edge);
1022
1023 PetscBool field_eval_flag = PETSC_FALSE;
1024 std::array<double, 3> field_eval_coords;
1025 boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> field_eval_data;
1026 PetscInt coords_dim = 3;
1027 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1028 field_eval_coords.data(), &coords_dim,
1029 &field_eval_flag);
1030
1031 auto scalar_field_ptr = boost::make_shared<VectorDouble>();
1032 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
1033 auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
1034
1035 if (field_eval_flag) {
1036 field_eval_data = m_field.getInterface<FieldEvaluatorInterface>()
1037 ->getData<VolumeElementForcesAndSourcesCore>();
1038 CHKERR m_field.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1039 field_eval_data, "ELASTIC");
1040
1041 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1042 auto no_rule = [](int, int, int) { return -1; };
1043
1044 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1045 field_eval_fe_ptr->getRuleHook = no_rule;
1046
1047 if (m_field.check_field("TEMP")) {
1048 field_eval_fe_ptr->getOpPtrVector().push_back(
1049 new OpCalculateScalarFieldValues("TEMP", scalar_field_ptr));
1050 }
1051 field_eval_fe_ptr->getOpPtrVector().push_back(
1052 new OpCalculateVectorFieldValues<3>("DISPLACEMENT", vector_field_ptr));
1053 field_eval_fe_ptr->getOpPtrVector().push_back(
1055 }
1056
1057 // Temperature field is defined on the mesh
1058 if (m_field.check_field("TEMP")) {
1059
1060 // Create thermal vector
1061 Vec F_thermal;
1062 CHKERR VecDuplicate(F, &F_thermal);
1063
1064 // Set up implementation for calculation of thermal stress vector. Look
1065 // how thermal stresses and vector is assembled in ThermalStressElement.
1066 CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1067 "DISPLACEMENT", "TEMP", F_thermal);
1068
1069 SeriesRecorder *recorder_ptr;
1070 CHKERR m_field.getInterface(recorder_ptr);
1071
1072 // Read time series and do thermo-elastic analysis, this is when time
1073 // dependent
1074 // temperature problem was run before on the mesh. It means that before
1075 // non-stationary
1076 // problem was solved for temperature and filed "TEMP" is stored for
1077 // subsequent time
1078 // steps in the recorder.
1079 if (recorder_ptr->check_series("THEMP_SERIES")) {
1080 // This is time dependent case, so loop of data series stored by tape
1081 // recorder.
1082 // Loop over time steps
1083 for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
1084 sit)) {
1085 MOFEM_LOG_C("ELASTIC", Sev::inform, "Process step %d",
1086 sit->get_step_number());
1087 // Load field data for this time step
1089 sit->get_step_number());
1090
1091 CHKERR VecZeroEntries(F_thermal);
1092 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1093 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1094
1095 // Calculate the right-hand side vector as result of thermal stresses.
1096 // It MetaNodalForces
1097 // that on "ELASTIC" element data structure the element implementation
1098 // in thermal_stress_elem
1099 // is executed.
1101 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1102
1103 // Assemble vector
1104 CHKERR VecAssemblyBegin(F_thermal);
1105 CHKERR VecAssemblyEnd(F_thermal);
1106 // Accumulate ghost dofs
1109
1110 // Calculate norm of vector and print values
1111 PetscReal nrm_F;
1112 CHKERR VecNorm(F, NORM_2, &nrm_F);
1113
1114 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1115 PetscReal nrm_F_thermal;
1116 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1117 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1118 nrm_F_thermal);
1119
1120 CHKERR VecScale(F_thermal, -1);
1121 // check this !!!
1122 CHKERR VecAXPY(F_thermal, 1, F);
1123
1124 // Set dirichlet boundary to thermal stresses vector
1125 dirichlet_bc_ptr->snes_x = D;
1126 dirichlet_bc_ptr->snes_f = F_thermal;
1127 CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1128
1129 // Solve problem
1130 CHKERR KSPSolve(solver, F_thermal, D);
1131 // Add boundary conditions vector
1132 CHKERR VecAXPY(D, 1., D0);
1133 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1134 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1135
1136 // Save data on the mesh
1137 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1138
1139 // Save data on mesh
1140 CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1141
1142 if (field_eval_flag) {
1144 ->evalFEAtThePoint3D(
1145 field_eval_coords.data(), 1e-12, "ELASTIC_PROB", "ELASTIC",
1146 field_eval_data, m_field.get_comm_rank(),
1147 m_field.get_comm_rank(), nullptr, MF_EXIST, QUIET);
1148 if (scalar_field_ptr->size()) {
1149 auto t_temp = getFTensor0FromVec(*scalar_field_ptr);
1150 MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
1151 << "Eval point TEMP: " << t_temp;
1152 }
1153 if (vector_field_ptr->size1()) {
1155 auto t_disp = getFTensor1FromMat<3>(*vector_field_ptr);
1156 MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
1157 << "Eval point DISPLACEMENT magnitude: "
1158 << sqrt(t_disp(i) * t_disp(i));
1159 }
1160 if (tensor_field_ptr->size1()) {
1162 auto t_disp_grad = getFTensor2FromMat<3, 3>(*tensor_field_ptr);
1163 MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
1165 }
1166
1168 }
1169
1170 // Post-process results
1171 if (is_post_proc_volume == PETSC_TRUE) {
1172 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1173 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1174 std::ostringstream o1;
1175 o1 << "out_" << sit->step_number << ".h5m";
1176 if (!test_nb)
1177 CHKERR post_proc.writeFile(o1.str().c_str());
1178 MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
1179 }
1180
1181 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1182 CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN",
1183 &post_proc_skin);
1184 std::ostringstream o1_skin;
1185 o1_skin << "out_skin_" << sit->step_number << ".h5m";
1186 if (!test_nb)
1187 CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1188 MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
1189 }
1190 } else {
1191
1192 // This is a case when stationary problem for temperature was solved.
1193 CHKERR VecZeroEntries(F_thermal);
1194 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1195 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1196
1197 // Calculate the right-hand side vector with thermal stresses
1199 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1200
1201 // Assemble vector
1202 CHKERR VecAssemblyBegin(F_thermal);
1203 CHKERR VecAssemblyEnd(F_thermal);
1204
1205 // Accumulate ghost dofs
1208
1209 // Calculate norm
1210 PetscReal nrm_F;
1211 CHKERR VecNorm(F, NORM_2, &nrm_F);
1212
1213 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1214 PetscReal nrm_F_thermal;
1215 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1216
1217 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1218 nrm_F_thermal);
1219
1220 // Add thermal stress vector and other forces vector
1221 CHKERR VecScale(F_thermal, -1);
1222 CHKERR VecAXPY(F_thermal, 1, F);
1223
1224 // Apply kinetic boundary conditions
1225 dirichlet_bc_ptr->snes_x = D;
1226 dirichlet_bc_ptr->snes_f = F_thermal;
1227 CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1228
1229 // Solve problem
1230 CHKERR KSPSolve(solver, F_thermal, D);
1231 CHKERR VecAXPY(D, 1., D0);
1232
1233 // Update ghost values for solution vector
1234 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1235 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1236 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1237
1238 // Save data on mesh
1239 if (is_post_proc_volume == PETSC_TRUE) {
1240 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1241 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1242 // Save results to file
1243 if (!test_nb)
1244 CHKERR post_proc.writeFile("out.h5m");
1245 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1246 }
1247
1248 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1249 CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1250 if (!test_nb)
1251 CHKERR post_proc_skin.writeFile("out_skin.h5m");
1252 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1253 }
1254
1255 // Destroy vector, no needed any more
1256 CHKERR VecDestroy(&F_thermal);
1257 } else {
1258 // Elastic analysis no temperature field
1259 // VecView(F, PETSC_VIEWER_STDOUT_WORLD);
1260 // Solve for vector D
1261 CHKERR KSPSolve(solver, F, D);
1262
1263 // VecView(D, PETSC_VIEWER_STDOUT_WORLD);
1264 // cerr << F;
1265
1266 // Add kinetic boundary conditions
1267 CHKERR VecAXPY(D, 1., D0);
1268 // Update ghost values
1269 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1270 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1271 // Save data from vector on mesh
1272 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1273 // Post-process results
1274 MOFEM_LOG("ELASTIC", Sev::inform) << "Post-process start ...";
1275 if (is_post_proc_volume == PETSC_TRUE) {
1276 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1277 }
1278 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
1279 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
1280 CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1281 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1282
1283 // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
1284 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1285 if (mesh_has_tets) {
1286 if (is_post_proc_volume == PETSC_TRUE) {
1287 if (!test_nb)
1288 CHKERR post_proc.writeFile("out.h5m");
1289 }
1290 if (!test_nb)
1291 CHKERR post_proc_skin.writeFile("out_skin.h5m");
1292 }
1293 if (mesh_has_prisms) {
1294 if (!test_nb)
1295 CHKERR prism_post_proc.writeFile("prism_out.h5m");
1296 }
1297 if (!edges_in_simple_rod.empty())
1298 if (!test_nb)
1299 CHKERR post_proc_edge.writeFile("out_edge.h5m");
1300 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1301 }
1302
1303 if (is_calculating_frequency == PETSC_TRUE) {
1304 // Calculate mode mass, m = u^T * M * u
1305 Vec u1;
1306 VecDuplicate(D, &u1);
1307 CHKERR MatMult(Mij, D, u1);
1308 double mode_mass;
1309 CHKERR VecDot(u1, D, &mode_mass);
1310 MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode mass %6.4e\n", mode_mass);
1311
1312 Vec v1;
1313 VecDuplicate(D, &v1);
1314 CHKERR MatMult(Aij, D, v1);
1315
1316 double mode_stiffness;
1317 CHKERR VecDot(v1, D, &mode_stiffness);
1318 MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode stiffness %6.4e\n",
1319 mode_stiffness);
1320
1321 double frequency;
1322 double pi = 3.14159265359;
1323 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1324 MOFEM_LOG_C("ELASTIC", Sev::inform, "Frequency %6.4e", frequency);
1325 }
1326
1327 // Calculate elastic energy
1328 auto calculate_strain_energy = [&]() {
1330
1331 SmartPetscObj<Vec> v_energy;
1332 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
1333 "MESH_NODE_POSITIONS", false, true,
1334 v_energy);
1335
1336 // Print elastic energy
1337 double energy;
1338 CHKERR VecSum(v_energy, &energy);
1339 MOFEM_LOG_C("ELASTIC", Sev::inform, "Elastic energy %6.4e", energy);
1340
1341 switch (test_nb) {
1342 case 1:
1343 if (fabs(energy - 17.129) > 1e-3)
1344 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1345 "atom test diverged!");
1346 break;
1347 case 2:
1348 if (fabs(energy - 5.6475e-03) > 1e-4)
1349 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1350 "atom test diverged!");
1351 break;
1352 case 3:
1353 if (fabs(energy - 7.4679e-03) > 1e-4)
1354 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1355 "atom test diverged!");
1356 break;
1357 case 4:
1358 if (fabs(energy - 2.4992e+00) > 1e-3)
1359 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1360 "atom test diverged!");
1361 break;
1362 // FIXME: Here are missing regersion tests
1363 case 8: {
1364 double min;
1365 CHKERR VecMin(D, PETSC_NULL, &min);
1366 constexpr double expected_val = 0.10001;
1367 if (fabs(min + expected_val) > 1e-10)
1368 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1369 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1370 } break;
1371 case 9: {
1372 if (fabs(energy - 4.7416e-04) > 1e-8)
1373 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1374 "atom test diverged!");
1375 }
1376 default:
1377 break;
1378 }
1379
1381 };
1382 CHKERR calculate_strain_energy();
1383
1384 MPI_Comm_free(&moab_comm_world);
1385 }
1387
1389
1390 return 0;
1391}
const std::string default_options
std::string param_file
Implementation of linear elastic material.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
int main()
Kronecker Delta class symmetric.
@ QUIET
Definition: definitions.h:208
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ MF_EXIST
Definition: definitions.h:100
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ NOBASE
Definition: definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: elasticity.cpp:94
static char help[]
Definition: elasticity.cpp:57
@ F
constexpr auto t_kd
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
virtual bool check_field(const std::string &name) const =0
check if field is in database
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"
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
virtual bool check_series(const std::string &name) const
check if series is in database
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
double D
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createKSP(MPI_Comm comm)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Body forces elements.
Definition: BodyForce.hpp:12
MyVolumeFE & getLoopFe()
Definition: BodyForce.hpp:23
MoFEMErrorCode addBlock(const std::string field_name, Vec F, int ms_id)
Definition: BodyForce.hpp:94
data for calculation inertia forces
Fluid pressure forces.
MyTriangleFE & getLoopFe()
MoFEMErrorCode setNeumannFluidPressureFiniteElementOperators(string field_name, Vec F, bool allow_negative_pressure=true, bool ho_geometry=false)
MoFEMErrorCode addNeumannFluidPressureBCElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
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
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
static MoFEMErrorCode setSimpleRodOperators(MoFEM::Interface &m_field, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_lhs_ptr, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of SimpleRod element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Managing BitRefLevels.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
virtual MPI_Comm & get_comm() const =0
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
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
Field evaluator interface.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Transform local reference derivatives of shape functions to global derivatives.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
data for calculation heat conductivity and heat capacity elements
Set data structures of MG pre-conditioner via approximation orders.
Postprocess on edge.
Postprocess on face.
Postprocess on prism.
Operator post-procesing stresses for Hook isotropic material.
Post processing.
int getRuleThroughThickness(int order)
int getRuleTrianglesOnly(int order)
Implentation of thermal stress element.
MyVolumeFE & getLoopThermalStressRhs()
MoFEMErrorCode addThermalStressElement(const std::string fe_name, const std::string field_name, const std::string thermal_field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
MoFEMErrorCode setThermalStressRhsOperators(string field_name, string thermal_field_name, Vec &F, int verb=0)
Set integration rule.
int operator()(int, int, int order) const
Definition: elasticity.cpp:80