v0.14.0
GriffithForceElement.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #ifndef __GRIFFITH_FORCE_ELEMENT_HPP__
16 #define __GRIFFITH_FORCE_ELEMENT_HPP__
17 
18 #ifndef WITH_ADOL_C
19 #error "MoFEM need to be compiled with ADOL-C"
20 #endif
21 
22 namespace FractureMechanics {
23 
24 static double calMax(double a, double b, double r) {
25  return (a + b + (1 / r) * pow(fabs(a - b), r)) / 2;
26 }
27 
28 static double diffCalMax_a(double a, double b, double r) {
29  double sgn = ((a - b) == 0) ? 0 : (((a - b) < 0) ? -1 : 1);
30  return (1 + pow(fabs(a - b), r - 1) * sgn * (+1)) / 2;
31 }
32 
33 /** \brief Implementation of Griffith element
34  */
35 struct GriffithForceElement {
36 
37  MoFEM::Interface &mField;
38 
39  struct MyTriangleFE : public FaceElementForcesAndSourcesCore {
40 
41  Mat B;
42  Vec F;
43 
44  MyTriangleFE(MoFEM::Interface &m_field)
45  : FaceElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
46  F(PETSC_NULL) {}
47 
48  int getRule(int order) { return 1; };
49  };
50 
51  boost::shared_ptr<MyTriangleFE> feRhsPtr;
52  boost::shared_ptr<MyTriangleFE> feLhsPtr;
53 
54  MyTriangleFE &feRhs;
55  MyTriangleFE &getLoopFeRhs() { return feRhs; }
56  MyTriangleFE &feLhs;
57  MyTriangleFE &getLoopFeLhs() { return feLhs; }
58 
59  GriffithForceElement(MoFEM::Interface &m_field)
60  : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
61  feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr),
62  feLhs(*feLhsPtr) {}
63 
64  struct CommonData {
65  VectorDouble griffithForce; ///< Griffith force at nodes
66  MatrixDouble tangentGriffithForce; ///< Tangent matrix
67  vector<double *>
68  tangentGriffithForceRowPtr; ///< Pointers to rows of tangent matrix
69  VectorDouble densityRho; ///< gc * rho^beta
70  MatrixDouble diffRho; ///< for lhs with heterogeneous gc
71  };
72  CommonData commonData; ///< Common data sgared between operators
73 
74  struct BlockData {
75  double gc; ///< Griffith energy
76  double E; ///< Young modulus
77  double r; ///< regularity parameter
78  Range frontNodes; ///< Nodes on crack front
79  BlockData() : r(1) {}
80  };
81  map<int, BlockData> blockData; ///< Shared data between finite elements
82 
83  static MoFEMErrorCode getElementOptions(BlockData &block_data) {
84  MoFEMFunctionBegin;
85 
86  MOFEM_LOG_CHANNEL("WORLD");
87  MOFEM_LOG_FUNCTION();
88  MOFEM_LOG_TAG("WORLD", "Griffith");
89 
90  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
91  "Get Griffith element options", "none");
92  CHKERRQ(ierr);
93  CHKERR PetscOptionsScalar("-griffith_gc", "release energy", "",
94  block_data.gc, &block_data.gc, PETSC_NULL);
95  CHKERR PetscOptionsScalar("-griffith_r", "release energy", "", block_data.r,
96  &block_data.r, PETSC_NULL);
97  CHKERR PetscOptionsScalar("-griffith_E", "stiffness", "", block_data.E,
98  &block_data.E, PETSC_NULL);
99  ierr = PetscOptionsEnd();
100  CHKERRQ(ierr);
101 
102  MOFEM_LOG_C("WORLD", Sev::inform, "### Input parameter: -griffith_E %6.4e",
103  block_data.E);
104  MOFEM_LOG_C("WORLD", Sev::inform, "### Input parameter: -griffith_r %6.4e",
105  block_data.r);
106 
107  MoFEMFunctionReturn(0);
108  }
109 
110  /**
111 
112  \f[
113  A = \| \mathbf{N} \| =
114  \left\|
115  \textrm{Spin}\left[
116  \frac{\partial\mathbf{N}\mathbf{X}}{\partial\xi}
117  \right]
118  \frac{\partial\mathbf{N}\mathbf{X}}{\partial\eta}
119  \right\|
120  \f]
121 
122  */
123  template <typename TYPE> struct AuxFunctions {
124 
125  ublas::vector<TYPE> referenceCoords;
126  ublas::vector<TYPE> currentCoords;
127  ublas::vector<TYPE> griffithForce;
128 
129  FTensor::Tensor1<TYPE, 3> t_Normal;
130  FTensor::Tensor1<TYPE, 3> t_Tangent1_ksi;
131  FTensor::Tensor1<TYPE, 3> t_Tangent2_eta;
132  TYPE currentArea;
133 
134  FTensor::Index<'i', 3> i;
135  FTensor::Index<'j', 3> j;
136  FTensor::Index<'k', 3> k;
137  FTensor::Number<0> N0;
138  FTensor::Number<1> N1;
139 
140  boost::function<FTensor::Tensor1<TYPE *, 3>(ublas::vector<TYPE> &v)>
141  get_ftensor_from_vec;
142 
143  boost::function<FTensor::Tensor1<FTensor::PackPtr<const double *, 2>, 2>(
144  const MatrixDouble &m)>
145  get_ftensor_from_mat_2d;
146 
147  AuxFunctions() {
148 
149  get_ftensor_from_vec = [](ublas::vector<TYPE> &v) {
150  return FTensor::Tensor1<TYPE *, 3>(&v(0), &v(1), &v(2), 3);
151  };
152 
153  get_ftensor_from_mat_2d = [](const MatrixDouble &m) {
154  return FTensor::Tensor1<FTensor::PackPtr<const double *, 2>, 2>(
155  &m(0, 0), &m(0, 1));
156  };
157  }
158 
159  MoFEMErrorCode calculateAreaAndNormal(const MatrixDouble &diff_n) {
160  MoFEMFunctionBeginHot;
161 
162  auto t_coords = get_ftensor_from_vec(currentCoords);
163  auto t_diff = get_ftensor_from_mat_2d(diff_n);
164  t_Tangent1_ksi(i) = (TYPE)0;
165  t_Tangent2_eta(i) = (TYPE)0;
166  for (int nn = 0; nn != 3; ++nn) {
167  t_Tangent1_ksi(i) += t_coords(i) * t_diff(N0);
168  t_Tangent2_eta(i) += t_coords(i) * t_diff(N1);
169  ++t_diff;
170  ++t_coords;
171  }
172 
173  t_Normal(j) =
174  FTensor::levi_civita(i, j, k) * t_Tangent1_ksi(k) * t_Tangent2_eta(i);
175  currentArea = 0.5 * sqrt(t_Normal(i) * t_Normal(i));
176 
177  MoFEMFunctionReturnHot(0);
178  }
179 
180  MoFEMErrorCode calculateMatA(const MatrixDouble &diff_n) {
181  MoFEMFunctionBeginHot;
182 
183  griffithForce.resize(9, false);
184  std::fill(griffithForce.begin(), griffithForce.end(), 0);
185  auto t_griffith = get_ftensor_from_vec(griffithForce);
186  auto t_diff = get_ftensor_from_mat_2d(diff_n);
187 
188  for (int dd = 0; dd != 3; dd++) {
189 
190  t_griffith(i) += FTensor::levi_civita(i, j, k) * t_Tangent1_ksi(k) *
191  t_diff(N1) * t_Normal(j) * 0.25 / currentArea;
192  t_griffith(i) -= FTensor::levi_civita(i, j, k) * t_Tangent2_eta(k) *
193  t_diff(N0) * t_Normal(j) * 0.25 / currentArea;
194  ++t_diff;
195  ++t_griffith;
196  }
197 
198  MoFEMFunctionReturnHot(0);
199  }
200 
201  MoFEMErrorCode calculateGriffithForce(const double gc, const double beta,
202  const MatrixDouble &diff_n) {
203  MoFEMFunctionBeginHot;
204 
205  CHKERR calculateAreaAndNormal(diff_n);
206  CHKERR calculateMatA(diff_n);
207 
208  for (int dd = 0; dd != 9; dd++)
209  griffithForce[dd] *= gc * beta;
210 
211  MoFEMFunctionReturnHot(0);
212  }
213  };
214 
215  struct AuxOp {
216 
217  int tAg;
218  BlockData &blockData;
219  CommonData &commonData;
220 
221  AuxOp(int tag, BlockData &block_data, CommonData &common_data)
222  : tAg(tag), blockData(block_data), commonData(common_data){};
223 
224  ublas::vector<int> rowIndices;
225  ublas::vector<int> rowIndicesLocal;
226 
227  VectorInt3 rowLambdaIndices;
228  VectorInt3 rowLambdaIndicesLocal;
229  VectorDouble3 lambdaAtNodes;
230  VectorDouble activeVariables;
231 
232  MoFEMErrorCode setIndices(DataForcesAndSourcesCore::EntData &data) {
233  MoFEMFunctionBegin;
234  const auto nb_dofs = data.getIndices().size();
235  rowIndices.resize(nb_dofs, false);
236  noalias(rowIndices) = data.getIndices();
237  rowIndicesLocal.resize(nb_dofs, false);
238  noalias(rowIndicesLocal) = data.getLocalIndices();
239  auto dit = data.getFieldDofs().begin();
240  auto hi_dit = data.getFieldDofs().end();
241  for (int ii = 0; dit != hi_dit; ++dit, ++ii) {
242  if (auto dof = (*dit)) {
243  if (blockData.frontNodes.find(dof->getEnt()) ==
244  blockData.frontNodes.end()) {
245 
246  const auto side = dof->getSideNumber();
247  const auto idx = dof->getEntDofIdx();
248  const auto rank = dof->getNbDofsOnEnt();
249  if (ii != side * rank + idx)
250  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong index");
251 
252  rowIndices[ii] = -1;
253  rowIndicesLocal[ii] = -1;
254  }
255  }
256  }
257  MoFEMFunctionReturn(0);
258  }
259 
260  MoFEMErrorCode
261  setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
262  DataForcesAndSourcesCore::EntData &data) {
263  MoFEMFunctionBeginHot;
264  int nb_dofs = data.getIndices().size();
265  activeVariables.resize(18, false);
266  for (int dd = 0; dd != nb_dofs; dd++) {
267  activeVariables[dd] = fe_ptr->getCoords()[dd];
268  }
269  for (int dd = 0; dd != nb_dofs; dd++) {
270  activeVariables[9 + dd] = data.getFieldData()[dd];
271  }
272  MoFEMFunctionReturnHot(0);
273  }
274 
275  MoFEMErrorCode
276  setLambdaNodes(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
277  const std::string &lambda_field_name) {
278  MoFEMFunctionBeginHot;
279  lambdaAtNodes.resize(3, false);
280  lambdaAtNodes.clear();
281  const auto bit_field_number =
282  fe_ptr->getFEMethod()->getFieldBitNumber(lambda_field_name);
283 
284  auto data_dofs = fe_ptr->getFEMethod()->getDataDofsPtr();
285  for (auto dit = data_dofs->get<Unique_mi_tag>().lower_bound(
286  FieldEntity::getLoBitNumberUId(bit_field_number));
287  dit != data_dofs->get<Unique_mi_tag>().upper_bound(
288  FieldEntity::getHiBitNumberUId(bit_field_number));
289  ++dit) {
290  if (dit->get()->getEntType() != MBVERTEX) {
291  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
292  "DOFs only on vertices");
293  }
294  int side = dit->get()->getSideNumber();
295  lambdaAtNodes[side] = dit->get()->getFieldData();
296  }
297  MoFEMFunctionReturnHot(0);
298  }
299 
300  MoFEMErrorCode
301  setLambdaIndices(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
302  const std::string &lambda_field_name) {
303  MoFEMFunctionBeginHot;
304  rowLambdaIndices.resize(3, false);
305  rowLambdaIndicesLocal.resize(3, false);
306  std::fill(rowLambdaIndices.begin(), rowLambdaIndices.end(), -1);
307  std::fill(rowLambdaIndicesLocal.begin(), rowLambdaIndicesLocal.end(), -1);
308  const auto bit_field_number =
309  fe_ptr->getFEMethod()->getFieldBitNumber(lambda_field_name);
310 
311  auto data_dofs = fe_ptr->getFEMethod()->getRowDofsPtr();
312  for (auto dit = data_dofs->get<Unique_mi_tag>().lower_bound(
313  FieldEntity::getLoBitNumberUId(bit_field_number));
314  dit != data_dofs->get<Unique_mi_tag>().upper_bound(
315  FieldEntity::getHiBitNumberUId(bit_field_number));
316  ++dit) {
317 
318  if (dit->get()->getEntType() != MBVERTEX)
319  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
320  "DOFs only on vertices");
321 
322  int side = dit->get()->getSideNumber();
323  rowLambdaIndices[side] = dit->get()->getPetscGlobalDofIdx();
324  rowLambdaIndicesLocal[side] = dit->get()->getPetscLocalDofIdx();
325  }
326  MoFEMFunctionReturnHot(0);
327  }
328  };
329 
330  struct OpRhs : public FaceElementForcesAndSourcesCore::UserDataOperator,
331  AuxOp {
332  OpRhs(int tag, BlockData &block_data, CommonData &common_data)
333  : FaceElementForcesAndSourcesCore::UserDataOperator(
334  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
335  AuxOp(tag, block_data, common_data) {}
336  MoFEMErrorCode doWork(int side, EntityType type,
337  DataForcesAndSourcesCore::EntData &data) {
338  MoFEMFunctionBegin;
339 
340  if (type != MBVERTEX)
341  MoFEMFunctionReturnHot(0);
342  CHKERR setIndices(data);
343  CHKERR setVariables(this, data);
344  int nb_dofs = data.getIndices().size();
345 
346  CHKERR VecSetValues(
347  getFEMethod()->snes_f, nb_dofs, &*rowIndices.data().begin(),
348  &*commonData.griffithForce.data().begin(), ADD_VALUES);
349  MoFEMFunctionReturn(0);
350  }
351  };
352 
353  /**
354  * @brief calculate griffith force vector
355  *
356  * \f[
357  \mathbf{f}_{\mathrm{m}}^{\mathrm{h}}=\frac{1}{2}\left(\tilde{\mathbf{A}}_{\Gamma}^{\mathrm{h}}\right)^{\mathrm{T}}
358  \mathbf{g}_{c}(\rho(\mathbf{X}))
359  * ]\f
360  *
361  * The matrix $\mathbf{A}_{\Gamma}^h$ comprising of direction vectors along
362  the crack front that are normal to the crack front and tangent to the crack
363  surface is defined as follows:
364  \f[
365  A^h_{\Gamma} =
366  \| \mathbf{N}(\tilde{\mathbf{X}}) \| =
367  \left\|
368  \epsilon_{ijk}
369  \frac{\partial \Phi^\alpha_p}{\partial \xi_i}
370  \frac{\partial \Phi^\beta_r}{\partial \xi_j}
371  \tilde{X}^\alpha_p
372  \tilde{X}^\beta_r
373  \right\|
374  \f]
375  *
376  *
377  */
378  struct OpCalculateGriffithForce
379  : public FaceElementForcesAndSourcesCore::UserDataOperator,
380  AuxOp {
381  OpCalculateGriffithForce(int tag, BlockData &block_data,
382  CommonData &common_data)
383  : FaceElementForcesAndSourcesCore::UserDataOperator(
384  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
385  AuxOp(tag, block_data, common_data) {}
386  AuxFunctions<double> auxFun;
387  MoFEMErrorCode doWork(int side, EntityType type,
388  DataForcesAndSourcesCore::EntData &data) {
389  MoFEMFunctionBegin;
390 
391  if (type != MBVERTEX)
392  MoFEMFunctionReturnHot(0);
393  CHKERR setIndices(data);
394  CHKERR setVariables(this, data);
395  int nb_dofs = data.getIndices().size();
396  int nb_gauss_pts = data.getN().size1();
397 
398  auxFun.referenceCoords.resize(nb_dofs, false);
399  auxFun.currentCoords.resize(nb_dofs, false);
400 
401  for (int dd = 0; dd != nb_dofs; dd++) {
402  auxFun.referenceCoords[dd] = getCoords()[dd];
403  auxFun.currentCoords[dd] = data.getFieldData()[dd];
404  }
405 
406  for (int gg = 0; gg != nb_gauss_pts; gg++) {
407  double val = getGaussPts()(2, gg) * 0.5;
408  CHKERR auxFun.calculateGriffithForce(blockData.gc, val,
409  data.getDiffN(gg));
410  }
411  commonData.griffithForce.resize(nb_dofs, false);
412  std::fill(commonData.griffithForce.begin(),
413  commonData.griffithForce.end(), 0);
414  std::copy(auxFun.griffithForce.begin(), auxFun.griffithForce.end(),
415  commonData.griffithForce.begin());
416 
417  MoFEMFunctionReturn(0);
418  }
419  };
420 
421  /**
422  * \brief Calculate density at the crack fron nodes and multiplity gc
423  *
424  * For heterogeneous case, gc is multiplied by a coefficient depending on
425  spatial distribution of density as follows:
426  *
427  \f[
428  g_c(\mathbf X) = g_c \cdot \rho(\mathbf X)^{\beta^{gc}}
429  \f]
430 
431  */
432  struct OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs
433  : public FaceElementForcesAndSourcesCore::UserDataOperator,
434  AuxOp {
435  const double betaGC;
436  std::string rhoTagName;
437  boost::shared_ptr<MWLSApprox> mwlsApprox;
438  MoFEM::Interface &mField;
439 
440  OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs(
441  string rho_tag_name, const double beta_gc,
442  boost::shared_ptr<MWLSApprox> mwls_approx, MoFEM::Interface &m_field,
443  int tag, BlockData &block_data, CommonData &common_data)
444  : FaceElementForcesAndSourcesCore::UserDataOperator(
445  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
446  rhoTagName(rho_tag_name), betaGC(beta_gc), mwlsApprox(mwls_approx),
447  mField(m_field), AuxOp(tag, block_data, common_data) {}
448 
449  MoFEMErrorCode doWork(int side, EntityType type,
450  DataForcesAndSourcesCore::EntData &row_data) {
451  MoFEMFunctionBegin;
452  FTensor::Index<'i', 3> i;
453  FTensor::Index<'j', 3> j;
454 
455  if (type != MBVERTEX)
456  MoFEMFunctionReturnHot(0);
457  const int row_nb_dofs = row_data.getIndices().size();
458  if (!row_nb_dofs)
459  MoFEMFunctionReturnHot(0);
460 
461  CHKERR setIndices(row_data);
462  CHKERR setVariables(this, row_data);
463 
464  auto &inv_AB_map = mwlsApprox->invABMap.at(getFEEntityHandle());
465  auto &influence_nodes_map =
466  mwlsApprox->influenceNodesMap.at(getFEEntityHandle());
467  auto &dm_nodes_map = mwlsApprox->dmNodesMap[this->getFEEntityHandle()];
468 
469  auto get_ftensor_from_vec = [](auto &v) {
470  return FTensor::Tensor1<double *, 3>(&v(0), &v(1), &v(2), 3);
471  };
472 
473  Tag th;
474  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(rhoTagName.c_str(), th);
475 
476  int nb_dofs = row_data.getFieldData().size();
477 
478  commonData.diffRho.resize(3, nb_dofs / 3, false);
479  std::fill(commonData.diffRho.data().begin(),
480  commonData.diffRho.data().end(), 0);
481 
482  auto t_coords = get_ftensor_from_vec(row_data.getFieldData());
483 
484  commonData.densityRho.resize(nb_dofs / 3, false);
485  std::fill(commonData.densityRho.begin(), commonData.densityRho.end(), 0);
486  auto rho = getFTensor0FromVec(commonData.densityRho);
487  VectorDouble gc_coeff(nb_dofs, false);
488  auto t_coeff = get_ftensor_from_vec(gc_coeff);
489  auto t_approx_base_point =
490  getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
491  const int nb_gauss_pts = row_data.getN().size1();
492 
493  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
494  for (int dd = 0; dd != nb_dofs / 3; ++dd) {
495 
496  // go only through nodes at the crack front
497  if (rowIndices[dd * 3] != -1) {
498  FTensor::Tensor1<double, 3> t_pos;
499  t_pos(i) = t_coords(i);
500 
501  for (int dd : {0, 1, 2})
502  mwlsApprox->approxPointCoords[dd] = t_approx_base_point(dd);
503 
504  mwlsApprox->getInfluenceNodes() = influence_nodes_map[gg];
505  mwlsApprox->getInvAB() = inv_AB_map[gg];
506  mwlsApprox->getShortenDm() = dm_nodes_map[gg];
507 
508  // approximate at MESH_NODE_POSITIONS
509  CHKERR mwlsApprox->evaluateApproxFun(&t_pos(0));
510  CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_pos(0), false);
511  CHKERR mwlsApprox->getTagData(th);
512  rho = mwlsApprox->getDataApprox()[0];
513  const auto &diff_vals = mwlsApprox->getDiffDataApprox();
514  for (int ii = 0; ii != 3; ++ii) {
515  commonData.diffRho(ii, dd) = diff_vals(ii, 0);
516  }
517  t_coeff(i) = pow(rho, betaGC);
518  }
519 
520  ++rho;
521  ++t_coords;
522  ++t_coeff;
523  }
524  }
525  // multiply griffith force with coefficients
526  commonData.griffithForce =
527  element_prod(commonData.griffithForce, gc_coeff);
528  MoFEMFunctionReturn(0);
529  }
530  };
531 
532  struct SaveGriffithForceOnTags {
533 
534  MoFEM::Interface &mField;
535 
536  Tag tH;
537  Range &tRis;
538  Tag *thCoordsPtr;
539 
540  SaveGriffithForceOnTags(MoFEM::Interface &m_field, Tag th, Range &tris,
541  Tag *th_coords = nullptr)
542  : mField(m_field), tH(th), tRis(tris), thCoordsPtr(th_coords) {}
543 
544  AuxFunctions<double> auxFun;
545 
546  MoFEMErrorCode operator()() {
547  MoFEMFunctionBegin;
548 
549  Range verts;
550  CHKERR mField.get_moab().get_connectivity(tRis, verts, true);
551  {
552  std::vector<double> data(verts.size() * 3, 0);
553  CHKERR mField.get_moab().tag_set_data(tH, verts, &*data.begin());
554  }
555  MatrixDouble diff_n(3, 2);
556  CHKERR ShapeDiffMBTRI(&*diff_n.data().begin());
557 
558  for (Range::iterator tit = tRis.begin(); tit != tRis.end(); ++tit) {
559  int num_nodes;
560  const EntityHandle *conn;
561  CHKERR mField.get_moab().get_connectivity(*tit, conn, num_nodes, true);
562  VectorDouble9 coords(9);
563 
564  int nb_dofs = num_nodes * 3;
565  auxFun.referenceCoords.resize(9, false);
566  auxFun.currentCoords.resize(9, false);
567 
568  if (thCoordsPtr) {
569  CHKERR mField.get_moab().tag_get_data(*thCoordsPtr, conn, num_nodes,
570  &*coords.begin());
571  } else {
572  CHKERR mField.get_moab().get_coords(conn, num_nodes,
573  &*coords.begin());
574  }
575  for (int dd = 0; dd != nb_dofs; dd++) {
576  auxFun.currentCoords[dd] = coords[dd];
577  }
578  if (thCoordsPtr)
579  CHKERR mField.get_moab().get_coords(conn, num_nodes,
580  &*coords.begin());
581  for (int dd = 0; dd != nb_dofs; dd++) {
582  auxFun.referenceCoords[dd] = coords[dd];
583  }
584 
585  CHKERR auxFun.calculateGriffithForce(1, 0.5, diff_n);
586 
587  double *tag_vals[3];
588  CHKERR mField.get_moab().tag_get_by_ptr(tH, conn, 3,
589  (const void **)tag_vals);
590  for (int nn = 0; nn != 3; ++nn) {
591  for (int dd = 0; dd != 3; ++dd) {
592  tag_vals[nn][dd] += auxFun.griffithForce[3 * nn + dd];
593  }
594  }
595  }
596 
597  MoFEMFunctionReturn(0);
598  }
599  };
600 
601  struct OpLhs : public FaceElementForcesAndSourcesCore::UserDataOperator,
602  AuxOp {
603 
604  OpLhs(int tag, BlockData &block_data, CommonData &common_data,
605  const double beta_gc = 0)
606  : FaceElementForcesAndSourcesCore::UserDataOperator(
607  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
608  UserDataOperator::OPROWCOL),
609  betaGC(beta_gc), AuxOp(tag, block_data, common_data) {}
610 
611  MatrixDouble mat;
612  const double betaGC;
613 
614  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
615  EntityType col_type,
616  DataForcesAndSourcesCore::EntData &row_data,
617  DataForcesAndSourcesCore::EntData &col_data) {
618  MoFEMFunctionBegin;
619 
620  if (row_type != MBVERTEX) {
621  MoFEMFunctionReturnHot(0);
622  }
623 
624  FTensor::Index<'i', 3> i;
625  FTensor::Index<'j', 3> j;
626  FTensor::Index<'l', 3> l;
627  FTensor::Index<'m', 3> m;
628  FTensor::Index<'k', 3> k;
629  FTensor::Number<0> N0;
630  FTensor::Number<1> N1;
631 
632  AuxFunctions<double> auxFun;
633  auxFun.currentCoords.resize(9, false);
634 
635  int nb_dofs = row_data.getFieldData().size();
636  for (int dd = 0; dd != nb_dofs; dd++)
637  auxFun.currentCoords[dd] = row_data.getFieldData()[dd];
638 
639  CHKERR setIndices(row_data);
640  CHKERR setVariables(this, row_data);
641  int nb_gauss_pts = row_data.getN().size1();
642  int nb_base_row = row_data.getFieldData().size() / 3;
643  int nb_base_col = col_data.getFieldData().size() / 3;
644  int row_nb_dofs = row_data.getIndices().size();
645 
646  mat.resize(9, 9, false);
647  mat.clear();
648  auto calculate_derivative = [&](FTensor::Tensor1<double *, 2> &t_diff) {
649  FTensor::Tensor2<double, 3, 3> t_d_n;
650  t_d_n(i, j) = FTensor::levi_civita(i, j, k) * auxFun.t_Tangent1_ksi(k) *
651  t_diff(N1);
652  t_d_n(i, j) -= FTensor::levi_civita(i, j, k) *
653  auxFun.t_Tangent2_eta(k) * t_diff(N0);
654 
655  return t_d_n;
656  };
657 
658  CHKERR auxFun.calculateAreaAndNormal(col_data.getDiffN());
659  auto t_normal = auxFun.t_Normal;
660  const double coefficient_1 = 0.5 * pow((t_normal(i) * t_normal(i)), -0.5);
661  const double coefficient_2 = 0.5 * pow((t_normal(i) * t_normal(i)), -1.5);
662 
663  // for homogeneous case, this is a bit ugly
664  auto &rho_v = commonData.densityRho;
665  if (rho_v.empty() || rho_v.size() != nb_base_row) {
666  rho_v.resize(nb_base_row, false);
667  }
668  auto rho = getFTensor0FromVec(rho_v);
669  for (int gg = 0; gg != nb_gauss_pts; gg++) {
670 
671  double val = getGaussPts()(2, gg) * 0.5;
672  auto t_row_diff = row_data.getFTensor1DiffN<2>(gg, 0);
673  for (int rrr = 0; rrr != nb_base_row; rrr++) {
674 
675  FTensor::Tensor2<double *, 3, 3> t_mat(
676  &mat(3 * rrr + 0, 0), &mat(3 * rrr + 0, 1), &mat(3 * rrr + 0, 2),
677  &mat(3 * rrr + 1, 0), &mat(3 * rrr + 1, 1), &mat(3 * rrr + 1, 2),
678  &mat(3 * rrr + 2, 0), &mat(3 * rrr + 2, 1), &mat(3 * rrr + 2, 2),
679  3);
680 
681  auto t_tan_row = calculate_derivative(t_row_diff);
682  auto t_col_diff = col_data.getFTensor1DiffN<2>(gg, 0);
683  const double coef = blockData.gc * val * pow(rho, betaGC);
684  for (int ccc = 0; ccc != nb_base_col; ccc++) {
685 
686  auto t_tan_col = calculate_derivative(t_col_diff);
687 
688  t_mat(i, j) -= coefficient_1 *
689  (t_tan_row(i, l) * t_tan_col(l, j) +
690  FTensor::levi_civita(i, j, k) * t_col_diff(N0) *
691  t_row_diff(N1) * t_normal(k) -
692  FTensor::levi_civita(i, j, k) * t_col_diff(N1) *
693  t_row_diff(N0) * t_normal(k));
694 
695  t_mat(i, j) += coefficient_2 * (t_tan_row(i, k) * t_normal(k)) *
696  (t_tan_col(l, j) * t_normal(l));
697 
698  t_mat(i, j) *= coef;
699  ++t_col_diff;
700  ++t_mat;
701  }
702  ++t_row_diff;
703  ++rho;
704  }
705  }
706 
707  int col_nb_dofs = col_data.getIndices().size();
708 
709  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
710  &*rowIndices.data().begin(), col_nb_dofs,
711  &*col_data.getIndices().data().begin(),
712  &*mat.data().begin(), ADD_VALUES);
713 
714  MoFEMFunctionReturn(0);
715  }
716  };
717 
718  struct OpHeterogeneousGcLhs
719  : public FaceElementForcesAndSourcesCore::UserDataOperator,
720  AuxOp {
721 
722  OpHeterogeneousGcLhs(int tag, BlockData &block_data,
723  CommonData &common_data, const double beta_gc)
724  : FaceElementForcesAndSourcesCore::UserDataOperator(
725  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
726  UserDataOperator::OPROWCOL),
727  betaGC(beta_gc), AuxOp(tag, block_data, common_data) {}
728 
729  MatrixDouble mat;
730  const double betaGC;
731 
732  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
733  EntityType col_type,
734  DataForcesAndSourcesCore::EntData &row_data,
735  DataForcesAndSourcesCore::EntData &col_data) {
736  MoFEMFunctionBegin;
737 
738  if (row_type != MBVERTEX) {
739  MoFEMFunctionReturnHot(0);
740  }
741 
742  FTensor::Index<'i', 3> i;
743  FTensor::Index<'j', 3> j;
744 
745  AuxFunctions<double> auxFun;
746  auxFun.currentCoords.resize(9, false);
747 
748  int nb_dofs = row_data.getFieldData().size();
749  for (int dd = 0; dd != nb_dofs; dd++)
750  auxFun.currentCoords[dd] = row_data.getFieldData()[dd];
751 
752  CHKERR setIndices(row_data);
753  CHKERR setVariables(this, row_data);
754  int nb_gauss_pts = row_data.getN().size1();
755  int nb_base_row = row_data.getFieldData().size() / 3;
756  int nb_base_col = col_data.getFieldData().size() / 3;
757 
758  int row_nb_dofs = row_data.getIndices().size();
759  mat.resize(9, 9, false);
760  mat.clear();
761 
762  const double &gc = blockData.gc;
763  auto rho = getFTensor0FromVec(commonData.densityRho);
764  auto get_ftensor_from_vec = [](auto &v) {
765  return FTensor::Tensor1<double *, 3>(&v(0), &v(1), &v(2), 3);
766  };
767 
768  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
769  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
770  &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
771  &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
772  &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
773  &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
774  &m(3 * r + 2, 3 * c + 2));
775  };
776 
777  for (int gg = 0; gg != nb_gauss_pts; gg++) {
778 
779  double val = getGaussPts()(2, gg) * 0.5;
780 
781  auxFun.calculateGriffithForce(1, 1, row_data.getDiffN(gg));
782 
783  auto t_griffith = get_ftensor_from_vec(auxFun.griffithForce);
784  auto t_diff_rho = getFTensor1FromMat<3>(commonData.diffRho);
785  for (int rrr = 0; rrr != nb_base_row; rrr++) {
786 
787  const double a = val * gc * pow(rho, betaGC - 1.) * betaGC;
788  FTensor::Tensor2<double, 3, 3> k;
789  k(i, j) = t_diff_rho(j) * t_griffith(i) * a;
790  // there is no dependence on the other nodes so row == col
791  auto tt_mat = get_tensor2(mat, rrr, rrr);
792  tt_mat(i, j) += k(i, j);
793 
794  ++t_diff_rho;
795  ++t_griffith;
796  ++rho;
797  }
798  }
799 
800  int col_nb_dofs = col_data.getIndices().size();
801 
802  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
803  &*rowIndices.data().begin(), col_nb_dofs,
804  &*col_data.getIndices().data().begin(),
805  &*mat.data().begin(), ADD_VALUES);
806 
807  MoFEMFunctionReturn(0);
808  }
809  };
810 
811  struct MyTriangleFEConstrainsDelta : public MyTriangleFE {
812 
813  const std::string lambdaFieldName;
814  bool ownDeltaVec;
815  SmartPetscObj<Vec> deltaVec;
816 
817  MyTriangleFEConstrainsDelta(MoFEM::Interface &m_field,
818  const std::string &lambda_field_name)
819  : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name) {}
820 
821  MoFEMErrorCode preProcess() {
822  MoFEMFunctionBegin;
823  Vec l;
824  CHKERR VecGhostGetLocalForm(deltaVec, &l);
825  double *a;
826  CHKERR VecGetArray(l, &a);
827  int n;
828  CHKERR VecGetLocalSize(l, &n);
829  std::fill(&a[0], &a[n], 0);
830  CHKERR VecRestoreArray(l, &a);
831  CHKERR VecGhostRestoreLocalForm(deltaVec, &l);
832  MoFEMFunctionReturn(0);
833  }
834 
835  MoFEMErrorCode postProcess() {
836  MoFEMFunctionBegin;
837  MOFEM_LOG_CHANNEL("WORLD");
838  MOFEM_LOG_FUNCTION();
839  MOFEM_LOG_TAG("WORLD", "Arc");
840 
841  CHKERR VecAssemblyBegin(deltaVec);
842  CHKERR VecAssemblyEnd(deltaVec);
843  CHKERR VecGhostUpdateBegin(deltaVec, ADD_VALUES, SCATTER_REVERSE);
844  CHKERR VecGhostUpdateEnd(deltaVec, ADD_VALUES, SCATTER_REVERSE);
845  CHKERR VecGhostUpdateBegin(deltaVec, INSERT_VALUES, SCATTER_FORWARD);
846  CHKERR VecGhostUpdateEnd(deltaVec, INSERT_VALUES, SCATTER_FORWARD);
847  double nrm;
848  CHKERR VecNorm(deltaVec, NORM_2, &nrm);
849  MOFEM_LOG_C("WORLD", Sev::noisy, "\tdelta vec nrm = %9.8e", nrm);
850  MoFEMFunctionReturn(0);
851  }
852  };
853 
854  struct OpConstrainsDelta
855  : public FaceElementForcesAndSourcesCore::UserDataOperator,
856  AuxOp {
857 
858  MoFEM::Interface &mField;
859  const std::string lambdaFieldName;
860  SmartPetscObj<Vec> &deltaVec;
861 
862  OpConstrainsDelta(MoFEM::Interface &m_field, int tag, BlockData &block_data,
863  CommonData &common_data,
864  const std::string &lambda_field_name,
865  SmartPetscObj<Vec> &delta_vec)
866  : FaceElementForcesAndSourcesCore::UserDataOperator(
867  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
868  AuxOp(tag, block_data, common_data), mField(m_field),
869  lambdaFieldName(lambda_field_name), deltaVec(delta_vec) {}
870 
871  AuxFunctions<double> auxFun;
872 
873  MoFEMErrorCode doWork(int side, EntityType type,
874  DataForcesAndSourcesCore::EntData &data) {
875  MoFEMFunctionBegin;
876  if (type != MBVERTEX) {
877  MoFEMFunctionReturnHot(0);
878  }
879  // get indices of material DOFs
880  CHKERR setIndices(data);
881  // get indices of Lagrange multiplier DOFs
882  CHKERR setLambdaIndices(this, lambdaFieldName);
883  // do calculations
884  auxFun.currentCoords.resize(9, false);
885  for (int dd = 0; dd != 9; ++dd)
886  auxFun.currentCoords[dd] = data.getFieldData()[dd];
887 
888  CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
889  data.getDiffN(0));
890 
891  // assemble area change
892  double delta[] = {0, 0, 0};
893  for (int nn = 0; nn != 3; ++nn) {
894  if (rowLambdaIndices[nn] == -1) {
895  if (rowIndices[3 * nn + 0] != -1 || rowIndices[3 * nn + 1] != -1 ||
896  rowIndices[3 * nn + 2] != -1) {
897  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
898  "Element has no nodes at crack front");
899  }
900  continue;
901  }
902  for (int dd = 0; dd != 3; ++dd) {
903  int idx = 3 * nn + dd;
904  delta[nn] += auxFun.griffithForce[idx] *
905  (data.getFieldData()[idx] - getCoords()[idx]);
906  }
907  }
908  // add nodes to vector of lambda indices
909  CHKERR VecSetValues(deltaVec, 3, &rowLambdaIndices[0], delta, ADD_VALUES);
910 
911  MoFEMFunctionReturn(0);
912  }
913  };
914 
915  struct MyTriangleFEConstrains : public MyTriangleFE {
916 
917  const std::string lambdaFieldName;
918  BlockData &blockData;
919  SmartPetscObj<Vec> deltaVec;
920 
921  MyTriangleFEConstrains(MoFEM::Interface &m_field,
922  const std::string &lambda_field_name,
923  BlockData &block_data, SmartPetscObj<Vec> &delta_vec)
924  : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name),
925  blockData(block_data), deltaVec(delta_vec) {}
926 
927  MoFEMErrorCode preProcess() {
928  MoFEMFunctionBegin;
929 
930  MOFEM_LOG_CHANNEL("WORLD");
931  MOFEM_LOG_FUNCTION();
932  MOFEM_LOG_TAG("WORLD", "GriffithConstrain");
933 
934  const double *delta;
935  CHKERR VecGetArrayRead(deltaVec, &delta);
936  double sum_of_delta = 0;
937  double sum_of_lambda = 0;
938  const auto bit_field_number =
939  mField.get_field_bit_number(lambdaFieldName);
940  switch (snes_ctx) {
941  case CTX_SNESSETFUNCTION: {
942  for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
943  bit_field_number, dit)) {
944  if (static_cast<int>(dit->get()->getPart()) ==
945  mField.get_comm_rank()) {
946  int local_idx = dit->get()->getPetscLocalDofIdx();
947  double lambda = dit->get()->getFieldData();
948  double W = lambda - blockData.E * delta[local_idx];
949  double val = lambda - calMax(W, 0, blockData.r);
950  int global_idx = dit->get()->getPetscGlobalDofIdx();
951  CHKERR VecSetValue(snes_f, global_idx, val, ADD_VALUES);
952  }
953  }
954  } break;
955  case CTX_SNESSETJACOBIAN: {
956  for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
957  bit_field_number, dit)) {
958  if (static_cast<int>(dit->get()->getPart()) ==
959  mField.get_comm_rank()) {
960  int local_idx = dit->get()->getPetscLocalDofIdx();
961  double lambda = dit->get()->getFieldData();
962  double W = lambda - blockData.E * delta[local_idx];
963  double diffW = 1;
964  double val = 1 - diffCalMax_a(W, 0, blockData.r) * diffW;
965  int global_idx = dit->get()->getPetscGlobalDofIdx();
966  MOFEM_LOG_C("WORLD", Sev::verbose,
967  "Constrains on node %lu diag = %+3.5e "
968  "delta = %+3.5e lambda = %+3.5e",
969  dit->get()->getEnt(), val, delta[local_idx], lambda);
970  CHKERR MatSetValue(snes_B, global_idx, global_idx, val, ADD_VALUES);
971  sum_of_delta += delta[local_idx];
972  sum_of_lambda += lambda;
973  }
974  }
975  MOFEM_LOG_C("WORLD", Sev::noisy,
976  "Sum of delta = %+6.4e Sum of lambda = %+6.4e",
977  sum_of_delta, sum_of_lambda);
978 
979  } break;
980  default:
981  break;
982  }
983  CHKERR VecRestoreArrayRead(deltaVec, &delta);
984  MoFEMFunctionReturn(0);
985  }
986 
987  MoFEMErrorCode postProcess() {
988  MoFEMFunctionBeginHot;
989  MoFEMFunctionReturnHot(0);
990  }
991  };
992 
993  struct OpConstrainsRhs
994  : public FaceElementForcesAndSourcesCore::UserDataOperator,
995  AuxOp {
996 
997  MoFEM::Interface &mField;
998  const std::string lambdaFieldName;
999 
1000  OpConstrainsRhs(MoFEM::Interface &m_field, int tag, BlockData &block_data,
1001  CommonData &common_data,
1002  const std::string &lambda_field_name)
1003  : FaceElementForcesAndSourcesCore::UserDataOperator(
1004  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
1005  AuxOp(tag, block_data, common_data), mField(m_field),
1006  lambdaFieldName(lambda_field_name) {}
1007 
1008  AuxFunctions<double> auxFun;
1009  VectorDouble9 nF;
1010  // VectorDouble3 nG;
1011 
1012  MoFEMErrorCode doWork(int side, EntityType type,
1013  DataForcesAndSourcesCore::EntData &data) {
1014  MoFEMFunctionBegin;
1015  if (type != MBVERTEX)
1016  MoFEMFunctionReturnHot(0);
1017 
1018  // get indices of material DOFs
1019  CHKERR setIndices(data);
1020  CHKERR setLambdaNodes(this, lambdaFieldName);
1021  CHKERR setLambdaIndices(this, lambdaFieldName);
1022 
1023  auxFun.currentCoords.resize(9, false);
1024  noalias(auxFun.currentCoords) = data.getFieldData();
1025 
1026  CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1027  data.getDiffN(0));
1028 
1029  // assemble local vector
1030  nF.resize(9, false);
1031  nF.clear();
1032  for (int nn = 0; nn != 3; nn++) {
1033  if (rowLambdaIndices[nn] != -1) {
1034  const auto lambda = lambdaAtNodes[nn];
1035  for (int dd = 0; dd != 3; dd++) {
1036  int idx = 3 * nn + dd;
1037  nF[idx] -= lambda * auxFun.griffithForce[idx];
1038  }
1039  }
1040  }
1041 
1042  // assemble the right hand vectors
1043  CHKERR VecSetValues(getFEMethod()->snes_f, 9, &*rowIndices.data().begin(),
1044  &nF[0], ADD_VALUES);
1045 
1046  MoFEMFunctionReturn(0);
1047  }
1048  };
1049 
1050  struct OpConstrainsLhs
1051  : public FaceElementForcesAndSourcesCore::UserDataOperator,
1052  AuxOp {
1053 
1054  MoFEM::Interface &mField;
1055  SmartPetscObj<Vec> &deltaVec;
1056 
1057  OpConstrainsLhs(MoFEM::Interface &m_field, int tag, BlockData &block_data,
1058  CommonData &common_data,
1059  const std::string &lambda_field_name,
1060  SmartPetscObj<Vec> &delta_vec)
1061  : FaceElementForcesAndSourcesCore::UserDataOperator(
1062  "MESH_NODE_POSITIONS", lambda_field_name,
1063  UserDataOperator::OPROWCOL,
1064  false // not symmetric operator
1065  ),
1066  AuxOp(tag, block_data, common_data), mField(m_field),
1067  deltaVec(delta_vec) {}
1068 
1069  MatrixDouble nG_dX;
1070  MatrixDouble nF_dLambda;
1071  MatrixDouble nF_dX;
1072 
1073  MatrixDouble jacDelta;
1074  ublas::vector<double *> jacDeltaPtr;
1075  VectorDouble d_dA;
1076  VectorDouble d_Delta;
1077 
1078  ublas::vector<adouble> a_Delta;
1079  AuxFunctions<adouble> a_auxFun;
1080  ublas::vector<adouble> a_dA;
1081  ublas::vector<adouble> a_lambdaAtNodes;
1082 
1083  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1084  EntityType col_type,
1085  DataForcesAndSourcesCore::EntData &row_data,
1086  DataForcesAndSourcesCore::EntData &col_data) {
1087  MoFEMFunctionBegin;
1088 
1089  if (row_type != MBVERTEX || col_type != MBVERTEX)
1090  MoFEMFunctionReturnHot(0);
1091 
1092  CHKERR setIndices(row_data);
1093  CHKERR setLambdaNodes(this, colFieldName);
1094  CHKERR setLambdaIndices(this, colFieldName);
1095 
1096  activeVariables.resize(18, false);
1097  a_auxFun.referenceCoords.resize(9, false);
1098  a_auxFun.currentCoords.resize(9, false);
1099 
1100  a_Delta.resize(3, false);
1101  d_Delta.resize(3, false);
1102  trace_on(tAg);
1103  for (int dd = 0; dd != 9; dd++) {
1104  double val = getCoords()[dd];
1105  a_auxFun.referenceCoords[dd] <<= val;
1106  activeVariables[dd] = val;
1107  }
1108  for (int dd = 0; dd != 9; dd++) {
1109  double val = row_data.getFieldData()[dd];
1110  a_auxFun.currentCoords[dd] <<= val;
1111  activeVariables[9 + dd] = val;
1112  }
1113 
1114  CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1115  row_data.getDiffN(0));
1116 
1117  // calculate area change
1118  std::fill(a_Delta.begin(), a_Delta.end(), 0);
1119  for (int nn = 0; nn != 3; nn++) {
1120  if (col_data.getIndices()[nn] == -1)
1121  continue;
1122  if (col_data.getIndices()[nn] != rowLambdaIndices[nn]) {
1123  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1124  "data inconsistency");
1125  }
1126  for (int dd = 0; dd != 3; dd++) {
1127  int idx = 3 * nn + dd;
1128  if (row_data.getFieldDofs()[idx]->getEnt() !=
1129  col_data.getFieldDofs()[nn]->getEnt()) {
1130  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1131  "data inconsistency");
1132  }
1133  a_Delta[nn] +=
1134  a_auxFun.griffithForce[idx] *
1135  (a_auxFun.currentCoords[idx] - a_auxFun.referenceCoords[idx]);
1136  }
1137  }
1138  for (int nn = 0; nn != 3; nn++) {
1139  a_Delta[nn] >>= d_Delta[nn];
1140  }
1141  trace_off();
1142 
1143  jacDelta.resize(3, 18, false);
1144  jacDeltaPtr.resize(3, false);
1145  for (int nn = 0; nn != 3; nn++) {
1146  jacDeltaPtr[nn] = &jacDelta(nn, 0);
1147  }
1148 
1149  {
1150  int r;
1151  // play recorder for jacobian
1152  r = ::jacobian(tAg, 3, 18, &activeVariables[0], &jacDeltaPtr[0]);
1153  if (r < 3) { // function is locally analytic
1154  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1155  "ADOL-C function evaluation with error r = %d", r);
1156  }
1157  }
1158 
1159  const double *delta;
1160  Vec delta_vec_local;
1161  CHKERR VecGhostGetLocalForm(deltaVec, &delta_vec_local);
1162  int vec_size;
1163  CHKERR VecGetSize(delta_vec_local, &vec_size);
1164  CHKERR VecGetArrayRead(delta_vec_local, &delta);
1165  nG_dX.resize(3, 9, false);
1166  nG_dX.clear();
1167  for (int nn = 0; nn != 3; ++nn) {
1168  int local_idx = col_data.getLocalIndices()[nn];
1169  if (local_idx != rowLambdaIndicesLocal[nn]) {
1170  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1171  "data inconsistency");
1172  }
1173  if (local_idx == -1)
1174  continue;
1175  if (vec_size < local_idx) {
1176  int g_vec_size;
1177  CHKERR VecGetLocalSize(deltaVec, &g_vec_size);
1178  SETERRQ3(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1179  "data inconsistency %d < %d (%d)", vec_size, local_idx,
1180  g_vec_size);
1181  }
1182  double W = col_data.getFieldData()[nn] - blockData.E * delta[local_idx];
1183  double diff_constrain = -diffCalMax_a(W, 0, blockData.r);
1184  for (int dd = 0; dd != 9; dd++) {
1185  double diffW = -blockData.E * jacDelta(nn, 9 + dd);
1186  nG_dX(nn, dd) = diff_constrain * diffW;
1187  }
1188  }
1189  CHKERR VecRestoreArrayRead(delta_vec_local, &delta);
1190  CHKERR VecGhostRestoreLocalForm(deltaVec, &delta_vec_local);
1191  CHKERR MatSetValues(getFEMethod()->snes_B, 3,
1192  &*col_data.getIndices().data().begin(), 9,
1193  &*row_data.getIndices().data().begin(),
1194  &*nG_dX.data().begin(), ADD_VALUES);
1195 
1196  // Dervatives of nF
1197 
1198  // Linearisation with adol-C
1199  activeVariables.resize(21, false);
1200  d_dA.resize(9, false);
1201  a_dA.resize(9, false);
1202  a_lambdaAtNodes.resize(3, false);
1203  a_auxFun.referenceCoords.resize(9, false);
1204 
1205  trace_on(tAg);
1206  for (int dd = 0; dd != 9; dd++) {
1207  double val = getCoords()[dd];
1208  a_auxFun.referenceCoords[dd] <<= val;
1209  activeVariables[dd] = val;
1210  }
1211  for (int dd = 0; dd != 9; dd++) {
1212  double val = row_data.getFieldData()[dd];
1213  a_auxFun.currentCoords[dd] <<= val;
1214  activeVariables[9 + dd] = val;
1215  }
1216  for (int nn = 0; nn != 3; nn++) {
1217  double val = col_data.getFieldData()[nn];
1218  a_lambdaAtNodes[nn] <<= val;
1219  activeVariables[18 + nn] = val;
1220  }
1221 
1222  CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1223  row_data.getDiffN(0));
1224 
1225  std::fill(a_dA.begin(), a_dA.end(), 0);
1226  for (int nn = 0; nn != 3; nn++) {
1227  if (col_data.getIndices()[nn] == -1)
1228  continue;
1229  for (int dd = 0; dd != 3; dd++) {
1230  int idx = 3 * nn + dd;
1231  a_dA[idx] -= a_lambdaAtNodes[nn] * a_auxFun.griffithForce[idx];
1232  }
1233  }
1234  for (int dd = 0; dd != 9; dd++) {
1235  a_dA[dd] >>= d_dA[dd];
1236  }
1237  trace_off();
1238 
1239  commonData.tangentGriffithForce.resize(9, 21, false);
1240  commonData.tangentGriffithForceRowPtr.resize(9);
1241  for (int dd = 0; dd != 9; ++dd) {
1242  commonData.tangentGriffithForceRowPtr[dd] =
1243  &commonData.tangentGriffithForce(dd, 0);
1244  }
1245 
1246  {
1247  int r;
1248  // play recorder for jacobian
1249  r = ::jacobian(tAg, 9, 21, &activeVariables[0],
1250  &commonData.tangentGriffithForceRowPtr[0]);
1251  if (r < 3) { // function is locally analytic
1252  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1253  "ADOL-C function evaluation with error r = %d", r);
1254  }
1255  }
1256 
1257  nF_dLambda.resize(9, 3);
1258  for (int rr = 0; rr != 9; rr++) {
1259  for (int cc = 0; cc != 3; cc++) {
1260  nF_dLambda(rr, cc) = commonData.tangentGriffithForce(rr, 18 + cc);
1261  }
1262  }
1263  CHKERR MatSetValues(getFEMethod()->snes_B, 9, &*rowIndices.data().begin(),
1264  3, &*col_data.getIndices().data().begin(),
1265  &*nF_dLambda.data().begin(), ADD_VALUES);
1266 
1267  nF_dX.resize(9, 9, false);
1268  for (int rr = 0; rr != 9; rr++) {
1269  for (int cc = 0; cc != 9; cc++) {
1270  nF_dX(rr, cc) = commonData.tangentGriffithForce(rr, 9 + cc);
1271  }
1272  }
1273  CHKERR MatSetValues(getFEMethod()->snes_B, 9, &*rowIndices.data().begin(),
1274  9, &*row_data.getIndices().data().begin(),
1275  &*nF_dX.data().begin(), ADD_VALUES);
1276 
1277  MoFEMFunctionReturn(0);
1278  }
1279  };
1280 
1281  /**
1282  * \brief arc-length element
1283  *
1284  * Calualte residual for arc length method amd vector db, which is derivative
1285  * of residual.
1286  *
1287  */
1288  struct FrontArcLengthControl : public FEMethod {
1289 
1290  int tAg;
1291  boost::shared_ptr<ArcLengthCtx> arcPtr;
1292  BlockData &blockData;
1293  CommonData &commonData;
1294  std::string lambdaFieldName;
1295 
1296  Vec ghostIntLambda;
1297  double intLambdaArray[2];
1298  double &intLambdaDelta;
1299  double &intLambdaLambda;
1300  double intLambda;
1301 
1302  boost::shared_ptr<MyTriangleFE> feIntLambda;
1303 
1304  /**
1305  * \brief assemble internal residual
1306  * @param block_data griffith block data
1307  * @param common_data griffith common data
1308  * @param lambda_field_name lagrange multipliers controling crack surface
1309  * area
1310  * @param arc_front reference to FrontArcLengthControl object
1311  */
1312  struct OpIntLambda : FaceElementForcesAndSourcesCore::UserDataOperator,
1313  GriffithForceElement::AuxOp {
1314 
1315  MoFEM::Interface &mField;
1316  std::string lambdaFieldName;
1317  FrontArcLengthControl &arcFront;
1318 
1319  OpIntLambda(MoFEM::Interface &m_field,
1320  GriffithForceElement::BlockData &block_data,
1321  GriffithForceElement::CommonData &common_data,
1322  const std::string &lambda_field_name,
1323  FrontArcLengthControl &arc_front)
1324  : FaceElementForcesAndSourcesCore::UserDataOperator(
1325  "MESH_NODE_POSITIONS", lambda_field_name,
1326  UserDataOperator::OPROW,
1327  false // not symmetric operator
1328  ),
1329  AuxOp(0, block_data, common_data), mField(m_field),
1330  lambdaFieldName(lambda_field_name), arcFront(arc_front) {}
1331 
1332  AuxFunctions<double> auxFun;
1333 
1334  MoFEMErrorCode doWork(int side, EntityType type,
1335  DataForcesAndSourcesCore::EntData &data) {
1336  MoFEMFunctionBegin;
1337  if (type != MBVERTEX) {
1338  MoFEMFunctionReturnHot(0);
1339  }
1340  // get indices of Lagrange multiplier DOFs
1341  CHKERR setLambdaIndices(this, lambdaFieldName);
1342  // get dofs increment
1343  const double *dx;
1344  CHKERR VecGetArrayRead(arcFront.arcPtr->dx, &dx);
1345  // get dofs at begining of load step
1346  const double *x0;
1347  CHKERR VecGetArrayRead(arcFront.arcPtr->x0, &x0);
1348  auxFun.referenceCoords.resize(9, false);
1349  auxFun.currentCoords.resize(9, false);
1350  for (int dd = 0; dd != 9; dd++) {
1351  int loc_idx = data.getLocalIndices()[dd];
1352  if (loc_idx != -1) {
1353  auxFun.referenceCoords[dd] = x0[loc_idx];
1354  auxFun.currentCoords[dd] = auxFun.referenceCoords[dd] + dx[loc_idx];
1355  } else {
1356  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
1357  }
1358  }
1359  CHKERR VecRestoreArrayRead(arcFront.arcPtr->dx, &dx);
1360  CHKERR VecRestoreArrayRead(arcFront.arcPtr->x0, &x0);
1361 
1362  // calculate crack area increment
1363  CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1364  data.getDiffN(0));
1365 
1366  // calculate area change
1367  double delta = 0;
1368  for (int nn = 0; nn != 3; nn++) {
1369  if (rowLambdaIndices[nn] == -1)
1370  continue;
1371  auto griffith_force =
1372  getVectorAdaptor(&(auxFun.griffithForce[3 * nn]), 3);
1373  // griffith_force /= norm_2(griffith_force);
1374  for (int dd = 0; dd != 3; dd++) {
1375  const int idx = 3 * nn + dd;
1376  delta += griffith_force[dd] *
1377  (auxFun.currentCoords[idx] - auxFun.referenceCoords[idx]);
1378  }
1379  }
1380  if (delta != delta) {
1381  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "nan value");
1382  }
1383  arcFront.intLambdaDelta += delta;
1384  MoFEMFunctionReturn(0);
1385  };
1386  };
1387 
1388  /**
1389  * \brief assemble internal residual derivative
1390  * @param block_data griffith block data
1391  * @param common_data griffith common data
1392  * @param lambda_field_name lagrange multipliers controling crack surface
1393  * area
1394  * @param arc_front reference to FrontArcLengthControl object
1395  */
1396  struct OpAssemble_db : public OpIntLambda {
1397  int tAg;
1398  OpAssemble_db(int tag, MoFEM::Interface &m_field,
1399  GriffithForceElement::BlockData &block_data,
1400  GriffithForceElement::CommonData &common_data,
1401  const std::string &lambda_field_name,
1402  FrontArcLengthControl &arc_front)
1403  : OpIntLambda(m_field, block_data, common_data, lambda_field_name,
1404  arc_front),
1405  tAg(tag) {}
1406 
1407  AuxFunctions<adouble> a_auxFun;
1408  VectorDouble gradDelta;
1409  VectorDouble activeVariables;
1410 
1411  MoFEMErrorCode doWork(int side, EntityType type,
1412  DataForcesAndSourcesCore::EntData &data) {
1413  MoFEMFunctionBeginHot;
1414  if (type != MBVERTEX) {
1415  MoFEMFunctionReturnHot(0);
1416  }
1417  CHKERR setIndices(data);
1418  // get indices of Lagrange multiplier DOFs
1419  CHKERR setLambdaIndices(this, lambdaFieldName);
1420  double d_delta;
1421  adouble a_delta;
1422  // get dofs increment
1423  const double *dx;
1424  CHKERR VecGetArrayRead(arcFront.arcPtr->dx, &dx);
1425  // get dofs at begining of load step
1426  const double *x0;
1427  CHKERR VecGetArrayRead(arcFront.arcPtr->x0, &x0);
1428  a_auxFun.referenceCoords.resize(9, false);
1429  a_auxFun.currentCoords.resize(9, false);
1430 
1431  activeVariables.resize(18, false);
1432  trace_on(tAg);
1433  for (int dd = 0; dd != 9; dd++) {
1434  int loc_idx = data.getLocalIndices()[dd];
1435  if (loc_idx != -1) {
1436  a_auxFun.referenceCoords[dd] <<= x0[loc_idx];
1437  activeVariables[dd] = x0[loc_idx];
1438  } else {
1439  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
1440  }
1441  }
1442  for (int dd = 0; dd != 9; dd++) {
1443  int loc_idx = data.getLocalIndices()[dd];
1444  if (loc_idx != -1) {
1445  a_auxFun.currentCoords[dd] <<= x0[dd] + dx[loc_idx];
1446  activeVariables[9 + dd] = x0[loc_idx] + dx[loc_idx];
1447  } else {
1448  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
1449  }
1450  }
1451  // calculate crack area increment
1452 
1453  CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1454  data.getDiffN(0));
1455  // calculate area change
1456  a_delta = 0;
1457  for (int nn = 0; nn != 3; nn++) {
1458  if (rowLambdaIndices[nn] == -1)
1459  continue;
1460  auto griffith_force =
1461  getVectorAdaptor(&(a_auxFun.griffithForce[3 * nn]), 3);
1462  // griffith_force /= sqrt(griffith_force[0] * griffith_force[0] +
1463  // griffith_force[1] * griffith_force[1] +
1464  // griffith_force[2] * griffith_force[2]);
1465  for (int dd = 0; dd != 3; dd++) {
1466  const int idx = 3 * nn + dd;
1467  a_delta += griffith_force[dd] * (a_auxFun.currentCoords[idx] -
1468  a_auxFun.referenceCoords[idx]);
1469  }
1470  }
1471  a_delta >>= d_delta;
1472  trace_off();
1473  CHKERR VecRestoreArrayRead(arcFront.arcPtr->dx, &dx);
1474  CHKERR VecRestoreArrayRead(arcFront.arcPtr->x0, &x0);
1475 
1476  gradDelta.resize(18, false);
1477  {
1478  int r;
1479  // play recorder for jacobian
1480  r = ::gradient(tAg, 18, &activeVariables[0], &gradDelta[0]);
1481  if (r < 3) { // function is locally analytic
1482  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1483  "ADOL-C function evaluation with error r = %d", r);
1484  }
1485  }
1486 
1487  for (int dd = 0; dd != 9; dd++) {
1488  double val = arcFront.arcPtr->alpha * gradDelta[9 + dd];
1489  CHKERR VecSetValue(arcFront.arcPtr->db, data.getIndices()[dd], val,
1490  ADD_VALUES);
1491  }
1492 
1493  MoFEMFunctionReturnHot(0);
1494  }
1495  };
1496 
1497  FrontArcLengthControl(int tag, GriffithForceElement::BlockData &block_data,
1498  GriffithForceElement::CommonData &common_data,
1499  const std::string &lambda_field_name,
1500  boost::shared_ptr<ArcLengthCtx> &arc_ptr)
1501  : tAg(tag), arcPtr(arc_ptr), blockData(block_data),
1502  commonData(common_data), lambdaFieldName(lambda_field_name),
1503  intLambdaDelta(intLambdaArray[0]),
1504  intLambdaLambda(intLambdaArray[1]) {
1505  int two[] = {0, 1};
1506  int nb_local = arcPtr->mField.get_comm_rank() == 0 ? 2 : 0;
1507  int nb_ghost = nb_local ? 0 : 2;
1508  ierr = VecCreateGhostWithArray(PETSC_COMM_WORLD, nb_local, 2, nb_ghost,
1509  two, intLambdaArray, &ghostIntLambda);
1510  CHKERRABORT(PETSC_COMM_WORLD, ierr);
1511  feIntLambda = boost::shared_ptr<MyTriangleFE>(
1512  new GriffithForceElement::MyTriangleFE(arcPtr->mField));
1513  }
1514  virtual ~FrontArcLengthControl() {
1515  ierr = VecDestroy(&ghostIntLambda);
1516  CHKERRABORT(PETSC_COMM_WORLD, ierr);
1517  }
1518 
1519  MoFEMErrorCode preProcess() {
1520  MoFEMFunctionBegin;
1521  switch (snes_ctx) {
1522  case CTX_SNESSETFUNCTION: {
1523  CHKERR VecAssemblyBegin(snes_f);
1524  CHKERR VecAssemblyEnd(snes_f);
1525  CHKERR calculateDxAndDlambda(snes_x);
1526  CHKERR calculateDb();
1527  } break;
1528  default:
1529  break;
1530  }
1531  MoFEMFunctionReturn(0);
1532  }
1533 
1534  /** \brief Calculate internal lambda
1535  */
1536  double calculateLambdaInt() {
1537  MoFEMFunctionBeginHot;
1538  return intLambda;
1539  MoFEMFunctionReturnHot(0);
1540  }
1541 
1542  MoFEMErrorCode operator()() {
1543  MoFEMFunctionBegin;
1544  MOFEM_LOG_CHANNEL("WORLD");
1545  MOFEM_LOG_FUNCTION();
1546  MOFEM_LOG_TAG("WORLD", "Arc");
1547 
1548  switch (snes_ctx) {
1549  case CTX_SNESSETFUNCTION: {
1550  arcPtr->res_lambda = calculateLambdaInt() - arcPtr->alpha * arcPtr->s;
1551  CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
1552  arcPtr->res_lambda, ADD_VALUES);
1553  } break;
1554  case CTX_SNESSETJACOBIAN: {
1555  arcPtr->dIag = 2 * arcPtr->dLambda * arcPtr->alpha *
1556  pow(arcPtr->beta, 2) * arcPtr->F_lambda2;
1557  MOFEM_LOG_C("WORLD", Sev::noisy, "diag = %9.8e", arcPtr->dIag);
1558  CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
1559  arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
1560  } break;
1561  default:
1562  break;
1563  }
1564  MoFEMFunctionReturn(0);
1565  }
1566 
1567  MoFEMErrorCode postProcess() {
1568  MoFEMFunctionBegin;
1569  switch (snes_ctx) {
1570  case CTX_SNESSETFUNCTION: {
1571  CHKERR VecAssemblyBegin(snes_f);
1572  CHKERR VecAssemblyEnd(snes_f);
1573  } break;
1574  case CTX_SNESSETJACOBIAN: {
1575  CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
1576  SCATTER_FORWARD);
1577  CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES,
1578  SCATTER_FORWARD);
1579  } break;
1580  default:
1581  break;
1582  }
1583  MoFEMFunctionReturn(0);
1584  }
1585 
1586  /** \brief Calculate db
1587  */
1588  MoFEMErrorCode calculateDb() {
1589  MoFEMFunctionBegin;
1590 
1591  MOFEM_LOG_CHANNEL("WORLD");
1592  MOFEM_LOG_FUNCTION();
1593  MOFEM_LOG_TAG("WORLD", "Arc");
1594 
1595  double alpha = arcPtr->alpha;
1596  double beta = arcPtr->beta;
1597  double d_lambda = arcPtr->dLambda;
1598 
1599  CHKERR VecZeroEntries(ghostIntLambda);
1600  CHKERR VecGhostUpdateBegin(ghostIntLambda, INSERT_VALUES,
1601  SCATTER_FORWARD);
1602  CHKERR VecGhostUpdateEnd(ghostIntLambda, INSERT_VALUES, SCATTER_FORWARD);
1603  CHKERR VecZeroEntries(arcPtr->db);
1604  CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1605  CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1606 
1607  // calculate internal lambda
1608  feIntLambda->getOpPtrVector().clear();
1609  feIntLambda->getOpPtrVector().push_back(new OpIntLambda(
1610  arcPtr->mField, blockData, commonData, lambdaFieldName, *this));
1611  feIntLambda->getOpPtrVector().push_back(new OpAssemble_db(
1612  tAg, arcPtr->mField, blockData, commonData, lambdaFieldName, *this));
1613  CHKERR arcPtr->mField.loop_finite_elements(
1614  problemPtr->getName(), "CRACKFRONT_AREA_ELEM", (*feIntLambda));
1615  const double *dx;
1616  const auto bit_field_number = getFieldBitNumber(lambdaFieldName);
1617  CHKERR VecGetArrayRead(arcPtr->dx, &dx);
1618  for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
1619  bit_field_number, dit)) {
1620  if (static_cast<int>(dit->get()->getPart()) !=
1621  arcPtr->mField.get_comm_rank())
1622  continue;
1623  int idx = dit->get()->getPetscLocalDofIdx();
1624  double val = dx[idx];
1625  CHKERR VecSetValue(ghostIntLambda, 1, val, ADD_VALUES);
1626  }
1627  CHKERR VecRestoreArrayRead(arcPtr->dx, &dx);
1628 
1629  CHKERR VecAssemblyBegin(ghostIntLambda);
1630  CHKERR VecAssemblyEnd(ghostIntLambda);
1631  CHKERR VecGhostUpdateBegin(ghostIntLambda, ADD_VALUES, SCATTER_REVERSE);
1632  CHKERR VecGhostUpdateEnd(ghostIntLambda, ADD_VALUES, SCATTER_REVERSE);
1633  CHKERR VecGhostUpdateBegin(ghostIntLambda, INSERT_VALUES,
1634  SCATTER_FORWARD);
1635  CHKERR VecGhostUpdateEnd(ghostIntLambda, INSERT_VALUES, SCATTER_FORWARD);
1636 
1637  CHKERR VecAssemblyBegin(arcPtr->db);
1638  CHKERR VecAssemblyEnd(arcPtr->db);
1639  CHKERR VecGhostUpdateBegin(arcPtr->db, ADD_VALUES, SCATTER_REVERSE);
1640  CHKERR VecGhostUpdateEnd(arcPtr->db, ADD_VALUES, SCATTER_REVERSE);
1641  CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1642  CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1643 
1644  intLambda = alpha * intLambdaDelta +
1645  alpha * beta * beta * d_lambda * d_lambda * arcPtr->F_lambda2;
1646 
1647  MOFEM_LOG_C(
1648  "WORLD", Sev::noisy,
1649  "\tintLambdaLambda = %9.8e intLambdaDelta = %9.8e intLambda = %9.8e",
1650  intLambdaLambda, intLambdaDelta, intLambda);
1651 
1652  double nrm;
1653  CHKERR VecNorm(arcPtr->db, NORM_2, &nrm);
1654  MOFEM_LOG_C("WORLD", Sev::noisy, "\tdb nrm = %9.8e", nrm);
1655 
1656  MoFEMFunctionReturn(0);
1657  }
1658 
1659  MoFEMErrorCode calculateDxAndDlambda(Vec x) {
1660  MoFEMFunctionBegin;
1661  MOFEM_LOG_CHANNEL("WORLD");
1662  MOFEM_LOG_FUNCTION();
1663  MOFEM_LOG_TAG("WORLD", "Arc");
1664 
1665  // Calculate dx
1666  CHKERR VecGhostUpdateBegin(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
1667  CHKERR VecGhostUpdateEnd(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
1668 
1669  Vec l_x, l_x0, l_dx;
1670  CHKERR VecGhostGetLocalForm(x, &l_x);
1671  CHKERR VecGhostGetLocalForm(arcPtr->x0, &l_x0);
1672  CHKERR VecGhostGetLocalForm(arcPtr->dx, &l_dx);
1673  {
1674  double *x_array, *x0_array, *dx_array;
1675  CHKERR VecGetArray(l_x, &x_array);
1676  CHKERR VecGetArray(l_x0, &x0_array);
1677  CHKERR VecGetArray(l_dx, &dx_array);
1678  int size =
1679  problemPtr->getNbLocalDofsRow() + problemPtr->getNbGhostDofsRow();
1680  for (int i = 0; i != size; ++i) {
1681  dx_array[i] = x_array[i] - x0_array[i];
1682  }
1683  CHKERR VecRestoreArray(l_x, &x_array);
1684  CHKERR VecRestoreArray(l_x0, &x0_array);
1685  CHKERR VecRestoreArray(l_dx, &dx_array);
1686  }
1687  CHKERR VecGhostRestoreLocalForm(x, &l_x);
1688  CHKERR VecGhostRestoreLocalForm(arcPtr->x0, &l_x0);
1689  CHKERR VecGhostRestoreLocalForm(arcPtr->dx, &l_dx);
1690 
1691  // Calculate dlambda
1692  if (arcPtr->getPetscLocalDofIdx() != -1) {
1693  double *array;
1694  CHKERR VecGetArray(arcPtr->dx, &array);
1695  arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
1696  array[arcPtr->getPetscLocalDofIdx()] = 0;
1697  CHKERR VecRestoreArray(arcPtr->dx, &array);
1698  }
1699  CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
1700  SCATTER_FORWARD);
1701  CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
1702  SCATTER_FORWARD);
1703 
1704  // Calculate dx2
1705  double x_nrm, x0_nrm;
1706  CHKERR VecNorm(x, NORM_2, &x_nrm);
1707  CHKERR VecNorm(arcPtr->x0, NORM_2, &x0_nrm);
1708  CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
1709  MOFEM_LOG_C("WORLD", Sev::noisy,
1710  "\tx norm = %6.4e x0 norm = %6.4e dx2 = %6.4e", x_nrm, x0_nrm,
1711  arcPtr->dx2);
1712  MoFEMFunctionReturn(0);
1713  }
1714  };
1715 };
1716 
1717 } // namespace FractureMechanics
1718 
1719 #endif //__GRIFFITH_FORCE_ELEMENT_HPP__
FractureMechanics::GriffithForceElement::AuxFunctions::calculateGriffithForce
MoFEMErrorCode calculateGriffithForce(const double gc, const double beta, const MatrixDouble &diff_n)
Definition: GriffithForceElement.hpp:201
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FractureMechanics::GriffithForceElement::FrontArcLengthControl::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: GriffithForceElement.hpp:1291
FractureMechanics::GriffithForceElement::BlockData
Definition: GriffithForceElement.hpp:74
FractureMechanics::GriffithForceElement::BlockData::r
double r
regularity parameter
Definition: GriffithForceElement.hpp:77
FractureMechanics::GriffithForceElement::getElementOptions
static MoFEMErrorCode getElementOptions(BlockData &block_data)
Definition: GriffithForceElement.hpp:83
FractureMechanics::GriffithForceElement::OpConstrainsRhs::OpConstrainsRhs
OpConstrainsRhs(MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data, const std::string &lambda_field_name)
Definition: GriffithForceElement.hpp:1000
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs::mwlsApprox
boost::shared_ptr< MWLSApprox > mwlsApprox
Definition: GriffithForceElement.hpp:437
FractureMechanics::GriffithForceElement::CommonData::tangentGriffithForceRowPtr
vector< double * > tangentGriffithForceRowPtr
Pointers to rows of tangent matrix.
Definition: GriffithForceElement.hpp:68
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:1644
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FractureMechanics::GriffithForceElement::AuxOp::AuxOp
AuxOp(int tag, BlockData &block_data, CommonData &common_data)
Definition: GriffithForceElement.hpp:221
FTensor::Tensor1< TYPE, 3 >
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
FractureMechanics::GriffithForceElement::OpConstrainsLhs::OpConstrainsLhs
OpConstrainsLhs(MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data, const std::string &lambda_field_name, SmartPetscObj< Vec > &delta_vec)
Definition: GriffithForceElement.hpp:1057
FractureMechanics::GriffithForceElement::OpConstrainsRhs::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:997
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:1334
FractureMechanics::GriffithForceElement::FrontArcLengthControl::~FrontArcLengthControl
virtual ~FrontArcLengthControl()
Definition: GriffithForceElement.hpp:1514
FractureMechanics::GriffithForceElement::AuxFunctions::k
FTensor::Index< 'k', 3 > k
Definition: GriffithForceElement.hpp:136
rho
double rho
Definition: plastic.cpp:140
FractureMechanics::GriffithForceElement::OpLhs::mat
MatrixDouble mat
Definition: GriffithForceElement.hpp:611
FractureMechanics::calMax
static double calMax(double a, double b, double r)
Definition: GriffithForceElement.hpp:24
FractureMechanics::GriffithForceElement::FrontArcLengthControl::lambdaFieldName
std::string lambdaFieldName
Definition: GriffithForceElement.hpp:1294
FractureMechanics::GriffithForceElement::OpConstrainsLhs::nF_dX
MatrixDouble nF_dX
Definition: GriffithForceElement.hpp:1071
FractureMechanics::GriffithForceElement::AuxOp::activeVariables
VectorDouble activeVariables
Definition: GriffithForceElement.hpp:230
FractureMechanics::GriffithForceElement::OpLhs::betaGC
const double betaGC
Definition: GriffithForceElement.hpp:612
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda::OpIntLambda
OpIntLambda(MoFEM::Interface &m_field, GriffithForceElement::BlockData &block_data, GriffithForceElement::CommonData &common_data, const std::string &lambda_field_name, FrontArcLengthControl &arc_front)
Definition: GriffithForceElement.hpp:1319
FractureMechanics::GriffithForceElement::AuxOp::setIndices
MoFEMErrorCode setIndices(DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:232
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FractureMechanics::GriffithForceElement::OpConstrainsLhs
Definition: GriffithForceElement.hpp:1050
FractureMechanics::GriffithForceElement::OpConstrainsRhs::nF
VectorDouble9 nF
Definition: GriffithForceElement.hpp:1009
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains::preProcess
MoFEMErrorCode preProcess()
Definition: GriffithForceElement.hpp:927
FractureMechanics::GriffithForceElement::AuxOp::rowLambdaIndicesLocal
VectorInt3 rowLambdaIndicesLocal
Definition: GriffithForceElement.hpp:228
FractureMechanics::GriffithForceElement::OpConstrainsLhs::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:1054
FractureMechanics::GriffithForceElement::AuxOp::tAg
int tAg
Definition: GriffithForceElement.hpp:217
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
FractureMechanics::GriffithForceElement::FrontArcLengthControl::blockData
BlockData & blockData
Definition: GriffithForceElement.hpp:1292
FractureMechanics::GriffithForceElement::AuxFunctions::i
FTensor::Index< 'i', 3 > i
Definition: GriffithForceElement.hpp:134
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
FractureMechanics::GriffithForceElement::AuxFunctions::N0
FTensor::Number< 0 > N0
Definition: GriffithForceElement.hpp:137
FractureMechanics::GriffithForceElement::OpConstrainsLhs::nF_dLambda
MatrixDouble nF_dLambda
Definition: GriffithForceElement.hpp:1070
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains::lambdaFieldName
const std::string lambdaFieldName
Definition: GriffithForceElement.hpp:917
FractureMechanics::GriffithForceElement::AuxOp::setLambdaNodes
MoFEMErrorCode setLambdaNodes(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, const std::string &lambda_field_name)
Definition: GriffithForceElement.hpp:276
FractureMechanics::GriffithForceElement::OpHeterogeneousGcLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: GriffithForceElement.hpp:732
FractureMechanics::GriffithForceElement::OpHeterogeneousGcLhs::OpHeterogeneousGcLhs
OpHeterogeneousGcLhs(int tag, BlockData &block_data, CommonData &common_data, const double beta_gc)
Definition: GriffithForceElement.hpp:722
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta::preProcess
MoFEMErrorCode preProcess()
Definition: GriffithForceElement.hpp:821
W
double W
Definition: free_surface.cpp:166
FractureMechanics::GriffithForceElement::OpConstrainsDelta::lambdaFieldName
const std::string lambdaFieldName
Definition: GriffithForceElement.hpp:859
FractureMechanics::GriffithForceElement::BlockData::frontNodes
Range frontNodes
Nodes on crack front.
Definition: GriffithForceElement.hpp:78
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::tRis
Range & tRis
Definition: GriffithForceElement.hpp:537
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta::ownDeltaVec
bool ownDeltaVec
Definition: GriffithForceElement.hpp:814
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db::activeVariables
VectorDouble activeVariables
Definition: GriffithForceElement.hpp:1409
FractureMechanics::GriffithForceElement::AuxFunctions::get_ftensor_from_mat_2d
boost::function< FTensor::Tensor1< FTensor::PackPtr< const double *, 2 >, 2 > const MatrixDouble &m)> get_ftensor_from_mat_2d
Definition: GriffithForceElement.hpp:145
FractureMechanics::GriffithForceElement::BlockData::gc
double gc
Griffith energy.
Definition: GriffithForceElement.hpp:75
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta::deltaVec
SmartPetscObj< Vec > deltaVec
Definition: GriffithForceElement.hpp:815
FractureMechanics::GriffithForceElement::FrontArcLengthControl::postProcess
MoFEMErrorCode postProcess()
Definition: GriffithForceElement.hpp:1567
FractureMechanics::GriffithForceElement::FrontArcLengthControl
arc-length element
Definition: GriffithForceElement.hpp:1288
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db::gradDelta
VectorDouble gradDelta
Definition: GriffithForceElement.hpp:1408
MoFEM::Types::VectorInt3
VectorBoundedArray< int, 3 > VectorInt3
Definition: Types.hpp:87
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db::OpAssemble_db
OpAssemble_db(int tag, MoFEM::Interface &m_field, GriffithForceElement::BlockData &block_data, GriffithForceElement::CommonData &common_data, const std::string &lambda_field_name, FrontArcLengthControl &arc_front)
Definition: GriffithForceElement.hpp:1398
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs
Calculate density at the crack fron nodes and multiplity gc.
Definition: GriffithForceElement.hpp:432
FractureMechanics::GriffithForceElement::AuxFunctions::currentCoords
ublas::vector< TYPE > currentCoords
Definition: GriffithForceElement.hpp:126
FractureMechanics::GriffithForceElement::feLhs
MyTriangleFE & feLhs
Definition: GriffithForceElement.hpp:56
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::auxFun
AuxFunctions< double > auxFun
Definition: GriffithForceElement.hpp:544
FractureMechanics::GriffithForceElement::OpConstrainsLhs::d_Delta
VectorDouble d_Delta
Definition: GriffithForceElement.hpp:1076
FractureMechanics::GriffithForceElement::OpRhs::OpRhs
OpRhs(int tag, BlockData &block_data, CommonData &common_data)
Definition: GriffithForceElement.hpp:332
MOFEM_LOG_FUNCTION
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:325
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
FractureMechanics::GriffithForceElement::OpConstrainsDelta::deltaVec
SmartPetscObj< Vec > & deltaVec
Definition: GriffithForceElement.hpp:860
FractureMechanics::GriffithForceElement::OpConstrainsDelta::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:873
FractureMechanics::GriffithForceElement::AuxOp::setVariables
MoFEMErrorCode setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:261
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
FractureMechanics::GriffithForceElement::BlockData::BlockData
BlockData()
Definition: GriffithForceElement.hpp:79
FractureMechanics::GriffithForceElement::OpHeterogeneousGcLhs
Definition: GriffithForceElement.hpp:718
FractureMechanics::GriffithForceElement::FrontArcLengthControl::calculateDb
MoFEMErrorCode calculateDb()
Calculate db.
Definition: GriffithForceElement.hpp:1588
FractureMechanics::GriffithForceElement::FrontArcLengthControl::tAg
int tAg
Definition: GriffithForceElement.hpp:1290
FEMethod
FractureMechanics::GriffithForceElement::GriffithForceElement
GriffithForceElement(MoFEM::Interface &m_field)
Definition: GriffithForceElement.hpp:59
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FractureMechanics::GriffithForceElement::AuxFunctions::referenceCoords
ublas::vector< TYPE > referenceCoords
Definition: GriffithForceElement.hpp:125
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db
assemble internal residual derivative
Definition: GriffithForceElement.hpp:1396
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
FractureMechanics::GriffithForceElement::OpCalculateGriffithForce::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:387
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:1411
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta::MyTriangleFEConstrainsDelta
MyTriangleFEConstrainsDelta(MoFEM::Interface &m_field, const std::string &lambda_field_name)
Definition: GriffithForceElement.hpp:817
_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_
#define _IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(PROBLEMPTR, FIELD_BIT_NUMBER, IT)
Definition: ProblemsMultiIndices.hpp:338
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::operator()
MoFEMErrorCode operator()()
Definition: GriffithForceElement.hpp:546
FractureMechanics::GriffithForceElement::MyTriangleFE::MyTriangleFE
MyTriangleFE(MoFEM::Interface &m_field)
Definition: GriffithForceElement.hpp:44
FractureMechanics::GriffithForceElement::AuxFunctions::N1
FTensor::Number< 1 > N1
Definition: GriffithForceElement.hpp:138
a
constexpr double a
Definition: approx_sphere.cpp:30
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:534
FractureMechanics::GriffithForceElement::AuxFunctions::j
FTensor::Index< 'j', 3 > j
Definition: GriffithForceElement.hpp:135
FractureMechanics::GriffithForceElement::OpConstrainsDelta::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:858
FractureMechanics::GriffithForceElement::OpConstrainsDelta::auxFun
AuxFunctions< double > auxFun
Definition: GriffithForceElement.hpp:871
convert.type
type
Definition: convert.py:64
FractureMechanics::GriffithForceElement::AuxFunctions::t_Tangent1_ksi
FTensor::Tensor1< TYPE, 3 > t_Tangent1_ksi
Definition: GriffithForceElement.hpp:130
FractureMechanics::GriffithForceElement::commonData
CommonData commonData
Common data sgared between operators.
Definition: GriffithForceElement.hpp:72
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
FractureMechanics::GriffithForceElement::FrontArcLengthControl::calculateDxAndDlambda
MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition: GriffithForceElement.hpp:1659
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains
Definition: GriffithForceElement.hpp:915
FractureMechanics::GriffithForceElement::FrontArcLengthControl::intLambda
double intLambda
Definition: GriffithForceElement.hpp:1300
FractureMechanics::GriffithForceElement::FrontArcLengthControl::feIntLambda
boost::shared_ptr< MyTriangleFE > feIntLambda
Definition: GriffithForceElement.hpp:1302
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::SaveGriffithForceOnTags
SaveGriffithForceOnTags(MoFEM::Interface &m_field, Tag th, Range &tris, Tag *th_coords=nullptr)
Definition: GriffithForceElement.hpp:540
FractureMechanics::GriffithForceElement::OpConstrainsRhs
Definition: GriffithForceElement.hpp:993
FractureMechanics::GriffithForceElement::AuxFunctions::get_ftensor_from_vec
boost::function< FTensor::Tensor1< TYPE *, 3 >ublas::vector< TYPE > &v)> get_ftensor_from_vec
Definition: GriffithForceElement.hpp:141
FractureMechanics::GriffithForceElement::AuxOp::setLambdaIndices
MoFEMErrorCode setLambdaIndices(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, const std::string &lambda_field_name)
Definition: GriffithForceElement.hpp:301
FractureMechanics::GriffithForceElement::AuxFunctions::t_Normal
FTensor::Tensor1< TYPE, 3 > t_Normal
Definition: GriffithForceElement.hpp:129
FractureMechanics::GriffithForceElement::FrontArcLengthControl::ghostIntLambda
Vec ghostIntLambda
Definition: GriffithForceElement.hpp:1296
FractureMechanics::GriffithForceElement::MyTriangleFE
Definition: GriffithForceElement.hpp:39
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs
OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs(string rho_tag_name, const double beta_gc, boost::shared_ptr< MWLSApprox > mwls_approx, MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data)
Definition: GriffithForceElement.hpp:440
FractureMechanics::GriffithForceElement::OpConstrainsLhs::a_dA
ublas::vector< adouble > a_dA
Definition: GriffithForceElement.hpp:1080
FractureMechanics::GriffithForceElement::OpCalculateGriffithForce::auxFun
AuxFunctions< double > auxFun
Definition: GriffithForceElement.hpp:386
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta::lambdaFieldName
const std::string lambdaFieldName
Definition: GriffithForceElement.hpp:813
FractureMechanics::GriffithForceElement::BlockData::E
double E
Young modulus.
Definition: GriffithForceElement.hpp:76
FractureMechanics::GriffithForceElement::OpConstrainsLhs::deltaVec
SmartPetscObj< Vec > & deltaVec
Definition: GriffithForceElement.hpp:1055
FractureMechanics::GriffithForceElement::AuxFunctions::t_Tangent2_eta
FTensor::Tensor1< TYPE, 3 > t_Tangent2_eta
Definition: GriffithForceElement.hpp:131
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda::auxFun
AuxFunctions< double > auxFun
Definition: GriffithForceElement.hpp:1332
FractureMechanics::GriffithForceElement::AuxOp::rowLambdaIndices
VectorInt3 rowLambdaIndices
Definition: GriffithForceElement.hpp:227
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains::MyTriangleFEConstrains
MyTriangleFEConstrains(MoFEM::Interface &m_field, const std::string &lambda_field_name, BlockData &block_data, SmartPetscObj< Vec > &delta_vec)
Definition: GriffithForceElement.hpp:921
FractureMechanics::GriffithForceElement::MyTriangleFE::getRule
int getRule(int order)
Definition: GriffithForceElement.hpp:48
FractureMechanics::GriffithForceElement::AuxFunctions
Definition: GriffithForceElement.hpp:123
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta
Definition: GriffithForceElement.hpp:811
FractureMechanics::GriffithForceElement::getLoopFeRhs
MyTriangleFE & getLoopFeRhs()
Definition: GriffithForceElement.hpp:55
FractureMechanics::GriffithForceElement::AuxOp::blockData
BlockData & blockData
Definition: GriffithForceElement.hpp:218
FractureMechanics::GriffithForceElement::OpConstrainsRhs::lambdaFieldName
const std::string lambdaFieldName
Definition: GriffithForceElement.hpp:998
FractureMechanics::GriffithForceElement::FrontArcLengthControl::preProcess
MoFEMErrorCode preProcess()
Definition: GriffithForceElement.hpp:1519
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FractureMechanics::GriffithForceElement::AuxFunctions::currentArea
TYPE currentArea
Definition: GriffithForceElement.hpp:132
FractureMechanics::GriffithForceElement::OpLhs::OpLhs
OpLhs(int tag, BlockData &block_data, CommonData &common_data, const double beta_gc=0)
Definition: GriffithForceElement.hpp:604
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FaceElementForcesAndSourcesCore
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:1315
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs::betaGC
const double betaGC
Definition: GriffithForceElement.hpp:435
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
FractureMechanics::GriffithForceElement::OpLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: GriffithForceElement.hpp:614
FractureMechanics::GriffithForceElement::OpHeterogeneousGcLhs::mat
MatrixDouble mat
Definition: GriffithForceElement.hpp:729
FractureMechanics::GriffithForceElement::OpHeterogeneousGcLhs::betaGC
const double betaGC
Definition: GriffithForceElement.hpp:730
convert.n
n
Definition: convert.py:82
FractureMechanics::GriffithForceElement::OpRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:336
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags
Definition: GriffithForceElement.hpp:532
FractureMechanics::GriffithForceElement::OpConstrainsDelta
Definition: GriffithForceElement.hpp:854
FractureMechanics::GriffithForceElement::OpConstrainsLhs::nG_dX
MatrixDouble nG_dX
Definition: GriffithForceElement.hpp:1069
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db::a_auxFun
AuxFunctions< adouble > a_auxFun
Definition: GriffithForceElement.hpp:1407
FractureMechanics::GriffithForceElement::CommonData::diffRho
MatrixDouble diffRho
for lhs with heterogeneous gc
Definition: GriffithForceElement.hpp:70
FractureMechanics::GriffithForceElement::OpConstrainsLhs::a_auxFun
AuxFunctions< adouble > a_auxFun
Definition: GriffithForceElement.hpp:1079
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
adouble
FractureMechanics::GriffithForceElement::OpCalculateGriffithForce::OpCalculateGriffithForce
OpCalculateGriffithForce(int tag, BlockData &block_data, CommonData &common_data)
Definition: GriffithForceElement.hpp:381
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
FractureMechanics::GriffithForceElement::OpConstrainsLhs::a_Delta
ublas::vector< adouble > a_Delta
Definition: GriffithForceElement.hpp:1078
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains::postProcess
MoFEMErrorCode postProcess()
Definition: GriffithForceElement.hpp:987
FractureMechanics::GriffithForceElement::AuxOp::lambdaAtNodes
VectorDouble3 lambdaAtNodes
Definition: GriffithForceElement.hpp:229
FractureMechanics::GriffithForceElement::AuxFunctions::griffithForce
ublas::vector< TYPE > griffithForce
Definition: GriffithForceElement.hpp:127
FractureMechanics::GriffithForceElement::OpCalculateGriffithForce
calculate griffith force vector
Definition: GriffithForceElement.hpp:378
FractureMechanics::GriffithForceElement::AuxOp::rowIndicesLocal
ublas::vector< int > rowIndicesLocal
Definition: GriffithForceElement.hpp:225
FractureMechanics::GriffithForceElement::AuxFunctions::calculateAreaAndNormal
MoFEMErrorCode calculateAreaAndNormal(const MatrixDouble &diff_n)
Definition: GriffithForceElement.hpp:159
FractureMechanics::GriffithForceElement::MyTriangleFEConstrainsDelta::postProcess
MoFEMErrorCode postProcess()
Definition: GriffithForceElement.hpp:835
FractureMechanics::GriffithForceElement::CommonData::densityRho
VectorDouble densityRho
gc * rho^beta
Definition: GriffithForceElement.hpp:69
FractureMechanics::GriffithForceElement::feLhsPtr
boost::shared_ptr< MyTriangleFE > feLhsPtr
Definition: GriffithForceElement.hpp:52
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FractureMechanics::GriffithForceElement::feRhsPtr
boost::shared_ptr< MyTriangleFE > feRhsPtr
Definition: GriffithForceElement.hpp:51
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda::arcFront
FrontArcLengthControl & arcFront
Definition: GriffithForceElement.hpp:1317
FractureMechanics::GriffithForceElement::getLoopFeLhs
MyTriangleFE & getLoopFeLhs()
Definition: GriffithForceElement.hpp:57
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
FractureMechanics::GriffithForceElement
Implementation of Griffith element.
Definition: GriffithForceElement.hpp:35
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
FractureMechanics::GriffithForceElement::FrontArcLengthControl::commonData
CommonData & commonData
Definition: GriffithForceElement.hpp:1293
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FractureMechanics::GriffithForceElement::blockData
map< int, BlockData > blockData
Shared data between finite elements.
Definition: GriffithForceElement.hpp:81
FractureMechanics::diffCalMax_a
static double diffCalMax_a(double a, double b, double r)
Definition: GriffithForceElement.hpp:28
FractureMechanics::GriffithForceElement::FrontArcLengthControl::FrontArcLengthControl
FrontArcLengthControl(int tag, GriffithForceElement::BlockData &block_data, GriffithForceElement::CommonData &common_data, const std::string &lambda_field_name, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
Definition: GriffithForceElement.hpp:1497
FractureMechanics::GriffithForceElement::OpConstrainsLhs::jacDelta
MatrixDouble jacDelta
Definition: GriffithForceElement.hpp:1073
FractureMechanics::GriffithForceElement::AuxOp::rowIndices
ublas::vector< int > rowIndices
Definition: GriffithForceElement.hpp:222
FractureMechanics::GriffithForceElement::OpConstrainsLhs::jacDeltaPtr
ublas::vector< double * > jacDeltaPtr
Definition: GriffithForceElement.hpp:1074
FractureMechanics::GriffithForceElement::FrontArcLengthControl::operator()
MoFEMErrorCode operator()()
Definition: GriffithForceElement.hpp:1542
FractureMechanics::GriffithForceElement::OpConstrainsLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: GriffithForceElement.hpp:1083
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains::deltaVec
SmartPetscObj< Vec > deltaVec
Definition: GriffithForceElement.hpp:919
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpAssemble_db::tAg
int tAg
Definition: GriffithForceElement.hpp:1397
FractureMechanics::GriffithForceElement::FrontArcLengthControl::intLambdaLambda
double & intLambdaLambda
Definition: GriffithForceElement.hpp:1299
FractureMechanics::GriffithForceElement::feRhs
MyTriangleFE & feRhs
Definition: GriffithForceElement.hpp:54
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
FractureMechanics::GriffithForceElement::FrontArcLengthControl::intLambdaDelta
double & intLambdaDelta
Definition: GriffithForceElement.hpp:1298
FractureMechanics::GriffithForceElement::FrontArcLengthControl::intLambdaArray
double intLambdaArray[2]
Definition: GriffithForceElement.hpp:1297
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FractureMechanics::GriffithForceElement::OpConstrainsLhs::a_lambdaAtNodes
ublas::vector< adouble > a_lambdaAtNodes
Definition: GriffithForceElement.hpp:1081
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
FractureMechanics::GriffithForceElement::OpConstrainsRhs::auxFun
AuxFunctions< double > auxFun
Definition: GriffithForceElement.hpp:1008
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda
assemble internal residual
Definition: GriffithForceElement.hpp:1312
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &row_data)
Definition: GriffithForceElement.hpp:449
FractureMechanics::GriffithForceElement::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:37
FractureMechanics::GriffithForceElement::OpLhs
Definition: GriffithForceElement.hpp:601
MoFEM::Types::VectorDouble9
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:96
FractureMechanics::GriffithForceElement::CommonData
Definition: GriffithForceElement.hpp:64
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
FractureMechanics::GriffithForceElement::MyTriangleFE::B
Mat B
Definition: GriffithForceElement.hpp:41
FractureMechanics::GriffithForceElement::OpConstrainsLhs::d_dA
VectorDouble d_dA
Definition: GriffithForceElement.hpp:1075
FractureMechanics::GriffithForceElement::MyTriangleFEConstrains::blockData
BlockData & blockData
Definition: GriffithForceElement.hpp:918
FractureMechanics::GriffithForceElement::AuxFunctions::calculateMatA
MoFEMErrorCode calculateMatA(const MatrixDouble &diff_n)
Definition: GriffithForceElement.hpp:180
FractureMechanics::GriffithForceElement::CommonData::griffithForce
VectorDouble griffithForce
Griffith force at nodes.
Definition: GriffithForceElement.hpp:65
FractureMechanics::GriffithForceElement::OpConstrainsDelta::OpConstrainsDelta
OpConstrainsDelta(MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data, const std::string &lambda_field_name, SmartPetscObj< Vec > &delta_vec)
Definition: GriffithForceElement.hpp:862
FractureMechanics::GriffithForceElement::CommonData::tangentGriffithForce
MatrixDouble tangentGriffithForce
Tangent matrix.
Definition: GriffithForceElement.hpp:66
FractureMechanics::GriffithForceElement::AuxFunctions::AuxFunctions
AuxFunctions()
Definition: GriffithForceElement.hpp:147
FractureMechanics::GriffithForceElement::OpConstrainsRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: GriffithForceElement.hpp:1012
FractureMechanics::GriffithForceElement::AuxOp::commonData
CommonData & commonData
Definition: GriffithForceElement.hpp:219
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
FractureMechanics
Definition: AnalyticalFun.hpp:15
FractureMechanics::GriffithForceElement::OpRhs
Definition: GriffithForceElement.hpp:330
FractureMechanics::GriffithForceElement::AuxOp
Definition: GriffithForceElement.hpp:215
FractureMechanics::GriffithForceElement::FrontArcLengthControl::OpIntLambda::lambdaFieldName
std::string lambdaFieldName
Definition: GriffithForceElement.hpp:1316
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::thCoordsPtr
Tag * thCoordsPtr
Definition: GriffithForceElement.hpp:538
FractureMechanics::GriffithForceElement::SaveGriffithForceOnTags::tH
Tag tH
Definition: GriffithForceElement.hpp:536
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs::mField
MoFEM::Interface & mField
Definition: GriffithForceElement.hpp:438
FractureMechanics::GriffithForceElement::FrontArcLengthControl::calculateLambdaInt
double calculateLambdaInt()
Calculate internal lambda.
Definition: GriffithForceElement.hpp:1536
FractureMechanics::GriffithForceElement::MyTriangleFE::F
Vec F
Definition: GriffithForceElement.hpp:42
FractureMechanics::GriffithForceElement::OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs::rhoTagName
std::string rhoTagName
Definition: GriffithForceElement.hpp:436
Generated by Doxygen 1.8.17 and hosted at University of Glasgow