v0.15.0
Loading...
Searching...
No Matches
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
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
41 FTensor::Index<'i', 3> i;
42 FTensor::Index<'j', 3> j;
43 FTensor::Index<'k', 3> k;
44
45 VolumeLengthQuality(VolumeLengthQualityType type, double alpha, double gamma)
47 tYpe(type), aLpha(alpha), gAmma(gamma) {
49 CHKERRABORT(PETSC_COMM_WORLD, ierr);
50 }
51
52 MoFEMErrorCode getMaterialOptions() {
54 PetscOptionsBegin(PETSC_COMM_WORLD, "",
55 "Get VolumeLengthQuality material options", "none");
56 CHKERR PetscOptionsEList(
57 "-volume_length_type", "Volume length quality type", "",
59 VolumeLengthQualityTypeNames[tYpe], &tYpe, PETSC_NULLPTR);
60 CHKERR PetscOptionsScalar("-volume_length_alpha",
61 "volume length alpha parameter", "", aLpha,
62 &aLpha, PETSC_NULLPTR);
63 CHKERR PetscOptionsScalar("-volume_length_gamma",
64 "volume length parameter (barrier)", "", gAmma,
65 &gAmma, PETSC_NULLPTR);
66 PetscOptionsEnd();
68 }
69
70 VectorDouble coordsEdges;
71 double lrmsSquared0;
72 VectorDouble deltaChi;
73
74 ublas::vector<TYPE> deltaX;
75 ublas::matrix<TYPE> Q, dXdChiT;
77
78 /** Get coordinates of edges using cannonical element numeration
79 */
80 MoFEMErrorCode getEdgesFromElemCoords() {
82 if (coordsEdges.empty()) {
83 coordsEdges.resize(6 * 2 * 3, false);
84 }
85 cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
86 &coordsEdges[0 * 3 * 2 + 0], 1);
87 cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
88 &coordsEdges[0 * 3 * 2 + 3], 1);
89 cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
90 &coordsEdges[1 * 3 * 2 + 0], 1);
91 cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
92 &coordsEdges[1 * 3 * 2 + 3], 1);
93 cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
94 &coordsEdges[2 * 3 * 2 + 0], 1);
95 cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
96 &coordsEdges[2 * 3 * 2 + 3], 1);
97 cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
98 &coordsEdges[3 * 3 * 2 + 0], 1);
99 cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
100 &coordsEdges[3 * 3 * 2 + 3], 1);
101 cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
102 &coordsEdges[4 * 3 * 2 + 0], 1);
103 cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
104 &coordsEdges[4 * 3 * 2 + 3], 1);
105 cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
106 &coordsEdges[5 * 3 * 2 + 0], 1);
107 cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
108 &coordsEdges[5 * 3 * 2 + 3], 1);
110 }
111
112 /** \brief Calculate mean element edge length
113
114 \f[
115 \Delta \boldsymbol\chi = \boldsymbol\chi^1 - \boldsymbol\chi^2
116 \f]
117
118 \f[
119 \Delta X = \mathbf{F} \Delta \boldsymbol\chi
120 \f]
121
122 \f[
123 l_\textrm{rms} = \sqrt{\frac{1}{6} \sum_{i=0}^6 l_i^2 } =
124 L_\textrm{rms}dl_\textrm{rms} \f]
125
126 */
127 MoFEMErrorCode calculateLrms() {
129 if (deltaChi.size() != 3) {
130 deltaChi.resize(3);
131 deltaX.resize(3);
132 dXdChiT.resize(3, 3);
133 }
134 lrmsSquared = 0;
135 lrmsSquared0 = 0;
136 std::fill(dXdChiT.data().begin(), dXdChiT.data().end(), 0);
137 for (int ee = 0; ee < 6; ee++) {
138 for (int dd = 0; dd < 3; dd++) {
139 deltaChi[dd] = coordsEdges[6 * ee + dd] - coordsEdges[6 * ee + 3 + dd];
140 }
141 std::fill(deltaX.begin(), deltaX.end(), 0);
142 for (int ii = 0; ii != 3; ++ii)
143 for (int jj = 0; jj != 3; ++jj)
144 deltaX(ii) += this->F(ii, jj) * deltaChi(jj);
145
146 for (int dd = 0; dd < 3; dd++) {
147 lrmsSquared += (1. / 6.) * deltaX[dd] * deltaX[dd];
148 lrmsSquared0 += (1. / 6.) * deltaChi[dd] * deltaChi[dd];
149 }
150 for (int ii = 0; ii != 3; ++ii)
151 for (int jj = 0; jj != 3; ++jj)
152 dXdChiT(ii, jj) += deltaX(ii) * deltaChi(jj);
153 }
155 }
156
157 /** \brief Calculate Q
158
159 \f[
160 \mathbf{Q} =
161 \mathbf{F}^{-\mathsf{T}}
162 -
163 \frac{1}{2}
164 \frac{1}{l^2_\textrm{rms}}
165 \sum_i^6
166 \Delta\mathbf{X}_i
167 \Delta\boldsymbol\chi_i^\mathsf{T}
168 \f]
169
170 */
171 MoFEMErrorCode calculateQ() {
173 Q.resize(3, 3, false);
174 auto t_Q = getFTensor2FromArray3by3(Q, FTensor::Number<0>(), 0);
175 auto t_invF = getFTensor2FromArray3by3(this->invF, FTensor::Number<0>(), 0);
176 auto t_dXdChiT = getFTensor2FromArray3by3(dXdChiT, FTensor::Number<0>(), 0);
177 t_Q(i, j) = t_invF(j, i) - 0.5 * t_dXdChiT(i, j) / lrmsSquared;
179 }
180
181 /** \brief Volume Length Quality
182
183 Based on:
184 Three‐dimensional brittle fracture: configurational‐force‐driven crack
185 propagation International Journal for Numerical Methods in Engineering 97
186 (7), 531-550
187
188 \f[
189 \mathcal{B}(a)=\frac{a}{2(1-\gamma)}-\ln{(a-\gamma)}
190 \f]
191
192 \f[
193 q = q_0 b,
194 \quad q_0 = 6\sqrt{2}\frac{V_0}{L^3_\textrm{rms}},
195 \quad b = \frac{\textrm{det}(\mathbf{F})}{\textrm{d}l^3_\textrm{rms}}
196 \f]
197
198 \f[
199 \mathbf{P} = \mathcal{B}(a)\mathbf{Q},
200 \f]
201 where \f$a\f$ depending on problem could be \f$q\f$ or \f$b\f$.
202
203 */
204 virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(
205 const NonlinearElasticElement::BlockData block_data,
206 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
208
210 detF = determinantTensor3by3(this->F);
211 if (this->invF.size1() != 3)
212 this->invF.resize(3, 3);
213
214 CHKERR invertTensor3by3(this->F, detF, this->invF);
217
218 double lrms03 = lrmsSquared0 * sqrt(lrmsSquared0);
219 b = detF / (lrmsSquared * sqrt(lrmsSquared) / lrms03);
220
221 currentVolume = detF * this->opPtr->getVolume();
222 q = 6. * sqrt(2.) * currentVolume / (lrmsSquared * sqrt(lrmsSquared));
223
224 if (this->P.size1() != 3)
225 this->P.resize(3, 3);
226
227 switch (tYpe) {
228 case QUALITY:
229 // Only use for testing, simple quality gradient
230 noalias(this->P) = q * Q;
231 break;
233 // This is used form mesh smoothing
234 tMp = q / (1.0 - gAmma) - 1.0 / (q - gAmma);
235 noalias(this->P) = tMp * Q;
236 break;
238 // Works well with Arbitrary Lagrangian Formulation
239 tMp = b / (1.0 - gAmma) - 1.0 / (b - gAmma);
240 noalias(this->P) = tMp * Q;
241 break;
243 // When scaled by volume works well with ALE and face flipping.
244 // Works well with smooth crack propagation
246 tMp *= b / (1.0 - gAmma) - 1.0 / (b - gAmma);
247 noalias(this->P) = tMp * Q;
248 break;
249 }
250
251 // Divide by volume, to make units as they should be
252 this->P *= aLpha / this->opPtr->getVolume();
253
255 }
256};
257
258#endif //__VOLUME_LENGTH_QUALITY_HPP__
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.
FTensor::Index< 'k', 3 > k
MoFEMErrorCode calculateQ()
Calculate Q.
ublas::matrix< TYPE > Q
ublas::vector< TYPE > deltaX
FTensor::Index< 'j', 3 > j
VolumeLengthQuality(VolumeLengthQualityType type, double alpha, double gamma)
FTensor::Index< 'i', 3 > i
MoFEMErrorCode getEdgesFromElemCoords()
MoFEMErrorCode getMaterialOptions()
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
VolumeLengthQualityType
@ LASTOP_VOLUMELENGTHQUALITYTYPE
@ BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME
@ BARRIER_AND_QUALITY
@ BARRIER_AND_CHANGE_QUALITY
static const char * VolumeLengthQualityTypeNames[]
VolumeLengthQualityType
@ LASTOP_VOLUMELENGTHQUALITYTYPE
@ BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME
@ BARRIER_AND_QUALITY
@ BARRIER_AND_CHANGE_QUALITY
static const char * VolumeLengthQualityTypeNames[]