167 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
173 if (row_type != MBVERTEX)
184 dot_W.resize(3,
false);
185 a_res.resize(3,
false);
186 g.resize(3, 3,
false);
187 G.resize(3, 3,
false);
188 h.resize(3, 3,
false);
189 H.resize(3, 3,
false);
190 invH.resize(3, 3,
false);
191 F.resize(3, 3,
false);
195 std::fill(
H.data().begin(),
H.data().end(), 0);
196 std::fill(
invH.data().begin(),
invH.data().end(), 0);
197 for (
int ii = 0; ii != 3; ii++) {
202 int nb_gauss_pts = row_data.
getN().size1();
207 const std::vector<VectorDouble> &dot_spacial_vel =
210 const std::vector<MatrixDouble> &spatial_positions_grad =
213 const std::vector<MatrixDouble> &spatial_velocities_grad =
216 const std::vector<VectorDouble> &meshpos_vel =
219 const std::vector<MatrixDouble> &mesh_positions_gradient =
222 int nb_active_vars = 0;
223 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
229 for (
int nn1 = 0; nn1 < 3; nn1++) {
231 a[nn1] <<= dot_spacial_vel[gg][nn1];
234 for (
int nn1 = 0; nn1 < 3; nn1++) {
235 for (
int nn2 = 0; nn2 < 3; nn2++) {
237 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
248 for (
int nn1 = 0; nn1 < 3; nn1++) {
249 for (
int nn2 = 0; nn2 < 3; nn2++) {
251 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
255 for (
int nn1 = 0; nn1 < 3; nn1++) {
257 dot_W(nn1) <<= meshpos_vel[gg][nn1];
260 for (
int nn1 = 0; nn1 < 3; nn1++) {
261 for (
int nn2 = 0; nn2 < 3; nn2++) {
263 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
289 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
290 t_a_res(
i) = t_a(
i) - t_a0(
i) + t_G(
i,
j) * t_dotW(
j);
296 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
297 t_a_res(
i) *= rho0 * detH;
302 t_a_res(
i) *= rho0 * detH;
308 for (
int rr = 0; rr < 3; rr++) {
309 a_res[rr] >>= res[rr];
315 active.resize(nb_active_vars);
317 for (
int nn1 = 0; nn1 < 3; nn1++) {
318 active[aa++] = dot_spacial_vel[gg][nn1];
320 for (
int nn1 = 0; nn1 < 3; nn1++) {
321 for (
int nn2 = 0; nn2 < 3; nn2++) {
323 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
331 for (
int nn1 = 0; nn1 < 3; nn1++) {
332 for (
int nn2 = 0; nn2 < 3; nn2++) {
333 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
336 for (
int nn1 = 0; nn1 < 3; nn1++) {
337 active[aa++] = meshpos_vel[gg][nn1];
339 for (
int nn1 = 0; nn1 < 3; nn1++) {
340 for (
int nn2 = 0; nn2 < 3; nn2++) {
341 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
351 r = ::function(
tAg, 3, nb_active_vars, &
active[0], &res[0]);
354 "ADOL-C function evaluation with error r = %d",
r);
357 double val = getVolume() * getGaussPts()(3, gg);
363 for (
int nn1 = 0; nn1 < 3; nn1++) {
368 r = jacobian(
tAg, 3, nb_active_vars, &
active[0],
372 "ADOL-C function evaluation with error");
374 double val = getVolume() * getGaussPts()(3, gg);