v0.16.0
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | List of all members
HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp > Struct Template Reference

#include "tutorials/vec-2_nonlinear_elasticity/src/HenckyOps.hpp"

Inheritance diagram for HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >:
[legend]
Collaboration diagram for HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >:
[legend]

Public Member Functions

 OpCalculateEigenValsImpl (const std::string field_name, boost::shared_ptr< CommonData > common_data)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 
 OpCalculateEigenValsImpl (const std::string field_name, boost::shared_ptr< CommonData > common_data)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 

Private Attributes

boost::shared_ptr< CommonData > commonDataPtr
 

Detailed Description

template<int DIM, typename DomainEleOp>
struct HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >

Definition at line 83 of file HenckyOps.hpp.

Constructor & Destructor Documentation

◆ OpCalculateEigenValsImpl() [1/2]

template<int DIM, typename DomainEleOp >
HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::OpCalculateEigenValsImpl ( const std::string  field_name,
boost::shared_ptr< CommonData >  common_data 
)
inline

Definition at line 85 of file HenckyOps.hpp.

87 : DomainEleOp(field_name, DomainEleOp::OPROW),
88 commonDataPtr(common_data) {
89 std::fill(&DomainEleOp::doEntities[MBEDGE],
90 &DomainEleOp::doEntities[MBMAXTYPE], false);
91 }
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Definition child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13

◆ OpCalculateEigenValsImpl() [2/2]

template<int DIM, typename DomainEleOp >
HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::OpCalculateEigenValsImpl ( const std::string  field_name,
boost::shared_ptr< CommonData >  common_data 
)
inline

Definition at line 156 of file HenckyOps.hpp.

158 : DomainEleOp(field_name, DomainEleOp::OPROW),
159 commonDataPtr(common_data) {
160 std::fill(&DomainEleOp::doEntities[MBEDGE],
161 &DomainEleOp::doEntities[MBMAXTYPE], false);
162 }

Member Function Documentation

◆ doWork() [1/2]

template<int DIM, typename DomainEleOp >
MoFEMErrorCode HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::doWork ( int  side,
EntityType  type,
EntData &  data 
)
inline

Definition at line 93 of file HenckyOps.hpp.

93 {
94 MoFEMFunctionBegin;
95
96 FTensor::Index<'i', DIM> i;
97 FTensor::Index<'j', DIM> j;
98 FTensor::Index<'k', DIM> k;
99
100 [[maybe_unused]] constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
101 // const size_t nb_gauss_pts = matGradPtr->size2();
102 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
103
104 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
105
106 commonDataPtr->matEigVal.resize(nb_gauss_pts, DIM, false);
107 commonDataPtr->matEigVec.resize(nb_gauss_pts, DIM * DIM, false);
108 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
109 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
110
111 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
112
113 FTensor::Tensor2<double, DIM, DIM> F;
114 FTensor::Tensor2_symmetric<double, DIM> C;
115 FTensor::Tensor1<double, DIM> eig;
116 FTensor::Tensor2<double, DIM, DIM> eigen_vec;
117
118 F(i, j) = t_grad(i, j) + t_kd(i, j);
119 C(i, j) = F(k, i) ^ F(k, j);
120
121 for (int ii = 0; ii != DIM; ii++)
122 for (int jj = 0; jj != DIM; jj++)
123 eigen_vec(ii, jj) = C(ii, jj);
124
125 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
126 for (auto ii = 0; ii != DIM; ++ii)
127 eig(ii) = std::max(eps, eig(ii));
128
129 // rare case when two eigen values are equal
130 auto nb_uniq =
131 getUniqNb<DIM>(getVectorAdaptor(&eig(0), DIM), FTensor::Number<DIM>{});
132 if constexpr (DIM == 3) {
133 if (nb_uniq == 2) {
134 CHKERR sortEigenVals<DIM>(
135 getVectorAdaptor(&eig(0), DIM),
136 getMatrixAdaptor(&eigen_vec(0, 0), DIM, DIM),
137 FTensor::Number<DIM>{});
138 }
139 }
140
141 t_eig_val(i) = eig(i);
142 t_eig_vec(i, j) = eigen_vec(i, j);
143
144#ifndef NDEBUG
145 auto nb_uniq_test = getUniqNb<DIM>(getVectorAdaptor(&t_eig_val(0), DIM),
146 FTensor::Number<DIM>{});
147 if (nb_uniq_test != nb_uniq) {
148 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
149 "Inconsistent number of unique eigen values %ld != %ld",
150 nb_uniq, nb_uniq_test);
151 }
152#endif
153
154 ++t_grad;
155 ++t_eig_val;
156 ++t_eig_vec;
157 }
158
159 MoFEMFunctionReturn(0);
160 }
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
F
@ F
Definition free_surface.cpp:630
t_kd
constexpr auto t_kd
Definition free_surface.cpp:161
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
HenckyOps::eps
const double eps
Definition HenckyOps.hpp:13
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:49
MoFEM::getMatrixAdaptor
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition Templates.hpp:75
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition Templates.hpp:2118

