v0.9.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 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifndef __VOLUME_LENGTH_QUALITY_HPP__
21 #define __VOLUME_LENGTH_QUALITY_HPP__
22 
23 #ifndef WITH_ADOL_C
24 #error "MoFEM need to be compiled with ADOL-C"
25 #endif
26 
33 };
34 
35 static const char *VolumeLengthQualityTypeNames[] = {
36  "QUALITY", "BARRIER_AND_QUALITY", "BARRIER_AND_CHANGE_QUALITY",
37  "BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME"};
38 
39 /** \brief Volume Length Quality
40  \ingroup nonlinear_elastic_elem
41 
42  */
43 template <typename TYPE>
46  TYPE> {
47 
48  // VolumeLengthQualityType tYpe;
49  int tYpe;
50  double aLpha;
51  double gAmma;
52 
53  VolumeLengthQuality(VolumeLengthQualityType type, double alpha, double gamma)
54  : NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<TYPE>(),
55  tYpe(type), aLpha(alpha), gAmma(gamma) {
57  CHKERRABORT(PETSC_COMM_WORLD, ierr);
58  }
59 
62  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
63  "Get VolumeLengthQuality material options",
64  "none");
65  CHKERR PetscOptionsEList(
66  "-volume_length_type", "Volume length quality type", "",
68  VolumeLengthQualityTypeNames[tYpe], &tYpe, PETSC_NULL);
69  CHKERR PetscOptionsScalar("-volume_length_alpha",
70  "volume length alpha parameter", "", aLpha,
71  &aLpha, PETSC_NULL);
72  CHKERR PetscOptionsScalar("-volume_length_gamma",
73  "volume length parameter (barrier)", "", gAmma,
74  &gAmma, PETSC_NULL);
75  ierr = PetscOptionsEnd();
77  }
78 
80  double lrmsSquared0;
82 
83  ublas::vector<TYPE> deltaX;
84  ublas::matrix<TYPE> Q, dXdChiT;
86 
87  /** Get coordinates of edges using cannonical element numeration
88  */
91  if (coordsEdges.empty()) {
92  coordsEdges.resize(6 * 2 * 3, false);
93  }
94  cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
95  &coordsEdges[0 * 3 * 2 + 0], 1);
96  cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
97  &coordsEdges[0 * 3 * 2 + 3], 1);
98  cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
99  &coordsEdges[1 * 3 * 2 + 0], 1);
100  cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
101  &coordsEdges[1 * 3 * 2 + 3], 1);
102  cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
103  &coordsEdges[2 * 3 * 2 + 0], 1);
104  cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
105  &coordsEdges[2 * 3 * 2 + 3], 1);
106  cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
107  &coordsEdges[3 * 3 * 2 + 0], 1);
108  cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
109  &coordsEdges[3 * 3 * 2 + 3], 1);
110  cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
111  &coordsEdges[4 * 3 * 2 + 0], 1);
112  cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
113  &coordsEdges[4 * 3 * 2 + 3], 1);
114  cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
115  &coordsEdges[5 * 3 * 2 + 0], 1);
116  cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
117  &coordsEdges[5 * 3 * 2 + 3], 1);
119  }
120 
121  /** \brief Calculate mean element edge length
122 
123  \f[
124  \Delta \boldsymbol\chi = \boldsymbol\chi^1 - \boldsymbol\chi^2
125  \f]
126 
127  \f[
128  \Delta X = \mathbf{F} \Delta \boldsymbol\chi
129  \f]
130 
131  \f[
132  l_\textrm{rms} = \sqrt{\frac{1}{6} \sum_{i=0}^6 l_i^2 } =
133  L_\textrm{rms}dl_\textrm{rms} \f]
134 
135  */
138  if (deltaChi.size() != 3) {
139  deltaChi.resize(3);
140  deltaX.resize(3);
141  dXdChiT.resize(3, 3);
142  }
143  lrmsSquared = 0;
144  lrmsSquared0 = 0;
145  dXdChiT.clear();
146  for (int ee = 0; ee < 6; ee++) {
147  for (int dd = 0; dd < 3; dd++) {
148  deltaChi[dd] = coordsEdges[6 * ee + dd] - coordsEdges[6 * ee + 3 + dd];
149  }
150  noalias(deltaX) = prod(this->F, deltaChi);
151  for (int dd = 0; dd < 3; dd++) {
152  lrmsSquared += (1. / 6.) * deltaX[dd] * deltaX[dd];
153  lrmsSquared0 += (1. / 6.) * deltaChi[dd] * deltaChi[dd];
154  }
155  noalias(dXdChiT) += outer_prod(deltaX, deltaChi);
156  }
158  }
159 
160  /** \brief Calculate Q
161 
162  \f[
163  \mathbf{Q} =
164  \mathbf{F}^{-\mathsf{T}}
165  -
166  \frac{1}{2}
167  \frac{1}{l^2_\textrm{rms}}
168  \sum_i^6
169  \Delta\mathbf{X}_i
170  \Delta\boldsymbol\chi_i^\mathsf{T}
171  \f]
172 
173  */
176  if (Q.size1() == 0) {
177  Q.resize(3, 3);
178  }
179  noalias(Q) = trans(this->invF) - 0.5 * dXdChiT / lrmsSquared;
181  }
182 
183  /** \brief Volume Length Quality
184 
185  Based on:
186  Three‐dimensional brittle fracture: configurational‐force‐driven crack
187  propagation International Journal for Numerical Methods in Engineering 97
188  (7), 531-550
189 
190  \f[
191  \mathcal{B}(a)=\frac{a}{2(1-\gamma)}-\ln{(a-\gamma)}
192  \f]
193 
194  \f[
195  q = q_0 b,
196  \quad q_0 = 6\sqrt{2}\frac{V_0}{L^3_\textrm{rms}},
197  \quad b = \frac{\textrm{det}(\mathbf{F})}{\textrm{d}l^3_\textrm{rms}}
198  \f]
199 
200  \f[
201  \mathbf{P} = \mathcal{B}(a)\mathbf{Q},
202  \f]
203  where \f$a\f$ depending on problem could be \f$q\f$ or \f$b\f$.
204 
205  */
207  const NonlinearElasticElement::BlockData block_data,
208  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
210 
212  CHKERR this->dEterminant(this->F, detF);
213  if (this->invF.size1() != 3)
214  this->invF.resize(3, 3);
215 
216  CHKERR this->iNvert(detF, this->F, this->invF);
218  CHKERR calculateQ();
219 
220  double lrms03 = lrmsSquared0 * sqrt(lrmsSquared0);
221  b = detF / (lrmsSquared * sqrt(lrmsSquared) / lrms03);
222 
223  currentVolume = detF * this->opPtr->getVolume();
224  q = 6. * sqrt(2.) * currentVolume / (lrmsSquared * sqrt(lrmsSquared));
225 
226  if (this->P.size1() != 3)
227  this->P.resize(3, 3);
228 
229  switch (tYpe) {
230  case QUALITY:
231  // Only use for testing, simple quality gradient
232  noalias(this->P) = q * Q;
233  break;
234  case BARRIER_AND_QUALITY:
235  // This is used form mesh smoothing
236  tMp = q / (1.0 - gAmma) - 1.0 / (q - gAmma);
237  noalias(this->P) = tMp * Q;
238  break;
240  // Works well with Arbitrary Lagrangian Formulation
241  tMp = b / (1.0 - gAmma) - 1.0 / (b - gAmma);
242  noalias(this->P) = tMp * Q;
243  break;
245  // When scaled by volume works well with ALE and face flipping.
246  // Works well with smooth crack propagation
247  tMp = currentVolume;
248  tMp *= b / (1.0 - gAmma) - 1.0 / (b - gAmma);
249  noalias(this->P) = tMp * Q;
250  break;
251  }
252 
253  // Divide by volume, to make units as they should be
254  this->P *= aLpha / this->opPtr->getVolume();
255 
257  }
258 };
259 
260 #endif //__VOLUME_LENGTH_QUALITY_HPP__
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > F
MoFEMErrorCode calculateLrms()
Calculate mean element edge length.
VolumeLengthQuality(VolumeLengthQualityType type, double alpha, double gamma)
MoFEMErrorCode getMaterialOptions()
MoFEMErrorCode dEterminant(ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &a, TYPE &det)
Calculate determinant of 3x3 matrix.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
Implementation of elastic (non-linear) St. Kirchhoff equation.
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > P
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
MoFEMErrorCode calculateQ()
Calculate Q.
ublas::matrix< TYPE > dXdChiT
Volume Length Quality.
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
VolumeLengthQualityType
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator * opPtr
pointer to finite element tetrahedral operator
ublas::matrix< TYPE > Q
ublas::vector< TYPE > deltaX
#define CHKERR
Inline error check.
Definition: definitions.h:601
MoFEMErrorCode getEdgesFromElemCoords()
MoFEMErrorCode iNvert(TYPE det, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &a, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &inv_a)
Calculate inverse of 3x3 matrix.
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Volume Length Quality.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
static const char * VolumeLengthQualityTypeNames[]
data for calculation het conductivity and heat capacity elements
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > invF
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...