v0.14.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
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 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
71 VectorDouble coordsEdges;
73 VectorDouble deltaChi;
74
75 ublas::vector<TYPE> deltaX;
76 ublas::matrix<TYPE> Q, dXdChiT;
78
79 /** Get coordinates of edges using cannonical element numeration
80 */
81 MoFEMErrorCode getEdgesFromElemCoords() {
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 */
128 MoFEMErrorCode calculateLrms() {
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 std::fill(dXdChiT.data().begin(), dXdChiT.data().end(), 0);
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 std::fill(deltaX.begin(), deltaX.end(), 0);
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 */
172 MoFEMErrorCode calculateQ() {
174 Q.resize(3, 3, false);
175 auto t_Q = getFTensor2FromArray3by3(Q, FTensor::Number<0>(), 0);
176 auto t_invF = getFTensor2FromArray3by3(this->invF, FTensor::Number<0>(), 0);
177 auto t_dXdChiT = getFTensor2FromArray3by3(dXdChiT, FTensor::Number<0>(), 0);
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 */
205 virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(
206 const NonlinearElasticElement::BlockData block_data,
207 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
209
211 detF = determinantTensor3by3(this->F);
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[]
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.
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