164 {
166
167 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
170 }
171
172
173 if (row_type != MBVERTEX)
175
177 if (nb_dofs == 0)
179
180 {
181
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);
192 }
193
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++) {
200 }
201
202 int nb_gauss_pts = row_data.
getN().size1();
206
207 const std::vector<VectorDouble> &dot_spacial_vel =
209
210 const std::vector<MatrixDouble> &spatial_positions_grad =
212
213 const std::vector<MatrixDouble> &spatial_velocities_grad =
215
216 const std::vector<VectorDouble> &meshpos_vel =
218
219 const std::vector<MatrixDouble> &mesh_positions_gradient =
221
222 int nb_active_vars = 0;
223 for (int gg = 0; gg < nb_gauss_pts; gg++) {
224
225 if (gg == 0) {
226
228
229 for (int nn1 = 0; nn1 < 3; nn1++) {
230
231 a[nn1] <<= dot_spacial_vel[gg][nn1];
232 nb_active_vars++;
233 }
234 for (int nn1 = 0; nn1 < 3; nn1++) {
235 for (int nn2 = 0; nn2 < 3; nn2++) {
236
237 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
239 if (nn1 == nn2) {
241 }
242 }
243 nb_active_vars++;
244 }
245 }
247 .size() > 0) {
248 for (int nn1 = 0; nn1 < 3; nn1++) {
249 for (int nn2 = 0; nn2 < 3; nn2++) {
250
251 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
252 nb_active_vars++;
253 }
254 }
255 for (int nn1 = 0; nn1 < 3; nn1++) {
256
257 dot_W(nn1) <<= meshpos_vel[gg][nn1];
258 nb_active_vars++;
259 }
260 for (int nn1 = 0; nn1 < 3; nn1++) {
261 for (int nn2 = 0; nn2 < 3; nn2++) {
262
263 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
264 nb_active_vars++;
265 }
266 }
267 }
268
271
272 auto t_a_res =
276 auto t_dotW =
283
285
288
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);
291
292
293
295
296 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
297 t_a_res(
i) *= rho0 * detH;
299
300 } else {
301
302 t_a_res(
i) *= rho0 * detH;
303 }
304
305
307 res.resize(3);
308 for (int rr = 0; rr < 3; rr++) {
309 a_res[rr] >>= res[rr];
310 }
311
312 trace_off();
313 }
314
315 active.resize(nb_active_vars);
316 int aa = 0;
317 for (int nn1 = 0; nn1 < 3; nn1++) {
318 active[aa++] = dot_spacial_vel[gg][nn1];
319 }
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;
324 } else {
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
326 }
327 }
328 }
330 0) {
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);
334 }
335 }
336 for (int nn1 = 0; nn1 < 3; nn1++) {
337 active[aa++] = meshpos_vel[gg][nn1];
338 }
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);
342 }
343 }
344 }
345
348 if (gg > 0) {
349 res.resize(3);
351 r = ::function(
tAg, 3, nb_active_vars, &
active[0], &res[0]);
352 if (r != 3) {
354 "ADOL-C function evaluation with error r = %d", r);
355 }
356 }
357 double val = getVolume() * getGaussPts()(3, gg);
358 res *= val;
359
360 } else {
363 for (int nn1 = 0; nn1 < 3; nn1++) {
366 }
368 r = jacobian(
tAg, 3, nb_active_vars, &
active[0],
370 if (r != 3) {
372 "ADOL-C function evaluation with error");
373 }
374 double val = getVolume() * getGaussPts()(3, gg);
376 }
377 }
378 }
379
381}
#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 ...
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.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
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.