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

Go to the source code of this file.

Classes

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

Typedefs

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

Functions

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

Variables

static char help []
 

Typedef Documentation

◆ BlockData

◆ MassBlockData

Definition at line 90 of file elasticity.cpp.

◆ VolUserDataOperator

Function Documentation

◆ main()

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

Definition at line 109 of file elasticity.cpp.

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

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.