v0.7.24
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 80 of file elasticity.cpp.

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