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

Go to the source code of this file.

Classes

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

Typedefs

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

Functions

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

Variables

static char help []
 

Typedef Documentation

◆ BlockData

◆ EdgeEle

Examples
elasticity.cpp.

Definition at line 105 of file elasticity.cpp.

◆ MassBlockData

Definition at line 86 of file elasticity.cpp.

◆ VolUserDataOperator

Function Documentation

◆ main()

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

Definition at line 107 of file elasticity.cpp.

107  {
108 
109  const string default_options = "-ksp_type gmres \n"
110  "-pc_type lu \n"
111  "-pc_factor_mat_solver_type mumps \n"
112  "-mat_mumps_icntl_20 0 \n"
113  "-ksp_monitor \n"
114  "-snes_type newtonls \n"
115  "-snes_linesearch_type basic \n"
116  "-snes_atol 1e-8 \n"
117  "-snes_rtol 1e-8 \n"
118  "-snes_monitor \n"
119  "-ts_monitor \n"
120  "-ts_type beuler \n";
121 
122  string param_file = "param_file.petsc";
123  if (!static_cast<bool>(ifstream(param_file))) {
124  std::ofstream file(param_file.c_str(), std::ios::ate);
125  if (file.is_open()) {
126  file << default_options;
127  file.close();
128  }
129  }
130 
131  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
132  auto core_log = logging::core::get();
133  core_log->add_sink(
134  LogManager::createSink(LogManager::getStrmWorld(), "ELASTIC"));
135  LogManager::setLog("ELASTIC");
136  MOFEM_LOG_TAG("ELASTIC", "elasticity")
137 
138  try {
139 
140  PetscBool flg_block_config, flg_file;
141  char mesh_file_name[255];
142  char block_config_file[255];
143  PetscInt test_nb = 0;
144  PetscInt order = 2;
145  PetscBool is_partitioned = PETSC_FALSE;
146  PetscBool is_calculating_frequency = PETSC_FALSE;
147  PetscBool is_post_proc_volume = PETSC_TRUE;
148 
149  // Select base
150  enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
151  const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
152  "jacobi"};
153  PetscInt choice_base_value = LOBATTO;
154 
155  // Read options from command line
156  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
157  CHKERR(ierr);
158  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
159  mesh_file_name, 255, &flg_file);
160 
161  CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
162  order, &order, PETSC_NULL);
163 
164  CHKERR PetscOptionsEList("-base", "approximation base", "", list_bases,
165  LASBASETOP, list_bases[choice_base_value],
166  &choice_base_value, PETSC_NULL);
167 
168  CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "", test_nb,
169  &test_nb, PETSC_NULL);
170 
171  CHKERR PetscOptionsBool("-my_is_partitioned",
172  "set if mesh is partitioned (this result that each "
173  "process keeps only one part of the mesh)",
174  "", is_partitioned, &is_partitioned, PETSC_NULL);
175 
176  CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
177  "", "block_conf.in", block_config_file, 255,
178  &flg_block_config);
179 
180  CHKERR PetscOptionsBool(
181  "-my_is_calculating_frequency", "set if frequency will be calculated",
182  "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
183 
184  CHKERR PetscOptionsBool("-my_is_post_proc_volume",
185  "if true post proc volume", "", is_post_proc_volume,
186  &is_post_proc_volume, PETSC_NULL);
187 
188  ierr = PetscOptionsEnd();
189  CHKERRG(ierr);
190 
191  // Throw error if file with mesh is not provided
192  if (flg_file != PETSC_TRUE) {
193  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
194  }
195 
196  // Create mesh database
197  moab::Core mb_instance;
198  moab::Interface &moab = mb_instance;
199 
200  // Create moab communicator
201  // Create separate MOAB communicator, it is duplicate of PETSc communicator.
202  // NOTE That this should eliminate potential communication problems between
203  // MOAB and PETSC functions.
204  MPI_Comm moab_comm_world;
205  MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
206  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
207  if (pcomm == NULL)
208  pcomm = new ParallelComm(&moab, moab_comm_world);
209 
210  // Read whole mesh or part of it if partitioned
211  if (is_partitioned == PETSC_TRUE) {
212  // This is a case of distributed mesh and algebra. In this case each
213  // processor keeps only one part of the problem.
214  const char *option;
215  option = "PARALLEL=READ_PART;"
216  "PARALLEL_RESOLVE_SHARED_ENTS;"
217  "PARTITION=PARALLEL_PARTITION;";
218  CHKERR moab.load_file(mesh_file_name, 0, option);
219  } else {
220  // If that case we have distributed algebra, i.e. assembly of vectors and
221  // matrices is in parallel, but whole mesh is stored on all processors.
222  // Solver and matrix scales well, however problem set-up of problem is
223  // not fully parallel.
224  const char *option;
225  option = "";
226  CHKERR moab.load_file(mesh_file_name, 0, option);
227  }
228 
229  // Create MoFEM database and link it to MoAB
230  MoFEM::Core core(moab);
231  MoFEM::Interface &m_field = core;
232 
233  // Print boundary conditions and material parameters
234  MeshsetsManager *meshsets_mng_ptr;
235  CHKERR m_field.getInterface(meshsets_mng_ptr);
236  CHKERR meshsets_mng_ptr->printDisplacementSet();
237  CHKERR meshsets_mng_ptr->printForceSet();
238  CHKERR meshsets_mng_ptr->printMaterialsSet();
239 
240  bool mesh_has_tets = false;
241  bool mesh_has_prisms = false;
242  int nb_tets = 0;
243  int nb_hexs = 0;
244  int nb_prisms = 0;
245 
246  CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets, true);
247  CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs, true);
248  CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms, true);
249 
250  mesh_has_tets = (nb_tets + nb_hexs) > 0;
251  mesh_has_prisms = nb_prisms > 0;
252 
253  // Set bit refinement level to all entities (we do not refine mesh in
254  // this example so all entities have the same bit refinement level)
255  BitRefLevel bit_level0;
256  bit_level0.set(0);
257  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
258  0, 3, bit_level0);
259 
260  // CHECK IF EDGE BLOCKSET EXIST AND IF IT IS ADD ALL ENTITIES FROM IT
261  // CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
262  // MESHSET_OF_EDGE_BLOCKSET, 1, bit_level0);
263 
265  if (bit->getName().compare(0, 3, "ROD") == 0) {
266  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
267  0, 1, bit_level0);
268  }
269  }
270 
271  // Declare approximation fields
273  switch (choice_base_value) {
274  case LEGENDRE:
276  break;
277  case LOBATTO:
278  base = AINSWORTH_LOBATTO_BASE;
279  break;
280  case BERNSTEIN_BEZIER:
282  break;
283  case JACOBI:
284  base = DEMKOWICZ_JACOBI_BASE;
285  break;
286  default:
287  SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
288  };
289  CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_DENSE,
290  MF_ZERO);
291 
292  // We can use higher oder geometry to define body
293  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, base, 3, MB_TAG_DENSE,
294  MF_ZERO);
295 
296  // Declare problem
297 
298  // Add entities (by tets) to the field (all entities in the mesh, root_set
299  // = 0 )
300  CHKERR m_field.add_ents_to_field_by_dim(0, 3, "DISPLACEMENT");
301  CHKERR m_field.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
302 
303  // Get all edges in the mesh
304  Range all_edges;
305  CHKERR m_field.get_moab().get_entities_by_type(0, MBEDGE, all_edges, true);
306 
307  // Get edges associated with simple rod
308  Range edges_in_simple_rod;
310  if (bit->getName().compare(0, 3, "ROD") == 0) {
311  Range edges;
312  CHKERR m_field.get_moab().get_entities_by_type(bit->getMeshset(),
313  MBEDGE, edges, true);
314  edges_in_simple_rod.merge(edges);
315  }
316  }
317 
318  CHKERR m_field.add_ents_to_field_by_type(edges_in_simple_rod, MBEDGE,
319  "DISPLACEMENT");
320 
321  // Set order of edge in rod to be 1
322  CHKERR m_field.set_field_order(edges_in_simple_rod, "DISPLACEMENT", 1);
323 
324  // Get remaining edges (not associated with simple rod) to set order
325  Range edges_to_set_order;
326  edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
327 
328  // Set approximation order.
329  // See Hierarchical Finite Element Bases on Unstructured Tetrahedral
330  // Meshes.
331  CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
332  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
333  CHKERR m_field.set_field_order(0, MBHEX, "DISPLACEMENT", order);
334  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
335  CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
336  CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
337  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
338 
340  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
341  else
342  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
343 
344  // Set order of approximation of geometry.
345  // Apply 2nd order only on skin (or in whole body)
346  auto setting_second_order_geometry = [&m_field]() {
348  // Setting geometry order everywhere
349  Range tets, edges;
350  CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
351  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
352  moab::Interface::UNION);
353 
354  // Setting 2nd geometry order only on skin
355  // Range tets, faces, edges;
356  // Skinner skin(&m_field.get_moab());
357  // CHKERR skin.find_skin(0,tets,false,faces);
358  // CHKERR m_field.get_moab().get_adjacencies(
359  // faces,1,false,edges,moab::Interface::UNION
360  // );
361  // CHKERR
362  // m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
363 
364  CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
365  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
366 
368  };
369  CHKERR setting_second_order_geometry();
370 
371  // Configure blocks by parsing config file. It allows setting
372  // approximation order for each block independently.
373  std::map<int, BlockOptionData> block_data;
374  auto setting_blocks_data_and_order_from_config_file =
375  [&m_field, &moab, &block_data, flg_block_config, block_config_file,
376  order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
378  if (flg_block_config) {
379  ifstream ini_file(block_config_file);
380  po::variables_map vm;
381  po::options_description config_file_options;
383  it)) {
384  std::ostringstream str_order;
385  str_order << "block_" << it->getMeshsetId()
386  << ".displacement_order";
387  config_file_options.add_options()(
388  str_order.str().c_str(),
389  po::value<int>(&block_data[it->getMeshsetId()].oRder)
390  ->default_value(order));
391 
392  std::ostringstream str_cond;
393  str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
394  config_file_options.add_options()(
395  str_cond.str().c_str(),
396  po::value<double>(&block_data[it->getMeshsetId()].yOung)
397  ->default_value(-1));
398 
399  std::ostringstream str_capa;
400  str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
401  config_file_options.add_options()(
402  str_capa.str().c_str(),
403  po::value<double>(&block_data[it->getMeshsetId()].pOisson)
404  ->default_value(-2));
405 
406  std::ostringstream str_init_temp;
407  str_init_temp << "block_" << it->getMeshsetId()
408  << ".initial_temperature";
409  config_file_options.add_options()(
410  str_init_temp.str().c_str(),
411  po::value<double>(&block_data[it->getMeshsetId()].initTemp)
412  ->default_value(0));
413  }
414  po::parsed_options parsed =
415  parse_config_file(ini_file, config_file_options, true);
416  store(parsed, vm);
417  po::notify(vm);
419  it)) {
420  if (block_data[it->getMeshsetId()].oRder == -1)
421  continue;
422  if (block_data[it->getMeshsetId()].oRder == order)
423  continue;
424  MOFEM_LOG_C("ELASTIC", Sev::inform, "Set block %d order to %d",
425  it->getMeshsetId(),
426  block_data[it->getMeshsetId()].oRder);
427  Range block_ents;
428  CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
429  true);
430  Range ents_to_set_order;
431  CHKERR moab.get_adjacencies(block_ents, 3, false,
432  ents_to_set_order,
433  moab::Interface::UNION);
434  ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
435  CHKERR moab.get_adjacencies(block_ents, 2, false,
436  ents_to_set_order,
437  moab::Interface::UNION);
438  CHKERR moab.get_adjacencies(block_ents, 1, false,
439  ents_to_set_order,
440  moab::Interface::UNION);
441  CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
442  ents_to_set_order);
443 
444  CHKERR m_field.set_field_order(
445  ents_to_set_order, "DISPLACEMENT",
446  block_data[it->getMeshsetId()].oRder);
447  }
448  std::vector<std::string> additional_parameters;
449  additional_parameters =
450  collect_unrecognized(parsed.options, po::include_positional);
451  for (std::vector<std::string>::iterator vit =
452  additional_parameters.begin();
453  vit != additional_parameters.end(); vit++) {
454  MOFEM_LOG_C("ELASTIC", Sev::warning, "Unrecognized option %s",
455  vit->c_str());
456  }
457  }
458 
459  // Update material parameters. Set material parameters block by
460  // block.
462  m_field, BLOCKSET | MAT_ELASTICSET, it)) {
463  const int id = it->getMeshsetId();
464  auto &bd = (*block_sets_ptr)[id];
465  if (block_data[id].yOung > 0)
466  bd.E = block_data[id].yOung;
467  if (block_data[id].pOisson >= -1)
468  bd.PoissonRatio = block_data[id].pOisson;
469  MOFEM_LOG_C("ELASTIC", Sev::inform, "Block %d", id);
470  MOFEM_LOG_C("ELASTIC", Sev::inform, "\tYoung modulus %3.4g", bd.E);
471  MOFEM_LOG_C("ELASTIC", Sev::inform, "\tPoisson ratio %3.4g",
472  bd.PoissonRatio);
473  }
474 
476  };
477 
478  // Add elastic element
479 
480  boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
481  boost::make_shared<std::map<int, HookeElement::BlockData>>();
482  CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
483  CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
484 
485  boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
486  boost::make_shared<std::map<int, MassBlockData>>();
487  CHKERR ConvectiveMassElement::setBlocks(m_field, mass_block_sets_ptr);
488 
489  auto fe_lhs_ptr =
490  boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
491  auto fe_rhs_ptr =
492  boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
493  fe_lhs_ptr->getRuleHook = VolRule();
494  fe_rhs_ptr->getRuleHook = VolRule();
495 
496  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_lhs_ptr, true, false, false,
497  false);
498  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_rhs_ptr, true, false, false,
499  false);
500 
501  boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
502  new PrismFE(m_field));
503  boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
504  new PrismFE(m_field));
505 
506  CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
507  "DISPLACEMENT",
508  "MESH_NODE_POSITIONS", false);
509 
510  auto add_skin_element_for_post_processing = [&]() {
512  Range elastic_element_ents;
514  "ELASTIC", 3, elastic_element_ents);
515  Skinner skin(&m_field.get_moab());
516  Range skin_faces; // skin faces from 3d ents
517  CHKERR skin.find_skin(0, elastic_element_ents, false, skin_faces);
518  Range proc_skin;
519  if (is_partitioned) {
520  CHKERR pcomm->filter_pstatus(skin_faces,
521  PSTATUS_SHARED | PSTATUS_MULTISHARED,
522  PSTATUS_NOT, -1, &proc_skin);
523  } else {
524  proc_skin = skin_faces;
525  }
526  CHKERR m_field.add_finite_element("POST_PROC_SKIN");
527  CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
528  "DISPLACEMENT");
529  CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
530  "DISPLACEMENT");
531  CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
532  "DISPLACEMENT");
534  "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
535  CHKERR m_field.add_ents_to_finite_element_by_dim(proc_skin, 2,
536  "POST_PROC_SKIN");
538  };
539  CHKERR add_skin_element_for_post_processing();
540 
541  auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
542  if (mesh_has_tets) {
543  CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
544  "DISPLACEMENT", "MESH_NODE_POSITIONS",
545  false, true, MBTET, data_at_pts);
546  }
547  if (mesh_has_prisms) {
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 =
557  FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
558  constexpr double alpha = 1;
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 
614  m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
615  "MESH_NODE_POSITIONS");
616 
617  // Add body force element. This is only declaration of element. not its
618  // implementation.
619  CHKERR m_field.add_finite_element("BODY_FORCE");
620  CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
621  "DISPLACEMENT");
622  CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
623  "DISPLACEMENT");
624  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
625  "DISPLACEMENT");
626  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
627  "MESH_NODE_POSITIONS");
628 
630  m_field, BLOCKSET | BODYFORCESSET, it)) {
631  Range tets;
632  CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
633  true);
634  CHKERR m_field.add_ents_to_finite_element_by_dim(tets, 3, "BODY_FORCE");
635  }
636  CHKERR m_field.build_finite_elements("BODY_FORCE");
637 
638  // Add Neumann forces, i.e. pressure or traction forces applied on body
639  // surface. This is only declaration not implementation.
640  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
641  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
642  CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
643 
644  // Add fluid pressure finite elements. This is special pressure on the
645  // surface from fluid, i.e. pressure which linearly change with the depth.
646  FluidPressure fluid_pressure_fe(m_field);
647  // This function only declare element. Element is implemented by operators
648  // in class FluidPressure.
649  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
650 
651  // Add elements for thermo elasticity if temperature field is defined.
652  ThermalStressElement thermal_stress_elem(m_field);
653  // Check if TEMP field exist, and then add element.
654  if (!m_field.check_field("TEMP")) {
655  bool add_temp_field = false;
657  if (block_data[it->getMeshsetId()].initTemp != 0) {
658  add_temp_field = true;
659  break;
660  }
661  }
662  if (add_temp_field) {
663  CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
664  MB_TAG_SPARSE, MF_ZERO);
665 
666  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
667  CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
668  }
669  }
670  if (m_field.check_field("TEMP")) {
671  CHKERR thermal_stress_elem.addThermalStressElement(
672  "ELASTIC", "DISPLACEMENT", "TEMP");
673  }
674 
675  // All is declared, at this point build fields first,
676  CHKERR m_field.build_fields();
677  // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
678  // Projection10NodeCoordsOnField
679  // do the trick.
680  Projection10NodeCoordsOnField ent_method_material(m_field,
681  "MESH_NODE_POSITIONS");
682  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
683  if (m_field.check_field("TEMP")) {
685  if (block_data[it->getMeshsetId()].initTemp != 0) {
686  MOFEM_LOG_C("ELASTIC", Sev::inform,
687  "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
688  block_data[it->getMeshsetId()].initTemp);
689  Range block_ents;
690  CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
691  Range vertices;
692  CHKERR moab.get_connectivity(block_ents, vertices, true);
693  CHKERR m_field.getInterface<FieldBlas>()->setField(
694  block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
695  "TEMP");
696  }
697  }
698  }
699 
700  // Build database for elements. Actual implementation of element is not need
701  // here, only elements has to be declared.
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
708 
709  // Create DM manager
710  auto dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
711  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
712  CHKERR DMSetFromOptions(dm);
713  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
714  // Add elements to DM manager
715  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
716  CHKERR DMMoFEMAddElement(dm, "SPRING");
717  CHKERR DMMoFEMAddElement(dm, "SIMPLE_ROD");
718  CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
719  CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
720  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
721  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
722  CHKERR DMMoFEMAddElement(dm, "POST_PROC_SKIN");
723  CHKERR DMSetUp(dm);
724 
725  // Create matrices & vectors. Note that native PETSc DM interface is used,
726  // but under the PETSc interface MoFEM implementation is running.
729  auto D = smartVectorDuplicate(F);
730  auto D0 = smartVectorDuplicate(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 = smartMatDuplicate(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);
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  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
931  CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
932  CHKERR PCSetFromOptions(pc);
933  } else {
934  // Operators are already set, do not use DM for doing that
935  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
936  }
937  }
938  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
939  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
940  // Set up solver
941  CHKERR KSPSetUp(solver);
942 
943  auto set_post_proc_skin = [&](auto &post_proc_skin) {
945  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false,
946  false);
947  CHKERR post_proc_skin.generateReferenceElementMesh();
948  CHKERR post_proc_skin.addFieldValuesPostProc("DISPLACEMENT");
949  CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
950  CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
951  "DISPLACEMENT", "ELASTIC", data_at_pts->hMat, true);
952  CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
953 
954  "MESH_NODE_POSITIONS", "ELASTIC", data_at_pts->HMat, false);
955  post_proc_skin.getOpPtrVector().push_back(
956  new HookeElement::OpPostProcHookeElement<
958  "DISPLACEMENT", data_at_pts, *block_sets_ptr,
959  post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
960  true));
962  };
963 
964  auto set_post_proc_tets = [&](auto &post_proc) {
966  // Add operators to the elements, starting with some generic operators
967  CHKERR post_proc.generateReferenceElementMesh();
968  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, true, false, false,
969  false);
970  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
971  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
972  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
973  // Add problem specific operator on element to post-process stresses
974  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
975  m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
976  "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
978  };
979 
980  auto set_post_proc_edge = [&](auto &post_proc_edge) {
982  CHKERR post_proc_edge.generateReferenceElementMesh();
983  CHKERR post_proc_edge.addFieldValuesPostProc("DISPLACEMENT");
985  };
986 
987  auto set_post_proc_prisms = [&](auto &prism_post_proc) {
989  CHKERR prism_post_proc.generateReferenceElementMesh();
990  boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
991  prism_post_proc.getOpPtrVector().push_back(
992  new OpCalculateInvJacForFatPrism(inv_jac_ptr));
993  prism_post_proc.getOpPtrVector().push_back(
994  new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
995  CHKERR prism_post_proc.addFieldValuesPostProc("DISPLACEMENT");
996  CHKERR prism_post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
997  CHKERR prism_post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
998  prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
999  m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
1000  "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
1002  };
1003 
1004  PostProcFaceOnRefinedMesh post_proc_skin(m_field);
1005  PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
1006  PostProcEdgeOnRefinedMesh post_proc_edge(m_field);
1007  PostProcVolumeOnRefinedMesh post_proc(m_field);
1008 
1009  CHKERR set_post_proc_skin(post_proc_skin);
1010  CHKERR set_post_proc_tets(post_proc);
1011  CHKERR set_post_proc_prisms(prism_post_proc);
1012  CHKERR set_post_proc_edge(post_proc_edge);
1013 
1014  // Temperature field is defined on the mesh
1015  if (m_field.check_field("TEMP")) {
1016 
1017  // Create thermal vector
1018  Vec F_thermal;
1019  CHKERR VecDuplicate(F, &F_thermal);
1020 
1021  // Set up implementation for calculation of thermal stress vector. Look
1022  // how thermal stresses and vector is assembled in ThermalStressElement.
1023  CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1024  "DISPLACEMENT", "TEMP", F_thermal);
1025 
1026  SeriesRecorder *recorder_ptr;
1027  CHKERR m_field.getInterface(recorder_ptr);
1028 
1029  // Read time series and do thermo-elastic analysis, this is when time
1030  // dependent
1031  // temperature problem was run before on the mesh. It means that before
1032  // non-stationary
1033  // problem was solved for temperature and filed "TEMP" is stored for
1034  // subsequent time
1035  // steps in the recorder.
1036  if (recorder_ptr->check_series("THEMP_SERIES")) {
1037  // This is time dependent case, so loop of data series stored by tape
1038  // recorder.
1039  // Loop over time steps
1040  for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
1041  sit)) {
1042  MOFEM_LOG_C("ELASTIC", Sev::inform, "Process step %d",
1043  sit->get_step_number());
1044  // Load field data for this time step
1045  CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
1046  sit->get_step_number());
1047 
1048  CHKERR VecZeroEntries(F_thermal);
1049  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1050  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1051 
1052  // Calculate the right-hand side vector as result of thermal stresses.
1053  // It MetaNodalForces
1054  // that on "ELASTIC" element data structure the element implementation
1055  // in thermal_stress_elem
1056  // is executed.
1058  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1059 
1060  // Assemble vector
1061  CHKERR VecAssemblyBegin(F_thermal);
1062  CHKERR VecAssemblyEnd(F_thermal);
1063  // Accumulate ghost dofs
1064  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1065  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1066 
1067  // Calculate norm of vector and print values
1068  PetscReal nrm_F;
1069  CHKERR VecNorm(F, NORM_2, &nrm_F);
1070 
1071  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1072  PetscReal nrm_F_thermal;
1073  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1074  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1075  nrm_F_thermal);
1076 
1077  CHKERR VecScale(F_thermal, -1);
1078  // check this !!!
1079  CHKERR VecAXPY(F_thermal, 1, F);
1080 
1081  // Set dirichlet boundary to thermal stresses vector
1082  dirichlet_bc_ptr->snes_x = D;
1083  dirichlet_bc_ptr->snes_f = F_thermal;
1084  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1085 
1086  // Solve problem
1087  CHKERR KSPSolve(solver, F_thermal, D);
1088  // Add boundary conditions vector
1089  CHKERR VecAXPY(D, 1., D0);
1090  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1091  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1092 
1093  // Save data on the mesh
1094  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1095 
1096  // Save data on mesh
1097  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1098 
1099 
1100  // Post-process results
1101  if (is_post_proc_volume == PETSC_TRUE) {
1102  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1103  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1104  std::ostringstream o1;
1105  o1 << "out_" << sit->step_number << ".h5m";
1106  if (!test_nb)
1107  CHKERR post_proc.writeFile(o1.str().c_str());
1108  MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
1109  }
1110 
1111  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1112  CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN",
1113  &post_proc_skin);
1114  std::ostringstream o1_skin;
1115  o1_skin << "out_skin" << sit->step_number << ".h5m";
1116  if (!test_nb)
1117  CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1118  MOFEM_LOG("POST_PROC_SKIN", Sev::inform) << "done ...";
1119  }
1120  } else {
1121 
1122  // This is a case when stationary problem for temperature was solved.
1123  CHKERR VecZeroEntries(F_thermal);
1124  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1125  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1126 
1127  // Calculate the right-hand side vector with thermal stresses
1129  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1130 
1131  // Assemble vector
1132  CHKERR VecAssemblyBegin(F_thermal);
1133  CHKERR VecAssemblyEnd(F_thermal);
1134 
1135  // Accumulate ghost dofs
1136  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1137  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1138 
1139  // Calculate norm
1140  PetscReal nrm_F;
1141  CHKERR VecNorm(F, NORM_2, &nrm_F);
1142 
1143  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1144  PetscReal nrm_F_thermal;
1145  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1146 
1147  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1148  nrm_F_thermal);
1149 
1150  // Add thermal stress vector and other forces vector
1151  CHKERR VecScale(F_thermal, -1);
1152  CHKERR VecAXPY(F_thermal, 1, F);
1153 
1154  // Apply kinetic boundary conditions
1155  dirichlet_bc_ptr->snes_x = D;
1156  dirichlet_bc_ptr->snes_f = F_thermal;
1157  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1158 
1159  // Solve problem
1160  CHKERR KSPSolve(solver, F_thermal, D);
1161  CHKERR VecAXPY(D, 1., D0);
1162 
1163  // Update ghost values for solution vector
1164  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1165  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1166  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1167 
1168  // Save data on mesh
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  // Save results to file
1173  if (!test_nb)
1174  CHKERR post_proc.writeFile("out.h5m");
1175  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1176  }
1177 
1178  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1179  CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1180  if (!test_nb)
1181  CHKERR post_proc_skin.writeFile("out_skin.h5m");
1182  MOFEM_LOG("POST_PROC_SKIN", Sev::inform) << "done";
1183  }
1184 
1185  // Destroy vector, no needed any more
1186  CHKERR VecDestroy(&F_thermal);
1187  } else {
1188  // Elastic analysis no temperature field
1189  // VecView(F, PETSC_VIEWER_STDOUT_WORLD);
1190  // Solve for vector D
1191  CHKERR KSPSolve(solver, F, D);
1192 
1193  // VecView(D, PETSC_VIEWER_STDOUT_WORLD);
1194  // cerr << F;
1195 
1196  // Add kinetic boundary conditions
1197  CHKERR VecAXPY(D, 1., D0);
1198  // Update ghost values
1199  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1200  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1201  // Save data from vector on mesh
1202  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1203  // Post-process results
1204  MOFEM_LOG("ELASTIC", Sev::inform) << "Post-process start ...";
1205  if (is_post_proc_volume == PETSC_TRUE) {
1206  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1207  }
1208  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
1209  CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
1210  CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1211  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1212 
1213  // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
1214  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1215  if (mesh_has_tets) {
1216  if (is_post_proc_volume == PETSC_TRUE) {
1217  if (!test_nb)
1218  CHKERR post_proc.writeFile("out.h5m");
1219  }
1220  if (!test_nb)
1221  CHKERR post_proc_skin.writeFile("out_skin.h5m");
1222  }
1223  if (mesh_has_prisms) {
1224  if (!test_nb)
1225  CHKERR prism_post_proc.writeFile("prism_out.h5m");
1226  }
1227  if (!edges_in_simple_rod.empty())
1228  if (!test_nb)
1229  CHKERR post_proc_edge.writeFile("out_edge.h5m");
1230  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1231  }
1232 
1233  if (is_calculating_frequency == PETSC_TRUE) {
1234  // Calculate mode mass, m = u^T * M * u
1235  Vec u1;
1236  VecDuplicate(D, &u1);
1237  CHKERR MatMult(Mij, D, u1);
1238  double mode_mass;
1239  CHKERR VecDot(u1, D, &mode_mass);
1240  MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode mass %6.4e\n", mode_mass);
1241 
1242  Vec v1;
1243  VecDuplicate(D, &v1);
1244  CHKERR MatMult(Aij, D, v1);
1245 
1246  double mode_stiffness;
1247  CHKERR VecDot(v1, D, &mode_stiffness);
1248  MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode stiffness %6.4e\n",
1249  mode_stiffness);
1250 
1251  double frequency;
1252  double pi = 3.14159265359;
1253  frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1254  MOFEM_LOG_C("ELASTIC", Sev::inform, "Frequency %6.4e", frequency);
1255  }
1256 
1257  // Calculate elastic energy
1258  auto calculate_strain_energy = [&]() {
1260 
1261  SmartPetscObj<Vec> v_energy;
1262  CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
1263  "MESH_NODE_POSITIONS", false, true,
1264  v_energy);
1265 
1266  // Print elastic energy
1267  double energy;
1268  CHKERR VecSum(v_energy, &energy);
1269  MOFEM_LOG_C("ELASTIC", Sev::inform, "Elastic energy %6.4e", energy);
1270 
1271  switch (test_nb) {
1272  case 1:
1273  if (fabs(energy - 17.129) > 1e-3)
1274  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1275  "atom test diverged!");
1276  break;
1277  case 2:
1278  if (fabs(energy - 5.6475e-03) > 1e-4)
1279  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1280  "atom test diverged!");
1281  break;
1282  case 3:
1283  if (fabs(energy - 7.4679e-03) > 1e-4)
1284  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1285  "atom test diverged!");
1286  break;
1287  case 4:
1288  if (fabs(energy - 2.4992e+00) > 1e-3)
1289  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1290  "atom test diverged!");
1291  break;
1292  // FIXME: Here are missing regersion tests
1293  case 8: {
1294  double min;
1295  CHKERR VecMin(D, PETSC_NULL, &min);
1296  constexpr double expected_val = 0.10001;
1297  if (fabs(min + expected_val) > 1e-10)
1298  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1299  "atom test diverged! %3.4e != %3.4e", min, expected_val);
1300  } break;
1301  case 9: {
1302  if (fabs(energy - 4.7416e-04) > 1e-8)
1303  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1304  "atom test diverged!");
1305  }
1306  default:
1307  break;
1308  }
1309 
1311  };
1312  CHKERR calculate_strain_energy();
1313 
1314  MPI_Comm_free(&moab_comm_world);
1315  }
1316  CATCH_ERRORS;
1317 
1319 
1320  return 0;
1321 }
const std::string default_options
std::string param_file
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
Kronecker Delta class symmetric.
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ NOBASE
Definition: definitions.h:72
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:175
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: elasticity.cpp:105
static char help[]
Definition: elasticity.cpp:68
constexpr auto t_kd
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
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:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:514
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:481
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1135
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1105
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:504
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
virtual bool check_series(const std::string &name) const
check if series is in database
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
static double const pi
Definition: clipper.cpp:53
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
auto createSmartDM
Creates smart DM object.
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
auto createKSP
CoreTmp< 0 > Core
Definition: Core.hpp:1096
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
Body forces elements.
Definition: BodyForce.hpp:24
Fluid pressure forces.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:109
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
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:140
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
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
static MoFEMErrorCode setSimpleRodOperators(MoFEM::Interface &m_field, boost::shared_ptr< EdgeElementForcesAndSourcesCoreBase > fe_simple_rod_lhs_ptr, boost::shared_ptr< EdgeElementForcesAndSourcesCoreBase > fe_simple_rod_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of SimpleRod element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Managing BitRefLevels.
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:31
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
Transform local reference derivatives of shape functions to global derivatives.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Set data structures of MG pre-conditioner via approximation orders.
Postprocess on edge.
Postprocess on face.
Postprocess on prism.
Operator post-procesing stresses for Hook isotropic material.
Post processing.
Implentation of thermal stress element.
Set integration rule.

Variable Documentation

◆ help

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

Definition at line 68 of file elasticity.cpp.