v0.14.0
ConstantArea.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 __CONSTANT_AREA_HPP__
16 #define __CONSTANT_AREA_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 /** \brief Constant area constrains
25  */
26 struct ConstantArea {
27 
28  MoFEM::Interface &mField;
29  ConstantArea(MoFEM::Interface &m_field)
30  : mField(m_field), commonData(m_field) {}
31 
32  struct CommonData {
33 
34  MoFEM::Interface &mField;
35  CommonData(MoFEM::Interface &m_field) : mField(m_field) {}
36 
37  VectorDouble crackAreaIncrement;
38  MatrixDouble tangentCrackAreaIncrement;
39  vector<double *> tangentCrackAreaIncrementRowPtr;
40 
41  VectorDouble crackFrontTangent;
42  MatrixDouble tangentCrackFrontTangent;
43  vector<double *> tangentCrackFrontTangentRowPtr;
44 
45  map<int, Vec> mapV;
46  bool setMapV;
47 
48  Tag thInterfaceSide;
49  };
50  CommonData commonData;
51 
52  struct MyTriangleFE : public FaceElementForcesAndSourcesCore {
53 
54  CommonData &commonData;
55  std::string fieldName;
56 
57  Mat petscB;
58  Vec petscF;
59  Mat petscQ;
60 
61  MyTriangleFE(MoFEM::Interface &m_field, CommonData &common_data,
62  std::string field_name = "LAMBDA_CRACKFRONT_AREA")
63  : FaceElementForcesAndSourcesCore(m_field), commonData(common_data),
64  fieldName(field_name), petscB(PETSC_NULL), petscF(PETSC_NULL),
65  petscQ(PETSC_NULL) {}
66  int getRule(int order) { return order; };
67 
68  PetscErrorCode preProcess() {
69  MoFEMFunctionBegin;
70  CHKERR FaceElementForcesAndSourcesCore::preProcess();
71 
72  if (petscB != PETSC_NULL) {
73  snes_B = petscB;
74  }
75 
76  if (petscF != PETSC_NULL) {
77  snes_f = petscF;
78  }
79 
80  switch (ts_ctx) {
81  case CTX_TSSETIFUNCTION: {
82  if (!petscF) {
83  snes_ctx = CTX_SNESSETFUNCTION;
84  snes_f = ts_F;
85  }
86  break;
87  }
88  case CTX_TSSETIJACOBIAN: {
89  if (!petscB) {
90  snes_ctx = CTX_SNESSETJACOBIAN;
91  snes_B = ts_B;
92  }
93  break;
94  }
95  default:
96  break;
97  }
98 
99  // Tag to detect which side is it
100  CHKERR mField.get_moab().tag_get_handle("INTERFACE_SIDE",
101  commonData.thInterfaceSide);
102 
103  // If petscQ matrix is defined appply projection to each row
104  if (petscQ != PETSC_NULL) {
105  if (petscB == PETSC_NULL) {
106  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
107  "Matrix B should be set");
108  }
109  const auto field_bit_number = mField.get_field_bit_number(fieldName);
110  for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr, field_bit_number,
111  dit)) {
112  CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
113  problemPtr->getName(), COL,
114  &commonData.mapV[dit->get()->getPetscGlobalDofIdx()]);
115  CHKERR VecZeroEntries(
116  commonData.mapV[dit->get()->getPetscGlobalDofIdx()]);
117  CHKERR VecGhostUpdateBegin(
118  commonData.mapV[dit->get()->getPetscGlobalDofIdx()],
119  INSERT_VALUES, SCATTER_FORWARD);
120  CHKERR VecGhostUpdateEnd(
121  commonData.mapV[dit->get()->getPetscGlobalDofIdx()],
122  INSERT_VALUES, SCATTER_FORWARD);
123  }
124  commonData.setMapV = true;
125  } else {
126  commonData.setMapV = false;
127  }
128 
129  MoFEMFunctionReturn(0);
130  }
131 
132  PetscErrorCode postProcess() {
133  MoFEMFunctionBegin;
134  CHKERR FaceElementForcesAndSourcesCore::postProcess();
135  if (petscQ != PETSC_NULL) {
136  if (petscB == PETSC_NULL) {
137  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
138  "Matrix B should be set");
139  }
140 
141  // Q vectors
142  void *void_ctx;
143  CHKERR MatShellGetContext(petscQ, &void_ctx);
144  ConstrainMatrixCtx *q_mat_ctx = (ConstrainMatrixCtx *)void_ctx;
145 
146  Vec Qv;
147  CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
148  problemPtr->getName(), COL, &Qv);
149  CHKERR MatAssemblyBegin(petscB, MAT_FLUSH_ASSEMBLY);
150  CHKERR MatAssemblyEnd(petscB, MAT_FLUSH_ASSEMBLY);
151  auto project_vectors_and_assemble = [this, Qv](auto dit) {
152  MoFEMFunctionBegin;
153  if (commonData.mapV.find(dit->get()->getPetscGlobalDofIdx()) ==
154  commonData.mapV.end()) {
155  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
156  "Vector on row not found");
157  }
158  CHKERR VecAssemblyBegin(
159  commonData.mapV[dit->get()->getPetscGlobalDofIdx()]);
160  CHKERR VecAssemblyEnd(
161  commonData.mapV[dit->get()->getPetscGlobalDofIdx()]);
162  CHKERR VecGhostUpdateBegin(
163  commonData.mapV[dit->get()->getPetscGlobalDofIdx()], ADD_VALUES,
164  SCATTER_REVERSE);
165  CHKERR VecGhostUpdateEnd(
166  commonData.mapV[dit->get()->getPetscGlobalDofIdx()], ADD_VALUES,
167  SCATTER_REVERSE);
168  CHKERR VecGhostUpdateBegin(
169  commonData.mapV[dit->get()->getPetscGlobalDofIdx()],
170  INSERT_VALUES, SCATTER_FORWARD);
171  CHKERR VecGhostUpdateEnd(
172  commonData.mapV[dit->get()->getPetscGlobalDofIdx()],
173  INSERT_VALUES, SCATTER_FORWARD);
174 
175  CHKERR MatMult(
176  petscQ, commonData.mapV[dit->get()->getPetscGlobalDofIdx()], Qv);
177  CHKERR VecGhostUpdateBegin(Qv, INSERT_VALUES, SCATTER_FORWARD);
178  CHKERR VecGhostUpdateEnd(Qv, INSERT_VALUES, SCATTER_FORWARD);
179 
180  if (0) {
181  auto save_vec = [](std::string name, Vec v) {
182  PetscViewer viewer;
183  PetscViewerASCIIOpen(PETSC_COMM_WORLD, name.c_str(), &viewer);
184  PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
185  VecView(v, viewer);
186  PetscViewerPopFormat(viewer);
187  PetscViewerDestroy(&viewer);
188  };
189  save_vec(
190  "V_" + boost::lexical_cast<std::string>(dit->get()->getEnt()) +
191  ".m",
192  commonData.mapV[dit->get()->getPetscGlobalDofIdx()]);
193  save_vec(
194  "Qv_" + boost::lexical_cast<std::string>(dit->get()->getEnt()) +
195  ".m",
196  Qv);
197  }
198 
199  if (dit->get()->getPart() == mField.get_comm_rank()) {
200  double *array;
201  CHKERR VecGetArray(Qv, &array);
202  vector<int> glob_idx;
203  vector<double> vals;
204  int row = dit->get()->getPetscGlobalDofIdx();
205  for (_IT_NUMEREDDOF_COL_BY_BITNUMBER_FOR_LOOP_(
206  problemPtr,
207  mField.get_field_bit_number("MESH_NODE_POSITIONS"),
208  diit)) {
209  int idx = diit->get()->getPetscGlobalDofIdx();
210  double val = array[diit->get()->getPetscGlobalDofIdx()];
211  glob_idx.push_back(idx);
212  vals.push_back(val);
213  }
214  CHKERR MatSetValues(petscB, 1, &row, glob_idx.size(),
215  &*glob_idx.begin(), &*vals.begin(), ADD_VALUES);
216  CHKERR VecRestoreArray(Qv, &array);
217  }
218 
219  MoFEMFunctionReturn(0);
220  };
221 
222  for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(
223  problemPtr, mField.get_field_bit_number(fieldName), dit)) {
224  CHKERR project_vectors_and_assemble(dit);
225  CHKERR VecDestroy(
226  &commonData.mapV[dit->get()->getPetscGlobalDofIdx()]);
227  }
228 
229  CHKERR VecDestroy(&Qv);
230  commonData.mapV.clear();
231  }
232 
233  MoFEMFunctionReturn(0);
234  }
235  };
236 
237  boost::shared_ptr<MyTriangleFE> feRhsPtr;
238  boost::shared_ptr<MyTriangleFE> feLhsPtr;
239 
240  struct OpAreaJacobian
241  : public FaceElementForcesAndSourcesCore::UserDataOperator {
242 
243  int tAg;
244  CommonData &commonData;
245  bool isTapeRecorded;
246 
247  OpAreaJacobian(int tag, CommonData &common_data)
248  : FaceElementForcesAndSourcesCore::UserDataOperator(
249  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
250  tAg(tag), commonData(common_data), isTapeRecorded(false) {
251  sYmm = false;
252  }
253 
254  /**
255 
256  \f[
257  A = \| \mathbf{N} \| =
258  \left\|
259  \textrm{Spin}\left[
260  \frac{\partial\mathbf{N}\mathbf{X}}{\partial\xi}
261  \right]
262  \frac{\partial\mathbf{N}\mathbf{X}}{\partial\eta}
263  \right\|
264  \f]
265 
266  */
267  template <typename TYPE> struct AuxFunctions {
268 
269  MatrixDouble N, NTN;
270  MatrixDouble Bksi;
271  MatrixDouble Beta;
272 
273  ublas::vector<TYPE> referenceCoords;
274 
275  ublas::vector<TYPE> referenceXdKsi;
276  ublas::vector<TYPE> referenceXdEta;
277  ublas::matrix<TYPE> referenceSpinKsi;
278  ublas::matrix<TYPE> referenceSpinEta;
279  ublas::vector<TYPE> referenceNormal;
280  TYPE referenceArea;
281 
282  ublas::vector<TYPE> currentCoords;
283  ublas::vector<TYPE> currentXdKsi;
284  ublas::vector<TYPE> currentXdEta;
285  ublas::matrix<TYPE> currentSpinKsi;
286  ublas::matrix<TYPE> currentSpinEta;
287  ublas::vector<TYPE> currentNormal;
288  TYPE currentArea;
289 
290  ublas::matrix<TYPE> A;
291  ublas::vector<TYPE> crackAreaIncrement;
292  ublas::vector<TYPE> crackFrontTangent;
293 
294  PetscErrorCode sPin(ublas::matrix<TYPE> &spin, ublas::vector<TYPE> &vec) {
295  MoFEMFunctionBeginHot;
296  spin.resize(3, 3, false);
297  std::fill(spin.data().begin(), spin.data().end(), 0);
298  spin(0, 1) = -vec[2];
299  spin(0, 2) = +vec[1];
300  spin(1, 0) = +vec[2];
301  spin(1, 2) = -vec[0];
302  spin(2, 0) = -vec[1];
303  spin(2, 1) = +vec[0];
304  MoFEMFunctionReturnHot(0);
305  }
306 
307  PetscErrorCode matrixB(int gg, DataForcesAndSourcesCore::EntData &data) {
308  MoFEMFunctionBeginHot;
309  Bksi.resize(3, 9, false);
310  Bksi.clear();
311  for (int ii = 0; ii < 3; ++ii) {
312  for (int jj = 0; jj < 3; ++jj) {
313  Bksi(jj, ii * 3 + jj) = data.getDiffN(gg)(ii, 0);
314  }
315  }
316  Beta.resize(3, 9, false);
317  Beta.clear();
318  for (int ii = 0; ii < 3; ++ii) {
319  for (int jj = 0; jj < 3; ++jj) {
320  Beta(jj, ii * 3 + jj) = data.getDiffN(gg)(ii, 1);
321  }
322  }
323  MoFEMFunctionReturnHot(0);
324  }
325 
326  PetscErrorCode dIffX() {
327  MoFEMFunctionBeginHot;
328  currentXdKsi.resize(3, false);
329  currentXdEta.resize(3, false);
330  std::fill(currentXdKsi.begin(), currentXdKsi.end(), 0);
331  std::fill(currentXdEta.begin(), currentXdEta.end(), 0);
332  for (auto ii = 0; ii != Bksi.size1(); ++ii)
333  for (auto jj = 0; jj != Bksi.size2(); ++jj)
334  currentXdKsi(ii) += Bksi(ii, jj) * currentCoords(jj);
335  for (auto ii = 0; ii != Beta.size1(); ++ii)
336  for (auto jj = 0; jj != Beta.size2(); ++jj)
337  currentXdEta(ii) += Beta(ii, jj) * currentCoords(jj);
338  // noalias(currentXdKsi) = prod(Bksi, currentCoords);
339  // noalias(currentXdEta) = prod(Beta, currentCoords);
340  MoFEMFunctionReturnHot(0);
341  }
342 
343  PetscErrorCode nOrmal() {
344  MoFEMFunctionBeginHot;
345 
346  CHKERR sPin(currentSpinKsi, currentXdKsi);
347  CHKERR sPin(currentSpinEta, currentXdEta);
348  currentNormal.resize(3, false);
349  std::fill(currentNormal.begin(), currentNormal.end(), 0);
350  for (auto ii = 0; ii != currentSpinKsi.size1(); ++ii)
351  for (auto jj = 0; jj != currentSpinKsi.size2(); ++jj)
352  currentNormal(ii) +=
353  0.5 * currentSpinKsi(ii, jj) * currentXdEta(jj);
354  // noalias(currentNormal) = 0.5 * prod(currentSpinKsi, currentXdEta);
355  currentArea = 0;
356  for (int dd = 0; dd != 3; ++dd) {
357  currentArea += currentNormal[dd] * currentNormal[dd];
358  }
359  currentArea = sqrt(currentArea);
360  MoFEMFunctionReturnHot(0);
361  }
362 
363  PetscErrorCode matrixA() {
364  MoFEMFunctionBeginHot;
365  A.resize(3, 9, false);
366  A.clear();
367  for (auto ii = 0; ii != A.size1(); ++ii)
368  for (auto jj = 0; jj != A.size2(); ++jj)
369  for (auto kk = 0; kk != currentSpinKsi.size2(); ++kk)
370  A(ii, jj) += 0.5 * (currentSpinKsi(ii, kk) * Beta(kk, jj) -
371  currentSpinEta(ii, kk) * Bksi(kk, jj));
372  // noalias(A) =
373  // 0.5 * (prod(currentSpinKsi, Beta) - prod(currentSpinEta, Bksi));
374  MoFEMFunctionReturnHot(0);
375  }
376 
377  PetscErrorCode calculateAreaIncrementDirection(double beta) {
378  MoFEMFunctionBeginHot;
379  for (int dd = 0; dd != 9; dd++) {
380  for (int ii = 0; ii != 3; ii++) {
381  crackAreaIncrement[dd] +=
382  beta * A(ii, dd) * currentNormal[ii] / currentArea;
383  }
384  }
385  MoFEMFunctionReturnHot(0);
386  }
387 
388  // FIXME Need to implement this with antisymmetric version need fix in
389  // FTensor
390 
391  FTensor::Tensor2<TYPE, 3, 3> tSpin;
392 
393  PetscErrorCode calculateFrontTangent(double beta) {
394  MoFEMFunctionBeginHot;
395  FTensor::Index<'i', 3> i;
396  FTensor::Index<'j', 3> j;
397  tSpin(i, j) = 0;
398  tSpin(0, 1) = -currentNormal[2];
399  tSpin(0, 2) = +currentNormal[1];
400  tSpin(1, 0) = +currentNormal[2];
401  tSpin(1, 2) = -currentNormal[0];
402  tSpin(2, 0) = -currentNormal[1];
403  tSpin(2, 1) = +currentNormal[0];
404  tSpin(i, j) /= beta * currentArea;
405  crackFrontTangent.resize(9, false);
406  FTensor::Tensor1<TYPE *, 3> t_crack_front_tangent(
407  &crackFrontTangent[0], &crackFrontTangent[1], &crackFrontTangent[2],
408  3);
409  FTensor::Tensor1<TYPE *, 3> t_area_increment(&crackAreaIncrement[0],
410  &crackAreaIncrement[1],
411  &crackAreaIncrement[2], 3);
412  for (int ii = 0; ii != 3; ++ii) {
413  t_crack_front_tangent(i) = tSpin(i, j) * t_area_increment(j);
414  ++t_area_increment;
415  ++t_crack_front_tangent;
416  }
417  MoFEMFunctionReturnHot(0);
418  }
419  };
420 
421  AuxFunctions<adouble> a_auxFun;
422 
423  PetscErrorCode doWork(int side, EntityType type,
424  DataForcesAndSourcesCore::EntData &data) {
425  MoFEMFunctionBegin;
426 
427  if (isTapeRecorded) {
428  MoFEMFunctionReturnHot(0);
429  }
430  if (type != MBVERTEX) {
431  MoFEMFunctionReturnHot(0);
432  }
433 
434  int nb_dofs = data.getFieldData().size();
435  int nb_gauss_pts = data.getN().size1();
436 
437  a_auxFun.referenceCoords.resize(9, false);
438  a_auxFun.currentCoords.resize(9, false);
439 
440  trace_on(tAg);
441 
442  for (int dd = 0; dd != nb_dofs; dd++) {
443  a_auxFun.referenceCoords[dd] <<= getCoords()[dd];
444  }
445  for (int dd = 0; dd != nb_dofs; dd++) {
446  a_auxFun.currentCoords[dd] <<= data.getFieldData()[dd];
447  }
448 
449  a_auxFun.crackAreaIncrement.resize(9, false);
450  std::fill(a_auxFun.crackAreaIncrement.begin(),
451  a_auxFun.crackAreaIncrement.end(), 0);
452 
453  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
454 
455  double val = getGaussPts()(2, gg) * 0.5;
456 
457  CHKERR a_auxFun.matrixB(gg, data);
458  CHKERR a_auxFun.dIffX();
459  CHKERR a_auxFun.nOrmal();
460  CHKERR a_auxFun.matrixA();
461  CHKERR a_auxFun.calculateAreaIncrementDirection(val);
462  }
463 
464  // cerr << "Griffith Force " << auxFun.griffithForce << endl;
465 
466  commonData.crackAreaIncrement.resize(nb_dofs, false);
467  for (int dd = 0; dd != nb_dofs; dd++) {
468  a_auxFun.crackAreaIncrement[dd] >>= commonData.crackAreaIncrement[dd];
469  }
470 
471  trace_off();
472 
473  isTapeRecorded = true;
474 
475  MoFEMFunctionReturn(0);
476  }
477  };
478 
479  struct AuxOp {
480 
481  int tAg;
482  CommonData &commonData;
483  AuxOp(int tag, CommonData &common_data)
484  : tAg(tag), commonData(common_data) {}
485 
486  ublas::vector<int> rowIndices;
487  VectorDouble activeVariables;
488 
489  PetscErrorCode
490  setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
491  int side, EntityType type,
492  DataForcesAndSourcesCore::EntData &data) {
493  MoFEMFunctionBegin;
494  int nb_dofs = data.getIndices().size();
495  activeVariables.resize(18);
496  for (int dd = 0; dd != nb_dofs; dd++) {
497  activeVariables[dd] = fe_ptr->getCoords()[dd];
498  }
499  for (int dd = 0; dd != nb_dofs; dd++) {
500  activeVariables[9 + dd] = data.getFieldData()[dd];
501  }
502  MoFEMFunctionReturn(0);
503  }
504  };
505 
506  struct OpAreaC : public FaceElementForcesAndSourcesCore::UserDataOperator,
507  AuxOp {
508 
509  OpAreaC(int tag, CommonData &common_data)
510  : FaceElementForcesAndSourcesCore::UserDataOperator(
511  "LAMBDA_CRACKFRONT_AREA", "MESH_NODE_POSITIONS",
512  UserDataOperator::OPROWCOL),
513  AuxOp(tag, common_data) {
514  sYmm = false;
515  }
516 
517  MatrixDouble cMat;
518 
519  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
520  EntityType col_type,
521  DataForcesAndSourcesCore::EntData &row_data,
522  DataForcesAndSourcesCore::EntData &col_data) {
523  MoFEMFunctionBegin;
524 
525  if (row_type != MBVERTEX) {
526  MoFEMFunctionReturnHot(0);
527  }
528  VectorInt &row_indices = row_data.getIndices();
529  const int nb_rows = row_indices.size();
530  if (!nb_rows)
531  MoFEMFunctionReturnHot(0);
532  VectorInt &col_indices = col_data.getIndices();
533  const int nb_cols = col_indices.size();
534  if (!nb_cols)
535  MoFEMFunctionReturnHot(0);
536 
537  const EntityHandle *conn = getConn();
538  for (int nn = 0; nn != 3; ++nn) {
539  if (row_indices[nn] != -1) {
540  if (row_data.getFieldDofs()[nn]->getEnt() != conn[nn]) {
541  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
542  "Data inconsistency");
543  }
544  }
545  for (int dd = 0; dd != 3; ++dd) {
546  if (col_indices[3 * nn + dd] != -1) {
547  if (col_data.getFieldDofs()[3 * nn + dd]->getEnt() != conn[nn]) {
548  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
549  "Data inconsistency");
550  }
551  }
552  }
553  }
554 
555  auto row_dofs = getFEMethod()->getRowDofsPtr();
556  for (auto it = row_dofs->begin(); it != row_dofs->end(); ++it) {
557  int side = it->get()->getSideNumber();
558  if (conn[side] != it->get()->getEnt()) {
559  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
560  "Data inconsistency");
561  }
562  }
563 
564  auto col_dofs = getFEMethod()->getColDofsPtr();
565  for (auto it = col_dofs->begin(); it != col_dofs->end(); ++it) {
566  int side = it->get()->getSideNumber();
567  if (conn[side] != it->get()->getEnt()) {
568  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
569  "Data inconsistency");
570  }
571  }
572 
573  cMat.resize(3, 9, false);
574  cMat.clear();
575 
576  CHKERR setVariables(this, col_side, col_type, col_data);
577 
578  int nb_dofs = col_data.getFieldData().size();
579  commonData.crackAreaIncrement.resize(nb_dofs, false);
580  int r;
581  // play recorder for values
582  r = ::function(tAg, nb_dofs, 18, &activeVariables[0],
583  &commonData.crackAreaIncrement[0]);
584  if (r < 3) {
585  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
586  "ADOL-C function evaluation with error r = %d", r);
587  }
588 
589  auto check_isnormal = [](double v) {
590  MoFEMFunctionBeginHot;
591  if (v != 0 && !std::isnormal(v)) {
592  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
593  "Not a number %3.4e", v);
594  }
595  MoFEMFunctionReturnHot(0);
596  };
597 
598  for (int rr = 0; rr != nb_rows; rr++) {
599  if (row_indices[rr] != -1) {
600  CHKERR check_isnormal(commonData.crackAreaIncrement[3 * rr + 0]);
601  CHKERR check_isnormal(commonData.crackAreaIncrement[3 * rr + 1]);
602  CHKERR check_isnormal(commonData.crackAreaIncrement[3 * rr + 2]);
603  cMat(rr, 3 * rr + 0) = commonData.crackAreaIncrement[3 * rr + 0];
604  cMat(rr, 3 * rr + 1) = commonData.crackAreaIncrement[3 * rr + 1];
605  cMat(rr, 3 * rr + 2) = commonData.crackAreaIncrement[3 * rr + 2];
606  }
607  }
608 
609  if (!commonData.setMapV) {
610  CHKERR MatSetValues(getFEMethod()->snes_B, nb_rows, &row_indices[0],
611  nb_cols, &col_indices[0], &cMat(0, 0), ADD_VALUES);
612  } else {
613  for (int rr = 0; rr != nb_rows; ++rr) {
614  if (row_indices[rr] != -1) {
615  if (commonData.mapV.find(row_indices[rr]) ==
616  commonData.mapV.end()) {
617  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
618  "Vector for given DOF in map not found");
619  }
620  for (int cc = 0; cc != nb_cols; ++cc) {
621  if (col_indices[cc] != -1) {
622  CHKERR check_isnormal(cMat(rr, cc));
623  }
624  }
625  CHKERR VecSetValues(commonData.mapV[row_indices[rr]], nb_cols,
626  &col_indices[0], &cMat(rr, 0), ADD_VALUES);
627  }
628  }
629  }
630 
631  MoFEMFunctionReturn(0);
632  }
633  };
634 
635  struct OpTangentJacobian : public OpAreaJacobian {
636 
637  OpTangentJacobian(int tag, CommonData &common_data)
638  : OpAreaJacobian(tag, common_data) {
639  sYmm = false;
640  }
641 
642  PetscErrorCode doWork(int side, EntityType type,
643  DataForcesAndSourcesCore::EntData &data) {
644  MoFEMFunctionBegin;
645 
646  if (isTapeRecorded) {
647  MoFEMFunctionReturnHot(0);
648  }
649  if (type != MBVERTEX) {
650  MoFEMFunctionReturnHot(0);
651  }
652 
653  int nb_dofs = data.getFieldData().size();
654  int nb_gauss_pts = data.getN().size1();
655 
656  a_auxFun.referenceCoords.resize(9, false);
657  a_auxFun.currentCoords.resize(9, false);
658 
659  trace_on(tAg);
660 
661  for (int dd = 0; dd != nb_dofs; dd++) {
662  a_auxFun.referenceCoords[dd] <<= getCoords()[dd];
663  }
664  for (int dd = 0; dd != nb_dofs; dd++) {
665  a_auxFun.currentCoords[dd] <<= data.getFieldData()[dd];
666  }
667 
668  a_auxFun.crackAreaIncrement.resize(9, false);
669  a_auxFun.crackFrontTangent.resize(9, false);
670  std::fill(a_auxFun.crackAreaIncrement.begin(),
671  a_auxFun.crackAreaIncrement.end(), 0);
672  std::fill(a_auxFun.crackFrontTangent.begin(),
673  a_auxFun.crackFrontTangent.end(), 0);
674 
675  for (int gg = 0; gg != nb_gauss_pts; gg++) {
676 
677  double val = getGaussPts()(2, gg) * 0.5;
678 
679  CHKERR a_auxFun.matrixB(gg, data);
680  CHKERR a_auxFun.dIffX();
681  CHKERR a_auxFun.nOrmal();
682  CHKERR a_auxFun.matrixA();
683  CHKERR a_auxFun.calculateAreaIncrementDirection(1);
684  CHKERR a_auxFun.calculateFrontTangent(val);
685  }
686 
687  commonData.crackFrontTangent.resize(nb_dofs, false);
688  for (int dd = 0; dd != nb_dofs; dd++) {
689  a_auxFun.crackFrontTangent[dd] >>= commonData.crackFrontTangent[dd];
690  }
691 
692  trace_off();
693 
694  isTapeRecorded = true;
695 
696  MoFEMFunctionReturn(0);
697  }
698  };
699 
700  struct OpTangentC : public FaceElementForcesAndSourcesCore::UserDataOperator,
701  AuxOp {
702 
703  OpTangentC(int tag, CommonData &common_data)
704  : FaceElementForcesAndSourcesCore::UserDataOperator(
705  "LAMBDA_CRACKFRONT_AREA_TANGENT", "MESH_NODE_POSITIONS",
706  UserDataOperator::OPROWCOL),
707  AuxOp(tag, common_data) {
708  sYmm = false;
709  }
710 
711  MatrixDouble cMat;
712 
713  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
714  EntityType col_type,
715  DataForcesAndSourcesCore::EntData &row_data,
716  DataForcesAndSourcesCore::EntData &col_data) {
717  MoFEMFunctionBegin;
718 
719  if (row_type != MBVERTEX) {
720  MoFEMFunctionReturnHot(0);
721  }
722 
723  VectorInt &row_indices = row_data.getIndices();
724  const int nb_rows = row_indices.size();
725  if (!nb_rows)
726  MoFEMFunctionReturnHot(0);
727  VectorInt &col_indices = col_data.getIndices();
728  const int nb_cols = col_indices.size();
729  if (!nb_cols)
730  MoFEMFunctionReturnHot(0);
731  cMat.resize(3, 9, false);
732  cMat.clear();
733 
734  CHKERR setVariables(this, col_side, col_type, col_data);
735 
736  int nb_dofs = col_data.getFieldData().size();
737  commonData.crackFrontTangent.resize(nb_dofs, false);
738  int r;
739  // play recorder for values
740  r = ::function(tAg, nb_dofs, 18, &activeVariables[0],
741  &commonData.crackFrontTangent[0]);
742  if (r < 3) {
743  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
744  "ADOL-C function evaluation with error r = %d", r);
745  }
746 
747  // Set face orientation
748  {
749  // Get orientation in respect to adjacent tet
750  const EntityHandle face = getFEEntityHandle();
751  const BitRefLevel bit = getFEMethod()->problemPtr->getBitRefLevel();
752  Range adj_side_elems;
753  CHKERR commonData.mField.getInterface<BitRefManager>()->getAdjacencies(
754  bit, &face, 1, 3, adj_side_elems);
755  adj_side_elems = adj_side_elems.subset_by_type(MBTET);
756  if (adj_side_elems.size() != 1) {
757  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
758  "expect 1 tet but is %u", adj_side_elems.size());
759  }
760  EntityHandle side_elem = *adj_side_elems.begin();
761  int side_number, sense, offset;
762  CHKERR commonData.mField.get_moab().side_number(
763  side_elem, face, side_number, sense, offset);
764  if (sense == -1) {
765  commonData.crackFrontTangent *= -1;
766  }
767  // Get orientaton of face in respect to face in the prism
768  int side;
769  CHKERR commonData.mField.get_moab().tag_get_data(
770  commonData.thInterfaceSide, &face, 1, &side);
771  if (side == 1) {
772  commonData.crackFrontTangent *= -1;
773  }
774  }
775 
776  auto check_isnormal = [](double v) {
777  MoFEMFunctionBeginHot;
778  if (v != 0 && !std::isnormal(v)) {
779  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Not a number");
780  }
781  MoFEMFunctionReturnHot(0);
782  };
783 
784  for (int rr = 0; rr != nb_rows; ++rr) {
785  if (row_indices[rr] != -1) {
786  // cerr << commonData.crackFrontTangent(3*rr+0) << " "
787  // << commonData.crackFrontTangent(3*rr+1) << " "
788  // << commonData.crackFrontTangent(3*rr+2) << endl;
789  CHKERR check_isnormal(commonData.crackFrontTangent[3 * rr + 0]);
790  CHKERR check_isnormal(commonData.crackFrontTangent[3 * rr + 1]);
791  CHKERR check_isnormal(commonData.crackFrontTangent[3 * rr + 2]);
792  cMat(rr, 3 * rr + 0) = commonData.crackFrontTangent[3 * rr + 0];
793  cMat(rr, 3 * rr + 1) = commonData.crackFrontTangent[3 * rr + 1];
794  cMat(rr, 3 * rr + 2) = commonData.crackFrontTangent[3 * rr + 2];
795  }
796  }
797 
798  if (!commonData.setMapV) {
799  CHKERR MatSetValues(getFEMethod()->snes_B, nb_rows, &row_indices[0],
800  nb_cols, &col_indices[0], &cMat(0, 0), ADD_VALUES);
801  } else {
802  for (int rr = 0; rr != nb_rows; ++rr) {
803  if (row_indices[rr] != -1) {
804  if (commonData.mapV.find(row_indices[rr]) ==
805  commonData.mapV.end()) {
806  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
807  "Vector for given DOF in map not found");
808  }
809  for (int cc = 0; cc != nb_cols; ++cc) {
810  if (col_indices[cc] != -1) {
811  CHKERR check_isnormal(cMat(rr, cc));
812  }
813  }
814  CHKERR VecSetValues(commonData.mapV[row_indices[rr]], nb_cols,
815  &col_indices[0], &cMat(rr, 0), ADD_VALUES);
816  }
817  }
818  }
819 
820  MoFEMFunctionReturn(0);
821  }
822  };
823 };
824 
825 } // namespace FractureMechanics
826 #endif //__CONSTANT_AREA_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
FractureMechanics::ConstantArea::CommonData
Definition: ConstantArea.hpp:32
FractureMechanics::ConstantArea::MyTriangleFE::postProcess
PetscErrorCode postProcess()
Definition: ConstantArea.hpp:132
_IT_NUMEREDDOF_COL_BY_BITNUMBER_FOR_LOOP_
#define _IT_NUMEREDDOF_COL_BY_BITNUMBER_FOR_LOOP_(PROBLEMPTR, FIELD_BIT_NUMBER, IT)
use with loops to iterate col DOFs
Definition: ProblemsMultiIndices.hpp:358
FractureMechanics::ConstantArea::feRhsPtr
boost::shared_ptr< MyTriangleFE > feRhsPtr
Definition: ConstantArea.hpp:237
FractureMechanics::ConstantArea::CommonData::setMapV
bool setMapV
Definition: ConstantArea.hpp:46
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::matrixA
PetscErrorCode matrixA()
Definition: ConstantArea.hpp:363
FractureMechanics::ConstantArea::CommonData::crackAreaIncrement
VectorDouble crackAreaIncrement
Definition: ConstantArea.hpp:37
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::dIffX
PetscErrorCode dIffX()
Definition: ConstantArea.hpp:326
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
FractureMechanics::ConstantArea::MyTriangleFE::petscQ
Mat petscQ
Definition: ConstantArea.hpp:59
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::Beta
MatrixDouble Beta
Definition: ConstantArea.hpp:271
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentSpinKsi
ublas::matrix< TYPE > currentSpinKsi
Definition: ConstantArea.hpp:285
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceCoords
ublas::vector< TYPE > referenceCoords
Definition: ConstantArea.hpp:273
FractureMechanics::ConstantArea::OpAreaC::OpAreaC
OpAreaC(int tag, CommonData &common_data)
Definition: ConstantArea.hpp:509
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceXdKsi
ublas::vector< TYPE > referenceXdKsi
Definition: ConstantArea.hpp:275
FractureMechanics::ConstantArea::OpAreaJacobian::OpAreaJacobian
OpAreaJacobian(int tag, CommonData &common_data)
Definition: ConstantArea.hpp:247
FractureMechanics::ConstantArea::OpTangentC
Definition: ConstantArea.hpp:700
FractureMechanics::ConstantArea::CommonData::tangentCrackFrontTangent
MatrixDouble tangentCrackFrontTangent
Definition: ConstantArea.hpp:42
FractureMechanics::ConstantArea
Constant area constrains.
Definition: ConstantArea.hpp:26
FractureMechanics::ConstantArea::MyTriangleFE::preProcess
PetscErrorCode preProcess()
Definition: ConstantArea.hpp:68
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceArea
TYPE referenceArea
Definition: ConstantArea.hpp:280
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
FractureMechanics::ConstantArea::CommonData::mField
MoFEM::Interface & mField
Definition: ConstantArea.hpp:34
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentXdEta
ublas::vector< TYPE > currentXdEta
Definition: ConstantArea.hpp:284
FractureMechanics::ConstantArea::CommonData::thInterfaceSide
Tag thInterfaceSide
Definition: ConstantArea.hpp:48
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions
Definition: ConstantArea.hpp:267
FractureMechanics::ConstantArea::OpAreaJacobian
Definition: ConstantArea.hpp:240
FractureMechanics::ConstantArea::CommonData::tangentCrackAreaIncrement
MatrixDouble tangentCrackAreaIncrement
Definition: ConstantArea.hpp:38
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::crackFrontTangent
ublas::vector< TYPE > crackFrontTangent
Definition: ConstantArea.hpp:292
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::N
MatrixDouble N
Definition: ConstantArea.hpp:269
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentArea
TYPE currentArea
Definition: ConstantArea.hpp:288
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
sdf.r
int r
Definition: sdf.py:8
FractureMechanics::ConstantArea::CommonData::mapV
map< int, Vec > mapV
Definition: ConstantArea.hpp:45
FTensor::Tensor2< TYPE, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FractureMechanics::ConstantArea::feLhsPtr
boost::shared_ptr< MyTriangleFE > feLhsPtr
Definition: ConstantArea.hpp:238
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::crackAreaIncrement
ublas::vector< TYPE > crackAreaIncrement
Definition: ConstantArea.hpp:291
FractureMechanics::ConstantArea::MyTriangleFE::getRule
int getRule(int order)
Definition: ConstantArea.hpp:66
FractureMechanics::ConstantArea::MyTriangleFE::MyTriangleFE
MyTriangleFE(MoFEM::Interface &m_field, CommonData &common_data, std::string field_name="LAMBDA_CRACKFRONT_AREA")
Definition: ConstantArea.hpp:61
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::calculateFrontTangent
PetscErrorCode calculateFrontTangent(double beta)
Definition: ConstantArea.hpp:393
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::matrixB
PetscErrorCode matrixB(int gg, DataForcesAndSourcesCore::EntData &data)
Definition: ConstantArea.hpp:307
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
_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::ConstantArea::OpTangentJacobian
Definition: ConstantArea.hpp:635
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::A
ublas::matrix< TYPE > A
Definition: ConstantArea.hpp:290
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceSpinEta
ublas::matrix< TYPE > referenceSpinEta
Definition: ConstantArea.hpp:278
FractureMechanics::ConstantArea::AuxOp::tAg
int tAg
Definition: ConstantArea.hpp:481
FractureMechanics::ConstantArea::MyTriangleFE::petscB
Mat petscB
Definition: ConstantArea.hpp:57
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
FractureMechanics::ConstantArea::CommonData::tangentCrackAreaIncrementRowPtr
vector< double * > tangentCrackAreaIncrementRowPtr
Definition: ConstantArea.hpp:39
FractureMechanics::ConstantArea::MyTriangleFE::petscF
Vec petscF
Definition: ConstantArea.hpp:58
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentCoords
ublas::vector< TYPE > currentCoords
Definition: ConstantArea.hpp:282
FractureMechanics::ConstantArea::OpTangentC::OpTangentC
OpTangentC(int tag, CommonData &common_data)
Definition: ConstantArea.hpp:703
COL
@ COL
Definition: definitions.h:123
FractureMechanics::ConstantArea::AuxOp::AuxOp
AuxOp(int tag, CommonData &common_data)
Definition: ConstantArea.hpp:483
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::NTN
MatrixDouble NTN
Definition: ConstantArea.hpp:269
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceXdEta
ublas::vector< TYPE > referenceXdEta
Definition: ConstantArea.hpp:276
FractureMechanics::ConstantArea::OpAreaJacobian::a_auxFun
AuxFunctions< adouble > a_auxFun
Definition: ConstantArea.hpp:421
FractureMechanics::ConstantArea::MyTriangleFE
Definition: ConstantArea.hpp:52
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
FractureMechanics::ConstantArea::CommonData::tangentCrackFrontTangentRowPtr
vector< double * > tangentCrackFrontTangentRowPtr
Definition: ConstantArea.hpp:43
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::nOrmal
PetscErrorCode nOrmal()
Definition: ConstantArea.hpp:343
FractureMechanics::ConstantArea::OpAreaJacobian::doWork
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: ConstantArea.hpp:423
FractureMechanics::ConstantArea::AuxOp::activeVariables
VectorDouble activeVariables
Definition: ConstantArea.hpp:487
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FractureMechanics::ConstantArea::OpAreaJacobian::commonData
CommonData & commonData
Definition: ConstantArea.hpp:244
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::tSpin
FTensor::Tensor2< TYPE, 3, 3 > tSpin
Definition: ConstantArea.hpp:391
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::sPin
PetscErrorCode sPin(ublas::matrix< TYPE > &spin, ublas::vector< TYPE > &vec)
Definition: ConstantArea.hpp:294
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FractureMechanics::ConstantArea::mField
MoFEM::Interface & mField
Definition: ConstantArea.hpp:28
FTensor::Index< 'i', 3 >
FractureMechanics::ConstantArea::MyTriangleFE::fieldName
std::string fieldName
Definition: ConstantArea.hpp:55
FractureMechanics::ConstantArea::AuxOp::commonData
CommonData & commonData
Definition: ConstantArea.hpp:482
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceSpinKsi
ublas::matrix< TYPE > referenceSpinKsi
Definition: ConstantArea.hpp:277
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::calculateAreaIncrementDirection
PetscErrorCode calculateAreaIncrementDirection(double beta)
Definition: ConstantArea.hpp:377
FractureMechanics::ConstantArea::MyTriangleFE::commonData
CommonData & commonData
Definition: ConstantArea.hpp:54
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentXdKsi
ublas::vector< TYPE > currentXdKsi
Definition: ConstantArea.hpp:283
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
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::ConstantArea::ConstantArea
ConstantArea(MoFEM::Interface &m_field)
Definition: ConstantArea.hpp:29
FractureMechanics::ConstantArea::CommonData::CommonData
CommonData(MoFEM::Interface &m_field)
Definition: ConstantArea.hpp:35
FractureMechanics::ConstantArea::OpAreaC::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: ConstantArea.hpp:519
FractureMechanics::ConstantArea::AuxOp::setVariables
PetscErrorCode setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: ConstantArea.hpp:490
FractureMechanics::ConstantArea::OpTangentJacobian::doWork
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: ConstantArea.hpp:642
FractureMechanics::ConstantArea::AuxOp::rowIndices
ublas::vector< int > rowIndices
Definition: ConstantArea.hpp:486
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentSpinEta
ublas::matrix< TYPE > currentSpinEta
Definition: ConstantArea.hpp:286
FractureMechanics::ConstantArea::OpTangentC::cMat
MatrixDouble cMat
Definition: ConstantArea.hpp:711
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FractureMechanics::ConstantArea::AuxOp
Definition: ConstantArea.hpp:479
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FractureMechanics::ConstantArea::OpAreaC::cMat
MatrixDouble cMat
Definition: ConstantArea.hpp:517
FractureMechanics::ConstantArea::commonData
CommonData commonData
Definition: ConstantArea.hpp:50
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
FractureMechanics::ConstantArea::OpTangentJacobian::OpTangentJacobian
OpTangentJacobian(int tag, CommonData &common_data)
Definition: ConstantArea.hpp:637
FractureMechanics::ConstantArea::OpAreaC
Definition: ConstantArea.hpp:506
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::referenceNormal
ublas::vector< TYPE > referenceNormal
Definition: ConstantArea.hpp:279
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FractureMechanics::ConstantArea::OpAreaJacobian::isTapeRecorded
bool isTapeRecorded
Definition: ConstantArea.hpp:245
FractureMechanics::ConstantArea::OpAreaJacobian::tAg
int tAg
Definition: ConstantArea.hpp:243
ConstrainMatrixCtx
structure for projection matrices
Definition: ConstrainMatrixCtx.hpp:16
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::currentNormal
ublas::vector< TYPE > currentNormal
Definition: ConstantArea.hpp:287
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
FractureMechanics
Definition: AnalyticalFun.hpp:15
FractureMechanics::ConstantArea::OpAreaJacobian::AuxFunctions::Bksi
MatrixDouble Bksi
Definition: ConstantArea.hpp:270
FractureMechanics::ConstantArea::OpTangentC::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: ConstantArea.hpp:713
FractureMechanics::ConstantArea::CommonData::crackFrontTangent
VectorDouble crackFrontTangent
Definition: ConstantArea.hpp:41
Generated by Doxygen 1.8.17 and hosted at University of Glasgow