v0.14.0
Public Member Functions | Public Attributes | Private Attributes | List of all members
FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT > Struct Template Reference

#include <users_modules/fracture_mechanics/src/CrackFrontElement.hpp>

Inheritance diagram for FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >:
[legend]
Collaboration diagram for FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >:
[legend]

Public Member Functions

MoFEMErrorCode getElementOptions ()
 
 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)
 
virtual ~CrackFrontSingularBase ()=default
 
MoFEMErrorCode operator() ()
 
MoFEMErrorCode getSpaceBaseAndOrderOnElement ()
 
template<typename E , typename B >
MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl (TwoType< E, B >)
 
template<typename E >
MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl (TwoType< E, VolumeElementForcesAndSourcesCore >)
 Partial specialization for volume elements. More...
 
MoFEMErrorCode calculateHOJacobian ()
 
template<typename E >
MoFEMErrorCode calculateHOJacobianImpl (TwoType< E, VolumeElementForcesAndSourcesCore >)
 Partial specialization for volume element. More...
 
int getRule (int order)
 
template<typename E , typename B >
int getRuleImpl (TwoType< E, B >, int order)
 
MoFEMErrorCode setGaussPts (int order)
 Generate specific user integration quadrature. More...
 
template<typename E , typename B >
int setGaussPtsImpl (TwoType< E, B >, int order)
 
int getRuleImpl (TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >, int order)
 
MoFEMErrorCode setGaussPtsImpl (TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >, int order)
 
int getRuleImpl (TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >, int order)
 
MoFEMErrorCode setGaussPtsImpl (TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >, int order)
 

Public Attributes

Range crackFrontNodes
 
Range crackFrontNodesEdges
 
Range crackFrontElements
 
MatrixDouble sJac
 
MatrixDouble invSJac
 
MatrixDouble singularDisp
 
MatrixDouble singularRefCoords
 
VectorDouble detS
 
bool singularElement
 
bool singularEdges [6]
 
bool singularNodes [4]
 
bool setSingularCoordinates
 
PetscBool addSingularity
 
PetscBool refineIntegration
 
int refinementLevels
 

Private Attributes

const bool & setSingularCoordinatesPriv
 
MatrixDouble refCoords
 
MatrixDouble refGaussCoords
 

Detailed Description

template<typename ELEMENT, typename BASE_ELEMENT>
struct FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >

Use idea from [10] to approximate singularity at crack tip. The idea is to shift mid noted to quarter edge length to crack tip on edges adjacent to crack front.

Definition at line 35 of file CrackFrontElement.hpp.

Constructor & Destructor Documentation

◆ CrackFrontSingularBase() [1/2]

template<typename ELEMENT , typename BASE_ELEMENT >
FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::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 
)
inline

Definition at line 53 of file CrackFrontElement.hpp.

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) {
64  ierr = getElementOptions();
65  CHKERRABORT(PETSC_COMM_WORLD, ierr);
66  }

◆ CrackFrontSingularBase() [2/2]

template<typename ELEMENT , typename BASE_ELEMENT >
FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::CrackFrontSingularBase ( MoFEM::Interface &  m_field)
inline

Definition at line 69 of file CrackFrontElement.hpp.

70  : ELEMENT(m_field), setSingularCoordinatesPriv(setSingularCoordinates) {}

◆ ~CrackFrontSingularBase()

template<typename ELEMENT , typename BASE_ELEMENT >
virtual FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::~CrackFrontSingularBase ( )
virtualdefault

Member Function Documentation

◆ calculateHOJacobian()

template<typename ELEMENT , typename BASE_ELEMENT >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::calculateHOJacobian ( )
inline

Definition at line 146 of file CrackFrontElement.hpp.

146  {
147  return calculateHOJacobianImpl(
148  TwoType<ELEMENT, VolumeElementForcesAndSourcesCore>());
149  }

◆ calculateHOJacobianImpl()

template<typename ELEMENT , typename BASE_ELEMENT >
template<typename E >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::calculateHOJacobianImpl ( TwoType< E, VolumeElementForcesAndSourcesCore >  )
inline

Partial specialization for volume element.

Definition at line 154 of file CrackFrontElement.hpp.

