v0.13.1
Classes | Typedefs | Functions | Variables
elasticity.cpp File Reference
#include <BasicFiniteElements.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes

struct  BlockOptionData
 
struct  VolRule
 Set integration rule. More...
 
struct  PrismFE
 

Typedefs

using BlockData = NonlinearElasticElement::BlockData
 
using MassBlockData = ConvectiveMassElement::BlockData
 
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
 
using EdgeEle = MoFEM::EdgeElementForcesAndSourcesCore
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help []
 

Typedef Documentation

◆ BlockData

◆ EdgeEle

Examples
elasticity.cpp.

Definition at line 94 of file elasticity.cpp.

◆ MassBlockData

Definition at line 75 of file elasticity.cpp.

◆ VolUserDataOperator

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
elasticity.cpp.

Definition at line 96 of file elasticity.cpp.

96 {
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();
122 core_log->add_sink(
123 LogManager::createSink(LogManager::getStrmWorld(), "ELASTIC"));
124 LogManager::setLog("ELASTIC");
125 MOFEM_LOG_TAG("ELASTIC", "elasticity")
126
127 try {
128
129 PetscBool flg_block_config, flg_file;
130 char mesh_file_name[255];
131 char block_config_file[255];
132 PetscInt test_nb = 0;
133 PetscInt order = 2;
134 PetscBool is_partitioned = PETSC_FALSE;
135 PetscBool is_calculating_frequency = PETSC_FALSE;
136 PetscBool is_post_proc_volume = PETSC_TRUE;
137
138 // Select base
139 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
140 const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
141 "jacobi"};
142 PetscInt choice_base_value = LOBATTO;
143
144 // Read options from command line
145 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
146 CHKERR(ierr);
147 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
148 mesh_file_name, 255, &flg_file);
149
150 CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
151 order, &order, PETSC_NULL);
152
153 CHKERR PetscOptionsEList("-base", "approximation base", "", list_bases,
154 LASBASETOP, list_bases[choice_base_value],
155 &choice_base_value, PETSC_NULL);
156
157 CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "", test_nb,
158 &test_nb, PETSC_NULL);
159
160 CHKERR PetscOptionsBool("-my_is_partitioned",
161 "set if mesh is partitioned (this result that each "
162 "process keeps only one part of the mesh)",
163 "", is_partitioned, &is_partitioned, PETSC_NULL);
164
165 CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
166 "", "block_conf.in", block_config_file, 255,
167 &flg_block_config);
168
169 CHKERR PetscOptionsBool(
170 "-my_is_calculating_frequency", "set if frequency will be calculated",
171 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
172
173 CHKERR PetscOptionsBool("-my_is_post_proc_volume",
174 "if true post proc volume", "", is_post_proc_volume,
175 &is_post_proc_volume, PETSC_NULL);
176
177 ierr = PetscOptionsEnd();
178 CHKERRG(ierr);
179
180 // Throw error if file with mesh is not provided
181 if (flg_file != PETSC_TRUE) {
182 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
183 }
184
185 // Create mesh database
186 moab::Core mb_instance;
187 moab::Interface &moab = mb_instance;
188
189 // Create moab communicator
190 // Create separate MOAB communicator, it is duplicate of PETSc communicator.
191 // NOTE That this should eliminate potential communication problems between
192 // MOAB and PETSC functions.
193 MPI_Comm moab_comm_world;
194 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
195 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
196 if (pcomm == NULL)
197 pcomm = new ParallelComm(&moab, moab_comm_world);
198
199 // Read whole mesh or part of it if partitioned
200 if (is_partitioned == PETSC_TRUE) {
201 // This is a case of distributed mesh and algebra. In this case each
202 // processor keeps only one part of the problem.
203 const char *option;
204 option = "PARALLEL=READ_PART;"
205 "PARALLEL_RESOLVE_SHARED_ENTS;"
206 "PARTITION=PARALLEL_PARTITION;";
207 CHKERR moab.load_file(mesh_file_name, 0, option);
208 } else {
209 // If that case we have distributed algebra, i.e. assembly of vectors and
210 // matrices is in parallel, but whole mesh is stored on all processors.
211 // Solver and matrix scales well, however problem set-up of problem is
212 // not fully parallel.
213 const char *option;
214 option = "";
215 CHKERR moab.load_file(mesh_file_name, 0, option);
216 }
217
218 // Create MoFEM database and link it to MoAB
219 MoFEM::Core core(moab);
220 MoFEM::Interface &m_field = core;
221
222 // Print boundary conditions and material parameters
223 MeshsetsManager *meshsets_mng_ptr;
224 CHKERR m_field.getInterface(meshsets_mng_ptr);
225 CHKERR meshsets_mng_ptr->printDisplacementSet();
226 CHKERR meshsets_mng_ptr->printForceSet();
227 CHKERR meshsets_mng_ptr->printMaterialsSet();
228
229 bool mesh_has_tets = false;
230 bool mesh_has_prisms = false;
231 int nb_tets = 0;
232 int nb_hexs = 0;
233 int nb_prisms = 0;
234
235 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets, true);
236 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs, true);
237 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms, true);
238
239 mesh_has_tets = (nb_tets + nb_hexs) > 0;
240 mesh_has_prisms = nb_prisms > 0;
241
242 // Set bit refinement level to all entities (we do not refine mesh in
243 // this example so all entities have the same bit refinement level)
244 BitRefLevel bit_level0;
245 bit_level0.set(0);
246 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
247 0, 3, bit_level0);
248
249 // CHECK IF EDGE BLOCKSET EXIST AND IF IT IS ADD ALL ENTITIES FROM IT
250 // CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
251 // MESHSET_OF_EDGE_BLOCKSET, 1, bit_level0);
252
254 if (bit->getName().compare(0, 3, "ROD") == 0) {
255 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
256 0, 1, bit_level0);
257 }
258 }
259
260 // Declare approximation fields
262 switch (choice_base_value) {
263 case LEGENDRE:
265 break;
266 case LOBATTO:
268 break;
269 case BERNSTEIN_BEZIER:
271 break;
272 case JACOBI:
274 break;
275 default:
276 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
277 };
278 CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_DENSE,
279 MF_ZERO);
280
281 // We can use higher oder geometry to define body
282 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, base, 3, MB_TAG_DENSE,
283 MF_ZERO);
284
285 // Declare problem
286
287 // Add entities (by tets) to the field (all entities in the mesh, root_set
288 // = 0 )
289 CHKERR m_field.add_ents_to_field_by_dim(0, 3, "DISPLACEMENT");
290 CHKERR m_field.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
291
292 // Get all edges in the mesh
293 Range all_edges;
294 CHKERR m_field.get_moab().get_entities_by_type(0, MBEDGE, all_edges, true);
295
296 // Get edges associated with simple rod
297 Range edges_in_simple_rod;
299 if (bit->getName().compare(0, 3, "ROD") == 0) {
300 Range edges;
301 CHKERR m_field.get_moab().get_entities_by_type(bit->getMeshset(),
302 MBEDGE, edges, true);
303 edges_in_simple_rod.merge(edges);
304 }
305 }
306
307 CHKERR m_field.add_ents_to_field_by_type(edges_in_simple_rod, MBEDGE,
308 "DISPLACEMENT");
309
310 // Set order of edge in rod to be 1
311 CHKERR m_field.set_field_order(edges_in_simple_rod, "DISPLACEMENT", 1);
312
313 // Get remaining edges (not associated with simple rod) to set order
314 Range edges_to_set_order;
315 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
316
317 // Set approximation order.
318 // See Hierarchical Finite Element Bases on Unstructured Tetrahedral
319 // Meshes.
320 CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
321 CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
322 CHKERR m_field.set_field_order(0, MBHEX, "DISPLACEMENT", order);
323 CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
324 CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
325 CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
326 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
327
329 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
330 else
331 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
332
333 // Set order of approximation of geometry.
334 // Apply 2nd order only on skin (or in whole body)
335 auto setting_second_order_geometry = [&m_field]() {
337 // Setting geometry order everywhere
338 Range tets, edges;
339 CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
340 CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
341 moab::Interface::UNION);
342
343 // Setting 2nd geometry order only on skin
344 // Range tets, faces, edges;
345 // Skinner skin(&m_field.get_moab());
346 // CHKERR skin.find_skin(0,tets,false,faces);
347 // CHKERR m_field.get_moab().get_adjacencies(
348 // faces,1,false,edges,moab::Interface::UNION
349 // );
350 // CHKERR
351 // m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
352
353 CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
354 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
355
357 };
358 CHKERR setting_second_order_geometry();
359
360 // Configure blocks by parsing config file. It allows setting
361 // approximation order for each block independently.
362 std::map<int, BlockOptionData> block_data;
363 auto setting_blocks_data_and_order_from_config_file =
364 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
365 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
367 if (flg_block_config) {
368 ifstream ini_file(block_config_file);
369 po::variables_map vm;
370 po::options_description config_file_options;
372 it)) {
373 std::ostringstream str_order;
374 str_order << "block_" << it->getMeshsetId()
375 << ".displacement_order";
376 config_file_options.add_options()(
377 str_order.str().c_str(),
378 po::value<int>(&block_data[it->getMeshsetId()].oRder)
379 ->default_value(order));
380
381 std::ostringstream str_cond;
382 str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
383 config_file_options.add_options()(
384 str_cond.str().c_str(),
385 po::value<double>(&block_data[it->getMeshsetId()].yOung)
386 ->default_value(-1));
387
388 std::ostringstream str_capa;
389 str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
390 config_file_options.add_options()(
391 str_capa.str().c_str(),
392 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
393 ->default_value(-2));
394
395 std::ostringstream str_init_temp;
396 str_init_temp << "block_" << it->getMeshsetId()
397 << ".initial_temperature";
398 config_file_options.add_options()(
399 str_init_temp.str().c_str(),
400 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
401 ->default_value(0));
402 }
403 po::parsed_options parsed =
404 parse_config_file(ini_file, config_file_options, true);
405 store(parsed, vm);
406 po::notify(vm);
408 it)) {
409 if (block_data[it->getMeshsetId()].oRder == -1)
410 continue;
411 if (block_data[it->getMeshsetId()].oRder == order)
412 continue;
413 MOFEM_LOG_C("ELASTIC", Sev::inform, "Set block %d order to %d",
414 it->getMeshsetId(),
415 block_data[it->getMeshsetId()].oRder);
416 Range block_ents;
417 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
418 true);
419 Range ents_to_set_order;
420 CHKERR moab.get_adjacencies(block_ents, 3, false,
421 ents_to_set_order,
422 moab::Interface::UNION);
423 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
424 CHKERR moab.get_adjacencies(block_ents, 2, false,
425 ents_to_set_order,
426 moab::Interface::UNION);
427 CHKERR moab.get_adjacencies(block_ents, 1, false,
428 ents_to_set_order,
429 moab::Interface::UNION);
430 CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
431 ents_to_set_order);
432
433 CHKERR m_field.set_field_order(
434 ents_to_set_order, "DISPLACEMENT",
435 block_data[it->getMeshsetId()].oRder);
436 }
437 std::vector<std::string> additional_parameters;
438 additional_parameters =
439 collect_unrecognized(parsed.options, po::include_positional);
440 for (std::vector<std::string>::iterator vit =
441 additional_parameters.begin();
442 vit != additional_parameters.end(); vit++) {
443 MOFEM_LOG_C("ELASTIC", Sev::warning, "Unrecognized option %s",
444 vit->c_str());
445 }
446 }
447
448 // Update material parameters. Set material parameters block by
449 // block.
451 m_field, BLOCKSET | MAT_ELASTICSET, it)) {
452 const int id = it->getMeshsetId();
453 auto &bd = (*block_sets_ptr)[id];
454 if (block_data[id].yOung > 0)
455 bd.E = block_data[id].yOung;
456 if (block_data[id].pOisson >= -1)
457 bd.PoissonRatio = block_data[id].pOisson;
458 MOFEM_LOG_C("ELASTIC", Sev::inform, "Block %d", id);
459 MOFEM_LOG_C("ELASTIC", Sev::inform, "\tYoung modulus %3.4g", bd.E);
460 MOFEM_LOG_C("ELASTIC", Sev::inform, "\tPoisson ratio %3.4g",
461 bd.PoissonRatio);
462 }
463
465 };
466
467 // Add elastic element
468
469 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
470 boost::make_shared<std::map<int, HookeElement::BlockData>>();
471 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
472 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
473
474 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
475 boost::make_shared<std::map<int, MassBlockData>>();
476 CHKERR ConvectiveMassElement::setBlocks(m_field, mass_block_sets_ptr);
477
478 auto fe_lhs_ptr =
479 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
480 auto fe_rhs_ptr =
481 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
482 fe_lhs_ptr->getRuleHook = VolRule();
483 fe_rhs_ptr->getRuleHook = VolRule();
484
485 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_lhs_ptr, true, false, false,
486 false);
487 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_rhs_ptr, true, false, false,
488 false);
489
490 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
491 new PrismFE(m_field));
492 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
493 new PrismFE(m_field));
494
495 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
496 "DISPLACEMENT",
497 "MESH_NODE_POSITIONS", false);
498
499 auto add_skin_element_for_post_processing = [&]() {
501 Range elastic_element_ents;
503 "ELASTIC", 3, elastic_element_ents);
504 Skinner skin(&m_field.get_moab());
505 Range skin_faces; // skin faces from 3d ents
506 CHKERR skin.find_skin(0, elastic_element_ents, false, skin_faces);
507 Range proc_skin;
508 if (is_partitioned) {
509 CHKERR pcomm->filter_pstatus(skin_faces,
510 PSTATUS_SHARED | PSTATUS_MULTISHARED,
511 PSTATUS_NOT, -1, &proc_skin);
512 } else {
513 proc_skin = skin_faces;
514 }
515 CHKERR m_field.add_finite_element("POST_PROC_SKIN");
516 CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
517 "DISPLACEMENT");
518 CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
519 "DISPLACEMENT");
520 CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
521 "DISPLACEMENT");
523 "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
524 CHKERR m_field.add_ents_to_finite_element_by_dim(proc_skin, 2,
525 "POST_PROC_SKIN");
527 };
528 CHKERR add_skin_element_for_post_processing();
529
530 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
531 if (mesh_has_tets) {
532 CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
533 "DISPLACEMENT", "MESH_NODE_POSITIONS",
534 false, true, MBTET, data_at_pts);
535 }
536 if (mesh_has_prisms) {
538 prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
539 "MESH_NODE_POSITIONS", false, true, MBPRISM, data_at_pts);
540 }
541
542 if (test_nb == 4) {
543
544 auto thermal_strain =
547 constexpr double alpha = 1;
551 t_thermal_strain(i, j) = alpha * t_coords(2) * t_kd(i, j);
552 return t_thermal_strain;
553 };
554
555 fe_rhs_ptr->getOpPtrVector().push_back(
556 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
557 "DISPLACEMENT", data_at_pts, thermal_strain));
558 }
559
560 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
562
563 for (auto &sit : *block_sets_ptr) {
564 for (auto &mit : *mass_block_sets_ptr) {
565 fe_mass_ptr->getOpPtrVector().push_back(
566 new HookeElement::OpCalculateMassMatrix("DISPLACEMENT",
567 "DISPLACEMENT", sit.second,
568 mit.second, data_at_pts));
569 }
570 }
571
572 // Add spring boundary condition applied on surfaces.
573 // This is only declaration not implementation.
574 CHKERR MetaSpringBC::addSpringElements(m_field, "DISPLACEMENT",
575 "MESH_NODE_POSITIONS");
576
577 // Implementation of spring element
578 // Create new instances of face elements for springs
579 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
581 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
583
584 CHKERR MetaSpringBC::setSpringOperators(m_field, fe_spring_lhs_ptr,
585 fe_spring_rhs_ptr, "DISPLACEMENT",
586 "MESH_NODE_POSITIONS");
587
588 // Add Simple Rod elements
589 // This is only declaration not implementation.
591 "MESH_NODE_POSITIONS");
592
593 // CHKERR m_field.add_ents_to_finite_element_by_type(edges_in_simple_rod,
594 // MBEDGE, "SIMPLE_ROD");
595
596 // Implementation of Simple Rod element
597 // Create new instances of edge elements for Simple Rod
598 boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(new EdgeEle(m_field));
599 boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(new EdgeEle(m_field));
600
601
603 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
604 "MESH_NODE_POSITIONS");
605
606 // Add body force element. This is only declaration of element. not its
607 // implementation.
608 CHKERR m_field.add_finite_element("BODY_FORCE");
609 CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
610 "DISPLACEMENT");
611 CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
612 "DISPLACEMENT");
613 CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
614 "DISPLACEMENT");
615 CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
616 "MESH_NODE_POSITIONS");
617
619 m_field, BLOCKSET | BODYFORCESSET, it)) {
620 Range tets;
621 CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
622 true);
623 CHKERR m_field.add_ents_to_finite_element_by_dim(tets, 3, "BODY_FORCE");
624 }
625 CHKERR m_field.build_finite_elements("BODY_FORCE");
626
627 // Add Neumann forces, i.e. pressure or traction forces applied on body
628 // surface. This is only declaration not implementation.
629 CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
630 CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
631 CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
632
633 // Add fluid pressure finite elements. This is special pressure on the
634 // surface from fluid, i.e. pressure which linearly change with the depth.
635 FluidPressure fluid_pressure_fe(m_field);
636 // This function only declare element. Element is implemented by operators
637 // in class FluidPressure.
638 fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
639
640 // Add elements for thermo elasticity if temperature field is defined.
641 ThermalStressElement thermal_stress_elem(m_field);
642 // Check if TEMP field exist, and then add element.
643 if (!m_field.check_field("TEMP")) {
644 bool add_temp_field = false;
646 if (block_data[it->getMeshsetId()].initTemp != 0) {
647 add_temp_field = true;
648 break;
649 }
650 }
651 if (add_temp_field) {
652 CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
653 MB_TAG_SPARSE, MF_ZERO);
654
655 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
656 CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
657 }
658 }
659 if (m_field.check_field("TEMP")) {
660 CHKERR thermal_stress_elem.addThermalStressElement(
661 "ELASTIC", "DISPLACEMENT", "TEMP");
662 }
663
664 // All is declared, at this point build fields first,
665 CHKERR m_field.build_fields();
666 // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
667 // Projection10NodeCoordsOnField
668 // do the trick.
669 Projection10NodeCoordsOnField ent_method_material(m_field,
670 "MESH_NODE_POSITIONS");
671 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
672 if (m_field.check_field("TEMP")) {
674 if (block_data[it->getMeshsetId()].initTemp != 0) {
675 MOFEM_LOG_C("ELASTIC", Sev::inform,
676 "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
677 block_data[it->getMeshsetId()].initTemp);
678 Range block_ents;
679 CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
680 Range vertices;
681 CHKERR moab.get_connectivity(block_ents, vertices, true);
682 CHKERR m_field.getInterface<FieldBlas>()->setField(
683 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
684 "TEMP");
685 }
686 }
687 }
688
689 // Build database for elements. Actual implementation of element is not need
690 // here, only elements has to be declared.
692 // Build adjacencies between elements and field entities
693 CHKERR m_field.build_adjacencies(bit_level0);
694
695 // Register MOFEM DM implementation in PETSc
697
698 // Create DM manager
699 auto dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
700 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
701 CHKERR DMSetFromOptions(dm);
702 CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
703 // Add elements to DM manager
704 CHKERR DMMoFEMAddElement(dm, "ELASTIC");
705 CHKERR DMMoFEMAddElement(dm, "SPRING");
706 CHKERR DMMoFEMAddElement(dm, "SIMPLE_ROD");
707 CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
708 CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
709 CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
710 CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
711 CHKERR DMMoFEMAddElement(dm, "POST_PROC_SKIN");
712 CHKERR DMSetUp(dm);
713
714 // Create matrices & vectors. Note that native PETSc DM interface is used,
715 // but under the PETSc interface MoFEM implementation is running.
718 auto D = smartVectorDuplicate(F);
719 auto D0 = smartVectorDuplicate(F);
722 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
723
724 // Initialise mass matrix
726 if (is_calculating_frequency == PETSC_TRUE) {
727 Mij = smartMatDuplicate(Aij, MAT_DO_NOT_COPY_VALUES);
728 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
729 // MatView(Mij, PETSC_VIEWER_STDOUT_SELF);
730 }
731
732 // Assign global matrix/vector contributed by springs
733 fe_spring_lhs_ptr->ksp_B = Aij;
734 fe_spring_rhs_ptr->ksp_f = F;
735
736 // Assign global matrix/vector contributed by Simple Rod
737 fe_simple_rod_lhs_ptr->ksp_B = Aij;
738 fe_simple_rod_rhs_ptr->ksp_f = F;
739
740 // Zero vectors and matrices
741 CHKERR VecZeroEntries(F);
742 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
743 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
744 CHKERR VecZeroEntries(D);
745 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
746 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
747 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
748 CHKERR MatZeroEntries(Aij);
749
750 // Below particular implementations of finite elements are used to assemble
751 // problem matrixes and vectors. Implementation of element does not change
752 // how element is declared.
753
754 // Assemble Aij and F. Define Dirichlet bc element, which sets constrains
755 // to MatrixDouble and the right hand side vector.
756
757 // if normally defined boundary conditions are not found,
758 // DirichletDisplacementBc will try to use DISPLACEMENT blockset. Two
759 // implementations are available, depending how BC is defined on mesh file.
760 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
761 m_field, "DISPLACEMENT", Aij, D0, F);
762
763 // That sets Dirichlet bc objects that problem is linear, i.e. no newton
764 // (SNES) solver is run for this problem.
765 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
766 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
767
768 // D0 vector will store initial displacements
769 CHKERR VecZeroEntries(D0);
770 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
771 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
772 CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
773 // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
774 // implementation
775 // of particular dirichlet_bc.
776 CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
777 // Set values from D0 on the field (on the mesh)
778
779 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
780 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
781 CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
782
783 // Calculate residual forces as result of applied kinematic constrains. Run
784 // implementation
785 // of particular finite element implementation. Look how
786 // NonlinearElasticElement is implemented,
787 // in that case. We run NonlinearElasticElement with hook material.
788 // Calculate right hand side vector
789 fe_rhs_ptr->snes_f = F;
790 prism_fe_rhs_ptr->snes_f = F;
791 MOFEM_LOG("ELASTIC", Sev::inform) << "Assemble external force vector ...";
792 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
793 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
794 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
795 // Assemble matrix
796 fe_lhs_ptr->snes_B = Aij;
797 prism_fe_lhs_ptr->snes_B = Aij;
798 MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate stiffness matrix ...";
799 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
800 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
801 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
802
803 // Assemble springs
804 CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_lhs_ptr);
805 CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_rhs_ptr);
806
807 // Assemble Simple Rod
808 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_lhs_ptr);
809 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_rhs_ptr);
810
811 if (is_calculating_frequency == PETSC_TRUE) {
812 // Assemble mass matrix
813 fe_mass_ptr->snes_B = Mij;
814 MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate mass matrix ...";
815 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_mass_ptr);
816 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
817 }
818
819 // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
820
821 // Assemble pressure and traction forces. Run particular implemented for do
822 // this, see
823 // MetaNeumannForces how this is implemented.
824 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
826 F, "DISPLACEMENT");
827
828 {
829 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
830 neumann_forces.begin();
831 for (; mit != neumann_forces.end(); mit++) {
832 CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
833 &mit->second->getLoopFe());
834 }
835 }
836 // Assemble forces applied to nodes, see implementation in NodalForce
837 boost::ptr_map<std::string, NodalForce> nodal_forces;
838 CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
839 "DISPLACEMENT");
840
841 {
842 boost::ptr_map<std::string, NodalForce>::iterator fit =
843 nodal_forces.begin();
844 for (; fit != nodal_forces.end(); fit++) {
845 CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
846 &fit->second->getLoopFe());
847 }
848 }
849 // Assemble edge forces
850 boost::ptr_map<std::string, EdgeForce> edge_forces;
851 CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
852 "DISPLACEMENT");
853 {
854 auto fit = edge_forces.begin();
855 for (; fit != edge_forces.end(); fit++) {
856 auto &fe = fit->second->getLoopFe();
857 CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(), &fe);
858 }
859 }
860 // Assemble body forces, implemented in BodyForceConstantField
861 BodyForceConstantField body_forces_methods(m_field);
863 m_field, BLOCKSET | BODYFORCESSET, it)) {
864 CHKERR body_forces_methods.addBlock("DISPLACEMENT", F,
865 it->getMeshsetId());
866 }
867 CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
868 &body_forces_methods.getLoopFe());
869 // Assemble fluid pressure forces
870 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fluid_pressure_fe.getLoopFe(),
871 false, false);
872 CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
873 "DISPLACEMENT", F, false, true);
874
875 CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
876 &fluid_pressure_fe.getLoopFe());
877 // Apply kinematic constrain to right hand side vector and matrix
878 CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
879
880 // Matrix View
881 PetscViewerPushFormat(
882 PETSC_VIEWER_STDOUT_SELF,
883 PETSC_VIEWER_ASCII_MATLAB); /// PETSC_VIEWER_ASCII_DENSE,
884 /// PETSC_VIEWER_ASCII_MATLAB
885 // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
886 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
887 // std::string wait;
888 // std::cin >> wait;
889
890 if (is_calculating_frequency == PETSC_TRUE) {
891 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
892 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
893 }
894
895 // Set matrix positive defined and symmetric for Cholesky and icc
896 // pre-conditioner
897
898 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
899 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
900 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
901 CHKERR VecAssemblyBegin(F);
902 CHKERR VecAssemblyEnd(F);
903 CHKERR VecScale(F, -1);
904
905 // Create solver
906 auto solver = createKSP(PETSC_COMM_WORLD);
907 CHKERR KSPSetDM(solver, dm);
908 CHKERR KSPSetFromOptions(solver);
909 CHKERR KSPSetOperators(solver, Aij, Aij);
910
911 // Setup multi-grid pre-conditioner if set from command line
912 {
913 // from PETSc example ex42.c
914 PetscBool same = PETSC_FALSE;
915 PC pc;
916 CHKERR KSPGetPC(solver, &pc);
917 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
918 if (same) {
919 PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
920 CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
921 CHKERR PCSetFromOptions(pc);
922 } else {
923 // Operators are already set, do not use DM for doing that
924 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
925 }
926 }
927 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
928 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
929 // Set up solver
930 CHKERR KSPSetUp(solver);
931
932 auto set_post_proc_skin = [&](auto &post_proc_skin) {
934 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false,
935 false);
936 CHKERR post_proc_skin.generateReferenceElementMesh();
937 CHKERR post_proc_skin.addFieldValuesPostProc("DISPLACEMENT");
938 CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
939 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
940 "DISPLACEMENT", "ELASTIC", data_at_pts->hMat, true);
941 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
942
943 "MESH_NODE_POSITIONS", "ELASTIC", data_at_pts->HMat, false);
944 post_proc_skin.getOpPtrVector().push_back(
945 new HookeElement::OpPostProcHookeElement<
947 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
948 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
949 true));
951 };
952
953 auto set_post_proc_tets = [&](auto &post_proc) {
955 // Add operators to the elements, starting with some generic operators
956 CHKERR post_proc.generateReferenceElementMesh();
957 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, true, false, false,
958 false);
959 CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
960 CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
961 CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
962 // Add problem specific operator on element to post-process stresses
963 post_proc.getOpPtrVector().push_back(new PostProcHookStress(
964 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
965 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
967 };
968
969 auto set_post_proc_edge = [&](auto &post_proc_edge) {
971 CHKERR post_proc_edge.generateReferenceElementMesh();
972 CHKERR post_proc_edge.addFieldValuesPostProc("DISPLACEMENT");
974 };
975
976 auto set_post_proc_prisms = [&](auto &prism_post_proc) {
978 CHKERR prism_post_proc.generateReferenceElementMesh();
979 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
980 prism_post_proc.getOpPtrVector().push_back(
981 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
982 prism_post_proc.getOpPtrVector().push_back(
983 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
984 CHKERR prism_post_proc.addFieldValuesPostProc("DISPLACEMENT");
985 CHKERR prism_post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
986 CHKERR prism_post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
987 prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
988 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
989 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
991 };
992
993 PostProcFaceOnRefinedMesh post_proc_skin(m_field);
994 PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
995 PostProcEdgeOnRefinedMesh post_proc_edge(m_field);
996 PostProcVolumeOnRefinedMesh post_proc(m_field);
997
998 CHKERR set_post_proc_skin(post_proc_skin);
999 CHKERR set_post_proc_tets(post_proc);
1000 CHKERR set_post_proc_prisms(prism_post_proc);
1001 CHKERR set_post_proc_edge(post_proc_edge);
1002
1003 // Temperature field is defined on the mesh
1004 if (m_field.check_field("TEMP")) {
1005
1006 // Create thermal vector
1007 Vec F_thermal;
1008 CHKERR VecDuplicate(F, &F_thermal);
1009
1010 // Set up implementation for calculation of thermal stress vector. Look
1011 // how thermal stresses and vector is assembled in ThermalStressElement.
1012 CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1013 "DISPLACEMENT", "TEMP", F_thermal);
1014
1015 SeriesRecorder *recorder_ptr;
1016 CHKERR m_field.getInterface(recorder_ptr);
1017
1018 // Read time series and do thermo-elastic analysis, this is when time
1019 // dependent
1020 // temperature problem was run before on the mesh. It means that before
1021 // non-stationary
1022 // problem was solved for temperature and filed "TEMP" is stored for
1023 // subsequent time
1024 // steps in the recorder.
1025 if (recorder_ptr->check_series("THEMP_SERIES")) {
1026 // This is time dependent case, so loop of data series stored by tape
1027 // recorder.
1028 // Loop over time steps
1029 for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
1030 sit)) {
1031 MOFEM_LOG_C("ELASTIC", Sev::inform, "Process step %d",
1032 sit->get_step_number());
1033 // Load field data for this time step
1034 CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
1035 sit->get_step_number());
1036
1037 CHKERR VecZeroEntries(F_thermal);
1038 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1039 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1040
1041 // Calculate the right-hand side vector as result of thermal stresses.
1042 // It MetaNodalForces
1043 // that on "ELASTIC" element data structure the element implementation
1044 // in thermal_stress_elem
1045 // is executed.
1047 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1048
1049 // Assemble vector
1050 CHKERR VecAssemblyBegin(F_thermal);
1051 CHKERR VecAssemblyEnd(F_thermal);
1052 // Accumulate ghost dofs
1053 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1054 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1055
1056 // Calculate norm of vector and print values
1057 PetscReal nrm_F;
1058 CHKERR VecNorm(F, NORM_2, &nrm_F);
1059
1060 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1061 PetscReal nrm_F_thermal;
1062 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1063 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1064 nrm_F_thermal);
1065
1066 CHKERR VecScale(F_thermal, -1);
1067 // check this !!!
1068 CHKERR VecAXPY(F_thermal, 1, F);
1069
1070 // Set dirichlet boundary to thermal stresses vector
1071 dirichlet_bc_ptr->snes_x = D;
1072 dirichlet_bc_ptr->snes_f = F_thermal;
1073 CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1074
1075 // Solve problem
1076 CHKERR KSPSolve(solver, F_thermal, D);
1077 // Add boundary conditions vector
1078 CHKERR VecAXPY(D, 1., D0);
1079 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1080 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1081
1082 // Save data on the mesh
1083 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1084
1085 // Save data on mesh
1086 CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1087
1088
1089 // Post-process results
1090 if (is_post_proc_volume == PETSC_TRUE) {
1091 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1092 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1093 std::ostringstream o1;
1094 o1 << "out_" << sit->step_number << ".h5m";
1095 if (!test_nb)
1096 CHKERR post_proc.writeFile(o1.str().c_str());
1097 MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
1098 }
1099
1100 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1101 CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN",
1102 &post_proc_skin);
1103 std::ostringstream o1_skin;
1104 o1_skin << "out_skin" << sit->step_number << ".h5m";
1105 if (!test_nb)
1106 CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1107 MOFEM_LOG("POST_PROC_SKIN", Sev::inform) << "done ...";
1108 }
1109 } else {
1110
1111 // This is a case when stationary problem for temperature was solved.
1112 CHKERR VecZeroEntries(F_thermal);
1113 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1114 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1115
1116 // Calculate the right-hand side vector with thermal stresses
1118 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1119
1120 // Assemble vector
1121 CHKERR VecAssemblyBegin(F_thermal);
1122 CHKERR VecAssemblyEnd(F_thermal);
1123
1124 // Accumulate ghost dofs
1125 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1126 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1127
1128 // Calculate norm
1129 PetscReal nrm_F;
1130 CHKERR VecNorm(F, NORM_2, &nrm_F);
1131
1132 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1133 PetscReal nrm_F_thermal;
1134 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1135
1136 MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1137 nrm_F_thermal);
1138
1139 // Add thermal stress vector and other forces vector
1140 CHKERR VecScale(F_thermal, -1);
1141 CHKERR VecAXPY(F_thermal, 1, F);
1142
1143 // Apply kinetic boundary conditions
1144 dirichlet_bc_ptr->snes_x = D;
1145 dirichlet_bc_ptr->snes_f = F_thermal;
1146 CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1147
1148 // Solve problem
1149 CHKERR KSPSolve(solver, F_thermal, D);
1150 CHKERR VecAXPY(D, 1., D0);
1151
1152 // Update ghost values for solution vector
1153 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1154 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1155 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1156
1157 // Save data on mesh
1158 if (is_post_proc_volume == PETSC_TRUE) {
1159 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1160 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1161 // Save results to file
1162 if (!test_nb)
1163 CHKERR post_proc.writeFile("out.h5m");
1164 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1165 }
1166
1167 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1168 CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1169 if (!test_nb)
1170 CHKERR post_proc_skin.writeFile("out_skin.h5m");
1171 MOFEM_LOG("POST_PROC_SKIN", Sev::inform) << "done";
1172 }
1173
1174 // Destroy vector, no needed any more
1175 CHKERR VecDestroy(&F_thermal);
1176 } else {
1177 // Elastic analysis no temperature field
1178 // VecView(F, PETSC_VIEWER_STDOUT_WORLD);
1179 // Solve for vector D
1180 CHKERR KSPSolve(solver, F, D);
1181
1182 // VecView(D, PETSC_VIEWER_STDOUT_WORLD);
1183 // cerr << F;
1184
1185 // Add kinetic boundary conditions
1186 CHKERR VecAXPY(D, 1., D0);
1187 // Update ghost values
1188 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1189 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1190 // Save data from vector on mesh
1191 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1192 // Post-process results
1193 MOFEM_LOG("ELASTIC", Sev::inform) << "Post-process start ...";
1194 if (is_post_proc_volume == PETSC_TRUE) {
1195 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1196 }
1197 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
1198 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
1199 CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1200 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1201
1202 // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
1203 MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1204 if (mesh_has_tets) {
1205 if (is_post_proc_volume == PETSC_TRUE) {
1206 if (!test_nb)
1207 CHKERR post_proc.writeFile("out.h5m");
1208 }
1209 if (!test_nb)
1210 CHKERR post_proc_skin.writeFile("out_skin.h5m");
1211 }
1212 if (mesh_has_prisms) {
1213 if (!test_nb)
1214 CHKERR prism_post_proc.writeFile("prism_out.h5m");
1215 }
1216 if (!edges_in_simple_rod.empty())
1217 if (!test_nb)
1218 CHKERR post_proc_edge.writeFile("out_edge.h5m");
1219 MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1220 }
1221
1222 if (is_calculating_frequency == PETSC_TRUE) {
1223 // Calculate mode mass, m = u^T * M * u
1224 Vec u1;
1225 VecDuplicate(D, &u1);
1226 CHKERR MatMult(Mij, D, u1);
1227 double mode_mass;
1228 CHKERR VecDot(u1, D, &mode_mass);
1229 MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode mass %6.4e\n", mode_mass);
1230
1231 Vec v1;
1232 VecDuplicate(D, &v1);
1233 CHKERR MatMult(Aij, D, v1);
1234
1235 double mode_stiffness;
1236 CHKERR VecDot(v1, D, &mode_stiffness);
1237 MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode stiffness %6.4e\n",
1238 mode_stiffness);
1239
1240 double frequency;
1241 double pi = 3.14159265359;
1242 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1243 MOFEM_LOG_C("ELASTIC", Sev::inform, "Frequency %6.4e", frequency);
1244 }
1245
1246 // Calculate elastic energy
1247 auto calculate_strain_energy = [&]() {
1249
1250 SmartPetscObj<Vec> v_energy;
1251 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
1252 "MESH_NODE_POSITIONS", false, true,
1253 v_energy);
1254
1255 // Print elastic energy
1256 double energy;
1257 CHKERR VecSum(v_energy, &energy);
1258 MOFEM_LOG_C("ELASTIC", Sev::inform, "Elastic energy %6.4e", energy);
1259
1260 switch (test_nb) {
1261 case 1:
1262 if (fabs(energy - 17.129) > 1e-3)
1263 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1264 "atom test diverged!");
1265 break;
1266 case 2:
1267 if (fabs(energy - 5.6475e-03) > 1e-4)
1268 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1269 "atom test diverged!");
1270 break;
1271 case 3:
1272 if (fabs(energy - 7.4679e-03) > 1e-4)
1273 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1274 "atom test diverged!");
1275 break;
1276 case 4:
1277 if (fabs(energy - 2.4992e+00) > 1e-3)
1278 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1279 "atom test diverged!");
1280 break;
1281 // FIXME: Here are missing regersion tests
1282 case 8: {
1283 double min;
1284 CHKERR VecMin(D, PETSC_NULL, &min);
1285 constexpr double expected_val = 0.10001;
1286 if (fabs(min + expected_val) > 1e-10)
1287 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1288 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1289 } break;
1290 case 9: {
1291 if (fabs(energy - 4.7416e-04) > 1e-8)
1292 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1293 "atom test diverged!");
1294 }
1295 default:
1296 break;
1297 }
1298
1300 };
1301 CHKERR calculate_strain_energy();
1302
1303 MPI_Comm_free(&moab_comm_world);
1304 }
1306
1308
1309 return 0;
1310}
const std::string default_options
std::string param_file
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
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
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: elasticity.cpp:94
static char help[]
Definition: elasticity.cpp:57
constexpr auto t_kd
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1070
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: DMMMoFEM.cpp:118
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:450
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:503
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1144
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1114
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: DMMMoFEM.cpp:493
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 add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to 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
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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_col(const std::string &fe_name, const std::string &name_row)=0
set field col 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
Add entities to field meshset.
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
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
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
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
static double const pi
Definition: clipper.cpp:53
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
auto createKSP(MPI_Comm comm)
DEPRECATED MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
DEPRECATED MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
Body forces elements.
Definition: BodyForce.hpp:12
Fluid pressure forces.
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
build adjacencies
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.
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
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.
Transform local reference derivatives of shape functions to global derivatives.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
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.
Implentation of thermal stress element.
Set integration rule.

Variable Documentation

◆ help

char help[]
static
Initial value:
= "-my_block_config set block data\n"
"-my_order approximation order\n"
"-my_is_partitioned set if mesh is partitioned\n"
"\n"
Examples
elasticity.cpp.

Definition at line 57 of file elasticity.cpp.