167 {
169
170 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
173 }
174
175
176 if (row_type != MBVERTEX)
178
180 if (nb_dofs == 0)
182
183 {
184
187 dot_W.resize(3,
false);
188 a_res.resize(3,
false);
189 g.resize(3, 3,
false);
190 G.resize(3, 3,
false);
191 h.resize(3, 3,
false);
192 H.resize(3, 3,
false);
193 invH.resize(3, 3,
false);
194 F.resize(3, 3,
false);
195 }
196
200 for (
int dd = 0;
dd < 3;
dd++) {
203 }
204
205 int nb_gauss_pts = row_data.
getN().size1();
209
210 const std::vector<VectorDouble> &dot_spacial_vel =
212
213 const std::vector<MatrixDouble> &spatial_positions_grad =
215
216 const std::vector<MatrixDouble> &spatial_velocities_grad =
218
219 const std::vector<VectorDouble> &meshpos_vel =
221
222 const std::vector<MatrixDouble> &mesh_positions_gradient =
224
225 int nb_active_vars = 0;
226 for (int gg = 0; gg < nb_gauss_pts; gg++) {
227
228 if (gg == 0) {
229
231
232 for (int nn1 = 0; nn1 < 3; nn1++) {
233
234 a[nn1] <<= dot_spacial_vel[gg][nn1];
235 nb_active_vars++;
236 }
237 for (int nn1 = 0; nn1 < 3; nn1++) {
238 for (int nn2 = 0; nn2 < 3; nn2++) {
239
240 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
242 if (nn1 == nn2) {
244 }
245 }
246 nb_active_vars++;
247 }
248 }
250 .size() > 0) {
251 for (int nn1 = 0; nn1 < 3; nn1++) {
252 for (int nn2 = 0; nn2 < 3; nn2++) {
253
254 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
255 nb_active_vars++;
256 }
257 }
258 for (int nn1 = 0; nn1 < 3; nn1++) {
259
260 dot_W(nn1) <<= meshpos_vel[gg][nn1];
261 nb_active_vars++;
262 }
263 for (int nn1 = 0; nn1 < 3; nn1++) {
264 for (int nn2 = 0; nn2 < 3; nn2++) {
265
266 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
267 nb_active_vars++;
268 }
269 }
270 }
271
274
275 auto t_a_res =
279 auto t_dotW =
286
288
291
292 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
293 t_a_res(
i) = t_a(
i) - t_a0(
i) + t_G(
i,
j) * t_dotW(
j);
294
296
297 t_F(
i,
j) = t_h(
i,
k)*t_invH(
k,
j);
298 t_a_res(
i) *= rho0 * detH;
300
301 } else {
302
303 t_a_res(
i) *= rho0 * detH;
304
305 }
306
307
309 res.resize(3);
310 for (int rr = 0; rr < 3; rr++) {
311 a_res[rr] >>= res[rr];
312 }
313
314 trace_off();
315 }
316
317 active.resize(nb_active_vars);
318 int aa = 0;
319 for (int nn1 = 0; nn1 < 3; nn1++) {
320 active[aa++] = dot_spacial_vel[gg][nn1];
321 }
322 for (int nn1 = 0; nn1 < 3; nn1++) {
323 for (int nn2 = 0; nn2 < 3; nn2++) {
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
326 } else {
327 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
328 }
329 }
330 }
332 0) {
333 for (int nn1 = 0; nn1 < 3; nn1++) {
334 for (int nn2 = 0; nn2 < 3; nn2++) {
335 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
336 }
337 }
338 for (int nn1 = 0; nn1 < 3; nn1++) {
339 active[aa++] = meshpos_vel[gg][nn1];
340 }
341 for (int nn1 = 0; nn1 < 3; nn1++) {
342 for (int nn2 = 0; nn2 < 3; nn2++) {
343 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
344 }
345 }
346 }
347
350 if (gg > 0) {
351 res.resize(3);
353 r = ::function(
tAg, 3, nb_active_vars, &
active[0], &res[0]);
356 "ADOL-C function evaluation with error r = %d",
r);
357 }
358 }
359 double val = getVolume() * getGaussPts()(3, gg);
360 res *= val;
361 } else {
364 for (int nn1 = 0; nn1 < 3; nn1++) {
367 }
369 r = jacobian(
tAg, 3, nb_active_vars, &
active[0],
373 "ADOL-C function evaluation with error");
374 }
375 double val = getVolume() * getGaussPts()(3, gg);
377 }
378 }
379 }
380
382}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_OPERATION_UNSUCCESSFUL
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
UBlasVector< double > VectorDouble
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
const double r
rate factor
Range tEts
elements in block set
double rho0
reference density
VectorDouble a0
constant acceleration
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< MatrixDouble > jacMass
std::vector< VectorDouble > valMass
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< std::vector< double * > > jacMassRowPtr
MatrixBoundedArray< adouble, 9 > invH
MatrixBoundedArray< adouble, 9 > G
FTensor::Index< 'k', 3 > k
FTensor::Index< 'i', 3 > i
MatrixBoundedArray< adouble, 9 > H
MatrixBoundedArray< adouble, 9 > g
MatrixBoundedArray< adouble, 9 > h
VectorBoundedArray< adouble, 3 > a_res
MatrixBoundedArray< adouble, 9 > F
VectorBoundedArray< adouble, 3 > dot_W
std::vector< double > active
FTensor::Index< 'j', 3 > j
VectorBoundedArray< adouble, 3 > a
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.