154  {
155  MoFEMFunctionBegin;
156 
157  if (setSingularCoordinatesPriv) {
158  singularElement = false;
159  MoFEMFunctionReturnHot(0);
160  }
161  if (addSingularity == PETSC_FALSE) {
162  singularElement = false;
163  MoFEMFunctionReturnHot(0);
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}};
170  FTensor::Index<'i', 3> i;
171  FTensor::Index<'j', 3> j;
172  FTensor::Index<'k', 3> k;
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  {
186  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_s_jac(
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);
202  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_diff_base_n(
203  diff_base_n, &diff_base_n[1], &diff_base_n[2]);
204  double diff_n[12];
205  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_diff_n(
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
260  FTensor::Tensor1<double, 3> t_dir;
261  t_dir(i) = this->tJac(i, j) * t_base_dir(j);
262  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_s_jac(
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]);
266  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
267  t_singular_displacement(&singularDisp(0, 0), &singularDisp(0, 1),
268  &singularDisp(0, 2));
269  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
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);
276  FTensor::Tensor1<double, 3> t_diff_n(
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  {
294  FTensor::Tensor2<double *, 3, 3> t_s_jac(
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);
301  FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
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 
319  MoFEMFunctionReturn(0);
320  }

◆ getElementOptions()

template<typename ELEMENT , typename BASE_ELEMENT >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getElementOptions

Definition at line 361 of file CrackFrontElement.hpp.

361  {
362  MoFEMFunctionBeginHot;
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);
379  MoFEMFunctionReturnHot(0);
380 }

◆ getRule()

template<typename ELEMENT , typename BASE_ELEMENT >
int FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getRule ( int  order)
inline
Returns
returning negative number -1, i.e. special user quadrature is generated

Definition at line 326 of file CrackFrontElement.hpp.

326  {
327  return getRuleImpl(TwoType<ELEMENT, BASE_ELEMENT>(), order);
328  }

◆ getRuleImpl() [1/3]

template<typename ELEMENT , typename BASE_ELEMENT >
template<typename E , typename B >
int FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getRuleImpl ( TwoType< E, B >  ,
int  order 
)
inline

Definition at line 331 of file CrackFrontElement.hpp.

331  {
332  return E::getRule(order);
333  }

◆ getRuleImpl() [2/3]

int FractureMechanics::CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >::getRuleImpl ( TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >  ,
int  order 
)

Definition at line 59 of file CrackFrontElement.cpp.

61  {
62  return -1;
63 }

◆ getRuleImpl() [3/3]

int FractureMechanics::CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >::getRuleImpl ( TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >  ,
int  order 
)

◆ getSpaceBaseAndOrderOnElement()

template<typename ELEMENT , typename BASE_ELEMENT >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getSpaceBaseAndOrderOnElement ( )
inline

Definition at line 103 of file CrackFrontElement.hpp.

103  {
104  return getSpaceBaseAndOrderOnElementImpl(TwoType<ELEMENT, BASE_ELEMENT>());
105  }

◆ getSpaceBaseAndOrderOnElementImpl() [1/2]

template<typename ELEMENT , typename BASE_ELEMENT >
template<typename E , typename B >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getSpaceBaseAndOrderOnElementImpl ( TwoType< E, B >  )
inline

Definition at line 108 of file CrackFrontElement.hpp.

108  {
109  return E::getSpaceBaseAndOrderOnElement();
110  }

◆ getSpaceBaseAndOrderOnElementImpl() [2/2]

template<typename ELEMENT , typename BASE_ELEMENT >
template<typename E >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getSpaceBaseAndOrderOnElementImpl ( TwoType< E, VolumeElementForcesAndSourcesCore >  )
inline

Partial specialization for volume elements.

Definition at line 114 of file CrackFrontElement.hpp.

115  {
116  MoFEMFunctionBeginHot;
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  }
143  MoFEMFunctionReturnHot(0);
144  }

◆ operator()()

template<typename ELEMENT , typename BASE_ELEMENT >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::operator() ( )
inline

Definition at line 74 of file CrackFrontElement.hpp.

74  {
75  MoFEMFunctionBegin;
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(MBTET);
85  CHKERR this->calculateVolumeAndJacobian();
86  CHKERR this->getSpaceBaseAndOrderOnElement();
87  CHKERR this->setIntegrationPts();
88  CHKERR this->calculateCoordinatesAtGaussPts();
89  CHKERR this->calHierarchicalBaseFunctionsOnElement();
90  CHKERR this->calculateHOJacobian();
91  CHKERR this->transformBaseFunctions();
92 
93  // Iterate over operators
94  CHKERR this->loopOverOperators();
95 
96  MoFEMFunctionReturn(0);
97  }

◆ setGaussPts()

template<typename ELEMENT , typename BASE_ELEMENT >
MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::setGaussPts ( int  order)
inline

Generate specific user integration quadrature.

Element is refined at the crack top and quadrature is set to each of sub-tets.

Definition at line 342 of file CrackFrontElement.hpp.

342  {
343  return setGaussPtsImpl(TwoType<ELEMENT, BASE_ELEMENT>(), order);
344  }

◆ setGaussPtsImpl() [1/3]

template<typename ELEMENT , typename BASE_ELEMENT >
template<typename E , typename B >
int FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::setGaussPtsImpl ( TwoType< E, B >  ,
int  order 
)
inline

Definition at line 348 of file CrackFrontElement.hpp.

348  {
349  return E::setGaussPts(order);
350  }

◆ setGaussPtsImpl() [2/3]

MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >::setGaussPtsImpl ( TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >  ,
int  order 
)

Definition at line 69 of file CrackFrontElement.cpp.

