v0.14.0
NonLinearElasticElement.cpp
Go to the documentation of this file.
1 /**
2  * \brief Operators and data structures for nonlinear elastic material
3  *
4  * Implementation of nonlinear elastic element.
5  */
6 
7 
8 
9 #include <MoFEM.hpp>
10 
11 using namespace MoFEM;
13 
14 #include <adolc/adolc.h>
16 
18  : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULL), F(PETSC_NULL),
19  addToRule(1) {
20 
21  auto create_vec = [&]() {
22  constexpr int ghosts[] = {0};
23  if (mField.get_comm_rank() == 0) {
24  return createVectorMPI(mField.get_comm(), 1, 1);
25  } else {
26  return createVectorMPI(mField.get_comm(), 0, 1);
27  }
28  };
29 
30  V = create_vec();
31 }
32 
34  return 2 * (order - 1) + addToRule;
35 };
36 
39 
40  CHKERR VolumeElementForcesAndSourcesCore::preProcess();
41 
42  if (A != PETSC_NULL) {
43  snes_B = A;
44  }
45 
46  if (F != PETSC_NULL) {
47  snes_f = F;
48  }
49 
50  switch (snes_ctx) {
51  case CTX_SNESNONE:
52  CHKERR VecZeroEntries(V);
53  break;
54  default:
55  break;
56  }
57 
59 }
60 
63 
64  const double *array;
65 
66  switch (snes_ctx) {
67  case CTX_SNESNONE:
68  CHKERR VecAssemblyBegin(V);
69  CHKERR VecAssemblyEnd(V);
70  CHKERR VecSum(V, &eNergy);
71  break;
72  default:
73  break;
74  }
75 
76  CHKERR VolumeElementForcesAndSourcesCore::postProcess();
77 
79 }
80 
82  short int tag)
83  : feRhs(m_field), feLhs(m_field), feEnergy(m_field), mField(m_field),
84  tAg(tag) {}
85 
87  const std::string field_name,
88  std::vector<VectorDouble> &values_at_gauss_pts,
89  std::vector<MatrixDouble> &gardient_at_gauss_pts)
92  valuesAtGaussPts(values_at_gauss_pts),
93  gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
94 
96  int side, EntityType type, EntitiesFieldData::EntData &data) {
98 
99  const int nb_dofs = data.getFieldData().size();
100  const int nb_base_functions = data.getN().size2();
101  if (nb_dofs == 0) {
103  }
104  const int nb_gauss_pts = data.getN().size1();
105  const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
106 
107  // initialize
108  if (type == zeroAtType) {
109  valuesAtGaussPts.resize(nb_gauss_pts);
110  gradientAtGaussPts.resize(nb_gauss_pts);
111  for (int gg = 0; gg != nb_gauss_pts; gg++) {
112  valuesAtGaussPts[gg].resize(rank, false);
113  gradientAtGaussPts[gg].resize(rank, 3, false);
114  }
115  for (int gg = 0; gg != nb_gauss_pts; gg++) {
116  valuesAtGaussPts[gg].clear();
117  gradientAtGaussPts[gg].clear();
118  }
119  }
120 
121  auto base_function = data.getFTensor0N();
122  auto diff_base_functions = data.getFTensor1DiffN<3>();
125 
126  if (rank == 1) {
127 
128  for (int gg = 0; gg != nb_gauss_pts; gg++) {
129  auto field_data = data.getFTensor0FieldData();
130  double &val = valuesAtGaussPts[gg][0];
131  FTensor::Tensor1<double *, 3> grad(&gradientAtGaussPts[gg](0, 0),
132  &gradientAtGaussPts[gg](0, 1),
133  &gradientAtGaussPts[gg](0, 2));
134  int bb = 0;
135  for (; bb != nb_dofs; bb++) {
136  val += base_function * field_data;
137  grad(i) += diff_base_functions(i) * field_data;
138  ++diff_base_functions;
139  ++base_function;
140  ++field_data;
141  }
142  for (; bb != nb_base_functions; bb++) {
143  ++diff_base_functions;
144  ++base_function;
145  }
146  }
147 
148  } else if (rank == 3) {
149 
150  for (int gg = 0; gg != nb_gauss_pts; gg++) {
151  auto field_data = data.getFTensor1FieldData<3>();
152  FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
153  &valuesAtGaussPts[gg][1],
154  &valuesAtGaussPts[gg][2]);
156  &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
157  &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
158  &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
159  &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
160  &gradientAtGaussPts[gg](2, 2));
161  int bb = 0;
162  for (; bb != nb_dofs / 3; bb++) {
163  values(i) += base_function * field_data(i);
164  gradient(i, j) += field_data(i) * diff_base_functions(j);
165  ++diff_base_functions;
166  ++base_function;
167  ++field_data;
168  }
169  for (; bb != nb_base_functions; bb++) {
170  ++diff_base_functions;
171  ++base_function;
172  }
173  }
174 
175  } else {
176  // FIXME: THat part is inefficient
177  VectorDouble &values = data.getFieldData();
178  for (int gg = 0; gg < nb_gauss_pts; gg++) {
179  VectorAdaptor N = data.getN(gg, nb_dofs / rank);
180  MatrixAdaptor diffN = data.getDiffN(gg, nb_dofs / rank);
181  for (int dd = 0; dd < nb_dofs / rank; dd++) {
182  for (int rr1 = 0; rr1 < rank; rr1++) {
183  valuesAtGaussPts[gg][rr1] += N[dd] * values[rank * dd + rr1];
184  for (int rr2 = 0; rr2 < 3; rr2++) {
185  gradientAtGaussPts[gg](rr1, rr2) +=
186  diffN(dd, rr2) * values[rank * dd + rr1];
187  }
188  }
189  }
190  }
191  }
192 
194 }
195 
197  const std::string field_name, CommonData &common_data)
198  : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
199  common_data.gradAtGaussPts[field_name]) {}
200 
203  BlockData &data, CommonData &common_data,
204  int tag, bool jacobian, bool ale,
205  bool field_disp)
207  field_name, UserDataOperator::OPROW),
208  dAta(data), commonData(common_data), tAg(tag), adlocReturnValue(0),
209  jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
210 
211 }
212 
215  const int gg) {
217 
218  CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
219  dAta, getNumeredEntFiniteElementPtr());
220 
221  if (aLe) {
222  auto &t_P = dAta.materialAdoublePtr->t_P;
223  auto &t_invH = dAta.materialAdoublePtr->t_invH;
224  t_P(i, j) = t_P(i, k) * t_invH(j, k);
225  t_P(i, j) *= dAta.materialAdoublePtr->detH;
226  }
227 
228  commonData.sTress[gg].resize(3, 3, false);
229  for (int dd1 = 0; dd1 < 3; dd1++) {
230  for (int dd2 = 0; dd2 < 3; dd2++) {
231  dAta.materialAdoublePtr->P(dd1, dd2) >>=
232  (commonData.sTress[gg])(dd1, dd2);
233  }
234  }
235 
237 }
238 
241  const int gg) {
243 
244  trace_on(tAg, 0);
245 
246  dAta.materialAdoublePtr->F.resize(3, 3, false);
247 
248  if (!aLe) {
249 
250  nbActiveVariables = 0;
251  for (int dd1 = 0; dd1 < 3; dd1++) {
252  for (int dd2 = 0; dd2 < 3; dd2++) {
253  dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
254  if (fieldDisp) {
255  if (dd1 == dd2) {
256  dAta.materialAdoublePtr->F(dd1, dd2) += 1;
257  }
258  }
259  nbActiveVariables++;
260  }
261  }
262 
263  } else {
264 
265  nbActiveVariables = 0;
266 
267  dAta.materialAdoublePtr->h.resize(3, 3, false);
268  for (int dd1 = 0; dd1 < 3; dd1++) {
269  for (int dd2 = 0; dd2 < 3; dd2++) {
270  dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
271  nbActiveVariables++;
272  }
273  }
274 
275  dAta.materialAdoublePtr->H.resize(3, 3, false);
276  for (int dd1 = 0; dd1 < 3; dd1++) {
277  for (int dd2 = 0; dd2 < 3; dd2++) {
278  dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
279  nbActiveVariables++;
280  }
281  }
282 
283  dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
284  dAta.materialAdoublePtr->invH.resize(3, 3, false);
285  CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
286  dAta.materialAdoublePtr->detH,
287  dAta.materialAdoublePtr->invH);
288 
289  auto &t_F = dAta.materialAdoublePtr->t_F;
290  auto &t_h = dAta.materialAdoublePtr->t_h;
291  auto &t_invH = dAta.materialAdoublePtr->t_invH;
292 
293  t_F(i, j) = t_h(i, k) * t_invH(k, j);
294 
295  }
296 
297  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
298  CHKERR calculateStress(gg);
299 
300  trace_off();
301 
303 }
304 
308 
309  int r;
310 
311  if (fUnction) {
312  commonData.sTress[gg].resize(3, 3, false);
313  // play recorder for values
314  r = ::function(tAg, 9, nbActiveVariables, &activeVariables[0],
315  &commonData.sTress[gg](0, 0));
316  if (r < adlocReturnValue) { // function is locally analytic
317  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
318  "ADOL-C function evaluation with error r = %d", r);
319  }
320  }
321 
322  if (jAcobian) {
323  commonData.jacStress[gg].resize(9, nbActiveVariables, false);
324  double *jac_ptr[] = {
325  &(commonData.jacStress[gg](0, 0)), &(commonData.jacStress[gg](1, 0)),
326  &(commonData.jacStress[gg](2, 0)), &(commonData.jacStress[gg](3, 0)),
327  &(commonData.jacStress[gg](4, 0)), &(commonData.jacStress[gg](5, 0)),
328  &(commonData.jacStress[gg](6, 0)), &(commonData.jacStress[gg](7, 0)),
329  &(commonData.jacStress[gg](8, 0))};
330  // play recorder for jacobians
331  r = jacobian(tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
332  if (r < adlocReturnValue) {
333  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
334  "ADOL-C function evaluation with error");
335  }
336  }
337 
339 }
340 
342  int row_side, EntityType row_type,
343  EntitiesFieldData::EntData &row_data) {
345 
346  // do it only once, no need to repeat this for edges,faces or tets
347  if (row_type != MBVERTEX)
349 
350  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
351  dAta.tEts.end()) {
353  }
354 
355  int nb_dofs = row_data.getFieldData().size();
356  if (nb_dofs == 0)
358  dAta.materialAdoublePtr->commonDataPtr = &commonData;
359  dAta.materialAdoublePtr->opPtr = this;
360 
361  int nb_gauss_pts = row_data.getN().size1();
362  commonData.sTress.resize(nb_gauss_pts);
363  commonData.jacStress.resize(nb_gauss_pts);
364 
366  if (aLe) {
368  }
369 
370  for (int gg = 0; gg != nb_gauss_pts; gg++) {
371 
372  dAta.materialAdoublePtr->gG = gg;
373 
374  // Record tag and calculate stress
375  if (recordTagForIntegrationPoint(gg)) {
376  CHKERR recordTag(gg);
377  }
378 
379  // Set active variables vector
380  if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
381  activeVariables.resize(nbActiveVariables, false);
382  if (!aLe) {
383  for (int dd1 = 0; dd1 < 3; dd1++) {
384  for (int dd2 = 0; dd2 < 3; dd2++) {
385  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
386  }
387  }
388  } else {
389  for (int dd1 = 0; dd1 < 3; dd1++) {
390  for (int dd2 = 0; dd2 < 3; dd2++) {
391  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
392  }
393  }
394  for (int dd1 = 0; dd1 < 3; dd1++) {
395  for (int dd2 = 0; dd2 < 3; dd2++) {
396  activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
397  }
398  }
399  }
400  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
401 
402  // Play tag and calculate stress or tangent
403  if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
404  CHKERR playTag(gg);
405  }
406  }
407  }
408 
410 }
411 
413  const std::string
414  field_name, ///< field name for spatial positions or displacements
415  BlockData &data, CommonData &common_data, int tag, bool gradient,
416  bool hessian, bool ale, bool field_disp)
418  field_name, UserDataOperator::OPROW),
419  dAta(data), commonData(common_data), tAg(tag), gRadient(gradient),
420  hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
421 
425  CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
426  dAta, getNumeredEntFiniteElementPtr());
427  dAta.materialAdoublePtr->eNergy >>= commonData.eNergy[gg];
429 }
430 
434 
435  trace_on(tAg, 0);
436 
437  if (!aLe) {
438 
439  nbActiveVariables = 0;
440  for (int dd1 = 0; dd1 < 3; dd1++) {
441  for (int dd2 = 0; dd2 < 3; dd2++) {
442  dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
443  if (fieldDisp) {
444  if (dd1 == dd2) {
445  dAta.materialAdoublePtr->F(dd1, dd2) += 1;
446  }
447  }
448  nbActiveVariables++;
449  }
450  }
451 
452  } else {
453 
454  nbActiveVariables = 0;
455 
456  dAta.materialAdoublePtr->h.resize(3, 3, false);
457  for (int dd1 = 0; dd1 < 3; dd1++) {
458  for (int dd2 = 0; dd2 < 3; dd2++) {
459  dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
460  nbActiveVariables++;
461  }
462  }
463 
464  dAta.materialAdoublePtr->H.resize(3, 3, false);
465  for (int dd1 = 0; dd1 < 3; dd1++) {
466  for (int dd2 = 0; dd2 < 3; dd2++) {
467  dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
468  nbActiveVariables++;
469  }
470  }
471 
472  dAta.materialAdoublePtr->detH = determinantTensor3by3(dAta.materialAdoublePtr->H);
473  dAta.materialAdoublePtr->invH.resize(3, 3, false);
474  CHKERR invertTensor3by3(dAta.materialAdoublePtr->H,
475  dAta.materialAdoublePtr->detH,
476  dAta.materialAdoublePtr->invH);
477 
478  auto &t_F = dAta.materialAdoublePtr->t_F;
479  auto &t_h = dAta.materialAdoublePtr->t_h;
480  auto &t_invH = dAta.materialAdoublePtr->t_invH;
481 
482  t_F(i, j) = t_h(i, k) * t_invH(k, j);
483 
484  }
485 
486  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
488 
489  trace_off();
490 
492 }
493 
497 
498  if (gRadient) {
499  commonData.jacEnergy[gg].resize(nbActiveVariables, false);
500  int r = ::gradient(tAg, nbActiveVariables, &activeVariables[0],
501  &commonData.jacEnergy[gg][0]);
502  if (r < 0) {
503  // That means that energy function is not smooth and derivative
504  // can not be calculated,
505  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
506  "ADOL-C function evaluation with error");
507  }
508  }
509 
510  if (hEssian) {
511  commonData.hessianEnergy[gg].resize(nbActiveVariables * nbActiveVariables,
512  false);
513  double *H[nbActiveVariables];
514  for (int n = 0; n != nbActiveVariables; n++) {
515  H[n] = &(commonData.hessianEnergy[gg][n * nbActiveVariables]);
516  }
517  int r = ::hessian(tAg, nbActiveVariables, &*activeVariables.begin(), H);
518  if (r < 0) {
519  // That means that energy function is not smooth and derivative
520  // can not be calculated,
521  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
522  "ADOL-C function evaluation with error");
523  }
524  }
525 
527 }
528 
530  int row_side, EntityType row_type,
531  EntitiesFieldData::EntData &row_data) {
533 
534  // do it only once, no need to repeat this for edges,faces or tets
535  if (row_type != MBVERTEX)
537 
538  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
539  dAta.tEts.end()) {
541  }
542 
543  int nb_dofs = row_data.getFieldData().size();
544  if (nb_dofs == 0)
546  dAta.materialAdoublePtr->commonDataPtr = &commonData;
547  dAta.materialAdoublePtr->opPtr = this;
548 
549  int nb_gauss_pts = row_data.getN().size1();
550  commonData.eNergy.resize(nb_gauss_pts);
551  commonData.jacEnergy.resize(nb_gauss_pts);
552 
554  if (aLe) {
556  }
557 
558  for (int gg = 0; gg != nb_gauss_pts; gg++) {
559 
560  dAta.materialAdoublePtr->gG = gg;
561 
562  // Record tag and calualte stress
563  if (recordTagForIntegrationPoint(gg)) {
564  CHKERR recordTag(gg);
565  }
566 
567  activeVariables.resize(nbActiveVariables, false);
568  if (!aLe) {
569  for (int dd1 = 0; dd1 < 3; dd1++) {
570  for (int dd2 = 0; dd2 < 3; dd2++) {
571  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
572  }
573  }
574  } else {
575  for (int dd1 = 0; dd1 < 3; dd1++) {
576  for (int dd2 = 0; dd2 < 3; dd2++) {
577  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
578  }
579  }
580  for (int dd1 = 0; dd1 < 3; dd1++) {
581  for (int dd2 = 0; dd2 < 3; dd2++) {
582  activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
583  }
584  }
585  }
586  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
587 
588  // Play tag and calculate stress or tangent
589  CHKERR playTag(gg);
590  }
591 
593 }
594 
596  const std::string field_name, BlockData &data, CommonData &common_data)
598  field_name, UserDataOperator::OPROW),
599  dAta(data), commonData(common_data), aLe(false) {}
600 
602  int row_side, EntityType row_type,
603  EntitiesFieldData::EntData &row_data) {
605 
606  int nb_dofs = row_data.getIndices().size();
607  int *indices_ptr = &row_data.getIndices()[0];
608  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
609  iNdices.resize(nb_dofs, false);
610  noalias(iNdices) = row_data.getIndices();
611  indices_ptr = &iNdices[0];
612  VectorDofs &dofs = row_data.getFieldDofs();
613  VectorDofs::iterator dit = dofs.begin();
614  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
615  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
616  dAta.forcesOnlyOnEntitiesRow.end()) {
617  iNdices[ii] = -1;
618  }
619  }
620  }
621  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
622  ADD_VALUES);
624 }
625 
627  int row_side, EntityType row_type,
628  EntitiesFieldData::EntData &row_data) {
630 
631  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
632  dAta.tEts.end()) {
634  }
635 
636  const int nb_dofs = row_data.getIndices().size();
637  if (nb_dofs == 0)
639  if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
640  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
641  }
642  const int nb_base_functions = row_data.getN().size2();
643  const int nb_gauss_pts = row_data.getN().size1();
644 
645  nf.resize(nb_dofs, false);
646  nf.clear();
647 
648  auto diff_base_functions = row_data.getFTensor1DiffN<3>();
651 
652  for (int gg = 0; gg != nb_gauss_pts; gg++) {
653  double val = getVolume() * getGaussPts()(3, gg);
654  MatrixDouble3by3 &stress = commonData.sTress[gg];
656  &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
657  &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
658  &stress(2, 2));
659  FTensor::Tensor1<double *, 3> rhs(&nf[0], &nf[1], &nf[2], 3);
660  int bb = 0;
661  for (; bb != nb_dofs / 3; bb++) {
662  rhs(i) += val * t3(i, j) * diff_base_functions(j);
663  ++rhs;
664  ++diff_base_functions;
665  }
666  for (; bb != nb_base_functions; bb++) {
667  ++diff_base_functions;
668  }
669  }
670 
671  CHKERR aSemble(row_side, row_type, row_data);
672 
674 }
675 
677  BlockData &data,
678  CommonData &common_data,
679  SmartPetscObj<Vec> ghost_vec,
680  bool field_disp)
682  field_name, UserDataOperator::OPROW),
683  dAta(data), commonData(common_data), ghostVec(ghost_vec, true),
684  fieldDisp(field_disp) {}
685 
687  int row_side, EntityType row_type,
688  EntitiesFieldData::EntData &row_data) {
690 
691  if (row_type != MBVERTEX)
693  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
694  dAta.tEts.end()) {
696  }
697 
698  std::vector<MatrixDouble> &F =
700  dAta.materialDoublePtr->F.resize(3, 3, false);
701 
702  double energy = 0;
703 
704  for (unsigned int gg = 0; gg != row_data.getN().size1(); ++gg) {
705  double val = getVolume() * getGaussPts()(3, gg);
706  noalias(dAta.materialDoublePtr->F) = F[gg];
707  if (fieldDisp) {
708  for (int dd = 0; dd < 3; dd++) {
709  dAta.materialDoublePtr->F(dd, dd) += 1;
710  }
711  }
712  int nb_active_variables = 0;
713  CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
714  CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
715  dAta, getNumeredEntFiniteElementPtr());
716  energy += val * dAta.materialDoublePtr->eNergy;
717  }
718 
719  CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
721 }
722 
724  const std::string vel_field, const std::string field_name, BlockData &data,
725  CommonData &common_data)
727  vel_field, field_name, UserDataOperator::OPROWCOL),
728  dAta(data), commonData(common_data), aLe(false) {}
729 
730 template <int S>
732  int gg, MatrixDouble &jac_stress,
733  MatrixDouble &jac) {
734  jac.clear();
738  int nb_col = col_data.getFieldData().size();
739  double *diff_ptr =
740  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
741  // First two indices 'i','j' derivatives of 1st Piola-stress, third index 'k'
742  // is displacement component
744  &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
745  &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
746  &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
747  &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
748  &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
749  &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
750  &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
751  &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
752  &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
753  &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
754  &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
755  &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
756  &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
757  &jac_stress(3 * 2 + 2, S + 2));
759  &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
760  &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
761  &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
762  &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
763  &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
764  &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
765  &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
766  &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
767  &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
768  &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
769  &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
770  &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
771  &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
772  &jac_stress(3 * 2 + 2, S + 5));
774  &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
775  &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
776  &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
777  &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
778  &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
779  &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
780  &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
781  &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
782  &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
783  &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
784  &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
785  &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
786  &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
787  &jac_stress(3 * 2 + 2, S + 8));
788  // Derivate of 1st Piola-stress multiplied by gradient of defamation for
789  // base function (dd) and displacement component (rr)
791  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
792  &jac(6, 0), &jac(7, 0), &jac(8, 0));
794  &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
795  &jac(6, 1), &jac(7, 1), &jac(8, 1));
797  &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
798  &jac(6, 2), &jac(7, 2), &jac(8, 2));
800  diff_ptr, &diff_ptr[1], &diff_ptr[2]);
801  for (int dd = 0; dd != nb_col / 3; ++dd) {
802  t2_1_0(i, j) += t3_1_0(i, j, k) * diff(k);
803  t2_1_1(i, j) += t3_1_1(i, j, k) * diff(k);
804  t2_1_2(i, j) += t3_1_2(i, j, k) * diff(k);
805  ++t2_1_0;
806  ++t2_1_1;
807  ++t2_1_2;
808  ++diff;
809  }
811 }
812 
814  EntitiesFieldData::EntData &col_data, int gg) {
815  return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
816 }
817 
819  int row_side, int col_side, EntityType row_type, EntityType col_type,
820  EntitiesFieldData::EntData &row_data,
821  EntitiesFieldData::EntData &col_data) {
823 
824  int nb_row = row_data.getIndices().size();
825  int nb_col = col_data.getIndices().size();
826 
827  int *row_indices_ptr = &row_data.getIndices()[0];
828  int *col_indices_ptr = &col_data.getIndices()[0];
829 
830  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
831  rowIndices.resize(nb_row, false);
832  noalias(rowIndices) = row_data.getIndices();
833  row_indices_ptr = &rowIndices[0];
834  VectorDofs &dofs = row_data.getFieldDofs();
835  VectorDofs::iterator dit = dofs.begin();
836  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
837  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
838  dAta.forcesOnlyOnEntitiesRow.end()) {
839  rowIndices[ii] = -1;
840  }
841  }
842  }
843 
844  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
845  colIndices.resize(nb_col, false);
846  noalias(colIndices) = col_data.getIndices();
847  col_indices_ptr = &colIndices[0];
848  VectorDofs &dofs = col_data.getFieldDofs();
849  VectorDofs::iterator dit = dofs.begin();
850  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
851  if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
852  dAta.forcesOnlyOnEntitiesCol.end()) {
853  colIndices[ii] = -1;
854  }
855  }
856  }
857 
858  CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
859  col_indices_ptr, &k(0, 0), ADD_VALUES);
860 
861  // is symmetric
862  if (row_side != col_side || row_type != col_type) {
863 
864  row_indices_ptr = &row_data.getIndices()[0];
865  col_indices_ptr = &col_data.getIndices()[0];
866 
867  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
868  rowIndices.resize(nb_row, false);
869  noalias(rowIndices) = row_data.getIndices();
870  row_indices_ptr = &rowIndices[0];
871  VectorDofs &dofs = row_data.getFieldDofs();
872  VectorDofs::iterator dit = dofs.begin();
873  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
874  if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
875  dAta.forcesOnlyOnEntitiesCol.end()) {
876  rowIndices[ii] = -1;
877  }
878  }
879  }
880 
881  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
882  colIndices.resize(nb_col, false);
883  noalias(colIndices) = col_data.getIndices();
884  col_indices_ptr = &colIndices[0];
885  VectorDofs &dofs = col_data.getFieldDofs();
886  VectorDofs::iterator dit = dofs.begin();
887  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
888  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
889  dAta.forcesOnlyOnEntitiesRow.end()) {
890  colIndices[ii] = -1;
891  }
892  }
893  }
894 
895  trans_k.resize(nb_col, nb_row, false);
896  noalias(trans_k) = trans(k);
897  CHKERR MatSetValues(getFEMethod()->snes_B, nb_col, col_indices_ptr, nb_row,
898  row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
899  }
900 
902 }
903 
905  int row_side, int col_side, EntityType row_type, EntityType col_type,
906  EntitiesFieldData::EntData &row_data,
907  EntitiesFieldData::EntData &col_data) {
909 
910  int nb_row = row_data.getIndices().size();
911  int nb_col = col_data.getIndices().size();
912  if (nb_row == 0)
914  if (nb_col == 0)
916 
917  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
918  dAta.tEts.end()) {
920  }
921 
922  // const int nb_base_functions = row_data.getN().size2();
923  const int nb_gauss_pts = row_data.getN().size1();
924 
928 
929  k.resize(nb_row, nb_col, false);
930  k.clear();
931  jac.resize(9, nb_col, false);
932 
933  for (int gg = 0; gg != nb_gauss_pts; gg++) {
934  CHKERR getJac(col_data, gg);
935  double val = getVolume() * getGaussPts()(3, gg);
937  &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
938  &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
939  &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
940  &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
941  &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
942  &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
943  &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
944  &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
945  &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
946  for (int cc = 0; cc != nb_col / 3; cc++) {
947  auto diff_base_functions = row_data.getFTensor1DiffN<3>(gg, 0);
949  &k(0, 3 * cc + 0), &k(0, 3 * cc + 1), &k(0, 3 * cc + 2),
950  &k(1, 3 * cc + 0), &k(1, 3 * cc + 1), &k(1, 3 * cc + 2),
951  &k(2, 3 * cc + 0), &k(2, 3 * cc + 1), &k(2, 3 * cc + 2), 3 * nb_col);
952  for (int rr = 0; rr != nb_row / 3; rr++) {
953  lhs(i, j) += val * t3_1(i, m, j) * diff_base_functions(m);
954  ++diff_base_functions;
955  ++lhs;
956  }
957  ++t3_1;
958  }
959  }
960 
961  CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
962 
964 }
965 
967  const std::string vel_field, const std::string field_name, BlockData &data,
968  CommonData &common_data)
969  : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {
970  sYmm = false;
971 }
972 
974  EntitiesFieldData::EntData &col_data, int gg) {
975  return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
976 }
977 
979  int row_side, int col_side, EntityType row_type, EntityType col_type,
980  EntitiesFieldData::EntData &row_data,
981  EntitiesFieldData::EntData &col_data) {
983 
984  int nb_row = row_data.getIndices().size();
985  int nb_col = col_data.getIndices().size();
986 
987  int *row_indices_ptr = &row_data.getIndices()[0];
988  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
989  rowIndices.resize(nb_row, false);
990  noalias(rowIndices) = row_data.getIndices();
991  row_indices_ptr = &rowIndices[0];
992  VectorDofs &dofs = row_data.getFieldDofs();
993  VectorDofs::iterator dit = dofs.begin();
994  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
995  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
996  dAta.forcesOnlyOnEntitiesRow.end()) {
997  rowIndices[ii] = -1;
998  }
999  }
1000  }
1001 
1002  int *col_indices_ptr = &col_data.getIndices()[0];
1003  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
1004  colIndices.resize(nb_col, false);
1005  noalias(colIndices) = col_data.getIndices();
1006  col_indices_ptr = &colIndices[0];
1007  VectorDofs &dofs = col_data.getFieldDofs();
1008  VectorDofs::iterator dit = dofs.begin();
1009  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1010  if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
1011  dAta.forcesOnlyOnEntitiesCol.end()) {
1012  colIndices[ii] = -1;
1013  }
1014  }
1015  }
1016 
1017  /*for(int dd1 = 0;dd1<k.size1();dd1++) {
1018  for(int dd2 = 0;dd2<k.size2();dd2++) {
1019  if(k(dd1,dd2)!=k(dd1,dd2)) {
1020  SETERRQ(PETSC_COMM_SELF,1,"Wrong result");
1021  }
1022  }
1023  }*/
1024 
1025  CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
1026  col_indices_ptr, &k(0, 0), ADD_VALUES);
1027 
1029 }
1030 
1032  const std::string field_name, BlockData &data, CommonData &common_data,
1033  int tag, bool jacobian, bool ale)
1034  : OpJacobianPiolaKirchhoffStress(field_name, data, common_data, tag,
1035  jacobian, ale, false) {}
1036 
1039  const int gg) {
1041 
1042  CHKERR dAta.materialAdoublePtr->calculatesIGma_EshelbyStress(
1043  dAta, getNumeredEntFiniteElementPtr());
1044  if (aLe) {
1045  auto &t_sIGma = dAta.materialAdoublePtr->t_sIGma;
1046  auto &t_invH = dAta.materialAdoublePtr->t_invH;
1047  t_sIGma(i, j) = t_sIGma(i, k) * t_invH(j, k);
1048  t_sIGma(i, j) *= dAta.materialAdoublePtr->detH;
1049 
1050  }
1051  commonData.sTress[gg].resize(3, 3, false);
1052  for (int dd1 = 0; dd1 < 3; dd1++) {
1053  for (int dd2 = 0; dd2 < 3; dd2++) {
1054  dAta.materialAdoublePtr->sIGma(dd1, dd2) >>=
1055  (commonData.sTress[gg])(dd1, dd2);
1056  }
1057  }
1058 
1060 }
1061 
1063  const std::string field_name, BlockData &data, CommonData &common_data)
1064  : OpRhsPiolaKirchhoff(field_name, data, common_data) {}
1065 
1067  const std::string vel_field, const std::string field_name, BlockData &data,
1068  CommonData &common_data)
1069  : OpLhsPiolaKirchhoff_dX(vel_field, field_name, data, common_data) {}
1070 
1072  EntitiesFieldData::EntData &col_data, int gg) {
1073  return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
1074 }
1075 
1077  const std::string vel_field, const std::string field_name, BlockData &data,
1078  CommonData &common_data)
1079  : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {}
1080 
1082  EntitiesFieldData::EntData &col_data, int gg) {
1083  return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
1084 }
1085 
1088  materialDoublePtr,
1090  materialAdoublePtr) {
1092 
1093  if (!materialDoublePtr) {
1095  "Pointer for materialDoublePtr not allocated");
1096  }
1097  if (!materialAdoublePtr) {
1099  "Pointer for materialAdoublePtr not allocated");
1100  }
1101 
1103  mField, BLOCKSET | MAT_ELASTICSET, it)) {
1104  Mat_Elastic mydata;
1105  CHKERR it->getAttributeDataStructure(mydata);
1106  int id = it->getMeshsetId();
1107  EntityHandle meshset = it->getMeshset();
1108  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1109  setOfBlocks[id].tEts, true);
1110  setOfBlocks[id].iD = id;
1111  setOfBlocks[id].E = mydata.data.Young;
1112  setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
1113  setOfBlocks[id].materialDoublePtr = materialDoublePtr;
1114  setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
1115  }
1116 
1118 }
1119 
1121  const std::string element_name,
1122  const std::string spatial_position_field_name,
1123  const std::string material_position_field_name, const bool ale) {
1125 
1126  CHKERR mField.add_finite_element(element_name, MF_ZERO);
1128  element_name, spatial_position_field_name);
1130  element_name, spatial_position_field_name);
1132  element_name, spatial_position_field_name);
1133  if (mField.check_field(material_position_field_name)) {
1134  if (ale) {
1136  element_name, material_position_field_name);
1138  element_name, material_position_field_name);
1139  }
1141  element_name, material_position_field_name);
1142  }
1143 
1144  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1145  for (; sit != setOfBlocks.end(); sit++) {
1146  CHKERR mField.add_ents_to_finite_element_by_type(sit->second.tEts, MBTET,
1147  element_name);
1148  }
1149 
1151 }
1152 
1154  const std::string spatial_position_field_name,
1155  const std::string material_position_field_name, const bool ale,
1156  const bool field_disp) {
1158 
1159  commonData.spatialPositions = spatial_position_field_name;
1160  commonData.meshPositions = material_position_field_name;
1161 
1162  // Rhs
1163  feRhs.getOpPtrVector().push_back(
1164  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1165  if (mField.check_field(material_position_field_name)) {
1167  material_position_field_name, commonData));
1168  }
1169  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1170  for (; sit != setOfBlocks.end(); sit++) {
1172  spatial_position_field_name, sit->second, commonData, tAg, false, ale,
1173  field_disp));
1174  feRhs.getOpPtrVector().push_back(new OpRhsPiolaKirchhoff(
1175  spatial_position_field_name, sit->second, commonData));
1176  }
1177 
1178  // Energy
1179  feEnergy.getOpPtrVector().push_back(
1180  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1181  if (mField.check_field(material_position_field_name)) {
1183  material_position_field_name, commonData));
1184  }
1185  sit = setOfBlocks.begin();
1186  for (; sit != setOfBlocks.end(); sit++) {
1187  feEnergy.getOpPtrVector().push_back(
1188  new OpEnergy(spatial_position_field_name, sit->second, commonData,
1189  feEnergy.V, field_disp));
1190  }
1191 
1192  // Lhs
1193  feLhs.getOpPtrVector().push_back(
1194  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1195  if (mField.check_field(material_position_field_name)) {
1197  material_position_field_name, commonData));
1198  }
1199  sit = setOfBlocks.begin();
1200  for (; sit != setOfBlocks.end(); sit++) {
1202  spatial_position_field_name, sit->second, commonData, tAg, true, ale,
1203  field_disp));
1205  spatial_position_field_name, spatial_position_field_name, sit->second,
1206  commonData));
1207  }
1208 
1210 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
NonlinearElasticElement::CommonData::eNergy
std::vector< double > eNergy
Definition: NonLinearElasticElement.hpp:117
NonlinearElasticElement::OpGetCommonDataAtGaussPts::OpGetCommonDataAtGaussPts
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:196
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
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
NonlinearElasticElement::OpGetDataAtGaussPts::OpGetDataAtGaussPts
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
Definition: NonLinearElasticElement.cpp:86
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
NonlinearElasticElement::mField
MoFEM::Interface & mField
Definition: NonLinearElasticElement.hpp:73
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< double >
NonlinearElasticElement::MyVolumeFE::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: NonLinearElasticElement.cpp:61
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
NonlinearElasticElement::CommonData::jacStress
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
Definition: NonLinearElasticElement.hpp:113
NonlinearElasticElement::OpJacobianEnergy::OpJacobianEnergy
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
Definition: NonLinearElasticElement.cpp:412
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
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
NonlinearElasticElement::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: NonLinearElasticElement.hpp:65
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
NonlinearElasticElement::tAg
short int tAg
Definition: NonLinearElasticElement.hpp:74
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
NonlinearElasticElement::OpLhsPiolaKirchhoff_dX::aSemble
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: NonLinearElasticElement.cpp:978
NonlinearElasticElement::OpJacobianEshelbyStress::OpJacobianEshelbyStress
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
Definition: NonLinearElasticElement.cpp:1031
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
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
NonlinearElasticElement::OpEnergy::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: NonLinearElasticElement.cpp:686
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1316
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: NonLinearElasticElement.cpp:904
NonlinearElasticElement::OpRhsPiolaKirchhoff::OpRhsPiolaKirchhoff
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:595
NonlinearElasticElement::MyVolumeFE::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: NonLinearElasticElement.cpp:37
NonlinearElasticElement::i
FTensor::Index< 'i', 3 > i
Definition: NonLinearElasticElement.hpp:123
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
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
NonlinearElasticElement::j
FTensor::Index< 'j', 3 > j
Definition: NonLinearElasticElement.hpp:124
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
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
NonlinearElasticElement::OpJacobianEshelbyStress::calculateStress
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition: NonLinearElasticElement.cpp:1038
NonlinearElasticElement::CommonData::spatialPositions
string spatialPositions
Definition: NonLinearElasticElement.hpp:109
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
NonlinearElasticElement::feLhs
MyVolumeFE feLhs
Definition: NonLinearElasticElement.hpp:67
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::OpLhsPiolaKirchhoff_dx
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:723
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
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::OpJacobianPiolaKirchhoffStress
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
Definition: NonLinearElasticElement.cpp:202
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
NonlinearElasticElement::CommonData::jacEnergy
std::vector< VectorDouble > jacEnergy
Definition: NonLinearElasticElement.hpp:118
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
NonlinearElasticElement::OpLhsEshelby_dx::OpLhsEshelby_dx
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:1066
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
convert.type
type
Definition: convert.py:64
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
NonlinearElasticElement::MyVolumeFE::getRule
int getRule(int order)
it is used to calculate nb. of Gauss integration points
Definition: NonLinearElasticElement.cpp:33
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::calculateStress
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition: NonLinearElasticElement.cpp:214
NonlinearElasticElement::CommonData::meshPositions
string meshPositions
Definition: NonLinearElasticElement.hpp:110
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
NonlinearElasticElement::MyVolumeFE::V
SmartPetscObj< Vec > V
Definition: NonLinearElasticElement.hpp:58
NonlinearElasticElement::CommonData::hessianEnergy
std::vector< VectorDouble > hessianEnergy
Definition: NonLinearElasticElement.hpp:119
NonlinearElasticElement::OpLhsPiolaKirchhoff_dX::OpLhsPiolaKirchhoff_dX
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:966
Projection10NodeCoordsOnField.hpp
get_jac
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
Definition: NonLinearElasticElement.cpp:731
NonlinearElasticElement::OpGetDataAtGaussPts::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
Definition: NonLinearElasticElement.cpp:95
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:155
NonlinearElasticElement::addElement
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition: NonLinearElasticElement.cpp:1120
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
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
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::playTag
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
Definition: NonLinearElasticElement.cpp:306
NonlinearElasticElement::OpLhsPiolaKirchhoff_dX
Definition: NonLinearElasticElement.hpp:598
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
NonlinearElasticElement::OpJacobianEnergy::playTag
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
Definition: NonLinearElasticElement.cpp:495
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::recordTag
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
Definition: NonLinearElasticElement.cpp:240
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
NonlinearElasticElement::OpRhsPiolaKirchhoff::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: NonLinearElasticElement.cpp:626
NonlinearElasticElement::setOperators
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Definition: NonLinearElasticElement.cpp:1153
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
NonlinearElasticElement::OpEnergy::OpEnergy
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > ghost_vec, bool field_disp)
Definition: NonLinearElasticElement.cpp:676
convert.n
n
Definition: convert.py:82
NonlinearElasticElement::OpJacobianEnergy::recordTag
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
Definition: NonLinearElasticElement.cpp:432
NonlinearElasticElement::commonData
CommonData commonData
Definition: NonLinearElasticElement.hpp:121
N
const int N
Definition: speed_test.cpp:3
NonlinearElasticElement::CommonData::sTress
std::vector< MatrixDouble3by3 > sTress
Definition: NonLinearElasticElement.hpp:111
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
NonlinearElasticElement::k
FTensor::Index< 'k', 3 > k
Definition: NonLinearElasticElement.hpp:125
MoFEM::EntitiesFieldData::EntData::getFTensor0FieldData
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
Definition: EntitiesFieldData.cpp:506
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
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
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
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
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
NonlinearElasticElement::OpGetDataAtGaussPts
Definition: NonLinearElasticElement.hpp:342
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
NonlinearElasticElement::OpEnergy
Definition: NonLinearElasticElement.hpp:540
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
NonlinearElasticElement::MyVolumeFE::MyVolumeFE
MyVolumeFE(MoFEM::Interface &m_field)
Definition: NonLinearElasticElement.cpp:17
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
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::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
NonlinearElasticElement::feEnergy
MyVolumeFE feEnergy
calculate elastic energy
Definition: NonLinearElasticElement.hpp:70
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Calculate stress or jacobian at gauss points.
Definition: NonLinearElasticElement.cpp:341
NonlinearElasticElement::OpJacobianEnergy::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: NonLinearElasticElement.cpp:529
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
NonlinearElasticElement::OpRhsEshelbyStress::OpRhsEshelbyStress
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:1062
MoFEM::Types::MatrixAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132
NonlinearElasticElement::OpRhsPiolaKirchhoff::aSemble
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: NonLinearElasticElement.cpp:601
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getJac
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
Definition: VolumeElementForcesAndSourcesCore.hpp:170
NonlinearElasticElement::OpJacobianEnergy::calculateEnergy
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition: NonLinearElasticElement.cpp:423
NonlinearElasticElement::NonlinearElasticElement
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
Definition: NonLinearElasticElement.cpp:81
MoFEM::SmartPetscObj< Vec >
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
NonlinearElasticElement::OpLhsEshelby_dX::OpLhsEshelby_dX
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:1076
NonlinearElasticElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition: NonLinearElasticElement.hpp:108
H
double H
Hardening.
Definition: plastic.cpp:124
NonlinearElasticElement::setBlocks
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
Definition: NonLinearElasticElement.cpp:1086
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
F
@ F
Definition: free_surface.cpp:394
NonlinearElasticElement::CommonData
common data used by volume elements
Definition: NonLinearElasticElement.hpp:105
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::aSemble
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: NonLinearElasticElement.cpp:818