v0.8.15
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_SPARSE, 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_SPARSE, 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.
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)
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:452
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:432
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:459
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:526
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:514
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:142
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:222
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:475
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:485
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:147
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:578
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:219
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:282
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:403
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:429
continuous field
Definition: definitions.h:168
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.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:845
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 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.
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.