3constexpr double eps = std::numeric_limits<float>::epsilon();
5auto f = [](
double v) {
return 0.5 * std::log(
static_cast<long double>(
v)); };
6auto d_f = [](
double v) {
return 0.5 /
static_cast<long double>(
v); };
8 return -0.5 / (
static_cast<long double>(
v) *
static_cast<long double>(
v));
11inline auto is_eq(
const double &
a,
const double &b) {
12 return std::abs(
a - b) <
eps;
16 std::array<double, DIM> tmp;
17 std::copy(ptr, &ptr[DIM], tmp.begin());
18 std::sort(tmp.begin(), tmp.end());
19 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(),
is_eq));
25 std::array<std::pair<double, size_t>, DIM> tmp;
26 auto is_eq_pair = [](
const auto &
a,
const auto &b) {
27 return a.first < b.first;
30 for (
size_t i = 0;
i != DIM; ++
i)
31 tmp[
i] = std::make_pair(eig(
i),
i);
32 std::sort(tmp.begin(), tmp.end(), is_eq_pair);
35 if (is_eq_pair(tmp[0], tmp[1])) {
46 eigen_vec(
i, 0), eigen_vec(
i, 1), eigen_vec(
i, 2),
48 eigen_vec(
j, 0), eigen_vec(
j, 1), eigen_vec(
j, 2),
50 eigen_vec(
k, 0), eigen_vec(
k, 1), eigen_vec(
k, 2)};
58 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
62struct CommonData :
public boost::enable_shared_from_this<CommonData> {
77 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
82 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
87 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
matLogC);
91 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
matTangent);
98 boost::shared_ptr<CommonData> common_data)
115 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(
commonDataPtr->matGradPtr));
118 commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts,
false);
119 auto t_eig_val = getFTensor1FromMat<DIM>(
commonDataPtr->matEigVal);
120 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(
commonDataPtr->matEigVec);
122 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
130 C(
i,
j) = F(
k,
i) ^ F(
k,
j);
132 for (
int ii = 0; ii != DIM; ii++)
133 for (
int jj = 0; jj != DIM; jj++)
134 eigen_vec(ii, jj) = C(ii, jj);
136 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
139 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
140 if (DIM == 3 && nb_uniq == 2)
141 sort_eigen_vals<DIM>(eig, eigen_vec);
143 t_eig_val(
i) = eig(
i);
144 t_eig_vec(
i,
j) = eigen_vec(
i,
j);
161 boost::shared_ptr<CommonData> common_data)
176 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
179 auto t_eig_val = getFTensor1FromMat<DIM>(
commonDataPtr->matEigVal);
180 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(
commonDataPtr->matEigVec);
182 auto t_logC = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matLogC);
184 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
187 auto nb_uniq = get_uniq_nb<DIM>(&(t_eig_val(0)));
191 eig(
i) = t_eig_val(
i);
192 eigen_vec(
i,
j) = t_eig_vec(
i,
j);
194 t_logC(
i,
j) = logC(
i,
j);
211 boost::shared_ptr<CommonData> common_data)
227 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
229 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->matLogCdC);
230 auto t_eig_val = getFTensor1FromMat<DIM>(
commonDataPtr->matEigVal);
231 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(
commonDataPtr->matEigVec);
233 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
237 eig(
i) = t_eig_val(
i);
238 eigen_vec(
i,
j) = t_eig_vec(
i,
j);
241 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
243 dlogC_dC(
i,
j,
k,
l) *= 2;
245 t_logC_dC(
i,
j,
k,
l) = dlogC_dC(
i,
j,
k,
l);
262 boost::shared_ptr<CommonData> common_data)
280 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
commonDataPtr->matDPtr);
281 auto t_logC = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matLogC);
282 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
284 auto t_T = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matHenckyStress);
286 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
287 t_T(
i,
j) = t_D(
i,
j,
k,
l) * t_logC(
k,
l);
303 boost::shared_ptr<CommonData> common_data,
304 boost::shared_ptr<MatrixDouble> mat_D_ptr,
305 const double scale = 1)
325 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
matDPtr);
326 auto t_logC = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matLogC);
327 auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*
matLogCPlastic);
328 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
330 auto t_T = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matHenckyStress);
332 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
333 t_T(
i,
j) = t_D(
i,
j,
k,
l) * (t_logC(
k,
l) - t_logCPlastic(
k,
l));
354 boost::shared_ptr<CommonData> common_data)
372 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
commonDataPtr->matDPtr);
373 auto t_logC = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matLogC);
374 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->matLogCdC);
375 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
376 commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts,
false);
378 auto t_P = getFTensor2FromMat<DIM, DIM>(
commonDataPtr->matFirstPiolaStress);
379 auto t_T = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matHenckyStress);
381 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matSecondPiolaStress);
382 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(
commonDataPtr->matGradPtr));
384 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
387 t_S(
k,
l) = t_T(
i,
j) * t_logC_dC(
i,
j,
k,
l);
388 t_P(
i,
l) = t_F(
i,
k) * t_S(
k,
l);
408 boost::shared_ptr<CommonData> common_data,
409 boost::shared_ptr<MatrixDouble> mat_D_ptr =
nullptr)
434 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
435 commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
437 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(
commonDataPtr->matTangent);
439 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
matDPtr);
440 auto t_eig_val = getFTensor1FromMat<DIM>(
commonDataPtr->matEigVal);
441 auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(
commonDataPtr->matEigVec);
442 auto t_T = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matHenckyStress);
444 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->matSecondPiolaStress);
445 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(
commonDataPtr->matGradPtr));
446 auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->matLogCdC);
448 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
456 eig(
i) = t_eig_val(
i);
457 eigen_vec(
i,
j) = t_eig_vec(
i,
j);
461 auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
471 P_D_P_plus_TL(
i,
j,
k,
l) =
473 (t_logC_dC(
i,
j,
o,
p) * t_D(
o,
p,
m,
n)) * t_logC_dC(
m,
n,
k,
l);
474 P_D_P_plus_TL(
i,
j,
k,
l) *= 0.5;
477 t_F(
i,
k) * (P_D_P_plus_TL(
k,
j,
o,
p) * dC_dF(
o,
p,
m,
n));
500 moab::Interface &post_proc_mesh,
501 std::vector<EntityHandle> &map_gauss_pts,
502 boost::shared_ptr<CommonData> common_data_ptr);
513 const std::string
field_name, moab::Interface &post_proc_mesh,
514 std::vector<EntityHandle> &map_gauss_pts,
515 boost::shared_ptr<CommonData> common_data_ptr)
517 mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
532 std::array<double, 9> def;
533 std::fill(def.begin(), def.end(), 0);
535 auto get_tag = [&](
const std::string name,
size_t size) {
537 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
538 MB_TAG_CREAT | MB_TAG_SPARSE,
543 MatrixDouble3by3 mat(3, 3);
545 auto set_matrix = [&](
auto &
t) -> MatrixDouble3by3 & {
553 auto set_tag = [&](
auto th,
auto gg, MatrixDouble3by3 &mat) {
554 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
555 &*mat.data().begin());
558 auto th_grad = get_tag(
"GRAD", 9);
559 auto th_stress = get_tag(
"STRESS", 9);
562 getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
563 auto t_grad = getFTensor2FromMat<DIM, DIM>(*commonDataPtr->matGradPtr);
568 for (
int gg = 0; gg != commonDataPtr->matGradPtr->size2(); ++gg) {
570 F(
i,
j) = t_grad(
i,
j) + kronecker_delta(
i,
j);
574 determinantTensor2by2(F, inv_t_detF);
576 determinantTensor3by3(F, inv_t_detF);
578 inv_t_detF = 1. / inv_t_detF;
580 cauchy_stress(
i,
j) = inv_t_detF * (t_piola(
i,
k) ^ F(
j,
k));
582 CHKERR set_tag(th_grad, gg, set_matrix(t_grad));
583 CHKERR set_tag(th_stress, gg, set_matrix(cauchy_stress));
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
auto get_uniq_nb(double *ptr)
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
auto is_eq(const double &a, const double &b)
constexpr double t
plate stiffness
constexpr auto field_name
auto getMatHenckyStress()
auto getMatFirstPiolaStress()
boost::shared_ptr< MatrixDouble > matLogCPlastic
MatrixDouble matFirstPiolaStress
boost::shared_ptr< MatrixDouble > matDPtr
MatrixDouble matSecondPiolaStress
MatrixDouble matHenckyStress
boost::shared_ptr< MatrixDouble > matGradPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateEigenVals(const std::string field_name, boost::shared_ptr< CommonData > common_data)
OpCalculateHenckyPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr, const double scale=1)
boost::shared_ptr< MatrixDouble > matDPtr
boost::shared_ptr< MatrixDouble > matLogCPlastic
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHenckyStress(const std::string field_name, boost::shared_ptr< CommonData > common_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateLogC_dC(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateLogC(const std::string field_name, boost::shared_ptr< CommonData > common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
OpCalculatePiolaStress(const std::string field_name, boost::shared_ptr< CommonData > common_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > matDPtr
OpHenckyTangent(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
OpPostProcHencky(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
boost::shared_ptr< CommonData > commonDataPtr
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element