v0.7.2
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.getInterface<BitRefManager>()->setBitRefLevelByDim(
169  0, 3, bit_level0);
170 
171  // Declare approximation fields
172  CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LOBATTO_BASE, 3,
173  MB_TAG_SPARSE, MF_ZERO);
174 
175  // We can use higher oder geometry to define body
176  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
177  3, MB_TAG_SPARSE, MF_ZERO);
178 
179  // Declare problem
180 
181  // Add entities (by tets) to the field ( all entities in the mesh, root_set
182  // = 0 )
183  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
184  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
185 
186  // Set apportion order.
187  // See Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes.
188  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
189  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
190  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
191  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
192 
193  // Set order of approximation of geometry.
194  // Apply 2nd order only on skin (or in in whole body)
195  {
196  Skinner skin(&m_field.get_moab());
197  Range faces, tets;
198  CHKERR m_field.get_moab().get_entities_by_type(0, MBTET, tets);
199  // CHKERR skin.find_skin(0,tets,false,faces); CHKERRG(rval);
200  Range edges;
201  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
202  moab::Interface::UNION);
203  // CHKERR m_field.get_moab().get_adjacencies(
204  // faces,1,false,edges,moab::Interface::UNION
205  // ); CHKERRG(rval);
206  // CHKERR m_field.synchronise_entities(edges);
207  CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
208 
209  }
210  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
211 
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, 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  }
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", 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
304  elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false, true);
305 
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, "Block %d set Young modulus %3.4g\n",
313  id, elastic.setOfBlocks[id].E);
314 
315  }
316  if (block_data[id].pOisson >= -1) {
317  elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
318  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Block %d set Poisson ratio %3.4g\n",
319  id, elastic.setOfBlocks[id].PoissonRatio);
320 
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", "DISPLACEMENT");
328  CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE", "DISPLACEMENT");
329  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
330  "DISPLACEMENT");
331  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
332  "MESH_NODE_POSITIONS");
333 
335  m_field, BLOCKSET | BODYFORCESSET, it)) {
336  Range tets;
337  m_field.get_moab().get_entities_by_type(it->meshset, MBTET, tets, true);
338  CHKERR
339  m_field.add_ents_to_finite_element_by_type(tets, MBTET, "BODY_FORCE");
340  }
341 
342  // Add Neumann forces, i.e. pressure or traction forces applied on body surface. This
343  // is only declaration not implementation.
344  CHKERR MetaNeummanForces::addNeumannBCElements(m_field, "DISPLACEMENT");
345  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
346  CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
347 
348  // Add fluid pressure finite elements. This is special pressure on the surface
349  // from fluid, i.e. pressure which linearly change with the depth.
350  FluidPressure fluid_pressure_fe(m_field);
351  // This function only declare element. Element is implemented by operators in
352  // class FluidPressure.
353  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
354 
355  // Add elements for thermo elasticity if temperature field is defined.
356  ThermalStressElement thermal_stress_elem(m_field);
357  // Check if TEMP field exist, and then add element.
358  if (!m_field.check_field("TEMP")) {
359  bool add_temp_field = false;
361  if (block_data[it->getMeshsetId()].initTemp != 0) {
362  add_temp_field = true;
363  break;
364  }
365  }
366  if (add_temp_field) {
367  CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
368  MB_TAG_SPARSE, MF_ZERO);
369 
370  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
371 
372  CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
373 
374  }
375  }
376  if (m_field.check_field("TEMP")) {
377  CHKERR thermal_stress_elem.addThermalStressElement("ELASTIC",
378  "DISPLACEMENT", "TEMP");
379 
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.set_field(block_data[it->getMeshsetId()].initTemp,
401  MBVERTEX, vertices, "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, but
431  // 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 to MatrixDouble
464  // 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 DISPLACEMENT blockset.
468  // To implementations available here, depending how BC is defined on mesh file.
469  if(!flag_cubit_disp){
470  dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(new
472  m_field,"DISPLACEMENT","DISPLACEMENT",Aij,D0,F
473  )
474  );
475  } else {
476  dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(
477  new DirichletDisplacementBc(m_field,"DISPLACEMENT",Aij,D0,F)
478  );
479  }
480  // That sets Dirichlet bc objects that problem is linear, i.e. no newton (SNES) solver
481  // is run for this problem.
482  dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
483  dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
484 
485  // D0 vector will store initial displacements
486  CHKERR VecZeroEntries(D0);
487  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
488  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
489  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
490  // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
491  // implementation
492  // of particular dirichlet_bc.
493  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
494  // Set values from D0 on the field (on the mesh)
495  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
496 
497  // Calculate residual forces as result of applied kinematic constrains. Run
498  // implementation
499  // of particular finite element implementation. Look how
500  // NonlinearElasticElement is implemented,
501  // in that case. We run NonlinearElasticElement with hook material.
502  // Calculate right hand side vector
503  elastic.getLoopFeRhs().snes_f = F;
504  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeRhs());
505  // Assemble matrix
506  elastic.getLoopFeLhs().snes_B = Aij;
507  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeLhs());
508 
509  // Assemble pressure and traction forces. Run particular implemented for do
510  // this, see
511  // MetaNeummanForces how this is implemented.
512  boost::ptr_map<std::string, NeummanForcesSurface> neumann_forces;
513  CHKERR MetaNeummanForces::setMomentumFluxOperators(m_field, neumann_forces, F,
514  "DISPLACEMENT");
515 
516  {
517  boost::ptr_map<std::string, NeummanForcesSurface>::iterator mit =
518  neumann_forces.begin();
519  for (; mit != neumann_forces.end(); mit++) {
520  CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
521  &mit->second->getLoopFe());
522 
523  }
524  }
525  // Assemble forces applied to nodes, see implementation in NodalForce
526  boost::ptr_map<std::string, NodalForce> nodal_forces;
527  CHKERR
528  MetaNodalForces::setOperators(m_field, nodal_forces, F, "DISPLACEMENT");
529 
530  {
531  boost::ptr_map<std::string, NodalForce>::iterator fit =
532  nodal_forces.begin();
533  for (; fit != nodal_forces.end(); fit++) {
534  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
535  &fit->second->getLoopFe());
536  }
537  }
538  // Assemble edge forces
539  boost::ptr_map<std::string, EdgeForce> edge_forces;
540  CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F, "DISPLACEMENT");
541  {
542  boost::ptr_map<std::string, EdgeForce>::iterator fit = 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, it->getMeshsetId());
553 
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 
567  // Matrix View
568  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
569  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
570  // std::string wait;
571  // std::cin >> wait;
572 
573  // Set matrix positive defined and symmetric for Cholesky and icc
574  // pre-conditioner
575  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
576  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
577  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
578  CHKERR VecAssemblyBegin(F);
579  CHKERR VecAssemblyEnd(F);
580  CHKERR VecScale(F, -1);
581 
582  // Create solver
583  KSP solver;
584  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
585  CHKERR KSPSetDM(solver, dm);
586  CHKERR KSPSetFromOptions(solver);
587  CHKERR KSPSetOperators(solver, Aij, Aij);
588  // Setup multi-grid pre-conditioner if set from command line
589  {
590  // from PETSc example ex42.c
591  PetscBool same = PETSC_FALSE;
592  PC pc;
593  CHKERR KSPGetPC(solver, &pc);
594  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
595  if (same) {
596  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
597  CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
598  CHKERR PCSetFromOptions(pc);
599  } else {
600  // Operators are already set, do not use DM for doing that
601  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
602  }
603  }
604  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
605  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
606  // Set up solver
607  CHKERR KSPSetUp(solver);
608 
609  // Set up post-processor. It is some generic implementation of finite 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 how
629  // thermal stresses and vector is assembled in ThermalStressElement.
630  CHKERR thermal_stress_elem.setThermalStressRhsOperators("DISPLACEMENT",
631  "TEMP", F_thermal);
632 
633 
634  SeriesRecorder *recorder_ptr;
635  CHKERR m_field.getInterface(recorder_ptr);
636 
637  // Read time series and do thermo-elastic analysis, this is when time
638  // dependent
639  // temperature problem was run before on the mesh. It means that before
640  // non-stationary
641  // problem was solved for temperature and filed "TEMP" is stored for
642  // subsequent time
643  // steps in the recorder.
644  if (recorder_ptr->check_series("THEMP_SERIES")) {
645  // This is time dependent case, so loop of data series stored by tape
646  // recorder.
647  // Loop over time steps
648  for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
649  sit)) {
650  PetscPrintf(PETSC_COMM_WORLD, "Process step %d\n",
651  sit->get_step_number());
652  // Load field data for this time step
653  CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
654  sit->get_step_number());
655 
656  CHKERR VecZeroEntries(F_thermal);
657  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
658  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
659 
660  // Calculate the right-hand side vector as result of thermal stresses.
661  // It MetaNodalForces
662  // that on "ELASTIC" element data structure the element implementation
663  // in thermal_stress_elem
664  // is executed.
666  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
667 
668  // Assemble vector
669  CHKERR VecAssemblyBegin(F_thermal);
670  CHKERR VecAssemblyEnd(F_thermal);
671  // Accumulate ghost dofs
672  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
673  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
674 
675  // Calculate norm of vector and print values
676  PetscReal nrm_F;
677  CHKERR VecNorm(F, NORM_2, &nrm_F);
678 
679  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
680  PetscReal nrm_F_thermal;
681  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
682  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
683  nrm_F_thermal);
684 
685  CHKERR VecScale(F_thermal, -1);
686  // check this !!!
687  CHKERR VecAXPY(F_thermal, 1, F);
688 
689  // Set dirichlet boundary to thermal stresses vector
690  dirichlet_bc_ptr->snes_x = D;
691  dirichlet_bc_ptr->snes_f = F_thermal;
692  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
693 
694  // Solve problem
695  CHKERR KSPSolve(solver, F_thermal, D);
696  // Add boundary conditions vector
697  CHKERR VecAXPY(D, 1., D0);
698  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
699  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
700 
701  // Save data on the mesh
702  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
703 
704  // Save data on mesh
705  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
706  // Post-process results
707  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
708 
709  std::ostringstream o1;
710  o1 << "out_" << sit->step_number << ".h5m";
711  CHKERR post_proc.writeFile(o1.str().c_str());
712 
713  }
714  } else {
715 
716  // This is a case when stationary problem for temperature was solved.
717  CHKERR VecZeroEntries(F_thermal);
718  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
719  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
720 
721  // Calculate the right-hand side vector with thermal stresses
723  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
724 
725  // Assemble vector
726  CHKERR VecAssemblyBegin(F_thermal);
727  CHKERR VecAssemblyEnd(F_thermal);
728 
729  // Accumulate ghost dofs
730  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
731  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
732 
733  // Calculate norm
734  PetscReal nrm_F;
735  CHKERR VecNorm(F, NORM_2, &nrm_F);
736 
737  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
738  PetscReal nrm_F_thermal;
739  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
740 
741  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n", nrm_F_thermal);
742 
743  // Add thermal stress vector and other forces vector
744  CHKERR VecScale(F_thermal, -1);
745  CHKERR VecAXPY(F_thermal, 1, F);
746 
747  // Apply kinetic boundary conditions
748  dirichlet_bc_ptr->snes_x = D;
749  dirichlet_bc_ptr->snes_f = F_thermal;
750  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
751 
752  // Solve problem
753  CHKERR KSPSolve(solver, F_thermal, D);
754  CHKERR VecAXPY(D, 1., D0);
755 
756  // Update ghost values for solution vector
757  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
758  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
759 
760  // Save data on mesh
761  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
762  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
763  // Save results to file
764  CHKERR post_proc.writeFile("out.h5m");
765 
766  }
767 
768  // Destroy vector, no needed any more
769  CHKERR VecDestroy(&F_thermal);
770 
771 
772  } else {
773  // Elastic analysis no temperature field
774  // Solve for vector D
775  CHKERR KSPSolve(solver, F, D);
776  // Add kinetic boundary conditions
777  CHKERR VecAXPY(D, 1., D0);
778  // Update ghost values
779  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
780  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
781  // Save data from vector on mesh
782  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
783  // Post-process results
784  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
785  // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
786  CHKERR post_proc.writeFile("out.h5m");
787 
788  }
789 
790  // Calculate elastic energy
791  elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
792  elastic.getLoopFeEnergy().eNergy = 0;
793  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &elastic.getLoopFeEnergy());
794 
795  // Print elastic energy
796  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
797  elastic.getLoopFeEnergy().eNergy);
798 
799  // Destroy matrices, vecors, solver and DM
800  CHKERR VecDestroy(&F);
801  CHKERR VecDestroy(&D);
802  CHKERR VecDestroy(&D0);
803  CHKERR MatDestroy(&Aij);
804  CHKERR KSPDestroy(&solver);
805  CHKERR DMDestroy(&dm);
806 
807  MPI_Comm_free(&moab_comm_world);
808  }
809  CATCH_ERRORS;
810 
812 
813  return 0;
814 }
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 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:469
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:449
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:558
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:531
static char help[]
Definition: elasticity.cpp:67
char mesh_file_name[255]
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 meshset.
#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.
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:141
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:238
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:492
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:502
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:610
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:142
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:316
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:462
continuous field
Definition: definitions.h:176
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:851
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.