v0.13.1
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 boolsetSingularCoordinatesPriv
 
MatrixDouble refCoords
 
MatrixDouble refGaussCoords
 

Detailed Description

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

Use idea from [9] 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 
)

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) {
65 CHKERRABORT(PETSC_COMM_WORLD, ierr);
66 }
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76

◆ CrackFrontSingularBase() [2/2]

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

◆ ~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 ( )

Definition at line 146 of file CrackFrontElement.hpp.

146 {
148 TwoType<ELEMENT, VolumeElementForcesAndSourcesCore>());
149 }
MoFEMErrorCode calculateHOJacobianImpl(TwoType< E, VolumeElementForcesAndSourcesCore >)
Partial specialization for volume element.

◆ calculateHOJacobianImpl()

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

Partial specialization for volume element.

Definition at line 154 of file CrackFrontElement.hpp.

154 {
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 }
@ 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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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
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
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
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
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199

◆ getElementOptions()

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

Definition at line 361 of file CrackFrontElement.hpp.

361 {
363 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Get singular element options",
364 "none");
365 CHKERRQ(ierr);
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}
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440

◆ getRule()

template<typename ELEMENT , typename BASE_ELEMENT >
int FractureMechanics::CrackFrontSingularBase< ELEMENT, BASE_ELEMENT >::getRule ( int  order)
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 }
int getRuleImpl(TwoType< E, B >, int order)

◆ 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 
)

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 
)

◆ getRuleImpl() [3/3]

Definition at line 58 of file CrackFrontElement.cpp.

61 {
62 return -1;
63}

◆ getSpaceBaseAndOrderOnElement()

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

Definition at line 103 of file CrackFrontElement.hpp.

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

◆ getSpaceBaseAndOrderOnElementImpl() [1/2]

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

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 >  )

Partial specialization for volume elements.

Definition at line 114 of file CrackFrontElement.hpp.

115 {
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 }

◆ operator()()

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

Definition at line 74 of file CrackFrontElement.hpp.

74 {
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 }
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31

◆ setGaussPts()

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

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 }
int setGaussPtsImpl(TwoType< E, B >, int order)

◆ 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 
)

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 
)

◆ setGaussPtsImpl() [3/3]

Definition at line 68 of file CrackFrontElement.cpp.

71 {
72
74
75 auto set_gauss_pts = [this]() {
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);
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
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)
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
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
261}
const int dim
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
double ShapeVolumeMBTET(double *diffN, const double *coords)
calculate TET volume
Definition: fem_tools.c:298
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static std::map< long int, MatrixDouble > mapRefCoords
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
CoreTmp< 0 > Core
Definition: Core.hpp:1086
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
Managing BitRefLevels.
Deprecated interface functions.
Mesh refinement interface.
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
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
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
int order
Definition: quad.h:28
int npoints
Definition: quad.h:29

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: