v0.9.0
Classes | Typedefs | Functions | Variables
#include <BasicFiniteElements.hpp>
#include <boost/program_options.hpp>
#include <Hooke.hpp>

Go to the source code of this file.

Classes

struct  BlockOptionData
 stores some optional data More...
 
struct  VolRule
 Set integration rule. More...
 
struct  PrismFE
 

Typedefs

using BlockData = NonlinearElasticElement::BlockData
 
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
 

Functions

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

Variables

static char help []
 

Typedef Documentation

◆ BlockData

◆ VolUserDataOperator

Definition at line 90 of file elasticity.cpp.

Function Documentation

◆ main()

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

Definition at line 108 of file elasticity.cpp.

108  {
109 
110  const string default_options = "-ksp_type gmres \n"
111  "-pc_type lu \n"
112  "-pc_factor_mat_solver_package mumps \n"
113  "-ksp_monitor \n"
114  "-snes_type newtonls \n"
115  "-snes_linesearch_type basic \n"
116  "-snes_atol 1e-8 \n"
117  "-snes_rtol 1e-8 \n"
118  "-snes_monitor \n"
119  "-ts_monitor \n"
120  "-ts_type beuler \n";
121 
122  string param_file = "param_file.petsc";
123  if (!static_cast<bool>(ifstream(param_file))) {
124  std::ofstream file(param_file.c_str(), std::ios::ate);
125  if (file.is_open()) {
126  file << default_options;
127  file.close();
128  }
129  }
130 
131  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
132 
133  try {
134 
135  PetscBool flg_block_config, flg_file;
136  char mesh_file_name[255];
137  char block_config_file[255];
138  PetscInt test_nb = 0;
139  PetscInt order = 2;
140  PetscBool is_partitioned = PETSC_FALSE;
141 
142  // Select base
143  enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
144  const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier"};
145  PetscInt choice_base_value = LOBATTO;
146 
147  // Read options from command line
148  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
149  CHKERR(ierr);
150  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
151  mesh_file_name, 255, &flg_file);
152 
153  CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
154  order, &order, PETSC_NULL);
155 
156  CHKERR PetscOptionsEList("-base", "approximation base", "", list_bases,
157  LASBASETOP, list_bases[choice_base_value],
158  &choice_base_value, PETSC_NULL);
159 
160  CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "",
161  test_nb, &test_nb, PETSC_NULL);
162 
163  CHKERR PetscOptionsBool("-my_is_partitioned",
164  "set if mesh is partitioned (this result that each "
165  "process keeps only one part of the mesh)",
166  "", is_partitioned, &is_partitioned, PETSC_NULL);
167 
168  CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
169  "", "block_conf.in", block_config_file, 255,
170  &flg_block_config);
171 
172  ierr = PetscOptionsEnd();
173  CHKERRG(ierr);
174 
175  // Throw error if file with mesh is not provided
176  if (flg_file != PETSC_TRUE) {
177  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
178  }
179 
180  // Create mesh database
181  moab::Core mb_instance;
182  moab::Interface &moab = mb_instance;
183 
184  // Create moab communicator
185  // Create separate MOAB communicator, it is duplicate of PETSc communicator.
186  // NOTE That this should eliminate potential communication problems between
187  // MOAB and PETSC functions.
188  MPI_Comm moab_comm_world;
189  MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
190  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
191  if (pcomm == NULL)
192  pcomm = new ParallelComm(&moab, moab_comm_world);
193 
194  // Read whole mesh or part of it if partitioned
195  if (is_partitioned == PETSC_TRUE) {
196  // This is a case of distributed mesh and algebra. In this case each
197  // processor keeps only one part of the problem.
198  const char *option;
199  option = "PARALLEL=READ_PART;"
200  "PARALLEL_RESOLVE_SHARED_ENTS;"
201  "PARTITION=PARALLEL_PARTITION;";
202  CHKERR moab.load_file(mesh_file_name, 0, option);
203  } else {
204  // If that case we have distributed algebra, i.e. assembly of vectors and
205  // matrices is in parallel, but whole mesh is stored on all processors.
206  // Solver and matrix scales well, however problem set-up of problem is
207  // not fully parallel.
208  const char *option;
209  option = "";
210  CHKERR moab.load_file(mesh_file_name, 0, option);
211  }
212 
213  // Create MoFEM database and link it to MoAB
214  MoFEM::Core core(moab);
215  MoFEM::Interface &m_field = core;
216 
217  // Print boundary conditions and material parameters
218  MeshsetsManager *meshsets_mng_ptr;
219  CHKERR m_field.getInterface(meshsets_mng_ptr);
220  CHKERR meshsets_mng_ptr->printDisplacementSet();
221  CHKERR meshsets_mng_ptr->printForceSet();
222  CHKERR meshsets_mng_ptr->printMaterialsSet();
223 
224  // Set bit refinement level to all entities (we do not refine mesh in this
225  // example so all entities have the same bit refinement level)
226  BitRefLevel bit_level0;
227  bit_level0.set(0);
228  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
229  0, 3, bit_level0);
230 
231  // Declare approximation fields
233  switch (choice_base_value) {
234  case LEGENDRE:
236  break;
237  case LOBATTO:
238  base = AINSWORTH_LOBATTO_BASE;
239  break;
240  case BERNSTEIN_BEZIER:
242  break;
243  default:
244  SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
245  };
246  CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3,
247  MB_TAG_DENSE, MF_ZERO);
248 
249  // We can use higher oder geometry to define body
250  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
251  3, MB_TAG_DENSE, MF_ZERO);
252 
253  // Declare problem
254 
255  // Add entities (by tets) to the field (all entities in the mesh, root_set
256  // = 0 )
257  CHKERR m_field.add_ents_to_field_by_dim(0, 3, "DISPLACEMENT");
258  CHKERR m_field.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
259 
260  // Set approximation order.
261  // See Hierarchical Finite Element Bases on Unstructured Tetrahedral
262  // Meshes.
263  CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
264  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
265  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
266  CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
267  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
269  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
270  else
271  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
272 
273  // Set order of approximation of geometry.
274  // Apply 2nd order only on skin (or in whole body)
275  auto setting_second_order_geometry = [&m_field]() {
277  // Setting geometry order everywhere
278  Range tets, edges;
279  CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
280  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
281  moab::Interface::UNION);
282 
283  // Setting 2nd geometry order only on skin
284  // Range tets, faces, edges;
285  // Skinner skin(&m_field.get_moab());
286  // CHKERR skin.find_skin(0,tets,false,faces);
287  // CHKERR m_field.get_moab().get_adjacencies(
288  // faces,1,false,edges,moab::Interface::UNION
289  // );
290  // CHKERR
291  // m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
292 
293  CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
294  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
295 
297  };
298  CHKERR setting_second_order_geometry();
299 
300  // Configure blocks by parsing config file. It allows setting
301  // approximation order for each block independently.
302  std::map<int, BlockOptionData> block_data;
303  auto setting_blocks_data_and_order_from_config_file =
304  [&m_field, &moab, &block_data, flg_block_config, block_config_file,
305  order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
307  if (flg_block_config) {
308  ifstream ini_file(block_config_file);
309  po::variables_map vm;
310  po::options_description config_file_options;
312  it)) {
313  std::ostringstream str_order;
314  str_order << "block_" << it->getMeshsetId()
315  << ".displacement_order";
316  config_file_options.add_options()(
317  str_order.str().c_str(),
318  po::value<int>(&block_data[it->getMeshsetId()].oRder)
319  ->default_value(order));
320 
321  std::ostringstream str_cond;
322  str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
323  config_file_options.add_options()(
324  str_cond.str().c_str(),
325  po::value<double>(&block_data[it->getMeshsetId()].yOung)
326  ->default_value(-1));
327 
328  std::ostringstream str_capa;
329  str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
330  config_file_options.add_options()(
331  str_capa.str().c_str(),
332  po::value<double>(&block_data[it->getMeshsetId()].pOisson)
333  ->default_value(-2));
334 
335  std::ostringstream str_init_temp;
336  str_init_temp << "block_" << it->getMeshsetId()
337  << ".initial_temperature";
338  config_file_options.add_options()(
339  str_init_temp.str().c_str(),
340  po::value<double>(&block_data[it->getMeshsetId()].initTemp)
341  ->default_value(0));
342  }
343  po::parsed_options parsed =
344  parse_config_file(ini_file, config_file_options, true);
345  store(parsed, vm);
346  po::notify(vm);
348  it)) {
349  if (block_data[it->getMeshsetId()].oRder == -1)
350  continue;
351  if (block_data[it->getMeshsetId()].oRder == order)
352  continue;
353  PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
354  it->getMeshsetId(),
355  block_data[it->getMeshsetId()].oRder);
356  Range block_ents;
357  CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
358  true);
359  Range ents_to_set_order;
360  CHKERR moab.get_adjacencies(block_ents, 3, false,
361  ents_to_set_order,
362  moab::Interface::UNION);
363  ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
364  CHKERR moab.get_adjacencies(block_ents, 2, false,
365  ents_to_set_order,
366  moab::Interface::UNION);
367  CHKERR moab.get_adjacencies(block_ents, 1, false,
368  ents_to_set_order,
369  moab::Interface::UNION);
370  CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
371  ents_to_set_order);
372 
373  CHKERR m_field.set_field_order(
374  ents_to_set_order, "DISPLACEMENT",
375  block_data[it->getMeshsetId()].oRder);
376  }
377  std::vector<std::string> additional_parameters;
378  additional_parameters =
379  collect_unrecognized(parsed.options, po::include_positional);
380  for (std::vector<std::string>::iterator vit =
381  additional_parameters.begin();
382  vit != additional_parameters.end(); vit++) {
383  CHKERR PetscPrintf(PETSC_COMM_WORLD,
384  "** WARNING Unrecognized option %s\n",
385  vit->c_str());
386  }
387  }
388 
389  // Update material parameters. Set material parameters block by
390  // block.
392  m_field, BLOCKSET | MAT_ELASTICSET, it)) {
393  const int id = it->getMeshsetId();
394  auto &bd = (*block_sets_ptr)[id];
395  if (block_data[id].yOung > 0)
396  bd.E = block_data[id].yOung;
397  if (block_data[id].pOisson >= -1)
398  bd.PoissonRatio = block_data[id].pOisson;
399  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Block %d\n", id);
400  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung modulus %3.4g\n",
401  bd.E);
402  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson ratio %3.4g\n",
403  bd.PoissonRatio);
404  }
405 
407  };
408 
409  // Add elastic element
410 
411  boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
412  boost::make_shared<std::map<int, BlockData>>();
413  CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
414  CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
415 
416  boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
417  new VolumeElementForcesAndSourcesCore(m_field));
418  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
419  new VolumeElementForcesAndSourcesCore(m_field));
420  fe_lhs_ptr->getRuleHook = VolRule();
421  fe_rhs_ptr->getRuleHook = VolRule();
422 
423  boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
424  new PrismFE(m_field));
425  boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
426  new PrismFE(m_field));
427 
428  CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
429  "DISPLACEMENT",
430  "MESH_NODE_POSITIONS", false);
431  CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
432  "DISPLACEMENT", "MESH_NODE_POSITIONS",
433  false, true);
435  prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
436  "MESH_NODE_POSITIONS", false, true, MBPRISM);
437 
438  // Add spring boundary condition applied on surfaces.
439  // This is only declaration not implementation.
440  CHKERR MetaSpringBC::addSpringElements(m_field, "DISPLACEMENT",
441  "MESH_NODE_POSITIONS");
442 
443  // Implementation of spring element
444  // Create new instances of face elements for springs
445  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
446  new FaceElementForcesAndSourcesCore(m_field));
447  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
448  new FaceElementForcesAndSourcesCore(m_field));
449 
450  CHKERR MetaSpringBC::setSpringOperators(m_field, fe_spring_lhs_ptr,
451  fe_spring_rhs_ptr, "DISPLACEMENT",
452  "MESH_NODE_POSITIONS");
453 
454  // Add body force element. This is only declaration of element. not its
455  // implementation.
456  CHKERR m_field.add_finite_element("BODY_FORCE");
457  CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
458  "DISPLACEMENT");
459  CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
460  "DISPLACEMENT");
461  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
462  "DISPLACEMENT");
463  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
464  "MESH_NODE_POSITIONS");
465 
467  m_field, BLOCKSET | BODYFORCESSET, it)) {
468  Range tets;
469  CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
470  true);
471  CHKERR m_field.add_ents_to_finite_element_by_dim(tets, 3, "BODY_FORCE");
472  }
473 
474  // Add Neumann forces, i.e. pressure or traction forces applied on body
475  // surface. This is only declaration not implementation.
476  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
477  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
478  CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
479 
480  // Add fluid pressure finite elements. This is special pressure on the
481  // surface from fluid, i.e. pressure which linearly change with the depth.
482  FluidPressure fluid_pressure_fe(m_field);
483  // This function only declare element. Element is implemented by operators
484  // in class FluidPressure.
485  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
486 
487  // Add elements for thermo elasticity if temperature field is defined.
488  ThermalStressElement thermal_stress_elem(m_field);
489  // Check if TEMP field exist, and then add element.
490  if (!m_field.check_field("TEMP")) {
491  bool add_temp_field = false;
493  if (block_data[it->getMeshsetId()].initTemp != 0) {
494  add_temp_field = true;
495  break;
496  }
497  }
498  if (add_temp_field) {
499  CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
500  MB_TAG_SPARSE, MF_ZERO);
501 
502  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
503 
504  CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
505  }
506  }
507  if (m_field.check_field("TEMP")) {
508  CHKERR thermal_stress_elem.addThermalStressElement(
509  "ELASTIC", "DISPLACEMENT", "TEMP");
510  }
511 
512  // All is declared, at this point build fields first,
513  CHKERR m_field.build_fields();
514  // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
515  // Projection10NodeCoordsOnField
516  // do the trick.
517  Projection10NodeCoordsOnField ent_method_material(m_field,
518  "MESH_NODE_POSITIONS");
519  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
520  if (m_field.check_field("TEMP")) {
522  if (block_data[it->getMeshsetId()].initTemp != 0) {
523  PetscPrintf(PETSC_COMM_WORLD, "Set block %d temperature to %3.2g\n",
524  it->getMeshsetId(),
525  block_data[it->getMeshsetId()].initTemp);
526  Range block_ents;
527  CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
528  Range vertices;
529  CHKERR moab.get_connectivity(block_ents, vertices, true);
530  CHKERR m_field.getInterface<FieldBlas>()->setField(
531  block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
532  "TEMP");
533  }
534  }
535  }
536 
537  // Build database for elements. Actual implementation of element is not need
538  // here, only elements has to be declared.
539  CHKERR m_field.build_finite_elements();
540  // Build adjacencies between elements and field entities
541  CHKERR m_field.build_adjacencies(bit_level0);
542 
543  // Register MOFEM DM implementation in PETSc
545 
546  // Create DM manager
547  DM dm;
548  CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
549  CHKERR DMSetType(dm, "MOFEM");
550  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
551  CHKERR DMSetFromOptions(dm);
552  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
553  // Add elements to DM manager
554  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
555  CHKERR DMMoFEMAddElement(dm, "SPRING");
556  CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
557  CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
558  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
559  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
560  CHKERR DMSetUp(dm);
561 
562  // Create matrices & vectors. Note that native PETSc DM interface is used,
563  // but under the PETSc interface MOFEM implementation is running.
564  Vec F, D, D0;
565  CHKERR DMCreateGlobalVector(dm, &F);
566  CHKERR VecDuplicate(F, &D);
567  CHKERR VecDuplicate(F, &D0);
568  Mat Aij;
569  CHKERR DMCreateMatrix(dm, &Aij);
570  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
571 
572  // Assign global matrix/vector contributed by springs
573  fe_spring_lhs_ptr->ksp_B = Aij;
574  fe_spring_rhs_ptr->ksp_f = F;
575 
576  // Zero vectors and matrices
577  CHKERR VecZeroEntries(F);
578  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
579  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
580  CHKERR VecZeroEntries(D);
581  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
582  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
583  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
584  CHKERR MatZeroEntries(Aij);
585 
586  // This controls how kinematic constrains are set, by blockset or nodeset.
587  // Cubit sets kinetic boundary conditions by blockset.
588  bool flag_cubit_disp = false;
590  m_field, NODESET | DISPLACEMENTSET, it)) {
591  flag_cubit_disp = true;
592  }
593 
594  // Below particular implementations of finite elements are used to assemble
595  // problem matrixes and vectors. Implementation of element does not change
596  // how element is declared.
597 
598  // Assemble Aij and F. Define Dirichlet bc element, which stets constrains
599  // to MatrixDouble and the right hand side vector.
600  boost::shared_ptr<FEMethod> dirichlet_bc_ptr;
601 
602  // if normally defined boundary conditions are not found, try to use
603  // DISPLACEMENT blockset. To implementations available here, depending how
604  // BC is defined on mesh file.
605  if (!flag_cubit_disp) {
606  dirichlet_bc_ptr =
607  boost::shared_ptr<FEMethod>(new DirichletSetFieldFromBlockWithFlags(
608  m_field, "DISPLACEMENT", "DISPLACEMENT", Aij, D0, F));
609  } else {
610  dirichlet_bc_ptr = boost::shared_ptr<FEMethod>(
611  new DirichletDisplacementBc(m_field, "DISPLACEMENT", Aij, D0, F));
612  }
613  // That sets Dirichlet bc objects that problem is linear, i.e. no newton
614  // (SNES) solver is run for this problem.
615  dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
616  dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
617 
618  // D0 vector will store initial displacements
619  CHKERR VecZeroEntries(D0);
620  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
621  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
622  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
623  // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
624  // implementation
625  // of particular dirichlet_bc.
626  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
627  // Set values from D0 on the field (on the mesh)
628 
629  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
630  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
631  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
632 
633  // Calculate residual forces as result of applied kinematic constrains. Run
634  // implementation
635  // of particular finite element implementation. Look how
636  // NonlinearElasticElement is implemented,
637  // in that case. We run NonlinearElasticElement with hook material.
638  // Calculate right hand side vector
639  fe_rhs_ptr->snes_f = F;
640  prism_fe_rhs_ptr->snes_f = F;
641  PetscPrintf(PETSC_COMM_WORLD, "Assemble external force vector ...");
642  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
643  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
644  PetscPrintf(PETSC_COMM_WORLD, " done\n");
645  // Assemble matrix
646  fe_lhs_ptr->snes_B = Aij;
647  prism_fe_lhs_ptr->snes_B = Aij;
648  PetscPrintf(PETSC_COMM_WORLD, "Calculate stiffness matrix ...");
649  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
650  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
651  PetscPrintf(PETSC_COMM_WORLD, " done\n");
652 
653  // Assemble springs
654  CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_lhs_ptr);
655  CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_rhs_ptr);
656 
657  // Assemble pressure and traction forces. Run particular implemented for do
658  // this, see
659  // MetaNeumannForces how this is implemented.
660  boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
661  CHKERR MetaNeumannForces::setMomentumFluxOperators(m_field, neumann_forces,
662  F, "DISPLACEMENT");
663 
664  {
665  boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
666  neumann_forces.begin();
667  for (; mit != neumann_forces.end(); mit++) {
668  CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
669  &mit->second->getLoopFe());
670  }
671  }
672  // Assemble forces applied to nodes, see implementation in NodalForce
673  boost::ptr_map<std::string, NodalForce> nodal_forces;
674  CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
675  "DISPLACEMENT");
676 
677  {
678  boost::ptr_map<std::string, NodalForce>::iterator fit =
679  nodal_forces.begin();
680  for (; fit != nodal_forces.end(); fit++) {
681  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
682  &fit->second->getLoopFe());
683  }
684  }
685  // Assemble edge forces
686  boost::ptr_map<std::string, EdgeForce> edge_forces;
687  CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
688  "DISPLACEMENT");
689  {
690  boost::ptr_map<std::string, EdgeForce>::iterator fit =
691  edge_forces.begin();
692  for (; fit != edge_forces.end(); fit++) {
693  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
694  &fit->second->getLoopFe());
695  }
696  }
697  // Assemble body forces, implemented in BodyForceConstantField
698  BodyForceConstantField body_forces_methods(m_field);
700  m_field, BLOCKSET | BODYFORCESSET, it)) {
701  CHKERR body_forces_methods.addBlock("DISPLACEMENT", F,
702  it->getMeshsetId());
703  }
704  CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
705  &body_forces_methods.getLoopFe());
706  // Assemble fluid pressure forces
707  CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
708  "DISPLACEMENT", F, false, true);
709 
710  CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
711  &fluid_pressure_fe.getLoopFe());
712  // Apply kinematic constrain to right hand side vector and matrix
713  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
714 
715  // Matrix View
716  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
717  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
718  // std::string wait;
719  // std::cin >> wait;
720 
721  // Set matrix positive defined and symmetric for Cholesky and icc
722  // pre-conditioner
723  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
724  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
725  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
726  CHKERR VecAssemblyBegin(F);
727  CHKERR VecAssemblyEnd(F);
728  CHKERR VecScale(F, -1);
729 
730  // Create solver
731  KSP solver;
732  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
733  CHKERR KSPSetDM(solver, dm);
734  CHKERR KSPSetFromOptions(solver);
735  CHKERR KSPSetOperators(solver, Aij, Aij);
736  // Setup multi-grid pre-conditioner if set from command line
737  {
738  // from PETSc example ex42.c
739  PetscBool same = PETSC_FALSE;
740  PC pc;
741  CHKERR KSPGetPC(solver, &pc);
742  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
743  if (same) {
744  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
745  CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
746  CHKERR PCSetFromOptions(pc);
747  } else {
748  // Operators are already set, do not use DM for doing that
749  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
750  }
751  }
752  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
753  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
754  // Set up solver
755  CHKERR KSPSetUp(solver);
756 
757  // Set up post-processor. It is some generic implementation of finite
758  // element.
759  PostProcVolumeOnRefinedMesh post_proc(m_field);
760  // Add operators to the elements, starting with some generic operators
761  CHKERR post_proc.generateReferenceElementMesh();
762  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
763  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
764  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
765  // Add problem specific operator on element to post-process stresses
766  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
767  m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "DISPLACEMENT",
768  post_proc.commonData, block_sets_ptr.get()));
769 
770  PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
771  CHKERR prism_post_proc.generateReferenceElementMesh();
772  CHKERR prism_post_proc.addFieldValuesPostProc("DISPLACEMENT");
773  CHKERR prism_post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
774  CHKERR prism_post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
775 
776  // Temperature field is defined on the mesh
777  if (m_field.check_field("TEMP")) {
778 
779  // Create thermal vector
780  Vec F_thermal;
781  CHKERR VecDuplicate(F, &F_thermal);
782 
783  // Set up implementation for calculation of thermal stress vector. Look
784  // how thermal stresses and vector is assembled in ThermalStressElement.
785  CHKERR thermal_stress_elem.setThermalStressRhsOperators(
786  "DISPLACEMENT", "TEMP", F_thermal);
787 
788  SeriesRecorder *recorder_ptr;
789  CHKERR m_field.getInterface(recorder_ptr);
790 
791  // Read time series and do thermo-elastic analysis, this is when time
792  // dependent
793  // temperature problem was run before on the mesh. It means that before
794  // non-stationary
795  // problem was solved for temperature and filed "TEMP" is stored for
796  // subsequent time
797  // steps in the recorder.
798  if (recorder_ptr->check_series("THEMP_SERIES")) {
799  // This is time dependent case, so loop of data series stored by tape
800  // recorder.
801  // Loop over time steps
802  for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
803  sit)) {
804  PetscPrintf(PETSC_COMM_WORLD, "Process step %d\n",
805  sit->get_step_number());
806  // Load field data for this time step
807  CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
808  sit->get_step_number());
809 
810  CHKERR VecZeroEntries(F_thermal);
811  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
812  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
813 
814  // Calculate the right-hand side vector as result of thermal stresses.
815  // It MetaNodalForces
816  // that on "ELASTIC" element data structure the element implementation
817  // in thermal_stress_elem
818  // is executed.
820  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
821 
822  // Assemble vector
823  CHKERR VecAssemblyBegin(F_thermal);
824  CHKERR VecAssemblyEnd(F_thermal);
825  // Accumulate ghost dofs
826  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
827  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
828 
829  // Calculate norm of vector and print values
830  PetscReal nrm_F;
831  CHKERR VecNorm(F, NORM_2, &nrm_F);
832 
833  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
834  PetscReal nrm_F_thermal;
835  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
836  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
837  nrm_F_thermal);
838 
839  CHKERR VecScale(F_thermal, -1);
840  // check this !!!
841  CHKERR VecAXPY(F_thermal, 1, F);
842 
843  // Set dirichlet boundary to thermal stresses vector
844  dirichlet_bc_ptr->snes_x = D;
845  dirichlet_bc_ptr->snes_f = F_thermal;
846  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
847 
848  // Solve problem
849  CHKERR KSPSolve(solver, F_thermal, D);
850  // Add boundary conditions vector
851  CHKERR VecAXPY(D, 1., D0);
852  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
853  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
854 
855  // Save data on the mesh
856  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
857 
858  // Save data on mesh
859  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
860  // Post-process results
861  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
862 
863  std::ostringstream o1;
864  o1 << "out_" << sit->step_number << ".h5m";
865  CHKERR post_proc.writeFile(o1.str().c_str());
866  }
867  } else {
868 
869  // This is a case when stationary problem for temperature was solved.
870  CHKERR VecZeroEntries(F_thermal);
871  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
872  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
873 
874  // Calculate the right-hand side vector with thermal stresses
876  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
877 
878  // Assemble vector
879  CHKERR VecAssemblyBegin(F_thermal);
880  CHKERR VecAssemblyEnd(F_thermal);
881 
882  // Accumulate ghost dofs
883  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
884  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
885 
886  // Calculate norm
887  PetscReal nrm_F;
888  CHKERR VecNorm(F, NORM_2, &nrm_F);
889 
890  PetscPrintf(PETSC_COMM_WORLD, "norm2 F = %6.4e\n", nrm_F);
891  PetscReal nrm_F_thermal;
892  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
893 
894  PetscPrintf(PETSC_COMM_WORLD, "norm2 F_thermal = %6.4e\n",
895  nrm_F_thermal);
896 
897  // Add thermal stress vector and other forces vector
898  CHKERR VecScale(F_thermal, -1);
899  CHKERR VecAXPY(F_thermal, 1, F);
900 
901  // Apply kinetic boundary conditions
902  dirichlet_bc_ptr->snes_x = D;
903  dirichlet_bc_ptr->snes_f = F_thermal;
904  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
905 
906  // Solve problem
907  CHKERR KSPSolve(solver, F_thermal, D);
908  CHKERR VecAXPY(D, 1., D0);
909 
910  // Update ghost values for solution vector
911  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
912  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
913 
914  // Save data on mesh
915  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
916  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
917  // Save results to file
918  PetscPrintf(PETSC_COMM_WORLD, "Write output file ...");
919  CHKERR post_proc.writeFile("out.h5m");
920  PetscPrintf(PETSC_COMM_WORLD, " done\n");
921  }
922 
923  // Destroy vector, no needed any more
924  CHKERR VecDestroy(&F_thermal);
925 
926  } else {
927  // Elastic analysis no temperature field
928  // Solve for vector D
929  CHKERR KSPSolve(solver, F, D);
930  // Add kinetic boundary conditions
931  CHKERR VecAXPY(D, 1., D0);
932  // Update ghost values
933  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
934  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
935  // Save data from vector on mesh
936  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
937  // Post-process results
938  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
939  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
940  // CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", &post_proc);
941  // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
942  PetscPrintf(PETSC_COMM_WORLD, "Write output file ..,");
943  CHKERR post_proc.writeFile("out.h5m");
944  CHKERR prism_post_proc.writeFile("prism_out.h5m");
945  PetscPrintf(PETSC_COMM_WORLD, " done\n");
946  }
947 
948  // Calculate elastic energy
949  auto calculate_strain_energy = [dm, &block_sets_ptr, test_nb]() {
951 
952  Vec v_energy;
953  CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
954  "MESH_NODE_POSITIONS", false, true,
955  &v_energy);
956 
957  // Print elastic energy
958  const double *eng_ptr;
959  CHKERR VecGetArrayRead(v_energy, &eng_ptr);
960  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n", *eng_ptr);
961 
962  switch (test_nb) {
963  case 1:
964  if (fabs(*eng_ptr - 17.129) > 1e-3)
965  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
966  "atom test diverged!");
967  break;
968  case 2:
969  if (fabs(*eng_ptr - 5.6475e-03) > 1e-3)
970  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
971  "atom test diverged!");
972  break;
973 
974  default:
975  break;
976  }
977 
978  CHKERR VecRestoreArrayRead(v_energy, &eng_ptr);
979  CHKERR VecDestroy(&v_energy);
981  };
982  CHKERR calculate_strain_energy();
983 
984  // Destroy matrices, vecors, solver and DM
985  CHKERR VecDestroy(&F);
986  CHKERR VecDestroy(&D);
987  CHKERR VecDestroy(&D0);
988  CHKERR MatDestroy(&Aij);
989  CHKERR KSPDestroy(&solver);
990  CHKERR DMDestroy(&dm);
991 
992  MPI_Comm_free(&moab_comm_world);
993  }
994  CATCH_ERRORS;
995 
997 
998  return 0;
999 }
Set integration rule.
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.
Implentation of thermal stress element.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
Managing BitRefLevels.
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.
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET)
virtual moab::Interface & get_moab()=0
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:434
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Interface for managing meshsets containing materials and boundary conditions.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:74
Core (interface) class.
Definition: Core.hpp:50
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
static char help[]
Definition: elasticity.cpp:72
char mesh_file_name[255]
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of spring element. Set operators to calculate LHS and RHS.
Projection of edge entities with one mid-node on hierarchical basis.
Postprocess on prism.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
Operator post-procesing stresses for Hook isotropic material.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:140
FieldApproximationBase
approximation base
Definition: definitions.h:144
virtual bool check_series(const std::string &name) const
check if series is in database
Body forces elements.
Definition: BodyForce.hpp:24
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Basic algebra on fields.
Definition: FieldBlas.hpp:34
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
Add boundary conditions form block set having 6 attributes.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:43
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
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:225
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
Set data structures of MG pre-conditioner via approximation orders.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:457
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, Vec *v_energy_ptr)
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:467
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:109
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
block name is "MAT_ELASTIC"
Definition: definitions.h:222
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: EdgeForce.hpp:109
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &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.
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 MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:971
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61

Variable Documentation

◆ help

char help[]
static
Initial value:
= "-my_block_config set block data\n"
"-my_order approximation order\n"
"-my_is_partitioned set if mesh is partitioned\n"
"\n"
Examples
elasticity.cpp.

Definition at line 72 of file elasticity.cpp.