v0.13.1
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;
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
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
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
155 &m(0, 0), &m(0, 1));
156 };
157 }
158
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
182
183 griffithForce.resize(9, false);
184 griffithForce.clear();
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
231
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
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
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
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
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) {}
339
340 if (type != MBVERTEX)
342 CHKERR setIndices(data);
343 CHKERR setVariables(this, data);
344 int nb_dofs = data.getIndices().size();
345
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 */
380 AuxOp {
382 CommonData &common_data)
384 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
385 AuxOp(tag, block_data, common_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);
413
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 */
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
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 commonData.diffRho.clear();
480 auto t_coords = get_ftensor_from_vec(row_data.getFieldData());
481
482 commonData.densityRho.resize(nb_dofs / 3, false);
483 commonData.densityRho.clear();
485 VectorDouble gc_coeff(nb_dofs, false);
486 auto t_coeff = get_ftensor_from_vec(gc_coeff);
487 auto t_approx_base_point =
488 getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
489 const int nb_gauss_pts = row_data.getN().size1();
490
491 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
492 for (int dd = 0; dd != nb_dofs / 3; ++dd) {
493
494 // go only through nodes at the crack front
495 if (rowIndices[dd * 3] != -1) {
497 t_pos(i) = t_coords(i);
498
499 for (int dd : {0, 1, 2})
500 mwlsApprox->approxPointCoords[dd] = t_approx_base_point(dd);
501
502 mwlsApprox->getInfluenceNodes() = influence_nodes_map[gg];
503 mwlsApprox->getInvAB() = inv_AB_map[gg];
504 mwlsApprox->getShortenDm() = dm_nodes_map[gg];
505
506 // approximate at MESH_NODE_POSITIONS
507 CHKERR mwlsApprox->evaluateApproxFun(&t_pos(0));
508 CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_pos(0), false);
509 CHKERR mwlsApprox->getTagData(th);
510 rho = mwlsApprox->getDataApprox()[0];
511 const auto &diff_vals = mwlsApprox->getDiffDataApprox();
512 for (int ii = 0; ii != 3; ++ii) {
513 commonData.diffRho(ii, dd) = diff_vals(ii, 0);
514 }
515 t_coeff(i) = pow(rho, betaGC);
516 }
517
518 ++rho;
519 ++t_coords;
520 ++t_coeff;
521 }
522 }
523 // multiply griffith force with coefficients
525 element_prod(commonData.griffithForce, gc_coeff);
527 }
528 };
529
531
533
534 Tag tH;
537
539 Tag *th_coords = nullptr)
540 : mField(m_field), tH(th), tRis(tris), thCoordsPtr(th_coords) {}
541
543
546
547 Range verts;
548 CHKERR mField.get_moab().get_connectivity(tRis, verts, true);
549 {
550 std::vector<double> data(verts.size() * 3, 0);
551 CHKERR mField.get_moab().tag_set_data(tH, verts, &*data.begin());
552 }
553 MatrixDouble diff_n(3, 2);
554 CHKERR ShapeDiffMBTRI(&*diff_n.data().begin());
555
556 for (Range::iterator tit = tRis.begin(); tit != tRis.end(); ++tit) {
557 int num_nodes;
558 const EntityHandle *conn;
559 CHKERR mField.get_moab().get_connectivity(*tit, conn, num_nodes, true);
560 VectorDouble9 coords(9);
561
562 int nb_dofs = num_nodes * 3;
563 auxFun.referenceCoords.resize(9, false);
564 auxFun.currentCoords.resize(9, false);
565
566 if (thCoordsPtr) {
567 CHKERR mField.get_moab().tag_get_data(*thCoordsPtr, conn, num_nodes,
568 &*coords.begin());
569 } else {
570 CHKERR mField.get_moab().get_coords(conn, num_nodes,
571 &*coords.begin());
572 }
573 for (int dd = 0; dd != nb_dofs; dd++) {
574 auxFun.currentCoords[dd] = coords[dd];
575 }
576 if (thCoordsPtr)
577 CHKERR mField.get_moab().get_coords(conn, num_nodes,
578 &*coords.begin());
579 for (int dd = 0; dd != nb_dofs; dd++) {
580 auxFun.referenceCoords[dd] = coords[dd];
581 }
582
583 CHKERR auxFun.calculateGriffithForce(1, 0.5, diff_n);
584
585 double *tag_vals[3];
586 CHKERR mField.get_moab().tag_get_by_ptr(tH, conn, 3,
587 (const void **)tag_vals);
588 for (int nn = 0; nn != 3; ++nn) {
589 for (int dd = 0; dd != 3; ++dd) {
590 tag_vals[nn][dd] += auxFun.griffithForce[3 * nn + dd];
591 }
592 }
593 }
594
596 }
597 };
598
600 AuxOp {
601
602 OpLhs(int tag, BlockData &block_data, CommonData &common_data,
603 const double beta_gc = 0)
605 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
606 UserDataOperator::OPROWCOL),
607 betaGC(beta_gc), AuxOp(tag, block_data, common_data) {}
608
610 const double betaGC;
611
612 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
613 EntityType col_type,
617
618 if (row_type != MBVERTEX) {
620 }
621
629
631 auxFun.currentCoords.resize(9, false);
632
633 int nb_dofs = row_data.getFieldData().size();
634 for (int dd = 0; dd != nb_dofs; dd++)
635 auxFun.currentCoords[dd] = row_data.getFieldData()[dd];
636
637 CHKERR setIndices(row_data);
638 CHKERR setVariables(this, row_data);
639 int nb_gauss_pts = row_data.getN().size1();
640 int nb_base_row = row_data.getFieldData().size() / 3;
641 int nb_base_col = col_data.getFieldData().size() / 3;
642 int row_nb_dofs = row_data.getIndices().size();
643
644 mat.resize(9, 9, false);
645 mat.clear();
646 auto calculate_derivative = [&](FTensor::Tensor1<double *, 2> &t_diff) {
648 t_d_n(i, j) = FTensor::levi_civita(i, j, k) * auxFun.t_Tangent1_ksi(k) *
649 t_diff(N1);
650 t_d_n(i, j) -= FTensor::levi_civita(i, j, k) *
651 auxFun.t_Tangent2_eta(k) * t_diff(N0);
652
653 return t_d_n;
654 };
655
656 CHKERR auxFun.calculateAreaAndNormal(col_data.getDiffN());
657 auto t_normal = auxFun.t_Normal;
658 const double coefficient_1 = 0.5 * pow((t_normal(i) * t_normal(i)), -0.5);
659 const double coefficient_2 = 0.5 * pow((t_normal(i) * t_normal(i)), -1.5);
660
661 // for homogeneous case, this is a bit ugly
662 auto &rho_v = commonData.densityRho;
663 if (rho_v.empty() || rho_v.size() != nb_base_row) {
664 rho_v.resize(nb_base_row, false);
665 }
666 auto rho = getFTensor0FromVec(rho_v);
667 for (int gg = 0; gg != nb_gauss_pts; gg++) {
668
669 double val = getGaussPts()(2, gg) * 0.5;
670 auto t_row_diff = row_data.getFTensor1DiffN<2>(gg, 0);
671 for (int rrr = 0; rrr != nb_base_row; rrr++) {
672
674 &mat(3 * rrr + 0, 0), &mat(3 * rrr + 0, 1), &mat(3 * rrr + 0, 2),
675 &mat(3 * rrr + 1, 0), &mat(3 * rrr + 1, 1), &mat(3 * rrr + 1, 2),
676 &mat(3 * rrr + 2, 0), &mat(3 * rrr + 2, 1), &mat(3 * rrr + 2, 2),
677 3);
678
679 auto t_tan_row = calculate_derivative(t_row_diff);
680 auto t_col_diff = col_data.getFTensor1DiffN<2>(gg, 0);
681 const double coef = blockData.gc * val * pow(rho, betaGC);
682 for (int ccc = 0; ccc != nb_base_col; ccc++) {
683
684 auto t_tan_col = calculate_derivative(t_col_diff);
685
686 t_mat(i, j) -= coefficient_1 *
687 (t_tan_row(i, l) * t_tan_col(l, j) +
688 FTensor::levi_civita(i, j, k) * t_col_diff(N0) *
689 t_row_diff(N1) * t_normal(k) -
690 FTensor::levi_civita(i, j, k) * t_col_diff(N1) *
691 t_row_diff(N0) * t_normal(k));
692
693 t_mat(i, j) += coefficient_2 * (t_tan_row(i, k) * t_normal(k)) *
694 (t_tan_col(l, j) * t_normal(l));
695
696 t_mat(i, j) *= coef;
697 ++t_col_diff;
698 ++t_mat;
699 }
700 ++t_row_diff;
701 ++rho;
702 }
703 }
704
705 int col_nb_dofs = col_data.getIndices().size();
706
707 CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
708 &*rowIndices.data().begin(), col_nb_dofs,
709 &*col_data.getIndices().data().begin(),
710 &*mat.data().begin(), ADD_VALUES);
711
713 }
714 };
715
718 AuxOp {
719
720 OpHeterogeneousGcLhs(int tag, BlockData &block_data,
721 CommonData &common_data, const double beta_gc)
723 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
724 UserDataOperator::OPROWCOL),
725 betaGC(beta_gc), AuxOp(tag, block_data, common_data) {}
726
728 const double betaGC;
729
730 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
731 EntityType col_type,
735
736 if (row_type != MBVERTEX) {
738 }
739
742
744 auxFun.currentCoords.resize(9, false);
745
746 int nb_dofs = row_data.getFieldData().size();
747 for (int dd = 0; dd != nb_dofs; dd++)
748 auxFun.currentCoords[dd] = row_data.getFieldData()[dd];
749
750 CHKERR setIndices(row_data);
751 CHKERR setVariables(this, row_data);
752 int nb_gauss_pts = row_data.getN().size1();
753 int nb_base_row = row_data.getFieldData().size() / 3;
754 int nb_base_col = col_data.getFieldData().size() / 3;
755
756 int row_nb_dofs = row_data.getIndices().size();
757 mat.resize(9, 9, false);
758 mat.clear();
759
760 const double &gc = blockData.gc;
762 auto get_ftensor_from_vec = [](auto &v) {
763 return FTensor::Tensor1<double *, 3>(&v(0), &v(1), &v(2), 3);
764 };
765
766 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
768 &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
769 &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
770 &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
771 &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
772 &m(3 * r + 2, 3 * c + 2));
773 };
774
775 for (int gg = 0; gg != nb_gauss_pts; gg++) {
776
777 double val = getGaussPts()(2, gg) * 0.5;
778
779 auxFun.calculateGriffithForce(1, 1, row_data.getDiffN(gg));
780
781 auto t_griffith = get_ftensor_from_vec(auxFun.griffithForce);
782 auto t_diff_rho = getFTensor1FromMat<3>(commonData.diffRho);
783 for (int rrr = 0; rrr != nb_base_row; rrr++) {
784
785 const double a = val * gc * pow(rho, betaGC - 1.) * betaGC;
787 k(i, j) = t_diff_rho(j) * t_griffith(i) * a;
788 // there is no dependence on the other nodes so row == col
789 auto tt_mat = get_tensor2(mat, rrr, rrr);
790 tt_mat(i, j) += k(i, j);
791
792 ++t_diff_rho;
793 ++t_griffith;
794 ++rho;
795 }
796 }
797
798 int col_nb_dofs = col_data.getIndices().size();
799
800 CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
801 &*rowIndices.data().begin(), col_nb_dofs,
802 &*col_data.getIndices().data().begin(),
803 &*mat.data().begin(), ADD_VALUES);
804
806 }
807 };
808
810
811 const std::string lambdaFieldName;
813 SmartPetscObj<Vec> deltaVec;
814
816 const std::string &lambda_field_name)
817 : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name) {}
818
821 Vec l;
822 CHKERR VecGhostGetLocalForm(deltaVec, &l);
823 double *a;
824 CHKERR VecGetArray(l, &a);
825 int n;
826 CHKERR VecGetLocalSize(l, &n);
827 std::fill(&a[0], &a[n], 0);
828 CHKERR VecRestoreArray(l, &a);
829 CHKERR VecGhostRestoreLocalForm(deltaVec, &l);
831 }
832
835 MOFEM_LOG_CHANNEL("WORLD");
837 MOFEM_LOG_TAG("WORLD", "Arc");
838
839 CHKERR VecAssemblyBegin(deltaVec);
840 CHKERR VecAssemblyEnd(deltaVec);
841 CHKERR VecGhostUpdateBegin(deltaVec, ADD_VALUES, SCATTER_REVERSE);
842 CHKERR VecGhostUpdateEnd(deltaVec, ADD_VALUES, SCATTER_REVERSE);
843 CHKERR VecGhostUpdateBegin(deltaVec, INSERT_VALUES, SCATTER_FORWARD);
844 CHKERR VecGhostUpdateEnd(deltaVec, INSERT_VALUES, SCATTER_FORWARD);
845 double nrm;
846 CHKERR VecNorm(deltaVec, NORM_2, &nrm);
847 MOFEM_LOG_C("WORLD", Sev::noisy, "\tdelta vec nrm = %9.8e", nrm);
849 }
850 };
851
854 AuxOp {
855
857 const std::string lambdaFieldName;
858 SmartPetscObj<Vec> &deltaVec;
859
860 OpConstrainsDelta(MoFEM::Interface &m_field, int tag, BlockData &block_data,
861 CommonData &common_data,
862 const std::string &lambda_field_name,
863 SmartPetscObj<Vec> &delta_vec)
865 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
866 AuxOp(tag, block_data, common_data), mField(m_field),
867 lambdaFieldName(lambda_field_name), deltaVec(delta_vec) {}
868
870
874 if (type != MBVERTEX) {
876 }
877 // get indices of material DOFs
878 CHKERR setIndices(data);
879 // get indices of Lagrange multiplier DOFs
881 // do calculations
882 auxFun.currentCoords.resize(9, false);
883 for (int dd = 0; dd != 9; ++dd)
884 auxFun.currentCoords[dd] = data.getFieldData()[dd];
885
886 CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
887 data.getDiffN(0));
888
889 // assemble area change
890 double delta[] = {0, 0, 0};
891 for (int nn = 0; nn != 3; ++nn) {
892 if (rowLambdaIndices[nn] == -1) {
893 if (rowIndices[3 * nn + 0] != -1 || rowIndices[3 * nn + 1] != -1 ||
894 rowIndices[3 * nn + 2] != -1) {
895 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
896 "Element has no nodes at crack front");
897 }
898 continue;
899 }
900 for (int dd = 0; dd != 3; ++dd) {
901 int idx = 3 * nn + dd;
902 delta[nn] += auxFun.griffithForce[idx] *
903 (data.getFieldData()[idx] - getCoords()[idx]);
904 }
905 }
906 // add nodes to vector of lambda indices
907 CHKERR VecSetValues(deltaVec, 3, &rowLambdaIndices[0], delta, ADD_VALUES);
908
910 }
911 };
912
914
915 const std::string lambdaFieldName;
917 SmartPetscObj<Vec> deltaVec;
918
920 const std::string &lambda_field_name,
921 BlockData &block_data, SmartPetscObj<Vec> &delta_vec)
922 : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name),
923 blockData(block_data), deltaVec(delta_vec) {}
924
927
928 MOFEM_LOG_CHANNEL("WORLD");
930 MOFEM_LOG_TAG("WORLD", "GriffithConstrain");
931
932 const double *delta;
933 CHKERR VecGetArrayRead(deltaVec, &delta);
934 double sum_of_delta = 0;
935 double sum_of_lambda = 0;
936 const auto bit_field_number =
938 switch (snes_ctx) {
939 case CTX_SNESSETFUNCTION: {
941 bit_field_number, dit)) {
942 if (static_cast<int>(dit->get()->getPart()) ==
944 int local_idx = dit->get()->getPetscLocalDofIdx();
945 double lambda = dit->get()->getFieldData();
946 double W = lambda - blockData.E * delta[local_idx];
947 double val = lambda - calMax(W, 0, blockData.r);
948 int global_idx = dit->get()->getPetscGlobalDofIdx();
949 CHKERR VecSetValue(snes_f, global_idx, val, ADD_VALUES);
950 }
951 }
952 } break;
953 case CTX_SNESSETJACOBIAN: {
955 bit_field_number, dit)) {
956 if (static_cast<int>(dit->get()->getPart()) ==
958 int local_idx = dit->get()->getPetscLocalDofIdx();
959 double lambda = dit->get()->getFieldData();
960 double W = lambda - blockData.E * delta[local_idx];
961 double diffW = 1;
962 double val = 1 - diffCalMax_a(W, 0, blockData.r) * diffW;
963 int global_idx = dit->get()->getPetscGlobalDofIdx();
964 MOFEM_LOG_C("WORLD", Sev::verbose,
965 "Constrains on node %lu diag = %+3.5e "
966 "delta = %+3.5e lambda = %+3.5e",
967 dit->get()->getEnt(), val, delta[local_idx], lambda);
968 CHKERR MatSetValue(snes_B, global_idx, global_idx, val, ADD_VALUES);
969 sum_of_delta += delta[local_idx];
970 sum_of_lambda += lambda;
971 }
972 }
973 MOFEM_LOG_C("WORLD", Sev::noisy,
974 "Sum of delta = %+6.4e Sum of lambda = %+6.4e",
975 sum_of_delta, sum_of_lambda);
976
977 } break;
978 default:
979 break;
980 }
981 CHKERR VecRestoreArrayRead(deltaVec, &delta);
983 }
984
988 }
989 };
990
993 AuxOp {
994
996 const std::string lambdaFieldName;
997
998 OpConstrainsRhs(MoFEM::Interface &m_field, int tag, BlockData &block_data,
999 CommonData &common_data,
1000 const std::string &lambda_field_name)
1002 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
1003 AuxOp(tag, block_data, common_data), mField(m_field),
1004 lambdaFieldName(lambda_field_name) {}
1005
1008 // VectorDouble3 nG;
1009
1013 if (type != MBVERTEX)
1015
1016 // get indices of material DOFs
1017 CHKERR setIndices(data);
1020
1021 auxFun.currentCoords.resize(9, false);
1022 noalias(auxFun.currentCoords) = data.getFieldData();
1023
1024 CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1025 data.getDiffN(0));
1026
1027 // assemble local vector
1028 nF.resize(9, false);
1029 nF.clear();
1030 for (int nn = 0; nn != 3; nn++) {
1031 if (rowLambdaIndices[nn] != -1) {
1032 const auto lambda = lambdaAtNodes[nn];
1033 for (int dd = 0; dd != 3; dd++) {
1034 int idx = 3 * nn + dd;
1035 nF[idx] -= lambda * auxFun.griffithForce[idx];
1036 }
1037 }
1038 }
1039
1040 // assemble the right hand vectors
1041 CHKERR VecSetValues(getFEMethod()->snes_f, 9, &*rowIndices.data().begin(),
1042 &nF[0], ADD_VALUES);
1043
1045 }
1046 };
1047
1050 AuxOp {
1051
1053 SmartPetscObj<Vec> &deltaVec;
1054
1055 OpConstrainsLhs(MoFEM::Interface &m_field, int tag, BlockData &block_data,
1056 CommonData &common_data,
1057 const std::string &lambda_field_name,
1058 SmartPetscObj<Vec> &delta_vec)
1060 "MESH_NODE_POSITIONS", lambda_field_name,
1061 UserDataOperator::OPROWCOL,
1062 false // not symmetric operator
1063 ),
1064 AuxOp(tag, block_data, common_data), mField(m_field),
1065 deltaVec(delta_vec) {}
1066
1070
1072 ublas::vector<double *> jacDeltaPtr;
1075
1076 ublas::vector<adouble> a_Delta;
1078 ublas::vector<adouble> a_dA;
1079 ublas::vector<adouble> a_lambdaAtNodes;
1080
1081 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1082 EntityType col_type,
1086
1087 if (row_type != MBVERTEX || col_type != MBVERTEX)
1089
1090 CHKERR setIndices(row_data);
1091 CHKERR setLambdaNodes(this, colFieldName);
1092 CHKERR setLambdaIndices(this, colFieldName);
1093
1094 activeVariables.resize(18, false);
1095 a_auxFun.referenceCoords.resize(9, false);
1096 a_auxFun.currentCoords.resize(9, false);
1097
1098 a_Delta.resize(3, false);
1099 d_Delta.resize(3, false);
1100 trace_on(tAg);
1101 for (int dd = 0; dd != 9; dd++) {
1102 double val = getCoords()[dd];
1103 a_auxFun.referenceCoords[dd] <<= val;
1104 activeVariables[dd] = val;
1105 }
1106 for (int dd = 0; dd != 9; dd++) {
1107 double val = row_data.getFieldData()[dd];
1108 a_auxFun.currentCoords[dd] <<= val;
1109 activeVariables[9 + dd] = val;
1110 }
1111
1112 CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1113 row_data.getDiffN(0));
1114
1115 // calculate area change
1116 a_Delta.clear();
1117 for (int nn = 0; nn != 3; nn++) {
1118 if (col_data.getIndices()[nn] == -1)
1119 continue;
1120 if (col_data.getIndices()[nn] != rowLambdaIndices[nn]) {
1121 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1122 "data inconsistency");
1123 }
1124 for (int dd = 0; dd != 3; dd++) {
1125 int idx = 3 * nn + dd;
1126 if (row_data.getFieldDofs()[idx]->getEnt() !=
1127 col_data.getFieldDofs()[nn]->getEnt()) {
1128 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1129 "data inconsistency");
1130 }
1131 a_Delta[nn] +=
1132 a_auxFun.griffithForce[idx] *
1133 (a_auxFun.currentCoords[idx] - a_auxFun.referenceCoords[idx]);
1134 }
1135 }
1136 for (int nn = 0; nn != 3; nn++) {
1137 a_Delta[nn] >>= d_Delta[nn];
1138 }
1139 trace_off();
1140
1141 jacDelta.resize(3, 18, false);
1142 jacDeltaPtr.resize(3, false);
1143 for (int nn = 0; nn != 3; nn++) {
1144 jacDeltaPtr[nn] = &jacDelta(nn, 0);
1145 }
1146
1147 {
1148 int r;
1149 // play recorder for jacobian
1150 r = ::jacobian(tAg, 3, 18, &activeVariables[0], &jacDeltaPtr[0]);
1151 if (r < 3) { // function is locally analytic
1152 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1153 "ADOL-C function evaluation with error r = %d", r);
1154 }
1155 }
1156
1157 const double *delta;
1158 Vec delta_vec_local;
1159 CHKERR VecGhostGetLocalForm(deltaVec, &delta_vec_local);
1160 int vec_size;
1161 CHKERR VecGetSize(delta_vec_local, &vec_size);
1162 CHKERR VecGetArrayRead(delta_vec_local, &delta);
1163 nG_dX.resize(3, 9, false);
1164 nG_dX.clear();
1165 for (int nn = 0; nn != 3; ++nn) {
1166 int local_idx = col_data.getLocalIndices()[nn];
1167 if (local_idx != rowLambdaIndicesLocal[nn]) {
1168 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1169 "data inconsistency");
1170 }
1171 if (local_idx == -1)
1172 continue;
1173 if (vec_size < local_idx) {
1174 int g_vec_size;
1175 CHKERR VecGetLocalSize(deltaVec, &g_vec_size);
1176 SETERRQ3(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1177 "data inconsistency %d < %d (%d)", vec_size, local_idx,
1178 g_vec_size);
1179 }
1180 double W = col_data.getFieldData()[nn] - blockData.E * delta[local_idx];
1181 double diff_constrain = -diffCalMax_a(W, 0, blockData.r);
1182 for (int dd = 0; dd != 9; dd++) {
1183 double diffW = -blockData.E * jacDelta(nn, 9 + dd);
1184 nG_dX(nn, dd) = diff_constrain * diffW;
1185 }
1186 }
1187 CHKERR VecRestoreArrayRead(delta_vec_local, &delta);
1188 CHKERR VecGhostRestoreLocalForm(deltaVec, &delta_vec_local);
1189 CHKERR MatSetValues(getFEMethod()->snes_B, 3,
1190 &*col_data.getIndices().data().begin(), 9,
1191 &*row_data.getIndices().data().begin(),
1192 &*nG_dX.data().begin(), ADD_VALUES);
1193
1194 // Dervatives of nF
1195
1196 // Linearisation with adol-C
1197 activeVariables.resize(21, false);
1198 d_dA.resize(9, false);
1199 a_dA.resize(9, false);
1200 a_lambdaAtNodes.resize(3, false);
1201 a_auxFun.referenceCoords.resize(9, false);
1202
1203 trace_on(tAg);
1204 for (int dd = 0; dd != 9; dd++) {
1205 double val = getCoords()[dd];
1206 a_auxFun.referenceCoords[dd] <<= val;
1207 activeVariables[dd] = val;
1208 }
1209 for (int dd = 0; dd != 9; dd++) {
1210 double val = row_data.getFieldData()[dd];
1211 a_auxFun.currentCoords[dd] <<= val;
1212 activeVariables[9 + dd] = val;
1213 }
1214 for (int nn = 0; nn != 3; nn++) {
1215 double val = col_data.getFieldData()[nn];
1216 a_lambdaAtNodes[nn] <<= val;
1217 activeVariables[18 + nn] = val;
1218 }
1219
1220 CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1221 row_data.getDiffN(0));
1222
1223 a_dA.clear();
1224 for (int nn = 0; nn != 3; nn++) {
1225 if (col_data.getIndices()[nn] == -1)
1226 continue;
1227 for (int dd = 0; dd != 3; dd++) {
1228 int idx = 3 * nn + dd;
1229 a_dA[idx] -= a_lambdaAtNodes[nn] * a_auxFun.griffithForce[idx];
1230 }
1231 }
1232 for (int dd = 0; dd != 9; dd++) {
1233 a_dA[dd] >>= d_dA[dd];
1234 }
1235 trace_off();
1236
1237 commonData.tangentGriffithForce.resize(9, 21, false);
1239 for (int dd = 0; dd != 9; ++dd) {
1242 }
1243
1244 {
1245 int r;
1246 // play recorder for jacobian
1247 r = ::jacobian(tAg, 9, 21, &activeVariables[0],
1249 if (r < 3) { // function is locally analytic
1250 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1251 "ADOL-C function evaluation with error r = %d", r);
1252 }
1253 }
1254
1255 nF_dLambda.resize(9, 3);
1256 for (int rr = 0; rr != 9; rr++) {
1257 for (int cc = 0; cc != 3; cc++) {
1258 nF_dLambda(rr, cc) = commonData.tangentGriffithForce(rr, 18 + cc);
1259 }
1260 }
1261 CHKERR MatSetValues(getFEMethod()->snes_B, 9, &*rowIndices.data().begin(),
1262 3, &*col_data.getIndices().data().begin(),
1263 &*nF_dLambda.data().begin(), ADD_VALUES);
1264
1265 nF_dX.resize(9, 9, false);
1266 for (int rr = 0; rr != 9; rr++) {
1267 for (int cc = 0; cc != 9; cc++) {
1268 nF_dX(rr, cc) = commonData.tangentGriffithForce(rr, 9 + cc);
1269 }
1270 }
1271 CHKERR MatSetValues(getFEMethod()->snes_B, 9, &*rowIndices.data().begin(),
1272 9, &*row_data.getIndices().data().begin(),
1273 &*nF_dX.data().begin(), ADD_VALUES);
1274
1276 }
1277 };
1278
1279 /**
1280 * \brief arc-length element
1281 *
1282 * Calualte residual for arc length method amd vector db, which is derivative
1283 * of residual.
1284 *
1285 */
1287
1288 int tAg;
1289 boost::shared_ptr<ArcLengthCtx> arcPtr;
1292 std::string lambdaFieldName;
1293
1299
1300 boost::shared_ptr<MyTriangleFE> feIntLambda;
1301
1302 /**
1303 * \brief assemble internal residual
1304 * @param block_data griffith block data
1305 * @param common_data griffith common data
1306 * @param lambda_field_name lagrange multipliers controling crack surface
1307 * area
1308 * @param arc_front reference to FrontArcLengthControl object
1309 */
1312
1314 std::string lambdaFieldName;
1316
1320 const std::string &lambda_field_name,
1321 FrontArcLengthControl &arc_front)
1323 "MESH_NODE_POSITIONS", lambda_field_name,
1324 UserDataOperator::OPROW,
1325 false // not symmetric operator
1326 ),
1327 AuxOp(0, block_data, common_data), mField(m_field),
1328 lambdaFieldName(lambda_field_name), arcFront(arc_front) {}
1329
1331
1335 if (type != MBVERTEX) {
1337 }
1338 // get indices of Lagrange multiplier DOFs
1340 // get dofs increment
1341 const double *dx;
1342 CHKERR VecGetArrayRead(arcFront.arcPtr->dx, &dx);
1343 // get dofs at begining of load step
1344 const double *x0;
1345 CHKERR VecGetArrayRead(arcFront.arcPtr->x0, &x0);
1346 auxFun.referenceCoords.resize(9, false);
1347 auxFun.currentCoords.resize(9, false);
1348 for (int dd = 0; dd != 9; dd++) {
1349 int loc_idx = data.getLocalIndices()[dd];
1350 if (loc_idx != -1) {
1351 auxFun.referenceCoords[dd] = x0[loc_idx];
1352 auxFun.currentCoords[dd] = auxFun.referenceCoords[dd] + dx[loc_idx];
1353 } else {
1354 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
1355 }
1356 }
1357 CHKERR VecRestoreArrayRead(arcFront.arcPtr->dx, &dx);
1358 CHKERR VecRestoreArrayRead(arcFront.arcPtr->x0, &x0);
1359
1360 // calculate crack area increment
1361 CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1362 data.getDiffN(0));
1363
1364 // calculate area change
1365 double delta = 0;
1366 for (int nn = 0; nn != 3; nn++) {
1367 if (rowLambdaIndices[nn] == -1)
1368 continue;
1369 auto griffith_force =
1370 getVectorAdaptor(&(auxFun.griffithForce[3 * nn]), 3);
1371 // griffith_force /= norm_2(griffith_force);
1372 for (int dd = 0; dd != 3; dd++) {
1373 const int idx = 3 * nn + dd;
1374 delta += griffith_force[dd] *
1376 }
1377 }
1378 if (delta != delta) {
1379 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "nan value");
1380 }
1383 };
1384 };
1385
1386 /**
1387 * \brief assemble internal residual derivative
1388 * @param block_data griffith block data
1389 * @param common_data griffith common data
1390 * @param lambda_field_name lagrange multipliers controling crack surface
1391 * area
1392 * @param arc_front reference to FrontArcLengthControl object
1393 */
1394 struct OpAssemble_db : public OpIntLambda {
1395 int tAg;
1399 const std::string &lambda_field_name,
1400 FrontArcLengthControl &arc_front)
1401 : OpIntLambda(m_field, block_data, common_data, lambda_field_name,
1402 arc_front),
1403 tAg(tag) {}
1404
1408
1412 if (type != MBVERTEX) {
1414 }
1415 CHKERR setIndices(data);
1416 // get indices of Lagrange multiplier DOFs
1418 double d_delta;
1419 adouble a_delta;
1420 // get dofs increment
1421 const double *dx;
1422 CHKERR VecGetArrayRead(arcFront.arcPtr->dx, &dx);
1423 // get dofs at begining of load step
1424 const double *x0;
1425 CHKERR VecGetArrayRead(arcFront.arcPtr->x0, &x0);
1426 a_auxFun.referenceCoords.resize(9, false);
1427 a_auxFun.currentCoords.resize(9, false);
1428
1429 activeVariables.resize(18, false);
1430 trace_on(tAg);
1431 for (int dd = 0; dd != 9; dd++) {
1432 int loc_idx = data.getLocalIndices()[dd];
1433 if (loc_idx != -1) {
1434 a_auxFun.referenceCoords[dd] <<= x0[loc_idx];
1435 activeVariables[dd] = x0[loc_idx];
1436 } else {
1437 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
1438 }
1439 }
1440 for (int dd = 0; dd != 9; dd++) {
1441 int loc_idx = data.getLocalIndices()[dd];
1442 if (loc_idx != -1) {
1443 a_auxFun.currentCoords[dd] <<= x0[dd] + dx[loc_idx];
1444 activeVariables[9 + dd] = x0[loc_idx] + dx[loc_idx];
1445 } else {
1446 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
1447 }
1448 }
1449 // calculate crack area increment
1450
1451 CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
1452 data.getDiffN(0));
1453 // calculate area change
1454 a_delta = 0;
1455 for (int nn = 0; nn != 3; nn++) {
1456 if (rowLambdaIndices[nn] == -1)
1457 continue;
1458 auto griffith_force =
1459 getVectorAdaptor(&(a_auxFun.griffithForce[3 * nn]), 3);
1460 // griffith_force /= sqrt(griffith_force[0] * griffith_force[0] +
1461 // griffith_force[1] * griffith_force[1] +
1462 // griffith_force[2] * griffith_force[2]);
1463 for (int dd = 0; dd != 3; dd++) {
1464 const int idx = 3 * nn + dd;
1465 a_delta += griffith_force[dd] * (a_auxFun.currentCoords[idx] -
1466 a_auxFun.referenceCoords[idx]);
1467 }
1468 }
1469 a_delta >>= d_delta;
1470 trace_off();
1471 CHKERR VecRestoreArrayRead(arcFront.arcPtr->dx, &dx);
1472 CHKERR VecRestoreArrayRead(arcFront.arcPtr->x0, &x0);
1473
1474 gradDelta.resize(18, false);
1475 {
1476 int r;
1477 // play recorder for jacobian
1478 r = ::gradient(tAg, 18, &activeVariables[0], &gradDelta[0]);
1479 if (r < 3) { // function is locally analytic
1480 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1481 "ADOL-C function evaluation with error r = %d", r);
1482 }
1483 }
1484
1485 for (int dd = 0; dd != 9; dd++) {
1486 double val = arcFront.arcPtr->alpha * gradDelta[9 + dd];
1487 CHKERR VecSetValue(arcFront.arcPtr->db, data.getIndices()[dd], val,
1488 ADD_VALUES);
1489 }
1490
1492 }
1493 };
1494
1497 const std::string &lambda_field_name,
1498 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
1499 : tAg(tag), arcPtr(arc_ptr), blockData(block_data),
1500 commonData(common_data), lambdaFieldName(lambda_field_name),
1503 int two[] = {0, 1};
1504 int nb_local = arcPtr->mField.get_comm_rank() == 0 ? 2 : 0;
1505 int nb_ghost = nb_local ? 0 : 2;
1506 ierr = VecCreateGhostWithArray(PETSC_COMM_WORLD, nb_local, 2, nb_ghost,
1508 CHKERRABORT(PETSC_COMM_WORLD, ierr);
1509 feIntLambda = boost::shared_ptr<MyTriangleFE>(
1511 }
1513 ierr = VecDestroy(&ghostIntLambda);
1514 CHKERRABORT(PETSC_COMM_WORLD, ierr);
1515 }
1516
1519 switch (snes_ctx) {
1520 case CTX_SNESSETFUNCTION: {
1521 CHKERR VecAssemblyBegin(snes_f);
1522 CHKERR VecAssemblyEnd(snes_f);
1525 } break;
1526 default:
1527 break;
1528 }
1530 }
1531
1532 /** \brief Calculate internal lambda
1533 */
1536 return intLambda;
1538 }
1539
1542 MOFEM_LOG_CHANNEL("WORLD");
1544 MOFEM_LOG_TAG("WORLD", "Arc");
1545
1546 switch (snes_ctx) {
1547 case CTX_SNESSETFUNCTION: {
1548 arcPtr->res_lambda = calculateLambdaInt() - arcPtr->alpha * arcPtr->s;
1549 CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
1550 arcPtr->res_lambda, ADD_VALUES);
1551 } break;
1552 case CTX_SNESSETJACOBIAN: {
1553 arcPtr->dIag = 2 * arcPtr->dLambda * arcPtr->alpha *
1554 pow(arcPtr->beta, 2) * arcPtr->F_lambda2;
1555 MOFEM_LOG_C("WORLD", Sev::noisy, "diag = %9.8e", arcPtr->dIag);
1556 CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
1557 arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
1558 } break;
1559 default:
1560 break;
1561 }
1563 }
1564
1567 switch (snes_ctx) {
1568 case CTX_SNESSETFUNCTION: {
1569 CHKERR VecAssemblyBegin(snes_f);
1570 CHKERR VecAssemblyEnd(snes_f);
1571 } break;
1572 case CTX_SNESSETJACOBIAN: {
1573 CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
1574 SCATTER_FORWARD);
1575 CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES,
1576 SCATTER_FORWARD);
1577 } break;
1578 default:
1579 break;
1580 }
1582 }
1583
1584 /** \brief Calculate db
1585 */
1588
1589 MOFEM_LOG_CHANNEL("WORLD");
1591 MOFEM_LOG_TAG("WORLD", "Arc");
1592
1593 double alpha = arcPtr->alpha;
1594 double beta = arcPtr->beta;
1595 double d_lambda = arcPtr->dLambda;
1596
1597 CHKERR VecZeroEntries(ghostIntLambda);
1598 CHKERR VecGhostUpdateBegin(ghostIntLambda, INSERT_VALUES,
1599 SCATTER_FORWARD);
1600 CHKERR VecGhostUpdateEnd(ghostIntLambda, INSERT_VALUES, SCATTER_FORWARD);
1601 CHKERR VecZeroEntries(arcPtr->db);
1602 CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1603 CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1604
1605 // calculate internal lambda
1606 feIntLambda->getOpPtrVector().clear();
1607 feIntLambda->getOpPtrVector().push_back(new OpIntLambda(
1608 arcPtr->mField, blockData, commonData, lambdaFieldName, *this));
1609 feIntLambda->getOpPtrVector().push_back(new OpAssemble_db(
1610 tAg, arcPtr->mField, blockData, commonData, lambdaFieldName, *this));
1611 CHKERR arcPtr->mField.loop_finite_elements(
1612 problemPtr->getName(), "CRACKFRONT_AREA_ELEM", (*feIntLambda));
1613 const double *dx;
1614 const auto bit_field_number = getFieldBitNumber(lambdaFieldName);
1615 CHKERR VecGetArrayRead(arcPtr->dx, &dx);
1617 bit_field_number, dit)) {
1618 if (static_cast<int>(dit->get()->getPart()) !=
1619 arcPtr->mField.get_comm_rank())
1620 continue;
1621 int idx = dit->get()->getPetscLocalDofIdx();
1622 double val = dx[idx];
1623 CHKERR VecSetValue(ghostIntLambda, 1, val, ADD_VALUES);
1624 }
1625 CHKERR VecRestoreArrayRead(arcPtr->dx, &dx);
1626
1627 CHKERR VecAssemblyBegin(ghostIntLambda);
1628 CHKERR VecAssemblyEnd(ghostIntLambda);
1629 CHKERR VecGhostUpdateBegin(ghostIntLambda, ADD_VALUES, SCATTER_REVERSE);
1630 CHKERR VecGhostUpdateEnd(ghostIntLambda, ADD_VALUES, SCATTER_REVERSE);
1631 CHKERR VecGhostUpdateBegin(ghostIntLambda, INSERT_VALUES,
1632 SCATTER_FORWARD);
1633 CHKERR VecGhostUpdateEnd(ghostIntLambda, INSERT_VALUES, SCATTER_FORWARD);
1634
1635 CHKERR VecAssemblyBegin(arcPtr->db);
1636 CHKERR VecAssemblyEnd(arcPtr->db);
1637 CHKERR VecGhostUpdateBegin(arcPtr->db, ADD_VALUES, SCATTER_REVERSE);
1638 CHKERR VecGhostUpdateEnd(arcPtr->db, ADD_VALUES, SCATTER_REVERSE);
1639 CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1640 CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
1641
1643 alpha * beta * beta * d_lambda * d_lambda * arcPtr->F_lambda2;
1644
1646 "WORLD", Sev::noisy,
1647 "\tintLambdaLambda = %9.8e intLambdaDelta = %9.8e intLambda = %9.8e",
1649
1650 double nrm;
1651 CHKERR VecNorm(arcPtr->db, NORM_2, &nrm);
1652 MOFEM_LOG_C("WORLD", Sev::noisy, "\tdb nrm = %9.8e", nrm);
1653
1655 }
1656
1659 MOFEM_LOG_CHANNEL("WORLD");
1661 MOFEM_LOG_TAG("WORLD", "Arc");
1662
1663 // Calculate dx
1664 CHKERR VecGhostUpdateBegin(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
1665 CHKERR VecGhostUpdateEnd(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
1666
1667 Vec l_x, l_x0, l_dx;
1668 CHKERR VecGhostGetLocalForm(x, &l_x);
1669 CHKERR VecGhostGetLocalForm(arcPtr->x0, &l_x0);
1670 CHKERR VecGhostGetLocalForm(arcPtr->dx, &l_dx);
1671 {
1672 double *x_array, *x0_array, *dx_array;
1673 CHKERR VecGetArray(l_x, &x_array);
1674 CHKERR VecGetArray(l_x0, &x0_array);
1675 CHKERR VecGetArray(l_dx, &dx_array);
1676 int size =
1677 problemPtr->getNbLocalDofsRow() + problemPtr->getNbGhostDofsRow();
1678 for (int i = 0; i != size; ++i) {
1679 dx_array[i] = x_array[i] - x0_array[i];
1680 }
1681 CHKERR VecRestoreArray(l_x, &x_array);
1682 CHKERR VecRestoreArray(l_x0, &x0_array);
1683 CHKERR VecRestoreArray(l_dx, &dx_array);
1684 }
1685 CHKERR VecGhostRestoreLocalForm(x, &l_x);
1686 CHKERR VecGhostRestoreLocalForm(arcPtr->x0, &l_x0);
1687 CHKERR VecGhostRestoreLocalForm(arcPtr->dx, &l_dx);
1688
1689 // Calculate dlambda
1690 if (arcPtr->getPetscLocalDofIdx() != -1) {
1691 double *array;
1692 CHKERR VecGetArray(arcPtr->dx, &array);
1693 arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
1694 array[arcPtr->getPetscLocalDofIdx()] = 0;
1695 CHKERR VecRestoreArray(arcPtr->dx, &array);
1696 }
1697 CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
1698 SCATTER_FORWARD);
1699 CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
1700 SCATTER_FORWARD);
1701
1702 // Calculate dx2
1703 double x_nrm, x0_nrm;
1704 CHKERR VecNorm(x, NORM_2, &x_nrm);
1705 CHKERR VecNorm(arcPtr->x0, NORM_2, &x0_nrm);
1706 CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
1707 MOFEM_LOG_C("WORLD", Sev::noisy,
1708 "\tx norm = %6.4e x0 norm = %6.4e dx2 = %6.4e", x_nrm, x0_nrm,
1709 arcPtr->dx2);
1711 }
1712 };
1713};
1714
1715} // namespace FractureMechanics
1716
1717#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:304
constexpr double a
EntitiesFieldData::EntData EntData
#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
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
constexpr double lambda
constexpr double W
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:318
#define _IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(PROBLEMPTR, FIELD_BIT_NUMBER, IT)
FTensor::Index< 'i', SPACE_DIM > i
@ TYPE
Definition: inflate.h:32
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
const FTensor::Tensor2< T, Dim, Dim > Vec
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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static double diffCalMax_a(double a, double b, double r)
static double calMax(double a, double b, double r)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
VectorBoundedArray< int, 3 > VectorInt3
Definition: Types.hpp:87
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:96
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
const double r
rate factor
double rho
Definition: plastic.cpp:98
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.