71  {
72 
73  MoFEMFunctionBeginHot;
74 
75  auto set_gauss_pts = [this]() {
76  MoFEMFunctionBegin;
77  gaussPts.swap(refGaussCoords);
78  const size_t nb_gauss_pts = gaussPts.size2();
79  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
80  double *shape_ptr =
81  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
82  CHKERR ShapeMBTET(shape_ptr, &gaussPts(0, 0), &gaussPts(1, 0),
83  &gaussPts(2, 0), nb_gauss_pts);
84  MoFEMFunctionReturn(0);
85  };
86 
87  std::bitset<4> singular_nodes_ids;
88  if (refineIntegration == PETSC_TRUE && singularElement) {
89 
90  for (int nn = 0; nn != 4; nn++) {
91  if (singularNodes[nn])
92  singular_nodes_ids.set(nn);
93  }
94 
95  auto it_map_ref_coords = mapRefCoords.find(singular_nodes_ids.to_ulong());
96  if (it_map_ref_coords != mapRefCoords.end()) {
97 
98  refGaussCoords = it_map_ref_coords->second;
99  CHKERR set_gauss_pts();
100 
101  MoFEMFunctionReturnHot(0);
102  }
103  }
104 
105  order = order ? order : 1;
106  int rule = (singularElement) ? 2 * (std::max(order, 3) - 1) : 2 * (order - 1);
107 
108  if (rule < QUAD_3D_TABLE_SIZE) {
109  if (QUAD_3D_TABLE[rule]->dim != 3) {
110  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "wrong dimension");
111  }
112  if (QUAD_3D_TABLE[rule]->order < rule) {
113  SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
114  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
115  }
116  const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
117  gaussPts.resize(4, nb_gauss_pts, false);
118  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4, &gaussPts(0, 0),
119  1);
120  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4, &gaussPts(1, 0),
121  1);
122  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4, &gaussPts(2, 0),
123  1);
124  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1, &gaussPts(3, 0),
125  1);
126  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
127  false);
128  double *shape_ptr =
129  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
130  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr, 1);
131  } else {
132  SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
133  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
134  }
135 
136  // Refine quadrature by splitting edges adjacent to crack front. On sub
137  // triangle place quadrature points.
138  if (!singularElement)
139  MoFEMFunctionReturnHot(0);
140 
141  // refine mesh on refine elements
142  if (refineIntegration == PETSC_TRUE) {
143 
144  const int max_level = refinementLevels;
145  EntityHandle tet;
146 
147  moab::Core moab_ref;
148  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
149  EntityHandle nodes[4];
150  for (int nn = 0; nn != 4; nn++)
151  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
152  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
153  MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
154  MoFEM::Interface &m_field_ref = mofem_ref_core;
155  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
156  0, 3, BitRefLevel().set(0));
157 
158  Range singular_nodes;
159  for (int nn = 0; nn != 4; nn++) {
160  if (singularNodes[nn]) {
161  EntityHandle ent;
162  CHKERR moab_ref.side_element(tet, 0, nn, ent);
163  singular_nodes.insert(ent);
164  }
165  }
166 
167  EntityHandle meshset;
168  if (singular_nodes.size() > 1) {
169  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
170  for (int ee = 0; ee != 6; ee++) {
171  if (singularEdges[ee]) {
172  EntityHandle ent;
173  CHKERR moab_ref.side_element(tet, 1, ee, ent);
174  CHKERR moab_ref.add_entities(meshset, &ent, 1);
175  }
176  }
177  }
178 
179  MeshRefinement *m_ref;
180  CHKERR m_field_ref.getInterface(m_ref);
181  for (int ll = 0; ll != max_level; ll++) {
182  // PetscPrintf(T::mField.get_comm(),"Refine Level %d\n",ll);
183  Range edges;
184  CHKERR m_field_ref.getInterface<BitRefManager>()
185  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
186  BitRefLevel().set(), MBEDGE, edges);
187  Range ref_edges;
188  CHKERR moab_ref.get_adjacencies(singular_nodes, 1, true, ref_edges,
189  moab::Interface::UNION);
190  ref_edges = intersect(ref_edges, edges);
191  if (singular_nodes.size() > 1) {
192  Range ents;
193  CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
194  ref_edges = intersect(ref_edges, ents);
195  }
196  Range tets;
197  CHKERR m_field_ref.getInterface<BitRefManager>()
198  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
199  BitRefLevel().set(), MBTET, tets);
200  // refine mesh
201  CHKERR m_ref->addVerticesInTheMiddleOfEdges(
202  ref_edges, BitRefLevel().set(ll + 1));
203  CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
204  if (singular_nodes.size() > 1) {
205  CHKERR m_field_ref.getInterface<BitRefManager>()
206  ->updateMeshsetByEntitiesChildren(
207  meshset, BitRefLevel().set(ll + 1), meshset, MBEDGE, true);
208  }
209  }
210 
211  Range tets;
212  CHKERR m_field_ref.getInterface<BitRefManager>()
213  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
214  BitRefLevel().set(), MBTET, tets);
215  refCoords.resize(tets.size(), 12, false);
216  int tt = 0;
217  for (Range::iterator tit = tets.begin(); tit != tets.end(); tit++, tt++) {
218  int num_nodes;
219  const EntityHandle *conn;
220  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
221  CHKERR moab_ref.get_coords(conn, num_nodes, &refCoords(tt, 0));
222  }
223 
224  if (0) {
225  EntityHandle out_meshset;
226  CHKERR moab_ref.create_meshset(MESHSET_SET, out_meshset);
227  CHKERR moab_ref.add_entities(out_meshset, tets);
228  std::string file_name =
229  "ref_tet_" + boost::lexical_cast<std::string>(nInTheLoop) + ".vtk";
230  CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &out_meshset, 1);
231  CHKERR moab_ref.delete_entities(&out_meshset, 1);
232  }
233 
234  // Set integration points
235  double diff_n[12];
236  CHKERR ShapeDiffMBTET(diff_n);
237  const size_t nb_gauss_pts = gaussPts.size2();
238  refGaussCoords.resize(4, nb_gauss_pts * refCoords.size1());
239  MatrixDouble &shape_n = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
240  int gg = 0;
241  for (size_t tt = 0; tt != refCoords.size1(); tt++) {
242  double *tet_coords = &refCoords(tt, 0);
243  double det = ShapeVolumeMBTET(diff_n, tet_coords);
244  det *= 6;
245  for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
246  for (int dd = 0; dd != 3; dd++) {
247  refGaussCoords(dd, gg) = shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
248  shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
249  shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
250  shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
251  }
252  refGaussCoords(3, gg) = gaussPts(3, ggg) * det;
253  }
254  }
255 
256  mapRefCoords[singular_nodes_ids.to_ulong()] = refGaussCoords;
257  CHKERR set_gauss_pts();
258  }
259 
260  MoFEMFunctionReturnHot(0);
261 }

◆ setGaussPtsImpl() [3/3]

MoFEMErrorCode FractureMechanics::CrackFrontSingularBase< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >::setGaussPtsImpl ( TwoType< NonlinearElasticElement::MyVolumeFE, VolumeElementForcesAndSourcesCore >  ,
int  order 
)

Member Data Documentation

◆ addSingularity

template<typename ELEMENT , typename BASE_ELEMENT >
PetscBool FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::addSingularity

Definition at line 99 of file CrackFrontElement.hpp.

◆ crackFrontElements

template<typename ELEMENT , typename BASE_ELEMENT >
Range FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::crackFrontElements

Definition at line 39 of file CrackFrontElement.hpp.

◆ crackFrontNodes

template<typename ELEMENT , typename BASE_ELEMENT >
Range FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::crackFrontNodes

Definition at line 37 of file CrackFrontElement.hpp.

◆ crackFrontNodesEdges

template<typename ELEMENT , typename BASE_ELEMENT >
Range FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::crackFrontNodesEdges

Definition at line 38 of file CrackFrontElement.hpp.

◆ detS

template<typename ELEMENT , typename BASE_ELEMENT >
VectorDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::detS

Definition at line 45 of file CrackFrontElement.hpp.

◆ invSJac

template<typename ELEMENT , typename BASE_ELEMENT >
MatrixDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::invSJac

Definition at line 42 of file CrackFrontElement.hpp.

◆ refCoords

template<typename ELEMENT , typename BASE_ELEMENT >
MatrixDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::refCoords
private

Definition at line 355 of file CrackFrontElement.hpp.

◆ refGaussCoords

template<typename ELEMENT , typename BASE_ELEMENT >
MatrixDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::refGaussCoords
private

Definition at line 356 of file CrackFrontElement.hpp.

◆ refineIntegration

template<typename ELEMENT , typename BASE_ELEMENT >
PetscBool FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::refineIntegration

Definition at line 100 of file CrackFrontElement.hpp.

◆ refinementLevels

template<typename ELEMENT , typename BASE_ELEMENT >
int FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::refinementLevels

Definition at line 101 of file CrackFrontElement.hpp.

◆ setSingularCoordinates

template<typename ELEMENT , typename BASE_ELEMENT >
bool FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::setSingularCoordinates

Definition at line 68 of file CrackFrontElement.hpp.

◆ setSingularCoordinatesPriv

template<typename ELEMENT , typename BASE_ELEMENT >
const bool& FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::setSingularCoordinatesPriv
private

Definition at line 353 of file CrackFrontElement.hpp.

◆ singularDisp

template<typename ELEMENT , typename BASE_ELEMENT >
MatrixDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::singularDisp

Definition at line 43 of file CrackFrontElement.hpp.

◆ singularEdges

template<typename ELEMENT , typename BASE_ELEMENT >
bool FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::singularEdges[6]

Definition at line 48 of file CrackFrontElement.hpp.

◆ singularElement

template<typename ELEMENT , typename BASE_ELEMENT >
bool FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::singularElement

Definition at line 47 of file CrackFrontElement.hpp.

◆ singularNodes

template<typename ELEMENT , typename BASE_ELEMENT >
bool FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::singularNodes[4]

Definition at line 49 of file CrackFrontElement.hpp.

◆ singularRefCoords

template<typename ELEMENT , typename BASE_ELEMENT >
MatrixDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::singularRefCoords

Definition at line 44 of file CrackFrontElement.hpp.

◆ sJac

template<typename ELEMENT , typename BASE_ELEMENT >
MatrixDouble FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::sJac

Definition at line 41 of file CrackFrontElement.hpp.

The documentation for this struct was generated from the following file:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
FractureMechanics::CrackFrontSingularBase::crackFrontNodes
Range crackFrontNodes
Definition: CrackFrontElement.hpp:37
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::MeshRefinement::addVerticesInTheMiddleOfEdges
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
Definition: MeshRefinement.cpp:42
FractureMechanics::CrackFrontSingularBase::setSingularCoordinates
bool setSingularCoordinates
Definition: CrackFrontElement.hpp:68
FractureMechanics::CrackFrontSingularBase::singularElement
bool singularElement
Definition: CrackFrontElement.hpp:47
FractureMechanics::CrackFrontSingularBase::refineIntegration
PetscBool refineIntegration
Definition: CrackFrontElement.hpp:100
FractureMechanics::CrackFrontSingularBase::refGaussCoords
MatrixDouble refGaussCoords
Definition: CrackFrontElement.hpp:356
NOBASE
@ NOBASE
Definition: definitions.h:59
ELEMENT
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
QUAD_::order
int order
Definition: quad.h:28
FractureMechanics::CrackFrontSingularBase::singularNodes
bool singularNodes[4]
Definition: CrackFrontElement.hpp:49
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
FractureMechanics::CrackFrontSingularBase::invSJac
MatrixDouble invSJac
Definition: CrackFrontElement.hpp:42
order
constexpr int order
Definition: dg_projection.cpp:18
FractureMechanics::CrackFrontSingularBase::calculateHOJacobian
MoFEMErrorCode calculateHOJacobian()
Definition: CrackFrontElement.hpp:146
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FractureMechanics::CrackFrontSingularBase::getSpaceBaseAndOrderOnElement
MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Definition: CrackFrontElement.hpp:103
MoFEM::invertTensor3by3
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:1559
MoFEM::MeshRefinement::refineTets
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
Definition: MeshRefinement.cpp:197
FractureMechanics::CrackFrontSingularBase::crackFrontNodesEdges
Range crackFrontNodesEdges
Definition: CrackFrontElement.hpp:38
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::CoreTmp
Definition: Core.hpp:36
FractureMechanics::CrackFrontSingularBase::getSpaceBaseAndOrderOnElementImpl
MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl(TwoType< E, B >)
Definition: CrackFrontElement.hpp:108
FractureMechanics::CrackFrontSingularBase::sJac
MatrixDouble sJac
Definition: CrackFrontElement.hpp:41
FractureMechanics::CrackFrontSingularBase::singularDisp
MatrixDouble singularDisp
Definition: CrackFrontElement.hpp:43
FractureMechanics::CrackFrontSingularBase::addSingularity
PetscBool addSingularity
Definition: CrackFrontElement.hpp:99
FractureMechanics::CrackFrontSingularBase::detS
VectorDouble detS
Definition: CrackFrontElement.hpp:45
QUAD_::npoints
int npoints
Definition: quad.h:29
ShapeVolumeMBTET
double ShapeVolumeMBTET(double *diffN, const double *coords)
calculate TET volume
Definition: fem_tools.c:298
FractureMechanics::CrackFrontSingularBase::getElementOptions
MoFEMErrorCode getElementOptions()
Definition: CrackFrontElement.hpp:361
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
FractureMechanics::CrackFrontSingularBase::getRuleImpl
int getRuleImpl(TwoType< E, B >, int order)
Definition: CrackFrontElement.hpp:331
FractureMechanics::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: CrackFrontElement.cpp:53
FTensor::Index< 'i', 3 >
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
Range
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
FractureMechanics::CrackFrontSingularBase::refCoords
MatrixDouble refCoords
Definition: CrackFrontElement.hpp:355
FractureMechanics::CrackFrontSingularBase::calculateHOJacobianImpl
MoFEMErrorCode calculateHOJacobianImpl(TwoType< E, VolumeElementForcesAndSourcesCore >)
Partial specialization for volume element.
Definition: CrackFrontElement.hpp:154
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FractureMechanics::CrackFrontSingularBase::setGaussPtsImpl
int setGaussPtsImpl(TwoType< E, B >, int order)
Definition: CrackFrontElement.hpp:348
FractureMechanics::CrackFrontSingularBase::refinementLevels
int refinementLevels
Definition: CrackFrontElement.hpp:101
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
FractureMechanics::CrackFrontSingularBase::crackFrontElements
Range crackFrontElements
Definition: CrackFrontElement.hpp:39
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FractureMechanics::CrackFrontSingularBase::singularEdges
bool singularEdges[6]
Definition: CrackFrontElement.hpp:48
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FractureMechanics::CrackFrontSingularBase::singularRefCoords
MatrixDouble singularRefCoords
Definition: CrackFrontElement.hpp:44
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FractureMechanics::CrackFrontSingularBase::setSingularCoordinatesPriv
const bool & setSingularCoordinatesPriv
Definition: CrackFrontElement.hpp:353
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
Generated by Doxygen 1.8.17 and hosted at University of Glasgow