143 {
145
149
151
152 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
153
154 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(
commonDataPtr->matGradPtr));
155
157 commonDataPtr->matEigVec.resize(nb_gauss_pts, DIM * DIM,
false);
158 auto t_eig_val = getFTensor1FromMat<DIM>(
commonDataPtr->matEigVal);
159 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(
commonDataPtr->matEigVec);
160
161 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
162
167
170
171 for (int ii = 0; ii != DIM; ii++)
172 for (int jj = 0; jj != DIM; jj++)
173 eigen_vec(ii, jj) = C(ii, jj);
174
176 for (auto ii = 0; ii != DIM; ++ii)
177 eig(ii) = std::max(
eps, eig(ii));
178
179
180 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
181 if constexpr (DIM == 3) {
182 if (nb_uniq == 2) {
183 sort_eigen_vals<DIM>(eig, eigen_vec);
184 }
185 }
186
187 t_eig_val(
i) = eig(
i);
188 t_eig_vec(
i,
j) = eigen_vec(
i,
j);
189
190#ifndef NDEBUG
191 auto nb_uniq_test = get_uniq_nb<DIM>(&t_eig_val(0));
192 if (nb_uniq_test != nb_uniq) {
194 "Inconsistent number of unique eigen values %ld != %ld",
195 nb_uniq, nb_uniq_test);
196 }
197#endif
198
199 ++t_grad;
200 ++t_eig_val;
201 ++t_eig_vec;
202 }
203
205 }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev