v0.5.86
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 86 of file elasticity.cpp.

86  {
87 
88  // Initialize PETSCc
89  PetscInitialize(&argc,&argv,(char *)0,help);
90 
91  try {
92 
93  PetscBool flg_block_config,flg_file;
94  char mesh_file_name[255];
95  char block_config_file[255];
96  PetscInt order = 2;
97  PetscBool is_partitioned = PETSC_FALSE;
98 
99  // Read options from command line
100  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Elastic Config","none"); CHKERRQ(ierr);
101  ierr = PetscOptionsString(
102  "-my_file",
103  "mesh file name","",
104  "mesh.h5m",mesh_file_name,255,&flg_file
105  ); CHKERRQ(ierr);
106  ierr = PetscOptionsInt(
107  "-my_order",
108  "default approximation order","",
109  2,&order,PETSC_NULL
110  ); CHKERRQ(ierr);
111  ierr = PetscOptionsBool(
112  "-my_is_partitioned",
113  "set if mesh is partitioned (this result that each process keeps only part of the mes","",
114  PETSC_FALSE,&is_partitioned,PETSC_NULL
115  ); CHKERRQ(ierr);
116  ierr = PetscOptionsString("-my_block_config",
117  "elastic configure file name","",
118  "block_conf.in",block_config_file,255,&flg_block_config
119  ); CHKERRQ(ierr);
120  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
121 
122  // Throw error if file with mesh is not give
123  if(flg_file != PETSC_TRUE) {
124  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
125  }
126 
127  // Create mesh databse
128  moab::Core mb_instance;
129  moab::Interface& moab = mb_instance;
130 
131  // Create moab communicator
132  // Create separate MOAB communicator, it is duplicate of PETSc communicator.
133  // NOTE That this should eliminate potential communication problems between
134  // MOAB and PETSC functions.
135  MPI_Comm moab_comm_world;
136  MPI_Comm_dup(PETSC_COMM_WORLD,&moab_comm_world);
137  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
138  if(pcomm == NULL) pcomm = new ParallelComm(&moab,moab_comm_world);
139 
140  // Read whole mesh or part of is if partitioned
141  if(is_partitioned == PETSC_TRUE) {
142  // This is a case of distributed mesh and algebra. In that case each processor
143  // keep only part of the problem.
144  const char *option;
145  option = "PARALLEL=READ_PART;"
146  "PARALLEL_RESOLVE_SHARED_ENTS;"
147  "PARTITION=PARALLEL_PARTITION;";
148  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
149  } else {
150  // If that case we have distributed algebra, i.e. assembly of vectors and
151  // matrices is in parallel, but whole mesh is stored on all processors.
152  // Solver and matrix scales well, however problem set-up of problem is
153  // not fully parallel.
154  const char *option;
155  option = "";
156  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
157  }
158 
159  // Create MoFEM databse and link it to MoAB
160  MoFEM::Core core(moab);
161  MoFEM::Interface& m_field = core;
162 
163  // Print boundary conditions and material parameters
164  MeshsetsManager *mmanager_ptr;
165  ierr = m_field.query_interface(mmanager_ptr); CHKERRQ(ierr);
166  ierr = mmanager_ptr->printDisplacementSet(); CHKERRQ(ierr);
167  ierr = mmanager_ptr->printForceSet(); CHKERRQ(ierr);
168  ierr = mmanager_ptr->printMaterialsSet(); CHKERRQ(ierr);
169 
170  // Set bit refinement level to all entities (we do not refine mesh in this example
171  // so all entities have the same bit refinement level)
172  BitRefLevel bit_level0;
173  bit_level0.set(0);
174  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
175 
176  // Declare approximation fields
177  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LOBBATO_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
178  // We can use higher oder geometry to define body
179  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
180 
181  // Declare problem
182 
183  // Add entities (by tets) to the field ( all entities in the mesh, root_set = 0 )
184  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
185  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
186 
187  // Set apportion order.
188  // See Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes.
189  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
190  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
191  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
192  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
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  rval = m_field.get_moab().get_entities_by_type(0,MBTET,tets); CHKERRQ_MOAB(rval);
200  // rval = skin.find_skin(0,tets,false,faces); CHKERRQ_MOAB(rval);
201  Range edges;
202  rval = m_field.get_moab().get_adjacencies(
203  tets,1,false,edges,moab::Interface::UNION
204  ); CHKERRQ_MOAB(rval);
205  // rval = m_field.get_moab().get_adjacencies(
206  // faces,1,false,edges,moab::Interface::UNION
207  // ); CHKERRQ_MOAB(rval);
208  // ierr = m_field.synchronise_entities(edges); CHKERRQ(ierr);
209  ierr = m_field.set_field_order(edges,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
210  }
211  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
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(),po::value<int>(&block_data[it->getMeshsetId()].oRder)->default_value(order));
227  std::ostringstream str_cond;
228  str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
229  config_file_options.add_options()
230  (str_cond.str().c_str(),po::value<double>(&block_data[it->getMeshsetId()].yOung)->default_value(-1));
231  std::ostringstream str_capa;
232  str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
233  config_file_options.add_options()
234  (str_capa.str().c_str(),po::value<double>(&block_data[it->getMeshsetId()].pOisson)->default_value(-2));
235  std::ostringstream str_init_temp;
236  str_init_temp << "block_" << it->getMeshsetId() << ".initial_temperature";
237  config_file_options.add_options()
238  (str_init_temp.str().c_str(),po::value<double>(&block_data[it->getMeshsetId()].initTemp)->default_value(0));
239  }
240  po::parsed_options parsed = parse_config_file(ini_file,config_file_options,true);
241  store(parsed,vm);
242  po::notify(vm);
244  if(block_data[it->getMeshsetId()].oRder == -1) continue;
245  if(block_data[it->getMeshsetId()].oRder == order) continue;
246  PetscPrintf(PETSC_COMM_WORLD,"Set block %d order to %d\n",it->getMeshsetId(),block_data[it->getMeshsetId()].oRder);
247  Range block_ents;
248  rval = moab.get_entities_by_handle(it->getMeshset(),block_ents,true); CHKERRQ_MOAB(rval);
249  Range ents_to_set_order;
250  rval = moab.get_adjacencies(block_ents,3,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
251  ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
252  rval = moab.get_adjacencies(block_ents,2,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
253  rval = moab.get_adjacencies(block_ents,1,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
254  ierr = m_field.synchronise_entities(ents_to_set_order); CHKERRQ(ierr);
255  ierr = m_field.set_field_order(ents_to_set_order,"DISPLACEMENT",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
256  }
257  std::vector<std::string> additional_parameters;
258  additional_parameters = collect_unrecognized(parsed.options,po::include_positional);
259  for(std::vector<std::string>::iterator vit = additional_parameters.begin();
260  vit!=additional_parameters.end();vit++) {
261  ierr = PetscPrintf(PETSC_COMM_WORLD,"** WARNING Unrecognized option %s\n",vit->c_str()); CHKERRQ(ierr);
262  }
263  } catch (const std::exception& ex) {
264  std::ostringstream ss;
265  ss << ex.what() << std::endl;
266  SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
267  }
268  }
269 
270  // Add elastic element
271  boost::shared_ptr<Hooke<adouble> > hooke_adouble_ptr(new Hooke<adouble>());
272  boost::shared_ptr<Hooke<double> > hooke_double_ptr(new Hooke<double>());
273  NonlinearElasticElement elastic(m_field,2);
274  ierr = elastic.setBlocks(hooke_double_ptr,hooke_adouble_ptr); CHKERRQ(ierr);
275  ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
276  ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
277 
278  // Update material parameters. Set material parameters block by block.
280  int id = it->getMeshsetId();
281  if(block_data[id].yOung>0) {
282  elastic.setOfBlocks[id].E = block_data[id].yOung;
283  ierr = PetscPrintf(
284  PETSC_COMM_WORLD,"Block %d set Young modulus %3.4g\n",id,elastic.setOfBlocks[id].E
285  ); CHKERRQ(ierr);
286  }
287  if(block_data[id].pOisson>=-1) {
288  elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
289  ierr = PetscPrintf(
290  PETSC_COMM_WORLD,"Block %d set Poisson ratio %3.4g\n",id,elastic.setOfBlocks[id].PoissonRatio
291  ); CHKERRQ(ierr);
292  }
293  }
294 
295  // Add body force element. This is only declaration of element. not its implementation.
296  ierr = m_field.add_finite_element("BODY_FORCE"); CHKERRQ(ierr);
297  ierr = m_field.modify_finite_element_add_field_row("BODY_FORCE","DISPLACEMENT"); CHKERRQ(ierr);
298  ierr = m_field.modify_finite_element_add_field_col("BODY_FORCE","DISPLACEMENT"); CHKERRQ(ierr);
299  ierr = m_field.modify_finite_element_add_field_data("BODY_FORCE","DISPLACEMENT"); CHKERRQ(ierr);
300  ierr = m_field.modify_finite_element_add_field_data("BODY_FORCE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
302  Range tets;
303  rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTET,tets,true); CHKERRQ_MOAB(rval);
304  ierr = m_field.add_ents_to_finite_element_by_type(tets,MBTET,"BODY_FORCE"); CHKERRQ(ierr);
305  }
306 
307  // Add Neumann forces, i.e. pressure or traction forces applied on body surface. This
308  // is only declaration not implementation.
309  ierr = MetaNeummanForces::addNeumannBCElements(m_field,"DISPLACEMENT"); CHKERRQ(ierr);
310  ierr = MetaNodalForces::addElement(m_field,"DISPLACEMENT"); CHKERRQ(ierr);
311  ierr = MetaEdgeForces::addElement(m_field,"DISPLACEMENT"); CHKERRQ(ierr);
312 
313  // Add fluid pressure finite elements. This is special pressure on the surface from
314  // fluid, i.e. preseeure which linearly change with the depth.
315  FluidPressure fluid_pressure_fe(m_field);
316  // This function only declare element. Element is implemented by operators in
317  // class FluidPressure.
318  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
319 
320  // Add elements for thermo elasticity if temperature field is defined.
321  ThermalStressElement thermal_stress_elem(m_field);
322  // Check if TEMP field exist, and then add element.
323  if(!m_field.check_field("TEMP")) {
324  bool add_temp_field = false;
326  if(block_data[it->getMeshsetId()].initTemp!=0) {
327  add_temp_field = true;
328  break;
329  }
330  }
331  if(add_temp_field) {
332  ierr = m_field.add_field("TEMP",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
333  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"TEMP"); CHKERRQ(ierr);
334  ierr = m_field.set_field_order(0,MBVERTEX,"TEMP",1); CHKERRQ(ierr);
335  }
336  }
337  if(m_field.check_field("TEMP")) {
338  ierr = thermal_stress_elem.addThermalSterssElement("ELASTIC","DISPLACEMENT","TEMP"); CHKERRQ(ierr);
339  }
340 
341  // All is declared, at that point build fields first,
342  ierr = m_field.build_fields(); CHKERRQ(ierr);
343  // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class Projection10NodeCoordsOnField
344  // do the trick.
345  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
346  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
347  if(m_field.check_field("TEMP")) {
349  if(block_data[it->getMeshsetId()].initTemp!=0) {
350  PetscPrintf(PETSC_COMM_WORLD,"Set block %d temperature to %3.2g\n",
351  it->getMeshsetId(),block_data[it->getMeshsetId()].initTemp);
352  Range block_ents;
353  rval = moab.get_entities_by_handle(it->meshset,block_ents,true); CHKERRQ_MOAB(rval);
354  Range vertices;
355  rval = moab.get_connectivity(block_ents,vertices,true); CHKERRQ_MOAB(rval);
356  ierr = m_field.set_field(block_data[it->getMeshsetId()].initTemp,MBVERTEX,vertices,"TEMP"); CHKERRQ(ierr);
357  }
358  }
359  }
360 
361  // Build database for elements. Actual implementation of element is not need here,
362  // only elements has to be declared.
363  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
364  // Build adjacencies between elements and field entities
365  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
366 
367  // Register MOFEM DM implementation in PETSc
369 
370  // Create DM manager
371  DM dm;
372  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
373  ierr = DMSetType(dm,"MOFEM");CHKERRQ(ierr);
374  ierr = DMMoFEMCreateMoFEM(dm,&m_field,"ELASTIC_PROB",bit_level0); CHKERRQ(ierr);
375  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
376  ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
377  // Add elements to DM manager
378  ierr = DMMoFEMAddElement(dm,"ELASTIC"); CHKERRQ(ierr);
379  ierr = DMMoFEMAddElement(dm,"BODY_FORCE"); CHKERRQ(ierr);
380  ierr = DMMoFEMAddElement(dm,"FLUID_PRESSURE_FE"); CHKERRQ(ierr);
381  ierr = DMMoFEMAddElement(dm,"FORCE_FE"); CHKERRQ(ierr);
382  ierr = DMMoFEMAddElement(dm,"PRESSURE_FE"); CHKERRQ(ierr);
383  ierr = DMSetUp(dm); CHKERRQ(ierr);
384 
385  // Create matrices & vectors. Note that native PETSc DM interface is used, but
386  // ubder the PETSc interface MOFEM implementation is running.
387  Vec F,D,D0;
388  ierr = DMCreateGlobalVector(dm,&F); CHKERRQ(ierr);
389  ierr = VecDuplicate(F,&D); CHKERRQ(ierr);
390  ierr = VecDuplicate(F,&D0); CHKERRQ(ierr);
391  Mat Aij;
392  ierr = DMCreateMatrix(dm,&Aij); CHKERRQ(ierr);
393  ierr = MatSetOption(Aij,MAT_SPD,PETSC_TRUE); CHKERRQ(ierr);
394 
395  // Zero vectors and matrices
396  ierr = VecZeroEntries(F); CHKERRQ(ierr);
397  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
398  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
399  ierr = VecZeroEntries(D); CHKERRQ(ierr);
400  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
401  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
402  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
403  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
404 
405  // This controls how kinematic constrains are set, by blockset or nodeset. Cubit
406  // sets kinetic booundary conditions by blockset.
407  bool flag_cubit_disp = false;
409  flag_cubit_disp = true;
410  }
411 
412  // Below particular implementations of finite elements are used to assemble
413  // problem matrixes and vectors. Implementation of element does not change
414  // how element is declared.
415 
416  // Assemble Aij and F. Define dirihlet_bc element, which stets constrains to MatrixDouble
417  // and the right hand side vector.
418  boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
419 
420  // if normally defined boundary conditions are not found, try to use DISPLACEMENT blockset.
421  // To implementations available here, depending how BC is defined on mesh file.
422  if(!flag_cubit_disp){
423  dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(new
425  m_field,"DISPLACEMENT","DISPLACEMENT",Aij,D0,F
426  )
427  );
428  } else {
429  dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
430  new DirichletDisplacementBc(m_field,"DISPLACEMENT",Aij,D0,F)
431  );
432  }
433  // That sets dirihlet_bc objects that problem is linear, i.e. no newton (SNES) solver
434  // is run for this problem.
435  dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
436  dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
437 
438  // D0 vector will store initial displacements
439  ierr = VecZeroEntries(D0); CHKERRQ(ierr);
440  ierr = VecGhostUpdateBegin(D0,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
441  ierr = VecGhostUpdateEnd(D0,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
442  ierr = DMoFEMMeshToLocalVector(dm,D0,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
443  // Run dirihlet_bc, from that on the mesh set values in vector D0. Run implementation
444  // of particular dirihlet_bc.
445  ierr = DMoFEMPreProcessFiniteElements(dm,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
446  // Set values from D0 on the field (on the mesh)
447  ierr = DMoFEMMeshToLocalVector(dm,D0,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
448 
449  // Calculate residual forces as result of applied kinematic constrains. Run implementation
450  // of particular finite element implementation. Look how NonlinearElasticElement is implemented,
451  // in that case. We run NonlinearElasticElement with hook material.
452  // Calculate right hand side vector
453  elastic.getLoopFeRhs().snes_f = F;
454  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&elastic.getLoopFeRhs()); CHKERRQ(ierr);
455  // Assemble matrix
456  elastic.getLoopFeLhs().snes_B = Aij;
457  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&elastic.getLoopFeLhs()); CHKERRQ(ierr);
458 
459  // Assemble pressure and traction forces. Run particular implemented for do this, see
460  // MetaNeummanForces how this is implemented.
461  boost::ptr_map<std::string,NeummanForcesSurface> neumann_forces;
462  ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,F,"DISPLACEMENT"); CHKERRQ(ierr);
463  {
464  boost::ptr_map<std::string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
465  for(;mit!=neumann_forces.end();mit++) {
466  ierr = DMoFEMLoopFiniteElements(dm,mit->first.c_str(),&mit->second->getLoopFe()); CHKERRQ(ierr);
467  }
468  }
469  // Assemble forces applied to nodes, see implementation in NodalForce
470  boost::ptr_map<std::string,NodalForce> nodal_forces;
471  ierr = MetaNodalForces::setOperators(m_field,nodal_forces,F,"DISPLACEMENT"); CHKERRQ(ierr);
472  {
473  boost::ptr_map<std::string,NodalForce>::iterator fit = nodal_forces.begin();
474  for(;fit!=nodal_forces.end();fit++) {
475  ierr = DMoFEMLoopFiniteElements(dm,fit->first.c_str(),&fit->second->getLoopFe()); CHKERRQ(ierr);
476  }
477  }
478  // Assemble edge forces
479  boost::ptr_map<std::string,EdgeForce> edge_forces;
480  ierr = MetaEdgeForces::setOperators(m_field,edge_forces,F,"DISPLACEMENT"); CHKERRQ(ierr);
481  {
482  boost::ptr_map<std::string,EdgeForce>::iterator fit = edge_forces.begin();
483  for(;fit!=edge_forces.end();fit++) {
484  ierr = DMoFEMLoopFiniteElements(dm,fit->first.c_str(),&fit->second->getLoopFe()); CHKERRQ(ierr);
485  }
486  }
487  // Assemble body forces, implemented in BodyFroceConstantField
488  BodyFroceConstantField body_forces_methods(m_field);
490  ierr = body_forces_methods.addBlock("DISPLACEMENT",F,it->getMeshsetId()); CHKERRQ(ierr);
491  }
492  ierr = DMoFEMLoopFiniteElements(dm,"BODY_FORCE",&body_forces_methods.getLoopFe()); CHKERRQ(ierr);
493  // Assemble fluid pressure forces
494  ierr = fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators("DISPLACEMENT",F,false,true); CHKERRQ(ierr);
495  ierr = DMoFEMLoopFiniteElements(dm,"FLUID_PRESSURE_FE",&fluid_pressure_fe.getLoopFe()); CHKERRQ(ierr);
496 
497  // Apply kinematic constrain to right hand side vector and matrix
498  ierr = DMoFEMPostProcessFiniteElements(dm,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
499 
500  //Matrix View
501  //MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
502  //MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
503  //std::string wait;
504  //std::cin >> wait;
505 
506  // Set matrix positive defined and symmetric for Cholesky and icc pre-conditioner
507  ierr = MatSetOption(Aij,MAT_SPD,PETSC_TRUE); CHKERRQ(ierr);
508  ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
509  ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
510  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
511  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
512  ierr = VecScale(F,-1); CHKERRQ(ierr);
513 
514  // Create solver
515  KSP solver;
516  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
517  ierr = KSPSetDM(solver,dm); CHKERRQ(ierr);
518  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
519  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
520  // Setup multi-grid pre-conditioner if set from command line
521  {
522  //from PETSc example ex42.c
523  PetscBool same = PETSC_FALSE;
524  PC pc;
525  ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
526  PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
527  if (same) {
528  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm,Aij,true);
529  ierr = PCMGSetUpViaApproxOrders(pc,&pc_ctx); CHKERRQ(ierr);
530  ierr = PCSetFromOptions(pc); CHKERRQ(ierr);
531  } else {
532  // Operators are already set, do not use DM for doing that
533  ierr = KSPSetDMActive(solver,PETSC_FALSE); CHKERRQ(ierr);
534  }
535  }
536  ierr = KSPSetInitialGuessKnoll(solver,PETSC_FALSE); CHKERRQ(ierr);
537  ierr = KSPSetInitialGuessNonzero(solver,PETSC_TRUE); CHKERRQ(ierr);
538  // Set up solver
539  ierr = KSPSetUp(solver); CHKERRQ(ierr);
540 
541  // Set up post-procesor. It is some generic implementation of finite element.
542  PostProcVolumeOnRefinedMesh post_proc(m_field);
543  // Add operators to the elements, strating with some generic
544  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
545  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
546  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
547  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
548  // Add problem specific operator on element to post-process stresses
549  post_proc.getOpPtrVector().push_back(
550  new PostPorcHookStress(
551  m_field,
552  post_proc.postProcMesh,
553  post_proc.mapGaussPts,
554  "DISPLACEMENT",
555  post_proc.commonData,
556  &elastic.setOfBlocks
557  )
558  );
559 
560  // Temperature field is defined on the mesh
561  if(m_field.check_field("TEMP")) {
562 
563  // Create thermal vector
564  Vec F_thermal;
565  ierr = VecDuplicate(F,&F_thermal); CHKERRQ(ierr);
566 
567  // Set up implementation for calculation of thermal stress vector. Look how
568  // thermal stresses and vector is assembled in ThermalStressElement.
569  ierr = thermal_stress_elem.setThermalStressRhsOperators("DISPLACEMENT","TEMP",F_thermal); CHKERRQ(ierr);
570 
571  SeriesRecorder *recorder_ptr;
572  ierr = m_field.query_interface(recorder_ptr); CHKERRQ(ierr);
573  // Read time series and do thermo-elastic analysis, this is when time dependent
574  // temperature problem was run before on the mesh. It means that before non-stationary
575  // problem was solved for temperature and filed "TEMP" is stored for subsequent time
576  // steps in the recorder.
577  if( recorder_ptr->check_series("THEMP_SERIES") ) {
578  // This is time dependent case, so loop of data series stored by tape recorder.
579  // Loop over time steps
580  for(_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr,"THEMP_SERIES",sit)) {
581  PetscPrintf(PETSC_COMM_WORLD,"Process step %d\n",sit->get_step_number());
582  // Load field data for this time step
583  ierr = recorder_ptr->load_series_data("THEMP_SERIES",sit->get_step_number()); CHKERRQ(ierr);
584  ierr = VecZeroEntries(F_thermal); CHKERRQ(ierr);
585  ierr = VecGhostUpdateBegin(F_thermal,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
586  ierr = VecGhostUpdateEnd(F_thermal,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
587  // Calculate the right-hand side vector as result of thermal stresses. It MetaNodalForces
588  // that on "ELASTIC" element data structure the element implementation in thermal_stress_elem
589  // is executed.
590  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&thermal_stress_elem.getLoopThermalStressRhs()); CHKERRQ(ierr);
591  // Assemble vector
592  ierr = VecAssemblyBegin(F_thermal); CHKERRQ(ierr);
593  ierr = VecAssemblyEnd(F_thermal); CHKERRQ(ierr);
594  // Accumulate ghost dofs
595  ierr = VecGhostUpdateBegin(F_thermal,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
596  ierr = VecGhostUpdateEnd(F_thermal,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
597 
598  // Calculate norm of vector and print values
599  PetscReal nrm_F;
600  ierr = VecNorm(F,NORM_2,&nrm_F); CHKERRQ(ierr);
601  PetscPrintf(PETSC_COMM_WORLD,"norm2 F = %6.4e\n",nrm_F);
602  PetscReal nrm_F_thremal;
603  ierr = VecNorm(F_thermal,NORM_2,&nrm_F_thremal); CHKERRQ(ierr);
604  PetscPrintf(PETSC_COMM_WORLD,"norm2 F_thernal = %6.4e\n",nrm_F_thremal);
605 
606  ierr = VecScale(F_thermal,-1); CHKERRQ(ierr); //check this !!!
607  ierr = VecAXPY(F_thermal,1,F); CHKERRQ(ierr);
608 
609  // Set dirihlet boundary to thermal stresses vector
610  dirihlet_bc_ptr->snes_x = D;
611  dirihlet_bc_ptr->snes_f = F_thermal;
612  ierr = DMoFEMPostProcessFiniteElements(dm,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
613 
614  // Solve problem
615  ierr = KSPSolve(solver,F_thermal,D); CHKERRQ(ierr);
616  // Add boundary conditions vector
617  ierr = VecAXPY(D,1.,D0); CHKERRQ(ierr);
618  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
619  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
620 
621  // Sava data on the mesh
622  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
623 
624  // Save data on mesh
625  ierr = DMoFEMPreProcessFiniteElements(dm,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
626  // Post-process results
627  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&post_proc); CHKERRQ(ierr);
628  std::ostringstream o1;
629  o1 << "out_" << sit->step_number << ".h5m";
630  ierr = post_proc.writeFile(o1.str().c_str()); CHKERRQ(ierr);
631  }
632  } else {
633 
634  // This is a case when stationary problem for temperature was solved.
635  ierr = VecZeroEntries(F_thermal); CHKERRQ(ierr);
636  ierr = VecGhostUpdateBegin(F_thermal,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
637  ierr = VecGhostUpdateEnd(F_thermal,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
638  // Calculate the right-hand side vector with thermal stresses
639  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&thermal_stress_elem.getLoopThermalStressRhs()); CHKERRQ(ierr);
640  // Assemble vector
641  ierr = VecAssemblyBegin(F_thermal); CHKERRQ(ierr);
642  ierr = VecAssemblyEnd(F_thermal); CHKERRQ(ierr);
643  // Accumulate ghost dofs
644  ierr = VecGhostUpdateBegin(F_thermal,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
645  ierr = VecGhostUpdateEnd(F_thermal,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
646 
647  // Calculate norm
648  PetscReal nrm_F;
649  ierr = VecNorm(F,NORM_2,&nrm_F); CHKERRQ(ierr);
650  PetscPrintf(PETSC_COMM_WORLD,"norm2 F = %6.4e\n",nrm_F);
651  PetscReal nrm_F_thremal;
652  ierr = VecNorm(F_thermal,NORM_2,&nrm_F_thremal); CHKERRQ(ierr);
653  PetscPrintf(PETSC_COMM_WORLD,"norm2 F_thernal = %6.4e\n",nrm_F_thremal);
654 
655  // Add thermal stress vector and other forces vector
656  ierr = VecScale(F_thermal,-1); CHKERRQ(ierr);
657  ierr = VecAXPY(F_thermal,1,F); CHKERRQ(ierr);
658 
659  // Apply kinetic boundary conditions
660  dirihlet_bc_ptr->snes_x = D;
661  dirihlet_bc_ptr->snes_f = F_thermal;
662  ierr = DMoFEMPostProcessFiniteElements(dm,dirihlet_bc_ptr.get()); CHKERRQ(ierr);
663 
664  // Solve problem
665  ierr = KSPSolve(solver,F_thermal,D); CHKERRQ(ierr);
666  ierr = VecAXPY(D,1.,D0); CHKERRQ(ierr);
667  // Update ghost values for solution vector
668  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
669  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
670 
671  // Save data on mesh
672  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
673  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&post_proc); CHKERRQ(ierr);
674  // Save results to file
675  ierr = post_proc.writeFile("out.h5m"); CHKERRQ(ierr);
676 
677  }
678 
679  // Destroy vector, no needed any more
680  ierr = VecDestroy(&F_thermal); CHKERRQ(ierr);
681 
682  } else {
683  // Elastic analysis no temperature field
684  // Solve for vector D
685  ierr = KSPSolve(solver,F,D); CHKERRQ(ierr);
686  // Add kinetic boundary conditions
687  ierr = VecAXPY(D,1.,D0); CHKERRQ(ierr);
688  // Update ghost values
689  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
690  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
691  // Save data from vector on mesh
692  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
693  // Porticoes results by running post_proc implementation of element
694  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&post_proc); CHKERRQ(ierr);
695  // Weite mesh in parallel (uisng h5m MOAB format, writing is in parallel)
696  ierr = post_proc.writeFile("out.h5m"); CHKERRQ(ierr);
697  }
698 
699  // Calculate elastic energy
700  elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
701  elastic.getLoopFeEnergy().eNergy = 0;
702  ierr = DMoFEMLoopFiniteElements(dm,"ELASTIC",&elastic.getLoopFeEnergy()); CHKERRQ(ierr);
703  // Print elastic energy
704  PetscPrintf(PETSC_COMM_WORLD,"Elastic energy %6.4e\n",elastic.getLoopFeEnergy().eNergy);
705 
706  // Destroy matrices, vecors, solver and DM
707  ierr = VecDestroy(&F); CHKERRQ(ierr);
708  ierr = VecDestroy(&D); CHKERRQ(ierr);
709  ierr = VecDestroy(&D0); CHKERRQ(ierr);
710  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
711  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
712  ierr = DMDestroy(&dm); CHKERRQ(ierr);
713 
714  MPI_Comm_free(&moab_comm_world);
715 
716  } catch (MoFEMException const &e) {
717  SETERRQ(PETSC_COMM_SELF,e.errorCode,e.errorMessage);
718  }
719 
720  PetscFinalize();
721 
722  return 0;
723 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:362
static PetscErrorCode 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.
#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...
static PetscErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:104
MoFEMErrorCodes errorCode
Definition: Common.hpp:105
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
Fluid pressure forces.
Implentation of thermal stress element.
virtual PetscErrorCode 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...
virtual PetscErrorCode build_adjacencies(const Range &ents, int verb=-1)=0
build adjacencies
Hook equation.
Definition: Hooke.hpp:31
static MoABErrorCode rval
Definition: Common.hpp:25
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:428
static PetscErrorCode 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
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:444
PetscErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
virtual PetscErrorCode 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.
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_maks=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:154
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:494
virtual moab::Interface & get_moab()=0
virtual PetscErrorCode 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.
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:473
virtual PetscErrorCode synchronise_entities(Range &ent, int verb=-1)=0
Interface for managing meshsets containing materials and boundary conditions.
virtual PetscErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual PetscErrorCode load_series_data(const std::string &serie_name, const int step_number)
Core (interface) classThis is the implementation of abstract MoFEM::Interface class. Similarly to the convention used in MoAB, we use small letters to name function of purely abstract classes. This is an exception used only here. For more details about naming functions see Coding practice.
Definition: Core.hpp:44
Body forces elements.
Definition: BodyForce.hpp:25
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
static char help[]
Definition: elasticity.cpp:67
Projection of edge entities with one mid-node on hierarchical basis.
char errorMessage[255]
Definition: Common.hpp:106
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:48
PetscErrorCode query_interface(IFace *&ptr) const
Definition: Interface.hpp:47
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:465
virtual bool check_series(const std::string &name) const
check if series is in database
CHKERRQ(ierr)
Add boundary conditions form block set having 6 attributes.
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:43
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:125
static PetscErrorCode 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.
Post processing.
block name is "BODY_FORCES"
Definition: definitions.h:225
Set data structures of MG pre-conditioner via approximation orders.
virtual PetscErrorCode build_finite_elements(int verb=-1)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
PetscErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
virtual PetscErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_cooficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=-1)=0
Add field.
virtual DEPRECATED PetscErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)=0
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
PetscErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:305
virtual PetscErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL)=0
add finite element
virtual PetscErrorCode build_fields(int verb=-1)=0
virtual PetscErrorCode 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
Like AINSWORTH_LEGENDRE_BASE but with Lobatto base instead Legendre .
Definition: definitions.h:126
virtual bool check_field(const std::string &name) const =0
check if field is in database
InterfaceThis interface is used by user to:
Definition: Interface.hpp:38
virtual DEPRECATED PetscErrorCode set_field(const double val, const EntityType type, const std::string &field_name)=0
use FieldBlas
continuous field
Definition: definitions.h:159
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:776
virtual PetscErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual PetscErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
static PetscErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:79
static PetscErrorCode 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
PetscErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Operator post-procesing stresses for Hook&#39;e isotropic material.
PetscErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
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.