waves
Basic FE playground
ElementMatrices.hpp
Go to the documentation of this file.
1 //#define ElementMatrices_DEBUG
2 
16 template <typename scalar, int element_type, int element_size>
18  size_t _numPrimalDPN,
19  size_t pool_size) : domain(_domain), numPrimalDPN(_numPrimalDPN)
20 {
21 
22  K = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A", pool_size, 3 * element_size, 3 * element_size);
23  L = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A", pool_size, element_size, element_size);
24  S = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A", pool_size, 3 * element_size, element_size);
25 
26  material = Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A", pool_size, 3 * element_size);
27 }
28 
42 template <typename scalar, int element_type, int element_size>
44 {
45 
46  tbox::Cache *cache;
47 
48  /* @todo could probably be refactored into cache = element->getVCache(); */
49  if (element_type == static_cast<int>(tbox::ELTYPE::HEX8))
50  {
51  tbox::CacheHex8 &cachehex8 = this->trilinosHex8GetCache();
52  cache = &cachehex8;
53  }
54  else if (element_type == static_cast<int>(tbox::ELTYPE::TETRA4))
55  {
56  tbox::CacheTetra4 &cachetetra4 = this->trilinosTetra4GetCache();
57  cache = &cachetetra4;
58  }
59 
60  std::vector<Eigen::Vector3d> gauss_p;
61  std::vector<double> gauss_w;
62  for (size_t a = 0; a < cache->getVGauss().getN(); ++a)
63  {
64  gauss_p.push_back(cache->getVGauss().getP(a));
65  gauss_w.push_back(cache->getVGauss().getW(a));
66  }
67 
69 
70  auto K_current = subview(K, i_thread, Kokkos::ALL(), Kokkos::ALL());
71  auto L_current = subview(L, i_thread, Kokkos::ALL(), Kokkos::ALL());
72  auto S_current = subview(S, i_thread, Kokkos::ALL(), Kokkos::ALL());
73 
74  Kokkos::deep_copy(K_current, 0);
75  Kokkos::deep_copy(L_current, 0);
76  Kokkos::deep_copy(S_current, 0);
77 
78  auto element_material = domain->elementsList->getElementMaterial(e);
79 
80  tVector<double, 3> factor1, factor2, dffi, dffj;
81  tVector<scalar, 3> factor3;
82  tMatrix<double, 3, 3> J, Jinv;
86  scalar E, nu, lambda, G, conductivity, beta;
87  tMatrix<scalar, 3, 3> K_conductivity;
89 
90  katoptron::Material<scalar> material = domain->materialsList->getMaterial(element_material, i_thread);
91  if (domain->random_field->isRandom)
92  {
93  double x = 0., y = 0., z = 0.;
94  for (auto j = 0; j < element_size; ++j)
95  {
96  x += domain->nodesList->nodes(domain->elementsList->getElementNode(e, j), 0) / element_size;
97  y += domain->nodesList->nodes(domain->elementsList->getElementNode(e, j), 1) / element_size;
98  z += domain->nodesList->nodes(domain->elementsList->getElementNode(e, j), 2) / element_size;
99  }
100 
101  G = domain->random_field->operator()(x, y, z);
102 #ifdef ElementMatrices_DEBUG
103  printf("G_0 = %f\n", EnsembleTraits<scalar>::coeff(G, 0));
104 #endif
105 
106  //E = material.getE();
107  nu = material.getNu();
108  lambda = 2. * G * nu / (1. - 2. * nu);
109  //lambda = (nu*E)/((1+nu)*(1-2*nu));
110  }
111  else
112  {
113  E = material.getE();
114  nu = material.getNu();
115  lambda = (nu * E) / ((1 + nu) * (1 - 2 * nu));
116  G = E / (2 * (1 + nu));
117  }
118 
119  // Machine epsilon
120 
121  double epsilon = 2e-16;
122  // Looping on gauss points
123  for (auto a = 0; a < gauss_p.size(); ++a)
124  {
125  for (auto j = 0; j < element_size; ++j)
126  if (element_size == 8)
127  ff(j) = (static_cast<tbox::CacheHex8 *>(cache)->getSf(a))(j);
128  else if (element_size == 4)
129  ff(j) = (static_cast<tbox::CacheTetra4 *>(cache)->getSf(a))(j);
130 
131  // H matrix
132  H.clean();
133 
134  for (auto i = 0; i < 3; ++i)
135  for (auto j = 0; j < 3; ++j)
136  for (auto k = 0; k < 3; ++k)
137  for (auto l = 0; l < 3; ++l)
138  {
139  if (i == j && k == l)
140  H(i * 3 + j, k * 3 + l) = lambda;
141  if (i == k && j == l)
142  H(i * 3 + j, k * 3 + l) = H(i * 3 + j, k * 3 + l) + G;
143  if (i == l && j == k)
144  H(i * 3 + j, k * 3 + l) = H(i * 3 + j, k * 3 + l) + G;
145  }
146 
147  conductivity = material.getK();
148  K_conductivity.clean();
149  for (auto i = 0; i < 3; ++i)
150  K_conductivity(i, i) = conductivity;
151 
152  beta = material.getBeta();
153  D.clean();
154  for (auto i = 0; i < 3; ++i)
155  D(i, i) = -3 * (lambda + (2. / 3.) * G) * beta;
156 
157  double detJ = this->buildJ(a, J, e, *(domain->elementsList), *(domain->nodesList));
158 
159 #ifdef ElementMatrices_DEBUG
160  printf("detJ = %f\n", detJ);
161 #endif
162 
163  Jinv = J.inv();
164 
165  for (auto i = 0; i < element_size; ++i)
166  {
167 
168  for (auto k = 0; k < 3; ++k)
169  if (element_size == 8)
170  dffi(k) = (static_cast<tbox::CacheHex8 *>(cache)->getDsf(a))(k, i);
171  else if (element_size == 4)
172  dffi(k) = (static_cast<tbox::CacheTetra4 *>(cache)->getDsf(a))(k, i);
173  factor1 = Jinv * dffi;
174 #ifdef ElementMatrices_DEBUG
175  printf("dffi = [%f, %f, %f]\n", dffi(0), dffi(1), dffi(2));
176  printf("factor1 = [%f, %f, %f]\n", factor1(0), factor1(1), factor1(2));
177 #endif
178  // K
179  for (auto j = i; j < element_size; ++j)
180  {
181  for (auto k = 0; k < 3; ++k)
182  if (element_size == 8)
183  dffj(k) = (static_cast<tbox::CacheHex8 *>(cache)->getDsf(a))(k, j);
184  else if (element_size == 4)
185  dffj(k) = (static_cast<tbox::CacheTetra4 *>(cache)->getDsf(a))(k, j);
186  factor2 = Jinv * dffj;
187 
188  TMP1.clean();
189  for (auto k = 0; k < 3; ++k)
190  for (auto m = 0; m < 3; ++m)
191  {
192  for (auto l = 0; l < 3; ++l)
193  for (auto n = 0; n < 3; ++n)
194  TMP1(l, n) = H(k * 3 + l, m * 3 + n);
195 
196  factor3 = TMP1 * factor2;
197 
198  scalar dotProduct = factor1.dotproduct(factor3);
199 #ifdef ElementMatrices_DEBUG
200  printf("dotProduct_0 = %f\n", EnsembleTraits<scalar>::coeff(dotProduct, 0));
201 #endif
202  K_current(i * 3 + k, j * 3 + m) += dotProduct * gauss_w[a] * detJ;
203  }
204  }
205 
206  // L
207  factor3 = K_conductivity * factor1;
208  for (auto j = i; j < element_size; ++j)
209  {
210  for (auto k = 0; k < 3; ++k)
211  if (element_size == 8)
212  dffj(k) = (static_cast<tbox::CacheHex8 *>(cache)->getDsf(a))(k, j);
213  else if (element_size == 4)
214  dffj(k) = (static_cast<tbox::CacheTetra4 *>(cache)->getDsf(a))(k, j);
215 
216  factor2 = Jinv * dffj;
217  scalar dotProduct = factor3.dotproduct(factor2);
218  L_current(i, j) += dotProduct * gauss_w[a] * detJ;
219  }
220 
221  // S
222  factor3 = D * factor1;
223 
224  for (auto j = 0; j < element_size; ++j)
225  for (auto k = 0; k < 3; ++k)
226  S_current(i * 3 + k, j) += factor3(k) * ff(j) * gauss_w[a] * detJ;
227  }
228  }
229 
230  for (auto i = 0; i < 3 * element_size; ++i)
231  for (auto j = i; j < 3 * element_size; ++j)
232  {
233  if (fabs(K_current(i, j)) < epsilon)
234  K_current(i, j) = 0;
235  if (i != j)
236  K_current(j, i) = K_current(i, j);
237  }
238 
239  for (auto i = 0; i < element_size; ++i)
240  for (auto j = i; j < element_size; ++j)
241  {
242  if (fabs(L_current(i, j)) < epsilon)
243  L_current(i, j) = 0;
244  if (i != j)
245  L_current(j, i) = L_current(i, j);
246  }
247 
248  for (auto i = 0; i < 3 * element_size; ++i)
249  for (auto j = 0; j < element_size; ++j)
250  if (fabs(S_current(i, j)) < epsilon)
251  S_current(i, j) = 0;
252 }
253 
263 template <typename scalar, int element_type, int element_size>
264 Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> ElementMatrices<scalar, element_type, element_size>::getK(size_t i_thread) const
265 {
266  Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> K_current = subview(K, i_thread, Kokkos::ALL(), Kokkos::ALL());
267  return K_current;
268 }
269 
279 template <typename scalar, int element_type, int element_size>
280 Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> ElementMatrices<scalar, element_type, element_size>::getL(size_t i_thread) const
281 {
282  Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> L_current = subview(L, i_thread, Kokkos::ALL(), Kokkos::ALL());
283  return L_current;
284 }
285 
295 template <typename scalar, int element_type, int element_size>
296 Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> ElementMatrices<scalar, element_type, element_size>::getS(size_t i_thread) const
297 {
298  Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> S_current = subview(S, i_thread, Kokkos::ALL(), Kokkos::ALL());
299  return S_current;
300 }
katoptron::Domain
Class which is used to store all the information related to the discretized domain:
Definition: Domain.h:26
katoptron::Material::getNu
scalar getNu()
Return .
Definition: Material.h:45
ElementMatrices::material
view_type_2D material
Definition: ElementMatrices.h:54
ElementMatrices::S
view_type_3D S
Definition: ElementMatrices.h:52
katoptron::Material
Class which includes all constitutive values of a given material.
Definition: Material.h:26
ElementMatrices::L
view_type_3D L
Definition: ElementMatrices.h:52
tMatrix::clean
void clean()
Definition: tMatrix.h:29
EnsembleTraits
Definition: EnsembleTraits.h:8
tVector::dotproduct
T dotproduct(tVector rhs)
Definition: tMatrix.h:188
ElementMatrices::getS
view_type_2D getS(size_t i_thread) const
Return the coupling element matrix computed by the thread i_thread.
Definition: ElementMatrices.hpp:296
ElementMatrices::getK
view_type_2D getK(size_t i_thread) const
Return the element stifness matrix computed by the thread i_thread.
Definition: ElementMatrices.hpp:264
ElementMatrices::ElementMatrices
ElementMatrices(Teuchos::RCP< katoptron::Domain< scalar >> domain, size_t numPrimalDPN, size_t pool_size)
ElementMatrices constructor.
Definition: ElementMatrices.hpp:17
tVector
Definition: tMatrix.h:7
tMatrix
Definition: tMatrix.h:10
katoptron::Material::getE
scalar getE()
Return .
Definition: Material.h:40
tMatrix::inv
tMatrix inv()
Definition: tMatrix.h:97
ElementMatrices::getL
view_type_2D getL(size_t i_thread) const
Return the thermal element matrix computed by the thread i_thread.
Definition: ElementMatrices.hpp:280
ElementMatrices::K
view_type_3D K
Definition: ElementMatrices.h:52
katoptron::Material::getK
scalar getK()
Return .
Definition: Material.h:50
katoptron::Material::getBeta
scalar getBeta()
Return .
Definition: Material.h:55
ElementMatrices::compute
void compute(int e, int i_thread)
Compute all the matrices and store them as member data. This function must be called before getK,...
Definition: ElementMatrices.hpp:43