v0.10.0
Classes | Typedefs | Functions | Variables
elasticity.cpp File Reference
User modules » NonLinear Elastic ElementUser modules » | NonLinear Elastic Element
#include <BasicFiniteElements.hpp>
#include <boost/program_options.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes

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

Typedefs

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

Functions

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

Variables

static char help []
 

Typedef Documentation

◆ BlockData

using BlockData = NonlinearElasticElement::BlockData
Examples
elasticity.cpp, testing_jacobian_of_hook_element.cpp, and testing_jacobian_of_hook_scaled_with_density_element.cpp.

Definition at line 89 of file elasticity.cpp.

◆ MassBlockData

using MassBlockData = ConvectiveMassElement::BlockData

Definition at line 90 of file elasticity.cpp.

◆ VolUserDataOperator

using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
Examples
testing_jacobian_of_hook_scaled_with_density_element.cpp.

Definition at line 91 of file elasticity.cpp.

Function Documentation

◆ main()

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

Definition at line 109 of file elasticity.cpp.

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

Variable Documentation

◆ help

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

Definition at line 72 of file elasticity.cpp.

Generated on Sun Nov 22 2020 11:26:51 for MoFEM by Doxygen 1.8.15 and hosted at University of Glasgow