v0.9.2
Classes | Typedefs | Functions | Variables
#include <BasicFiniteElements.hpp>
#include <boost/program_options.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes

struct  BlockOptionData
 stores some optional data More...
 
struct  VolRule
 Set integration rule. More...
 
struct  PrismFE
 

Typedefs

using BlockData = NonlinearElasticElement::BlockData
 
using MassBlockData = ConvectiveMassElement::BlockData
 
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
 

Functions

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

Variables

static char help []
 

Typedef Documentation

◆ BlockData

◆ MassBlockData

Definition at line 90 of file elasticity.cpp.

◆ VolUserDataOperator

Function Documentation

◆ main()

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

Definition at line 109 of file elasticity.cpp.

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

Variable Documentation

◆ help

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

Definition at line 72 of file elasticity.cpp.