v0.6.10
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 79 of file elasticity.cpp.

79  {
80 
81  // Initialize PETSCc
82  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
83 
84  try {
85 
86  PetscBool flg_block_config, flg_file;
87  char mesh_file_name[255];
88  char block_config_file[255];
89  PetscInt order = 2;
90  PetscBool is_partitioned = PETSC_FALSE;
91 
92  // Read options from command line
93  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
94  CHKERR(ierr);
95  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
96  mesh_file_name, 255, &flg_file);
97 
98  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
99  &order, PETSC_NULL);
100 
101  CHKERR PetscOptionsBool("-my_is_partitioned",
102  "set if mesh is partitioned (this result that each "
103  "process keeps only part of the mes",
104  "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
105 
106  CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
107  "", "block_conf.in", block_config_file, 255,
108  &flg_block_config);
109 
110  ierr = PetscOptionsEnd();
111  CHKERRG(ierr);
112 
113  // Throw error if file with mesh is not give
114  if (flg_file != PETSC_TRUE) {
115  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
116  }
117 
118  // Create mesh database
119  moab::Core mb_instance;
120  moab::Interface &moab = mb_instance;
121 
122  // Create moab communicator
123  // Create separate MOAB communicator, it is duplicate of PETSc communicator.
124  // NOTE That this should eliminate potential communication problems between
125  // MOAB and PETSC functions.
126  MPI_Comm moab_comm_world;
127  MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
128  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
129  if (pcomm == NULL)
130  pcomm = new ParallelComm(&moab, moab_comm_world);
131 
132  // Read whole mesh or part of is if partitioned
133  if (is_partitioned == PETSC_TRUE) {
134  // This is a case of distributed mesh and algebra. In that case each
135  // processor
136  // keep only part of the problem.
137  const char *option;
138  option = "PARALLEL=READ_PART;"
139  "PARALLEL_RESOLVE_SHARED_ENTS;"
140  "PARTITION=PARALLEL_PARTITION;";
141  CHKERR moab.load_file(mesh_file_name, 0, option);
142  } else {
143  // If that case we have distributed algebra, i.e. assembly of vectors and
144  // matrices is in parallel, but whole mesh is stored on all processors.
145  // Solver and matrix scales well, however problem set-up of problem is
146  // not fully parallel.
147  const char *option;
148  option = "";
149  CHKERR moab.load_file(mesh_file_name, 0, option);
150  }
151 
152  // Create MoFEM database and link it to MoAB
153  MoFEM::Core core(moab);
154  MoFEM::Interface &m_field = core;
155 
156  // Print boundary conditions and material parameters
157  MeshsetsManager *meshsets_mng_ptr;
158  CHKERR m_field.getInterface(meshsets_mng_ptr);
159  CHKERR meshsets_mng_ptr->printDisplacementSet();
160  CHKERR meshsets_mng_ptr->printForceSet();
161  CHKERR meshsets_mng_ptr->printMaterialsSet();
162 
163  // Set bit refinement level to all entities (we do not refine mesh in this
164  // example
165  // so all entities have the same bit refinement level)
166  BitRefLevel bit_level0;
167  bit_level0.set(0);
168  CHKERR m_field.seed_ref_level_3D(0, bit_level0);
169 
170  // Declare approximation fields
171  CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LOBATTO_BASE, 3,
172  MB_TAG_SPARSE, MF_ZERO);
173 
174  // We can use higher oder geometry to define body
175  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
176  3, MB_TAG_SPARSE, MF_ZERO);
177 
178  // Declare problem
179 
180  // Add entities (by tets) to the field ( all entities in the mesh, root_set
181  // = 0 )
182  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
183  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
184 
185  // Set apportion order.
186  // See Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes.
187  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
188  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
189  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
190  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
191 
192  // Set order of approximation of geometry.
193  // Apply 2nd order only on skin (or in in whole body)
194  {
195  Skinner skin(&m_field.get_moab());
196  Range faces, tets;
197  CHKERR m_field.get_moab().get_entities_by_type(0, MBTET, tets);
198  // CHKERR skin.find_skin(0,tets,false,faces); CHKERRG(rval);
199  Range edges;
200  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
201  moab::Interface::UNION);
202  // CHKERR m_field.get_moab().get_adjacencies(
203  // faces,1,false,edges,moab::Interface::UNION
204  // ); CHKERRG(rval);
205  // CHKERR m_field.synchronise_entities(edges);
206  CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
207 
208  }
209  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
210 
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, true);
264  CHKERRG(rval);
265  Range ents_to_set_order;
266  CHKERR moab.get_adjacencies(block_ents, 3, false, ents_to_set_order,
267  moab::Interface::UNION);
268  ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
269  CHKERR moab.get_adjacencies(block_ents, 2, false, ents_to_set_order,
270  moab::Interface::UNION);
271  CHKERR moab.get_adjacencies(block_ents, 1, false, ents_to_set_order,
272  moab::Interface::UNION);
273  CHKERR m_field.synchronise_entities(ents_to_set_order);
274 
275  CHKERR m_field.set_field_order(ents_to_set_order, "DISPLACEMENT",
276  block_data[it->getMeshsetId()].oRder);
277 
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", vit->c_str());
287  }
288  } catch (const std::exception &ex) {
289  std::ostringstream ss;
290  ss << ex.what() << std::endl;
291  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
292  }
293  }
294 
295  // Add elastic element
296  boost::shared_ptr<Hooke<adouble> > hooke_adouble_ptr(new Hooke<adouble>());
297  boost::shared_ptr<Hooke<double> > hooke_double_ptr(new Hooke<double>());
298  NonlinearElasticElement elastic(m_field, 2);
299  CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
300 
301  CHKERR elastic.addElement("ELASTIC", "DISPLACEMENT");
302  CHKERR
303  elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false, true);
304 
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, "Block %d set Young modulus %3.4g\n",
312  id, elastic.setOfBlocks[id].E);
313 
314  }
315  if (block_data[id].pOisson >= -1) {
316  elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
317  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Block %d set Poisson ratio %3.4g\n",
318  id, elastic.setOfBlocks[id].PoissonRatio);
319 
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
327  m_field.modify_finite_element_add_field_row("BODY_FORCE", "DISPLACEMENT");
328  CHKERR
329  m_field.modify_finite_element_add_field_col("BODY_FORCE", "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 surface. This
344  // 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 surface
350  // 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 in
353  // 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  }
377  if (m_field.check_field("TEMP")) {
378  CHKERR thermal_stress_elem.addThermalStressElement("ELASTIC",
379  "DISPLACEMENT", "TEMP");
380 
381  }
382 
383  // All is declared, at that point build fields first,
384  CHKERR m_field.build_fields();
385  // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
386  // Projection10NodeCoordsOnField
387  // do the trick.
388  Projection10NodeCoordsOnField ent_method_material(m_field,
389  "MESH_NODE_POSITIONS");
390  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
391  if (m_field.check_field("TEMP")) {
393  if (block_data[it->getMeshsetId()].initTemp != 0) {
394  PetscPrintf(PETSC_COMM_WORLD, "Set block %d temperature to %3.2g\n",
395  it->getMeshsetId(),
396  block_data[it->getMeshsetId()].initTemp);
397  Range block_ents;
398  CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
399  Range vertices;
400  CHKERR moab.get_connectivity(block_ents, vertices, true);
401  CHKERR m_field.set_field(block_data[it->getMeshsetId()].initTemp,
402  MBVERTEX, vertices, "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, but
432  // 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 to MatrixDouble
465  // 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 DISPLACEMENT blockset.
469  // To implementations available here, depending how BC is defined on mesh file.
470  if(!flag_cubit_disp){
471  dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(new
473  m_field,"DISPLACEMENT","DISPLACEMENT",Aij,D0,F
474  )
475  );
476  } else {
477  dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(
478  new DirichletDisplacementBc(m_field,"DISPLACEMENT",Aij,D0,F)
479  );
480  }
481  // That sets Dirichlet bc objects that problem is linear, i.e. no newton (SNES) solver
482  // is run for this problem.
483  dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
484  dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
485 
486  // D0 vector will store initial displacements
487  CHKERR VecZeroEntries(D0);
488  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
489  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
490  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
491  // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
492  // implementation
493  // of particular dirichlet_bc.
494  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
495  // Set values from D0 on the field (on the mesh)
496  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
497 
498  // Calculate residual forces as result of applied kinematic constrains. Run
499  // implementation
500  // of particular finite element implementation. Look how
501  // NonlinearElasticElement is implemented,
502  // in that case. We run NonlinearElasticElement with hook material.
503  // Calculate right hand side vector
504  elastic.getLoopFeRhs().snes_f = F;
505  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeRhs());
506  // Assemble matrix
507  elastic.getLoopFeLhs().snes_B = Aij;
508  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeLhs());
509 
510  // Assemble pressure and traction forces. Run particular implemented for do
511  // this, see
512  // MetaNeummanForces how this is implemented.
513  boost::ptr_map<std::string, NeummanForcesSurface> neumann_forces;
514  CHKERR MetaNeummanForces::setMomentumFluxOperators(m_field, neumann_forces, F,
515  "DISPLACEMENT");
516 
517  {
518  boost::ptr_map<std::string, NeummanForcesSurface>::iterator mit =
519  neumann_forces.begin();
520  for (; mit != neumann_forces.end(); mit++) {
521  CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
522  &mit->second->getLoopFe());
523 
524  }
525  }
526  // Assemble forces applied to nodes, see implementation in NodalForce
527  boost::ptr_map<std::string, NodalForce> nodal_forces;
528  CHKERR
529  MetaNodalForces::setOperators(m_field, nodal_forces, F, "DISPLACEMENT");
530 
531  {
532  boost::ptr_map<std::string, NodalForce>::iterator fit =
533  nodal_forces.begin();
534  for (; fit != nodal_forces.end(); fit++) {
535  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
536  &fit->second->getLoopFe());
537  }
538  }
539  // Assemble edge forces
540  boost::ptr_map<std::string, EdgeForce> edge_forces;
541  CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F, "DISPLACEMENT");
542  {
543  boost::ptr_map<std::string, EdgeForce>::iterator fit = edge_forces.begin();
544  for (; fit != edge_forces.end(); fit++) {
545  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
546  &fit->second->getLoopFe());
547  }
548  }
549  // Assemble body forces, implemented in BodyForceConstantField
550  BodyForceConstantField body_forces_methods(m_field);
552  m_field, BLOCKSET | BODYFORCESSET, it)) {
553  CHKERR body_forces_methods.addBlock("DISPLACEMENT", F, it->getMeshsetId());
554 
555  }
556  CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
557  &body_forces_methods.getLoopFe());
558  // Assemble fluid pressure forces
559  CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
560  "DISPLACEMENT", F, false, true);
561 
562  CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
563  &fluid_pressure_fe.getLoopFe());
564  // Apply kinematic constrain to right hand side vector and matrix
565  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
566 
567 
568  // Matrix View
569  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
570  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
571  // std::string wait;
572  // std::cin >> wait;
573 
574  // Set matrix positive defined and symmetric for Cholesky and icc
575  // pre-conditioner
576  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
577  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
578  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
579  CHKERR VecAssemblyBegin(F);
580  CHKERR VecAssemblyEnd(F);
581  CHKERR VecScale(F, -1);
582 
583  // Create solver
584  KSP solver;
585  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
586  CHKERR KSPSetDM(solver, dm);
587  CHKERR KSPSetFromOptions(solver);
588  CHKERR KSPSetOperators(solver, Aij, Aij);
589  // Setup multi-grid pre-conditioner if set from command line
590  {
591  // from PETSc example ex42.c
592  PetscBool same = PETSC_FALSE;
593  PC pc;
594  CHKERR KSPGetPC(solver, &pc);
595  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
596  if (same) {
597  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
598  CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
599  CHKERR PCSetFromOptions(pc);
600  } else {
601  // Operators are already set, do not use DM for doing that
602  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
603  }
604  }
605  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
606  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
607  // Set up solver
608  CHKERR KSPSetUp(solver);
609 
610  // Set up post-processor. It is some generic implementation of finite element.
611  PostProcVolumeOnRefinedMesh post_proc(m_field);
612  // Add operators to the elements, starting with some generic operators
613  CHKERR post_proc.generateReferenceElementMesh();
614  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
615  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
616  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
617  // Add problem specific operator on element to post-process stresses
618  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
619  m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "DISPLACEMENT",
620  post_proc.commonData, &elastic.setOfBlocks));
621 
622  // Temperature field is defined on the mesh
623  if (m_field.check_field("TEMP")) {
624 
625  // Create thermal vector
626  Vec F_thermal;
627  CHKERR VecDuplicate(F, &F_thermal);
628 
629  // Set up implementation for calculation of thermal stress vector. Look how
630  // thermal stresses and vector is assembled in ThermalStressElement.
631  CHKERR thermal_stress_elem.setThermalStressRhsOperators("DISPLACEMENT",
632  "TEMP", F_thermal);
633 
634 
635  SeriesRecorder *recorder_ptr;
636  CHKERR m_field.getInterface(recorder_ptr);
637 
638  // Read time series and do thermo-elastic analysis, this is when time
639  // dependent
640  // temperature problem was run before on the mesh. It means that before
641  // non-stationary
642  // problem was solved for temperature and filed "TEMP" is stored for
643  // subsequent time
644  // steps in the recorder.
645  if (recorder_ptr->check_series("THEMP_SERIES")) {
646  // This is time dependent case, so loop of data series stored by tape
647  // recorder.
648  // Loop over time steps
649  for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
650  sit)) {
651  PetscPrintf(PETSC_COMM_WORLD, "Process step %d\n",
652  sit->get_step_number());
653  // Load field data for this time step
654  CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
655  sit->get_step_number());
656 
657  CHKERR VecZeroEntries(F_thermal);
658  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
659  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
660 
661  // Calculate the right-hand side vector as result of thermal stresses.
662  // It MetaNodalForces
663  // that on "ELASTIC" element data structure the element implementation
664  // in thermal_stress_elem
665  // is executed.
667  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
668 
669  // Assemble vector
670  CHKERR VecAssemblyBegin(F_thermal);
671  CHKERR VecAssemblyEnd(F_thermal);
672  // Accumulate ghost dofs
673  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
674  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
675 
676  // Calculate norm of vector and print values
677  PetscReal nrm_F;
678  CHKERR VecNorm(F, NORM_2, &nrm_F);
679 
680  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
681  PetscReal nrm_F_thermal;
682  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
683  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
684  nrm_F_thermal);
685 
686  CHKERR VecScale(F_thermal, -1);
687  // check this !!!
688  CHKERR VecAXPY(F_thermal, 1, F);
689 
690  // Set dirichlet boundary to thermal stresses vector
691  dirichlet_bc_ptr->snes_x = D;
692  dirichlet_bc_ptr->snes_f = F_thermal;
693  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
694 
695  // Solve problem
696  CHKERR KSPSolve(solver, F_thermal, D);
697  // Add boundary conditions vector
698  CHKERR VecAXPY(D, 1., D0);
699  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
700  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
701 
702  // Save data on the mesh
703  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
704 
705  // Save data on mesh
706  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
707  // Post-process results
708  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
709 
710  std::ostringstream o1;
711  o1 << "out_" << sit->step_number << ".h5m";
712  CHKERR post_proc.writeFile(o1.str().c_str());
713 
714  }
715  } else {
716 
717  // This is a case when stationary problem for temperature was solved.
718  CHKERR VecZeroEntries(F_thermal);
719  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
720  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
721 
722  // Calculate the right-hand side vector with thermal stresses
724  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
725 
726  // Assemble vector
727  CHKERR VecAssemblyBegin(F_thermal);
728  CHKERR VecAssemblyEnd(F_thermal);
729 
730  // Accumulate ghost dofs
731  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
732  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
733 
734  // Calculate norm
735  PetscReal nrm_F;
736  CHKERR VecNorm(F, NORM_2, &nrm_F);
737 
738  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
739  PetscReal nrm_F_thermal;
740  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
741 
742  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n", nrm_F_thermal);
743 
744  // Add thermal stress vector and other forces vector
745  CHKERR VecScale(F_thermal, -1);
746  CHKERR VecAXPY(F_thermal, 1, F);
747 
748  // Apply kinetic boundary conditions
749  dirichlet_bc_ptr->snes_x = D;
750  dirichlet_bc_ptr->snes_f = F_thermal;
751  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
752 
753  // Solve problem
754  CHKERR KSPSolve(solver, F_thermal, D);
755  CHKERR VecAXPY(D, 1., D0);
756 
757  // Update ghost values for solution vector
758  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
759  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
760 
761  // Save data on mesh
762  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
763  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
764  // Save results to file
765  CHKERR post_proc.writeFile("out.h5m");
766 
767  }
768 
769  // Destroy vector, no needed any more
770  CHKERR VecDestroy(&F_thermal);
771 
772 
773  } else {
774  // Elastic analysis no temperature field
775  // Solve for vector D
776  CHKERR KSPSolve(solver, F, D);
777  // Add kinetic boundary conditions
778  CHKERR VecAXPY(D, 1., D0);
779  // Update ghost values
780  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
781  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
782  // Save data from vector on mesh
783  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
784  // Post-process results
785  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
786  // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
787  CHKERR post_proc.writeFile("out.h5m");
788 
789  }
790 
791  // Calculate elastic energy
792  elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
793  elastic.getLoopFeEnergy().eNergy = 0;
794  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeEnergy());
795 
796  // Print elastic energy
797  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
798  elastic.getLoopFeEnergy().eNergy);
799 
800  // Destroy matrices, vecors, solver and DM
801  CHKERR VecDestroy(&F);
802  CHKERR VecDestroy(&D);
803  CHKERR VecDestroy(&D0);
804  CHKERR MatDestroy(&Aij);
805  CHKERR KSPDestroy(&solver);
806  CHKERR DMDestroy(&dm);
807 
808  MPI_Comm_free(&moab_comm_world);
809  }
810  CATCH_ERRORS;
811 
813 
814  return 0;
815 }
#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 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:440
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
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:424
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
Interface for managing meshsets containing materials and boundary conditions.
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:490
static char help[]
Definition: elasticity.cpp:67
Projection of edge entities with one mid-node on hierarchical basis.
virtual MoFEMErrorCode build_finite_elements(int verb=-1)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
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
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=-1)=0
Add field.
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
Add boundary conditions form block set having 6 attributes.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=-1)=0
Add entities to field meshsetThe lower dimension entities are added depending on the space type...
#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
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=-1)=0
Set order approximation of the entities in the field.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:140
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:237
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:461
DEPRECATED MoFEMErrorCode set_field(const double val, const EntityType type, const std::string &field_name)
use FieldBlas
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:469
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:150
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
#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:563
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:141
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:315
virtual MoFEMErrorCode synchronise_entities(Range &ent, int verb=-1)=0
virtual MoFEMErrorCode build_fields(int verb=-1)=0
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, EntMethod &method, int lower_rank, int upper_rank, int verb=-1)=0
Make a loop over entities.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL)=0
add finite element
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 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:415
continuous field
Definition: definitions.h:175
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=-1)=0
build adjacencies
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:772
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 67 of file elasticity.cpp.