v0.8.13
Classes | 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
 

Functions

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

Variables

static char help []
 

Function Documentation

◆ main()

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

Definition at line 81 of file elasticity.cpp.

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

Variable Documentation

◆ help

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

Definition at line 70 of file elasticity.cpp.