v0.13.1
CrackFrontElement.hpp
Go to the documentation of this file.
1/** \file CrackFrontElement.hpp
2 */
3
4/* This file is part of MoFEM.
5 * MoFEM is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU Lesser General Public License as published by the
7 * Free Software Foundation, either version 3 of the License, or (at your
8 * option) any later version.
9 *
10 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13 * License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
17
18#ifndef __CRACK_FRONT_ELEMENT_HPP__
19#define __CRACK_FRONT_ELEMENT_HPP__
20
21namespace FractureMechanics {
22
24
25template <typename E, typename B> class TwoType {};
26
27/**
28 *
29 * Use idea from \cite barsoum1976use to approximate singularity at crack tip.
30 * The idea is to shift mid noted to quarter edge length to crack tip on edges
31 * adjacent to crack front.
32 *
33 */
34template <typename ELEMENT, typename BASE_ELEMENT>
36
40
46
50
52
54 const bool &set_singular_coordinates,
55 const Range &crack_front_nodes,
56 const Range &crack_front_nodes_edges,
57 const Range &crack_front_elements,
58 PetscBool add_singularity)
59 : ELEMENT(m_field), crackFrontNodes(crack_front_nodes),
60 crackFrontNodesEdges(crack_front_nodes_edges),
61 crackFrontElements(crack_front_elements),
62 addSingularity(add_singularity),
63 setSingularCoordinatesPriv(set_singular_coordinates) {
65 CHKERRABORT(PETSC_COMM_WORLD, ierr);
66 }
67
71
72 virtual ~CrackFrontSingularBase() = default;
73
76
77 if (this->numeredEntFiniteElementPtr->getEntType() != MBTET)
78 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79 "Element type not implemented");
80
81 this->getElementPolynomialBase() =
82 boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
83
84 CHKERR this->createDataOnElement();
85 CHKERR this->calculateVolumeAndJacobian();
87 CHKERR this->setIntegrationPts();
88 CHKERR this->calculateCoordinatesAtGaussPts();
89 CHKERR this->calHierarchicalBaseFunctionsOnElement();
91 CHKERR this->transformBaseFunctions();
92
93 // Iterate over operators
94 CHKERR this->loopOverOperators();
95
97 }
98
99 PetscBool addSingularity;
102
105 }
106
107 template <typename E, typename B>
109 return E::getSpaceBaseAndOrderOnElement();
110 }
111
112 /// Partial specialization for volume elements
113 template <typename E>
117 CHKERR E::getSpaceBaseAndOrderOnElement();
118 // singularElement = false;
119 // MoFEMFunctionReturnHot(0);
120 // Determine if this element is singular
121 EntityHandle ent = this->numeredEntFiniteElementPtr->getEnt();
122 if (crackFrontElements.find(ent) != crackFrontElements.end()) {
123 singularElement = true;
124 for (int nn = 0; nn != 4; nn++) {
125 if (crackFrontNodes.find(this->conn[nn]) != crackFrontNodes.end()) {
126 singularNodes[nn] = true;
127 } else {
128 singularNodes[nn] = false;
129 }
130 }
131 for (int ee = 0; ee != 6; ee++) {
132 EntityHandle edge;
133 CHKERR this->mField.get_moab().side_element(ent, 1, ee, edge);
134 if (crackFrontNodesEdges.find(edge) != crackFrontNodesEdges.end()) {
135 singularEdges[ee] = true;
136 } else {
137 singularEdges[ee] = false;
138 }
139 }
140 } else {
141 singularElement = false;
142 }
144 }
145
149 }
150
151 /// Partial specialization for volume element
152 template <typename E>
156
158 singularElement = false;
160 }
161 if (addSingularity == PETSC_FALSE) {
162 singularElement = false;
164 }
165
166 if (singularElement) {
167
168 const int edge_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
169 {0, 3}, {1, 3}, {2, 3}};
173
174 singularDisp.resize(this->coordsAtGaussPts.size1(),
175 this->coordsAtGaussPts.size2());
176 singularDisp.clear();
177 singularRefCoords.resize(this->coordsAtGaussPts.size1(),
178 this->coordsAtGaussPts.size2());
179 singularRefCoords.clear();
180
181 const size_t nb_gauss_pts = this->gaussPts.size2();
182 sJac.resize(nb_gauss_pts, 9, false);
183 double *s_jac_ptr = &sJac(0, 0);
184 // Set jacobian to 1 on diagonal
185 {
187 s_jac_ptr, &s_jac_ptr[1], &s_jac_ptr[2], &s_jac_ptr[3],
188 &s_jac_ptr[4], &s_jac_ptr[5], &s_jac_ptr[6], &s_jac_ptr[7],
189 &s_jac_ptr[8]);
190 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
191 t_s_jac(i, j) = 0;
192 t_s_jac(0, 0) = 1;
193 t_s_jac(1, 1) = 1;
194 t_s_jac(2, 2) = 1;
195 ++t_s_jac;
196 }
197 }
198
199 // get base fuctions direvatives
200 double diff_base_n[12];
201 CHKERR ShapeDiffMBTET(diff_base_n);
203 diff_base_n, &diff_base_n[1], &diff_base_n[2]);
204 double diff_n[12];
206 diff_n, &diff_n[1], &diff_n[2]);
207 for (int nn = 0; nn != 4; nn++) {
208 t_diff_n(i) = this->tInvJac(j, i) * t_diff_base_n(j);
209 ++t_diff_n;
210 ++t_diff_base_n;
211 }
212 MatrixDouble &shape_n =
213 this->dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
214
215 const double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
216
217 // Calculate singular jacobian
218 for (int ee = 0; ee != 6; ee++) {
219 if (!singularEdges[ee])
220 continue;
221 // i0 is node at crack front, i1 is node inside body
222 const int i0 = edge_nodes[ee][0];
223 const int i1 = edge_nodes[ee][1];
224 // Only one of the nodes at edge have to be at crack front
225 if ((singularNodes[i0] && singularNodes[i1]) ||
226 !(singularNodes[i0] || singularNodes[i1])) {
227 PetscPrintf(PETSC_COMM_SELF, "Warning singular edge on both ends\n");
228 // EntityHandle out_meshset;
229 // CHKERR this->mField.get_moab().create_meshset(
230 // MESHSET_SET, out_meshset);
231 // CHKERR this->mField.get_moab().add_entities(out_meshset,
232 // crackFrontNodesEdges);
233 // CHKERR this->mField.get_moab().write_file(
234 // "debug_error_crack_front_nodes_edges.vtk", "VTK", "",
235 // &out_meshset, 1);
236 // CHKERR this->mField.get_moab().delete_entities(&out_meshset, 1);
237 // CHKERR this->mField.get_moab().create_meshset(
238 // MESHSET_SET, out_meshset);
239 // CHKERR this->mField.get_moab().add_entities(out_meshset,
240 // crackFrontElements);
241 // CHKERR this->mField.get_moab().write_file(
242 // "debug_error_crack_elements.vtk", "VTK", "", &out_meshset, 1);
243 // CHKERR this->mField.get_moab().delete_entities(&out_meshset, 1);
244 continue;
245 // SETERRQ(
246 // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
247 // "Huston we have problem, only one on the edge can be
248 // singular");
249 }
250 // edge directions oriented from crack front node towards inside
252 base_coords[3 * i1 + 0] - base_coords[3 * i0 + 0],
253 base_coords[3 * i1 + 1] - base_coords[3 * i0 + 1],
254 base_coords[3 * i1 + 2] - base_coords[3 * i0 + 2]);
255 // set orientation
256 if (singularNodes[edge_nodes[ee][0]]) {
257 t_base_dir(i) *= -1;
258 }
259 // push direction to current configuration
261 t_dir(i) = this->tJac(i, j) * t_base_dir(j);
263 s_jac_ptr, &s_jac_ptr[1], &s_jac_ptr[2], &s_jac_ptr[3],
264 &s_jac_ptr[4], &s_jac_ptr[5], &s_jac_ptr[6], &s_jac_ptr[7],
265 &s_jac_ptr[8]);
267 t_singular_displacement(&singularDisp(0, 0), &singularDisp(0, 1),
268 &singularDisp(0, 2));
270 t_singular_ref_displacement(&singularRefCoords(0, 0),
271 &singularRefCoords(0, 1),
272 &singularRefCoords(0, 2));
273 // Calculate jacobian
274 for (size_t gg = 0; gg != nb_gauss_pts; gg++) {
275 double t_n = shape_n(gg, i0) * shape_n(gg, i1);
277 shape_n(gg, i0) * diff_n[3 * i1 + 0] +
278 shape_n(gg, i1) * diff_n[3 * i0 + 0],
279 shape_n(gg, i0) * diff_n[3 * i1 + 1] +
280 shape_n(gg, i1) * diff_n[3 * i0 + 1],
281 shape_n(gg, i0) * diff_n[3 * i1 + 2] +
282 shape_n(gg, i1) * diff_n[3 * i0 + 2]);
283 t_s_jac(i, j) += t_dir(i) * t_diff_n(j);
284 t_singular_displacement(i) += t_dir(i) * t_n;
285 t_singular_ref_displacement(i) += t_base_dir(i) * t_n;
286 ++t_singular_displacement;
287 ++t_singular_ref_displacement;
288 ++t_s_jac;
289 }
290 }
291
292 // invert singular jacobian and set integration points
293 {
295 s_jac_ptr, &s_jac_ptr[1], &s_jac_ptr[2], &s_jac_ptr[3],
296 &s_jac_ptr[4], &s_jac_ptr[5], &s_jac_ptr[6], &s_jac_ptr[7],
297 &s_jac_ptr[8], 9);
298 // Allocate memory for singular inverse jacobian and setup tensor.
299 invSJac.resize(nb_gauss_pts, 9, false);
300 double *inv_s_jac_ptr = &invSJac(0, 0);
302 inv_s_jac_ptr, &inv_s_jac_ptr[1], &inv_s_jac_ptr[2],
303 &inv_s_jac_ptr[3], &inv_s_jac_ptr[4], &inv_s_jac_ptr[5],
304 &inv_s_jac_ptr[6], &inv_s_jac_ptr[7], &inv_s_jac_ptr[8], 9);
305 // Sum of all determinist has to be 1
306 // Calculate integrations weights
307 detS.resize(nb_gauss_pts, false);
308 for (size_t gg = 0; gg != nb_gauss_pts; gg++) {
309 double det;
310 CHKERR determinantTensor3by3(t_s_jac, det);
311 CHKERR invertTensor3by3(t_s_jac, det, t_inv_s_jac);
312 detS[gg] = det;
313 ++t_inv_s_jac;
314 ++t_s_jac;
315 }
316 }
317 }
318
320 }
321
322 /**
323 * @return returning negative number -1, i.e. special user quadrature is
324 * generated
325 */
326 int getRule(int order) {
328 }
329
330 // Generic specialization
331 template <typename E, typename B> int getRuleImpl(TwoType<E, B>, int order) {
332 return E::getRule(order);
333 }
334
335 /**
336 * \brief Generate specific user integration quadrature
337 *
338 * Element is refined at the crack top and quadrature is set to each of
339 * sub-tets.
340 *
341 */
344 }
345
346 // Generic specialization
347 template <typename E, typename B>
349 return E::setGaussPts(order);
350 }
351
352private:
354
357};
358
359template <typename ELEMENT, typename BASE_ELEMENT>
363 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Get singular element options",
364 "none");
365 CHKERRQ(ierr);
366 refinementLevels = 3;
367 ierr = PetscOptionsInt("-se_number_of_refinement_levels",
368 "approximation geometry order", "", refinementLevels,
369 &refinementLevels, PETSC_NULL);
370 CHKERRQ(ierr);
371 refineIntegration = PETSC_TRUE;
372 ierr =
373 PetscOptionsBool("-se_refined_integration",
374 "if set element is subdivided to generate quadrature",
375 "", refineIntegration, &refineIntegration, PETSC_NULL);
376 CHKERRQ(ierr);
377 ierr = PetscOptionsEnd();
378 CHKERRQ(ierr);
380}
381
382// specialization
383template <>
384template <>
386 VolumeElementForcesAndSourcesCore>::
388 VolumeElementForcesAndSourcesCore>,
389 int order);
390
391// specialization
392template <>
393template <>
395 VolumeElementForcesAndSourcesCore>::
397 VolumeElementForcesAndSourcesCore>,
398 int order);
399
401 VolumeElementForcesAndSourcesCore>
403
406
407 const Range tRis;
408 boost::ptr_vector<MethodForForceScaling> &methodsOp;
409 boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
411
412 CrackFrontSingularBase<VolumeElementForcesAndSourcesCoreOnSide,
413 VolumeElementForcesAndSourcesCore>
415 // boost::shared_ptr<MatrixDouble> hG; ///< spatial deformation gradient
416 // boost::shared_ptr<MatrixDouble> HG; ///< spatial deformation gradient
417 VectorDouble sNrm; ///< Length of the normal vector
419
421 MoFEM::Interface &m_field, const bool &set_singular_coordinates,
422 const Range &crack_front_nodes, const Range &crack_front_nodes_edges,
423 const Range &crack_front_elements, PetscBool add_singularity,
424 const std::string field_name, const Range tris,
425 boost::ptr_vector<MethodForForceScaling> &methods_op,
426 boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
427 &analytical_force_op);
428
429 VectorDouble Nf; //< Local force vector
430
433};
434
435struct FirendVolumeOnSide : public VolumeElementForcesAndSourcesCoreOnSide {
436 using VolumeElementForcesAndSourcesCoreOnSide::
437 VolumeElementForcesAndSourcesCoreOnSide;
440 using VolumeElementForcesAndSourcesCoreOnSide::
441 calHierarchicalBaseFunctionsOnElement;
442};
443
446
447 const Range tRis;
448 boost::ptr_vector<MethodForForceScaling> &methodsOp;
449 boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
451
454 boost::shared_ptr<MatrixDouble> hG; ///< spatial deformation gradient
455 boost::shared_ptr<MatrixDouble> HG; ///< spatial deformation gradient
456 VectorDouble sNrm; ///< Length of the normal vector
459
461
463 MoFEM::Interface &m_field, const bool &set_singular_coordinates,
464 const Range &crack_front_nodes, const Range &crack_front_nodes_edges,
465 const Range &crack_front_elements, PetscBool add_singularity,
466 const std::string field_name, const Range tris,
467 boost::ptr_vector<MethodForForceScaling> &methods_op,
468 boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
469 &analytical_force_op,
470 Range *forces_only_on_entities_row = NULL);
471
472 VectorDouble Nf; //< Local force vector
473
476};
477
479 : public OpCalculateVectorFieldGradient<3, 3> {
480
483
485 const std::string field_name, boost::shared_ptr<MatrixDouble> data_at_pts,
486 bool &singular_element, MatrixDouble &inv_s_jac)
487 : OpCalculateVectorFieldGradient<3, 3>(field_name, data_at_pts),
488 singularElement(singular_element), invSJac(inv_s_jac) {}
489
492};
493
496
499
501 const std::string field_name,
502 std::vector<VectorDouble> &values_at_gauss_pts,
503 std::vector<MatrixDouble> &gradient_at_gauss_pts, bool &singular_element,
504 MatrixDouble &inv_s_jac)
506 field_name, values_at_gauss_pts, gradient_at_gauss_pts),
507 singularElement(singular_element), invSJac(inv_s_jac) {}
510};
511
515 const std::string field_name,
516 NonlinearElasticElement::CommonData &common_data, bool &singular_element,
517 MatrixDouble &inv_s_jac)
519 common_data.dataAtGaussPts[field_name],
520 common_data.gradAtGaussPts[field_name],
521 singular_element, inv_s_jac) {}
522};
523
529 std::vector<EntityHandle> &mapGaussPts;
530 boost::shared_ptr<MatrixDouble> H;
531 OpPostProcDisplacements(bool &singular_element,
532 MatrixDouble &singular_ref_coords,
533 moab::Interface &post_proc_mesh,
534 std::vector<EntityHandle> &map_gauss_pts,
535 boost::shared_ptr<MatrixDouble> &H)
536 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
537 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
538 singularElement(singular_element), singularDisp(singular_ref_coords),
539 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), H(H) {}
542};
543
550 const bool applyDet;
552 bool &singular_element, VectorDouble &det_s, MatrixDouble &inv_s_jac,
554 bool apply_det = true)
555 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(H1),
556 singularElement(singular_element), detS(det_s), invSJac(inv_s_jac),
557 bAse(base), applyDet(apply_det) {}
561};
562
563/**
564* \brief Calculate explicit derivative of energy
565
566\f[
567\left(
568\frac{\partial \Psi}{\partial \mathbf{X}}
569\right)_\textrm{expl.}=
570\left.\left(
571\frac{\partial \Psi}{\partial \mathbf{X}}
572\right)\right|_{\mathbf{u}=const.,\mathbf{F}=const.}
573\f]
574]
575
576*/
579 HookeElement::DataAtIntegrationPts &commonData;
580 boost::shared_ptr<VectorDouble> rhoAtGaussPts;
581 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPts;
583 const double rHo0;
584 const double boneN;
586 ublas::vector<int> iNdices;
587
589 HookeElement::DataAtIntegrationPts &common_data,
590 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
591 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
592 Range tets_in_block, const double rho_0, const double bone_n,
593 Range *forces_only_on_entities_row = NULL);
595 MoFEMErrorCode doWork(int row_side, EntityType row_type,
597};
598/**
599* \brief Calculate explicit derivative of energy
600
601\f[
602\left(
603\frac{\partial \Psi}{\partial \mathbf{X}}
604\right)_\textrm{expl.}=
605\left.\left(
606\frac{\partial \Psi}{\partial \mathbf{X}}
607\right)\right|_{\mathbf{u}=const.,\mathbf{F}=const.}
608\f]
609]
610
611*/
614 HookeElement::DataAtIntegrationPts &commonData;
615 boost::shared_ptr<VectorDouble> rhoAtGaussPts;
616 boost::shared_ptr<MatrixDouble> diffRhoAtGaussPts;
617 boost::shared_ptr<MatrixDouble> diffDiffRhoAtGaussPts;
620 const double rHo0;
621 const double boneN;
623 ublas::vector<int> iNdices;
626 HookeElement::DataAtIntegrationPts &common_data,
627 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
628 boost::shared_ptr<MatrixDouble> diff_rho_at_gauss_pts,
629 boost::shared_ptr<MatrixDouble> diff_diff_rho_at_gauss_pts,
630 MatrixDouble &singular_displ, Range tets_in_block, const double rho_0,
631 const double bone_n, Range *forces_only_on_entities_row = NULL);
632 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
633 EntityType col_type,
636};
637
640 HookeElement::DataAtIntegrationPts &commonData;
641 boost::shared_ptr<VectorDouble> rhoAtGaussPts;
642 boost::shared_ptr<MatrixDouble> diffRhoAtGaussPts;
643 boost::shared_ptr<MatrixDouble> diffDiffRhoAtGaussPts;
646 const double rHo0;
647 const double boneN;
649 ublas::vector<int> iNdices;
652 HookeElement::DataAtIntegrationPts &common_data,
653 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
654 boost::shared_ptr<MatrixDouble> diff_rho_at_gauss_pts,
655 boost::shared_ptr<MatrixDouble> diff_diff_rho_at_gauss_pts,
656 MatrixDouble &singular_displ, Range tets_in_block, const double rho_0,
657 const double bone_n, Range *forces_only_on_entities_row = NULL);
658 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
659 EntityType col_type,
662};
663
666
667 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
668 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
669 boost::shared_ptr<MatrixDouble> singularDisplacement;
670
671 const double rhoN;
672 const double rHo0;
673
674 // Finite element stiffness sub-matrix K_ij
678
679 boost::shared_ptr<HookeElement::DataAtIntegrationPts> dataAtPts;
680
683
684 int nbRows; ///< number of dofs on rows
685 int nbCols; ///< number if dof on column
686 int nbIntegrationPts; ///< number of integration points
687 bool isDiag; ///< true if this block is on diagonal
688
690 const std::string row_field, const std::string col_field,
691 boost::shared_ptr<HookeElement::DataAtIntegrationPts> &data_at_pts,
692 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
693 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, const double rho_n,
694 const double rho_0,
695 boost::shared_ptr<MatrixDouble> singular_displacement);
696
697protected:
698 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
699 EntityType col_type,
702
705
708};
709
712
713 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
714 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
715 boost::shared_ptr<MatrixDouble> singularDisplacement;
716
717 const double rhoN;
718 const double rHo0;
719
720 // Finite element stiffness sub-matrix K_ij
724
725 boost::shared_ptr<HookeElement::DataAtIntegrationPts> dataAtPts;
726
729
730 int nbRows; ///< number of dofs on rows
731 int nbCols; ///< number if dof on column
732 int nbIntegrationPts; ///< number of integration points
733 bool isDiag; ///< true if this block is on diagonal
734
736 const std::string row_field, const std::string col_field,
737 boost::shared_ptr<HookeElement::DataAtIntegrationPts> &data_at_pts,
738 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
739 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, const double rho_n,
740 const double rho_0,
741 boost::shared_ptr<MatrixDouble> singular_displacement);
742
743protected:
744 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
745 EntityType col_type,
748
751
754};
755
756/**
757 * \brief Op to generate artificial density field
758 *
759 */
761
762 boost::shared_ptr<MatrixDouble> matCoordsPtr;
763 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
764 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
765 boost::shared_ptr<MatrixDouble> rhoGradGradAtGaussPtsPtr;
766 boost::shared_ptr<CrackFrontElement> feSingularPtr;
768 boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
770 const std::string row_field,
771 boost::shared_ptr<MatrixDouble> mat_coords_ptr,
772 boost::shared_ptr<VectorDouble> density_at_pts,
773 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr,
774 boost::shared_ptr<MatrixDouble> rho_grad_grad_at_gauss_pts_ptr,
775 boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
776 MatrixDouble &singular_disp,
777 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr)
778 : HookeElement::VolUserDataOperator(row_field, OPROW),
779 matCoordsPtr(mat_coords_ptr), rhoAtGaussPtsPtr(density_at_pts),
780 rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts_ptr),
781 rhoGradGradAtGaussPtsPtr(rho_grad_grad_at_gauss_pts_ptr),
782 feSingularPtr(fe_singular_ptr), singularDisp(singular_disp),
783 matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr) {}
784
785 MoFEMErrorCode doWork(int row_side, EntityType row_type,
786 HookeElement::EntData &row_data);
787};
788
789/**
790 * \brief Mark crack surfaces on skin
791 *
792 */
795
797 std::vector<EntityHandle> &mapGaussPts;
798 const string tagName;
801
803 vector<EntityHandle> &map_gauss_pts, Range &range_to_tag,
804 const string &tag_name, int tag_value = 1)
806 "MESH_NODE_POSITIONS", UserDataOperator::OPCOL),
807 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
808 rangeToTag(range_to_tag), tagName(tag_name), tagValue(tag_value) {}
809
810 MoFEMErrorCode doWork(int row_side, EntityType row_type,
812};
813
814} // namespace FractureMechanics
815
816#endif // __CRACK_FRONT_ELEMENT_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ 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
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
Definition: elasticity.cpp:76
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore > CrackFrontElement
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1210
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199
constexpr auto field_name
MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl(TwoType< E, VolumeElementForcesAndSourcesCore >)
Partial specialization for volume elements.
MoFEMErrorCode setGaussPts(int order)
Generate specific user integration quadrature.
int setGaussPtsImpl(TwoType< E, B >, int order)
MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl(TwoType< E, B >)
MoFEMErrorCode calculateHOJacobianImpl(TwoType< E, VolumeElementForcesAndSourcesCore >)
Partial specialization for volume element.
int getRuleImpl(TwoType< E, B >, int order)
CrackFrontSingularBase(MoFEM::Interface &m_field, const bool &set_singular_coordinates, const Range &crack_front_nodes, const Range &crack_front_nodes_edges, const Range &crack_front_elements, PetscBool add_singularity)
CrackFrontSingularBase(MoFEM::Interface &m_field)
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator UserDataOperator
OpAleLhsWithDensitySingularElement_dX_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< HookeElement::DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0, boost::shared_ptr< MatrixDouble > singular_displacement)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< HookeElement::DataAtIntegrationPts > dataAtPts
MoFEMErrorCode aSsemble(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)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< HookeElement::DataAtIntegrationPts > dataAtPts
OpAleLhsWithDensitySingularElement_dx_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< HookeElement::DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0, boost::shared_ptr< MatrixDouble > singular_displacement)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAnalyticalMaterialTraction(MoFEM::Interface &m_field, const bool &set_singular_coordinates, const Range &crack_front_nodes, const Range &crack_front_nodes_edges, const Range &crack_front_elements, PetscBool add_singularity, const std::string field_name, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > &analytical_force_op, Range *forces_only_on_entities_row=NULL)
boost::shared_ptr< MatrixDouble > hG
spatial deformation gradient
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
VectorDouble sNrm
Length of the normal vector.
boost::shared_ptr< MatrixDouble > HG
spatial deformation gradient
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CrackFrontSingularBase< FirendVolumeOnSide, VolumeElementForcesAndSourcesCore > volSideFe
boost::ptr_vector< MethodForForceScaling > & methodsOp
CrackFrontSingularBase< VolumeElementForcesAndSourcesCoreOnSide, VolumeElementForcesAndSourcesCore > volSideFe
VectorDouble sNrm
Length of the normal vector.
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
OpAnalyticalSpatialTraction(MoFEM::Interface &m_field, const bool &set_singular_coordinates, const Range &crack_front_nodes, const Range &crack_front_nodes_edges, const Range &crack_front_elements, PetscBool add_singularity, const std::string field_name, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > &analytical_force_op)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetCrackFrontCommonDataAtGaussPts(const std::string field_name, NonlinearElasticElement::CommonData &common_data, bool &singular_element, MatrixDouble &inv_s_jac)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetCrackFrontDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts, bool &singular_element, MatrixDouble &inv_s_jac)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetCrackFrontDataGradientAtGaussPts(const std::string field_name, boost::shared_ptr< MatrixDouble > data_at_pts, bool &singular_element, MatrixDouble &inv_s_jac)
Op to generate artificial density field.
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
boost::shared_ptr< MatrixDouble > matGradPosAtPtsPtr
boost::shared_ptr< CrackFrontElement > feSingularPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > matCoordsPtr
boost::shared_ptr< MatrixDouble > rhoGradGradAtGaussPtsPtr
OpGetDensityFieldForTesting(const std::string row_field, boost::shared_ptr< MatrixDouble > mat_coords_ptr, boost::shared_ptr< VectorDouble > density_at_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts_ptr, boost::shared_ptr< MatrixDouble > rho_grad_grad_at_gauss_pts_ptr, boost::shared_ptr< CrackFrontElement > fe_singular_ptr, MatrixDouble &singular_disp, boost::shared_ptr< MatrixDouble > mat_grad_pos_at_pts_ptr)
OpLhsBoneExplicitDerivariveWithHooke_dX(HookeElement::DataAtIntegrationPts &common_data, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_diff_rho_at_gauss_pts, MatrixDouble &singular_displ, Range tets_in_block, const double rho_0, const double bone_n, Range *forces_only_on_entities_row=NULL)
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)
OpLhsBoneExplicitDerivariveWithHooke_dx(HookeElement::DataAtIntegrationPts &common_data, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_diff_rho_at_gauss_pts, MatrixDouble &singular_displ, Range tets_in_block, const double rho_0, const double bone_n, Range *forces_only_on_entities_row=NULL)
OpPostProcDisplacements(bool &singular_element, MatrixDouble &singular_ref_coords, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< MatrixDouble > &H)
moab::Interface & postProcMesh
MatrixDouble & singularDisp
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
bool & singularElement
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< MatrixDouble > H
OpRhsBoneExplicitDerivariveWithHooke(HookeElement::DataAtIntegrationPts &common_data, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, Range tets_in_block, const double rho_0, const double bone_n, Range *forces_only_on_entities_row=NULL)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpSetTagRangeOnSkin(moab::Interface &post_proc_mesh, vector< EntityHandle > &map_gauss_pts, Range &range_to_tag, const string &tag_name, int tag_value=1)
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpTransfromSingularBaseFunctions(bool &singular_element, VectorDouble &det_s, MatrixDouble &inv_s_jac, FieldApproximationBase base=AINSWORTH_LOBATTO_BASE, bool apply_det=true)
Deprecated interface functions.
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
common data used by volume elements
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
structure grouping operators and data used for calculation of nonlinear elastic element