v0.14.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 94 of file elasticity.cpp.

◆ MassBlockData

Definition at line 75 of file elasticity.cpp.

◆ VolUserDataOperator

Function Documentation

◆ main()

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

PETSC_VIEWER_ASCII_DENSE, PETSC_VIEWER_ASCII_MATLAB

Examples
elasticity.cpp.

Definition at line 96 of file elasticity.cpp.

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

Variable Documentation

◆ help

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

Definition at line 57 of file elasticity.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:789
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MetaNodalForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
MoFEM::SeriesRecorder::check_series
virtual bool check_series(const std::string &name) const
check if series is in database
Definition: SeriesRecorder.cpp:301
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::SeriesRecorder
Definition: SeriesRecorder.hpp:25
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
PCMGSetUpViaApproxOrdersCtx
Set data structures of MG pre-conditioner via approximation orders.
Definition: PCMGSetUpViaApproxOrders.hpp:190
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
VolRule
Set integration rule.
Definition: simple_interface.cpp:88
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
PostProcFatPrismOnRefinedMesh
Postprocess on prism.
Definition: PostProcOnRefMesh.hpp:973
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
ClipperLib::pi
static const double pi
Definition: clipper.cpp:89
calculateEnergy
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)
MoFEM::CoreInterface::add_ents_to_field_by_dim
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3150
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MetaNodalForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:301
PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:727
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MetaSimpleRodElement::setSimpleRodOperators
static MoFEMErrorCode setSimpleRodOperators(MoFEM::Interface &m_field, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_lhs_ptr, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of SimpleRod element. Set operators to calculate LHS and RHS.
Definition: SimpleRodElement.cpp:298
BodyForceConstantField
Body forces elements.
Definition: BodyForce.hpp:12
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3177
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: elasticity.cpp:94
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
BODYFORCESSET
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
FluidPressure
Fluid pressure forces.
Definition: FluidPressure.hpp:20
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::SeriesRecorder::load_series_data
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
Definition: SeriesRecorder.cpp:309
_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
Definition: SeriesRecorder.hpp:205
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
PrismFE
Definition: prism_elements_from_surface.cpp:62
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
PostProcEdgeOnRefinedMesh
Postprocess on edge.
Definition: PostProcOnRefMesh.hpp:1148
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:234
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:114
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:322
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
FTensor::Index< 'i', 3 >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::CoreInterface::get_finite_element_entities_by_dimension
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
MetaEdgeForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
Range
ConvectiveMassElement::setBlocks
MoFEMErrorCode setBlocks()
Definition: ConvectiveMassElement.cpp:1761
DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:375
MetaSpringBC::setSpringOperators
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
Definition: SpringElement.cpp:1178
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:287
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
MetaSimpleRodElement::addSimpleRodElements
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
Definition: SimpleRodElement.cpp:270
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
addElasticElement
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)
Definition: HookeElement.cpp:533
setOperators
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)
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::SmartPetscObj< Vec >
help
static char help[]
Definition: elasticity.cpp:57
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MF_EXIST
@ MF_EXIST
Definition: definitions.h:100
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
ThermalStressElement
Implentation of thermal stress element.
Definition: ThermalStressElement.hpp:15
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::LogManager
Log manager is used to build and partition problems.
Definition: LogManager.hpp:26
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
F
@ F
Definition: free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MetaEdgeForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97