1 #ifndef KATOPTRON_RESULTSDEF_HPP
2 #define KATOPTRON_RESULTSDEF_HPP
8 #ifdef WAVES_USE_TBOXVTK
13 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal>
16 Teuchos::RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal>> xPost,
17 Teuchos::RCP<katoptron::Map> map,
25 for (
auto e : pbl.
msh->elems)
27 if (e->type() == tbox::ELTYPE::HEX8 || e->type() == tbox::ELTYPE::TETRA4)
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);
42 const size_t ensemble_size = ET::size;
44 typedef decltype(X) X_View;
45 typename Kokkos::FlatArrayType<X_View>::type X_flat = X;
47 size_t numPrimalDPN = map->numPrimalDPN;
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;
58 LocalOrdinal index = 0;
60 for (LocalOrdinal i = 0; i < map->mapElems->getNodeNumElements(); ++i)
62 if (pbl.
msh->elems[i]->type() == tbox::ELTYPE::HEX8 || pbl.
msh->elems[i]->type() == tbox::ELTYPE::TETRA4)
64 GlobalOrdinal global_i = map->mapElems->getGlobalElement(i);
65 elements_GO_to_LO[global_i] = index;
70 for (LocalOrdinal i = 0; i < map->mapNodesWO->getNodeNumElements(); ++i)
72 GlobalOrdinal global_i = map->mapNodesWO->getGlobalElement(i);
73 nodesWO_LO_to_GO[i] = global_i;
74 nodesWO_GO_to_LO[global_i] = i;
77 LocalOrdinal numNodes = nodesWO_LO_to_GO.
size();
78 LocalOrdinal numMechanicalDualDOFs = map->lagrange_to_dof.size();
79 LocalOrdinal numPrimalDOFs = numPrimalDPN * numNodes;
81 std::map<LocalOrdinal, GlobalOrdinal> elems_LO_to_GO;
82 std::map<GlobalOrdinal, LocalOrdinal> elems_GO_to_LO;
84 for (GlobalOrdinal i = 0; i < pbl.
msh->elems.size(); ++i)
86 if (pbl.
msh->elems[i]->type() == tbox::ELTYPE::HEX8 || pbl.
msh->elems[i]->type() == tbox::ELTYPE::TETRA4)
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;
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);
101 if (tmpcontactElementID.size() > 1)
103 std::sort(tmpcontactElementID.begin(), tmpcontactElementID.end());
105 if (tmpcontactElementID.size() >= 1)
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]);
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);
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);
128 std::vector<Eigen::MatrixXd> stress(pbl.
msh->elems.size());
129 std::vector<Eigen::MatrixXd> strain(pbl.
msh->elems.size());
131 std::vector<double> sigma_vm(pbl.
msh->elems.size());
133 std::vector<double> G_vector(pbl.
msh->elems.size());
135 std::vector<Eigen::Vector3d> displacement(numNodes);
137 std::vector<double> displacement_u(numNodes);
138 std::vector<double> T(numNodes);
139 std::vector<double> pressure(numMechanicalDualDOFs);
141 std::vector<double> h_fluxes(numMechanicalDualDOFs);
142 std::vector<double> t_x(numMechanicalDualDOFs);
143 std::vector<double> t_y(numMechanicalDualDOFs);
145 for (
auto &k : stress)
146 k = Eigen::Matrix3d::Zero();
148 for (
auto &k : strain)
149 k = Eigen::Matrix3d::Zero();
151 for (
auto &k : stress_at_nodes)
152 k = Eigen::Matrix3d::Zero();
154 for (
auto &k : strain_at_nodes)
155 k = Eigen::Matrix3d::Zero();
159 for (
size_t ell = 0; ell < ensemble_size; ++ell)
162 tbox::Results *results =
new tbox::Results();
163 tbox::Results *results_contact =
new tbox::Results();
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;
174 if (numPrimalDPN == 4)
175 results->scalars_at_nodes[
"T"] = &T;
177 results->tensors_at_elems[
"sigma"] = &stress;
178 results->tensors_at_elems[
"epsilon"] = &strain;
180 results->scalars_at_elems[
"von Mises"] = &sigma_vm;
182 results->scalars_at_elems[
"G"] = &G_vector;
184 results->vectors_at_nodes[
"displacement"] = &displacement;
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)
193 results_contact->scalars_at_nodes[
"tangential stress 1"] = &t_x;
194 results_contact->scalars_at_nodes[
"tangential stress 2"] = &t_y;
198 for (LocalOrdinal i = 0; i < numNodes; ++i)
201 if (numPrimalDPN == 4)
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];
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];
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];
227 for (LocalOrdinal i = 0; i < numMechanicalDualDOFs; ++i)
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];
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)
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];
245 for (
auto m : pbl.
media)
247 tbox::Tag *current_tag = m->tag;
249 double E = ET::coeff(m->E, ell);
250 double nu = ET::coeff(m->nu, ell);
252 double lambda = (nu * E) / ((1 + nu) * (1 - 2 * nu));
253 double G = E / (2 * (1 + nu));
255 for (LocalOrdinal i_elem = 0; i_elem < current_tag->elems.size(); ++i_elem)
257 if (random_field->isRandom)
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)
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;
267 G = ET::coeff(random_field->operator()(x, y, z), ell);
269 Eigen::MatrixXd H = Eigen::MatrixXd::Zero(9, 9);
270 Eigen::MatrixXd D = Eigen::MatrixXd::Zero(3, 3);
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)
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;
284 if (numPrimalDPN == 4)
286 double kappa = lambda + (2. / 3.) * G;
287 double beta = ET::coeff(m->beta, ell);
289 for (
auto i = 0; i < 3; ++i)
290 D(i, i) = -3. * kappa * beta;
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);
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());
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);
305 if (e->type() == tbox::ELTYPE::HEX8 || e->type() == tbox::ELTYPE::TETRA4)
307 for (
size_t i = 0; i < e->nodes.size(); ++i)
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];
319 strain_stress_x_y_z(e, Strain, Sigma, H, D, x_tmp, y_tmp, z_tmp, T_tmp, numPrimalDPN == 4);
321 strain[elems_LO_to_GO[e->no - 1]] = Strain;
322 stress[elems_LO_to_GO[e->no - 1]] = Sigma;
324 G_vector[elems_LO_to_GO[e->no - 1]] = G;
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);
330 for (
size_t i = 0; i < e->nodes.size(); ++i)
332 LocalOrdinal node_i = nodesWO_GO_to_LO[e->nodes[i]->row];
334 ++adjacent_elems[node_i];
336 strain_at_nodes[node_i] += Strain_at_node[i];
337 stress_at_nodes[node_i] += Sigma_at_node[i];
342 for (LocalOrdinal i = 0; i < numNodes; ++i)
344 strain_at_nodes[i] /= adjacent_elems[i];
345 stress_at_nodes[i] /= adjacent_elems[i];
347 Eigen::MatrixXd Sigma = stress_at_nodes[i];
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);
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);
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);
362 delete results_contact;
366 #endif //KATOPTRON_RESULTSDEF_HPP