16 template <
typename scalar,
int element_type,
int element_size>
19 size_t pool_size) : domain(_domain), numPrimalDPN(_numPrimalDPN)
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);
26 material = Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>(
"A", pool_size, 3 * element_size);
42 template <
typename scalar,
int element_type,
int element_size>
49 if (element_type ==
static_cast<int>(tbox::ElType::HEX8))
51 tbox::CacheHex8 &cachehex8 = this->trilinosHex8GetCache();
54 else if (element_type ==
static_cast<int>(tbox::ElType::TETRA4))
56 tbox::CacheTetra4 &cachetetra4 = this->trilinosTetra4GetCache();
60 std::vector<Eigen::Vector3d> gauss_p;
61 std::vector<double> gauss_w;
62 for (
size_t a = 0; a < cache->getVGauss().getN(); ++a)
64 gauss_p.push_back(cache->getVGauss().getP(a));
65 gauss_w.push_back(cache->getVGauss().getW(a));
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());
74 Kokkos::deep_copy(K_current, 0);
75 Kokkos::deep_copy(L_current, 0);
76 Kokkos::deep_copy(S_current, 0);
78 auto element_material = domain->elementsList->getElementMaterial(e);
86 scalar E, nu, lambda, G, conductivity, beta;
91 if (domain->random_field->isRandom)
93 double x = 0., y = 0., z = 0.;
94 for (
auto j = 0; j < element_size; ++j)
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;
101 G = domain->random_field->operator()(x, y, z);
102 #ifdef ElementMatrices_DEBUG
107 nu = material.
getNu();
108 lambda = 2. * G * nu / (1. - 2. * nu);
114 nu = material.
getNu();
115 lambda = (nu * E) / ((1 + nu) * (1 - 2 * nu));
116 G = E / (2 * (1 + nu));
121 double epsilon = 2e-16;
123 for (
auto a = 0; a < gauss_p.size(); ++a)
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);
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)
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;
147 conductivity = material.
getK();
148 K_conductivity.
clean();
149 for (
auto i = 0; i < 3; ++i)
150 K_conductivity(i, i) = conductivity;
154 for (
auto i = 0; i < 3; ++i)
155 D(i, i) = -3 * (lambda + (2. / 3.) * G) * beta;
157 double detJ = this->buildJ(a, J, e, *(domain->elementsList), *(domain->nodesList));
159 #ifdef ElementMatrices_DEBUG
160 printf(
"detJ = %f\n", detJ);
165 for (
auto i = 0; i < element_size; ++i)
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));
179 for (
auto j = i; j < element_size; ++j)
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;
189 for (
auto k = 0; k < 3; ++k)
190 for (
auto m = 0; m < 3; ++m)
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);
196 factor3 = TMP1 * factor2;
198 scalar dotProduct = factor1.
dotproduct(factor3);
199 #ifdef ElementMatrices_DEBUG
202 K_current(i * 3 + k, j * 3 + m) += dotProduct * gauss_w[a] * detJ;
207 factor3 = K_conductivity * factor1;
208 for (
auto j = i; j < element_size; ++j)
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);
216 factor2 = Jinv * dffj;
217 scalar dotProduct = factor3.
dotproduct(factor2);
218 L_current(i, j) += dotProduct * gauss_w[a] * detJ;
222 factor3 = D * factor1;
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;
230 for (
auto i = 0; i < 3 * element_size; ++i)
231 for (
auto j = i; j < 3 * element_size; ++j)
233 if (fabs(K_current(i, j)) < epsilon)
236 K_current(j, i) = K_current(i, j);
239 for (
auto i = 0; i < element_size; ++i)
240 for (
auto j = i; j < element_size; ++j)
242 if (fabs(L_current(i, j)) < epsilon)
245 L_current(j, i) = L_current(i, j);
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)
263 template <
typename scalar,
int element_type,
int element_size>
266 Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> K_current = subview(K, i_thread, Kokkos::ALL(), Kokkos::ALL());
279 template <
typename scalar,
int element_type,
int element_size>
282 Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> L_current = subview(L, i_thread, Kokkos::ALL(), Kokkos::ALL());
295 template <
typename scalar,
int element_type,
int element_size>
298 Kokkos::View<scalar **, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> S_current = subview(S, i_thread, Kokkos::ALL(), Kokkos::ALL());