v0.13.1
VolumeLengthQuality.hpp
Go to the documentation of this file.
1/** \file VolumeLengthQuality.hpp
2 * \ingroup nonlinear_elastic_elem
3 * \brief Implementation of Volume-Length-Quality measure with barrier
4 */
5
6
7
8#ifndef __VOLUME_LENGTH_QUALITY_HPP__
9#define __VOLUME_LENGTH_QUALITY_HPP__
10
11#ifndef WITH_ADOL_C
12#error "MoFEM need to be compiled with ADOL-C"
13#endif
14
21};
22
23static const char *VolumeLengthQualityTypeNames[] = {
24 "QUALITY", "BARRIER_AND_QUALITY", "BARRIER_AND_CHANGE_QUALITY",
25 "BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME"};
26
27/** \brief Volume Length Quality
28 \ingroup nonlinear_elastic_elem
29
30 */
31template <typename TYPE>
34 TYPE> {
35
36 // VolumeLengthQualityType tYpe;
37 int tYpe;
38 double aLpha;
39 double gAmma;
40
44
47 tYpe(type), aLpha(alpha), gAmma(gamma) {
49 CHKERRABORT(PETSC_COMM_WORLD, ierr);
50 }
51
54 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
55 "Get VolumeLengthQuality material options",
56 "none");
57 CHKERR PetscOptionsEList(
58 "-volume_length_type", "Volume length quality type", "",
61 CHKERR PetscOptionsScalar("-volume_length_alpha",
62 "volume length alpha parameter", "", aLpha,
63 &aLpha, PETSC_NULL);
64 CHKERR PetscOptionsScalar("-volume_length_gamma",
65 "volume length parameter (barrier)", "", gAmma,
66 &gAmma, PETSC_NULL);
67 ierr = PetscOptionsEnd();
69 }
70
74
75 ublas::vector<TYPE> deltaX;
76 ublas::matrix<TYPE> Q, dXdChiT;
78
79 /** Get coordinates of edges using cannonical element numeration
80 */
83 if (coordsEdges.empty()) {
84 coordsEdges.resize(6 * 2 * 3, false);
85 }
86 cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
87 &coordsEdges[0 * 3 * 2 + 0], 1);
88 cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
89 &coordsEdges[0 * 3 * 2 + 3], 1);
90 cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
91 &coordsEdges[1 * 3 * 2 + 0], 1);
92 cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
93 &coordsEdges[1 * 3 * 2 + 3], 1);
94 cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
95 &coordsEdges[2 * 3 * 2 + 0], 1);
96 cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
97 &coordsEdges[2 * 3 * 2 + 3], 1);
98 cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
99 &coordsEdges[3 * 3 * 2 + 0], 1);
100 cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
101 &coordsEdges[3 * 3 * 2 + 3], 1);
102 cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
103 &coordsEdges[4 * 3 * 2 + 0], 1);
104 cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
105 &coordsEdges[4 * 3 * 2 + 3], 1);
106 cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
107 &coordsEdges[5 * 3 * 2 + 0], 1);
108 cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
109 &coordsEdges[5 * 3 * 2 + 3], 1);
111 }
112
113 /** \brief Calculate mean element edge length
114
115 \f[
116 \Delta \boldsymbol\chi = \boldsymbol\chi^1 - \boldsymbol\chi^2
117 \f]
118
119 \f[
120 \Delta X = \mathbf{F} \Delta \boldsymbol\chi
121 \f]
122
123 \f[
124 l_\textrm{rms} = \sqrt{\frac{1}{6} \sum_{i=0}^6 l_i^2 } =
125 L_\textrm{rms}dl_\textrm{rms} \f]
126
127 */
130 if (deltaChi.size() != 3) {
131 deltaChi.resize(3);
132 deltaX.resize(3);
133 dXdChiT.resize(3, 3);
134 }
135 lrmsSquared = 0;
136 lrmsSquared0 = 0;
137 dXdChiT.clear();
138 for (int ee = 0; ee < 6; ee++) {
139 for (int dd = 0; dd < 3; dd++) {
140 deltaChi[dd] = coordsEdges[6 * ee + dd] - coordsEdges[6 * ee + 3 + dd];
141 }
142 deltaX.clear();
143 for (int ii = 0; ii != 3; ++ii)
144 for (int jj = 0; jj != 3; ++jj)
145 deltaX(ii) += this->F(ii, jj) * deltaChi(jj);
146
147 for (int dd = 0; dd < 3; dd++) {
148 lrmsSquared += (1. / 6.) * deltaX[dd] * deltaX[dd];
149 lrmsSquared0 += (1. / 6.) * deltaChi[dd] * deltaChi[dd];
150 }
151 for (int ii = 0; ii != 3; ++ii)
152 for (int jj = 0; jj != 3; ++jj)
153 dXdChiT(ii, jj) += deltaX(ii) * deltaChi(jj);
154 }
156 }
157
158 /** \brief Calculate Q
159
160 \f[
161 \mathbf{Q} =
162 \mathbf{F}^{-\mathsf{T}}
163 -
164 \frac{1}{2}
165 \frac{1}{l^2_\textrm{rms}}
166 \sum_i^6
167 \Delta\mathbf{X}_i
168 \Delta\boldsymbol\chi_i^\mathsf{T}
169 \f]
170
171 */
174 Q.resize(3, 3, false);
178 t_Q(i, j) = t_invF(j, i) - 0.5 * t_dXdChiT(i, j) / lrmsSquared;
180 }
181
182 /** \brief Volume Length Quality
183
184 Based on:
185 Three‐dimensional brittle fracture: configurational‐force‐driven crack
186 propagation International Journal for Numerical Methods in Engineering 97
187 (7), 531-550
188
189 \f[
190 \mathcal{B}(a)=\frac{a}{2(1-\gamma)}-\ln{(a-\gamma)}
191 \f]
192
193 \f[
194 q = q_0 b,
195 \quad q_0 = 6\sqrt{2}\frac{V_0}{L^3_\textrm{rms}},
196 \quad b = \frac{\textrm{det}(\mathbf{F})}{\textrm{d}l^3_\textrm{rms}}
197 \f]
198
199 \f[
200 \mathbf{P} = \mathcal{B}(a)\mathbf{Q},
201 \f]
202 where \f$a\f$ depending on problem could be \f$q\f$ or \f$b\f$.
203
204 */
206 const NonlinearElasticElement::BlockData block_data,
207 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
209
212 if (this->invF.size1() != 3)
213 this->invF.resize(3, 3);
214
215 CHKERR invertTensor3by3(this->F, detF, this->invF);
218
219 double lrms03 = lrmsSquared0 * sqrt(lrmsSquared0);
220 b = detF / (lrmsSquared * sqrt(lrmsSquared) / lrms03);
221
222 currentVolume = detF * this->opPtr->getVolume();
223 q = 6. * sqrt(2.) * currentVolume / (lrmsSquared * sqrt(lrmsSquared));
224
225 if (this->P.size1() != 3)
226 this->P.resize(3, 3);
227
228 switch (tYpe) {
229 case QUALITY:
230 // Only use for testing, simple quality gradient
231 noalias(this->P) = q * Q;
232 break;
234 // This is used form mesh smoothing
235 tMp = q / (1.0 - gAmma) - 1.0 / (q - gAmma);
236 noalias(this->P) = tMp * Q;
237 break;
239 // Works well with Arbitrary Lagrangian Formulation
240 tMp = b / (1.0 - gAmma) - 1.0 / (b - gAmma);
241 noalias(this->P) = tMp * Q;
242 break;
244 // When scaled by volume works well with ALE and face flipping.
245 // Works well with smooth crack propagation
247 tMp *= b / (1.0 - gAmma) - 1.0 / (b - gAmma);
248 noalias(this->P) = tMp * Q;
249 break;
250 }
251
252 // Divide by volume, to make units as they should be
253 this->P *= aLpha / this->opPtr->getVolume();
254
256 }
257};
258
259#endif //__VOLUME_LENGTH_QUALITY_HPP__
VolumeLengthQualityType
@ LASTOP_VOLUMELENGTHQUALITYTYPE
@ BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME
@ BARRIER_AND_QUALITY
@ BARRIER_AND_CHANGE_QUALITY
static const char * VolumeLengthQualityTypeNames[]
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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 MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:972
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:1204
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
data for calculation heat conductivity and heat capacity elements
Implementation of elastic (non-linear) St. Kirchhoff equation.
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_invF
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator * opPtr
pointer to finite element tetrahedral operator
structure grouping operators and data used for calculation of nonlinear elastic element
Volume Length Quality.
MoFEMErrorCode calculateQ()
Calculate Q.
ublas::matrix< TYPE > Q
VolumeLengthQuality(VolumeLengthQualityType type, double alpha, double gamma)
FTensor::Index< 'i', 3 > i
MoFEMErrorCode getEdgesFromElemCoords()
MoFEMErrorCode getMaterialOptions()
FTensor::Index< 'j', 3 > j
MoFEMErrorCode calculateLrms()
Calculate mean element edge length.
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Volume Length Quality.
ublas::matrix< TYPE > dXdChiT
ublas::vector< TYPE > deltaX
FTensor::Index< 'k', 3 > k