◆ doWork() [2/2]

template<int DIM, typename DomainEleOp >
MoFEMErrorCode HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::doWork ( int  side,
EntityType  type,
EntData &  data 
)
inline

Definition at line 164 of file HenckyOps.hpp.

164 {
165 MoFEMFunctionBegin;
166
167 FTensor::Index<'i', DIM> i;
168 FTensor::Index<'j', DIM> j;
169 FTensor::Index<'k', DIM> k;
170
171 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
172 // const size_t nb_gauss_pts = matGradPtr->size2();
173 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
174
175 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
176
177 commonDataPtr->matEigVal.resize(DIM, nb_gauss_pts, false);
178 commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts, false);
179 auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
180 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
181
182 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
183
184 FTensor::Tensor2<double, DIM, DIM> F;
185 FTensor::Tensor2_symmetric<double, DIM> C;
186 FTensor::Tensor1<double, DIM> eig;
187 FTensor::Tensor2<double, DIM, DIM> eigen_vec;
188
189 F(i, j) = t_grad(i, j) + t_kd(i, j);
190 C(i, j) = F(k, i) ^ F(k, j);
191
192 for (int ii = 0; ii != DIM; ii++)
193 for (int jj = 0; jj != DIM; jj++)
194 eigen_vec(ii, jj) = C(ii, jj);
195
196 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
197 for (auto ii = 0; ii != DIM; ++ii)
198 eig(ii) = std::max(std::numeric_limits<double>::epsilon(), eig(ii));
199
200 // rare case when two eigen values are equal
201 auto nb_uniq = getUniqNb(eig);
202 if constexpr (DIM == 3) {
203 if (nb_uniq == 2) {
204 sortEigenVals(eig, eigen_vec);
205 }
206 }
207
208 t_eig_val(i) = eig(i);
209 t_eig_vec(i, j) = eigen_vec(i, j);
210
211 ++t_grad;
212 ++t_eig_val;
213 ++t_eig_vec;
214 }
215
216 MoFEMFunctionReturn(0);
217 }
MoFEM::sortEigenVals
auto sortEigenVals(VectorShallowArrayAdaptor< T > eig, MatrixShallowArrayAdaptor< T > eigen_vec, FTensor::Number< DIM >)
Definition Templates.hpp:2058
MoFEM::getUniqNb
auto getUniqNb(VectorShallowArrayAdaptor< T > vec, FTensor::Number< DIM >)
Get the Uniq Nb objects.
Definition Templates.hpp:2044

Member Data Documentation

◆ commonDataPtr

template<int DIM, typename DomainEleOp >
boost::shared_ptr< CommonData > HenckyOps::OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
private

Definition at line 163 of file HenckyOps.hpp.

The documentation for this struct was generated from the following files:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow