v0.14.0
elasticity.cpp
Go to the documentation of this file.
1 /** \file elasticity.cpp
2  * \ingroup nonlinear_elastic_elem
3  * \example elasticity.cpp
4 
5  The example shows how to solve the linear elastic problem. An example can read
6  file with temperature field, then thermal stresses are included.
7 
8  What example can do:
9  - take into account temperature field, i.e. calculate thermal stresses and
10 deformation
11  - stationary and time depend field is considered
12  - take into account gravitational body forces
13  - take in account fluid pressure
14  - can work with higher order geometry definition
15  - works on distributed meshes
16  - multi-grid solver where each grid level is approximation order level
17  - each mesh block can have different material parameters and approximation
18 order
19 
20 See example how code can be used \cite jordi:2017,
21  \image html SquelaDamExampleByJordi.png "Example what you can do with this
22 code. Analysis of the arch dam of Susqueda, located in Catalonia (Spain)"
23 width=800px
24 
25  This is an example of application code; it does not show how elements are
26 implemented. Example presents how to:
27  - read mesh
28  - set-up problem
29  - run finite elements on the problem
30  - assemble matrices and vectors
31  - solve the linear system of equations
32  - save results
33 
34 
35  If you like to see how to implement finite elements, material, are other parts
36 of the code, look here;
37  - Hooke material, see \ref Hooke
38  - Thermal-stress assembly, see \ref ThermalElement
39  - Body forces element, see \ref BodyForceConstantField
40  - Fluid pressure element, see \ref FluidPressure
41  - The general implementation of an element for arbitrary Lagrangian-Eulerian
42  formulation for a nonlinear elastic problem is here \ref
43  NonlinearElasticElement. Here we limit ourselves to Hooke equation and fix
44  mesh, so the problem becomes linear. Not that elastic element is implemented
45  with automatic differentiation.
46 
47 */
48 
49 
50 
51 #include <BasicFiniteElements.hpp>
52 using namespace MoFEM;
53 
54 #include <Hooke.hpp>
55 using namespace boost::numeric;
56 
57 static char help[] = "-my_block_config set block data\n"
58  "-my_order approximation order\n"
59  "-my_is_partitioned set if mesh is partitioned\n"
60  "\n";
61 
62 struct BlockOptionData {
63  int oRder;
64  int iD;
65  double yOung;
66  double pOisson;
67  double initTemp;
68 
69  Range tRis;
70 
71  BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2), initTemp(0) {}
72 };
73 
74 using BlockData = NonlinearElasticElement::BlockData;
75 using MassBlockData = ConvectiveMassElement::BlockData;
76 using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator;
77 
78 /// Set integration rule
79 struct VolRule {
80  int operator()(int, int, int order) const { return 2 * order; }
81 };
82 
83 struct PrismFE : public FatPrismElementForcesAndSourcesCore {
84 
85  using FatPrismElementForcesAndSourcesCore::
86  FatPrismElementForcesAndSourcesCore;
87  int getRuleTrianglesOnly(int order);
88  int getRuleThroughThickness(int order);
89 };
90 
91 int PrismFE::getRuleTrianglesOnly(int order) { return 2 * order; };
92 int PrismFE::getRuleThroughThickness(int order) { return 2 * order; };
93 
94 using EdgeEle = MoFEM::EdgeElementForcesAndSourcesCore;
95 
96 int main(int argc, char *argv[]) {
97 
98  const string default_options = "-ksp_type gmres \n"
99  "-pc_type lu \n"
100  "-pc_factor_mat_solver_type mumps \n"
101  "-mat_mumps_icntl_20 0 \n"
102  "-ksp_monitor \n"
103  "-snes_type newtonls \n"
104  "-snes_linesearch_type basic \n"
105  "-snes_atol 1e-8 \n"
106  "-snes_rtol 1e-8 \n"
107  "-snes_monitor \n"
108  "-ts_monitor \n"
109  "-ts_type beuler \n";
110 
111  string param_file = "param_file.petsc";
112  if (!static_cast<bool>(ifstream(param_file))) {
113  std::ofstream file(param_file.c_str(), std::ios::ate);
114  if (file.is_open()) {
115  file << default_options;
116  file.close();
117  }
118  }
119 
120  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
121  auto core_log = logging::core::get();
122  core_log->add_sink(
123  LogManager::createSink(LogManager::getStrmWorld(), "ELASTIC"));
124  LogManager::setLog("ELASTIC");
125  MOFEM_LOG_TAG("ELASTIC", "elasticity")
126 
127  core_log->add_sink(
128  LogManager::createSink(LogManager::getStrmSync(), "ELASTIC_SYNC"));
129  LogManager::setLog("ELASTIC_SYNC");
130  MOFEM_LOG_TAG("ELASTIC_SYNC", "elastic_sync");
131 
132  try {
133 
134  PetscBool flg_block_config, flg_file;
135  char mesh_file_name[255];
136  char block_config_file[255];
137  PetscInt test_nb = 0;
138  PetscInt order = 2;
139  PetscBool is_partitioned = PETSC_FALSE;
140  PetscBool is_calculating_frequency = PETSC_FALSE;
141  PetscBool is_post_proc_volume = PETSC_TRUE;
142 
143  // Select base
144  enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
145  const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
146  "jacobi"};
147  PetscInt choice_base_value = LOBATTO;
148 
149  // Read options from command line
150  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
151  CHKERR(ierr);
152  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
153  mesh_file_name, 255, &flg_file);
154 
155  CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
156  order, &order, PETSC_NULL);
157 
158  CHKERR PetscOptionsEList("-base", "approximation base", "", list_bases,
159  LASBASETOP, list_bases[choice_base_value],
160  &choice_base_value, PETSC_NULL);
161 
162  CHKERR PetscOptionsInt("-is_atom_test", "ctest number", "", test_nb,
163  &test_nb, PETSC_NULL);
164 
165  CHKERR PetscOptionsBool("-my_is_partitioned",
166  "set if mesh is partitioned (this result that each "
167  "process keeps only one part of the mesh)",
168  "", is_partitioned, &is_partitioned, PETSC_NULL);
169 
170  CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
171  "", "block_conf.in", block_config_file, 255,
172  &flg_block_config);
173 
174  CHKERR PetscOptionsBool(
175  "-my_is_calculating_frequency", "set if frequency will be calculated",
176  "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
177 
178  CHKERR PetscOptionsBool("-my_is_post_proc_volume",
179  "if true post proc volume", "", is_post_proc_volume,
180  &is_post_proc_volume, PETSC_NULL);
181 
182  ierr = PetscOptionsEnd();
183  CHKERRG(ierr);
184 
185  // Throw error if file with mesh is not provided
186  if (flg_file != PETSC_TRUE) {
187  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
188  }
189 
190  // Create mesh database
191  moab::Core mb_instance;
192  moab::Interface &moab = mb_instance;
193 
194  // Create moab communicator
195  // Create separate MOAB communicator, it is duplicate of PETSc communicator.
196  // NOTE That this should eliminate potential communication problems between
197  // MOAB and PETSC functions.
198  MPI_Comm moab_comm_world;
199  MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
200  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
201  if (pcomm == NULL)
202  pcomm = new ParallelComm(&moab, moab_comm_world);
203 
204  // Read whole mesh or part of it if partitioned
205  if (is_partitioned == PETSC_TRUE) {
206  // This is a case of distributed mesh and algebra. In this case each
207  // processor keeps only one part of the problem.
208  const char *option;
209  option = "PARALLEL=READ_PART;"
210  "PARALLEL_RESOLVE_SHARED_ENTS;"
211  "PARTITION=PARALLEL_PARTITION;";
212  CHKERR moab.load_file(mesh_file_name, 0, option);
213  } else {
214  // If that case we have distributed algebra, i.e. assembly of vectors and
215  // matrices is in parallel, but whole mesh is stored on all processors.
216  // Solver and matrix scales well, however problem set-up of problem is
217  // not fully parallel.
218  const char *option;
219  option = "";
220  CHKERR moab.load_file(mesh_file_name, 0, option);
221  }
222 
223  // Create MoFEM database and link it to MoAB
224  MoFEM::Core core(moab);
225  MoFEM::Interface &m_field = core;
226 
227  // Print boundary conditions and material parameters
228  MeshsetsManager *meshsets_mng_ptr;
229  CHKERR m_field.getInterface(meshsets_mng_ptr);
230  CHKERR meshsets_mng_ptr->printDisplacementSet();
231  CHKERR meshsets_mng_ptr->printForceSet();
232  CHKERR meshsets_mng_ptr->printMaterialsSet();
233 
234  bool mesh_has_tets = false;
235  bool mesh_has_prisms = false;
236  int nb_tets = 0;
237  int nb_hexs = 0;
238  int nb_prisms = 0;
239 
240  CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets, true);
241  CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs, true);
242  CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms, true);
243 
244  mesh_has_tets = (nb_tets + nb_hexs) > 0;
245  mesh_has_prisms = nb_prisms > 0;
246 
247  // Set bit refinement level to all entities (we do not refine mesh in
248  // this example so all entities have the same bit refinement level)
249  BitRefLevel bit_level0;
250  bit_level0.set(0);
251  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
252  0, 3, bit_level0);
253 
254  // CHECK IF EDGE BLOCKSET EXIST AND IF IT IS ADD ALL ENTITIES FROM IT
255  // CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
256  // MESHSET_OF_EDGE_BLOCKSET, 1, bit_level0);
257 
258  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
259  if (bit->getName().compare(0, 3, "ROD") == 0) {
260  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
261  0, 1, bit_level0);
262  }
263  }
264 
265  // Declare approximation fields
266  FieldApproximationBase base = NOBASE;
267  switch (choice_base_value) {
268  case LEGENDRE:
269  base = AINSWORTH_LEGENDRE_BASE;
270  break;
271  case LOBATTO:
272  base = AINSWORTH_LOBATTO_BASE;
273  break;
274  case BERNSTEIN_BEZIER:
275  base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
276  break;
277  case JACOBI:
278  base = DEMKOWICZ_JACOBI_BASE;
279  break;
280  default:
281  SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
282  };
283  CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_DENSE,
284  MF_ZERO);
285 
286  // We can use higher oder geometry to define body
287  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, base, 3, MB_TAG_DENSE,
288  MF_ZERO);
289 
290  // Declare problem
291 
292  // Add entities (by tets) to the field (all entities in the mesh, root_set
293  // = 0 )
294  CHKERR m_field.add_ents_to_field_by_dim(0, 3, "DISPLACEMENT");
295  CHKERR m_field.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
296 
297  // Get all edges in the mesh
298  Range all_edges;
299  CHKERR m_field.get_moab().get_entities_by_type(0, MBEDGE, all_edges, true);
300 
301  // Get edges associated with simple rod
302  Range edges_in_simple_rod;
303  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
304  if (bit->getName().compare(0, 3, "ROD") == 0) {
305  Range edges;
306  CHKERR m_field.get_moab().get_entities_by_type(bit->getMeshset(),
307  MBEDGE, edges, true);
308  edges_in_simple_rod.merge(edges);
309  }
310  }
311 
312  CHKERR m_field.add_ents_to_field_by_type(edges_in_simple_rod, MBEDGE,
313  "DISPLACEMENT");
314 
315  // Set order of edge in rod to be 1
316  CHKERR m_field.set_field_order(edges_in_simple_rod, "DISPLACEMENT", 1);
317 
318  // Get remaining edges (not associated with simple rod) to set order
319  Range edges_to_set_order;
320  edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
321 
322  // Set approximation order.
323  // See Hierarchical Finite Element Bases on Unstructured Tetrahedral
324  // Meshes.
325  CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
326  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
327  CHKERR m_field.set_field_order(0, MBHEX, "DISPLACEMENT", order);
328  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
329  CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
330  CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
331  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
332 
333  if (base == AINSWORTH_BERNSTEIN_BEZIER_BASE)
334  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
335  else
336  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
337 
338  // Set order of approximation of geometry.
339  // Apply 2nd order only on skin (or in whole body)
340  auto setting_second_order_geometry = [&m_field]() {
341  MoFEMFunctionBegin;
342  // Setting geometry order everywhere
343  Range tets, edges;
344  CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
345  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
346  moab::Interface::UNION);
347 
348  // Setting 2nd geometry order only on skin
349  // Range tets, faces, edges;
350  // Skinner skin(&m_field.get_moab());
351  // CHKERR skin.find_skin(0,tets,false,faces);
352  // CHKERR m_field.get_moab().get_adjacencies(
353  // faces,1,false,edges,moab::Interface::UNION
354  // );
355  // CHKERR
356  // m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
357 
358  CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
359  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
360 
361  MoFEMFunctionReturn(0);
362  };
363  CHKERR setting_second_order_geometry();
364 
365  // Configure blocks by parsing config file. It allows setting
366  // approximation order for each block independently.
367  std::map<int, BlockOptionData> block_data;
368  auto setting_blocks_data_and_order_from_config_file =
369  [&m_field, &moab, &block_data, flg_block_config, block_config_file,
370  order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
371  MoFEMFunctionBegin;
372  if (flg_block_config) {
373  ifstream ini_file(block_config_file);
374  po::variables_map vm;
375  po::options_description config_file_options;
376  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET,
377  it)) {
378  std::ostringstream str_order;
379  str_order << "block_" << it->getMeshsetId()
380  << ".displacement_order";
381  config_file_options.add_options()(
382  str_order.str().c_str(),
383  po::value<int>(&block_data[it->getMeshsetId()].oRder)
384  ->default_value(order));
385 
386  std::ostringstream str_cond;
387  str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
388  config_file_options.add_options()(
389  str_cond.str().c_str(),
390  po::value<double>(&block_data[it->getMeshsetId()].yOung)
391  ->default_value(-1));
392 
393  std::ostringstream str_capa;
394  str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
395  config_file_options.add_options()(
396  str_capa.str().c_str(),
397  po::value<double>(&block_data[it->getMeshsetId()].pOisson)
398  ->default_value(-2));
399 
400  std::ostringstream str_init_temp;
401  str_init_temp << "block_" << it->getMeshsetId()
402  << ".initial_temperature";
403  config_file_options.add_options()(
404  str_init_temp.str().c_str(),
405  po::value<double>(&block_data[it->getMeshsetId()].initTemp)
406  ->default_value(0));
407  }
408  po::parsed_options parsed =
409  parse_config_file(ini_file, config_file_options, true);
410  store(parsed, vm);
411  po::notify(vm);
412  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET,
413  it)) {
414  if (block_data[it->getMeshsetId()].oRder == -1)
415  continue;
416  if (block_data[it->getMeshsetId()].oRder == order)
417  continue;
418  MOFEM_LOG_C("ELASTIC", Sev::inform, "Set block %d order to %d",
419  it->getMeshsetId(),
420  block_data[it->getMeshsetId()].oRder);
421  Range block_ents;
422  CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
423  true);
424  Range ents_to_set_order;
425  CHKERR moab.get_adjacencies(block_ents, 3, false,
426  ents_to_set_order,
427  moab::Interface::UNION);
428  ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
429  CHKERR moab.get_adjacencies(block_ents, 2, false,
430  ents_to_set_order,
431  moab::Interface::UNION);
432  CHKERR moab.get_adjacencies(block_ents, 1, false,
433  ents_to_set_order,
434  moab::Interface::UNION);
435  CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
436  ents_to_set_order);
437 
438  CHKERR m_field.set_field_order(
439  ents_to_set_order, "DISPLACEMENT",
440  block_data[it->getMeshsetId()].oRder);
441  }
442  std::vector<std::string> additional_parameters;
443  additional_parameters =
444  collect_unrecognized(parsed.options, po::include_positional);
445  for (std::vector<std::string>::iterator vit =
446  additional_parameters.begin();
447  vit != additional_parameters.end(); vit++) {
448  MOFEM_LOG_C("ELASTIC", Sev::warning, "Unrecognized option %s",
449  vit->c_str());
450  }
451  }
452 
453  // Update material parameters. Set material parameters block by
454  // block.
455  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
456  m_field, BLOCKSET | MAT_ELASTICSET, it)) {
457  const int id = it->getMeshsetId();
458  auto &bd = (*block_sets_ptr)[id];
459  if (block_data[id].yOung > 0)
460  bd.E = block_data[id].yOung;
461  if (block_data[id].pOisson >= -1)
462  bd.PoissonRatio = block_data[id].pOisson;
463  MOFEM_LOG_C("ELASTIC", Sev::inform, "Block %d", id);
464  MOFEM_LOG_C("ELASTIC", Sev::inform, "\tYoung modulus %3.4g", bd.E);
465  MOFEM_LOG_C("ELASTIC", Sev::inform, "\tPoisson ratio %3.4g",
466  bd.PoissonRatio);
467  }
468 
469  MoFEMFunctionReturn(0);
470  };
471 
472  // Add elastic element
473 
474  boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
475  boost::make_shared<std::map<int, HookeElement::BlockData>>();
476  CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
477  CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
478 
479  boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
480  boost::make_shared<std::map<int, MassBlockData>>();
481  CHKERR ConvectiveMassElement::setBlocks(m_field, mass_block_sets_ptr);
482 
483  auto fe_lhs_ptr =
484  boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
485  auto fe_rhs_ptr =
486  boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
487  fe_lhs_ptr->getRuleHook = VolRule();
488  fe_rhs_ptr->getRuleHook = VolRule();
489 
490  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_lhs_ptr, true, false, false,
491  false);
492  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_rhs_ptr, true, false, false,
493  false);
494 
495  boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
496  new PrismFE(m_field));
497  boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
498  new PrismFE(m_field));
499 
500  CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
501  "DISPLACEMENT",
502  "MESH_NODE_POSITIONS", false);
503 
504  auto add_skin_element_for_post_processing = [&]() {
505  MoFEMFunctionBegin;
506  Range elastic_element_ents;
507  CHKERR m_field.get_finite_element_entities_by_dimension(
508  "ELASTIC", 3, elastic_element_ents);
509  Skinner skin(&m_field.get_moab());
510  Range skin_faces; // skin faces from 3d ents
511  CHKERR skin.find_skin(0, elastic_element_ents, false, skin_faces);
512  Range proc_skin;
513  if (is_partitioned) {
514  CHKERR pcomm->filter_pstatus(skin_faces,
515  PSTATUS_SHARED | PSTATUS_MULTISHARED,
516  PSTATUS_NOT, -1, &proc_skin);
517  } else {
518  proc_skin = skin_faces;
519  }
520  CHKERR m_field.add_finite_element("POST_PROC_SKIN");
521  CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
522  "DISPLACEMENT");
523  CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
524  "DISPLACEMENT");
525  CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
526  "DISPLACEMENT");
527  if (m_field.check_field("TEMP")) {
528  // CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
529  // "TEMP");
530  // CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
531  // "TEMP");
532  CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
533  "TEMP");
534  }
535  CHKERR m_field.modify_finite_element_add_field_data(
536  "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
537  CHKERR m_field.add_ents_to_finite_element_by_dim(proc_skin, 2,
538  "POST_PROC_SKIN");
539  MoFEMFunctionReturn(0);
540  };
541  CHKERR add_skin_element_for_post_processing();
542 
543  auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
544  if (mesh_has_tets) {
545  CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
546  "DISPLACEMENT", "MESH_NODE_POSITIONS",
547  false, true, MBTET, data_at_pts);
548  }
549  if (mesh_has_prisms) {
550  CHKERR HookeElement::setOperators(
551  prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
552  "MESH_NODE_POSITIONS", false, true, MBPRISM, data_at_pts);
553  }
554 
555  if (test_nb == 4) {
556 
557  auto thermal_strain =
558  [](FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords) {
559  FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
560  constexpr double alpha = 1;
561  FTensor::Index<'i', 3> i;
562  FTensor::Index<'k', 3> j;
563  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
564  t_thermal_strain(i, j) = alpha * t_coords(2) * t_kd(i, j);
565  return t_thermal_strain;
566  };
567 
568  fe_rhs_ptr->getOpPtrVector().push_back(
569  new HookeElement::OpAnalyticalInternalStrain_dx<0>(
570  "DISPLACEMENT", data_at_pts, thermal_strain));
571  }
572 
573  boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
574  new VolumeElementForcesAndSourcesCore(m_field));
575 
576  for (auto &sit : *block_sets_ptr) {
577  for (auto &mit : *mass_block_sets_ptr) {
578  fe_mass_ptr->getOpPtrVector().push_back(
579  new HookeElement::OpCalculateMassMatrix("DISPLACEMENT",
580  "DISPLACEMENT", sit.second,
581  mit.second, data_at_pts));
582  }
583  }
584 
585  // Add spring boundary condition applied on surfaces.
586  // This is only declaration not implementation.
587  CHKERR MetaSpringBC::addSpringElements(m_field, "DISPLACEMENT",
588  "MESH_NODE_POSITIONS");
589 
590  // Implementation of spring element
591  // Create new instances of face elements for springs
592  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
593  new FaceElementForcesAndSourcesCore(m_field));
594  boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
595  new FaceElementForcesAndSourcesCore(m_field));
596 
597  CHKERR MetaSpringBC::setSpringOperators(m_field, fe_spring_lhs_ptr,
598  fe_spring_rhs_ptr, "DISPLACEMENT",
599  "MESH_NODE_POSITIONS");
600 
601  // Add Simple Rod elements
602  // This is only declaration not implementation.
603  CHKERR MetaSimpleRodElement::addSimpleRodElements(m_field, "DISPLACEMENT",
604  "MESH_NODE_POSITIONS");
605 
606  // CHKERR m_field.add_ents_to_finite_element_by_type(edges_in_simple_rod,
607  // MBEDGE, "SIMPLE_ROD");
608 
609  // Implementation of Simple Rod element
610  // Create new instances of edge elements for Simple Rod
611  boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(new EdgeEle(m_field));
612  boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(new EdgeEle(m_field));
613 
614 
615  CHKERR MetaSimpleRodElement::setSimpleRodOperators(
616  m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
617  "MESH_NODE_POSITIONS");
618 
619  // Add body force element. This is only declaration of element. not its
620  // implementation.
621  CHKERR m_field.add_finite_element("BODY_FORCE");
622  CHKERR m_field.modify_finite_element_add_field_row("BODY_FORCE",
623  "DISPLACEMENT");
624  CHKERR m_field.modify_finite_element_add_field_col("BODY_FORCE",
625  "DISPLACEMENT");
626  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
627  "DISPLACEMENT");
628  CHKERR m_field.modify_finite_element_add_field_data("BODY_FORCE",
629  "MESH_NODE_POSITIONS");
630 
631  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
632  m_field, BLOCKSET | BODYFORCESSET, it)) {
633  Range tets;
634  CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
635  true);
636  CHKERR m_field.add_ents_to_finite_element_by_dim(tets, 3, "BODY_FORCE");
637  }
638  CHKERR m_field.build_finite_elements("BODY_FORCE");
639 
640  // Add Neumann forces, i.e. pressure or traction forces applied on body
641  // surface. This is only declaration not implementation.
642  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
643  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
644  CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
645 
646  // Add fluid pressure finite elements. This is special pressure on the
647  // surface from fluid, i.e. pressure which linearly change with the depth.
648  FluidPressure fluid_pressure_fe(m_field);
649  // This function only declare element. Element is implemented by operators
650  // in class FluidPressure.
651  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
652 
653  // Add elements for thermo elasticity if temperature field is defined.
654  ThermalStressElement thermal_stress_elem(m_field);
655  // Check if TEMP field exist, and then add element.
656  if (!m_field.check_field("TEMP")) {
657  bool add_temp_field = false;
658  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
659  if (block_data[it->getMeshsetId()].initTemp != 0) {
660  add_temp_field = true;
661  break;
662  }
663  }
664  if (add_temp_field) {
665  CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
666  MB_TAG_SPARSE, MF_ZERO);
667 
668  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
669  CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
670  }
671  }
672  if (m_field.check_field("TEMP")) {
673  CHKERR thermal_stress_elem.addThermalStressElement(
674  "ELASTIC", "DISPLACEMENT", "TEMP");
675  }
676 
677  // All is declared, at this point build fields first,
678  CHKERR m_field.build_fields();
679  // If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
680  // Projection10NodeCoordsOnField
681  // do the trick.
682  Projection10NodeCoordsOnField ent_method_material(m_field,
683  "MESH_NODE_POSITIONS");
684  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
685  if (m_field.check_field("TEMP")) {
686  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
687  if (block_data[it->getMeshsetId()].initTemp != 0) {
688  MOFEM_LOG_C("ELASTIC", Sev::inform,
689  "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
690  block_data[it->getMeshsetId()].initTemp);
691  Range block_ents;
692  CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
693  Range vertices;
694  CHKERR moab.get_connectivity(block_ents, vertices, true);
695  CHKERR m_field.getInterface<FieldBlas>()->setField(
696  block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
697  "TEMP");
698  }
699  }
700  }
701 
702  // Build database for elements. Actual implementation of element is not need
703  // here, only elements has to be declared.
704  CHKERR m_field.build_finite_elements();
705  // Build adjacencies between elements and field entities
706  CHKERR m_field.build_adjacencies(bit_level0);
707 
708  // Register MOFEM DM implementation in PETSc
709  CHKERR DMRegister_MGViaApproxOrders("MOFEM");
710 
711  // Create DM manager
712  auto dm = createDM(PETSC_COMM_WORLD, "MOFEM");
713  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
714  CHKERR DMSetFromOptions(dm);
715  CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
716  // Add elements to DM manager
717  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
718  CHKERR DMMoFEMAddElement(dm, "SPRING");
719  CHKERR DMMoFEMAddElement(dm, "SIMPLE_ROD");
720  CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
721  CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
722  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
723  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
724  CHKERR DMMoFEMAddElement(dm, "POST_PROC_SKIN");
725  CHKERR DMSetUp(dm);
726 
727  // Create matrices & vectors. Note that native PETSc DM interface is used,
728  // but under the PETSc interface MoFEM implementation is running.
729  SmartPetscObj<Vec> F;
730  CHKERR DMCreateGlobalVector_MoFEM(dm, F);
731  auto D = vectorDuplicate(F);
732  auto D0 = vectorDuplicate(F);
733  SmartPetscObj<Mat> Aij;
734  CHKERR DMCreateMatrix_MoFEM(dm, Aij);
735  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
736 
737  // Initialise mass matrix
738  SmartPetscObj<Mat> Mij;
739  if (is_calculating_frequency == PETSC_TRUE) {
740  Mij = matDuplicate(Aij, MAT_DO_NOT_COPY_VALUES);
741  CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
742  // MatView(Mij, PETSC_VIEWER_STDOUT_SELF);
743  }
744 
745  // Assign global matrix/vector contributed by springs
746  fe_spring_lhs_ptr->ksp_B = Aij;
747  fe_spring_rhs_ptr->ksp_f = F;
748 
749  // Assign global matrix/vector contributed by Simple Rod
750  fe_simple_rod_lhs_ptr->ksp_B = Aij;
751  fe_simple_rod_rhs_ptr->ksp_f = F;
752 
753  // Zero vectors and matrices
754  CHKERR VecZeroEntries(F);
755  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
756  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
757  CHKERR VecZeroEntries(D);
758  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
759  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
760  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
761  CHKERR MatZeroEntries(Aij);
762 
763  // Below particular implementations of finite elements are used to assemble
764  // problem matrixes and vectors. Implementation of element does not change
765  // how element is declared.
766 
767  // Assemble Aij and F. Define Dirichlet bc element, which sets constrains
768  // to MatrixDouble and the right hand side vector.
769 
770  // if normally defined boundary conditions are not found,
771  // DirichletDisplacementBc will try to use DISPLACEMENT blockset. Two
772  // implementations are available, depending how BC is defined on mesh file.
773  auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
774  m_field, "DISPLACEMENT", Aij, D0, F);
775 
776  // That sets Dirichlet bc objects that problem is linear, i.e. no newton
777  // (SNES) solver is run for this problem.
778  dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
779  dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
780 
781  // D0 vector will store initial displacements
782  CHKERR VecZeroEntries(D0);
783  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
784  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
785  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
786  // Run dirichlet_bc, from that on the mesh set values in vector D0. Run
787  // implementation
788  // of particular dirichlet_bc.
789  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
790  // Set values from D0 on the field (on the mesh)
791 
792  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
793  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
794  CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
795 
796  // Calculate residual forces as result of applied kinematic constrains. Run
797  // implementation
798  // of particular finite element implementation. Look how
799  // NonlinearElasticElement is implemented,
800  // in that case. We run NonlinearElasticElement with hook material.
801  // Calculate right hand side vector
802  fe_rhs_ptr->snes_f = F;
803  prism_fe_rhs_ptr->snes_f = F;
804  MOFEM_LOG("ELASTIC", Sev::inform) << "Assemble external force vector ...";
805  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
806  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
807  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
808  // Assemble matrix
809  fe_lhs_ptr->snes_B = Aij;
810  prism_fe_lhs_ptr->snes_B = Aij;
811  MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate stiffness matrix ...";
812  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
813  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
814  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
815 
816  // Assemble springs
817  CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_lhs_ptr);
818  CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_rhs_ptr);
819 
820  // Assemble Simple Rod
821  CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_lhs_ptr);
822  CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_rhs_ptr);
823 
824  if (is_calculating_frequency == PETSC_TRUE) {
825  // Assemble mass matrix
826  fe_mass_ptr->snes_B = Mij;
827  MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate mass matrix ...";
828  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_mass_ptr);
829  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
830  }
831 
832  // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
833 
834  // Assemble pressure and traction forces. Run particular implemented for do
835  // this, see
836  // MetaNeumannForces how this is implemented.
837  boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
838  CHKERR MetaNeumannForces::setMomentumFluxOperators(m_field, neumann_forces,
839  F, "DISPLACEMENT");
840 
841  {
842  boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
843  neumann_forces.begin();
844  for (; mit != neumann_forces.end(); mit++) {
845  CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
846  &mit->second->getLoopFe());
847  }
848  }
849  // Assemble forces applied to nodes, see implementation in NodalForce
850  boost::ptr_map<std::string, NodalForce> nodal_forces;
851  CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
852  "DISPLACEMENT");
853 
854  {
855  boost::ptr_map<std::string, NodalForce>::iterator fit =
856  nodal_forces.begin();
857  for (; fit != nodal_forces.end(); fit++) {
858  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
859  &fit->second->getLoopFe());
860  }
861  }
862  // Assemble edge forces
863  boost::ptr_map<std::string, EdgeForce> edge_forces;
864  CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
865  "DISPLACEMENT");
866  {
867  auto fit = edge_forces.begin();
868  for (; fit != edge_forces.end(); fit++) {
869  auto &fe = fit->second->getLoopFe();
870  CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(), &fe);
871  }
872  }
873  // Assemble body forces, implemented in BodyForceConstantField
874  BodyForceConstantField body_forces_methods(m_field);
875  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
876  m_field, BLOCKSET | BODYFORCESSET, it)) {
877  CHKERR body_forces_methods.addBlock("DISPLACEMENT", F,
878  it->getMeshsetId());
879  }
880  CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
881  &body_forces_methods.getLoopFe());
882  // Assemble fluid pressure forces
883  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fluid_pressure_fe.getLoopFe(),
884  false, false);
885  CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
886  "DISPLACEMENT", F, false, true);
887 
888  CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
889  &fluid_pressure_fe.getLoopFe());
890  // Apply kinematic constrain to right hand side vector and matrix
891  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
892 
893  // Matrix View
894  PetscViewerPushFormat(
895  PETSC_VIEWER_STDOUT_SELF,
896  PETSC_VIEWER_ASCII_MATLAB); /// PETSC_VIEWER_ASCII_DENSE,
897  /// PETSC_VIEWER_ASCII_MATLAB
898  // MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
899  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
900  // std::string wait;
901  // std::cin >> wait;
902 
903  if (is_calculating_frequency == PETSC_TRUE) {
904  CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
905  CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
906  }
907 
908  // Set matrix positive defined and symmetric for Cholesky and icc
909  // pre-conditioner
910 
911  CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
912  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
913  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
914  CHKERR VecAssemblyBegin(F);
915  CHKERR VecAssemblyEnd(F);
916  CHKERR VecScale(F, -1);
917 
918  // Create solver
919  auto solver = createKSP(PETSC_COMM_WORLD);
920  CHKERR KSPSetDM(solver, dm);
921  CHKERR KSPSetFromOptions(solver);
922  CHKERR KSPSetOperators(solver, Aij, Aij);
923 
924  // Setup multi-grid pre-conditioner if set from command line
925  {
926  // from PETSc example ex42.c
927  PetscBool same = PETSC_FALSE;
928  PC pc;
929  CHKERR KSPGetPC(solver, &pc);
930  PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
931  if (same) {
932  PCMGSetUpViaApproxOrdersCtx pc_ctx(dm, Aij, true);
933  CHKERR PCMGSetUpViaApproxOrders(pc, &pc_ctx);
934  CHKERR PCSetFromOptions(pc);
935  } else {
936  // Operators are already set, do not use DM for doing that
937  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
938  }
939  }
940  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
941  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
942  // Set up solver
943  CHKERR KSPSetUp(solver);
944 
945  auto set_post_proc_skin = [&](auto &post_proc_skin) {
946  MoFEMFunctionBegin;
947  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false,
948  false);
949  CHKERR post_proc_skin.generateReferenceElementMesh();
950  CHKERR post_proc_skin.addFieldValuesPostProc("DISPLACEMENT");
951  CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
952  CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
953  "DISPLACEMENT", "ELASTIC", data_at_pts->hMat, true);
954  CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
955  "MESH_NODE_POSITIONS", "ELASTIC", data_at_pts->HMat, false);
956  if (m_field.check_field("TEMP")) {
957  CHKERR post_proc_skin.addFieldValuesPostProc("TEMP");
958  CHKERR post_proc_skin.addFieldValuesGradientPostProc("TEMP");
959  }
960  post_proc_skin.getOpPtrVector().push_back(
961  new HookeElement::OpPostProcHookeElement<
962  FaceElementForcesAndSourcesCore>(
963  "DISPLACEMENT", data_at_pts, *block_sets_ptr,
964  post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
965  true));
966  MoFEMFunctionReturn(0);
967  };
968 
969  auto set_post_proc_tets = [&](auto &post_proc) {
970  MoFEMFunctionBegin;
971  // Add operators to the elements, starting with some generic operators
972  CHKERR post_proc.generateReferenceElementMesh();
973  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, true, false, false,
974  false);
975  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
976  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
977  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
978  if (m_field.check_field("TEMP")) {
979  CHKERR post_proc.addFieldValuesPostProc("TEMP");
980  CHKERR post_proc.addFieldValuesGradientPostProc("TEMP");
981  }
982  // Add problem specific operator on element to post-process stresses
983  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
984  m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
985  "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
986  MoFEMFunctionReturn(0);
987  };
988 
989  auto set_post_proc_edge = [&](auto &post_proc_edge) {
990  MoFEMFunctionBegin;
991  CHKERR post_proc_edge.generateReferenceElementMesh();
992  CHKERR post_proc_edge.addFieldValuesPostProc("DISPLACEMENT");
993  MoFEMFunctionReturn(0);
994  };
995 
996  auto set_post_proc_prisms = [&](auto &prism_post_proc) {
997  MoFEMFunctionBegin;
998  CHKERR prism_post_proc.generateReferenceElementMesh();
999  boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
1000  prism_post_proc.getOpPtrVector().push_back(
1001  new OpCalculateInvJacForFatPrism(inv_jac_ptr));
1002  prism_post_proc.getOpPtrVector().push_back(
1003  new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
1004  CHKERR prism_post_proc.addFieldValuesPostProc("DISPLACEMENT");
1005  CHKERR prism_post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1006  CHKERR prism_post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
1007  prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
1008  m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
1009  "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
1010  MoFEMFunctionReturn(0);
1011  };
1012 
1013  PostProcFaceOnRefinedMesh post_proc_skin(m_field);
1014  PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
1015  PostProcEdgeOnRefinedMesh post_proc_edge(m_field);
1016  PostProcVolumeOnRefinedMesh post_proc(m_field);
1017 
1018  CHKERR set_post_proc_skin(post_proc_skin);
1019  CHKERR set_post_proc_tets(post_proc);
1020  CHKERR set_post_proc_prisms(prism_post_proc);
1021  CHKERR set_post_proc_edge(post_proc_edge);
1022 
1023  PetscBool field_eval_flag = PETSC_FALSE;
1024  std::array<double, 3> field_eval_coords;
1025  boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> field_eval_data;
1026  PetscInt coords_dim = 3;
1027  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1028  field_eval_coords.data(), &coords_dim,
1029  &field_eval_flag);
1030 
1031  auto scalar_field_ptr = boost::make_shared<VectorDouble>();
1032  auto vector_field_ptr = boost::make_shared<MatrixDouble>();
1033  auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
1034 
1035  if (field_eval_flag) {
1036  field_eval_data = m_field.getInterface<FieldEvaluatorInterface>()
1037  ->getData<VolumeElementForcesAndSourcesCore>();
1038  CHKERR m_field.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1039  field_eval_data, "ELASTIC");
1040 
1041  field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1042  auto no_rule = [](int, int, int) { return -1; };
1043 
1044  auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1045  field_eval_fe_ptr->getRuleHook = no_rule;
1046 
1047  if (m_field.check_field("TEMP")) {
1048  field_eval_fe_ptr->getOpPtrVector().push_back(
1049  new OpCalculateScalarFieldValues("TEMP", scalar_field_ptr));
1050  }
1051  field_eval_fe_ptr->getOpPtrVector().push_back(
1052  new OpCalculateVectorFieldValues<3>("DISPLACEMENT", vector_field_ptr));
1053  field_eval_fe_ptr->getOpPtrVector().push_back(
1054  new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENT", tensor_field_ptr));
1055  }
1056 
1057  // Temperature field is defined on the mesh
1058  if (m_field.check_field("TEMP")) {
1059 
1060  // Create thermal vector
1061  Vec F_thermal;
1062  CHKERR VecDuplicate(F, &F_thermal);
1063 
1064  // Set up implementation for calculation of thermal stress vector. Look
1065  // how thermal stresses and vector is assembled in ThermalStressElement.
1066  CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1067  "DISPLACEMENT", "TEMP", F_thermal);
1068 
1069  SeriesRecorder *recorder_ptr;
1070  CHKERR m_field.getInterface(recorder_ptr);
1071 
1072  // Read time series and do thermo-elastic analysis, this is when time
1073  // dependent
1074  // temperature problem was run before on the mesh. It means that before
1075  // non-stationary
1076  // problem was solved for temperature and filed "TEMP" is stored for
1077  // subsequent time
1078  // steps in the recorder.
1079  if (recorder_ptr->check_series("THEMP_SERIES")) {
1080  // This is time dependent case, so loop of data series stored by tape
1081  // recorder.
1082  // Loop over time steps
1083  for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
1084  sit)) {
1085  MOFEM_LOG_C("ELASTIC", Sev::inform, "Process step %d",
1086  sit->get_step_number());
1087  // Load field data for this time step
1088  CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
1089  sit->get_step_number());
1090 
1091  CHKERR VecZeroEntries(F_thermal);
1092  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1093  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1094 
1095  // Calculate the right-hand side vector as result of thermal stresses.
1096  // It MetaNodalForces
1097  // that on "ELASTIC" element data structure the element implementation
1098  // in thermal_stress_elem
1099  // is executed.
1100  CHKERR DMoFEMLoopFiniteElements(
1101  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1102 
1103  // Assemble vector
1104  CHKERR VecAssemblyBegin(F_thermal);
1105  CHKERR VecAssemblyEnd(F_thermal);
1106  // Accumulate ghost dofs
1107  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1108  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1109 
1110  // Calculate norm of vector and print values
1111  PetscReal nrm_F;
1112  CHKERR VecNorm(F, NORM_2, &nrm_F);
1113 
1114  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1115  PetscReal nrm_F_thermal;
1116  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1117  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1118  nrm_F_thermal);
1119 
1120  CHKERR VecScale(F_thermal, -1);
1121  // check this !!!
1122  CHKERR VecAXPY(F_thermal, 1, F);
1123 
1124  // Set dirichlet boundary to thermal stresses vector
1125  dirichlet_bc_ptr->snes_x = D;
1126  dirichlet_bc_ptr->snes_f = F_thermal;
1127  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1128 
1129  // Solve problem
1130  CHKERR KSPSolve(solver, F_thermal, D);
1131  // Add boundary conditions vector
1132  CHKERR VecAXPY(D, 1., D0);
1133  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1134  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1135 
1136  // Save data on the mesh
1137  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1138 
1139  // Save data on mesh
1140  CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1141 
1142  if (field_eval_flag) {
1143  CHKERR m_field.getInterface<FieldEvaluatorInterface>()
1144  ->evalFEAtThePoint3D(
1145  field_eval_coords.data(), 1e-12, "ELASTIC_PROB", "ELASTIC",
1146  field_eval_data, m_field.get_comm_rank(),
1147  m_field.get_comm_rank(), nullptr, MF_EXIST, QUIET);
1148  if (scalar_field_ptr->size()) {
1149  auto t_temp = getFTensor0FromVec(*scalar_field_ptr);
1150  MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
1151  << "Eval point TEMP: " << t_temp;
1152  }
1153  if (vector_field_ptr->size1()) {
1154  FTensor::Index<'i', 3> i;
1155  auto t_disp = getFTensor1FromMat<3>(*vector_field_ptr);
1156  MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
1157  << "Eval point DISPLACEMENT magnitude: "
1158  << sqrt(t_disp(i) * t_disp(i));
1159  }
1160  if (tensor_field_ptr->size1()) {
1161  FTensor::Index<'i', 3> i;
1162  auto t_disp_grad = getFTensor2FromMat<3, 3>(*tensor_field_ptr);
1163  MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
1164  << "Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(i, i);
1165  }
1166 
1167  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
1168  }
1169 
1170  // Post-process results
1171  if (is_post_proc_volume == PETSC_TRUE) {
1172  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1173  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1174  std::ostringstream o1;
1175  o1 << "out_" << sit->step_number << ".h5m";
1176  if (!test_nb)
1177  CHKERR post_proc.writeFile(o1.str().c_str());
1178  MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
1179  }
1180 
1181  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1182  CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN",
1183  &post_proc_skin);
1184  std::ostringstream o1_skin;
1185  o1_skin << "out_skin_" << sit->step_number << ".h5m";
1186  if (!test_nb)
1187  CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1188  MOFEM_LOG("ELASTIC", Sev::inform) << "done ...";
1189  }
1190  } else {
1191 
1192  // This is a case when stationary problem for temperature was solved.
1193  CHKERR VecZeroEntries(F_thermal);
1194  CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1195  CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1196 
1197  // Calculate the right-hand side vector with thermal stresses
1198  CHKERR DMoFEMLoopFiniteElements(
1199  dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1200 
1201  // Assemble vector
1202  CHKERR VecAssemblyBegin(F_thermal);
1203  CHKERR VecAssemblyEnd(F_thermal);
1204 
1205  // Accumulate ghost dofs
1206  CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1207  CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1208 
1209  // Calculate norm
1210  PetscReal nrm_F;
1211  CHKERR VecNorm(F, NORM_2, &nrm_F);
1212 
1213  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F = %6.4e", nrm_F);
1214  PetscReal nrm_F_thermal;
1215  CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1216 
1217  MOFEM_LOG_C("ELASTIC", Sev::inform, "norm2 F_thermal = %6.4e",
1218  nrm_F_thermal);
1219 
1220  // Add thermal stress vector and other forces vector
1221  CHKERR VecScale(F_thermal, -1);
1222  CHKERR VecAXPY(F_thermal, 1, F);
1223 
1224  // Apply kinetic boundary conditions
1225  dirichlet_bc_ptr->snes_x = D;
1226  dirichlet_bc_ptr->snes_f = F_thermal;
1227  CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1228 
1229  // Solve problem
1230  CHKERR KSPSolve(solver, F_thermal, D);
1231  CHKERR VecAXPY(D, 1., D0);
1232 
1233  // Update ghost values for solution vector
1234  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1235  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1236  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1237 
1238  // Save data on mesh
1239  if (is_post_proc_volume == PETSC_TRUE) {
1240  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1241  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1242  // Save results to file
1243  if (!test_nb)
1244  CHKERR post_proc.writeFile("out.h5m");
1245  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1246  }
1247 
1248  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
1249  CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1250  if (!test_nb)
1251  CHKERR post_proc_skin.writeFile("out_skin.h5m");
1252  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1253  }
1254 
1255  // Destroy vector, no needed any more
1256  CHKERR VecDestroy(&F_thermal);
1257  } else {
1258  // Elastic analysis no temperature field
1259  // VecView(F, PETSC_VIEWER_STDOUT_WORLD);
1260  // Solve for vector D
1261  CHKERR KSPSolve(solver, F, D);
1262 
1263  // VecView(D, PETSC_VIEWER_STDOUT_WORLD);
1264  // cerr << F;
1265 
1266  // Add kinetic boundary conditions
1267  CHKERR VecAXPY(D, 1., D0);
1268  // Update ghost values
1269  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1270  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1271  // Save data from vector on mesh
1272  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1273  // Post-process results
1274  MOFEM_LOG("ELASTIC", Sev::inform) << "Post-process start ...";
1275  if (is_post_proc_volume == PETSC_TRUE) {
1276  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
1277  }
1278  CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
1279  CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
1280  CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
1281  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1282 
1283  // Write mesh in parallel (using h5m MOAB format, writing is in parallel)
1284  MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
1285  if (mesh_has_tets) {
1286  if (is_post_proc_volume == PETSC_TRUE) {
1287  if (!test_nb)
1288  CHKERR post_proc.writeFile("out.h5m");
1289  }
1290  if (!test_nb)
1291  CHKERR post_proc_skin.writeFile("out_skin.h5m");
1292  }
1293  if (mesh_has_prisms) {
1294  if (!test_nb)
1295  CHKERR prism_post_proc.writeFile("prism_out.h5m");
1296  }
1297  if (!edges_in_simple_rod.empty())
1298  if (!test_nb)
1299  CHKERR post_proc_edge.writeFile("out_edge.h5m");
1300  MOFEM_LOG("ELASTIC", Sev::inform) << "done";
1301  }
1302 
1303  if (is_calculating_frequency == PETSC_TRUE) {
1304  // Calculate mode mass, m = u^T * M * u
1305  Vec u1;
1306  VecDuplicate(D, &u1);
1307  CHKERR MatMult(Mij, D, u1);
1308  double mode_mass;
1309  CHKERR VecDot(u1, D, &mode_mass);
1310  MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode mass %6.4e\n", mode_mass);
1311 
1312  Vec v1;
1313  VecDuplicate(D, &v1);
1314  CHKERR MatMult(Aij, D, v1);
1315 
1316  double mode_stiffness;
1317  CHKERR VecDot(v1, D, &mode_stiffness);
1318  MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode stiffness %6.4e\n",
1319  mode_stiffness);
1320 
1321  double frequency;
1322  double pi = 3.14159265359;
1323  frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1324  MOFEM_LOG_C("ELASTIC", Sev::inform, "Frequency %6.4e", frequency);
1325  }
1326 
1327  // Calculate elastic energy
1328  auto calculate_strain_energy = [&]() {
1329  MoFEMFunctionBegin;
1330 
1331  SmartPetscObj<Vec> v_energy;
1332  CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
1333  "MESH_NODE_POSITIONS", false, true,
1334  v_energy);
1335 
1336  // Print elastic energy
1337  double energy;
1338  CHKERR VecSum(v_energy, &energy);
1339  MOFEM_LOG_C("ELASTIC", Sev::inform, "Elastic energy %6.4e", energy);
1340 
1341  switch (test_nb) {
1342  case 1:
1343  if (fabs(energy - 17.129) > 1e-3)
1344  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1345  "atom test diverged!");
1346  break;
1347  case 2:
1348  if (fabs(energy - 5.6475e-03) > 1e-4)
1349  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1350  "atom test diverged!");
1351  break;
1352  case 3:
1353  if (fabs(energy - 7.4679e-03) > 1e-4)
1354  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1355  "atom test diverged!");
1356  break;
1357  case 4:
1358  if (fabs(energy - 2.4992e+00) > 1e-3)
1359  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1360  "atom test diverged!");
1361  break;
1362  // FIXME: Here are missing regersion tests
1363  case 8: {
1364  double min;
1365  CHKERR VecMin(D, PETSC_NULL, &min);
1366  constexpr double expected_val = 0.10001;
1367  if (fabs(min + expected_val) > 1e-10)
1368  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1369  "atom test diverged! %3.4e != %3.4e", min, expected_val);
1370  } break;
1371  case 9: {
1372  if (fabs(energy - 4.7416e-04) > 1e-8)
1373  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1374  "atom test diverged!");
1375  }
1376  default:
1377  break;
1378  }
1379 
1380  MoFEMFunctionReturn(0);
1381  };
1382  CHKERR calculate_strain_energy();
1383 
1384  MPI_Comm_free(&moab_comm_world);
1385  }
1386  CATCH_ERRORS;
1387 
1388  MoFEM::Core::Finalize();
1389 
1390  return 0;
1391 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:789
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
FluidPressure::addNeumannFluidPressureBCElements
MoFEMErrorCode addNeumannFluidPressureBCElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: FluidPressure.cpp:67
MetaNeumannForces::addNeumannBCElements
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.
Definition: SurfacePressure.cpp:1974
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
BodyForceConstantField::addBlock
MoFEMErrorCode addBlock(const std::string field_name, Vec F, int ms_id)
Definition: BodyForce.hpp:94
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MetaNodalForces::addElement
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:92
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
MoFEM::SeriesRecorder::check_series
virtual bool check_series(const std::string &name) const
check if series is in database
Definition: SeriesRecorder.cpp:301
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::SeriesRecorder
Definition: SeriesRecorder.hpp:25
PrismFE::getRuleTrianglesOnly
int getRuleTrianglesOnly(int order)
Definition: prism_elements_from_surface.cpp:436
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
PCMGSetUpViaApproxOrdersCtx
Set data structures of MG pre-conditioner via approximation orders.
Definition: PCMGSetUpViaApproxOrders.hpp:190
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ThermalStressElement::addThermalStressElement
MoFEMErrorCode addThermalStressElement(const std::string fe_name, const std::string field_name, const std::string thermal_field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: ThermalStressElement.hpp:155
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
ThermalStressElement::setThermalStressRhsOperators
MoFEMErrorCode setThermalStressRhsOperators(string field_name, string thermal_field_name, Vec &F, int verb=0)
Definition: ThermalStressElement.hpp:207
VolRule
Set integration rule.
Definition: simple_interface.cpp:88
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
BlockOptionData::tRis
Range tRis
Definition: elasticity.cpp:69
PostProcFatPrismOnRefinedMesh
Postprocess on prism.
Definition: PostProcOnRefMesh.hpp:973
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:257
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
BodyForceConstantField::getLoopFe
MyVolumeFE & getLoopFe()
Definition: BodyForce.hpp:23
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
BlockOptionData::initTemp
double initTemp
Definition: elasticity.cpp:67
ClipperLib::pi
static const double pi
Definition: clipper.cpp:89
calculateEnergy
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, SmartPetscObj< Vec > &v_energy_ptr)
MoFEM::CoreInterface::add_ents_to_field_by_dim
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.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3150
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MetaNodalForces::setOperators
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:128
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1171
MoFEM::FatPrismElementForcesAndSourcesCore::FatPrismElementForcesAndSourcesCore
FatPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: FatPrismElementForcesAndSourcesCore.cpp:11
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1201
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:301
PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:728
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MetaSimpleRodElement::setSimpleRodOperators
static MoFEMErrorCode setSimpleRodOperators(MoFEM::Interface &m_field, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_lhs_ptr, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of SimpleRod element. Set operators to calculate LHS and RHS.
Definition: SimpleRodElement.cpp:298
BodyForceConstantField
Body forces elements.
Definition: BodyForce.hpp:12
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
PrismFE::getRuleThroughThickness
int getRuleThroughThickness(int order)
Definition: prism_elements_from_surface.cpp:437
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3177
MetaNeumannForces::setMomentumFluxOperators
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.
Definition: SurfacePressure.cpp:2069
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:137
EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: elasticity.cpp:94
VolRule::operator()
int operator()(int, int, int order) const
Definition: elasticity.cpp:80
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
Hooke.hpp
Implementation of linear elastic material.
BODYFORCESSET
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
MoFEM::CoreInterface::add_field
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.
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
FluidPressure
Fluid pressure forces.
Definition: FluidPressure.hpp:20
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcOnRefMesh.hpp:253
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
Definition: LoopMethods.hpp:145
MoFEM::SeriesRecorder::load_series_data
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
Definition: SeriesRecorder.cpp:309
_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
Definition: SeriesRecorder.hpp:205
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition: ConvectiveMassElement.hpp:116
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:550
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
PrismFE
Definition: prism_elements_from_surface.cpp:62
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:560
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
PostProcEdgeOnRefinedMesh
Postprocess on edge.
Definition: PostProcOnRefMesh.hpp:1148
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:230
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::FatPrismElementForcesAndSourcesCore
FatPrism finite element.
Definition: FatPrismElementForcesAndSourcesCore.hpp:31
MoFEM::DMMoFEMCreateMoFEM
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: DMMoFEM.cpp:118
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:322
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
BlockOptionData::yOung
double yOung
Definition: elasticity.cpp:65
FTensor::Index< 'i', 3 >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::CoreInterface::get_finite_element_entities_by_dimension
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
MetaEdgeForces::addElement
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:62
Range
ConvectiveMassElement::setBlocks
MoFEMErrorCode setBlocks()
Definition: ConvectiveMassElement.cpp:1761
DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:376
MetaSpringBC::setSpringOperators
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", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
Definition: SpringElement.cpp:1178
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FluidPressure::getLoopFe
MyTriangleFE & getLoopFe()
Definition: FluidPressure.hpp:35
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:287
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
FluidPressure::setNeumannFluidPressureFiniteElementOperators
MoFEMErrorCode setNeumannFluidPressureFiniteElementOperators(string field_name, Vec F, bool allow_negative_pressure=true, bool ho_geometry=false)
Definition: FluidPressure.cpp:131
MetaSimpleRodElement::addSimpleRodElements
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
Definition: SimpleRodElement.cpp:270
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
BlockOptionData::oRder
int oRder
Definition: elasticity.cpp:63
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BlockData
Definition: MeshsetsManager.cpp:752
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
addElasticElement
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)
Definition: HookeElement.cpp:533
main
int main(int argc, char *argv[])
Definition: elasticity.cpp:96
setOperators
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, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
BlockOptionData
Definition: elasticity.cpp:62
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
BlockOptionData::BlockOptionData
BlockOptionData()
Definition: elasticity.cpp:71
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
ThermalStressElement::getLoopThermalStressRhs
MyVolumeFE & getLoopThermalStressRhs()
Definition: ThermalStressElement.hpp:24
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::CoreInterface::set_field_order
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.
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Vec >
help
static char help[]
Definition: elasticity.cpp:57
BlockOptionData::iD
int iD
Definition: elasticity.cpp:64
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MF_EXIST
@ MF_EXIST
Definition: definitions.h:100
convert.int
int
Definition: convert.py:64
BlockOptionData::pOisson
double pOisson
Definition: elasticity.cpp:66
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
ThermalStressElement
Implentation of thermal stress element.
Definition: ThermalStressElement.hpp:15
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
BlockData
NonlinearElasticElement::BlockData BlockData
Definition: elasticity.cpp:74
F
@ F
Definition: free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
MetaEdgeForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
Generated by Doxygen 1.8.17 and hosted at University of Glasgow