waves
Basic FE playground
ResultsDef.hpp
Go to the documentation of this file.
1 #ifndef KATOPTRON_RESULTSDEF_HPP
2 #define KATOPTRON_RESULTSDEF_HPP
3 
4 #include "ResultsDecl.hpp"
5 #include "StressComputation.h"
6 
7 /* @todo should be removed */
8 #ifdef WAVES_USE_TBOXVTK
9 #include "tboxVtk.h"
10 #include "wVtkExport_KIM2CLEAN.h"
11 #endif
12 
13 template <class Scalar, class LocalOrdinal, class GlobalOrdinal>
15  const katoptron::Problem &pbl,
16  Teuchos::RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal>> xPost,
17  Teuchos::RCP<katoptron::Map> map,
18  Teuchos::RCP<RandomField<Scalar, Kokkos::DefaultExecutionSpace>> random_field,
19  std::string name,
20  size_t ast,
21  bool final_result,
22  bool write_txt)
23 {
24  /* @todo [AC] temporary fix for initializing Mem in katoptron */
25  for (auto e : pbl.msh->elems)
26  {
27  if (e->type() == tbox::ElType::HEX8 || e->type() == tbox::ElType::TETRA4) // 3D
28  e->initValues(true);
29  //else // 2D
30  // e->initValues(false);
31  }
32  /* end temporary fix*/
33 
34  // Extraction of types and ensemble format
35 
36  xPost->template sync<Kokkos::HostSpace>();
37  auto x_2d = xPost->template getLocalView<Kokkos::HostSpace>();
38  auto X = Kokkos::subview(x_2d, Kokkos::ALL(), 0);
39 
40  typedef EnsembleTraits<Scalar> ET;
41 
42  const size_t ensemble_size = ET::size;
43 
44  typedef decltype(X) X_View;
45  typename Kokkos::FlatArrayType<X_View>::type X_flat = X;
46 
47  size_t numPrimalDPN = map->numPrimalDPN;
48 
49  // Construction of the maps
50 
51  std::map<LocalOrdinal, GlobalOrdinal> nodesWO_LO_to_GO;
52  std::map<GlobalOrdinal, LocalOrdinal> nodesWO_GO_to_LO;
53  std::map<GlobalOrdinal, LocalOrdinal> elements_GO_to_LO;
54  std::vector<GlobalOrdinal> contactElementID;
55  std::map<LocalOrdinal, GlobalOrdinal> lagrange_LO_to_nodes_GO;
56  std::map<GlobalOrdinal, LocalOrdinal> nodes_GO_to_lagrange_LO;
57 
58  LocalOrdinal index = 0;
59 
60  for (LocalOrdinal i = 0; i < map->mapElems->getNodeNumElements(); ++i)
61  {
62  if (pbl.msh->elems[i]->type() == tbox::ElType::HEX8 || pbl.msh->elems[i]->type() == tbox::ElType::TETRA4)
63  {
64  GlobalOrdinal global_i = map->mapElems->getGlobalElement(i);
65  elements_GO_to_LO[global_i] = index;
66  ++index;
67  }
68  }
69 
70  for (LocalOrdinal i = 0; i < map->mapNodesWO->getNodeNumElements(); ++i)
71  {
72  GlobalOrdinal global_i = map->mapNodesWO->getGlobalElement(i);
73  nodesWO_LO_to_GO[i] = global_i;
74  nodesWO_GO_to_LO[global_i] = i;
75  }
76 
77  LocalOrdinal numNodes = nodesWO_LO_to_GO.size();
78  LocalOrdinal numMechanicalDualDOFs = map->lagrange_to_dof.size();
79  LocalOrdinal numPrimalDOFs = numPrimalDPN * numNodes;
80 
81  std::map<LocalOrdinal, GlobalOrdinal> elems_LO_to_GO;
82  std::map<GlobalOrdinal, LocalOrdinal> elems_GO_to_LO;
83 
84  for (GlobalOrdinal i = 0; i < pbl.msh->elems.size(); ++i)
85  {
86  if (pbl.msh->elems[i]->type() == tbox::ElType::HEX8 || pbl.msh->elems[i]->type() == tbox::ElType::TETRA4)
87  {
88  LocalOrdinal local_i = pbl.msh->elems[i]->no - 1;
89  elems_LO_to_GO[local_i] = i;
90  elems_GO_to_LO[i] = local_i;
91  }
92  }
93 
94  if (pbl.Contacts.size() >= 1)
95  {
96  std::vector<GlobalOrdinal> tmpcontactElementID;
97  for (size_t i = 0; i < pbl.Contacts.size(); ++i)
98  for (LocalOrdinal j = 0; j < pbl.Contacts[i]->slave_elems.size(); ++j)
99  tmpcontactElementID.push_back(pbl.Contacts[i]->slave_elems[j]->no - 1);
100 
101  if (tmpcontactElementID.size() > 1)
102  {
103  std::sort(tmpcontactElementID.begin(), tmpcontactElementID.end());
104 
105  if (tmpcontactElementID.size() >= 1)
106  {
107  contactElementID.push_back(tmpcontactElementID[0]);
108  if (tmpcontactElementID.size() > 1)
109  for (LocalOrdinal i = 1; i < tmpcontactElementID.size(); ++i)
110  if (tmpcontactElementID[i] != tmpcontactElementID[i - 1])
111  contactElementID.push_back(tmpcontactElementID[i]);
112  }
113  }
114  }
115 
116  // Values in which we are interested
117 
118  std::vector<double> displacement_x(numNodes);
119  std::vector<double> displacement_y(numNodes);
120  std::vector<double> displacement_z(numNodes);
121  std::vector<double> shear_modulus_per_node(numNodes);
122 
123  std::vector<double> sigma_vm_nodes(numNodes);
124  std::vector<size_t> adjacent_elems(numNodes);
125  std::vector<Eigen::MatrixXd> stress_at_nodes(numNodes);
126  std::vector<Eigen::MatrixXd> strain_at_nodes(numNodes);
127 
128  std::vector<Eigen::MatrixXd> stress(pbl.msh->elems.size());
129  std::vector<Eigen::MatrixXd> strain(pbl.msh->elems.size());
130 
131  std::vector<double> sigma_vm(pbl.msh->elems.size());
132 
133  std::vector<double> G_vector(pbl.msh->elems.size());
134 
135  std::vector<Eigen::Vector3d> displacement(numNodes);
136 
137  std::vector<double> displacement_u(numNodes);
138  std::vector<double> T(numNodes);
139  std::vector<double> pressure(numMechanicalDualDOFs);
140 
141  std::vector<double> h_fluxes(numMechanicalDualDOFs);
142  std::vector<double> t_x(numMechanicalDualDOFs);
143  std::vector<double> t_y(numMechanicalDualDOFs);
144 
145  for (auto &k : stress)
146  k = Eigen::Matrix3d::Zero();
147 
148  for (auto &k : strain)
149  k = Eigen::Matrix3d::Zero();
150 
151  for (auto &k : stress_at_nodes)
152  k = Eigen::Matrix3d::Zero();
153 
154  for (auto &k : strain_at_nodes)
155  k = Eigen::Matrix3d::Zero();
156 
157  // Fill the result vectors
158 
159  for (size_t ell = 0; ell < ensemble_size; ++ell)
160  {
161 
162  tbox::Results *results = new tbox::Results();
163  tbox::Results *results_contact = new tbox::Results();
164 
165  results->scalars_at_nodes["x"] = &displacement_x;
166  results->scalars_at_nodes["y"] = &displacement_y;
167  results->scalars_at_nodes["z"] = &displacement_z;
168  results->scalars_at_nodes["u"] = &displacement_u;
169  results->scalars_at_nodes["G"] = &shear_modulus_per_node;
170  results->scalars_at_nodes["von Mises"] = &sigma_vm_nodes;
171  results->tensors_at_nodes["sigma"] = &stress_at_nodes;
172  results->tensors_at_nodes["epsilon"] = &strain_at_nodes;
173 
174  if (numPrimalDPN == 4)
175  results->scalars_at_nodes["T"] = &T;
176 
177  results->tensors_at_elems["sigma"] = &stress;
178  results->tensors_at_elems["epsilon"] = &strain;
179 
180  results->scalars_at_elems["von Mises"] = &sigma_vm;
181 
182  results->scalars_at_elems["G"] = &G_vector;
183 
184  results->vectors_at_nodes["displacement"] = &displacement;
185 
186  if (pbl.Contacts.size() > 0)
187  {
188  results_contact->scalars_at_nodes["p"] = &pressure;
189  if (numPrimalDPN == 4)
190  results_contact->scalars_at_nodes["heat fluxes"] = &h_fluxes;
191  if (map->numDualDPN == 3 || map->numDualDPN == 4)
192  {
193  results_contact->scalars_at_nodes["tangential stress 1"] = &t_x;
194  results_contact->scalars_at_nodes["tangential stress 2"] = &t_y;
195  }
196  }
197 
198  for (LocalOrdinal i = 0; i < numNodes; ++i)
199  {
200  double dx, dy, dz;
201  if (numPrimalDPN == 4)
202  {
203  dx = X_flat[ensemble_size * (numNodes + 3 * i) + ell];
204  dy = X_flat[ensemble_size * (numNodes + 3 * i + 1) + ell];
205  dz = X_flat[ensemble_size * (numNodes + 3 * i + 2) + ell];
206  }
207  else
208  {
209  dx = X_flat[ensemble_size * (3 * i) + ell];
210  dy = X_flat[ensemble_size * (3 * i + 1) + ell];
211  dz = X_flat[ensemble_size * (3 * i + 2) + ell];
212  }
213 
214  displacement_x[i] = dx;
215  displacement_y[i] = dy;
216  displacement_z[i] = dz;
217  displacement_u[i] = sqrt(dx * dx + dy * dy + dz * dz);
218  strain_at_nodes[i].setZero();
219  stress_at_nodes[i].setZero();
220  shear_modulus_per_node[i] = ET::coeff(random_field->operator()(pbl.msh->nodes[nodesWO_LO_to_GO[i]]->pos(0), pbl.msh->nodes[nodesWO_LO_to_GO[i]]->pos(1), pbl.msh->nodes[nodesWO_LO_to_GO[i]]->pos(2)), ell);
221  displacement[i] = Eigen::Vector3d(dx, dy, dz);
222  if (numPrimalDPN == 4 || numPrimalDPN == 1)
223  T[i] = X_flat[ensemble_size * i + ell];
224  }
225  if (pbl.Contacts.size() > 0)
226  {
227  for (LocalOrdinal i = 0; i < numMechanicalDualDOFs; ++i)
228  {
229  int i_node = map->lm_to_dof_global[i];
230  lagrange_LO_to_nodes_GO[i] = i_node;
231  nodes_GO_to_lagrange_LO[i_node] = i;
232  pressure[i] = X_flat[ensemble_size * (i * map->numDualDPN + numPrimalDOFs) + ell];
233  }
234  if (map->numDualDPN == 2 || map->numDualDPN == 4)
235  for (LocalOrdinal i = 0; i < numMechanicalDualDOFs; ++i)
236  h_fluxes[i] = X_flat[ensemble_size * (i * map->numDualDPN + map->numDualDPN - 1 + numPrimalDOFs) + ell];
237  if (map->numDualDPN == 3 || map->numDualDPN == 4)
238  for (LocalOrdinal i = 0; i < numMechanicalDualDOFs; ++i)
239  {
240  t_x[i] = X_flat[ensemble_size * (i * map->numDualDPN + 1 + numPrimalDOFs) + ell];
241  t_y[i] = X_flat[ensemble_size * (i * map->numDualDPN + 2 + numPrimalDOFs) + ell];
242  }
243  }
244 
245  for (auto m : pbl.media)
246  {
247  tbox::Tag *current_tag = m->tag;
248 
249  double E = ET::coeff(m->E, ell);
250  double nu = ET::coeff(m->nu, ell);
251 
252  double lambda = (nu * E) / ((1 + nu) * (1 - 2 * nu));
253  double G = E / (2 * (1 + nu));
254 
255  for (LocalOrdinal i_elem = 0; i_elem < current_tag->elems.size(); ++i_elem)
256  {
257  if (random_field->isRandom)
258  {
259  double x = 0., y = 0., z = 0.;
260  size_t element_size = current_tag->elems[i_elem]->nodes.size();
261  for (size_t j = 0; j < element_size; ++j)
262  {
263  x += current_tag->elems[i_elem]->nodes[j]->pos(0) / element_size;
264  y += current_tag->elems[i_elem]->nodes[j]->pos(1) / element_size;
265  z += current_tag->elems[i_elem]->nodes[j]->pos(2) / element_size;
266  }
267  G = ET::coeff(random_field->operator()(x, y, z), ell);
268  }
269  Eigen::MatrixXd H = Eigen::MatrixXd::Zero(9, 9);
270  Eigen::MatrixXd D = Eigen::MatrixXd::Zero(3, 3);
271 
272  for (auto i = 0; i < 3; ++i)
273  for (auto j = 0; j < 3; ++j)
274  for (auto k = 0; k < 3; ++k)
275  for (auto l = 0; l < 3; ++l)
276  {
277  if (i == j && k == l)
278  H(i * 3 + j, k * 3 + l) = lambda;
279  if (i == k && j == l)
280  H(i * 3 + j, k * 3 + l) = H(i * 3 + j, k * 3 + l) + G;
281  if (i == l && j == k)
282  H(i * 3 + j, k * 3 + l) = H(i * 3 + j, k * 3 + l) + G;
283  }
284  if (numPrimalDPN == 4)
285  {
286  double kappa = lambda + (2. / 3.) * G;
287  double beta = ET::coeff(m->beta, ell);
288 
289  for (auto i = 0; i < 3; ++i)
290  D(i, i) = -3. * kappa * beta;
291  }
292  std::vector<double> x_tmp(8), y_tmp(8), z_tmp(8), T_tmp(8);
293  Eigen::MatrixXd Strain(3, 3);
294  Eigen::MatrixXd Sigma(3, 3);
295 
296  auto e = current_tag->elems[i_elem];
297  std::vector<Eigen::MatrixXd> Strain_at_node(e->nodes.size());
298  std::vector<Eigen::MatrixXd> Sigma_at_node(e->nodes.size());
299 
300  for (auto &k : Strain_at_node)
301  k = Eigen::MatrixXd::Zero(3, 3);
302  for (auto &k : Sigma_at_node)
303  k = Eigen::MatrixXd::Zero(3, 3);
304 
305  if (e->type() == tbox::ElType::HEX8 || e->type() == tbox::ElType::TETRA4)
306  {
307  for (size_t i = 0; i < e->nodes.size(); ++i)
308  {
309  int node_i = nodesWO_GO_to_LO[e->nodes[i]->row];
310  x_tmp[i] = displacement_x[node_i];
311  y_tmp[i] = displacement_y[node_i];
312  z_tmp[i] = displacement_z[node_i];
313  if (numPrimalDPN == 4)
314  T_tmp[i] = T[node_i];
315  else
316  T_tmp[i] = 0.;
317  }
318 
319  strain_stress_x_y_z(e, Strain, Sigma, H, D, x_tmp, y_tmp, z_tmp, T_tmp, numPrimalDPN == 4);
320 
321  strain[elems_LO_to_GO[e->no - 1]] = Strain;
322  stress[elems_LO_to_GO[e->no - 1]] = Sigma;
323 
324  G_vector[elems_LO_to_GO[e->no - 1]] = G;
325 
326  sigma_vm[elems_LO_to_GO[e->no - 1]] = sqrt((pow(Sigma(0, 0) - Sigma(1, 1), 2.) + pow(Sigma(1, 1) - Sigma(2, 2), 2.) + pow(Sigma(2, 2) - Sigma(0, 0), 2.) + 6 * (pow(Sigma(0, 1), 2.) + pow(Sigma(1, 2), 2.) + pow(Sigma(2, 0), 2.))) / 2);
327 
328  strain_stress_x_y_z_at_node(e, Strain_at_node, Sigma_at_node, H, D, x_tmp, y_tmp, z_tmp, T_tmp, numPrimalDPN == 4);
329 
330  for (size_t i = 0; i < e->nodes.size(); ++i)
331  {
332  LocalOrdinal node_i = nodesWO_GO_to_LO[e->nodes[i]->row];
333  if (ell == 0)
334  ++adjacent_elems[node_i]; // Increases the number of adjacent elements by one only for the first sample
335 
336  strain_at_nodes[node_i] += Strain_at_node[i];
337  stress_at_nodes[node_i] += Sigma_at_node[i];
338  }
339  }
340  }
341  }
342  for (LocalOrdinal i = 0; i < numNodes; ++i)
343  {
344  strain_at_nodes[i] /= adjacent_elems[i];
345  stress_at_nodes[i] /= adjacent_elems[i];
346 
347  Eigen::MatrixXd Sigma = stress_at_nodes[i];
348 
349  sigma_vm_nodes[i] = sqrt((pow(Sigma(0, 0) - Sigma(1, 1), 2.) + pow(Sigma(1, 1) - Sigma(2, 2), 2.) + pow(Sigma(2, 2) - Sigma(0, 0), 2.) + 6 * (pow(Sigma(0, 1), 2.) + pow(Sigma(1, 2), 2.) + pow(Sigma(2, 0), 2.))) / 2);
350  }
351 
352  /* @todo temporary fix to output Kim's data */
353  // All the post-processing related to katoptron should be done here,
354  // then tbox::MshExport should be used, either by passing tbox::GmshExport or tboxVtk::VtkExport
355 
356  tboxVtk::VtkExport_KIM2CLEAN vtkWriter(pbl.msh);
357  vtkWriter.save_MPI(pbl.msh->name + "_ast" + std::to_string(ast) + "_s" + std::to_string(ell) + ".pvtu", *results, nodesWO_LO_to_GO, nodesWO_GO_to_LO, elements_GO_to_LO);
358  if (pbl.Contacts.size() >= 1)
359  vtkWriter.save_MPI(pbl.msh->name + "_ast" + std::to_string(ast) + "_s" + std::to_string(ell) + "_contact.vtu", *results_contact, contactElementID, lagrange_LO_to_nodes_GO, nodes_GO_to_lagrange_LO);
360 
361  delete results;
362  delete results_contact;
363  }
364 }
365 
366 #endif //KATOPTRON_RESULTSDEF_HPP
katoptron::Problem::media
std::vector< Medium * > media
Materials.
Definition: wProblem.h:26
tboxVtk::VtkExport_KIM2CLEAN::save_MPI
void save_MPI(std::string const &fname, tbox::Results const &r, std::map< int, int > nodeWO_LO_to_GO, std::map< int, int > nodeWO_GO_to_LO, std::map< int, int > elements_GO_to_LO) const
writeResultsVTK
void writeResultsVTK(const katoptron::Problem &pbl, Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal >> xPost, Teuchos::RCP< katoptron::Map > map, Teuchos::RCP< RandomField< Scalar, Kokkos::DefaultExecutionSpace >> random_field, std::string name, size_t ast, bool final_result, bool write_txt)
Function to write results on disk using VTK formats.
Definition: ResultsDef.hpp:14
ResultsDecl.hpp
wVtkExport_KIM2CLEAN.h
RandomField
Definition: wRandomField.h:10
StressComputation.h
EnsembleTraits< Scalar >
tboxVtk::VtkExport_KIM2CLEAN
Save routines of Kim.
Definition: wVtkExport_KIM2CLEAN.h:45
strain_stress_x_y_z_at_node
void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector< Eigen::MatrixXd > &Strain_at_node, std::vector< Eigen::MatrixXd > &Stress_at_node, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &T, bool thermal=false)
Definition: StressComputation.h:81
katoptron::Problem
Class which is used to specify in Python the thermomechanical to solve.
Definition: wProblem.h:19
katoptron::Problem::Contacts
std::vector< Contact * > Contacts
Definition: wProblem.h:29
EnsembleTraits::size
static const int size
Definition: EnsembleTraits.h:10
katoptron::Problem::msh
std::shared_ptr< tbox::MshData > msh
Mesh structure.
Definition: wProblem.h:22
strain_stress_x_y_z
void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &T, bool thermal=false)
Definition: StressComputation.h:55