v0.14.0
Loading...
Searching...
No Matches
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
22namespace FractureMechanics {
23
24static double calMax(double a, double b, double r) {
25 return (a + b + (1 / r) * pow(fabs(a - b), r)) / 2;
26}
27
28static 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 */
36
38
40
41 Mat B;
42 Vec F;
43
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
58
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) {
85
86 MOFEM_LOG_CHANNEL("WORLD");
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
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
133
139
140 boost::function<FTensor::Tensor1<TYPE *, 3>(ublas::vector<TYPE> &v)>
142
143 boost::function<FTensor::Tensor1<FTensor::PackPtr<const double *, 2>, 2>(
144 const MatrixDouble &m)>
146
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) {
155 &m(0, 0), &m(0, 1));
156 };
157 }
158
159 MoFEMErrorCode calculateAreaAndNormal(const MatrixDouble &diff_n) {
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) =
175 currentArea = 0.5 * sqrt(t_Normal(i) * t_Normal(i));
176
178 }
179
180 MoFEMErrorCode calculateMatA(const MatrixDouble &diff_n) {
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
199 }
200
201 MoFEMErrorCode calculateGriffithForce(const double gc, const double beta,
202 const MatrixDouble &diff_n) {
204
206 CHKERR calculateMatA(diff_n);
207
208 for (int dd = 0; dd != 9; dd++)
209 griffithForce[dd] *= gc * beta;
210
212 }
213 };
214
215 struct AuxOp {
216
217 int tAg;
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
229 VectorDouble3 lambdaAtNodes;
230 VectorDouble activeVariables;
231
232 MoFEMErrorCode setIndices(DataForcesAndSourcesCore::EntData &data) {
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 }
258 }
259
260 MoFEMErrorCode
261 setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
262 DataForcesAndSourcesCore::EntData &data) {
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 }
273 }
274
275 MoFEMErrorCode
276 setLambdaNodes(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
277 const std::string &lambda_field_name) {
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 }
298 }
299
300 MoFEMErrorCode
301 setLambdaIndices(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
302 const std::string &lambda_field_name) {
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 }
327 }
328 };
329
330 struct OpRhs : public FaceElementForcesAndSourcesCore::UserDataOperator,
331 AuxOp {
332 OpRhs(int tag, BlockData &block_data, CommonData &common_data)
334 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
335 AuxOp(tag, block_data, common_data) {}
336 MoFEMErrorCode doWork(int side, EntityType type,
337 DataForcesAndSourcesCore::EntData &data) {
339
340 if (type != MBVERTEX)
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);
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 */
379 : public FaceElementForcesAndSourcesCore::UserDataOperator,
380 AuxOp {
382 CommonData &common_data)
384 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
385 AuxOp(tag, block_data, common_data) {}
387 MoFEMErrorCode doWork(int side, EntityType type,
388 DataForcesAndSourcesCore::EntData &data) {
390
391 if (type != MBVERTEX)
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;
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
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 */
433 : public FaceElementForcesAndSourcesCore::UserDataOperator,
434 AuxOp {
435 const double betaGC;
436 std::string rhoTagName;
437 boost::shared_ptr<MWLSApprox> mwlsApprox;
439
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)
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) {
454
455 if (type != MBVERTEX)
457 const int row_nb_dofs = row_data.getIndices().size();
458 if (!row_nb_dofs)
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) {
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
527 element_prod(commonData.griffithForce, gc_coeff);
529 }
530 };
531
533
535
536 Tag tH;
539
541 Tag *th_coords = nullptr)
542 : mField(m_field), tH(th), tRis(tris), thCoordsPtr(th_coords) {}
543
545
546 MoFEMErrorCode operator()() {
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
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)
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) {
619
620 if (row_type != MBVERTEX) {
622 }
623
631
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) {
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
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
715 }
716 };
717
719 : public FaceElementForcesAndSourcesCore::UserDataOperator,
720 AuxOp {
721
722 OpHeterogeneousGcLhs(int tag, BlockData &block_data,
723 CommonData &common_data, const double beta_gc)
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) {
737
738 if (row_type != MBVERTEX) {
740 }
741
744
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) {
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;
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
808 }
809 };
810
812
813 const std::string lambdaFieldName;
815 SmartPetscObj<Vec> deltaVec;
816
818 const std::string &lambda_field_name)
819 : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name) {}
820
821 MoFEMErrorCode preProcess() {
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);
833 }
834
835 MoFEMErrorCode postProcess() {
837 MOFEM_LOG_CHANNEL("WORLD");
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);
851 }
852 };
853
855 : public FaceElementForcesAndSourcesCore::UserDataOperator,
856 AuxOp {
857
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)
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
872
873 MoFEMErrorCode doWork(int side, EntityType type,
874 DataForcesAndSourcesCore::EntData &data) {
876 if (type != MBVERTEX) {
878 }
879 // get indices of material DOFs
880 CHKERR setIndices(data);
881 // get indices of Lagrange multiplier DOFs
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
912 }
913 };
914
916
917 const std::string lambdaFieldName;
919 SmartPetscObj<Vec> deltaVec;
920
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() {
929
930 MOFEM_LOG_CHANNEL("WORLD");
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 =
940 switch (snes_ctx) {
941 case CTX_SNESSETFUNCTION: {
943 bit_field_number, dit)) {
944 if (static_cast<int>(dit->get()->getPart()) ==
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: {
957 bit_field_number, dit)) {
958 if (static_cast<int>(dit->get()->getPart()) ==
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);
985 }
986
987 MoFEMErrorCode postProcess() {
990 }
991 };
992
994 : public FaceElementForcesAndSourcesCore::UserDataOperator,
995 AuxOp {
996
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)
1004 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
1005 AuxOp(tag, block_data, common_data), mField(m_field),
1006 lambdaFieldName(lambda_field_name) {}
1007
1009 VectorDouble9 nF;
1010 // VectorDouble3 nG;
1011
1012 MoFEMErrorCode doWork(int side, EntityType type,
1013 DataForcesAndSourcesCore::EntData &data) {
1015 if (type != MBVERTEX)
1017
1018 // get indices of material DOFs
1019 CHKERR setIndices(data);
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
1047 }
1048 };
1049
1051 : public FaceElementForcesAndSourcesCore::UserDataOperator,
1052 AuxOp {
1053
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)
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;
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) {
1088
1089 if (row_type != MBVERTEX || col_type != MBVERTEX)
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);
1241 for (int dd = 0; dd != 9; ++dd) {
1244 }
1245
1246 {
1247 int r;
1248 // play recorder for jacobian
1249 r = ::jacobian(tAg, 9, 21, &activeVariables[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
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 */
1289
1290 int tAg;
1291 boost::shared_ptr<ArcLengthCtx> arcPtr;
1294 std::string lambdaFieldName;
1295
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,
1314
1316 std::string lambdaFieldName;
1318
1322 const std::string &lambda_field_name,
1323 FrontArcLengthControl &arc_front)
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
1333
1334 MoFEMErrorCode doWork(int side, EntityType type,
1335 DataForcesAndSourcesCore::EntData &data) {
1337 if (type != MBVERTEX) {
1339 }
1340 // get indices of Lagrange multiplier DOFs
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] *
1378 }
1379 }
1380 if (delta != delta) {
1381 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "nan value");
1382 }
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;
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
1408 VectorDouble gradDelta;
1409 VectorDouble activeVariables;
1410
1411 MoFEMErrorCode doWork(int side, EntityType type,
1412 DataForcesAndSourcesCore::EntData &data) {
1414 if (type != MBVERTEX) {
1416 }
1417 CHKERR setIndices(data);
1418 // get indices of Lagrange multiplier DOFs
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
1494 }
1495 };
1496
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),
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,
1510 CHKERRABORT(PETSC_COMM_WORLD, ierr);
1511 feIntLambda = boost::shared_ptr<MyTriangleFE>(
1513 }
1515 ierr = VecDestroy(&ghostIntLambda);
1516 CHKERRABORT(PETSC_COMM_WORLD, ierr);
1517 }
1518
1519 MoFEMErrorCode preProcess() {
1521 switch (snes_ctx) {
1522 case CTX_SNESSETFUNCTION: {
1523 CHKERR VecAssemblyBegin(snes_f);
1524 CHKERR VecAssemblyEnd(snes_f);
1527 } break;
1528 default:
1529 break;
1530 }
1532 }
1533
1534 /** \brief Calculate internal lambda
1535 */
1538 return intLambda;
1540 }
1541
1542 MoFEMErrorCode operator()() {
1544 MOFEM_LOG_CHANNEL("WORLD");
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 }
1565 }
1566
1567 MoFEMErrorCode postProcess() {
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 }
1584 }
1585
1586 /** \brief Calculate db
1587 */
1588 MoFEMErrorCode calculateDb() {
1590
1591 MOFEM_LOG_CHANNEL("WORLD");
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);
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
1648 "WORLD", Sev::noisy,
1649 "\tintLambdaLambda = %9.8e intLambdaDelta = %9.8e intLambda = %9.8e",
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
1657 }
1658
1659 MoFEMErrorCode calculateDxAndDlambda(Vec x) {
1661 MOFEM_LOG_CHANNEL("WORLD");
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);
1713 }
1714 };
1715};
1716
1717} // namespace FractureMechanics
1718
1719#endif //__GRIFFITH_FORCE_ELEMENT_HPP__
static Number< 1 > N1
static Number< 0 > N0
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
constexpr double a
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
double W
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:325
#define _IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(PROBLEMPTR, FIELD_BIT_NUMBER, IT)
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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
static double diffCalMax_a(double a, double b, double r)
static double calMax(double a, double b, double r)
double rho
Definition: plastic.cpp:191
static constexpr double delta
boost::function< FTensor::Tensor1< TYPE *, 3 >(ublas::vector< TYPE > &v)> get_ftensor_from_vec
MoFEMErrorCode calculateMatA(const MatrixDouble &diff_n)
MoFEMErrorCode calculateGriffithForce(const double gc, const double beta, const MatrixDouble &diff_n)
MoFEMErrorCode calculateAreaAndNormal(const MatrixDouble &diff_n)
boost::function< FTensor::Tensor1< FTensor::PackPtr< const double *, 2 >, 2 >(const MatrixDouble &m)> get_ftensor_from_mat_2d
MoFEMErrorCode setLambdaNodes(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, const std::string &lambda_field_name)
AuxOp(int tag, BlockData &block_data, CommonData &common_data)
MoFEMErrorCode setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode setLambdaIndices(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr, const std::string &lambda_field_name)
MoFEMErrorCode setIndices(DataForcesAndSourcesCore::EntData &data)
vector< double * > tangentGriffithForceRowPtr
Pointers to rows of tangent matrix.
MatrixDouble diffRho
for lhs with heterogeneous gc
VectorDouble griffithForce
Griffith force at nodes.
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)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpIntLambda(MoFEM::Interface &m_field, GriffithForceElement::BlockData &block_data, GriffithForceElement::CommonData &common_data, const std::string &lambda_field_name, FrontArcLengthControl &arc_front)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
FrontArcLengthControl(int tag, GriffithForceElement::BlockData &block_data, GriffithForceElement::CommonData &common_data, const std::string &lambda_field_name, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
MyTriangleFEConstrainsDelta(MoFEM::Interface &m_field, const std::string &lambda_field_name)
MyTriangleFEConstrains(MoFEM::Interface &m_field, const std::string &lambda_field_name, BlockData &block_data, SmartPetscObj< Vec > &delta_vec)
OpCalculateGriffithForce(int tag, BlockData &block_data, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &row_data)
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)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpConstrainsDelta(MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data, const std::string &lambda_field_name, SmartPetscObj< Vec > &delta_vec)
OpConstrainsLhs(MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data, const std::string &lambda_field_name, SmartPetscObj< Vec > &delta_vec)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpConstrainsRhs(MoFEM::Interface &m_field, int tag, BlockData &block_data, CommonData &common_data, const std::string &lambda_field_name)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpHeterogeneousGcLhs(int tag, BlockData &block_data, CommonData &common_data, const double beta_gc)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpLhs(int tag, BlockData &block_data, CommonData &common_data, const double beta_gc=0)
OpRhs(int tag, BlockData &block_data, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
SaveGriffithForceOnTags(MoFEM::Interface &m_field, Tag th, Range &tris, Tag *th_coords=nullptr)
Implementation of Griffith element.
map< int, BlockData > blockData
Shared data between finite elements.
boost::shared_ptr< MyTriangleFE > feRhsPtr
boost::shared_ptr< MyTriangleFE > feLhsPtr
static MoFEMErrorCode getElementOptions(BlockData &block_data)
CommonData commonData
Common data sgared between operators.
GriffithForceElement(MoFEM::Interface &m_field)
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
Deprecated interface functions.