v0.12.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 
21 namespace FractureMechanics {
22 
24 
25 template <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  */
34 template <typename ELEMENT, typename BASE_ELEMENT>
36 
40 
46 
48  bool singularEdges[6];
49  bool singularNodes[4];
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;
100  PetscBool refineIntegration;
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
251  FTensor::Tensor1<double, 3> t_base_dir(
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 
352 private:
354 
357 };
358 
359 template <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
383 template <>
384 template <>
389  int order);
390 
391 // specialization
392 template <>
393 template <>
398  int order);
399 
403 
406 
407  const Range tRis;
408  boost::ptr_vector<MethodForForceScaling> &methodsOp;
409  boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
411 
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 
431  MoFEMErrorCode doWork(int side, EntityType type,
433 };
434 
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 
474  MoFEMErrorCode doWork(int side, EntityType type,
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 
490  MoFEMErrorCode doWork(int side, EntityType type,
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) {}
508  MoFEMErrorCode doWork(int side, EntityType type,
510 };
511 
515  const std::string field_name,
516  NonlinearElasticElement::CommonData &common_data, bool &singular_element,
517  MatrixDouble &inv_s_jac)
518  : OpGetCrackFrontDataAtGaussPts(field_name,
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)
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) {}
540  MoFEMErrorCode doWork(int side, EntityType type,
542 };
543 
550  const bool applyDet;
552  bool &singular_element, VectorDouble &det_s, MatrixDouble &inv_s_jac,
554  bool apply_det = true)
556  singularElement(singular_element), detS(det_s), invSJac(inv_s_jac),
557  bAse(base), applyDet(apply_det) {}
559  MoFEMErrorCode doWork(int side, EntityType type,
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;
582  Range tetsInBlock;
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;
619  Range tetsInBlock;
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;
645  Range tetsInBlock;
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 
697 protected:
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 
743 protected:
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;
799  Range rangeToTag;
800  int tagValue;
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
FieldApproximationBase
approximation base
Definition: definitions.h:69
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:73
@ NOBASE
Definition: definitions.h:70
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:440
@ H1
continuous field
Definition: definitions.h:96
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:339
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:42
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:409
#define CHKERR
Inline error check.
Definition: definitions.h:528
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:433
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
Definition: elasticity.cpp:91
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
VolumeElementForcesAndSourcesCoreOnSideSwitch< 0 > VolumeElementForcesAndSourcesCoreOnSide
Volume element used to integrate on skeleton.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore > CrackFrontElement
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:990
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1949
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:1009
DataForcesAndSourcesCore::EntData EntData
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.
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