Implementation of operators.
#include <boost/math/constants/constants.hpp>
double f(
const double v) {
return exp(
v); }
double d_f(
const double v) {
return exp(
v); }
double dd_f(
const double v) {
return exp(
v); }
template <class T> struct TensorTypeExtractor {
typedef typename std::remove_pointer<T>::type
Type;
};
template <
class T,
int I>
struct TensorTypeExtractor<
FTensor::PackPtr<T, I>> {
typedef typename std::remove_pointer<T>::type
Type;
};
template <typename T>
&t_omega) {
sqrt(t_omega(
i) * t_omega(
i));
constexpr double tol = 1e-18;
if (std::abs(angle) <
tol) {
return t_R;
}
t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
2. * ss_2 * ss_2 / (angle * angle);
t_R(
i,
j) +=
a * t_Omega(
i,
j);
t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
return t_R;
}
template <typename T>
&t_omega) {
sqrt(t_omega(
i) * t_omega(
i));
constexpr double tol = 1e-18;
if (std::abs(angle) <
tol) {
t_diff_R(
i,
j,
k) = FTensor::levi_civita<double>(
i,
l,
k) * t_R(
l,
j);
return t_diff_R;
}
t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
t_diff_R(
i,
j,
k) =
a * FTensor::levi_civita<double>(
i,
j,
k);
b * (t_Omega(
i,
l) * FTensor::levi_civita<double>(
l,
j,
k) +
FTensor::levi_civita<double>(
i,
l,
k) * t_Omega(
l,
j));
(angle * cc - ss) / angle2;
(2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
t_diff_R(
i,
j,
k) += diff_a * t_Omega(
i,
j) * (t_omega(
k) / angle);
diff_b * t_Omega(
i,
l) * t_Omega(
l,
j) * (t_omega(
k) / angle);
return t_diff_R;
}
template <typename T1, typename T2, typename T3>
T3 &eig_vec) {
for (int ii = 0; ii != 3; ii++)
for (int jj = 0; jj != 3; jj++)
eig_vec(ii, jj) = X(ii, jj);
int lda = 3;
int lwork = (3 + 2) * 3;
std::array<
double, (3 + 2) * 3> work;
if (
lapack_dsyev(
'V',
'U',
n, &(eig_vec(0, 0)), lda, &eig(0), work.data(),
lwork) > 0)
"The algorithm failed to compute eigenvalues.");
}
inline bool is_eq(
const double &
a,
const double &b) {
constexpr double eps = 1e-8;
}
return distance(&ec(0), unique(&ec(0), &ec(0) + 3,
is_eq));
}
template <typename T1, typename T2>
if (
is_eq(eig(0), eig(1))) {
eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
} else {
eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
}
}
if (type != MBTET)
int nb_integration_pts = data.
getN().size1();
dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts,
false);
dataAtPts->streachTensorAtPts.resize(6, nb_integration_pts,
false);
dataAtPts->diffStreachTensorAtPts.resize(36, nb_integration_pts,
false);
dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
dataAtPts->nbUniq.resize(nb_integration_pts,
false);
auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_log_u =
getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachTensorAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_diff_R(
i,
j,
k) = t0_diff(
i,
j,
k);
t_eigen_vals(
i) = eig(
i);
t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
if (nbUniq[gg] == 2)
auto t_diff_u_tmp =
for (int ii = 0; ii != 3; ++ii) {
for (int jj = ii; jj != 3; ++jj) {
for (int kk = 0; kk != 3; ++kk) {
for (int ll = 0; ll < kk; ++ll) {
t_diff_u_tmp(ii, jj, kk, ll) *= 2;
}
}
}
}
t_u(
i,
j) = t_u_tmp(
i,
j);
t_diff_u(
i,
j,
k,
l) = t_diff_u_tmp(
i,
j,
k,
l);
t_h(
i,
j) = t_R(
i,
k) * t_u(
k,
j);
++t_h;
++t_R;
++t_diff_R;
++t_log_u;
++t_u;
++t_diff_u;
++t_eigen_vals;
++t_eigen_vecs;
++t_omega;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_div_P = getFTensor1FromMat<3>(
dataAtPts->divPAtPts);
auto t_s_dot_w = getFTensor1FromMat<3>(
dataAtPts->wDotAtPts);
dataAtPts->wDotDotAtPts.size2() != nb_integration_pts) {
dataAtPts->wDotDotAtPts.resize(3, nb_integration_pts);
}
auto t_s_dot_dot_w = getFTensor1FromMat<3>(
dataAtPts->wDotDotAtPts);
int nb_base_functions = data.
getN().size2();
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(
nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) +=
a * t_row_base_fun * t_div_P(
i);
t_nf(
i) -=
a * t_row_base_fun *
alphaW * t_s_dot_w(
i);
t_nf(
i) -=
a * t_row_base_fun *
alphaRho * t_s_dot_dot_w(
i);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_div_P;
++t_s_dot_w;
++t_s_dot_dot_w;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
int nb_base_functions = data.
getN().size2();
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(
nF);
t_PRT(
i,
m) = t_approx_P(
i,
j) * t_R(
m,
j);
t_leviPRT(
k) = levi_civita(
i,
m,
k) * t_PRT(
i,
m);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
k) +=
a * t_row_base_fun * t_leviPRT(
k);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_approx_P;
++t_R;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_dot_log_u =
getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachDotTensorAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
auto get_ftensor2 = [](
auto &
v) {
&
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
};
int nb_base_functions = data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor2(
nF);
t_RTP(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
t_deltaP(
i,
j) = t_RTP(
i,
j) - t_P(
i,
j);
t_viscous_stress(
i,
j) =
alphaU * t_dot_log_u(
i,
j);
a * (t_diff_u(
k,
l,
i,
j) * (t_deltaP(
k,
l) - t_viscous_stress(
k,
l)));
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
t_nf(
i,
j) += t_row_base_fun * t1(
i,
j);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_approx_P;
++t_P;
++t_R;
++t_dot_log_u;
++t_diff_u;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(
nF);
(levi_civita(
i,
m,
k) * t_omega_dot(
k)) * t_R(
m,
j);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum_omega(
i,
j);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_u;
++t_R;
++t_omega_dot;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
int nb_base_functions = data.
getN().size2() / 9;
auto get_ftensor0 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor0(
nF);
t_residuum_omega(
i,
j) = (levi_civita(
i,
m,
k) * t_omega_dot(
k)) * t_R(
m,
j);
int bb = 0;
for (; bb != nb_dofs; ++bb) {
t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum_omega(
i,
j);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
}
++t_w;
++t_R;
++t_u;
++t_omega_dot;
}
}
int nb_integration_pts = data.
getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_s_w = getFTensor1FromMat<3>(
dataAtPts->wAtPts);
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_nf = get_ftensor1(
nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
double div_row_base = t_row_diff_base_fun(
i,
i);
t_nf(
i) +=
a * div_row_base * t_s_w(
i);
++t_nf;
++t_row_diff_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_diff_base_fun;
}
++t_w;
++t_s_w;
}
}
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = data.
getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
t_bc_disp(
i) *= getFEMethod()->ts_t;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_bc_res(
i) = t_bc_disp(
i);
auto t_nf = get_ftensor1(
nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_res(
i) * 0.5;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
}
}
}
}
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = data.
getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.
getN().size2() / 3;
auto get_ftensor1 = [](
auto &
v) {
};
double theta = bc.theta;
theta *= getFEMethod()->ts_t;
const double a = sqrt(t_normal(
i) * t_normal(
i));
t_omega(
i) = t_normal(
i) * (theta /
a);
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_delta(
i) = t_center(
i) - t_coords(
i);
t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
auto t_nf = get_ftensor1(
nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_nf(
i) -= t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_coords;
}
}
}
}
int operator()(int p_row, int p_col, int p_data) const {
return 2 * (p_data + 1);
}
};
if (
ts_ctx == CTX_TSSETIFUNCTION) {
fe.getOpPtrVector().push_back(
new OpTractionBc(std::string("P") , *this));
fe.ts_t = ts_t;
fe);
}
}
auto t_normal = getFTensor1Normal();
const double nrm2 = sqrt(t_normal(
i) * t_normal(
i));
t_unit_normal(
i) = t_normal(
i) / nrm2;
int nb_integration_pts = data.
getN().size1();
int nb_base_functions = data.
getN().size2() / 3;
double ts_t = getFEMethod()->ts_t;
auto integrate_matrix = [&]() {
auto t_w = getFTensor0IntegrationWeight();
matPP.resize(nb_dofs / 3, nb_dofs / 3,
false);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != nb_dofs / 3; ++rr) {
const double a = t_row_base_fun(
i) * t_unit_normal(
i);
for (int cc = 0; cc != nb_dofs / 3; ++cc) {
const double b = t_col_base_fun(
i) * t_unit_normal(
i);
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
};
auto integrate_rhs = [&](auto &bc) {
auto t_w = getFTensor0IntegrationWeight();
vecPv.resize(3, nb_dofs / 3,
false);
double ts_t = getFEMethod()->ts_t;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != nb_dofs / 3; ++rr) {
const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
for (int dd = 0; dd != 3; ++dd)
if (bc.flags[dd])
vecPv(dd, rr) +=
t * bc.vals[dd];
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
};
auto integrate_rhs_cook = [&](auto &bc) {
vecPv.resize(3, nb_dofs / 3,
false);
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto calc_tau = [](double y) {
y -= 44;
y /= (60 - 44);
return -y * (y - 1) / 0.25;
};
const double tau = calc_tau(t_coords(1));
int rr = 0;
for (; rr != nb_dofs / 3; ++rr) {
const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
for (int dd = 0; dd != 3; ++dd)
if (bc.flags[dd])
vecPv(dd, rr) +=
t * tau * bc.vals[dd];
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_coords;
}
};
if (bc.faces.find(fe_ent) != bc.faces.end()) {
if (nb_dofs) {
if (bc.blockName.compare(20, 4, "COOK") == 0)
CHKERR integrate_rhs_cook(bc);
else
const auto info =
nb_dofs / 3, &*
vecPv.data().begin(), nb_dofs / 3);
if (info > 0)
"The leading minor of order %i is "
"not positive; definite;\nthe "
"solution could not be computed",
info);
for (int dd = 0; dd != 3; ++dd)
if (bc.flags[dd])
for (int rr = 0; rr != nb_dofs / 3; ++rr)
}
}
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor1(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
double div_col_base = t_col_diff_base_fun(
i,
i);
t_m(
i) +=
a * t_row_base_fun * div_col_base;
++t_m;
++t_col_diff_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
if (
alphaW < std::numeric_limits<double>::epsilon() &&
alphaRho < std::numeric_limits<double>::epsilon())
const int nb_integration_pts = row_data.
getN().size1();
const int row_nb_dofs = row_data.
getIndices().size();
&
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.
getN().size2();
double ts_scale =
alphaW * getTSa();
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon())
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a =
v * t_w * ts_scale;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor1(
K, 3 * rr, 0);
for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
const double b =
a * t_row_base_fun * t_col_base_fun;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
&
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
&
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
&
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor3(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_RTP_dP(
i,
j,
k) = t_R(
k,
i) * t_col_base_fun(
j) * t_row_base_fun;
for (int kk = 0; kk != 3; ++kk)
for (int ii = 0; ii != 3; ++ii)
for (int jj = ii; jj != 3; ++jj) {
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn)
t_m(ii, jj, kk) +=
a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
}
++t_col_base_fun;
++t_m;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_R;
++t_diff_u;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor2(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs; ++cc) {
a * t_row_base_fun * t_R(
k,
i) * t_col_base_fun(
k,
j);
t_m(
i,
j) += t_diff_u(
m,
n,
i,
j) * t_RTP_dBubble(
m,
n);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_R;
++t_diff_u;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
&
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
&
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
&
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
&
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
&
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
&
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
&
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
&
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
&
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
auto t_dot_log_u =
getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachDotTensorAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
t_one4_symmetric(
i,
j,
k,
l) = 0;
for (auto ii : {0, 1, 2})
for (auto jj : {0, 1, 2})
for (auto kk : {0, 1, 2})
for (auto ll : {0, 1, 2}) {
if (ll != kk)
t_one4_symmetric(ii, jj, kk, ll) =
t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
else
t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
}
int row_nb_base_functions = row_data.
getN().size2();
for (auto ii : {0, 1, 2})
t_one(ii, ii) = 1;
const double ts_a = getTSa();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_RTP(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
t_deltaP(
i,
j) = t_RTP(
i,
j) - t_P(
i,
j);
t_dot_u(
i,
j) = t_u(
i,
k) * t_dot_log_u(
k,
j);
t_stress_diff(
i,
j,
k,
l) = 0;
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int kk = 0; kk != 3; ++kk)
for (int ll = 0; ll != 3; ++ll)
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn)
t_stress_diff(ii, jj, mm, nn) -=
t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
t_viscous_stress(
i,
j) =
alphaU * t_dot_log_u(
i,
j);
t_viscous_stress_diff(
i,
j,
k,
l) = t_one4_symmetric(
i,
j,
k,
l);
t_viscous_stress_diff(
i,
j,
k,
l) *= (ts_a *
alphaU);
t_deltaP2(
i,
j) = t_deltaP(
i,
j) - t_viscous_stress(
i,
j);
t_eigen_vals, t_eigen_vecs,
f,
d_f,
dd_f, t_deltaP2, nbUniq[gg]);
for (int ii = 0; ii != 3; ++ii) {
for (int jj = 0; jj < ii; ++jj) {
for (int kk = 0; kk != 3; ++kk) {
for (int ll = 0; ll != kk; ++ll) {
t_diff2_uP2(ii, jj, kk, ll) *= 2;
}
}
}
}
for (int ii = 0; ii != 3; ++ii) {
for (int jj = ii; jj != 3; ++jj) {
for (int kk = 0; kk != 3; ++kk) {
for (int ll = 0; ll < kk; ++ll) {
t_diff2_uP2(ii, jj, kk, ll) *= 2;
}
}
}
}
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor4(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
const double b =
a * t_row_base_fun * t_col_base_fun;
for (int ii = 0; ii != 3; ++ii)
for (int jj = ii; jj != 3; ++jj)
for (int kk = 0; kk != 3; ++kk)
for (int ll = kk; ll != 3; ++ll)
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn)
t_m(ii, jj, kk, ll) +=
b * t_diff_u(mm, nn, ii, jj) *
(t_stress_diff(mm, nn, kk, ll) -
t_viscous_stress_diff(mm, nn, kk, ll));
t_m(
i,
j,
k,
l) += b * t_diff2_uP2(
i,
j,
k,
l);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_P_dh0;
++t_P_dh1;
++t_P_dh2;
++t_R;
++t_approx_P;
++t_P;
++t_dot_log_u;
++t_u;
++t_diff_u;
++t_eigen_vals;
++t_eigen_vecs;
}
}
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
&
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
&
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
&
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_diff_u =
getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
int row_nb_base_functions = row_data.
getN().size2();
int nb_integration_pts = row_data.
getN().size1();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_RT_domega_P(
i,
j,
m) = t_diff_R(
k,
i,
m) * t_P(
k,
j);
for (int ii = 0; ii != 3; ++ii)
for (int jj = ii; jj != 3; ++jj)
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn)
for (int kk = 0; kk != 3; ++kk)
t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor3(
K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
double v =
a * t_row_base_fun * t_col_base_fun;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_P;
++t_diff_R;
++t_diff_u;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(
K, 3 * rr, 0);
(
a * t_row_base_fun) * (levi_civita(
i,
m,
k) * t_R(
m,
j));
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
k,
i) += t_levi_R(
i,
j,
k) * t_col_base_fun(
j);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_R;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor1(
K, 3 * rr, 0);
(
a * t_row_base_fun) * (levi_civita(
i,
m,
k) * t_R(
m,
j));
for (int cc = 0; cc != col_nb_dofs; ++cc) {
t_m(
k) += t_levi_r(
i,
j,
k) * t_col_base_fun(
i,
j);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_R;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_PRT(
i,
m) = t_approx_P(
i,
j) * t_R(
m,
j);
levi_civita(
i,
j,
n) * (t_approx_P(
i,
k) * t_diff_R(
j,
k,
l));
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
const double b =
a * t_row_base_fun * t_col_base_fun;
t_m(
k,
m) += b * t_leviPRT_domega(
k,
m);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_approx_P;
++t_R;
++t_diff_R;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
&
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
&
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
int row_nb_base_functions = row_data.
getN().size2() / 3;
const double ts_a = getTSa();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
t_PRT(
i,
k) = t_row_base_fun(
j) * (t_diff_R(
i,
l,
k) * t_u(
l,
j));
t_levi1(
i,
k) = (levi_civita(
i,
m,
n) * t_omega_dot(
n)) *
(t_row_base_fun(
j) * t_diff_R(
m,
j,
k));
t_levi2(
i,
k) = levi_civita(
i,
m,
k) * (t_row_base_fun(
j) * t_R(
m,
j));
auto t_m = get_ftensor2(
K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
i,
j) += (
a * t_col_base_fun) * t_PRT(
i,
j);
t_m(
i,
j) += (
a * t_col_base_fun) * t_levi1(
i,
j);
t_m(
i,
j) += ((
a * ts_a) * t_col_base_fun) * t_levi2(
i,
j);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_R;
++t_diff_R;
++t_u;
++t_omega_dot;
}
}
int nb_integration_pts = row_data.
getN().size1();
&
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
};
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
int row_nb_base_functions = row_data.
getN().size2() / 9;
const double ts_a = getTSa();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
t_PRT(
k) = t_row_base_fun(
i,
j) * (t_diff_R(
i,
l,
k) * t_u(
l,
j));
t_levi_omega(
i,
m) = levi_civita(
i,
m,
n) * t_omega_dot(
n);
t_levi_omegaDiffR(
i,
j,
k) = t_levi_omega(
i,
m) * t_diff_R(
m,
j,
k);
t_levi_omegaR(
i,
j,
k) = levi_civita(
i,
m,
k) * t_R(
m,
j);
t_levi1(
k) = t_row_base_fun(
i,
j) * t_levi_omegaDiffR(
i,
j,
k);
t_levi2(
k) = t_row_base_fun(
i,
j) * t_levi_omegaR(
i,
j,
k);
auto t_m = get_ftensor2(
K, rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(
j) += (
a * t_col_base_fun) * t_PRT(
j);
t_m(
j) += (
a * t_col_base_fun) * t_levi1(
j);
t_m(
j) += ((
a * ts_a) * t_col_base_fun) * t_levi2(
j);
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_R;
++t_diff_R;
++t_u;
++t_omega_dot;
}
}
if (row_type != MBTET)
}
const int nb_integration_pts = row_data.
getN().size1();
const int row_nb_dofs = row_data.
getIndices().size();
&
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
);
};
auto t_w = getFTensor0IntegrationWeight();
m.resize(row_nb_dofs, row_nb_dofs,
false);
int row_nb_base_functions = row_data.
getN().size2();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor1(
m, 3 * rr, 0);
for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
const double b =
a * t_row_base_fun * t_col_base_fun;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
if (row_type != MBTET)
if (auto ep_fe_vol_ptr =
dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
getFEMethod())) {
auto set_data = [&](auto fe) {
aoSuu = fe->aoSuu;
aoSBB = fe->aoSBubble;
aoSOO = fe->aoSOmega;
aoSww = fe->aoSw;
Suu = fe->Suu;
SBB = fe->SBubble;
SOO = fe->SOmega;
Sww = fe->Sw;
};
set_data(ep_fe_vol_ptr);
if (Suu) {
auto &row_name, auto &col_name, auto row_side,
auto col_side, auto row_type, auto col_type) {
return add_bmc.get<0>().find(boost::make_tuple(
row_name, col_name, row_type, col_type, row_side, col_side));
};
auto &row_name, auto &col_name, auto row_side,
auto col_side, auto row_type, auto col_type,
const auto &
m,
const auto &row_ind,
const auto &col_ind) {
auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
row_type, col_type);
if (it != add_bmc.get<0>().end()) {
it->setInd(row_ind, col_ind);
it->setSetAtElement();
return it;
} else {
auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
row_name, col_name, row_type, col_type, row_side, col_side,
m,
row_ind, col_ind));
return p_it.first;
}
};
auto &row_name, auto &col_name, auto row_side,
auto col_side, auto row_type, auto col_type,
const auto &
m,
const auto &row_ind,
const auto &col_ind) {
auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
row_type, col_type);
if (it != add_bmc.get<0>().end()) {
it->setInd(row_ind, col_ind);
it->setSetAtElement();
return it;
} else {
auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
row_name, col_name, row_type, col_type, row_side, col_side,
m,
row_ind, col_ind));
return p_it.first;
}
};
auto assemble_block = [&](
auto &
bit, Mat S) {
CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
&*cind.begin(), &*
m.data().begin(), ADD_VALUES);
};
const int nb =
m.size1();
inv.resize(nb, nb, false);
inv.clear();
for (int cc = 0; cc != nb; ++cc)
inv(cc, cc) = -1;
const auto info =
lapack_dsysv(
'L', nb, nb, &*
m.data().begin(), nb,
&*
iPIV.begin(), &*inv.data().begin(), nb,
};
std::string field, auto &inv) {
bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
if (
bit != bmc.get<1>().end()) {
} else
"%s matrix not found", field.c_str());
};
const int nb =
m.size1();
trans_inv.resize(nb, nb, false);
trans_inv.clear();
for (
int c = 0;
c != nb; ++
c)
const auto info =
&*trans_inv.data().begin(), nb);
if (info != 0)
"Can not invert matrix info = %d", info);
inv.resize(nb, nb, false);
noalias(inv) = trans(trans_inv);
};
auto invert_nonsymm_schur =
auto &inv,
const bool debug =
false) {
auto bit = bmc.get<1>().find(
boost::make_tuple(field, field, MBTET, MBTET));
if (
bit != bmc.get<1>().end()) {
std::cerr << prod(
m, inv) << endl;
std::cerr << endl;
}
} else
"%s matrix not found", field.c_str());
};
auto create_block_schur =
for (
auto &
bit : add_bmc) {
}
if (
bit.setAtElement &&
bit.rowField != field &&
if (ao) {
CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
}
auto it = set_block(add_bmc,
bit.rowField,
bit.colField,
bit.colType,
m, rind, cind);
}
}
for (auto &bit_col : bmc) {
if (bit_col.setAtElement && bit_col.rowField == field &&
bit_col.colField != field) {
invMat.resize(inv.size1(), cm.size2(),
false);
noalias(
invMat) = prod(inv, cm);
if (ao)
CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
for (auto &bit_row : bmc) {
if (bit_row.setAtElement && bit_row.rowField != field &&
bit_row.colField == field) {
K.resize(rind.size(), cind.size(),
false);
if (ao)
CHKERR AOApplicationToPetsc(ao, rind.size(),
&*rind.begin());
auto it = add_block(add_bmc, bit_row.rowField,
bit_col.colField, bit_row.rowSide,
bit_col.colSide, bit_row.rowType,
bit_col.colType,
K, rind, cind);
}
}
}
}
};
auto assemble_schur =
for (
auto &
bit : add_bmc) {
}
for (
auto &
bit : add_bmc) {
std::cerr <<
"assemble: " <<
bit.rowField <<
" "
std::cerr <<
bit.M << endl;
}
}
std::cerr << std::endl;
}
};
auto precondition_schur =
for (
auto &
bit : add_bmc) {
}
if (
bit.rowField != field ||
bit.colField != field)
auto it =
set_block(add_bmc,
bit.rowField,
bit.colField,
}
}
auto bit = bmc.get<1>().find(
boost::make_tuple(field, field, MBTET, MBTET));
if (
bit->setAtElement &&
bit != bmc.get<1>().end()) {
auto it =
set_block(add_bmc,
bit->rowField,
bit->colField,
bit->rowSide,
} else {
auto row_it = bmc.get<3>().lower_bound(field);
for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
if (row_it->setAtElement) {
auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
diag_mat, row_it->rowInd, row_it->rowInd);
break;
}
}
if (row_it == bmc.get<3>().end())
"row field not found %s", field.c_str());
}
};
if (SBB) {
if (SOO) {
} else {
}
if (Sww) {
}
}
}
}
}
}
if (type != MBVERTEX)
auto create_tag = [this](const std::string tag_name, const int size) {
double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
Tag th;
th, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
return th;
};
Tag th_w = create_tag("SpatialDisplacement", 3);
Tag th_omega = create_tag("Omega", 3);
Tag th_approxP = create_tag("Piola1Stress", 9);
Tag th_sigma = create_tag("CauchyStress", 9);
Tag th_log_u = create_tag("LogSpatialStreach", 9);
Tag th_u = create_tag("SpatialStreach", 9);
Tag th_h = create_tag("h", 9);
Tag th_X = create_tag("X", 3);
Tag th_detF = create_tag("detF", 1);
Tag th_angular_momentum = create_tag("AngularMomentum", 3);
Tag th_u_eig_vec = create_tag("SpatialStreachEigenVec", 9);
Tag th_u_eig_vals = create_tag("SpatialStreachEigenVals", 3);
int nb_gauss_pts = data.
getN().size1();
if (
mapGaussPts.size() != (
unsigned int)nb_gauss_pts) {
"Nb. of integration points is not equal to number points on "
"post-processing mesh");
}
auto t_w = getFTensor1FromMat<3>(
dataAtPts->wAtPts);
auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
auto t_log_u =
getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachTensorAtPts);
auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
};
&
m(0, 0), &
m(0, 1), &
m(0, 2),
&
m(1, 0), &
m(1, 1), &
m(1, 2),
&
m(2, 0), &
m(2, 1), &
m(2, 2));
auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
};
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
CHKERR save_vec_tag(th_w, t_w, gg);
CHKERR save_vec_tag(th_X, t_coords, gg);
CHKERR save_vec_tag(th_omega, t_omega, gg);
CHKERR save_mat_tag(th_h, t_h, gg);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_u_tmp(ii, jj) = t_u(ii, jj);
CHKERR save_mat_tag(th_u, t_u_tmp, gg);
CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
const double jac = determinantTensor3by3(t_h);
t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
t_PhT(
i,
k) = t_approx_P(
i,
j) * t_R(
k,
j);
t_leviPRT(
k) = levi_civita(
i,
l,
k) * t_PhT(
i,
l);
&t_leviPRT(0));
auto get_eiegn_vector_symmetric = [&](auto &t_u) {
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_m(ii, jj) = t_u(ii, jj);
int n = 3, lda = 3, info, lwork = -1;
double wkopt;
&(eigen_values.data()[0]), &wkopt, lwork);
if (info != 0)
"is something wrong with lapack_dsyev info = %d", info);
lwork = (int)wkopt;
double work[lwork];
&(eigen_values.data()[0]), work, lwork);
if (info != 0)
"is something wrong with lapack_dsyev info = %d", info);
&*eigen_values.data().begin());
};
CHKERR get_eiegn_vector_symmetric(t_u);
++t_w;
++t_h;
++t_log_u;
++t_u;
++t_omega;
++t_R;
++t_approx_P;
++t_coords;
}
}
DataForcesAndSourcesCore::EntData &data) {
if (type == MBTET) {
int nb_integration_pts = data.getN().size1();
auto t_w = getFTensor0IntegrationWeight();
auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = t_w *
v;
(*energy) +=
a * t_P(
i,
J) * t_h(
i,
J);
++t_w;
++t_P;
++t_h;
}
}
}
DataForcesAndSourcesCore::EntData &data) {
if (type == MBTET) {
}
}
}
Eshelbian plasticity interface.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
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)
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
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)
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
MoFEMErrorCode get_eigen_val_and_proj_lapack(T1 &X, T2 &eig, T3 &eig_vec)
bool is_eq(const double &a, const double &b)
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
void sort_eigen_vals(T1 &eig, T2 &eigen_vec)
double dd_f(const double v)
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
double d_f(const double v)
Tensors class implemented by Walter Landry.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< int > VectorInt
implementation of Data Operators for Forces and Sources
constexpr double t
plate stiffness
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainor
MoFEMErrorCode preProcess()
boost::shared_ptr< TractionBcVec > bcData
MoFEM::Interface & mField
MatrixDouble K
local tangent matrix
VectorDouble nF
local right hand side vector
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< BcDispVec > bcDispPtr
MoFEMErrorCode integrate(EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(EntData &data)
boost::shared_ptr< BcRotVec > bcRotPtr
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data)
MoFEMErrorCode assemble(int side, EntityType type, EntData &data)
map< std::string, MatrixDouble > invBlockMat
map< std::string, DataAtIntegrationPts::BlockMatContainor > blockMat
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Set integration rule to boundary elements.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
intrusive_ptr for managing petsc objects