v0.14.0
ConvectiveMassElement.cpp
Go to the documentation of this file.
1 /** \file ConvectiveMassElement.cpp
2  * \brief Operators and data structures for mass and convective mass element
3  * \ingroup convective_mass_elem
4  *
5  */
6 
7 #include <MoFEM.hpp>
8 using namespace MoFEM;
9 
11 
12 #include <adolc/adolc.h>
14 #include <DirichletBC.hpp>
17 
18 #ifndef WITH_ADOL_C
19 #error "MoFEM need to be compiled with ADOL-C"
20 #endif
21 
23  : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULL), F(PETSC_NULL) {
24  meshPositionsFieldName = "NoNE";
25 
26  auto create_vec = [&]() {
27  constexpr int ghosts[] = {0};
28  if (mField.get_comm_rank() == 0) {
29  return createVectorMPI(mField.get_comm(), 1, 1);
30  } else {
31  return createVectorMPI(mField.get_comm(), 0, 1);
32  }
33  };
34 
35  V = create_vec();
36 }
37 
39 
42 
43  CHKERR VolumeElementForcesAndSourcesCore::preProcess();
44 
45  switch (ts_ctx) {
46  case CTX_TSNONE:
47  CHKERR VecZeroEntries(V);
48  break;
49  default:
50  break;
51  }
52 
54 }
55 
58 
59  CHKERR VolumeElementForcesAndSourcesCore::postProcess();
60 
61  const double *array;
62  switch (ts_ctx) {
63  case CTX_TSNONE:
64  CHKERR VecAssemblyBegin(V);
65  CHKERR VecAssemblyEnd(V);
66  CHKERR VecSum(V, &eNergy);
67  break;
68  default:
69  break;
70  }
71 
73 }
74 
76  short int tag)
77  : feMassRhs(m_field), feMassLhs(m_field), feMassAuxLhs(m_field),
78  feVelRhs(m_field), feVelLhs(m_field), feTRhs(m_field), feTLhs(m_field),
79  feEnergy(m_field), mField(m_field), tAg(tag) {}
80 
82  const std::string field_name,
83  std::vector<VectorDouble> &values_at_gauss_pts,
84  std::vector<MatrixDouble> &gardient_at_gauss_pts)
87  valuesAtGaussPts(values_at_gauss_pts),
88  gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
89 
91  int side, EntityType type, EntitiesFieldData::EntData &data) {
93 
94  int nb_dofs = data.getFieldData().size();
95  if (nb_dofs == 0) {
97  }
98  int nb_gauss_pts = data.getN().size1();
99  int nb_base_functions = data.getN().size2();
100 
101  // initialize
102  // VectorDouble& values = data.getFieldData();
103  valuesAtGaussPts.resize(nb_gauss_pts);
104  gradientAtGaussPts.resize(nb_gauss_pts);
105  for (int gg = 0; gg < nb_gauss_pts; gg++) {
106  valuesAtGaussPts[gg].resize(3);
107  gradientAtGaussPts[gg].resize(3, 3);
108  }
109 
110  if (type == zeroAtType) {
111  for (int gg = 0; gg < nb_gauss_pts; gg++) {
112  valuesAtGaussPts[gg].clear();
113  gradientAtGaussPts[gg].clear();
114  }
115  }
116 
117  auto base_function = data.getFTensor0N();
118  auto diff_base_functions = data.getFTensor1DiffN<3>();
121 
122  for (int gg = 0; gg != nb_gauss_pts; gg++) {
123  auto field_data = data.getFTensor1FieldData<3>();
124  FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
125  &valuesAtGaussPts[gg][1],
126  &valuesAtGaussPts[gg][2]);
128  &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
129  &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
130  &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
131  &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
132  &gradientAtGaussPts[gg](2, 2));
133  int bb = 0;
134  for (; bb != nb_dofs / 3; bb++) {
135  values(i) += base_function * field_data(i);
136  gradient(i, j) += field_data(i) * diff_base_functions(j);
137  ++diff_base_functions;
138  ++base_function;
139  ++field_data;
140  }
141  for (; bb != nb_base_functions; bb++) {
142  ++diff_base_functions;
143  ++base_function;
144  }
145  }
147 }
148 
150  const std::string field_name, CommonData &common_data)
151  : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
152  common_data.gradAtGaussPts[field_name]) {}
153 
155  const std::string field_name, BlockData &data, CommonData &common_data,
156  boost::ptr_vector<MethodForForceScaling> &methods_op, int tag,
157  bool jacobian)
160  dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
161  lInear(commonData.lInear), fieldDisp(false), methodsOp(methods_op) {}
162 
164  int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
166 
167  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
168  dAta.tEts.end()) {
170  }
171 
172  // do it only once, no need to repeat this for edges,faces or tets
173  if (row_type != MBVERTEX)
175 
176  int nb_dofs = row_data.getIndices().size();
177  if (nb_dofs == 0)
179 
180  {
181 
182  if (a.size() != 3) {
183  a.resize(3, false);
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 
194  std::fill(dot_W.begin(), dot_W.end(), 0);
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++) {
198  H(ii, ii) = 1;
199  invH(ii, ii) = 1;
200  }
201 
202  int nb_gauss_pts = row_data.getN().size1();
203  commonData.valMass.resize(nb_gauss_pts);
204  commonData.jacMassRowPtr.resize(nb_gauss_pts);
205  commonData.jacMass.resize(nb_gauss_pts);
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 
227  trace_on(tAg);
228 
229  for (int nn1 = 0; nn1 < 3; nn1++) { // 0
230  // commonData.dataAtGaussPts["DOT_"+commonData.spatialVelocities]
231  a[nn1] <<= dot_spacial_vel[gg][nn1];
232  nb_active_vars++;
233  }
234  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
235  for (int nn2 = 0; nn2 < 3; nn2++) {
236  // commonData.gradAtGaussPts[commonData.spatialPositions][gg]
237  h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
238  if (fieldDisp) {
239  if (nn1 == nn2) {
240  h(nn1, nn2) += 1;
241  }
242  }
243  nb_active_vars++;
244  }
245  }
247  .size() > 0) {
248  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
249  for (int nn2 = 0; nn2 < 3; nn2++) {
250  // commonData.gradAtGaussPts[commonData.spatialVelocities]
251  g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
252  nb_active_vars++;
253  }
254  }
255  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
256  // commonData.dataAtGaussPts["DOT_"+commonData.meshPositions]
257  dot_W(nn1) <<= meshpos_vel[gg][nn1];
258  nb_active_vars++;
259  }
260  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
261  for (int nn2 = 0; nn2 < 3; nn2++) {
262  // commonData.gradAtGaussPts[commonData.meshPositions][gg]
263  H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
264  nb_active_vars++;
265  }
266  }
267  }
268 
269  auto a0 = dAta.a0;
271 
272  auto t_a_res =
273  FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
274  auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
275  auto t_a0 = FTensor::Tensor1<double *, 3>{&a0[0], &a0[1], &a0[2]};
276  auto t_dotW =
277  FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
280  auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
283 
284  const double rho0 = dAta.rho0;
285 
287  CHKERR invertTensor3by3(H, detH, invH);
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  // FIXME: there is error somewhere for nonlinear case
293  // test dam example with -is_linear 0
294  if (!lInear) {
295 
296  t_F(i, j) = t_h(i, k) * t_invH(k, j);
297  t_a_res(i) *= rho0 * detH;
298  t_a_res(i) *= determinantTensor3by3(t_F);
299 
300  } else {
301 
302  t_a_res(i) *= rho0 * detH;
303  }
304 
305  // dependant
306  VectorDouble &res = commonData.valMass[gg];
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++) { // 0
318  active[aa++] = dot_spacial_vel[gg][nn1];
319  }
320  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
321  for (int nn2 = 0; nn2 < 3; nn2++) {
322  if (fieldDisp && nn1 == 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  }
329  if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
330  0) {
331  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
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++) { // 3+9+9=21
337  active[aa++] = meshpos_vel[gg][nn1];
338  }
339  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
340  for (int nn2 = 0; nn2 < 3; nn2++) {
341  active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
342  }
343  }
344  }
345 
346  if (!jAcobian) {
347  VectorDouble &res = commonData.valMass[gg];
348  if (gg > 0) {
349  res.resize(3);
350  int r;
351  r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
352  if (r != 3) { // function is locally analytic
353  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
354  "ADOL-C function evaluation with error r = %d", r);
355  }
356  }
357  double val = getVolume() * getGaussPts()(3, gg);
358  res *= val;
359  // cout << "my res " << res << endl;
360  } else {
361  commonData.jacMassRowPtr[gg].resize(3);
362  commonData.jacMass[gg].resize(3, nb_active_vars);
363  for (int nn1 = 0; nn1 < 3; nn1++) {
364  (commonData.jacMassRowPtr[gg])[nn1] =
365  &(commonData.jacMass[gg](nn1, 0));
366  }
367  int r;
368  r = jacobian(tAg, 3, nb_active_vars, &active[0],
369  &(commonData.jacMassRowPtr[gg])[0]);
370  if (r != 3) {
371  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
372  "ADOL-C function evaluation with error");
373  }
374  double val = getVolume() * getGaussPts()(3, gg);
375  commonData.jacMass[gg] *= val;
376  }
377  }
378  }
379 
381 }
382 
384  BlockData &data,
385  CommonData &common_data)
388  dAta(data), commonData(common_data) {}
389 
391 ConvectiveMassElement::OpMassRhs::doWork(int row_side, EntityType row_type,
392  EntitiesFieldData::EntData &row_data) {
394 
395  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
396  dAta.tEts.end()) {
398  }
399  if (row_data.getIndices().size() == 0)
401  int nb_dofs = row_data.getIndices().size();
402 
403  auto base = row_data.getFTensor0N();
404  int nb_base_functions = row_data.getN().size2();
405 
406  {
407 
408  nf.resize(nb_dofs);
409  nf.clear();
410 
412 
413  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
414  FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
416  &commonData.valMass[gg][1],
417  &commonData.valMass[gg][2]);
418  int dd = 0;
419  for (; dd < nb_dofs / 3; dd++) {
420  t_nf(i) += base * res(i);
421  ++base;
422  ++t_nf;
423  }
424  for (; dd != nb_base_functions; dd++) {
425  ++base;
426  }
427  }
428 
429  if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
430  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
431  }
432  CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, &row_data.getIndices()[0],
433  &nf[0], ADD_VALUES);
434  }
435 
437 }
438 
440  const std::string vel_field, const std::string field_name, BlockData &data,
441  CommonData &common_data, Range *forcesonlyonentities_ptr)
443  vel_field, field_name,
445  dAta(data), commonData(common_data) {
446  sYmm = false;
447  if (forcesonlyonentities_ptr != NULL) {
448  forcesOnlyOnEntities = *forcesonlyonentities_ptr;
449  }
450 }
451 
453  EntitiesFieldData::EntData &col_data, int gg) {
455  int nb_col = col_data.getIndices().size();
456  jac.clear();
457  if (!nb_col)
462  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
463  &jac(1, 0), &jac(1, 1), &jac(1, 2),
464  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
466  &commonData.jacMass[gg](0, 0), &commonData.jacMass[gg](0, 1),
467  &commonData.jacMass[gg](0, 2), &commonData.jacMass[gg](1, 0),
468  &commonData.jacMass[gg](1, 1), &commonData.jacMass[gg](1, 2),
469  &commonData.jacMass[gg](2, 0), &commonData.jacMass[gg](2, 1),
470  &commonData.jacMass[gg](2, 2));
471  double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
472  FTensor::Tensor0<double *> base(base_ptr, 1);
473  if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
474  0) {
475  for (int dd = 0; dd < nb_col / 3; dd++) {
476  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
477  ++base;
478  ++t_jac;
479  }
480  } else {
481  const int s = 3 + 9;
483  // T* d000, T* d001, T* d002,
484  // T* d010, T* d011, T* d012,
485  // T* d020, T* d021, T* d022,
486  // T* d100, T* d101, T* d102,
487  // T* d110, T* d111, T* d112,
488  // T* d120, T* d121, T* d122,
489  // T* d200, T* d201, T* d202,
490  // T* d210, T* d211, T* d212,
491  // T* d220, T* d221, T* d222,
492  &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
493  &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
494  &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
495  &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
496  &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
497  &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
498  &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
499  &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
500  &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
501  &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
502  &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
503  &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
504  &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
505  &commonData.jacMass[gg](2, s + 8));
506 
507  double *diff_ptr =
508  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
509  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
510  for (int dd = 0; dd < nb_col / 3; dd++) {
511  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
512  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
513  ++base;
514  ++diff;
515  ++t_jac;
516  }
517  }
519 }
520 
522  int row_side, int col_side, EntityType row_type, EntityType col_type,
523  EntitiesFieldData::EntData &row_data,
524  EntitiesFieldData::EntData &col_data) {
526 
527  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
528  dAta.tEts.end()) {
530  }
531 
532  int nb_row = row_data.getIndices().size();
533  int nb_col = col_data.getIndices().size();
534  if (nb_row == 0)
536  if (nb_col == 0)
538 
539  auto base = row_data.getFTensor0N();
540  int nb_base_functions = row_data.getN().size2();
541 
542  {
543 
544  k.resize(nb_row, nb_col);
545  k.clear();
546  jac.resize(3, nb_col);
547 
548  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
549 
550  try {
551  CHKERR getJac(col_data, gg);
552  } catch (const std::exception &ex) {
553  std::ostringstream ss;
554  ss << "throw in method: " << ex.what() << std::endl;
555  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
556  }
557 
560 
561  {
562  int dd1 = 0;
563  // integrate element stiffness matrix
564  for (; dd1 < nb_row / 3; dd1++) {
566  &jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
567  &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
568  for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
570  &k(3 * dd1 + 0, 3 * dd2 + 0), &k(3 * dd1 + 0, 3 * dd2 + 1),
571  &k(3 * dd1 + 0, 3 * dd2 + 2), &k(3 * dd1 + 1, 3 * dd2 + 0),
572  &k(3 * dd1 + 1, 3 * dd2 + 1), &k(3 * dd1 + 1, 3 * dd2 + 2),
573  &k(3 * dd1 + 2, 3 * dd2 + 0), &k(3 * dd1 + 2, 3 * dd2 + 1),
574  &k(3 * dd1 + 2, 3 * dd2 + 2));
575  t_k(i, j) += base * t_jac(i, j);
576  ++t_jac;
577  }
578  ++base;
579  // for(int rr1 = 0;rr1<3;rr1++) {
580  // for(int dd2 = 0;dd2<nb_col;dd2++) {
581  // k(3*dd1+rr1,dd2) += row_data.getN()(gg,dd1)*jac(rr1,dd2);
582  // }
583  // }
584  }
585  for (; dd1 != nb_base_functions; dd1++) {
586  ++base;
587  }
588  }
589  }
590 
591  if (!forcesOnlyOnEntities.empty()) {
592  VectorInt indices = row_data.getIndices();
593  VectorDofs &dofs = row_data.getFieldDofs();
594  VectorDofs::iterator dit = dofs.begin();
595  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
596  if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
597  forcesOnlyOnEntities.end()) {
598  indices[ii] = -1;
599  }
600  }
601  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, &indices[0], nb_col,
602  &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
603  } else {
604  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
605  &row_data.getIndices()[0], nb_col,
606  &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
607  }
608  }
610 }
611 
613  const std::string field_name, const std::string col_field, BlockData &data,
614  CommonData &common_data)
615  : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
616 
618  EntitiesFieldData::EntData &col_data, int gg) {
623  int nb_col = col_data.getIndices().size();
624  jac.clear();
625  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
626  &jac(1, 0), &jac(1, 1), &jac(1, 2),
627  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
628  const int s = 3;
630  // T* d000, T* d001, T* d002,
631  // T* d010, T* d011, T* d012,
632  // T* d020, T* d021, T* d022,
633  // T* d100, T* d101, T* d102,
634  // T* d110, T* d111, T* d112,
635  // T* d120, T* d121, T* d122,
636  // T* d200, T* d201, T* d202,
637  // T* d210, T* d211, T* d212,
638  // T* d220, T* d221, T* d222,
639  &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
640  &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
641  &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
642  &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
643  &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
644  &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
645  &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
646  &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
647  &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
648  &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
649  &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
650  &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
651  &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
652  &commonData.jacMass[gg](2, s + 8));
653  double *diff_ptr =
654  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
655  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
656  for (int dd = 0; dd < nb_col / 3; dd++) {
657  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
658  ++diff;
659  ++t_jac;
660  }
662 }
663 
665  const std::string field_name, const std::string col_field, BlockData &data,
666  CommonData &common_data)
667  : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
668 
670  EntitiesFieldData::EntData &col_data, int gg) {
672  int nb_col = col_data.getIndices().size();
673  jac.clear();
674  double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
675  FTensor::Tensor0<double *> base(base_ptr, 1);
676  double *diff_ptr =
677  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
678  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
679  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
680  &jac(1, 0), &jac(1, 1), &jac(1, 2),
681  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
682  const int u = 3 + 9 + 9;
684  &commonData.jacMass[gg](0, u + 0), &commonData.jacMass[gg](0, u + 1),
685  &commonData.jacMass[gg](0, u + 2), &commonData.jacMass[gg](1, u + 0),
686  &commonData.jacMass[gg](1, u + 1), &commonData.jacMass[gg](1, u + 2),
687  &commonData.jacMass[gg](2, u + 0), &commonData.jacMass[gg](2, u + 1),
688  &commonData.jacMass[gg](2, u + 2));
689  const int s = 3 + 9 + 9 + 3;
691  // T* d000, T* d001, T* d002,
692  // T* d010, T* d011, T* d012,
693  // T* d020, T* d021, T* d022,
694  // T* d100, T* d101, T* d102,
695  // T* d110, T* d111, T* d112,
696  // T* d120, T* d121, T* d122,
697  // T* d200, T* d201, T* d202,
698  // T* d210, T* d211, T* d212,
699  // T* d220, T* d221, T* d222,
700  &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
701  &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
702  &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
703  &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
704  &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
705  &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
706  &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
707  &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
708  &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
709  &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
710  &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
711  &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
712  &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
713  &commonData.jacMass[gg](2, s + 8));
717  for (int dd = 0; dd < nb_col / 3; dd++) {
718  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
719  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
720  ++base_ptr;
721  ++diff_ptr;
722  ++t_jac;
723  }
725 }
726 
728  BlockData &data,
729  CommonData &common_data,
733  dAta(data), commonData(common_data), V(v, true),
734  lInear(commonData.lInear) {}
735 
737 ConvectiveMassElement::OpEnergy::doWork(int row_side, EntityType row_type,
738  EntitiesFieldData::EntData &row_data) {
740 
741  if (row_type != MBVERTEX) {
743  }
744  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
745  dAta.tEts.end()) {
747  }
748 
749  {
750  double energy = 0;
751  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
752  double val = getVolume() * getGaussPts()(3, gg);
753  double rho0 = dAta.rho0;
754  double rho;
755  if (lInear) {
756  rho = rho0;
757  } else {
758  h.resize(3, 3);
759  noalias(h) =
762  .size() > 0) {
763  H.resize(3, 3);
764  noalias(H) =
766  auto detH = determinantTensor3by3(H);
767  invH.resize(3, 3);
768  CHKERR invertTensor3by3(H, detH, invH);
769  F.resize(3, 3);
770  noalias(F) = prod(h, invH);
771  } else {
772  F.resize(3, 3);
773  noalias(F) = h;
774  }
775  double detF = determinantTensor3by3(F);
776  rho = detF * rho0;
777  }
778  v.resize(3);
780  energy += 0.5 * (rho * val) * inner_prod(v, v);
781  }
782  CHKERR VecSetValue(V, 0, energy, ADD_VALUES);
783  }
784 
786 }
787 
789  const std::string field_name, BlockData &data, CommonData &common_data,
790  int tag, bool jacobian)
793  dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
794  fieldDisp(false) {}
795 
797  int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
799 
800  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
801  dAta.tEts.end()) {
803  }
804 
805  // do it only once, no need to repeat this for edges,faces or tets
806  if (row_type != MBVERTEX)
808 
809  int nb_dofs = row_data.getIndices().size();
810  if (nb_dofs == 0)
812 
813  {
814 
815  v.resize(3);
816  dot_w.resize(3);
817  h.resize(3, 3);
818  h.clear();
819  F.resize(3, 3);
820  dot_W.resize(3);
821  dot_W.clear();
822  H.resize(3, 3);
823  H.clear();
824  invH.resize(3, 3);
825  invH.clear();
826  dot_u.resize(3);
827  for (int dd = 0; dd < 3; dd++) {
828  H(dd, dd) = 1;
829  invH(dd, dd) = 1;
830  }
831 
832  a_res.resize(3);
833  int nb_gauss_pts = row_data.getN().size1();
834  commonData.valVel.resize(nb_gauss_pts);
835  commonData.jacVelRowPtr.resize(nb_gauss_pts);
836  commonData.jacVel.resize(nb_gauss_pts);
837 
838  int nb_active_vars = 0;
839  for (int gg = 0; gg < nb_gauss_pts; gg++) {
840 
841  if (gg == 0) {
842 
843  trace_on(tAg);
844 
845  for (int nn1 = 0; nn1 < 3; nn1++) { // 0
846  v[nn1] <<=
848  nb_active_vars++;
849  }
850  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
851  dot_w[nn1] <<=
853  [gg][nn1];
854  nb_active_vars++;
855  }
857  .size() > 0) {
858  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
859  for (int nn2 = 0; nn2 < 3; nn2++) {
860  h(nn1, nn2) <<=
862  nn1, nn2);
863  if (fieldDisp) {
864  if (nn1 == nn2) {
865  h(nn1, nn2) += 1;
866  }
867  }
868  nb_active_vars++;
869  }
870  }
871  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
872  dot_W[nn1] <<=
873  commonData
874  .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
875  nb_active_vars++;
876  }
877  }
879  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
880  for (int nn2 = 0; nn2 < 3; nn2++) {
881  H(nn1, nn2) <<=
883  nn2);
884  nb_active_vars++;
885  }
886  }
887  }
888  detH = 1;
889 
893 
897  auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
898  auto t_dot_u =
899  FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
900  auto t_dot_w =
901  FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
902  auto t_dot_W =
903  FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
904  auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
905  auto t_a_res =
906  FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
907 
909  detH = determinantTensor3by3(H);
910  CHKERR invertTensor3by3(H, detH, invH);
911  t_F(i, j) = t_h(i, k) * t_invH(k, j);
912  } else {
913  t_F(i, j) = t_h(i, j);
914  }
915 
916  t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
917  t_a_res(i) = t_v(i) - t_dot_u(i);
918  t_a_res(i) *= detH;
919 
920  // dependant
921  VectorDouble &res = commonData.valVel[gg];
922  res.resize(3);
923  for (int rr = 0; rr < 3; rr++) {
924  a_res[rr] >>= res[rr];
925  }
926  trace_off();
927  }
928 
929  active.resize(nb_active_vars);
930  int aa = 0;
931  for (int nn1 = 0; nn1 < 3; nn1++) {
932  active[aa++] =
934  }
935  for (int nn1 = 0; nn1 < 3; nn1++) {
936  active[aa++] =
937  commonData
938  .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
939  }
940  if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
941  0) {
942  for (int nn1 = 0; nn1 < 3; nn1++) {
943  for (int nn2 = 0; nn2 < 3; nn2++) {
944  if (fieldDisp && nn1 == nn2) {
945  active[aa++] =
947  nn1, nn2) +
948  1;
949  } else {
950  active[aa++] =
952  nn1, nn2);
953  }
954  }
955  }
956  for (int nn1 = 0; nn1 < 3; nn1++) {
957  active[aa++] =
958  commonData
959  .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
960  }
961  }
963  for (int nn1 = 0; nn1 < 3; nn1++) {
964  for (int nn2 = 0; nn2 < 3; nn2++) {
965  active[aa++] =
967  nn2);
968  }
969  }
970  }
971 
972  if (!jAcobian) {
973  VectorDouble &res = commonData.valVel[gg];
974  if (gg > 0) {
975  res.resize(3);
976  int r;
977  r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
978  if (r != 3) {
979  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
980  "ADOL-C function evaluation with error");
981  }
982  }
983  double val = getVolume() * getGaussPts()(3, gg);
984  res *= val;
985  } else {
986  commonData.jacVelRowPtr[gg].resize(3);
987  commonData.jacVel[gg].resize(3, nb_active_vars);
988  for (int nn1 = 0; nn1 < 3; nn1++) {
989  (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
990  }
991  int r;
992  r = jacobian(tAg, 3, nb_active_vars, &active[0],
993  &(commonData.jacVelRowPtr[gg])[0]);
994  if (r != 3) {
995  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
996  "ADOL-C function evaluation with error");
997  }
998  double val = getVolume() * getGaussPts()(3, gg);
999  commonData.jacVel[gg] *= val;
1000  // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1001  }
1002  }
1003  }
1005 }
1006 
1008  const std::string field_name, BlockData &data, CommonData &common_data)
1011  dAta(data), commonData(common_data) {}
1012 
1014  int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1016 
1017  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1018  dAta.tEts.end()) {
1020  }
1021  int nb_dofs = row_data.getIndices().size();
1022  if (nb_dofs == 0)
1024 
1025  auto base = row_data.getFTensor0N();
1026  int nb_base_functions = row_data.getN().size2();
1028 
1029  {
1030 
1031  nf.resize(nb_dofs);
1032  nf.clear();
1033 
1034  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1035  FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1037  &commonData.valVel[gg][1],
1038  &commonData.valVel[gg][2]);
1039  int dd = 0;
1040  for (; dd < nb_dofs / 3; dd++) {
1041  t_nf(i) += base * res(i);
1042  ++base;
1043  ++t_nf;
1044  }
1045  for (; dd != nb_base_functions; dd++) {
1046  ++base;
1047  }
1048  }
1049 
1050  if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1051  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1052  }
1053  CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1054  &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1055  }
1057 }
1058 
1060  const std::string vel_field, const std::string field_name, BlockData &data,
1061  CommonData &common_data)
1062  : OpMassLhs_dM_dv(vel_field, field_name, data, common_data) {}
1063 
1065  EntitiesFieldData::EntData &col_data, int gg) {
1067  int nb_col = col_data.getIndices().size();
1068  jac.clear();
1069  if (!nb_col)
1071  double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1072  FTensor::Tensor0<double *> base(base_ptr, 1);
1073  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1074  &jac(1, 0), &jac(1, 1), &jac(1, 2),
1075  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1077  &commonData.jacVel[gg](0, 0), &commonData.jacVel[gg](0, 1),
1078  &commonData.jacVel[gg](0, 2), &commonData.jacVel[gg](1, 0),
1079  &commonData.jacVel[gg](1, 1), &commonData.jacVel[gg](1, 2),
1080  &commonData.jacVel[gg](2, 0), &commonData.jacVel[gg](2, 1),
1081  &commonData.jacVel[gg](2, 2));
1084  for (int dd = 0; dd < nb_col / 3; dd++) {
1085  t_jac(i, j) += t_mass1(i, j) * base;
1086  ++base_ptr;
1087  ++t_jac;
1088  }
1089 
1091 }
1092 
1094  const std::string vel_field, const std::string field_name, BlockData &data,
1095  CommonData &common_data)
1096  : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1097 
1099  EntitiesFieldData::EntData &col_data, int gg) {
1101  int nb_col = col_data.getIndices().size();
1102  jac.clear();
1103  if (!nb_col)
1105  double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1106  FTensor::Tensor0<double *> base(base_ptr, 1);
1107  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1108  &jac(1, 0), &jac(1, 1), &jac(1, 2),
1109  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1110  const int u = 3;
1112  &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1113  &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1114  &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1115  &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1116  &commonData.jacVel[gg](2, u + 2));
1120  if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
1121  0) {
1122 
1123  for (int dd = 0; dd < nb_col / 3; dd++) {
1124  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1125  ++base_ptr;
1126  ++t_jac;
1127  }
1128  } else {
1129  double *diff_ptr =
1130  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1131  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1132  const int s = 3 + 3;
1134  // T* d000, T* d001, T* d002,
1135  // T* d010, T* d011, T* d012,
1136  // T* d020, T* d021, T* d022,
1137  // T* d100, T* d101, T* d102,
1138  // T* d110, T* d111, T* d112,
1139  // T* d120, T* d121, T* d122,
1140  // T* d200, T* d201, T* d202,
1141  // T* d210, T* d211, T* d212,
1142  // T* d220, T* d221, T* d222,
1143  &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1144  &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1145  &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1146  &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1147  &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1148  &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1149  &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1150  &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1151  &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1152  &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1153  &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1154  &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1155  &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1156  &commonData.jacVel[gg](2, s + 8));
1157  for (int dd = 0; dd < nb_col / 3; dd++) {
1158  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1159  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1160  ++base_ptr;
1161  ++diff_ptr;
1162  ++t_jac;
1163  }
1164  }
1166 }
1167 
1169  const std::string vel_field, const std::string field_name, BlockData &data,
1170  CommonData &common_data)
1171  : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1172 
1174  EntitiesFieldData::EntData &col_data, int gg) {
1176  int nb_col = col_data.getIndices().size();
1177  jac.clear();
1178  if (!nb_col)
1180  double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1181  FTensor::Tensor0<double *> base(base_ptr, 1);
1182  double *diff_ptr =
1183  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1184  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1185  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1186  &jac(1, 0), &jac(1, 1), &jac(1, 2),
1187  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1188  const int u = 3 + 3 + 9;
1190  &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1191  &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1192  &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1193  &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1194  &commonData.jacVel[gg](2, u + 2));
1195  const int s = 3 + 3 + 9 + 3;
1197  // T* d000, T* d001, T* d002,
1198  // T* d010, T* d011, T* d012,
1199  // T* d020, T* d021, T* d022,
1200  // T* d100, T* d101, T* d102,
1201  // T* d110, T* d111, T* d112,
1202  // T* d120, T* d121, T* d122,
1203  // T* d200, T* d201, T* d202,
1204  // T* d210, T* d211, T* d212,
1205  // T* d220, T* d221, T* d222,
1206  &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1207  &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1208  &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1209  &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1210  &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1211  &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1212  &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1213  &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1214  &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1215  &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1216  &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1217  &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1218  &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1219  &commonData.jacVel[gg](2, s + 8));
1223  for (int dd = 0; dd < nb_col / 3; dd++) {
1224  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1225  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1226  ++base_ptr;
1227  ++diff_ptr;
1228  ++t_jac;
1229  }
1231 }
1232 
1235  BlockData &data,
1236  CommonData &common_data, int tag,
1237  bool jacobian)
1240  dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
1241  fieldDisp(false) {}
1242 
1245  int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1247 
1248  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1249  dAta.tEts.end()) {
1251  }
1252 
1253  // do it only once, no need to repeat this for edges,faces or tets
1254  if (row_type != MBVERTEX)
1256 
1257  int nb_dofs = row_data.getIndices().size();
1258  if (nb_dofs == 0)
1260 
1261  try {
1262 
1263  a.resize(3);
1264  v.resize(3);
1265  g.resize(3, 3);
1266  G.resize(3, 3);
1267  h.resize(3, 3);
1268  F.resize(3, 3);
1269  H.resize(3, 3);
1270  H.clear();
1271  invH.resize(3, 3);
1272  invH.clear();
1273  for (int dd = 0; dd < 3; dd++) {
1274  H(dd, dd) = 1;
1275  invH(dd, dd) = 1;
1276  }
1277 
1278  int nb_gauss_pts = row_data.getN().size1();
1279  commonData.valT.resize(nb_gauss_pts);
1280  commonData.jacTRowPtr.resize(nb_gauss_pts);
1281  commonData.jacT.resize(nb_gauss_pts);
1282 
1283  int nb_active_vars = 0;
1284  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1285 
1286  if (gg == 0) {
1287 
1288  trace_on(tAg);
1289 
1290  for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1291  a[nn1] <<=
1293  [gg][nn1];
1294  nb_active_vars++;
1295  }
1296 
1297  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1298  v[nn1] <<=
1300  nb_active_vars++;
1301  }
1302  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1303  for (int nn2 = 0; nn2 < 3; nn2++) {
1304  g(nn1, nn2) <<=
1306  nn1, nn2);
1307  nb_active_vars++;
1308  }
1309  }
1310  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1311  for (int nn2 = 0; nn2 < 3; nn2++) {
1312  h(nn1, nn2) <<=
1314  nn2);
1315  nb_active_vars++;
1316  if (fieldDisp) {
1317  if (nn1 == nn2) {
1318  h(nn1, nn2) += 1;
1319  }
1320  }
1321  }
1322  }
1324  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1325  for (int nn2 = 0; nn2 < 3; nn2++) {
1326  H(nn1, nn2) <<=
1328  nn2);
1329  nb_active_vars++;
1330  }
1331  }
1332  }
1333  adouble detH;
1334  detH = 1;
1336  detH = determinantTensor3by3(H);
1337  }
1338  CHKERR invertTensor3by3(H, detH, invH);
1339 
1343 
1344  a_T.resize(3);
1345 
1346  auto t_h = getFTensor2FromArray3by3(h, FTensor::Number<0>(), 0);
1347  auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
1348  auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
1349  auto t_g = getFTensor2FromArray3by3(g, FTensor::Number<0>(), 0);
1350  auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
1351 
1352  auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
1353  auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
1354  auto t_a_T = FTensor::Tensor1<adouble *, 3>{&a_T[0], &a_T[1], &a_T[2]};
1355 
1356  t_F(i, j) = t_h(i, k) * t_invH(k, j);
1357  t_G(i, j) = t_g(i, k) * t_invH(k, j);
1358  t_a_T(i) = t_F(k, i) * t_a(k) + t_G(k, i) * t_v(k);
1359  const auto rho0 = dAta.rho0;
1360  t_a_T(i) = -rho0 * detH;
1361 
1362  commonData.valT[gg].resize(3);
1363  for (int nn = 0; nn < 3; nn++) {
1364  a_T[nn] >>= (commonData.valT[gg])[nn];
1365  }
1366  trace_off();
1367  }
1368 
1369  active.resize(nb_active_vars);
1370  int aa = 0;
1371  for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1372  active[aa++] =
1373  commonData
1374  .dataAtGaussPts["DOT_" + commonData.spatialVelocities][gg][nn1];
1375  }
1376 
1377  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1378  active[aa++] =
1380  }
1381  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1382  for (int nn2 = 0; nn2 < 3; nn2++) {
1383  active[aa++] =
1385  nn2);
1386  }
1387  }
1388  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1389  for (int nn2 = 0; nn2 < 3; nn2++) {
1390  if (fieldDisp && nn1 == nn2) {
1391  active[aa++] =
1393  nn1, nn2) +
1394  1;
1395  } else {
1396  active[aa++] =
1398  nn2);
1399  }
1400  }
1401  }
1403  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1404  for (int nn2 = 0; nn2 < 3; nn2++) {
1405  active[aa++] =
1407  nn2);
1408  }
1409  }
1410  }
1411 
1412  if (!jAcobian) {
1413  VectorDouble &res = commonData.valT[gg];
1414  if (gg > 0) {
1415  res.resize(3);
1416  int r;
1417  r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
1418  if (r != 3) { // function is locally analytic
1419  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1420  "ADOL-C function evaluation with error r = %d", r);
1421  }
1422  }
1423  double val = getVolume() * getGaussPts()(3, gg);
1424  res *= val;
1425  } else {
1426  commonData.jacTRowPtr[gg].resize(3);
1427  commonData.jacT[gg].resize(3, nb_active_vars);
1428  for (int nn1 = 0; nn1 < 3; nn1++) {
1429  (commonData.jacTRowPtr[gg])[nn1] = &(commonData.jacT[gg](nn1, 0));
1430  }
1431  int r;
1432  r = jacobian(tAg, 3, nb_active_vars, &active[0],
1433  &(commonData.jacTRowPtr[gg])[0]);
1434  if (r != 3) {
1435  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1436  "ADOL-C function evaluation with error");
1437  }
1438  double val = getVolume() * getGaussPts()(3, gg);
1439  commonData.jacT[gg] *= val;
1440  }
1441  }
1442 
1443  } catch (const std::exception &ex) {
1444  std::ostringstream ss;
1445  ss << "throw in method: " << ex.what() << std::endl;
1446  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1447  }
1448 
1450 }
1451 
1454  BlockData &data,
1455  CommonData &common_data,
1456  Range *forcesonlyonentities_ptr)
1459  dAta(data), commonData(common_data) {
1460  if (forcesonlyonentities_ptr != NULL) {
1461  forcesOnlyOnEntities = *forcesonlyonentities_ptr;
1462  }
1463 }
1464 
1467  int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1469 
1470  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1471  dAta.tEts.end()) {
1473  }
1474  int nb_dofs = row_data.getIndices().size();
1475  if (nb_dofs == 0)
1477 
1478  try {
1479 
1480  nf.resize(nb_dofs);
1481  nf.clear();
1482 
1483  auto base = row_data.getFTensor0N();
1484  int nb_base_functions = row_data.getN().size2();
1486 
1487  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1488  FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1490  &commonData.valT[gg][1],
1491  &commonData.valT[gg][2]);
1492  int dd = 0;
1493  for (; dd < nb_dofs / 3; dd++) {
1494  t_nf(i) += base * res(i);
1495  ++base;
1496  ++t_nf;
1497  }
1498  for (; dd != nb_base_functions; dd++) {
1499  ++base;
1500  }
1501  }
1502 
1503  if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1504  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1505  }
1506  if (!forcesOnlyOnEntities.empty()) {
1507  VectorInt indices = row_data.getIndices();
1508  VectorDofs &dofs = row_data.getFieldDofs();
1509  VectorDofs::iterator dit = dofs.begin();
1510  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1511  if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1512  forcesOnlyOnEntities.end()) {
1513  // std::cerr << **dit << std::endl;
1514  indices[ii] = -1;
1515  }
1516  }
1517  // std::cerr << indices << std::endl;
1518  CHKERR VecSetValues(getFEMethod()->ts_F, indices.size(), &indices[0],
1519  &nf[0], ADD_VALUES);
1520  } else {
1521  CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1522  &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1523  }
1524 
1525  } catch (const std::exception &ex) {
1526  std::ostringstream ss;
1527  ss << "throw in method: " << ex.what() << std::endl;
1528  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1529  }
1530 
1532 }
1533 
1535  OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field,
1536  const std::string field_name,
1537  BlockData &data,
1538  CommonData &common_data,
1539  Range *forcesonlyonentities_ptr)
1541  vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1542 
1545  EntitiesFieldData::EntData &col_data, int gg) {
1547  int nb_col = col_data.getIndices().size();
1548  jac.clear();
1549  double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1550  FTensor::Tensor0<double *> base(base_ptr, 1);
1551  double *diff_ptr =
1552  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1553  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1554  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1555  &jac(1, 0), &jac(1, 1), &jac(1, 2),
1556  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1557  const int u = 3;
1559  &commonData.jacT[gg](0, u + 0), &commonData.jacT[gg](0, u + 1),
1560  &commonData.jacT[gg](0, u + 2), &commonData.jacT[gg](1, u + 0),
1561  &commonData.jacT[gg](1, u + 1), &commonData.jacT[gg](1, u + 2),
1562  &commonData.jacT[gg](2, u + 0), &commonData.jacT[gg](2, u + 1),
1563  &commonData.jacT[gg](2, u + 2));
1564  const int s = 3 + 3;
1566  // T* d000, T* d001, T* d002,
1567  // T* d010, T* d011, T* d012,
1568  // T* d020, T* d021, T* d022,
1569  // T* d100, T* d101, T* d102,
1570  // T* d110, T* d111, T* d112,
1571  // T* d120, T* d121, T* d122,
1572  // T* d200, T* d201, T* d202,
1573  // T* d210, T* d211, T* d212,
1574  // T* d220, T* d221, T* d222,
1575  &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1576  &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1577  &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1578  &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1579  &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1580  &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1581  &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1582  &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1583  &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1584  &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1585  &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1586  &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1587  &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1588  &commonData.jacT[gg](2, s + 8));
1592  for (int dd = 0; dd < nb_col / 3; dd++) {
1593  t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1594  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1595  ++base_ptr;
1596  ++diff_ptr;
1597  ++t_jac;
1598  }
1600 }
1601 
1603  OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field,
1604  const std::string field_name,
1605  BlockData &data,
1606  CommonData &common_data,
1607  Range *forcesonlyonentities_ptr)
1609  vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1610 
1613  EntitiesFieldData::EntData &col_data, int gg) {
1615  int nb_col = col_data.getIndices().size();
1616  jac.clear();
1617  double *diff_ptr =
1618  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1619  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1620  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1621  &jac(1, 0), &jac(1, 1), &jac(1, 2),
1622  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1623  const int s = 3 + 3 + 9;
1625  // T* d000, T* d001, T* d002,
1626  // T* d010, T* d011, T* d012,
1627  // T* d020, T* d021, T* d022,
1628  // T* d100, T* d101, T* d102,
1629  // T* d110, T* d111, T* d112,
1630  // T* d120, T* d121, T* d122,
1631  // T* d200, T* d201, T* d202,
1632  // T* d210, T* d211, T* d212,
1633  // T* d220, T* d221, T* d222,
1634  &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1635  &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1636  &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1637  &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1638  &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1639  &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1640  &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1641  &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1642  &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1643  &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1644  &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1645  &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1646  &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1647  &commonData.jacT[gg](2, s + 8));
1651  for (int dd = 0; dd < nb_col / 3; dd++) {
1652  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1653  ++diff_ptr;
1654  ++t_jac;
1655  }
1657 }
1658 
1660  OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field,
1661  const std::string field_name,
1662  BlockData &data,
1663  CommonData &common_data,
1664  Range *forcesonlyonentities_ptr)
1666  vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1667 
1670  EntitiesFieldData::EntData &col_data, int gg) {
1672  int nb_col = col_data.getIndices().size();
1673  jac.clear();
1674  double *diff_ptr =
1675  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1676  FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1677  FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1678  &jac(1, 0), &jac(1, 1), &jac(1, 2),
1679  &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1680  const int s = 3 + 3 + 9 + 9;
1682  // T* d000, T* d001, T* d002,
1683  // T* d010, T* d011, T* d012,
1684  // T* d020, T* d021, T* d022,
1685  // T* d100, T* d101, T* d102,
1686  // T* d110, T* d111, T* d112,
1687  // T* d120, T* d121, T* d122,
1688  // T* d200, T* d201, T* d202,
1689  // T* d210, T* d211, T* d212,
1690  // T* d220, T* d221, T* d222,
1691  &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1692  &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1693  &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1694  &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1695  &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1696  &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1697  &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1698  &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1699  &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1700  &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1701  &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1702  &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1703  &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1704  &commonData.jacT[gg](2, s + 8));
1708  for (int dd = 0; dd < nb_col / 3; dd++) {
1709  t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1710  ++diff_ptr;
1711  ++t_jac;
1712  }
1714 }
1715 
1717  MoFEM::Interface &m_field, TS _ts, const std::string velocity_field,
1718  const std::string spatial_position_field)
1719  : mField(m_field), tS(_ts), velocityField(velocity_field),
1720  spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1721 
1724 
1725  switch (ts_ctx) {
1726  case CTX_TSSETIFUNCTION: {
1727  snes_f = ts_F;
1728  // FIXME: This global scattering because Kuu problem and Dynamic problem
1729  // not share partitions. Both problem should use the same partitioning to
1730  // resolve this problem.
1731  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
1732  problemPtr, COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1733  CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1734  problemPtr, velocityField, "DOT_" + velocityField, COL, ts_u_t,
1735  INSERT_VALUES, SCATTER_REVERSE);
1736  CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1737  problemPtr, spatialPositionField, "DOT_" + spatialPositionField, COL,
1738  ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1739  break;
1740  }
1741  case CTX_TSSETIJACOBIAN: {
1742  snes_B = ts_B;
1743  break;
1744  }
1745  default:
1746  break;
1747  }
1748 
1750 }
1751 
1754  //
1755  // SNES snes;
1756  // CHKERR TSGetSNES(tS,&snes);
1757  // CHKERR SNESSetLagJacobian(snes,jacobianLag);
1759 }
1760 
1763 
1764  Range added_tets;
1766  mField, BLOCKSET | BODYFORCESSET, it)) {
1767  int id = it->getMeshsetId();
1768  EntityHandle meshset = it->getMeshset();
1769  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1770  setOfBlocks[id].tEts, true);
1771  added_tets.merge(setOfBlocks[id].tEts);
1772  Block_BodyForces mydata;
1773  CHKERR it->getAttributeDataStructure(mydata);
1774  setOfBlocks[id].rho0 = mydata.data.density;
1775  setOfBlocks[id].a0.resize(3);
1776  setOfBlocks[id].a0[0] = mydata.data.acceleration_x;
1777  setOfBlocks[id].a0[1] = mydata.data.acceleration_y;
1778  setOfBlocks[id].a0[2] = mydata.data.acceleration_z;
1779  // std::cerr << setOfBlocks[id].tEts << std::endl;
1780  }
1781 
1783  mField, BLOCKSET | MAT_ELASTICSET, it)) {
1784  Mat_Elastic mydata;
1785  CHKERR it->getAttributeDataStructure(mydata);
1786  if (mydata.data.User1 == 0)
1787  continue;
1788  Range tets;
1789  EntityHandle meshset = it->getMeshset();
1790  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets, true);
1791  tets = subtract(tets, added_tets);
1792  if (tets.empty())
1793  continue;
1794  int id = it->getMeshsetId();
1795  setOfBlocks[-id].tEts = tets;
1796  setOfBlocks[-id].rho0 = mydata.data.User1;
1797  setOfBlocks[-id].a0.resize(3);
1798  setOfBlocks[-id].a0[0] = mydata.data.User2;
1799  setOfBlocks[-id].a0[1] = mydata.data.User3;
1800  setOfBlocks[-id].a0[2] = mydata.data.User4;
1801  // std::cerr << setOfBlocks[id].tEts << std::endl;
1802  }
1803 
1805 }
1806 
1808  MoFEM::Interface &m_field,
1809  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1811 
1812  if (!block_sets_ptr)
1813  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1814  "Pointer to block of sets is null");
1815 
1817  m_field, BLOCKSET | BODYFORCESSET, it)) {
1818  Block_BodyForces mydata;
1819  CHKERR it->getAttributeDataStructure(mydata);
1820  int id = it->getMeshsetId();
1821  auto &block_data = (*block_sets_ptr)[id];
1822  EntityHandle meshset = it->getMeshset();
1823  CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
1824  block_data.tEts, true);
1825  block_data.rho0 = mydata.data.density;
1826  block_data.a0.resize(3);
1827  block_data.a0[0] = mydata.data.acceleration_x;
1828  block_data.a0[1] = mydata.data.acceleration_y;
1829  block_data.a0[2] = mydata.data.acceleration_z;
1830  }
1831 
1833 }
1834 
1836  string element_name, string velocity_field_name,
1837  string spatial_position_field_name, string material_position_field_name,
1838  bool ale, BitRefLevel bit) {
1840 
1841  //
1842 
1843  CHKERR mField.add_finite_element(element_name, MF_ZERO);
1845  velocity_field_name);
1847  velocity_field_name);
1849  velocity_field_name);
1851  element_name, spatial_position_field_name);
1853  element_name, spatial_position_field_name);
1855  element_name, spatial_position_field_name);
1856  if (mField.check_field(material_position_field_name)) {
1857  if (ale) {
1859  element_name, material_position_field_name);
1861  element_name, material_position_field_name);
1863  element_name, "DOT_" + material_position_field_name);
1864  }
1866  element_name, material_position_field_name);
1867  }
1869  element_name, "DOT_" + velocity_field_name);
1871  element_name, "DOT_" + spatial_position_field_name);
1872 
1873  Range tets;
1874  if (bit.any()) {
1875  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1876  bit, BitRefLevel().set(), MBTET, tets);
1877  }
1878 
1879  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1880  for (; sit != setOfBlocks.end(); sit++) {
1881  Range add_tets = sit->second.tEts;
1882  if (!tets.empty()) {
1883  add_tets = intersect(add_tets, tets);
1884  }
1886  element_name);
1887  }
1888 
1890 }
1891 
1893  string element_name, string velocity_field_name,
1894  string spatial_position_field_name, string material_position_field_name,
1895  bool ale, BitRefLevel bit) {
1897 
1898  //
1899 
1900  CHKERR mField.add_finite_element(element_name, MF_ZERO);
1902  velocity_field_name);
1904  velocity_field_name);
1906  velocity_field_name);
1908  element_name, spatial_position_field_name);
1910  element_name, spatial_position_field_name);
1911  if (mField.check_field(material_position_field_name)) {
1912  if (ale) {
1914  element_name, material_position_field_name);
1916  element_name, "DOT_" + material_position_field_name);
1917  }
1919  element_name, material_position_field_name);
1920  }
1922  element_name, "DOT_" + velocity_field_name);
1924  element_name, "DOT_" + spatial_position_field_name);
1925 
1926  Range tets;
1927  if (bit.any()) {
1928  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1929  bit, BitRefLevel().set(), MBTET, tets);
1930  }
1931 
1932  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1933  for (; sit != setOfBlocks.end(); sit++) {
1934  Range add_tets = sit->second.tEts;
1935  if (!tets.empty()) {
1936  add_tets = intersect(add_tets, tets);
1937  }
1939  element_name);
1940  }
1941 
1943 }
1944 
1946  string element_name, string velocity_field_name,
1947  string spatial_position_field_name, string material_position_field_name,
1948  bool ale, BitRefLevel bit, Range *intersected) {
1950 
1951  //
1952 
1953  CHKERR mField.add_finite_element(element_name, MF_ZERO);
1955  velocity_field_name);
1957  velocity_field_name);
1959  element_name, spatial_position_field_name);
1961  element_name, spatial_position_field_name);
1962  if (mField.check_field(material_position_field_name)) {
1963  if (ale) {
1965  element_name, material_position_field_name);
1967  element_name, material_position_field_name);
1969  element_name, "DOT_" + material_position_field_name);
1970  }
1972  element_name, material_position_field_name);
1973  }
1975  element_name, "DOT_" + velocity_field_name);
1977  element_name, "DOT_" + spatial_position_field_name);
1978 
1979  Range tets;
1980  if (bit.any()) {
1981  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1982  bit, BitRefLevel().set(), MBTET, tets);
1983  }
1984  if (intersected != NULL) {
1985  if (tets.empty()) {
1986  tets = *intersected;
1987  } else {
1988  tets = intersect(*intersected, tets);
1989  }
1990  }
1991 
1992  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1993  for (; sit != setOfBlocks.end(); sit++) {
1994  Range add_tets = sit->second.tEts;
1995  if (!tets.empty()) {
1996  add_tets = intersect(add_tets, tets);
1997  }
1999  element_name);
2000  }
2001 
2003 }
2004 
2006  string velocity_field_name, string spatial_position_field_name,
2007  string material_position_field_name, bool ale, bool linear) {
2009 
2010  commonData.spatialPositions = spatial_position_field_name;
2011  commonData.meshPositions = material_position_field_name;
2012  commonData.spatialVelocities = velocity_field_name;
2013  commonData.lInear = linear;
2014 
2015  // Rhs
2016  feMassRhs.getOpPtrVector().push_back(
2017  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2018  feMassRhs.getOpPtrVector().push_back(
2019  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2020  feMassRhs.getOpPtrVector().push_back(
2021  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2022  feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2023  "DOT_" + spatial_position_field_name, commonData));
2024  if (mField.check_field(material_position_field_name)) {
2025  feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2026  material_position_field_name, commonData));
2027  if (ale) {
2028  feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2029  "DOT_" + material_position_field_name, commonData));
2030  } else {
2031  feMassRhs.meshPositionsFieldName = material_position_field_name;
2032  }
2033  }
2034  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2035  for (; sit != setOfBlocks.end(); sit++) {
2036  feMassRhs.getOpPtrVector().push_back(
2037  new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2038  methodsOp, tAg, false));
2039  feMassRhs.getOpPtrVector().push_back(
2040  new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2041  }
2042 
2043  // Lhs
2044  feMassLhs.getOpPtrVector().push_back(
2045  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2046  feMassLhs.getOpPtrVector().push_back(
2047  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2048  feMassLhs.getOpPtrVector().push_back(
2049  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2050  feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2051  "DOT_" + spatial_position_field_name, commonData));
2052  if (mField.check_field(material_position_field_name)) {
2053  feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2054  material_position_field_name, commonData));
2055  if (ale) {
2056  feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2057  "DOT_" + material_position_field_name, commonData));
2058  } else {
2059  feMassLhs.meshPositionsFieldName = material_position_field_name;
2060  }
2061  }
2062  sit = setOfBlocks.begin();
2063  for (; sit != setOfBlocks.end(); sit++) {
2064  feMassLhs.getOpPtrVector().push_back(
2065  new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2066  methodsOp, tAg, true));
2067  feMassLhs.getOpPtrVector().push_back(
2068  new OpMassLhs_dM_dv(spatial_position_field_name, velocity_field_name,
2069  sit->second, commonData));
2070  feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2071  spatial_position_field_name, spatial_position_field_name, sit->second,
2072  commonData));
2073  if (mField.check_field(material_position_field_name)) {
2074  if (ale) {
2075  feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dX(
2076  spatial_position_field_name, material_position_field_name,
2077  sit->second, commonData));
2078  } else {
2079  feMassLhs.meshPositionsFieldName = material_position_field_name;
2080  }
2081  }
2082  }
2083 
2084  // Energy
2085  feEnergy.getOpPtrVector().push_back(
2086  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2087  feEnergy.getOpPtrVector().push_back(
2088  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2089  if (mField.check_field(material_position_field_name)) {
2090  feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2091  material_position_field_name, commonData));
2092  feEnergy.meshPositionsFieldName = material_position_field_name;
2093  }
2094  sit = setOfBlocks.begin();
2095  for (; sit != setOfBlocks.end(); sit++) {
2096  feEnergy.getOpPtrVector().push_back(new OpEnergy(
2097  spatial_position_field_name, sit->second, commonData, feEnergy.V));
2098  }
2099 
2101 }
2102 
2104  string velocity_field_name, string spatial_position_field_name,
2105  string material_position_field_name, bool ale) {
2107 
2108  commonData.spatialPositions = spatial_position_field_name;
2109  commonData.meshPositions = material_position_field_name;
2110  commonData.spatialVelocities = velocity_field_name;
2111 
2112  // Rhs
2113  feVelRhs.getOpPtrVector().push_back(
2114  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2115  feVelRhs.getOpPtrVector().push_back(
2116  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2117  feVelRhs.getOpPtrVector().push_back(
2118  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2119  if (mField.check_field(material_position_field_name)) {
2120  feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2121  "DOT_" + spatial_position_field_name, commonData));
2122  feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2123  material_position_field_name, commonData));
2124  if (ale) {
2125  feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2126  "DOT_" + material_position_field_name, commonData));
2127  } else {
2128  feVelRhs.meshPositionsFieldName = material_position_field_name;
2129  }
2130  }
2131  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2132  for (; sit != setOfBlocks.end(); sit++) {
2133  feVelRhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2134  velocity_field_name, sit->second, commonData, tAg, false));
2135  feVelRhs.getOpPtrVector().push_back(
2136  new OpVelocityRhs(velocity_field_name, sit->second, commonData));
2137  }
2138 
2139  // Lhs
2140  feVelLhs.getOpPtrVector().push_back(
2141  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2142  feVelLhs.getOpPtrVector().push_back(
2143  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2144  feVelLhs.getOpPtrVector().push_back(
2145  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2146  if (mField.check_field(material_position_field_name)) {
2147  feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2148  "DOT_" + spatial_position_field_name, commonData));
2149  feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2150  material_position_field_name, commonData));
2151  if (ale) {
2152  feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2153  "DOT_" + material_position_field_name, commonData));
2154  } else {
2155  feVelLhs.meshPositionsFieldName = material_position_field_name;
2156  }
2157  }
2158  sit = setOfBlocks.begin();
2159  for (; sit != setOfBlocks.end(); sit++) {
2160  feVelLhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2161  velocity_field_name, sit->second, commonData, tAg));
2162  feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dv(
2163  velocity_field_name, velocity_field_name, sit->second, commonData));
2164  feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dx(
2165  velocity_field_name, spatial_position_field_name, sit->second,
2166  commonData));
2167  if (mField.check_field(material_position_field_name)) {
2168  if (ale) {
2169  feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dX(
2170  velocity_field_name, material_position_field_name, sit->second,
2171  commonData));
2172  } else {
2173  feVelLhs.meshPositionsFieldName = material_position_field_name;
2174  }
2175  }
2176  }
2177 
2179 }
2180 
2182  string velocity_field_name, string spatial_position_field_name,
2183  string material_position_field_name, Range *forces_on_entities_ptr) {
2185 
2186  commonData.spatialPositions = spatial_position_field_name;
2187  commonData.meshPositions = material_position_field_name;
2188  commonData.spatialVelocities = velocity_field_name;
2189 
2190  // Rhs
2191  feTRhs.getOpPtrVector().push_back(
2192  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2193  feTRhs.getOpPtrVector().push_back(
2194  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2195  feTRhs.getOpPtrVector().push_back(
2196  new OpGetCommonDataAtGaussPts(material_position_field_name, commonData));
2197  feTRhs.getOpPtrVector().push_back(
2198  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2199 
2200  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2201  for (; sit != setOfBlocks.end(); sit++) {
2202  feTRhs.getOpPtrVector().push_back(
2204  material_position_field_name, sit->second, commonData, tAg, false));
2205  feTRhs.getOpPtrVector().push_back(new OpEshelbyDynamicMaterialMomentumRhs(
2206  material_position_field_name, sit->second, commonData,
2207  forces_on_entities_ptr));
2208  }
2209 
2210  // Lhs
2211  feTLhs.getOpPtrVector().push_back(
2212  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2213  feTLhs.getOpPtrVector().push_back(
2214  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2215  feTLhs.getOpPtrVector().push_back(
2216  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2217  if (mField.check_field(material_position_field_name)) {
2218  feTLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2219  material_position_field_name, commonData));
2220  }
2221  sit = setOfBlocks.begin();
2222  for (; sit != setOfBlocks.end(); sit++) {
2223  feTLhs.getOpPtrVector().push_back(
2225  material_position_field_name, sit->second, commonData, tAg));
2226  feTLhs.getOpPtrVector().push_back(
2228  material_position_field_name, velocity_field_name, sit->second,
2229  commonData, forces_on_entities_ptr));
2230  feTLhs.getOpPtrVector().push_back(
2232  material_position_field_name, spatial_position_field_name,
2233  sit->second, commonData, forces_on_entities_ptr));
2234  feTLhs.getOpPtrVector().push_back(
2236  material_position_field_name, material_position_field_name,
2237  sit->second, commonData, forces_on_entities_ptr));
2238  }
2239 
2241 }
2242 
2244  string velocity_field_name, string spatial_position_field_name,
2245  string material_position_field_name, bool linear) {
2247 
2248  commonData.spatialPositions = spatial_position_field_name;
2249  commonData.meshPositions = material_position_field_name;
2250  commonData.spatialVelocities = velocity_field_name;
2251  commonData.lInear = linear;
2252 
2253  // Rhs
2254  feMassRhs.getOpPtrVector().push_back(
2255  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2256  feMassRhs.getOpPtrVector().push_back(
2257  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2258  feMassRhs.getOpPtrVector().push_back(
2259  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2260  if (mField.check_field(material_position_field_name)) {
2261  feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2262  material_position_field_name, commonData));
2263  feMassRhs.meshPositionsFieldName = material_position_field_name;
2264  }
2265  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2266  for (; sit != setOfBlocks.end(); sit++) {
2267  feMassRhs.getOpPtrVector().push_back(
2268  new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2269  methodsOp, tAg, false));
2270  feMassRhs.getOpPtrVector().push_back(
2271  new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2272  }
2273 
2274  // Lhs
2275  feMassLhs.getOpPtrVector().push_back(
2276  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2277  feMassLhs.getOpPtrVector().push_back(
2278  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2279  feMassLhs.getOpPtrVector().push_back(
2280  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2281  if (mField.check_field(material_position_field_name)) {
2282  feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2283  material_position_field_name, commonData));
2284  feMassLhs.meshPositionsFieldName = material_position_field_name;
2285  }
2286  sit = setOfBlocks.begin();
2287  for (; sit != setOfBlocks.end(); sit++) {
2288  feMassLhs.getOpPtrVector().push_back(
2289  new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2290  methodsOp, tAg, true));
2291  feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dv(
2292  spatial_position_field_name, spatial_position_field_name, sit->second,
2293  commonData));
2294  if (mField.check_field(material_position_field_name)) {
2295  feMassLhs.meshPositionsFieldName = material_position_field_name;
2296  }
2297  }
2298 
2299  // Aux Lhs
2300  feMassAuxLhs.getOpPtrVector().push_back(
2301  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2302  feMassAuxLhs.getOpPtrVector().push_back(
2303  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2304  feMassAuxLhs.getOpPtrVector().push_back(
2305  new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2306  if (mField.check_field(material_position_field_name)) {
2307  feMassAuxLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2308  material_position_field_name, commonData));
2309  feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2310  }
2311  sit = setOfBlocks.begin();
2312  for (; sit != setOfBlocks.end(); sit++) {
2313  feMassAuxLhs.getOpPtrVector().push_back(
2314  new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2315  methodsOp, tAg, true));
2316  feMassAuxLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2317  spatial_position_field_name, spatial_position_field_name, sit->second,
2318  commonData));
2319  if (mField.check_field(material_position_field_name)) {
2320  feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2321  }
2322  }
2323 
2324  // Energy E=0.5*rho*v*v
2325  feEnergy.getOpPtrVector().push_back(
2326  new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2327  feEnergy.getOpPtrVector().push_back(
2328  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2329  if (mField.check_field(material_position_field_name)) {
2330  feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2331  material_position_field_name, commonData));
2332  feEnergy.meshPositionsFieldName = material_position_field_name;
2333  }
2334  sit = setOfBlocks.begin();
2335  for (; sit != setOfBlocks.end(); sit++) {
2336  feEnergy.getOpPtrVector().push_back(new OpEnergy(
2337  spatial_position_field_name, sit->second, commonData, feEnergy.V));
2338  }
2339 
2341 }
2342 
2345  if (iNitialized) {
2346 
2347  CHKERR dEstroy();
2348  CHKERRABORT(PETSC_COMM_WORLD, ierr);
2349  }
2350 }
2351 
2354  if (!iNitialized) {
2355 
2356 #if PETSC_VERSION_GE(3, 5, 3)
2357  CHKERR MatCreateVecs(K, &u, &Ku);
2358  CHKERR MatCreateVecs(M, &v, &Mv);
2359 #else
2360  CHKERR MatGetVecs(K, &u, &Ku);
2361  CHKERR MatGetVecs(M, &v, &Mv);
2362 #endif
2363  CHKERR MatDuplicate(K, MAT_SHARE_NONZERO_PATTERN, &barK);
2364  iNitialized = true;
2365  }
2367 }
2368 
2371  if (iNitialized) {
2372 
2373  CHKERR VecDestroy(&u);
2374  CHKERR VecDestroy(&Ku);
2375  CHKERR VecDestroy(&v);
2376  CHKERR VecDestroy(&Mv);
2377  CHKERR MatDestroy(&barK);
2378  iNitialized = false;
2379  }
2381 }
2382 
2385 
2386  if (!initPC) {
2387  MPI_Comm comm;
2388  CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2389  CHKERR PCCreate(comm, &pC);
2390  initPC = true;
2391  }
2393 }
2394 
2397 
2398  if (initPC) {
2399  CHKERR PCDestroy(&pC);
2400  initPC = false;
2401  }
2403 }
2404 
2406  MoFEM::Interface &m_field)
2407  : mField(m_field) {}
2408 
2411 
2412  if (ts_ctx != CTX_TSSETIFUNCTION) {
2413  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2414  "It is used to residual of velocities");
2415  }
2416  if (!shellMatCtx->iNitialized) {
2417  CHKERR shellMatCtx->iNit();
2418  }
2419  // Note velocities calculate from displacements are stroed in shellMatCtx->u
2420  CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2421  INSERT_VALUES, SCATTER_FORWARD);
2422  CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2423  INSERT_VALUES, SCATTER_FORWARD);
2424  CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2425  INSERT_VALUES, SCATTER_FORWARD);
2426  CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2427  INSERT_VALUES, SCATTER_FORWARD);
2428  CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2429  CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2430  ADD_VALUES, SCATTER_REVERSE);
2431  CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2432  SCATTER_REVERSE);
2433  // VecView(shellMatCtx->v,PETSC_VIEWER_STDOUT_WORLD);
2434 
2436 }
2437 
2438 #ifdef __DIRICHLET_HPP__
2439 
2440 ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2441  MoFEM::Interface &m_field)
2442  : mField(m_field) {}
2443 
2444 MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2446 
2447  if (ts_ctx != CTX_TSSETIJACOBIAN) {
2448  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2449  "It is used to calculate shell matrix only");
2450  }
2451 
2452  shellMatCtx->ts_a = ts_a;
2453  DirichletBcPtr->copyTs(*((TSMethod *)this)); // copy context for TSMethod
2454 
2455  DirichletBcPtr->dIag = 1;
2456  DirichletBcPtr->ts_B = shellMatCtx->K;
2457  CHKERR MatZeroEntries(shellMatCtx->K);
2458  CHKERR mField.problem_basic_method_preProcess(problemName, *DirichletBcPtr);
2459  LoopsToDoType::iterator itk = loopK.begin();
2460  for (; itk != loopK.end(); itk++) {
2461  itk->second->copyTs(*((TSMethod *)this));
2462  itk->second->ts_B = shellMatCtx->K;
2463  CHKERR mField.loop_finite_elements(problemName, itk->first, *itk->second);
2464  }
2465  LoopsToDoType::iterator itam = loopAuxM.begin();
2466  for (; itam != loopAuxM.end(); itam++) {
2467  itam->second->copyTs(*((TSMethod *)this));
2468  itam->second->ts_B = shellMatCtx->K;
2469  CHKERR mField.loop_finite_elements(problemName, itam->first, *itam->second);
2470  }
2471  CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2472  CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2473  CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2474 
2475  DirichletBcPtr->dIag = 0;
2476  DirichletBcPtr->ts_B = shellMatCtx->M;
2477  CHKERR MatZeroEntries(shellMatCtx->M);
2478  // CHKERR mField.problem_basic_method_preProcess(problemName,*DirichletBcPtr);
2479  LoopsToDoType::iterator itm = loopM.begin();
2480  for (; itm != loopM.end(); itm++) {
2481  itm->second->copyTs(*((TSMethod *)this));
2482  itm->second->ts_B = shellMatCtx->M;
2483  CHKERR mField.loop_finite_elements(problemName, itm->first, *itm->second);
2484  }
2485  CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2486  CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2487  CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2488 
2489  // barK
2490  CHKERR MatZeroEntries(shellMatCtx->barK);
2491  CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2492  CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2493  CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2494  CHKERR MatAssemblyEnd(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2495 
2496  // Matrix View
2497  // MatView(shellMatCtx->barK,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
2498  // std::string wait;
2499  // std::cin >> wait;
2500 
2502 }
2503 
2504 #endif //__DIRICHLET_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::Block_BodyForces
Body force data structure.
Definition: MaterialBlocks.hpp:313
ConvectiveMassElement::addEshelbyDynamicMaterialMomentum
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
Definition: ConvectiveMassElement.cpp:1945
MethodForForceScaling.hpp
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
g
constexpr double g
Definition: shallow_wave.cpp:63
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dv
Definition: ConvectiveMassElement.hpp:389
ConvectiveMassElement::OpMassRhs::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:391
MoFEM::CoreInterface::problem_basic_method_postProcess
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
ConvectiveMassElement::OpMassRhs
Definition: ConvectiveMassElement.hpp:207
ConvectiveMassElement::OpVelocityLhs_dV_dv
Definition: ConvectiveMassElement.hpp:317
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
ConvectiveMassElement::OpVelocityJacobian::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:796
ConvectiveMassElement::MatShellCtx::iNit
MoFEMErrorCode iNit()
Definition: ConvectiveMassElement.cpp:2352
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
ConvectiveMassElement::setShellMatrixMassOperators
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
Definition: ConvectiveMassElement.cpp:2243
EntityHandle
ConvectiveMassElement::OpMassRhs::OpMassRhs
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:383
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumRhs::forcesOnlyOnEntities
Range forcesOnlyOnEntities
Definition: ConvectiveMassElement.hpp:376
ConvectiveMassElement::OpMassLhs_dM_dv
Definition: ConvectiveMassElement.hpp:222
rho
double rho
Definition: plastic.cpp:140
ConvectiveMassElement::feVelLhs
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
Definition: ConvectiveMassElement.hpp:79
ConvectiveMassElement::setVelocityOperators
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
Definition: ConvectiveMassElement.cpp:2103
ConvectiveMassElement::OpMassLhs_dM_dx::OpMassLhs_dM_dx
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:612
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
ConvectiveMassElement::OpEnergy
Definition: ConvectiveMassElement.hpp:261
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
ConvectiveMassElement::OpVelocityRhs
Definition: ConvectiveMassElement.hpp:301
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
ConvectiveMassElement::CommonData::jacVelRowPtr
std::vector< std::vector< double * > > jacVelRowPtr
Definition: ConvectiveMassElement.hpp:139
ConvectiveMassElement::CommonData::jacT
std::vector< MatrixDouble > jacT
Definition: ConvectiveMassElement.hpp:146
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
ConvectiveMassElement::OpVelocityLhs_dV_dx::OpVelocityLhs_dV_dx
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:1093
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
ConvectiveMassElement::feEnergy
MyVolumeFE feEnergy
calculate kinetic energy
Definition: ConvectiveMassElement.hpp:87
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1316
ConvectiveMassElement::CommonData::spatialVelocities
string spatialVelocities
Definition: ConvectiveMassElement.hpp:137
ConvectiveMassElement::setConvectiveMassOperators
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
Definition: ConvectiveMassElement.cpp:2005
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
ConvectiveMassElement::CommonData::dataAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
Definition: ConvectiveMassElement.hpp:133
ConvectiveMassElement::CommonData::jacTRowPtr
std::vector< std::vector< double * > > jacTRowPtr
Definition: ConvectiveMassElement.hpp:145
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double *, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ConvectiveMassElement::OpMassLhs_dM_dX
Definition: ConvectiveMassElement.hpp:253
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dX::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:1669
ConvectiveMassElement.hpp
Operators and data structures for mass and convective mass element.
ConvectiveMassElement::addVelocityElement
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
Definition: ConvectiveMassElement.cpp:1892
MoFEM::invertTensor3by3
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.
Definition: Templates.hpp:1566
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
ConvectiveMassElement
structure grouping operators and data used for calculation of mass (convective) element
Definition: ConvectiveMassElement.hpp:28
ConvectiveMassElement::OpVelocityLhs_dV_dx::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:1098
ConvectiveMassElement::MatShellCtx::MatShellCtx
MatShellCtx()
Definition: ConvectiveMassElement.cpp:2343
ConvectiveMassElement::CommonData::jacMassRowPtr
std::vector< std::vector< double * > > jacMassRowPtr
Definition: ConvectiveMassElement.hpp:142
ConvectiveMassElement::PCShellCtx::iNit
MoFEMErrorCode iNit()
Definition: ConvectiveMassElement.cpp:2383
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
ConvectiveMassElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition: ConvectiveMassElement.hpp:134
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEM::CoreInterface::problem_basic_method_preProcess
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
ConvectiveMassElement::OpGetDataAtGaussPts::OpGetDataAtGaussPts
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
Definition: ConvectiveMassElement.cpp:81
ConvectiveMassElement::MatShellCtx::~MatShellCtx
virtual ~MatShellCtx()
Definition: ConvectiveMassElement.cpp:2344
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
ConvectiveMassElement::addConvectiveMassElement
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
Definition: ConvectiveMassElement.cpp:1835
ConvectiveMassElement::OpVelocityRhs::OpVelocityRhs
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:1007
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
ConvectiveMassElement::OpVelocityLhs_dV_dX
Definition: ConvectiveMassElement.hpp:337
ConvectiveMassElement::OpMassLhs_dM_dx
Definition: ConvectiveMassElement.hpp:245
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
ConvectiveMassElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: ConvectiveMassElement.hpp:122
ConvectiveMassElement::OpEnergy::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:737
ConvectiveMassElement::setKinematicEshelbyOperators
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
Definition: ConvectiveMassElement.cpp:2181
ConvectiveMassElement::feVelRhs
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
Definition: ConvectiveMassElement.hpp:77
a
constexpr double a
Definition: approx_sphere.cpp:30
ConvectiveMassElement::OpGetCommonDataAtGaussPts
Definition: ConvectiveMassElement.hpp:170
ConvectiveMassElement::UpdateAndControl::preProcess
MoFEMErrorCode preProcess()
Scatter values from t_u_dt on the fields.
Definition: ConvectiveMassElement.cpp:1722
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
ConvectiveMassElement::OpVelocityLhs_dV_dv::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:1064
BODYFORCESSET
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
ConvectiveMassElement::ShellResidualElement::preProcess
MoFEMErrorCode preProcess()
Calculate inconsistency between approximation of velocities and velocities calculated from displaceme...
Definition: ConvectiveMassElement.cpp:2409
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
ConvectiveMassElement::OpMassLhs_dM_dX::getJac
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:669
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumRhs::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:1466
ConvectiveMassElement::OpEnergy::OpEnergy
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > v)
Definition: ConvectiveMassElement.cpp:727
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
ConvectiveMassElement::OpGetDataAtGaussPts::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
Definition: ConvectiveMassElement.cpp:90
ConvectiveMassElement::OpVelocityLhs_dV_dx
Definition: ConvectiveMassElement.hpp:327
h
double h
Definition: photon_diffusion.cpp:60
COL
@ COL
Definition: definitions.h:123
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
ConvectiveMassElement::OpMassLhs_dM_dv::forcesOnlyOnEntities
Range forcesOnlyOnEntities
Definition: ConvectiveMassElement.hpp:228
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
ConvectiveMassElement::feMassAuxLhs
MyVolumeFE feMassAuxLhs
Definition: ConvectiveMassElement.hpp:71
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dv::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:1544
ConvectiveMassElement::OpMassLhs_dM_dv::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ConvectiveMassElement.cpp:521
ConvectiveMassElement::OpVelocityJacobian::OpVelocityJacobian
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
Definition: ConvectiveMassElement.cpp:788
ConvectiveMassElement::OpMassLhs_dM_dX::OpMassLhs_dM_dX
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:664
Projection10NodeCoordsOnField.hpp
ConvectiveMassElement::OpVelocityJacobian
Definition: ConvectiveMassElement.hpp:279
ConvectiveMassElement::CommonData::jacMass
std::vector< MatrixDouble > jacMass
Definition: ConvectiveMassElement.hpp:143
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition: ConvectiveMassElement.hpp:116
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::OpEshelbyDynamicMaterialMomentumJacobian
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
Definition: ConvectiveMassElement.cpp:1234
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:155
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
ConvectiveMassElement::OpVelocityLhs_dV_dX::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:1173
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
ConvectiveMassElement::CommonData::jacVel
std::vector< MatrixDouble > jacVel
Definition: ConvectiveMassElement.hpp:140
ConvectiveMassElement::CommonData::spatialPositions
string spatialPositions
Definition: ConvectiveMassElement.hpp:135
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ConvectiveMassElement::OpMassLhs_dM_dv::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:452
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumRhs::OpEshelbyDynamicMaterialMomentumRhs
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
Definition: ConvectiveMassElement.cpp:1453
ConvectiveMassElement::OpMassLhs_dM_dx::getJac
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:617
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dx::OpEshelbyDynamicMaterialMomentumLhs_dx
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
Definition: ConvectiveMassElement.cpp:1603
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dX::OpEshelbyDynamicMaterialMomentumLhs_dX
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
Definition: ConvectiveMassElement.cpp:1660
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
ConvectiveMassElement::OpMassJacobian::OpMassJacobian
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
Definition: ConvectiveMassElement.cpp:154
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dX
Definition: ConvectiveMassElement.hpp:414
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumRhs
Definition: ConvectiveMassElement.hpp:370
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
ConvectiveMassElement::CommonData::lInear
bool lInear
Definition: ConvectiveMassElement.hpp:129
FTensor::Index< 'i', 3 >
ConvectiveMassElement::CommonData::valT
std::vector< VectorDouble > valT
Definition: ConvectiveMassElement.hpp:144
ConvectiveMassElement::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: ConvectiveMassElement.hpp:150
ConvectiveMassElement::MyVolumeFE::getRule
int getRule(int order)
it is used to calculate nb. of Gauss integration points
Definition: ConvectiveMassElement.cpp:38
ConvectiveMassElement::UpdateAndControl::postProcess
MoFEMErrorCode postProcess()
Definition: ConvectiveMassElement.cpp:1752
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::getFTensor2FromArray3by3
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1320
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
MoFEM::Block_BodyForces::data
_data_ data
Definition: MaterialBlocks.hpp:328
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:1244
ConvectiveMassElement::MyVolumeFE::V
SmartPetscObj< Vec > V
Definition: ConvectiveMassElement.hpp:55
Range
ConvectiveMassElement::setBlocks
MoFEMErrorCode setBlocks()
Definition: ConvectiveMassElement.cpp:1761
ConvectiveMassElement::commonData
CommonData commonData
Definition: ConvectiveMassElement.hpp:148
adouble
FTensor::dd
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)
Definition: ddTensor0.hpp:33
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
ConvectiveMassElement::OpVelocityLhs_dV_dX::OpVelocityLhs_dV_dX
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:1168
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
ConvectiveMassElement::feMassLhs
MyVolumeFE feMassLhs
Definition: ConvectiveMassElement.hpp:66
MoFEM::EntitiesFieldData::EntData::getFTensor1FieldData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coefficients.
Definition: EntitiesFieldData.hpp:1275
FTensor::Tensor0
Definition: Tensor0.hpp:16
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dx
Definition: ConvectiveMassElement.hpp:401
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
ConvectiveMassElement::ShellResidualElement::ShellResidualElement
ShellResidualElement(MoFEM::Interface &m_field)
Definition: ConvectiveMassElement.cpp:2405
ConvectiveMassElement::CommonData
common data used by volume elements
Definition: ConvectiveMassElement.hpp:127
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dv::OpEshelbyDynamicMaterialMomentumLhs_dv
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
Definition: ConvectiveMassElement.cpp:1535
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumLhs_dx::getJac
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Definition: ConvectiveMassElement.cpp:1612
G
constexpr IntegrationType G
Definition: level_set.cpp:33
ConvectiveMassElement::OpMassLhs_dM_dv::OpMassLhs_dM_dv
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
Definition: ConvectiveMassElement.cpp:439
MoFEM::TSMethod
data structure for TS (time stepping) context
Definition: LoopMethods.hpp:134
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
ConvectiveMassElement::OpMassJacobian::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:163
ConvectiveMassElement::OpGetCommonDataAtGaussPts::OpGetCommonDataAtGaussPts
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:149
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ConvectiveMassElement::CommonData::valMass
std::vector< VectorDouble > valMass
Definition: ConvectiveMassElement.hpp:141
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
ConvectiveMassElement::OpVelocityRhs::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: ConvectiveMassElement.cpp:1013
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
a0
constexpr double a0
Definition: hcurl_check_approx_in_2d.cpp:37
ConvectiveMassElement::feTLhs
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
Definition: ConvectiveMassElement.hpp:84
ConvectiveMassElement::OpMassJacobian
Definition: ConvectiveMassElement.hpp:177
ConvectiveMassElement::CommonData::meshPositions
string meshPositions
Definition: ConvectiveMassElement.hpp:136
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
ConvectiveMassElement::MyVolumeFE::postProcess
MoFEMErrorCode postProcess()
Definition: ConvectiveMassElement.cpp:56
ConvectiveMassElement::OpVelocityLhs_dV_dv::OpVelocityLhs_dV_dv
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: ConvectiveMassElement.cpp:1059
ConvectiveMassElement::feMassRhs
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
Definition: ConvectiveMassElement.hpp:62
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
ConvectiveMassElement::ConvectiveMassElement
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
Definition: ConvectiveMassElement.cpp:75
ConvectiveMassElement::MyVolumeFE::preProcess
MoFEMErrorCode preProcess()
Definition: ConvectiveMassElement.cpp:40
MoFEM::SmartPetscObj< Vec >
ConvectiveMassElement::MyVolumeFE::MyVolumeFE
MyVolumeFE(MoFEM::Interface &m_field)
Definition: ConvectiveMassElement.cpp:22
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
ConvectiveMassElement::mField
MoFEM::Interface & mField
Definition: ConvectiveMassElement.hpp:108
ConvectiveMassElement::feTRhs
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
Definition: ConvectiveMassElement.hpp:82
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
DirichletBC.hpp
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
ConvectiveMassElement::MatShellCtx::dEstroy
MoFEMErrorCode dEstroy()
Definition: ConvectiveMassElement.cpp:2369
H
double H
Hardening.
Definition: plastic.cpp:124
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
ConvectiveMassElement::PCShellCtx::dEstroy
MoFEMErrorCode dEstroy()
Definition: ConvectiveMassElement.cpp:2395
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian
Definition: ConvectiveMassElement.hpp:347
F
@ F
Definition: free_surface.cpp:394
ConvectiveMassElement::tAg
short int tAg
Definition: ConvectiveMassElement.hpp:109
ConvectiveMassElement::CommonData::valVel
std::vector< VectorDouble > valVel
Definition: ConvectiveMassElement.hpp:138
ConvectiveMassElement::OpGetDataAtGaussPts
Definition: ConvectiveMassElement.hpp:152
ConvectiveMassElement::UpdateAndControl::UpdateAndControl
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
Definition: ConvectiveMassElement.cpp:1716