v0.8.19
Classes | Typedefs | Functions | Variables
elasticity.cpp File Reference
#include <BasicFiniteElements.hpp>
#include <boost/program_options.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes

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

Typedefs

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

Functions

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

Variables

static char help []
 

Typedef Documentation

◆ BlockData

◆ VolUserDataOperator

Definition at line 85 of file elasticity.cpp.

Function Documentation

◆ main()

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

Definition at line 92 of file elasticity.cpp.

92  {
93 
94  // Initialize PETSCc
95  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
96 
97  try {
98 
99  PetscBool flg_block_config, flg_file;
100  char mesh_file_name[255];
101  char block_config_file[255];
102  PetscInt order = 2;
103  PetscBool is_partitioned = PETSC_FALSE;
104 
105  // Read options from command line
106  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
107  CHKERR(ierr);
108  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
109  mesh_file_name, 255, &flg_file);
110 
111  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
112  &order, PETSC_NULL);
113 
114  CHKERR PetscOptionsBool("-my_is_partitioned",
115  "set if mesh is partitioned (this result that each "
116  "process keeps only part of the mes",
117  "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
118 
119  CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
120  "", "block_conf.in", block_config_file, 255,
121  &flg_block_config);
122 
123  ierr = PetscOptionsEnd();
124  CHKERRG(ierr);
125 
126  // Throw error if file with mesh is not give
127  if (flg_file != PETSC_TRUE) {
128  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
129  }
130 
131  // Create mesh database
132  moab::Core mb_instance;
133  moab::Interface &moab = mb_instance;
134 
135  // Create moab communicator
136  // Create separate MOAB communicator, it is duplicate of PETSc communicator.
137  // NOTE That this should eliminate potential communication problems between
138  // MOAB and PETSC functions.
139  MPI_Comm moab_comm_world;
140  MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
141  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
142  if (pcomm == NULL)
143  pcomm = new ParallelComm(&moab, moab_comm_world);
144 
145  // Read whole mesh or part of is if partitioned
146  if (is_partitioned == PETSC_TRUE) {
147  // This is a case of distributed mesh and algebra. In that case each
148  // processor
149  // keep only part of the problem.
150  const char *option;
151  option = "PARALLEL=READ_PART;"
152  "PARALLEL_RESOLVE_SHARED_ENTS;"
153  "PARTITION=PARALLEL_PARTITION;";
154  CHKERR moab.load_file(mesh_file_name, 0, option);
155  } else {
156  // If that case we have distributed algebra, i.e. assembly of vectors and
157  // matrices is in parallel, but whole mesh is stored on all processors.
158  // Solver and matrix scales well, however problem set-up of problem is
159  // not fully parallel.
160  const char *option;
161  option = "";
162  CHKERR moab.load_file(mesh_file_name, 0, option);
163  }
164 
165  // Create MoFEM database and link it to MoAB
166  MoFEM::Core core(moab);
167  MoFEM::Interface &m_field = core;
168 
169  // Print boundary conditions and material parameters
170  MeshsetsManager *meshsets_mng_ptr;
171  CHKERR m_field.getInterface(meshsets_mng_ptr);
172  CHKERR meshsets_mng_ptr->printDisplacementSet();
173  CHKERR meshsets_mng_ptr->printForceSet();
174  CHKERR meshsets_mng_ptr->printMaterialsSet();
175 
176  // Set bit refinement level to all entities (we do not refine mesh in this
177  // example
178  // so all entities have the same bit refinement level)
179  BitRefLevel bit_level0;
180  bit_level0.set(0);
181  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
182  0, 3, bit_level0);
183 
184  // Declare approximation fields
185  CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LOBATTO_BASE, 3,
186  MB_TAG_DENSE, MF_ZERO);
187 
188  // We can use higher oder geometry to define body
189  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
190  3, MB_TAG_DENSE, MF_ZERO);
191 
192  // Declare problem
193 
194  // Add entities (by tets) to the field ( all entities in the mesh, root_set
195  // = 0 )
196  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
197  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
198 
199  // Set apportion order.
200  // See Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes.
201  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
202  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
203  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
204  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
205 
206  // Set order of approximation of geometry.
207  // Apply 2nd order only on skin (or in in whole body)
208  auto setting_second_order_geometry = [&m_field]() {
210  // Setting geometry order everywhere
211  Range tets, edges;
212  CHKERR m_field.get_moab().get_entities_by_type(0, MBTET, tets);
213  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
214  moab::Interface::UNION);
215 
216  // Setting 2nd geometry order only on skin
217  // Range tets, faces, edges;
218  // Skinner skin(&m_field.get_moab());
219  // CHKERR skin.find_skin(0,tets,false,faces);
220  // CHKERR m_field.get_moab().get_adjacencies(
221  // faces,1,false,edges,moab::Interface::UNION
222  // );
223  // CHKERR m_field.synchronise_entities(edges);
224 
225  CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
226  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
228  };
229  CHKERR setting_second_order_geometry();
230 
231  // Configure blocks by parsing config file. It allows setting approximation
232  // order for each block independently.
233  std::map<int, BlockOptionData> block_data;
234  auto setting_blocks_data_and_order_from_config_file =
235  [&m_field, &moab, &block_data, flg_block_config, block_config_file,
236  order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
238  if (flg_block_config) {
239  ifstream ini_file(block_config_file);
240  po::variables_map vm;
241  po::options_description config_file_options;
243  it)) {
244  std::ostringstream str_order;
245  str_order << "block_" << it->getMeshsetId()
246  << ".displacement_order";
247  config_file_options.add_options()(
248  str_order.str().c_str(),
249  po::value<int>(&block_data[it->getMeshsetId()].oRder)
250  ->default_value(order));
251 
252  std::ostringstream str_cond;
253  str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
254  config_file_options.add_options()(
255  str_cond.str().c_str(),
256  po::value<double>(&block_data[it->getMeshsetId()].yOung)
257  ->default_value(-1));
258 
259  std::ostringstream str_capa;
260  str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
261  config_file_options.add_options()(
262  str_capa.str().c_str(),
263  po::value<double>(&block_data[it->getMeshsetId()].pOisson)
264  ->default_value(-2));
265 
266  std::ostringstream str_init_temp;
267  str_init_temp << "block_" << it->getMeshsetId()
268  << ".initial_temperature";
269  config_file_options.add_options()(
270  str_init_temp.str().c_str(),
271  po::value<double>(&block_data[it->getMeshsetId()].initTemp)
272  ->default_value(0));
273  }
274  po::parsed_options parsed =
275  parse_config_file(ini_file, config_file_options, true);
276  store(parsed, vm);
277  po::notify(vm);
279  it)) {
280  if (block_data[it->getMeshsetId()].oRder == -1)
281  continue;
282  if (block_data[it->getMeshsetId()].oRder == order)
283  continue;
284  PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
285  it->getMeshsetId(),
286  block_data[it->getMeshsetId()].oRder);
287  Range block_ents;
288  CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
289  true);
290  Range ents_to_set_order;
291  CHKERR moab.get_adjacencies(block_ents, 3, false,
292  ents_to_set_order,
293  moab::Interface::UNION);
294  ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
295  CHKERR moab.get_adjacencies(block_ents, 2, false,
296  ents_to_set_order,
297  moab::Interface::UNION);
298  CHKERR moab.get_adjacencies(block_ents, 1, false,
299  ents_to_set_order,
300  moab::Interface::UNION);
301  CHKERR m_field.synchronise_entities(ents_to_set_order);
302 
303  CHKERR m_field.set_field_order(
304  ents_to_set_order, "DISPLACEMENT",
305  block_data[it->getMeshsetId()].oRder);
306  }
307  std::vector<std::string> additional_parameters;
308  additional_parameters =
309  collect_unrecognized(parsed.options, po::include_positional);
310  for (std::vector<std::string>::iterator vit =
311  additional_parameters.begin();
312  vit != additional_parameters.end(); vit++) {
313  CHKERR PetscPrintf(PETSC_COMM_WORLD,
314  "** WARNING Unrecognized option %s\n",
315  vit->c_str());
316  }
317  }
318 
319  // Update material parameters. Set material parameters block by block.
321  m_field, BLOCKSET | MAT_ELASTICSET, it)) {
322  const int id = it->getMeshsetId();
323  auto &bd = (*block_sets_ptr)[id];
324  if (block_data[id].yOung > 0)
325  bd.E = block_data[id].yOung;
326  if (block_data[id].pOisson >= -1)
327  bd.PoissonRatio = block_data[id].pOisson;
328  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Block %d\n", id);
329  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung modulus %3.4g\n",
330  bd.E);
331  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson ratio %3.4g\n",
332  bd.PoissonRatio);
333  }
334 
336  };
337 
338  // Add elastic element
339 
340  boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
341  boost::make_shared<std::map<int, BlockData>>();
342  CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
343  CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
344 
345  boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
346  new VolumeElementForcesAndSourcesCore(m_field));
347  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
348  new VolumeElementForcesAndSourcesCore(m_field));
349  fe_lhs_ptr->getRuleHook = VolRule();
350  fe_rhs_ptr->getRuleHook = VolRule();
351 
352  CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
353  "DISPLACEMENT",
354  "MESH_NODE_POSITIONS", false);
355  CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
356  "DISPLACEMENT", "MESH_NODE_POSITIONS",
357  false, true);
358 
359  // Add body force element. This is only declaration of element. not its
360  // implementation.
361  CHKERR m_field.add_finite_element("BODY_FORCE");
362  CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
363  "DISPLACEMENT");
364  CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
365  "DISPLACEMENT");
366  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
367  "DISPLACEMENT");
368  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
369  "MESH_NODE_POSITIONS");
370 
372  m_field, BLOCKSET | BODYFORCESSET, it)) {
373  Range tets;
374  m_field.get_moab().get_entities_by_type(it->meshset, MBTET, tets, true);
375  CHKERR
376  m_field.add_ents_to_finite_element_by_type(tets, MBTET, "BODY_FORCE");
377  }
378 
379  // Add Neumann forces, i.e. pressure or traction forces applied on body
380  // surface. This is only declaration not implementation.
381  CHKERR MetaNeummanForces::addNeumannBCElements(m_field, "DISPLACEMENT");
382  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
383  CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
384 
385  // Add fluid pressure finite elements. This is special pressure on the
386  // surface from fluid, i.e. pressure which linearly change with the depth.
387  FluidPressure fluid_pressure_fe(m_field);
388  // This function only declare element. Element is implemented by operators
389  // in class FluidPressure.
390  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
391 
392  // Add elements for thermo elasticity if temperature field is defined.
393  ThermalStressElement thermal_stress_elem(m_field);
394  // Check if TEMP field exist, and then add element.
395  if (!m_field.check_field("TEMP")) {
396  bool add_temp_field = false;
398  if (block_data[it->getMeshsetId()].initTemp != 0) {
399  add_temp_field = true;
400  break;
401  }
402  }
403  if (add_temp_field) {
404  CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
405  MB_TAG_SPARSE, MF_ZERO);
406 
407  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
408 
409  CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
410  }
411  }
412  if (m_field.check_field("TEMP")) {
413  CHKERR thermal_stress_elem.addThermalStressElement(
414  "ELASTIC", "DISPLACEMENT", "TEMP");
415  }
416 
417  // All is declared, at that point build fields first,
418  CHKERR m_field.build_fields();
419  // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
420  // Projection10NodeCoordsOnField
421  // do the trick.
422  Projection10NodeCoordsOnField ent_method_material(m_field,
423  "MESH_NODE_POSITIONS");
424  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
425  if (m_field.check_field("TEMP")) {
427  if (block_data[it->getMeshsetId()].initTemp != 0) {
428  PetscPrintf(PETSC_COMM_WORLD, "Set block %d temperature to %3.2g\n",
429  it->getMeshsetId(),
430  block_data[it->getMeshsetId()].initTemp);
431  Range block_ents;
432  CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
433  Range vertices;
434  CHKERR moab.get_connectivity(block_ents, vertices, true);
435  CHKERR m_field.getInterface<FieldBlas>()->setField(
436  block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
437  "TEMP");
438  }
439  }
440  }
441 
442  // Build database for elements. Actual implementation of element is not need
443  // here, only elements has to be declared.
444  CHKERR m_field.build_finite_elements();
445  // Build adjacencies between elements and field entities
446  CHKERR m_field.build_adjacencies(bit_level0);
447 
448  // Register MOFEM DM implementation in PETSc
450 
451  // Create DM manager
452  DM dm;
453  CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
454  CHKERR DMSetType(dm, "MOFEM");
455  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
456  CHKERR DMSetFromOptions(dm);
457  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
458  // Add elements to DM manager
459  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
460  CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
461  CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
462  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
463  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
464  CHKERR DMSetUp(dm);
465 
466  // Create matrices & vectors. Note that native PETSc DM interface is used,
467  // but under the PETSc interface MOFEM implementation is running.
468  Vec F, D, D0;
469  CHKERR DMCreateGlobalVector(dm, &F);
470  CHKERR VecDuplicate(F, &D);
471  CHKERR VecDuplicate(F, &D0);
472  Mat Aij;
473  CHKERR DMCreateMatrix(dm, &Aij);
474  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
475 
476  // Zero vectors and matrices
477  CHKERR VecZeroEntries(F);
478  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
479  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
480  CHKERR VecZeroEntries(D);
481  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
482  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
483  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
484  CHKERR MatZeroEntries(Aij);
485 
486  // This controls how kinematic constrains are set, by blockset or nodeset.
487  // Cubit
488  // sets kinetic boundary conditions by blockset.
489  bool flag_cubit_disp = false;
491  m_field, NODESET | DISPLACEMENTSET, it)) {
492  flag_cubit_disp = true;
493  }
494 
495  // Below particular implementations of finite elements are used to assemble
496  // problem matrixes and vectors. Implementation of element does not change
497  // how element is declared.
498 
499  // Assemble Aij and F. Define Dirichlet bc element, which stets constrains
500  // to MatrixDouble and the right hand side vector.
501  boost::shared_ptr<FEMethod> dirichlet_bc_ptr;
502 
503  // if normally defined boundary conditions are not found, try to use
504  // DISPLACEMENT blockset. To implementations available here, depending how
505  // BC is defined on mesh file.
506  if (!flag_cubit_disp) {
507  dirichlet_bc_ptr =
508  boost::shared_ptr<FEMethod>(new DirichletSetFieldFromBlockWithFlags(
509  m_field, "DISPLACEMENT", "DISPLACEMENT", Aij, D0, F));
510  } else {
511  dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(
512  new DirichletDisplacementBc(m_field, "DISPLACEMENT", Aij, D0, F));
513  }
514  // That sets Dirichlet bc objects that problem is linear, i.e. no newton
515  // (SNES) solver is run for this problem.
516  dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
517  dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
518 
519  // D0 vector will store initial displacements
520  CHKERR VecZeroEntries(D0);
521  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
522  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
523  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
524  // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
525  // implementation
526  // of particular dirichlet_bc.
527  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
528  // Set values from D0 on the field (on the mesh)
529  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
530 
531  // Calculate residual forces as result of applied kinematic constrains. Run
532  // implementation
533  // of particular finite element implementation. Look how
534  // NonlinearElasticElement is implemented,
535  // in that case. We run NonlinearElasticElement with hook material.
536  // Calculate right hand side vector
537  fe_rhs_ptr->snes_f = F;
538  PetscPrintf(PETSC_COMM_WORLD, "Assemble external force vector ...");
539  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
540  PetscPrintf(PETSC_COMM_WORLD, " done\n");
541  // Assemble matrix
542  fe_lhs_ptr->snes_B = Aij;
543  PetscPrintf(PETSC_COMM_WORLD, "Calculate stiffness matrix ...");
544  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
545  PetscPrintf(PETSC_COMM_WORLD, " done\n");
546 
547  // Assemble pressure and traction forces. Run particular implemented for do
548  // this, see
549  // MetaNeummanForces how this is implemented.
550  boost::ptr_map<std::string, NeummanForcesSurface> neumann_forces;
551  CHKERR MetaNeummanForces::setMomentumFluxOperators(m_field, neumann_forces,
552  F, "DISPLACEMENT");
553 
554  {
555  boost::ptr_map<std::string, NeummanForcesSurface>::iterator mit =
556  neumann_forces.begin();
557  for (; mit != neumann_forces.end(); mit++) {
558  CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
559  &mit->second->getLoopFe());
560  }
561  }
562  // Assemble forces applied to nodes, see implementation in NodalForce
563  boost::ptr_map<std::string, NodalForce> nodal_forces;
564  CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
565  "DISPLACEMENT");
566 
567  {
568  boost::ptr_map<std::string, NodalForce>::iterator fit =
569  nodal_forces.begin();
570  for (; fit != nodal_forces.end(); fit++) {
571  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
572  &fit->second->getLoopFe());
573  }
574  }
575  // Assemble edge forces
576  boost::ptr_map<std::string, EdgeForce> edge_forces;
577  CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
578  "DISPLACEMENT");
579  {
580  boost::ptr_map<std::string, EdgeForce>::iterator fit =
581  edge_forces.begin();
582  for (; fit != edge_forces.end(); fit++) {
583  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
584  &fit->second->getLoopFe());
585  }
586  }
587  // Assemble body forces, implemented in BodyForceConstantField
588  BodyForceConstantField body_forces_methods(m_field);
590  m_field, BLOCKSET | BODYFORCESSET, it)) {
591  CHKERR body_forces_methods.addBlock("DISPLACEMENT", F,
592  it->getMeshsetId());
593  }
594  CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
595  &body_forces_methods.getLoopFe());
596  // Assemble fluid pressure forces
597  CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
598  "DISPLACEMENT", F, false, true);
599 
600  CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
601  &fluid_pressure_fe.getLoopFe());
602  // Apply kinematic constrain to right hand side vector and matrix
603  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
604 
605  // Matrix View
606  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
607  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
608  // std::string wait;
609  // std::cin >> wait;
610 
611  // Set matrix positive defined and symmetric for Cholesky and icc
612  // pre-conditioner
613  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
614  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
615  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
616  CHKERR VecAssemblyBegin(F);
617  CHKERR VecAssemblyEnd(F);
618  CHKERR VecScale(F, -1);
619 
620  // Create solver
621  KSP solver;
622  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
623  CHKERR KSPSetDM(solver, dm);
624  CHKERR KSPSetFromOptions(solver);
625  CHKERR KSPSetOperators(solver, Aij, Aij);
626  // Setup multi-grid pre-conditioner if set from command line
627  {
628  // from PETSc example ex42.c
629  PetscBool same = PETSC_FALSE;
630  PC pc;
631  CHKERR KSPGetPC(solver, &pc);
632  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
633  if (same) {
634  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
635  CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
636  CHKERR PCSetFromOptions(pc);
637  } else {
638  // Operators are already set, do not use DM for doing that
639  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
640  }
641  }
642  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
643  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
644  // Set up solver
645  CHKERR KSPSetUp(solver);
646 
647  // Set up post-processor. It is some generic implementation of finite
648  // element.
649  PostProcVolumeOnRefinedMesh post_proc(m_field);
650  // Add operators to the elements, starting with some generic operators
651  CHKERR post_proc.generateReferenceElementMesh();
652  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
653  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
654  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
655  // Add problem specific operator on element to post-process stresses
656  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
657  m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "DISPLACEMENT",
658  post_proc.commonData, block_sets_ptr.get()));
659 
660  // Temperature field is defined on the mesh
661  if (m_field.check_field("TEMP")) {
662 
663  // Create thermal vector
664  Vec F_thermal;
665  CHKERR VecDuplicate(F, &F_thermal);
666 
667  // Set up implementation for calculation of thermal stress vector. Look
668  // how thermal stresses and vector is assembled in ThermalStressElement.
669  CHKERR thermal_stress_elem.setThermalStressRhsOperators(
670  "DISPLACEMENT", "TEMP", F_thermal);
671 
672  SeriesRecorder *recorder_ptr;
673  CHKERR m_field.getInterface(recorder_ptr);
674 
675  // Read time series and do thermo-elastic analysis, this is when time
676  // dependent
677  // temperature problem was run before on the mesh. It means that before
678  // non-stationary
679  // problem was solved for temperature and filed "TEMP" is stored for
680  // subsequent time
681  // steps in the recorder.
682  if (recorder_ptr->check_series("THEMP_SERIES")) {
683  // This is time dependent case, so loop of data series stored by tape
684  // recorder.
685  // Loop over time steps
686  for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
687  sit)) {
688  PetscPrintf(PETSC_COMM_WORLD, "Process step %d\n",
689  sit->get_step_number());
690  // Load field data for this time step
691  CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
692  sit->get_step_number());
693 
694  CHKERR VecZeroEntries(F_thermal);
695  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
696  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
697 
698  // Calculate the right-hand side vector as result of thermal stresses.
699  // It MetaNodalForces
700  // that on "ELASTIC" element data structure the element implementation
701  // in thermal_stress_elem
702  // is executed.
704  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
705 
706  // Assemble vector
707  CHKERR VecAssemblyBegin(F_thermal);
708  CHKERR VecAssemblyEnd(F_thermal);
709  // Accumulate ghost dofs
710  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
711  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
712 
713  // Calculate norm of vector and print values
714  PetscReal nrm_F;
715  CHKERR VecNorm(F, NORM_2, &nrm_F);
716 
717  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
718  PetscReal nrm_F_thermal;
719  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
720  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
721  nrm_F_thermal);
722 
723  CHKERR VecScale(F_thermal, -1);
724  // check this !!!
725  CHKERR VecAXPY(F_thermal, 1, F);
726 
727  // Set dirichlet boundary to thermal stresses vector
728  dirichlet_bc_ptr->snes_x = D;
729  dirichlet_bc_ptr->snes_f = F_thermal;
730  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
731 
732  // Solve problem
733  CHKERR KSPSolve(solver, F_thermal, D);
734  // Add boundary conditions vector
735  CHKERR VecAXPY(D, 1., D0);
736  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
737  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
738 
739  // Save data on the mesh
740  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
741 
742  // Save data on mesh
743  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
744  // Post-process results
745  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
746 
747  std::ostringstream o1;
748  o1 << "out_" << sit->step_number << ".h5m";
749  CHKERR post_proc.writeFile(o1.str().c_str());
750  }
751  } else {
752 
753  // This is a case when stationary problem for temperature was solved.
754  CHKERR VecZeroEntries(F_thermal);
755  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
756  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
757 
758  // Calculate the right-hand side vector with thermal stresses
760  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
761 
762  // Assemble vector
763  CHKERR VecAssemblyBegin(F_thermal);
764  CHKERR VecAssemblyEnd(F_thermal);
765 
766  // Accumulate ghost dofs
767  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
768  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
769 
770  // Calculate norm
771  PetscReal nrm_F;
772  CHKERR VecNorm(F, NORM_2, &nrm_F);
773 
774  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
775  PetscReal nrm_F_thermal;
776  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
777 
778  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
779  nrm_F_thermal);
780 
781  // Add thermal stress vector and other forces vector
782  CHKERR VecScale(F_thermal, -1);
783  CHKERR VecAXPY(F_thermal, 1, F);
784 
785  // Apply kinetic boundary conditions
786  dirichlet_bc_ptr->snes_x = D;
787  dirichlet_bc_ptr->snes_f = F_thermal;
788  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
789 
790  // Solve problem
791  CHKERR KSPSolve(solver, F_thermal, D);
792  CHKERR VecAXPY(D, 1., D0);
793 
794  // Update ghost values for solution vector
795  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
796  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
797 
798  // Save data on mesh
799  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
800  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
801  // Save results to file
802  PetscPrintf(PETSC_COMM_WORLD, "Write output file ...");
803  CHKERR post_proc.writeFile("out.h5m");
804  PetscPrintf(PETSC_COMM_WORLD, " done\n");
805  }
806 
807  // Destroy vector, no needed any more
808  CHKERR VecDestroy(&F_thermal);
809 
810  } else {
811  // Elastic analysis no temperature field
812  // Solve for vector D
813  CHKERR KSPSolve(solver, F, D);
814  // Add kinetic boundary conditions
815  CHKERR VecAXPY(D, 1., D0);
816  // Update ghost values
817  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
818  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
819  // Save data from vector on mesh
820  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
821  // Post-process results
822  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
823  // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
824  PetscPrintf(PETSC_COMM_WORLD, "Write output file ..,");
825  CHKERR post_proc.writeFile("out.h5m");
826  PetscPrintf(PETSC_COMM_WORLD, " done\n");
827  }
828 
829  // Calculate elastic energy
830  auto calculate_strain_energy = [dm, &block_sets_ptr]() {
832 
833  Vec v_energy;
835  dm, block_sets_ptr, "DISPLACEMENT", "MESH_NODE_POSITIONS", false,
836  true, &v_energy);
837 
838  // Print elastic energy
839  const double *eng_ptr;
840  CHKERR VecGetArrayRead(v_energy, &eng_ptr);
841  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n", *eng_ptr);
842  CHKERR VecRestoreArrayRead(v_energy, &eng_ptr);
843 
844  CHKERR VecDestroy(&v_energy);
845 
847  };
848  CHKERR calculate_strain_energy();
849 
850  // Destroy matrices, vecors, solver and DM
851  CHKERR VecDestroy(&F);
852  CHKERR VecDestroy(&D);
853  CHKERR VecDestroy(&D0);
854  CHKERR MatDestroy(&Aij);
855  CHKERR KSPDestroy(&solver);
856  CHKERR DMDestroy(&dm);
857 
858  MPI_Comm_free(&moab_comm_world);
859  }
860  CATCH_ERRORS;
861 
863 
864  return 0;
865 }
Set integration rule.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeummanForcesSurface > &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 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.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Implentation of thermal stress element.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode 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:453
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:433
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)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
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)
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:542
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
Volume finite elementUser is implementing own operator at Gauss point level, by class derived from Vo...
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:515
static char help[]
Definition: elasticity.cpp:73
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
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.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
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
virtual bool check_series(const std::string &name) const
check if series is in database
Body forces elements.
Definition: BodyForce.hpp:25
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.
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:43
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:143
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
Post processing.
block name is "BODY_FORCES"
Definition: definitions.h:223
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:476
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:486
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:148
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
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:594
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:220
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:283
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode synchronise_entities(Range &ent, 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
virtual bool check_field(const std::string &name) const =0
check if field is in database
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:405
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:431
continuous field
Definition: definitions.h:169
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:846
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:61

Variable Documentation

◆ help

char help[]
static
Initial value:
= "-my_block_config set block data\n"
"\n"
Examples:
elasticity.cpp.

Definition at line 73 of file elasticity.cpp.