v0.14.0
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 
23 static 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  */
31 template <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 
54  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
55  "Get VolumeLengthQuality material options",
56  "none");
57  CHKERR PetscOptionsEList(
58  "-volume_length_type", "Volume length quality type", "",
60  VolumeLengthQualityTypeNames[tYpe], &tYpe, PETSC_NULL);
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 
72  double lrmsSquared0;
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  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  */
174  Q.resize(3, 3, false);
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  */
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);
217  CHKERR calculateQ();
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;
233  case BARRIER_AND_QUALITY:
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
246  tMp = currentVolume;
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__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
VolumeLengthQuality::aLpha
double aLpha
Definition: VolumeLengthQuality.hpp:38
VolumeLengthQuality::deltaX
ublas::vector< TYPE > deltaX
Definition: VolumeLengthQuality.hpp:75
VolumeLengthQuality::i
FTensor::Index< 'i', 3 > i
Definition: VolumeLengthQuality.hpp:41
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI::FunctionsToCalculatePiolaKirchhoffI
FunctionsToCalculatePiolaKirchhoffI()
Definition: NonLinearElasticElement.hpp:136
QUALITY
@ QUALITY
Definition: VolumeLengthQuality.hpp:16
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI::opPtr
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator * opPtr
pointer to finite element tetrahedral operator
Definition: NonLinearElasticElement.hpp:158
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI
Implementation of elastic (non-linear) St. Kirchhoff equation.
Definition: NonLinearElasticElement.hpp:79
VolumeLengthQuality::coordsEdges
VectorDouble coordsEdges
Definition: VolumeLengthQuality.hpp:71
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VolumeLengthQuality::gAmma
double gAmma
Definition: VolumeLengthQuality.hpp:39
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getCoords
VectorDouble & getCoords()
nodal coordinates
Definition: VolumeElementForcesAndSourcesCore.hpp:180
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
VolumeLengthQuality::tMp
TYPE tMp
Definition: VolumeLengthQuality.hpp:77
VolumeLengthQualityType
VolumeLengthQualityType
Definition: VolumeLengthQuality.hpp:15
VolumeLengthQuality::VolumeLengthQuality
VolumeLengthQuality(VolumeLengthQualityType type, double alpha, double gamma)
Definition: VolumeLengthQuality.hpp:45
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:1588
VolumeLengthQuality::lrmsSquared
TYPE lrmsSquared
Definition: VolumeLengthQuality.hpp:77
LASTOP_VOLUMELENGTHQUALITYTYPE
@ LASTOP_VOLUMELENGTHQUALITYTYPE
Definition: VolumeLengthQuality.hpp:20
BARRIER_AND_CHANGE_QUALITY
@ BARRIER_AND_CHANGE_QUALITY
Definition: VolumeLengthQuality.hpp:18
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI::F
MatrixBoundedArray< TYPE, 9 > F
Definition: NonLinearElasticElement.hpp:147
convert.type
type
Definition: convert.py:64
VolumeLengthQuality::getEdgesFromElemCoords
MoFEMErrorCode getEdgesFromElemCoords()
Definition: VolumeLengthQuality.hpp:81
VolumeLengthQuality::detF
TYPE detF
Definition: VolumeLengthQuality.hpp:77
VolumeLengthQuality
Volume Length Quality.
Definition: VolumeLengthQuality.hpp:32
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
VolumeLengthQuality::Q
ublas::matrix< TYPE > Q
Definition: VolumeLengthQuality.hpp:76
FTensor::Index< 'i', 3 >
VolumeLengthQuality::lrmsSquared0
double lrmsSquared0
Definition: VolumeLengthQuality.hpp:72
VolumeLengthQuality::k
FTensor::Index< 'k', 3 > k
Definition: VolumeLengthQuality.hpp:43
VolumeLengthQuality::calculateP_PiolaKirchhoffI
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Volume Length Quality.
Definition: VolumeLengthQuality.hpp:205
MoFEM::getFTensor2FromArray3by3
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1342
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
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
VolumeLengthQuality::tYpe
int tYpe
Definition: VolumeLengthQuality.hpp:37
VolumeLengthQuality::calculateQ
MoFEMErrorCode calculateQ()
Calculate Q.
Definition: VolumeLengthQuality.hpp:172
VolumeLengthQuality::getMaterialOptions
MoFEMErrorCode getMaterialOptions()
Definition: VolumeLengthQuality.hpp:52
VolumeLengthQuality::b
TYPE b
Definition: VolumeLengthQuality.hpp:77
VolumeLengthQuality::dXdChiT
ublas::matrix< TYPE > dXdChiT
Definition: VolumeLengthQuality.hpp:76
VolumeLengthQuality::deltaChi
VectorDouble deltaChi
Definition: VolumeLengthQuality.hpp:73
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
VolumeLengthQuality::j
FTensor::Index< 'j', 3 > j
Definition: VolumeLengthQuality.hpp:42
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
VolumeLengthQualityTypeNames
static const char * VolumeLengthQualityTypeNames[]
Definition: VolumeLengthQuality.hpp:23
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
BARRIER_AND_QUALITY
@ BARRIER_AND_QUALITY
Definition: VolumeLengthQuality.hpp:17
VolumeLengthQuality::calculateLrms
MoFEMErrorCode calculateLrms()
Calculate mean element edge length.
Definition: VolumeLengthQuality.hpp:128
VolumeLengthQuality::q
TYPE q
Definition: VolumeLengthQuality.hpp:77
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI::invF
MatrixBoundedArray< TYPE, 9 > invF
Definition: NonLinearElasticElement.hpp:147
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME
@ BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME
Definition: VolumeLengthQuality.hpp:19
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI::P
MatrixBoundedArray< TYPE, 9 > P
Definition: NonLinearElasticElement.hpp:147
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI::t_invF
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_invF
Definition: NonLinearElasticElement.hpp:150
VolumeLengthQuality::currentVolume
TYPE currentVolume
Definition: VolumeLengthQuality.hpp:77