Volume Length Quality. More...

#include <users_modules/basic_finite_elements/nonlinear_elastic_materials/src/VolumeLengthQuality.hpp>

Inheritance diagram for VolumeLengthQuality< TYPE >:

Collaboration diagram for VolumeLengthQuality< TYPE >:

Public Member Functions
	VolumeLengthQuality (VolumeLengthQualityType type, double alpha, double gamma)

MoFEMErrorCode	getMaterialOptions ()

MoFEMErrorCode	getEdgesFromElemCoords ()

MoFEMErrorCode	calculateLrms ()
	Calculate mean element edge length. More...

MoFEMErrorCode	calculateQ ()
	Calculate Q. More...

virtual MoFEMErrorCode	calculateP_PiolaKirchhoffI (const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
	Volume Length Quality. More...

Public Member Functions inherited from NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< TYPE >
	FunctionsToCalculatePiolaKirchhoffI ()

virtual	~FunctionsToCalculatePiolaKirchhoffI ()=default

MoFEMErrorCode	calculateC_CauchyDeformationTensor ()

MoFEMErrorCode	calculateE_GreenStrain ()

MoFEMErrorCode	calculateS_PiolaKirchhoffII ()

virtual MoFEMErrorCode	calculateCauchyStress (const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
	Function overload to implement user material. More...

virtual MoFEMErrorCode	setUserActiveVariables (int &nb_active_variables)
	add additional active variables More...

virtual MoFEMErrorCode	setUserActiveVariables (VectorDouble &activeVariables)
	Add additional independent variables More complex physical models depend on gradient of defamation and some additional variables. For example can depend on temperature. This function adds additional independent variables to the model. More...

virtual MoFEMErrorCode	calculateElasticEnergy (const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
	Calculate elastic energy density. More...

virtual MoFEMErrorCode	calculatesIGma_EshelbyStress (const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
	Calculate Eshelby stress. More...

virtual MoFEMErrorCode	getDataOnPostProcessor (std::map< std::string, std::vector< VectorDouble >> &field_map, std::map< std::string, std::vector< MatrixDouble >> &grad_map)
	Do operations when pre-process. More...

Public Attributes
int	tYpe

double	aLpha

double	gAmma

FTensor::Index< 'i', 3 >	i

FTensor::Index< 'j', 3 >	j

FTensor::Index< 'k', 3 >	k

VectorDouble	coordsEdges

double	lrmsSquared0

VectorDouble	deltaChi

ublas::vector< TYPE >	deltaX

ublas::matrix< TYPE >	Q

ublas::matrix< TYPE >	dXdChiT

TYPE	lrmsSquared

TYPE	q

TYPE	b

TYPE	detF

TYPE	currentVolume

TYPE	tMp

Public Attributes inherited from NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< TYPE >
FTensor::Index< 'i', 3 >	i

FTensor::Index< 'j', 3 >	j

FTensor::Index< 'k', 3 >	k

double	lambda

double	mu

MatrixBoundedArray< TYPE, 9 >	F

MatrixBoundedArray< TYPE, 9 >	C

MatrixBoundedArray< TYPE, 9 >	E

MatrixBoundedArray< TYPE, 9 >	S

MatrixBoundedArray< TYPE, 9 >	invF

MatrixBoundedArray< TYPE, 9 >	P

MatrixBoundedArray< TYPE, 9 >	sIGma

MatrixBoundedArray< TYPE, 9 >	h

MatrixBoundedArray< TYPE, 9 >	H

MatrixBoundedArray< TYPE, 9 >	invH

MatrixBoundedArray< TYPE, 9 >	sigmaCauchy

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_F

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_C

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_E

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_S

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_invF

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_P

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_sIGma

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_h

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_H

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_invH

FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 >	t_sigmaCauchy

TYPE	J

TYPE	eNergy

TYPE	detH

TYPE	detF

int	gG
	Gauss point number. More...

CommonData *	commonDataPtr

MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator *	opPtr
	pointer to finite element tetrahedral operator More...

Additional Inherited Members
Static Protected Member Functions inherited from NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< TYPE >
static auto	resizeAndSet (MatrixBoundedArray< TYPE, 9 > &m)

Detailed Description

template<typename TYPE>
struct VolumeLengthQuality< TYPE >

Volume Length Quality.

Examples: mesh_smoothing.cpp.

Definition at line 32 of file VolumeLengthQuality.hpp.

Constructor & Destructor Documentation

◆ VolumeLengthQuality()

template<typename TYPE >

VolumeLengthQuality< TYPE >::VolumeLengthQuality	(	VolumeLengthQualityType	type,
		double	alpha,
		double	gamma
	)

inline

Definition at line 45 of file VolumeLengthQuality.hpp.

       : NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<TYPE>(),
         tYpe(type), aLpha(alpha), gAmma(gamma) {
     ierr = getMaterialOptions();
     CHKERRABORT(PETSC_COMM_WORLD, ierr);
   }

Member Function Documentation

◆ calculateLrms()

template<typename TYPE >

MoFEMErrorCode VolumeLengthQuality< TYPE >::calculateLrms ( )

inline

Calculate mean element edge length.

\[ \Delta \boldsymbol\chi = \boldsymbol\chi^1 - \boldsymbol\chi^2 \]

\[ \Delta X = \mathbf{F} \Delta \boldsymbol\chi \]

\[ l_\textrm{rms} = \sqrt{\frac{1}{6} \sum_{i=0}^6 l_i^2 } = L_\textrm{rms}dl_\textrm{rms} \]

Definition at line 128 of file VolumeLengthQuality.hpp.

                                  {
     MoFEMFunctionBegin;
     if (deltaChi.size() != 3) {
       deltaChi.resize(3);
       deltaX.resize(3);
       dXdChiT.resize(3, 3);
     }
     lrmsSquared = 0;
     lrmsSquared0 = 0;
     std::fill(dXdChiT.data().begin(), dXdChiT.data().end(), 0);
     for (int ee = 0; ee < 6; ee++) {
       for (int dd = 0; dd < 3; dd++) {
         deltaChi[dd] = coordsEdges[6 * ee + dd] - coordsEdges[6 * ee + 3 + dd];
       }
       std::fill(deltaX.begin(), deltaX.end(), 0);
       for (int ii = 0; ii != 3; ++ii)
         for (int jj = 0; jj != 3; ++jj)
           deltaX(ii) += this->F(ii, jj) * deltaChi(jj);
  
       for (int dd = 0; dd < 3; dd++) {
         lrmsSquared += (1. / 6.) * deltaX[dd] * deltaX[dd];
         lrmsSquared0 += (1. / 6.) * deltaChi[dd] * deltaChi[dd];
       }
       for (int ii = 0; ii != 3; ++ii)
         for (int jj = 0; jj != 3; ++jj)
           dXdChiT(ii, jj) += deltaX(ii) * deltaChi(jj);
     }
     MoFEMFunctionReturn(0);
   }

◆ calculateP_PiolaKirchhoffI()

template<typename TYPE >

virtual MoFEMErrorCode VolumeLengthQuality< TYPE >::calculateP_PiolaKirchhoffI	(	const NonlinearElasticElement::BlockData	block_data,
		boost::shared_ptr< const NumeredEntFiniteElement >	fe_ptr
	)

inlinevirtual

Volume Length Quality.

Based on: Three‐dimensional brittle fracture: configurational‐force‐driven crack propagation International Journal for Numerical Methods in Engineering 97 (7), 531-550

\[ \mathcal{B}(a)=\frac{a}{2(1-\gamma)}-\ln{(a-\gamma)} \]

\[ q = q_0 b, \quad q_0 = 6\sqrt{2}\frac{V_0}{L^3_\textrm{rms}}, \quad b = \frac{\textrm{det}(\mathbf{F})}{\textrm{d}l^3_\textrm{rms}} \]

\[ \mathbf{P} = \mathcal{B}(a)\mathbf{Q}, \]

where \(a\) depending on problem could be \(q\) or \(b\).

Reimplemented from NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< TYPE >.

Definition at line 205 of file VolumeLengthQuality.hpp.

                                                              {
     MoFEMFunctionBegin;
  
     CHKERR getEdgesFromElemCoords();
     detF = determinantTensor3by3(this->F);
     if (this->invF.size1() != 3)
       this->invF.resize(3, 3);
  
     CHKERR invertTensor3by3(this->F, detF, this->invF);
     CHKERR calculateLrms();
     CHKERR calculateQ();
  
     double lrms03 = lrmsSquared0 * sqrt(lrmsSquared0);
     b = detF / (lrmsSquared * sqrt(lrmsSquared) / lrms03);
  
     currentVolume = detF * this->opPtr->getVolume();
     q = 6. * sqrt(2.) * currentVolume / (lrmsSquared * sqrt(lrmsSquared));
  
     if (this->P.size1() != 3)
       this->P.resize(3, 3);
  
     switch (tYpe) {
     case QUALITY:
       // Only use for testing, simple quality gradient
       noalias(this->P) = q * Q;
       break;
     case BARRIER_AND_QUALITY:
       // This is used form mesh smoothing
       tMp = q / (1.0 - gAmma) - 1.0 / (q - gAmma);
       noalias(this->P) = tMp * Q;
       break;
     case BARRIER_AND_CHANGE_QUALITY:
       // Works well with Arbitrary Lagrangian Formulation
       tMp = b / (1.0 - gAmma) - 1.0 / (b - gAmma);
       noalias(this->P) = tMp * Q;
       break;
     case BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME:
       // When scaled by volume works well with ALE and face flipping.
       // Works well with smooth crack propagation
       tMp = currentVolume;
       tMp *= b / (1.0 - gAmma) - 1.0 / (b - gAmma);
       noalias(this->P) = tMp * Q;
       break;
     }
  
     // Divide by volume, to make units as they should be
     this->P *= aLpha / this->opPtr->getVolume();
  
     MoFEMFunctionReturn(0);
   }

◆ calculateQ()

template<typename TYPE >

MoFEMErrorCode VolumeLengthQuality< TYPE >::calculateQ ( )

inline

Calculate Q.

\[ \mathbf{Q} = \mathbf{F}^{-\mathsf{T}} - \frac{1}{2} \frac{1}{l^2_\textrm{rms}} \sum_i^6 \Delta\mathbf{X}_i \Delta\boldsymbol\chi_i^\mathsf{T} \]

Definition at line 172 of file VolumeLengthQuality.hpp.

                               {
     MoFEMFunctionBegin;
     Q.resize(3, 3, false);
     auto t_Q = getFTensor2FromArray3by3(Q, FTensor::Number<0>(), 0);
     auto t_invF = getFTensor2FromArray3by3(this->invF, FTensor::Number<0>(), 0);
     auto t_dXdChiT = getFTensor2FromArray3by3(dXdChiT, FTensor::Number<0>(), 0);
     t_Q(i, j) = t_invF(j, i) - 0.5 * t_dXdChiT(i, j) / lrmsSquared;
     MoFEMFunctionReturn(0);
   }

◆ getEdgesFromElemCoords()

template<typename TYPE >

MoFEMErrorCode VolumeLengthQuality< TYPE >::getEdgesFromElemCoords ( )

inline

Get coordinates of edges using cannonical element numeration

Definition at line 81 of file VolumeLengthQuality.hpp.

                                           {
     MoFEMFunctionBegin;
     if (coordsEdges.empty()) {
       coordsEdges.resize(6 * 2 * 3, false);
     }
     cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
                 &coordsEdges[0 * 3 * 2 + 0], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
                 &coordsEdges[0 * 3 * 2 + 3], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
                 &coordsEdges[1 * 3 * 2 + 0], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
                 &coordsEdges[1 * 3 * 2 + 3], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[0 * 3], 1,
                 &coordsEdges[2 * 3 * 2 + 0], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
                 &coordsEdges[2 * 3 * 2 + 3], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
                 &coordsEdges[3 * 3 * 2 + 0], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
                 &coordsEdges[3 * 3 * 2 + 3], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[1 * 3], 1,
                 &coordsEdges[4 * 3 * 2 + 0], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
                 &coordsEdges[4 * 3 * 2 + 3], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[2 * 3], 1,
                 &coordsEdges[5 * 3 * 2 + 0], 1);
     cblas_dcopy(3, &this->opPtr->getCoords()[3 * 3], 1,
                 &coordsEdges[5 * 3 * 2 + 3], 1);
     MoFEMFunctionReturn(0);
   }

◆ getMaterialOptions()

template<typename TYPE >

MoFEMErrorCode VolumeLengthQuality< TYPE >::getMaterialOptions ( )

inline

Definition at line 52 of file VolumeLengthQuality.hpp.

                                       {
     MoFEMFunctionBeginHot;
     CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
                              "Get VolumeLengthQuality material options",
                              "none");
     CHKERR PetscOptionsEList(
         "-volume_length_type", "Volume length quality type", "",
         VolumeLengthQualityTypeNames, LASTOP_VOLUMELENGTHQUALITYTYPE,
         VolumeLengthQualityTypeNames[tYpe], &tYpe, PETSC_NULL);
     CHKERR PetscOptionsScalar("-volume_length_alpha",
                               "volume length alpha parameter", "", aLpha,
                               &aLpha, PETSC_NULL);
     CHKERR PetscOptionsScalar("-volume_length_gamma",
                               "volume length parameter (barrier)", "", gAmma,
                               &gAmma, PETSC_NULL);
     ierr = PetscOptionsEnd();
     MoFEMFunctionReturnHot(0);
   }