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