v0.15.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
elasticity.cpp File Reference
User modules » NonLinear Elastic Element
#include <BasicFiniteElements.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes

struct  BlockOptionData
 
struct  VolRule
 Set integration rule to volume elements. 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

using BlockData = NonlinearElasticElement::BlockData

Definition at line 74 of file elasticity.cpp.

◆ EdgeEle

using EdgeEle = MoFEM::EdgeElementForcesAndSourcesCore

Definition at line 94 of file elasticity.cpp.

◆ MassBlockData

using MassBlockData = ConvectiveMassElement::BlockData

Definition at line 75 of file elasticity.cpp.

◆ VolUserDataOperator

using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator

Definition at line 76 of file elasticity.cpp.

Function Documentation

◆ main()

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

PETSC_VIEWER_ASCII_DENSE, PETSC_VIEWER_ASCII_MATLAB

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 core_log->add_sink(
128 LogManager::createSink(LogManager::getStrmSync(), "ELASTIC_SYNC"));
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
256 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
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
264 FieldApproximationBase base = NOBASE;
265 switch (choice_base_value) {
266 case LEGENDRE:
267 base = AINSWORTH_LEGENDRE_BASE;
268 break;
269 case LOBATTO:
270 base = AINSWORTH_LOBATTO_BASE;
271 break;
272 case BERNSTEIN_BEZIER:
273 base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
274 break;
275 case JACOBI:
276 base = DEMKOWICZ_JACOBI_BASE;
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;
301 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
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
331 if (base == AINSWORTH_BERNSTEIN_BEZIER_BASE)
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]() {
339 MoFEMFunctionBegin;
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
359 MoFEMFunctionReturn(0);
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) {
369 MoFEMFunctionBegin;
370 if (flg_block_config) {
371 ifstream ini_file(block_config_file);
372 po::variables_map vm;
373 po::options_description config_file_options;
374 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET,
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);
410 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET,
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.
453 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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
467 MoFEMFunctionReturn(0);
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 = [&]() {
503 MoFEMFunctionBegin;
504 Range elastic_element_ents;
505 CHKERR m_field.get_finite_element_entities_by_dimension(
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 }
533 CHKERR m_field.modify_finite_element_add_field_data(
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");
537 MoFEMFunctionReturn(0);
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 =
556 [](FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords) {
557 FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
558 constexpr double alpha = 1;
559 FTensor::Index<'i', 3> i;
560 FTensor::Index<'k', 3> j;
561 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
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(
572 new VolumeElementForcesAndSourcesCore(m_field));
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(
591 new FaceElementForcesAndSourcesCore(m_field));
592 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
593 new FaceElementForcesAndSourcesCore(m_field));
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.
601 CHKERR MetaSimpleRodElement::addSimpleRodElements(m_field, "DISPLACEMENT",
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
613 CHKERR MetaSimpleRodElement::setSimpleRodOperators(
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
629 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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;
656 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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")) {
684 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
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.
702 CHKERR m_field.build_finite_elements();
703 // Build adjacencies between elements and field entities
704 CHKERR m_field.build_adjacencies(bit_level0);
705
706 // Register MOFEM DM implementation in PETSc
707 CHKERR DMRegister_MGViaApproxOrders("MOFEM");
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.
727 SmartPetscObj<Vec> F;
728 CHKERR DMCreateGlobalVector_MoFEM(dm, F);
729 auto D = vectorDuplicate(F);
730 auto D0 = vectorDuplicate(F);
731 SmartPetscObj<Mat> Aij;
732 CHKERR DMCreateMatrix_MoFEM(dm, Aij);
733 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
734
735 // Initialise mass matrix
736 SmartPetscObj<Mat> Mij;
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;
836 CHKERR MetaNeumannForces::setMomentumFluxOperators(m_field, 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);
873 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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);
883 CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
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) {
930 CHKERR PCMGSetUpViaApproxOrders(
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) {
944 MoFEMFunctionBegin;
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<
960 FaceElementForcesAndSourcesCore>(
961 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
962 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
963 true));
964 MoFEMFunctionReturn(0);
965 };
966
967 auto set_post_proc_tets = [&](auto &post_proc) {
968 MoFEMFunctionBegin;
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()));
984 MoFEMFunctionReturn(0);
985 };
986
987 auto set_post_proc_edge = [&](auto &post_proc_edge) {
988 MoFEMFunctionBegin;
989 CHKERR post_proc_edge.generateReferenceElementMesh();
990 CHKERR post_proc_edge.addFieldValuesPostProc("DISPLACEMENT");
991 MoFEMFunctionReturn(0);
992 };
993
994 auto set_post_proc_prisms = [&](auto &prism_post_proc) {
995 MoFEMFunctionBegin;
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()));
1008 MoFEMFunctionReturn(0);
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.
1098 CHKERR DMoFEMLoopFiniteElements(
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) {
1141 CHKERR m_field.getInterface<FieldEvaluatorInterface>()
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
1165 MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
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
1196 CHKERR DMoFEMLoopFiniteElements(
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 = [&]() {
1327 MoFEMFunctionBegin;
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
1378 MoFEMFunctionReturn(0);
1379 };
1380 CHKERR calculate_strain_energy();
1381
1382 MPI_Comm_free(&moab_comm_world);
1383 }
1384 CATCH_ERRORS;
1385
1386 MoFEM::Core::Finalize();
1387
1388 return 0;
1389}
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition LogManager.hpp:353
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:319
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition Kronecker_Delta.hpp:49
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
QUIET
@ QUIET
Definition definitions.h:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MF_ZERO
@ MF_ZERO
Definition definitions.h:111
MF_EXIST
@ MF_EXIST
Definition definitions.h:113
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
NOBASE
@ NOBASE
Definition definitions.h:59
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
H1
@ H1
continuous field
Definition definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
BODYFORCESSET
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition definitions.h:175
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition definitions.h:172
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
order
constexpr int order
Definition dg_projection.cpp:19
EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition elasticity.cpp:94
help
static char help[]
Definition elasticity.cpp:57
F
@ F
Definition free_surface.cpp:624
t_kd
constexpr auto t_kd
Definition free_surface.cpp:163
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
MoFEM::DMMoFEMCreateMoFEM
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
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
MoFEM::DMoFEMMeshToLocalVector
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
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1187
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition PCMGSetUpViaApproxOrders.cpp:302
MoFEM::DMoFEMLoopFiniteElements
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
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition DMMoFEM.cpp:1157
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::add_ents_to_field_by_dim
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.
MoFEM::CoreInterface::get_finite_element_entities_by_dimension
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
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:392
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
Print meshsets with force boundary conditions.
Definition MeshsetsManager.cpp:315
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition MeshsetsManager.hpp:72
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
Print meshsets with material properties.
Definition MeshsetsManager.cpp:336
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
Print meshsets with displacement boundary conditions.
Definition MeshsetsManager.cpp:301
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition MeshsetsManager.hpp:94
MoFEM::SeriesRecorder::load_series_data
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
Definition SeriesRecorder.cpp:309
MoFEM::SeriesRecorder::check_series
virtual bool check_series(const std::string &name) const
check if series is in database
Definition SeriesRecorder.cpp:301
_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
Definition SeriesRecorder.hpp:205
bit
auto bit
set bit
Definition hanging_node_approx.cpp:75
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
D
double D
Definition initial_diffusion.cpp:44
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition PetscSmartObj.hpp:269
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition HODataOperators.hpp:1011
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::createPCMGSetUpViaApproxOrdersCtx
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
Definition PCMGSetUpViaApproxOrders.cpp:653
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition DeprecatedPetsc.hpp:192
MoFEM::PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
Definition PCMGSetUpViaApproxOrders.cpp:657
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:135
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition PetscSmartObj.hpp:236
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition PetscSmartObj.hpp:143
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition HODataOperators.hpp:986
BodyForceConstantField
Body forces elements.
Definition BodyForce.hpp:12
FluidPressure
Fluid pressure forces.
Definition FluidPressure.hpp:20
MetaEdgeForces::setOperators
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
MetaEdgeForces::addElement
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
MetaNeumannForces::addNeumannBCElements
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.
Definition SurfacePressure.cpp:1974
MetaNeumannForces::setMomentumFluxOperators
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.
Definition SurfacePressure.cpp:2069
MetaNodalForces::setOperators
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
MetaNodalForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition NodalForce.hpp:92
MetaSimpleRodElement::addSimpleRodElements
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
Definition SimpleRodElement.cpp:269
MetaSimpleRodElement::setSimpleRodOperators
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.
Definition SimpleRodElement.cpp:297
MetaSpringBC::setSpringOperators
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
Definition SpringElement.cpp:1178
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition SpringElement.cpp:1127
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CommInterface
Managing BitRefLevels.
Definition CommInterface.hpp:21
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::FieldBlas
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition FieldEvaluator.hpp:21
MoFEM::LogManager
Log manager is used to build and partition problems.
Definition LogManager.hpp:26
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition LogManager.cpp:293
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition LogManager.cpp:339
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:208
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition UserDataOperators.hpp:3983
MoFEM::OpCalculateScalarFieldValues
Specialization for double precision scalar field values calculation.
Definition UserDataOperators.hpp:204
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1745
MoFEM::OpCalculateVectorFieldValues
Specialization for double precision vector field values calculation.
Definition UserDataOperators.hpp:626
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition UserDataOperators.hpp:4010
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition PostProcHookStresses.hpp:41
ThermalStressElement
Implentation of thermal stress element.
Definition ThermalStressElement.hpp:15
VolRule
Set integration rule to volume elements.
Definition simple_interface.cpp:88

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"

Definition at line 57 of file elasticity.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow