v0.10.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  * This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #include <MoFEM.hpp>
21 using namespace MoFEM;
23 
24 #include <boost/numeric/ublas/vector_proxy.hpp>
25 #include <boost/numeric/ublas/matrix.hpp>
26 #include <boost/numeric/ublas/matrix_proxy.hpp>
27 #include <boost/numeric/ublas/vector.hpp>
28 
29 #include <adolc/adolc.h>
31 
33  : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULL), F(PETSC_NULL),
34  addToRule(1) {
35  int ghosts[] = {0};
36  if (mField.get_comm_rank() == 0) {
37  ierr = VecCreateGhost(mField.get_comm(), 1, 1, 0, ghosts, &V);
38  } else {
39  ierr = VecCreateGhost(mField.get_comm(), 0, 1, 1, ghosts, &V);
40  }
41  CHKERRABORT(PETSC_COMM_SELF, ierr);
42 }
43 
45  ierr = VecDestroy(&V);
46  CHKERRABORT(PETSC_COMM_SELF, ierr);
47 }
48 
50  return 2 * (order - 1) + addToRule;
51 };
52 
55 
56  CHKERR VolumeElementForcesAndSourcesCore::preProcess();
57 
58  if (A != PETSC_NULL) {
59  snes_B = A;
60  }
61 
62  if (F != PETSC_NULL) {
63  snes_f = F;
64  }
65 
66  switch (snes_ctx) {
67  case CTX_SNESNONE:
68  CHKERR VecZeroEntries(V);
69  CHKERR VecGhostUpdateBegin(V, INSERT_VALUES, SCATTER_FORWARD);
70  CHKERR VecGhostUpdateEnd(V, INSERT_VALUES, SCATTER_FORWARD);
71  break;
72  default:
73  break;
74  }
75 
77 }
78 
81 
82  double *array;
83 
84  switch (snes_ctx) {
85  case CTX_SNESNONE:
86  CHKERR VecAssemblyBegin(V);
87  CHKERR VecAssemblyEnd(V);
88  CHKERR VecGhostUpdateBegin(V, ADD_VALUES, SCATTER_REVERSE);
89  CHKERR VecGhostUpdateEnd(V, ADD_VALUES, SCATTER_REVERSE);
90  CHKERR VecGhostUpdateBegin(V, INSERT_VALUES, SCATTER_FORWARD);
91  CHKERR VecGhostUpdateEnd(V, INSERT_VALUES, SCATTER_FORWARD);
92  CHKERR VecGetArray(V, &array);
93  eNergy = array[0];
94  CHKERR VecRestoreArray(V, &array);
95  break;
96  default:
97  break;
98  }
99 
100  CHKERR VolumeElementForcesAndSourcesCore::postProcess();
101 
103 }
104 
106  short int tag)
107  : feRhs(m_field), feLhs(m_field), feEnergy(m_field), mField(m_field),
108  tAg(tag) {}
109 
111  const std::string field_name,
112  std::vector<VectorDouble> &values_at_gauss_pts,
113  std::vector<MatrixDouble> &gardient_at_gauss_pts)
115  field_name, UserDataOperator::OPROW),
116  valuesAtGaussPts(values_at_gauss_pts),
117  gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
118 
120  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
122 
123  const int nb_dofs = data.getFieldData().size();
124  const int nb_base_functions = data.getN().size2();
125  if (nb_dofs == 0) {
127  }
128  const int nb_gauss_pts = data.getN().size1();
129  const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
130 
131  // initialize
132  if (type == zeroAtType) {
133  valuesAtGaussPts.resize(nb_gauss_pts);
134  gradientAtGaussPts.resize(nb_gauss_pts);
135  for (int gg = 0; gg != nb_gauss_pts; gg++) {
136  valuesAtGaussPts[gg].resize(rank, false);
137  gradientAtGaussPts[gg].resize(rank, 3, false);
138  }
139  for (int gg = 0; gg != nb_gauss_pts; gg++) {
140  valuesAtGaussPts[gg].clear();
141  gradientAtGaussPts[gg].clear();
142  }
143  }
144 
145  auto base_function = data.getFTensor0N();
146  auto diff_base_functions = data.getFTensor1DiffN<3>();
149 
150  if (rank == 1) {
151 
152  for (int gg = 0; gg != nb_gauss_pts; gg++) {
153  auto field_data = data.getFTensor0FieldData();
154  double &val = valuesAtGaussPts[gg][0];
155  FTensor::Tensor1<double *, 3> grad(&gradientAtGaussPts[gg](0, 0),
156  &gradientAtGaussPts[gg](0, 1),
157  &gradientAtGaussPts[gg](0, 2));
158  int bb = 0;
159  for (; bb != nb_dofs; bb++) {
160  val += base_function * field_data;
161  grad(i) += diff_base_functions(i) * field_data;
162  ++diff_base_functions;
163  ++base_function;
164  ++field_data;
165  }
166  for (; bb != nb_base_functions; bb++) {
167  ++diff_base_functions;
168  ++base_function;
169  }
170  }
171 
172  } else if (rank == 3) {
173 
174  for (int gg = 0; gg != nb_gauss_pts; gg++) {
175  auto field_data = data.getFTensor1FieldData<3>();
176  FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
177  &valuesAtGaussPts[gg][1],
178  &valuesAtGaussPts[gg][2]);
180  &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
181  &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
182  &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
183  &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
184  &gradientAtGaussPts[gg](2, 2));
185  int bb = 0;
186  for (; bb != nb_dofs / 3; bb++) {
187  values(i) += base_function * field_data(i);
188  gradient(i, j) += field_data(i) * diff_base_functions(j);
189  ++diff_base_functions;
190  ++base_function;
191  ++field_data;
192  }
193  for (; bb != nb_base_functions; bb++) {
194  ++diff_base_functions;
195  ++base_function;
196  }
197  }
198 
199  } else {
200  // FIXME: THat part is inefficient
201  VectorDouble &values = data.getFieldData();
202  for (int gg = 0; gg < nb_gauss_pts; gg++) {
203  VectorAdaptor N = data.getN(gg, nb_dofs / rank);
204  MatrixAdaptor diffN = data.getDiffN(gg, nb_dofs / rank);
205  for (int dd = 0; dd < nb_dofs / rank; dd++) {
206  for (int rr1 = 0; rr1 < rank; rr1++) {
207  valuesAtGaussPts[gg][rr1] += N[dd] * values[rank * dd + rr1];
208  for (int rr2 = 0; rr2 < 3; rr2++) {
209  gradientAtGaussPts[gg](rr1, rr2) +=
210  diffN(dd, rr2) * values[rank * dd + rr1];
211  }
212  }
213  }
214  }
215  }
216 
218 }
219 
221  const std::string field_name, CommonData &common_data)
222  : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
223  common_data.gradAtGaussPts[field_name]) {}
224 
226  OpJacobianPiolaKirchhoffStress(const std::string field_name,
227  BlockData &data, CommonData &common_data,
228  int tag, bool jacobian, bool ale,
229  bool field_disp)
231  field_name, UserDataOperator::OPROW),
232  dAta(data), commonData(common_data), tAg(tag), adlocReturnValue(0),
233  jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
234 
235 }
236 
239  const int gg) {
241 
242  ierr = dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
243  dAta, getNumeredEntFiniteElementPtr());
244  CHKERRG(ierr);
245  if (aLe) {
246  dAta.materialAdoublePtr->P =
247  dAta.materialAdoublePtr->detH *
248  prod(dAta.materialAdoublePtr->P, trans(dAta.materialAdoublePtr->invH));
249  }
250  commonData.sTress[gg].resize(3, 3, false);
251  for (int dd1 = 0; dd1 < 3; dd1++) {
252  for (int dd2 = 0; dd2 < 3; dd2++) {
253  dAta.materialAdoublePtr->P(dd1, dd2) >>=
254  (commonData.sTress[gg])(dd1, dd2);
255  }
256  }
257 
259 }
260 
263  const int gg) {
265 
266  trace_on(tAg, 0);
267 
268  dAta.materialAdoublePtr->F.resize(3, 3, false);
269 
270  if (!aLe) {
271 
272  nbActiveVariables = 0;
273  for (int dd1 = 0; dd1 < 3; dd1++) {
274  for (int dd2 = 0; dd2 < 3; dd2++) {
275  dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
276  if (fieldDisp) {
277  if (dd1 == dd2) {
278  dAta.materialAdoublePtr->F(dd1, dd2) += 1;
279  }
280  }
281  nbActiveVariables++;
282  }
283  }
284 
285  } else {
286 
287  nbActiveVariables = 0;
288 
289  dAta.materialAdoublePtr->h.resize(3, 3, false);
290  for (int dd1 = 0; dd1 < 3; dd1++) {
291  for (int dd2 = 0; dd2 < 3; dd2++) {
292  dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
293  nbActiveVariables++;
294  }
295  }
296 
297  dAta.materialAdoublePtr->H.resize(3, 3, false);
298  for (int dd1 = 0; dd1 < 3; dd1++) {
299  for (int dd2 = 0; dd2 < 3; dd2++) {
300  dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
301  nbActiveVariables++;
302  }
303  }
304 
305  CHKERR dAta.materialAdoublePtr->dEterminant(dAta.materialAdoublePtr->H,
306  dAta.materialAdoublePtr->detH);
307  dAta.materialAdoublePtr->invH.resize(3, 3, false);
308  CHKERR dAta.materialAdoublePtr->iNvert(dAta.materialAdoublePtr->detH,
309  dAta.materialAdoublePtr->H,
310  dAta.materialAdoublePtr->invH);
311  noalias(dAta.materialAdoublePtr->F) =
312  prod(dAta.materialAdoublePtr->h, dAta.materialAdoublePtr->invH);
313  }
314 
315  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
316  CHKERR calculateStress(gg);
317 
318  trace_off();
319 
321 }
322 
326 
327  int r;
328 
329  if (fUnction) {
330  commonData.sTress[gg].resize(3, 3, false);
331  // play recorder for values
332  r = ::function(tAg, 9, nbActiveVariables, &activeVariables[0],
333  &commonData.sTress[gg](0, 0));
334  if (r < adlocReturnValue) { // function is locally analytic
335  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
336  "ADOL-C function evaluation with error r = %d", r);
337  }
338  }
339 
340  if (jAcobian) {
341  commonData.jacStress[gg].resize(9, nbActiveVariables, false);
342  double *jac_ptr[] = {
343  &(commonData.jacStress[gg](0, 0)), &(commonData.jacStress[gg](1, 0)),
344  &(commonData.jacStress[gg](2, 0)), &(commonData.jacStress[gg](3, 0)),
345  &(commonData.jacStress[gg](4, 0)), &(commonData.jacStress[gg](5, 0)),
346  &(commonData.jacStress[gg](6, 0)), &(commonData.jacStress[gg](7, 0)),
347  &(commonData.jacStress[gg](8, 0))};
348  // play recorder for jacobians
349  r = jacobian(tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
350  if (r < adlocReturnValue) {
351  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
352  "ADOL-C function evaluation with error");
353  }
354  }
355 
357 }
358 
360  int row_side, EntityType row_type,
363 
364  // do it only once, no need to repeat this for edges,faces or tets
365  if (row_type != MBVERTEX)
367 
368  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
369  dAta.tEts.end()) {
371  }
372 
373  int nb_dofs = row_data.getFieldData().size();
374  if (nb_dofs == 0)
376  dAta.materialAdoublePtr->commonDataPtr = &commonData;
377  dAta.materialAdoublePtr->opPtr = this;
378 
379  int nb_gauss_pts = row_data.getN().size1();
380  commonData.sTress.resize(nb_gauss_pts);
381  commonData.jacStress.resize(nb_gauss_pts);
382 
384  if (aLe) {
386  }
387 
388  for (int gg = 0; gg != nb_gauss_pts; gg++) {
389 
390  dAta.materialAdoublePtr->gG = gg;
391 
392  // Record tag and calculate stress
393  if (recordTagForIntegrationPoint(gg)) {
394  CHKERR recordTag(gg);
395  }
396 
397  // Set active variables vector
398  if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
399  activeVariables.resize(nbActiveVariables, false);
400  if (!aLe) {
401  for (int dd1 = 0; dd1 < 3; dd1++) {
402  for (int dd2 = 0; dd2 < 3; dd2++) {
403  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
404  }
405  }
406  } else {
407  for (int dd1 = 0; dd1 < 3; dd1++) {
408  for (int dd2 = 0; dd2 < 3; dd2++) {
409  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
410  }
411  }
412  for (int dd1 = 0; dd1 < 3; dd1++) {
413  for (int dd2 = 0; dd2 < 3; dd2++) {
414  activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
415  }
416  }
417  }
418  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
419 
420  // Play tag and calculate stress or tangent
421  if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
422  CHKERR playTag(gg);
423  }
424  }
425  }
426 
428 }
429 
431  const std::string
432  field_name, ///< field name for spatial positions or displacements
433  BlockData &data, CommonData &common_data, int tag, bool gradient,
434  bool hessian, bool ale, bool field_disp)
436  field_name, UserDataOperator::OPROW),
437  dAta(data), commonData(common_data), tAg(tag), gRadient(gradient),
438  hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
439 
443  CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
444  dAta, getNumeredEntFiniteElementPtr());
445  dAta.materialAdoublePtr->eNergy >>= commonData.eNergy[gg];
447 }
448 
452 
453  trace_on(tAg, 0);
454 
455  dAta.materialAdoublePtr->F.resize(3, 3, false);
456 
457  if (!aLe) {
458 
459  nbActiveVariables = 0;
460  for (int dd1 = 0; dd1 < 3; dd1++) {
461  for (int dd2 = 0; dd2 < 3; dd2++) {
462  dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
463  if (fieldDisp) {
464  if (dd1 == dd2) {
465  dAta.materialAdoublePtr->F(dd1, dd2) += 1;
466  }
467  }
468  nbActiveVariables++;
469  }
470  }
471 
472  } else {
473 
474  nbActiveVariables = 0;
475 
476  dAta.materialAdoublePtr->h.resize(3, 3, false);
477  for (int dd1 = 0; dd1 < 3; dd1++) {
478  for (int dd2 = 0; dd2 < 3; dd2++) {
479  dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
480  nbActiveVariables++;
481  }
482  }
483 
484  dAta.materialAdoublePtr->H.resize(3, 3, false);
485  for (int dd1 = 0; dd1 < 3; dd1++) {
486  for (int dd2 = 0; dd2 < 3; dd2++) {
487  dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
488  nbActiveVariables++;
489  }
490  }
491 
492  CHKERR dAta.materialAdoublePtr->dEterminant(dAta.materialAdoublePtr->H,
493  dAta.materialAdoublePtr->detH);
494  dAta.materialAdoublePtr->invH.resize(3, 3, false);
495  CHKERR dAta.materialAdoublePtr->iNvert(dAta.materialAdoublePtr->detH,
496  dAta.materialAdoublePtr->H,
497  dAta.materialAdoublePtr->invH);
498  noalias(dAta.materialAdoublePtr->F) =
499  prod(dAta.materialAdoublePtr->h, dAta.materialAdoublePtr->invH);
500  }
501 
502  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
504 
505  trace_off();
506 
508 }
509 
513 
514  if (gRadient) {
515  commonData.jacEnergy[gg].resize(nbActiveVariables, false);
516  int r = ::gradient(tAg, nbActiveVariables, &activeVariables[0],
517  &commonData.jacEnergy[gg][0]);
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 
526  if (hEssian) {
527  commonData.hessianEnergy[gg].resize(nbActiveVariables * nbActiveVariables,
528  false);
529  double *H[nbActiveVariables];
530  for (int n = 0; n != nbActiveVariables; n++) {
531  H[n] = &(commonData.hessianEnergy[gg][n * nbActiveVariables]);
532  }
533  int r = ::hessian(tAg, nbActiveVariables, &*activeVariables.begin(), H);
534  if (r < 0) {
535  // That means that energy function is not smooth and derivative
536  // can not be calculated,
537  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
538  "ADOL-C function evaluation with error");
539  }
540  }
541 
543 }
544 
546  int row_side, EntityType row_type,
549 
550  // do it only once, no need to repeat this for edges,faces or tets
551  if (row_type != MBVERTEX)
553 
554  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
555  dAta.tEts.end()) {
557  }
558 
559  int nb_dofs = row_data.getFieldData().size();
560  if (nb_dofs == 0)
562  dAta.materialAdoublePtr->commonDataPtr = &commonData;
563  dAta.materialAdoublePtr->opPtr = this;
564 
565  int nb_gauss_pts = row_data.getN().size1();
566  commonData.eNergy.resize(nb_gauss_pts);
567  commonData.jacEnergy.resize(nb_gauss_pts);
568 
570  if (aLe) {
572  }
573 
574  for (int gg = 0; gg != nb_gauss_pts; gg++) {
575 
576  dAta.materialAdoublePtr->gG = gg;
577 
578  // Record tag and calualte stress
579  if (recordTagForIntegrationPoint(gg)) {
580  CHKERR recordTag(gg);
581  }
582 
583  activeVariables.resize(nbActiveVariables, false);
584  if (!aLe) {
585  for (int dd1 = 0; dd1 < 3; dd1++) {
586  for (int dd2 = 0; dd2 < 3; dd2++) {
587  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
588  }
589  }
590  } else {
591  for (int dd1 = 0; dd1 < 3; dd1++) {
592  for (int dd2 = 0; dd2 < 3; dd2++) {
593  activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
594  }
595  }
596  for (int dd1 = 0; dd1 < 3; dd1++) {
597  for (int dd2 = 0; dd2 < 3; dd2++) {
598  activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
599  }
600  }
601  }
602  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
603 
604  // Play tag and calculate stress or tangent
605  CHKERR playTag(gg);
606  }
607 
609 }
610 
612  const std::string field_name, BlockData &data, CommonData &common_data)
614  field_name, UserDataOperator::OPROW),
615  dAta(data), commonData(common_data), aLe(false) {}
616 
618  int row_side, EntityType row_type,
621 
622  int nb_dofs = row_data.getIndices().size();
623  int *indices_ptr = &row_data.getIndices()[0];
624  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
625  iNdices.resize(nb_dofs, false);
626  noalias(iNdices) = row_data.getIndices();
627  indices_ptr = &iNdices[0];
628  VectorDofs &dofs = row_data.getFieldDofs();
629  VectorDofs::iterator dit = dofs.begin();
630  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
631  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
632  dAta.forcesOnlyOnEntitiesRow.end()) {
633  iNdices[ii] = -1;
634  }
635  }
636  }
637  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
638  ADD_VALUES);
640 }
641 
643  int row_side, EntityType row_type,
646 
647  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
648  dAta.tEts.end()) {
650  }
651 
652  const int nb_dofs = row_data.getIndices().size();
653  if (nb_dofs == 0)
655  if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
656  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
657  }
658  const int nb_base_functions = row_data.getN().size2();
659  const int nb_gauss_pts = row_data.getN().size1();
660 
661  nf.resize(nb_dofs, false);
662  nf.clear();
663 
664  FTensor::Tensor1<double *, 3> diff_base_functions =
665  row_data.getFTensor1DiffN<3>();
668 
669  for (int gg = 0; gg != nb_gauss_pts; gg++) {
670  double val = getVolume() * getGaussPts()(3, gg);
671  if ((!aLe) && getHoGaussPtsDetJac().size() > 0) {
672  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
673  }
674  MatrixDouble3by3 &stress = commonData.sTress[gg];
676  &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
677  &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
678  &stress(2, 2));
679  FTensor::Tensor1<double *, 3> rhs(&nf[0], &nf[1], &nf[2], 3);
680  int bb = 0;
681  for (; bb != nb_dofs / 3; bb++) {
682  rhs(i) += val * t3(i, j) * diff_base_functions(j);
683  ++rhs;
684  ++diff_base_functions;
685  }
686  for (; bb != nb_base_functions; bb++) {
687  ++diff_base_functions;
688  }
689  }
690 
691  CHKERR aSemble(row_side, row_type, row_data);
692 
694 }
695 
696 NonlinearElasticElement::OpEnergy::OpEnergy(const std::string field_name,
697  BlockData &data,
698  CommonData &common_data, Vec ghost_vec,
699  bool field_disp)
701  field_name, UserDataOperator::OPROW),
702  dAta(data), commonData(common_data), ghostVec(ghost_vec), fieldDisp(field_disp) {
703  ierr = PetscObjectReference((PetscObject)ghostVec);
704  CHKERRABORT(PETSC_COMM_SELF, ierr);
705 }
706 
708  ierr = VecDestroy(&ghostVec);
709  CHKERRABORT(PETSC_COMM_SELF, ierr);
710 }
711 
713  int row_side, EntityType row_type,
716 
717  if (row_type != MBVERTEX)
719  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
720  dAta.tEts.end()) {
722  }
723 
724  std::vector<MatrixDouble> &F =
726  dAta.materialDoublePtr->F.resize(3, 3, false);
727 
728  double *energy_ptr;
729  CHKERR VecGetArray(ghostVec, &energy_ptr);
730 
731  for (unsigned int gg = 0; gg != row_data.getN().size1(); ++gg) {
732  double val = getVolume() * getGaussPts()(3, gg);
733  if (getHoGaussPtsDetJac().size() > 0) {
734  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
735  }
736 
737  noalias(dAta.materialDoublePtr->F) = F[gg];
738  if (fieldDisp) {
739  for (int dd = 0; dd < 3; dd++) {
740  dAta.materialDoublePtr->F(dd, dd) += 1;
741  }
742  }
743 
744  int nb_active_variables = 0;
745  CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
746  CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
747  dAta, getNumeredEntFiniteElementPtr());
748  energy_ptr[0] += val * dAta.materialDoublePtr->eNergy;
749  }
750  CHKERR VecRestoreArray(ghostVec, &energy_ptr);
751 
753 }
754 
756  const std::string vel_field, const std::string field_name, BlockData &data,
757  CommonData &common_data)
759  vel_field, field_name, UserDataOperator::OPROWCOL),
760  dAta(data), commonData(common_data), aLe(false) {}
761 
762 template <int S>
764  int gg, MatrixDouble &jac_stress,
765  MatrixDouble &jac) {
766  jac.clear();
770  int nb_col = col_data.getFieldData().size();
771  double *diff_ptr =
772  const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
773  // First two indices 'i','j' derivatives of 1st Piola-stress, third index 'k'
774  // is displacement component
776  &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
777  &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
778  &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
779  &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
780  &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
781  &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
782  &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
783  &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
784  &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
785  &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
786  &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
787  &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
788  &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
789  &jac_stress(3 * 2 + 2, S + 2));
791  &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
792  &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
793  &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
794  &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
795  &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
796  &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
797  &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
798  &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
799  &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
800  &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
801  &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
802  &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
803  &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
804  &jac_stress(3 * 2 + 2, S + 5));
806  &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
807  &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
808  &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
809  &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
810  &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
811  &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
812  &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
813  &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
814  &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
815  &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
816  &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
817  &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
818  &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
819  &jac_stress(3 * 2 + 2, S + 8));
820  // Derivate of 1st Piola-stress multiplied by gradient of defamation for
821  // base function (dd) and displacement component (rr)
823  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
824  &jac(6, 0), &jac(7, 0), &jac(8, 0));
826  &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
827  &jac(6, 1), &jac(7, 1), &jac(8, 1));
829  &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
830  &jac(6, 2), &jac(7, 2), &jac(8, 2));
832  diff_ptr, &diff_ptr[1], &diff_ptr[2]);
833  for (int dd = 0; dd != nb_col / 3; ++dd) {
834  t2_1_0(i, j) += t3_1_0(i, j, k) * diff(k);
835  t2_1_1(i, j) += t3_1_1(i, j, k) * diff(k);
836  t2_1_2(i, j) += t3_1_2(i, j, k) * diff(k);
837  ++t2_1_0;
838  ++t2_1_1;
839  ++t2_1_2;
840  ++diff;
841  }
843 }
844 
846  DataForcesAndSourcesCore::EntData &col_data, int gg) {
847  return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
848 }
849 
851  int row_side, int col_side, EntityType row_type, EntityType col_type,
855 
856  int nb_row = row_data.getIndices().size();
857  int nb_col = col_data.getIndices().size();
858 
859  int *row_indices_ptr = &row_data.getIndices()[0];
860  int *col_indices_ptr = &col_data.getIndices()[0];
861 
862  /*for(int dd1 = 0;dd1<k.size1();dd1++) {
863  for(int dd2 = 0;dd2<k.size2();dd2++) {
864  if(k(dd1,dd2)!=k(dd1,dd2)) {
865  SETERRQ(PETSC_COMM_SELF,1,"Wrong result");
866  }
867  }
868  }*/
869 
870  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
871  rowIndices.resize(nb_row, false);
872  noalias(rowIndices) = row_data.getIndices();
873  row_indices_ptr = &rowIndices[0];
874  VectorDofs &dofs = row_data.getFieldDofs();
875  VectorDofs::iterator dit = dofs.begin();
876  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
877  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
878  dAta.forcesOnlyOnEntitiesRow.end()) {
879  rowIndices[ii] = -1;
880  }
881  }
882  }
883 
884  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
885  colIndices.resize(nb_col, false);
886  noalias(colIndices) = col_data.getIndices();
887  col_indices_ptr = &colIndices[0];
888  VectorDofs &dofs = col_data.getFieldDofs();
889  VectorDofs::iterator dit = dofs.begin();
890  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
891  if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
892  dAta.forcesOnlyOnEntitiesCol.end()) {
893  colIndices[ii] = -1;
894  }
895  }
896  }
897 
898  CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
899  col_indices_ptr, &k(0, 0), ADD_VALUES);
900 
901  // is symmetric
902  if (row_side != col_side || row_type != col_type) {
903 
904  row_indices_ptr = &row_data.getIndices()[0];
905  col_indices_ptr = &col_data.getIndices()[0];
906 
907  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
908  rowIndices.resize(nb_row, false);
909  noalias(rowIndices) = row_data.getIndices();
910  row_indices_ptr = &rowIndices[0];
911  VectorDofs &dofs = row_data.getFieldDofs();
912  VectorDofs::iterator dit = dofs.begin();
913  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
914  if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
915  dAta.forcesOnlyOnEntitiesCol.end()) {
916  rowIndices[ii] = -1;
917  }
918  }
919  }
920 
921  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
922  colIndices.resize(nb_col, false);
923  noalias(colIndices) = col_data.getIndices();
924  col_indices_ptr = &colIndices[0];
925  VectorDofs &dofs = col_data.getFieldDofs();
926  VectorDofs::iterator dit = dofs.begin();
927  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
928  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
929  dAta.forcesOnlyOnEntitiesRow.end()) {
930  colIndices[ii] = -1;
931  }
932  }
933  }
934 
935  trans_k.resize(nb_col, nb_row, false);
936  noalias(trans_k) = trans(k);
937  CHKERR MatSetValues(getFEMethod()->snes_B, nb_col, col_indices_ptr, nb_row,
938  row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
939  }
940 
942 }
943 
945  int row_side, int col_side, EntityType row_type, EntityType col_type,
949 
950  int nb_row = row_data.getIndices().size();
951  int nb_col = col_data.getIndices().size();
952  if (nb_row == 0)
954  if (nb_col == 0)
956 
957  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
958  dAta.tEts.end()) {
960  }
961 
962  // const int nb_base_functions = row_data.getN().size2();
963  const int nb_gauss_pts = row_data.getN().size1();
964 
968 
969  k.resize(nb_row, nb_col, false);
970  k.clear();
971  jac.resize(9, nb_col, false);
972 
973  for (int gg = 0; gg != nb_gauss_pts; gg++) {
974  CHKERR getJac(col_data, gg);
975  double val = getVolume() * getGaussPts()(3, gg);
976  if ((!aLe) && (getHoGaussPtsDetJac().size() > 0)) {
977  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
978  }
980  &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
981  &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
982  &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
983  &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
984  &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
985  &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
986  &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
987  &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
988  &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
989  for (int cc = 0; cc != nb_col / 3; cc++) {
990  FTensor::Tensor1<double *, 3> diff_base_functions =
991  row_data.getFTensor1DiffN<3>(gg, 0);
993  &k(0, 3 * cc + 0), &k(0, 3 * cc + 1), &k(0, 3 * cc + 2),
994  &k(1, 3 * cc + 0), &k(1, 3 * cc + 1), &k(1, 3 * cc + 2),
995  &k(2, 3 * cc + 0), &k(2, 3 * cc + 1), &k(2, 3 * cc + 2), 3 * nb_col);
996  for (int rr = 0; rr != nb_row / 3; rr++) {
997  lhs(i, j) += val * t3_1(i, m, j) * diff_base_functions(m);
998  ++diff_base_functions;
999  ++lhs;
1000  }
1001  ++t3_1;
1002  }
1003  }
1004 
1005  CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
1006 
1008 }
1009 
1011  const std::string vel_field, const std::string field_name, BlockData &data,
1012  CommonData &common_data)
1013  : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {
1014  sYmm = false;
1015 }
1016 
1018  DataForcesAndSourcesCore::EntData &col_data, int gg) {
1019  return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
1020 }
1021 
1023  int row_side, int col_side, EntityType row_type, EntityType col_type,
1027 
1028  int nb_row = row_data.getIndices().size();
1029  int nb_col = col_data.getIndices().size();
1030 
1031  int *row_indices_ptr = &row_data.getIndices()[0];
1032  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
1033  rowIndices.resize(nb_row, false);
1034  noalias(rowIndices) = row_data.getIndices();
1035  row_indices_ptr = &rowIndices[0];
1036  VectorDofs &dofs = row_data.getFieldDofs();
1037  VectorDofs::iterator dit = dofs.begin();
1038  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1039  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1040  dAta.forcesOnlyOnEntitiesRow.end()) {
1041  rowIndices[ii] = -1;
1042  }
1043  }
1044  }
1045 
1046  int *col_indices_ptr = &col_data.getIndices()[0];
1047  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
1048  colIndices.resize(nb_col, false);
1049  noalias(colIndices) = col_data.getIndices();
1050  col_indices_ptr = &colIndices[0];
1051  VectorDofs &dofs = col_data.getFieldDofs();
1052  VectorDofs::iterator dit = dofs.begin();
1053  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1054  if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
1055  dAta.forcesOnlyOnEntitiesCol.end()) {
1056  colIndices[ii] = -1;
1057  }
1058  }
1059  }
1060 
1061  /*for(int dd1 = 0;dd1<k.size1();dd1++) {
1062  for(int dd2 = 0;dd2<k.size2();dd2++) {
1063  if(k(dd1,dd2)!=k(dd1,dd2)) {
1064  SETERRQ(PETSC_COMM_SELF,1,"Wrong result");
1065  }
1066  }
1067  }*/
1068 
1069  CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr, nb_col,
1070  col_indices_ptr, &k(0, 0), ADD_VALUES);
1071 
1073 }
1074 
1076  const std::string field_name, BlockData &data, CommonData &common_data,
1077  int tag, bool jacobian, bool ale)
1078  : OpJacobianPiolaKirchhoffStress(field_name, data, common_data, tag,
1079  jacobian, ale, false) {}
1080 
1083  const int gg) {
1085 
1086  CHKERR dAta.materialAdoublePtr->calculateSiGma_EshelbyStress(
1087  dAta, getNumeredEntFiniteElementPtr());
1088  if (aLe) {
1089  dAta.materialAdoublePtr->SiGma = dAta.materialAdoublePtr->detH *
1090  prod(dAta.materialAdoublePtr->SiGma,
1091  trans(dAta.materialAdoublePtr->invH));
1092  }
1093  commonData.sTress[gg].resize(3, 3, false);
1094  for (int dd1 = 0; dd1 < 3; dd1++) {
1095  for (int dd2 = 0; dd2 < 3; dd2++) {
1096  dAta.materialAdoublePtr->SiGma(dd1, dd2) >>=
1097  (commonData.sTress[gg])(dd1, dd2);
1098  }
1099  }
1100 
1102 }
1103 
1105  const std::string field_name, BlockData &data, CommonData &common_data)
1106  : OpRhsPiolaKirchhoff(field_name, data, common_data) {}
1107 
1109  const std::string vel_field, const std::string field_name, BlockData &data,
1110  CommonData &common_data)
1111  : OpLhsPiolaKirchhoff_dX(vel_field, field_name, data, common_data) {}
1112 
1114  DataForcesAndSourcesCore::EntData &col_data, int gg) {
1115  return get_jac<0>(col_data, gg, commonData.jacStress[gg], jac);
1116 }
1117 
1119  const std::string vel_field, const std::string field_name, BlockData &data,
1120  CommonData &common_data)
1121  : OpLhsPiolaKirchhoff_dx(vel_field, field_name, data, common_data) {}
1122 
1124  DataForcesAndSourcesCore::EntData &col_data, int gg) {
1125  return get_jac<9>(col_data, gg, commonData.jacStress[gg], jac);
1126 }
1127 
1130  materialDoublePtr,
1132  materialAdoublePtr) {
1134 
1135  if (!materialDoublePtr) {
1137  "Pointer for materialDoublePtr not allocated");
1138  }
1139  if (!materialAdoublePtr) {
1141  "Pointer for materialAdoublePtr not allocated");
1142  }
1143 
1145  mField, BLOCKSET | MAT_ELASTICSET, it)) {
1146  Mat_Elastic mydata;
1147  CHKERR it->getAttributeDataStructure(mydata);
1148  int id = it->getMeshsetId();
1149  EntityHandle meshset = it->getMeshset();
1150  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1151  setOfBlocks[id].tEts, true);
1152  setOfBlocks[id].iD = id;
1153  setOfBlocks[id].E = mydata.data.Young;
1154  setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
1155  setOfBlocks[id].materialDoublePtr = materialDoublePtr;
1156  setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
1157  }
1158 
1160 }
1161 
1163  const std::string element_name,
1164  const std::string spatial_position_field_name,
1165  const std::string material_position_field_name, const bool ale) {
1167 
1168  CHKERR mField.add_finite_element(element_name, MF_ZERO);
1170  element_name, spatial_position_field_name);
1172  element_name, spatial_position_field_name);
1174  element_name, spatial_position_field_name);
1175  if (mField.check_field(material_position_field_name)) {
1176  if (ale) {
1178  element_name, material_position_field_name);
1180  element_name, material_position_field_name);
1181  }
1183  element_name, material_position_field_name);
1184  }
1185 
1186  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1187  for (; sit != setOfBlocks.end(); sit++) {
1188  CHKERR mField.add_ents_to_finite_element_by_type(sit->second.tEts, MBTET,
1189  element_name);
1190  }
1191 
1193 }
1194 
1196  const std::string spatial_position_field_name,
1197  const std::string material_position_field_name, const bool ale,
1198  const bool field_disp) {
1200 
1201  commonData.spatialPositions = spatial_position_field_name;
1202  commonData.meshPositions = material_position_field_name;
1203 
1204  // Rhs
1205  feRhs.getOpPtrVector().push_back(
1206  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1207  if (mField.check_field(material_position_field_name)) {
1209  material_position_field_name, commonData));
1210  }
1211  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1212  for (; sit != setOfBlocks.end(); sit++) {
1214  spatial_position_field_name, sit->second, commonData, tAg, false, ale,
1215  field_disp));
1216  feRhs.getOpPtrVector().push_back(new OpRhsPiolaKirchhoff(
1217  spatial_position_field_name, sit->second, commonData));
1218  }
1219 
1220  // Energy
1221  feEnergy.getOpPtrVector().push_back(
1222  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1223  if (mField.check_field(material_position_field_name)) {
1225  material_position_field_name, commonData));
1226  }
1227  sit = setOfBlocks.begin();
1228  for (; sit != setOfBlocks.end(); sit++) {
1229  feEnergy.getOpPtrVector().push_back(
1230  new OpEnergy(spatial_position_field_name, sit->second, commonData,
1231  feEnergy.V, field_disp));
1232  }
1233 
1234  // Lhs
1235  feLhs.getOpPtrVector().push_back(
1236  new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
1237  if (mField.check_field(material_position_field_name)) {
1239  material_position_field_name, commonData));
1240  }
1241  sit = setOfBlocks.begin();
1242  for (; sit != setOfBlocks.end(); sit++) {
1244  spatial_position_field_name, sit->second, commonData, tAg, true, ale,
1245  field_disp));
1247  spatial_position_field_name, spatial_position_field_name, sit->second,
1248  commonData));
1249  }
1250 
1252 }
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.
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
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
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
Deprecated interface functions.
virtual moab::Interface & get_moab()=0
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.
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
std::vector< MatrixDouble3by3 > sTress
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
Elastic material data structure.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
MyVolumeFE feEnergy
calculate elastic energy
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
constexpr double H
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
MoFEMErrorCode postProcess()
function is run at the end of loop
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
static Index< 'n', 3 > n
static MoFEMErrorCode get_jac(DataForcesAndSourcesCore::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
static Index< 'm', 3 > m
bool sYmm
If true assume that matrix is symmetric structure.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
virtual int get_comm_rank() const =0
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
constexpr int order
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
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
static Index< 'i', 3 > i
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
common data used by volume elements
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Calculate stress or jacobian at gauss points.
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:126
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, Vec ghost_vec, bool field_disp)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
MoFEMErrorCode preProcess()
function is run at the beginning of loop
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
static Index< 'j', 3 > j
Operators and data structures for non-linear elastic analysis.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:109
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
int getRule(int order)
it is used to calculate nb. of Gauss integration points
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
const double r
rate factor
block name is "MAT_ELASTIC"
Definition: definitions.h:228
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)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
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, Vec *v_energy_ptr)
static Index< 'k', 3 > k
virtual bool check_field(const std::string &name) const =0
check if field is in database
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating deformation gradient
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
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
const VectorDouble & getFieldData() const
get dofs values
data for calculation heat conductivity and heat capacity elements
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
const int N
Definition: speed_test.cpp:3
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)