v0.8.23
EshelbianOperators.cpp

Implementation of operators

/**
* \file EshelbianOperators.cpp
* \example EshelbianOperators.cpp
*
* \brief Implementation of operators
*/
#include <MoFEM.hpp>
using namespace MoFEM;
#include <boost/math/constants/constants.hpp>
#include <lapack_wrap.h>
MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
EntityType type,
EntData &data) {
if (type != MBTET)
int nb_integration_pts = data.getN().size1();
dataAtPts->hAtPts->resize(9, nb_integration_pts, false);
dataAtPts->hDeltaAtPts->resize(9, nb_integration_pts, false);
dataAtPts->xDetAtPts->resize(nb_integration_pts, false);
dataAtPts->xInvGradAtPts->resize(9, nb_integration_pts, false);
dataAtPts->rotMatAtPts->resize(9, nb_integration_pts, false);
dataAtPts->diffRotMatAtPts->resize(27, nb_integration_pts, false);
auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
auto t_x_inv_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xInvGradAtPts));
auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
auto t_omega = getFTensor1FromMat<3>(*(dataAtPts->rotAxisAtPts));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
CHKERR determinantTensor3by3(t_x_grad, t_x_det);
CHKERR invertTensor3by3(t_x_grad, t_x_det, t_x_inv_grad);
auto t0_diff = getDiffRotationFormVector(t_omega);
auto t0 = getRotationFormVector(t_omega);
t_diff_R(i, j, k) = t0_diff(i, j, k);
t_R(i, j) = t0(i, j);
t_delta_h(i, j) = t_R(i, k) * t_u(k, j) + t_R(i, j);
t_h(i, j) = t_delta_h(i, k) * t_x_grad(k, j);
++t_h;
++t_delta_h;
++t_x_grad;
++t_x_det;
++t_x_inv_grad;
++t_R;
++t_diff_R;
++t_u;
++t_omega;
}
}
MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_div_P = getFTensor1FromMat<3>(*(dataAtPts->divPAtPts));
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_div_P;
}
}
MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor1(nF);
// FTensor::Tensor2<double, 3, 3> t_PhT;
// t_PhT(i, m) = (t_x_grad(j, l) * t_approx_P(i, l)) * t_delta_h(m, j);
t_PRT(i, m) = (t_x_grad(j, l) * t_approx_P(i, l)) * 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_x_grad;
++t_delta_h;
++t_R;
}
}
MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
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_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
auto t_dot_u =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachDotTensorAtPts));
auto get_ftensor2 = [](auto &v) {
&v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
};
int nb_base_functions = data.getN().size2();
auto t_row_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor2(nF);
t_xP(i, j) = t_x_grad(j, k) * t_P(i, k);
t_RTxP(i, j) = t_R(k, i) * t_x_grad(j, l) * t_approx_P(k, l);
int bb = 0;
for (; bb != nb_dofs / 6; ++bb) {
for (int ii : {0, 1, 2})
for (int jj : {0, 1, 2}) {
t_nf(ii, jj) += a * t_row_base_fun * (t_RTxP(ii, jj) - t_xP(ii, jj));
t_nf(ii, jj) -= a * t_row_base_fun * alpha_u * t_dot_u(ii, jj);
}
++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_x_grad;
++t_dot_u;
}
}
MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int nb_base_functions = data.getN().size2() / 3;
auto t_row_base_fun = data.getFTensor1N<3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor1(nF);
int bb = 0;
for (; bb != nb_dofs / 3; ++bb) {
t_pushed_base(i) = t_x_grad(i, k) * t_row_base_fun(k);
t_nf(i) += a * t_pushed_base(j) * t_delta_h(i, j);
t_nf(i) -= a * t_pushed_base(i);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
}
++t_w;
++t_delta_h;
++t_x_grad;
}
}
MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int nb_base_functions = data.getN().size2() / 9;
auto t_row_base_fun = data.getFTensor2N<3, 3>();
auto get_ftensor0 = [](auto &v) {
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
auto t_nf = get_ftensor0(nF);
int bb = 0;
for (; bb != nb_dofs; ++bb) {
t_pushed_base(i, j) = t_x_grad(j, k) * t_row_base_fun(i, k);
t_nf += a * t_pushed_base(i, j) * t_delta_h(i, j);
t_nf -= a * t_pushed_base(i, i);
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb) {
++t_row_base_fun;
}
++t_w;
++t_delta_h;
++t_x_grad;
}
}
MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_s_w = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
int nb_base_functions = data.getN().size2() / 3;
auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
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;
}
}
MoFEMErrorCode OpDispBc::integrate(EntData &data) {
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
// interate over all boundary data
for (auto &bc : (*bcDispPtr)) {
// check if finite element entity is part of boundary condition
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.getN().size2() / 3;
auto t_row_base_fun = data.getFTensor1N<3>();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*(dataAtPts->xAtPts));
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
// get bc data
FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
t_bc_disp(i) *= getFEMethod()->ts_t;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_bc_res(i) = t_bc_disp(i) - t_x(i) + t_coords(i);
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_bc_res(i) * 0.5;
++t_nf;
++t_row_base_fun;
}
for (; bb != nb_base_functions; ++bb)
++t_row_base_fun;
++t_w;
++t_coords;
++t_x;
}
}
}
}
MoFEMErrorCode OpDispBc_dx::integrate(EntData &row_data, EntData &col_data) {
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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));
};
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
// interate over all boundary data
for (auto &bc : (*bcDispPtr)) {
// check if finite element entity is part of boundary condition
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = row_data.getN().size2() / 3;
auto t_row_base_fun = row_data.getFTensor1N<3>();
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);
const double val = t_w * (t_row_base_fun(i) * t_normal(i)) * 0.5;
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
const double set_val = val * t_col_base_fun;
t_m(0, 0) += set_val;
t_m(1, 1) += set_val;
t_m(2, 2) += set_val;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
}
}
MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
// interate over all boundary data
for (auto &bc : (*bcRotPtr)) {
// check if finite element entity is part of boundary condition
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_dofs = data.getIndices().size();
int nb_integration_pts = data.getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
int nb_base_functions = data.getN().size2() / 3;
auto t_row_base_fun = data.getFTensor1N<3>();
auto get_ftensor1 = [](auto &v) {
&v[2]);
};
// get bc data
FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
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_R = getRotationFormVector(t_omega);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*(dataAtPts->xAtPts));
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);
t_disp(i) -= t_x(i) - t_coords(i);
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;
++t_x;
}
}
}
}
MoFEMErrorCode OpRotationBc_dx::integrate(EntData &row_data,
EntData &col_data) {
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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));
};
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
// interate over all boundary data
for (auto &bc : (*bcRotPtr)) {
// check if finite element entity is part of boundary condition
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
int row_nb_base_functions = row_data.getN().size2() / 3;
auto t_row_base_fun = row_data.getFTensor1N<3>();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
// get bc data
for (int gg = 0; gg != nb_integration_pts; ++gg) {
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
const double val = t_w * (t_row_base_fun(i) * t_normal(i)) * 0.5;
auto t_m = get_ftensor2(K, 3 * rr, 0);
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
const double set_val = val * t_col_base_fun;
t_m(0, 0) += set_val;
t_m(1, 1) += set_val;
t_m(2, 2) += set_val;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
}
}
MoFEMErrorCode FeTractionBc::preProcess() {
struct FaceRule {
int operator()(int p_row, int p_col, int p_data) const {
return 2 * (p_data + 1);
}
};
if (ts_ctx == CTX_TSSETIFUNCTION) {
// Loop boundary elements with traction boundary conditions
fe.getOpPtrVector().push_back(
new OpTractionBc(std::string("P")/* + "_RT"*/, *this));
fe.getRuleHook = FaceRule();
fe.ts_t = ts_t;
CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
fe);
}
}
MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
auto t_normal = getFTensor1Normal();
const double nrm2 = sqrt(t_normal(i) * t_normal(i));
t_unit_normal(i) = t_normal(i) / nrm2;
int nb_dofs = data.getFieldData().size();
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);
matPP.clear();
auto t_row_base_fun = data.getFTensor1N<3>();
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);
auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
for (int cc = 0; cc != nb_dofs / 3; ++cc) {
const double b = t_col_base_fun(i) * t_unit_normal(i);
matPP(rr, cc) += t_w * a * b;
++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);
vecPv.clear();
auto t_row_base_fun = data.getFTensor1N<3>();
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);
vecPv.clear();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base_fun = data.getFTensor1N<3>();
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;
}
};
// get entity of face
EntityHandle fe_ent = getFEEntityHandle();
for (auto &bc : *(bcFe.bcData)) {
if (bc.faces.find(fe_ent) != bc.faces.end()) {
int nb_dofs = data.getFieldData().size();
if (nb_dofs) {
CHKERR integrate_matrix();
if (bc.blockName.compare(20, 4, "COOK") == 0)
CHKERR integrate_rhs_cook(bc);
else
CHKERR integrate_rhs(bc);
const auto info =
lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
if (info > 0)
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"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)
data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
// int nb_integration_pts = data.getN().size1();
// auto t_row_base_fun = data.getFTensor1N<3>();
// auto t_coords = getFTensor1CoordsAtGaussPts();
// for (int gg = 0; gg != nb_integration_pts; ++gg) {
// FTensor::Tensor1<double, 3> t_v;
// t_v(i) = 0;
// int rr = 0;
// for (; rr != nb_dofs / 3; ++rr) {
// for (int dd = 0; dd != 3; ++dd) {
// const double v = data.getFieldDofs()[3 * rr +
// dd]->getFieldData(); t_v(dd) += t_row_base_fun(i) *
// t_unit_normal(i) * v;
// }
// ++t_row_base_fun;
// }
// for (; rr != nb_base_functions; ++rr)
// ++t_row_base_fun;
// ++t_coords;
// cerr << t_v << endl;
// }
}
}
}
}
MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
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;
}
}
MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
auto t_m = get_ftensor3(K, 6 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_pushed_base(i) = t_x_grad(i, k) * t_col_base_fun(k);
t_RTP_dP(i, j, k) = t_R(k, i) * t_pushed_base(j) * t_row_base_fun;
for (int kk = 0; kk != 3; ++kk)
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj) {
t_m(ii, jj, kk) += a * t_RTP_dP(ii, jj, 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_x_grad;
}
}
MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
&m(r + 5, c));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_m = get_ftensor2(K, 6 * rr, 0);
auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
for (int cc = 0; cc != col_nb_dofs; ++cc) {
t_pushed_base(i, j) = t_x_grad(j, k) * t_col_base_fun(i, k);
t_RTP_dBubble(i, j) = t_R(k, i) * t_pushed_base(k, j) * t_row_base_fun;
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj) {
t_m(ii, jj) += a * t_RTP_dBubble(ii, jj);
}
++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_x_grad;
}
}
MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_P_dh0 = getFTensor3FromMat(*(dataAtPts->P_dh0));
auto t_P_dh1 = getFTensor3FromMat(*(dataAtPts->P_dh1));
auto t_P_dh2 = getFTensor3FromMat(*(dataAtPts->P_dh2));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
const double ts_a = getFEMethod()->ts_a;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
t_P_dh(i, j, N0, k) = t_x_grad(j, m) * t_P_dh0(i, m, k);
t_P_dh(i, j, N1, k) = t_x_grad(j, m) * t_P_dh1(i, m, k);
t_P_dh(i, j, N2, k) = t_x_grad(j, m) * t_P_dh2(i, m, k);
t_pushed_P_du(i, j, k, l) = t_P_dh(i, j, k, n) * t_x_grad(l, n);
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
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, 1, 2})
for (int jj : {0, 1, 2})
for (int kk : {0, 1, 2})
for (int ll : {0, 1, 2})
t_m(ii, jj, kk, ll) -= b * t_pushed_P_du(ii, jj, kk, ll);
for (int ii : {0, 1, 2})
for (int jj : {0, 1, 2})
t_m(ii, jj, ii, jj) -= b * alpha_u * ts_a;
++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_x_grad;
}
}
MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
EntData &col_data) {
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
int nb_integration_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
t_P(i, j) = t_x_grad(j, k) * t_approx_P(i, k);
t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
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;
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int kk = 0; kk != 3; ++kk)
t_m(ii, jj, kk) += v * t_RT_domega_P(ii, jj, kk);
++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_diff_R;
++t_x_grad;
}
}
MoFEMErrorCode OpSpatialPhysical_du_dx::integrate(EntData &row_data,
EntData &col_data) {
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->PAtPts));
auto t_P_dh0 = getFTensor3FromMat(*(dataAtPts->P_dh0));
auto t_P_dh1 = getFTensor3FromMat(*(dataAtPts->P_dh1));
auto t_P_dh2 = getFTensor3FromMat(*(dataAtPts->P_dh2));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
int nb_integration_pts = row_data.getN().size1();
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
t_delta_P(i, j) = t_R(l, i) * t_approx_P(l, j) - t_P(i, j);
t_P_dh(i, j, N0, l) = t_x_grad(j, m) * t_P_dh0(i, m, l);
t_P_dh(i, j, N1, l) = t_x_grad(j, m) * t_P_dh1(i, m, l);
t_P_dh(i, j, N2, l) = t_x_grad(j, m) * t_P_dh2(i, m, l);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
t_delta_h(ii, jj) = t_u(ii, jj);
for (int ii = 0; ii != 3; ++ii)
t_delta_h(ii, ii) += 1;
t_P_intermidiate_dh(i, j, k, l) = t_P_dh(i, j, m, l) * t_delta_h(m, k);
int rr = 0;
for (; rr != row_nb_dofs / 6; ++rr) {
auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
auto t_m = get_ftensor3(K, 6 * rr, 0);
const double b = a * t_row_base_fun;
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
FTensor::Tensor1<double, 3> t_dgrad_x_deltaP;
t_dgrad_x_deltaP(i) = b * t_col_diff_base_fun(k) * t_delta_P(i, k);
t_col_P_delta_dh(i, j, k) =
b * (t_P_intermidiate_dh(i, j, k, l) * t_col_diff_base_fun(l));
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj) {
t_m(ii, jj, jj) += t_dgrad_x_deltaP(ii);
for (int kk = 0; kk != 3; ++kk)
t_m(ii, jj, kk) -= t_col_P_delta_dh(ii, jj, kk);
}
++t_m;
++t_col_diff_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_approx_P;
++t_P;
++t_P_dh0;
++t_P_dh1;
++t_P_dh2;
++t_u;
++t_x_grad;
++t_R;
++t_w;
}
}
MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
auto t_m = get_ftensor2(K, 3 * rr, 0);
// t_levi_h(i, j, k) =
// (a * t_row_base_fun) * (levi_civita(i, m, k) * t_delta_h(m, j));
t_levi_h(i, j, k) =
(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_pushed_base(i) = t_x_grad(i, k) * t_col_base_fun(k);
t_m(k, i) += t_levi_h(i, j, k) * t_pushed_base(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_delta_h;
++t_R;
++t_x_grad;
}
}
MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
auto t_m = get_ftensor1(K, 3 * rr, 0);
// t_levi_h(i, j, k) =
// (a * t_row_base_fun) * (levi_civita(i, m, k) * t_delta_h(m, j));
t_levi_h(i, j, k) =
(a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
for (int cc = 0; cc != col_nb_dofs; ++cc) {
t_pushed_base(i, j) = t_x_grad(j, k) * t_col_base_fun(i, k);
t_m(k) += t_levi_h(i, j, k) * t_pushed_base(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_delta_h;
++t_R;
++t_x_grad;
}
}
MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
t_pushed_P(i, j) = t_x_grad(j, k) * t_approx_P(i, k);
// FTensor::Tensor2<double, 3, 3> t_leviPhT_domega;
// t_leviPhT_domega(n, l) =
// levi_civita(i, j, n) *
// (t_pushed_P(i, k) *
// (t_diff_R(j, m, l) * t_u(m, k) + t_diff_R(j, k, l))
// );
t_leviPRT_domega(n, l) =
levi_civita(i, j, n) * (t_pushed_P(i, k) * t_diff_R(j, k, l));
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
auto t_m = get_ftensor2(K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(k, m) +=
(a * t_row_base_fun * t_col_base_fun) * 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_u;
++t_diff_R;
++t_x_grad;
}
}
MoFEMErrorCode OpSpatialRotation_domega_du::integrate(EntData &row_data,
EntData &col_data) {
// int nb_integration_pts = row_data.getN().size1();
// int row_nb_dofs = row_data.getIndices().size();
// int col_nb_dofs = col_data.getIndices().size();
// auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
// return FTensor::Christof<FTensor::PackPtr<double *, 6>, 3, 3>(
// &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)
// );
// };
// FTensor::Index<'i', 3> i;
// FTensor::Index<'j', 3> j;
// FTensor::Index<'k', 3> k;
// FTensor::Index<'l', 3> l;
// FTensor::Index<'m', 3> m;
// FTensor::Index<'n', 3> n;
// auto v = getVolume();
// auto t_w = getFTensor0IntegrationWeight();
// auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
// auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
// auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
// auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
// auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
// int row_nb_base_functions = row_data.getN().size2();
// auto t_row_base_fun = row_data.getFTensor0N();
// for (int gg = 0; gg != nb_integration_pts; ++gg) {
// double a = v * t_w;
// FTensor::Tensor2<double, 3, 3> t_pushed_P;
// t_pushed_P(i, j) = t_x_grad(j, k) * t_approx_P(i, k);
// FTensor::Tensor4<double, 3, 3, 3, 3> t_pushed_PhT_du;
// t_pushed_PhT_du(i, j, k, l) = t_pushed_P(i, l) * t_R(j, k);
// FTensor::Tensor3<double, 3, 3, 3> t_levi;
// t_levi(i, j, k) = levi_civita(i, j, k);
// FTensor::Tensor3<double, 3, 3, 3> t_leviPhT_du;
// t_leviPhT_du(n, k, l) = 0;
// for (int ii : {0, 1, 2})
// for (int jj : {0, 1, 2})
// for (int nn : {0, 1, 2})
// for (int kk : {0, 1, 2})
// for (int ll : {0, 1, 2})
// t_leviPhT_du(nn, kk, ll) += t_levi(ii, jj, nn) *
// t_pushed_PhT_du(ii, jj, kk, ll);
// int rr = 0;
// for (; rr != row_nb_dofs / 3; ++rr) {
// auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
// auto t_m = get_ftensor3(K, 3 * 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, 1, 2})
// for (int jj : {0, 1, 2})
// for (int kk : {0, 1, 2})
// t_m(ii, jj, kk) += b * t_leviPhT_du(ii, jj, kk);
// ++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_u;
// ++t_x_grad;
// }
}
MoFEMErrorCode OpSpatialRotation_domega_dx::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
// FTensor::Tensor3<double, 3, 3, 3> t_levihT;
// t_levihT(i, l, k) = levi_civita(i, j, k) * t_delta_h(j, l);
t_leviRT(i, l, k) = levi_civita(i, j, k) * t_R(j, l);
auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
auto t_m = get_ftensor2(K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(k, l) +=
a * t_row_base_fun *
((t_col_diff_base_fun(m) * t_approx_P(i, m)) * t_leviRT(i, l, k));
++t_col_diff_base_fun;
++t_m;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_approx_P;
++t_delta_h;
++t_R;
}
}
MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2() / 3;
auto t_row_base_fun = row_data.getFTensor1N<3>();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
t_pushed_base(i) = t_x_grad(i, k) * t_row_base_fun(k);
t_PRT(i, k) = t_pushed_base(j) *
(t_diff_R(i, l, k) * t_u(l, j) + t_diff_R(i, j, k));
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
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;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_diff_R;
++t_u;
++t_x_grad;
}
}
OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
int row_nb_base_functions = row_data.getN().size2() / 9;
auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
t_pushed_base(i, j) = t_x_grad(j, k) * t_row_base_fun(i, k);
t_PRT(k) = t_pushed_base(i, j) *
(t_diff_R(i, l, k) * t_u(l, j) + t_diff_R(i, j, k));
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
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;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_diff_R;
++t_x_grad;
++t_u;
}
}
MoFEMErrorCode OpSpatialConsistency_dP_dx::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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 v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
int row_nb_base_functions = row_data.getN().size2() / 3;
auto t_row_base_fun = row_data.getFTensor1N<3>();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
auto t_m = get_ftensor2(K, 3 * rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
const double b = a * t_col_diff_base_fun(l) * t_row_base_fun(l);
t_m(i, j) += b * t_delta_h(i, j);
t_m(0, 0) -= b;
t_m(1, 1) -= b;
t_m(2, 2) -= b;
++t_m;
++t_col_diff_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_delta_h;
}
}
MoFEMErrorCode OpSpatialConsistency_dBubble_dx::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
int row_nb_base_functions = row_data.getN().size2() / 9;
auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs; ++rr) {
auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
auto t_m = get_ftensor2(K, rr, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
FTensor::Tensor1<double, 3> t_pushed_base_dx;
t_pushed_base_dx(i) = a * t_col_diff_base_fun(l) * t_row_base_fun(i, l);
t_m(j) += t_pushed_base_dx(i) * t_delta_h(i, j);
t_m(i) -= t_pushed_base_dx(i);
++t_m;
++t_col_diff_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_delta_h;
}
}
MoFEMErrorCode OpSpatialPrj::integrate(EntData &row_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
auto get_ftensor1 = [](VectorDouble &v, const int r) {
&v[r + 0], &v[r + 1], &v[r + 2]);
};
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
auto t_spatial_disp = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
auto t_rhs = get_ftensor1(nF, 0);
double a = v * t_w * t_x_det;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
t_rhs(i) += a * t_row_base_fun * (-t_spatial_disp(i));
++t_row_base_fun;
++t_rhs;
}
for (; rr < row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_spatial_disp;
++t_w;
++t_x_det;
}
}
MoFEMErrorCode OpSpatialPrj_dx_dx::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
auto t_x_inv_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xInvGradAtPts));
auto t_spatial_disp = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const double ts_a = getFEMethod()->ts_a;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w * t_x_det;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(K, 3 * rr, 0);
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
t_val(i) = a * t_row_base_fun * (-t_spatial_disp(i));
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
t_m(i, j) += t_val(i) * t_x_inv_grad(k, j) * t_col_diff_base_fun(k);
++t_m;
++t_col_base_fun;
++t_col_diff_base_fun;
}
++t_row_base_fun;
}
for (; rr < row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_spatial_disp;
++t_w;
++t_x_det;
++t_x_inv_grad;
}
}
MoFEMErrorCode OpSpatialPrj_dx_dw::integrate(EntData &row_data,
EntData &col_data) {
int nb_integration_pts = row_data.getN().size1();
int row_nb_dofs = row_data.getIndices().size();
int col_nb_dofs = col_data.getIndices().size();
auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
&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)
);
};
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w * t_x_det;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_m = get_ftensor2(K, 3 * rr, 0);
auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
const double val = a * t_row_base_fun * t_col_base_fun;
t_m(0, 0) -= val;
t_m(1, 1) -= val;
t_m(2, 2) -= val;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr < row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
++t_x_det;
}
};
MoFEMErrorCode OpSpatialSchurBegin::assemble(int row_side, EntityType row_type,
EntData &data) {
if (row_type != MBTET)
auto &bmc = dataAtPts->blockMatContainor;
for (auto &bit : bmc)
bit.unSetAtElement();
// bmc.clear();
}
MoFEMErrorCode OpSpatialPreconditionMass::integrate(EntData &row_data) {
const int nb_integration_pts = row_data.getN().size1();
const int row_nb_dofs = row_data.getIndices().size();
auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
);
};
auto v = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto &m = *mPtr;
m.resize(row_nb_dofs, row_nb_dofs, false);
m.clear();
int row_nb_base_functions = row_data.getN().size2();
auto t_row_base_fun = row_data.getFTensor0N();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double a = v * t_w;
int rr = 0;
for (; rr != row_nb_dofs / 3; ++rr) {
auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
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(i) += b;
++t_m;
++t_col_base_fun;
}
++t_row_base_fun;
}
for (; rr != row_nb_base_functions; ++rr)
++t_row_base_fun;
++t_w;
}
}
MoFEMErrorCode OpSpatialSchurEnd::assemble(int row_side, EntityType row_type,
EntData &data) {
if (row_type != MBTET)
if (auto ep_fe_vol_ptr =
dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
getFEMethod())) {
AO aoSuu = PETSC_NULL;
AO aoSBB = PETSC_NULL;
AO aoSOO = PETSC_NULL;
AO aoSww = PETSC_NULL;
Mat Suu = PETSC_NULL;
Mat SBB = PETSC_NULL;
Mat SOO = PETSC_NULL;
Mat Sww = PETSC_NULL;
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 find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
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 set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
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->setMat(m);
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 add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
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->addMat(m);
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) {
const VectorInt &rind = bit.rowInd;
const VectorInt &cind = bit.colInd;
const MatrixDouble &m = bit.M;
CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
&*cind.begin(), &*m.data().begin(), ADD_VALUES);
};
auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
const int nb = m.size1();
inv.resize(nb, nb, false);
inv.clear();
for (int cc = 0; cc != nb; ++cc)
inv(cc, cc) = -1;
iPIV.resize(nb, false);
lapackWork.resize(nb * nb, false);
const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
&*iPIV.begin(), &*inv.data().begin(), nb,
&*lapackWork.begin(), nb * nb);
if (info != 0)
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"Can not invert matrix info = %d", info);
};
auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
std::string field, auto &inv) {
auto bit =
bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
if (bit != bmc.get<1>().end()) {
auto &m = *const_cast<MatrixDouble *>(&(bit->M));
CHKERR invert_symm_mat(m, inv);
} else
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"%s matrix not found", field.c_str());
};
auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
const int nb = m.size1();
MatrixDouble trans_m = trans(m);
MatrixDouble trans_inv;
trans_inv.resize(nb, nb, false);
trans_inv.clear();
for (int c = 0; c != nb; ++c)
trans_inv(c, c) = -1;
iPIV.resize(nb, false);
const auto info =
lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
&*trans_inv.data().begin(), nb);
if (info != 0)
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"Can not invert matrix info = %d", info);
inv.resize(nb, nb, false);
noalias(inv) = trans(trans_inv);
};
auto invert_nonsymm_schur =
[&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
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()) {
auto &m = *const_cast<MatrixDouble *>(&(bit->M));
CHKERR invert_nonsymm_mat(m, inv);
if (debug) {
std::cerr << prod(m, inv) << endl;
std::cerr << endl;
}
} else
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"%s matrix not found", field.c_str());
};
auto create_block_schur =
[&](DataAtIntegrationPts::BlockMatContainor &bmc,
DataAtIntegrationPts::BlockMatContainor &add_bmc,
std::string field, AO ao, MatrixDouble &inv) {
for (auto &bit : add_bmc) {
bit.unSetAtElement();
bit.clearMat();
}
for (auto &bit : bmc) {
if (bit.setAtElement && bit.rowField != field &&
bit.colField != field) {
VectorInt rind = bit.rowInd;
VectorInt cind = bit.colInd;
const MatrixDouble &m = bit.M;
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.rowSide, bit.colSide, bit.rowType,
bit.colType, m, rind, cind);
}
}
for (auto &bit_col : bmc) {
if (bit_col.setAtElement && bit_col.rowField == field &&
bit_col.colField != field) {
const MatrixDouble &cm = bit_col.M;
VectorInt cind = bit_col.colInd;
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) {
const MatrixDouble &rm = bit_row.M;
VectorInt rind = bit_row.rowInd;
K.resize(rind.size(), cind.size(), false);
noalias(K) = prod(rm, invMat);
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 =
[&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
bool debug = false) {
for (auto &bit : add_bmc) {
if (bit.setAtElement)
CHKERR assemble_block(bit, S);
}
if (debug) {
for (auto &bit : add_bmc) {
if (bit.setAtElement) {
std::cerr << "assemble: " << bit.rowField << " "
<< bit.colField << endl;
std::cerr << bit.M << endl;
}
}
std::cerr << std::endl;
}
};
auto precondition_schur =
[&](DataAtIntegrationPts::BlockMatContainor &bmc,
DataAtIntegrationPts::BlockMatContainor &add_bmc,
const std::string field, const MatrixDouble &diag_mat,
const double eps) {
for (auto &bit : add_bmc) {
bit.unSetAtElement();
bit.clearMat();
}
for (auto &bit : bmc) {
if (bit.setAtElement) {
if (bit.rowField != field || bit.colField != field)
auto it =
set_block(add_bmc, bit.rowField, bit.colField,
bit.rowSide, bit.colSide, bit.rowType,
bit.colType, bit.M, bit.rowInd, bit.colInd);
}
}
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,
bit->colSide, bit->rowType, bit->colType, bit->M,
bit->rowInd, bit->colInd);
MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
m += eps * diag_mat;
} 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);
MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
m *= eps;
break;
}
}
if (row_it == bmc.get<3>().end())
SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"row field not found %s", field.c_str());
}
};
CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
invBlockMat["uu"]);
CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
"u", aoSuu, invBlockMat["uu"]);
CHKERR assemble_schur(blockMat["uu"], Suu);
if (SBB) {
CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
aoSBB, invBlockMat["BB"]);
CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
*dataAtPts->ooMatPtr, eps);
CHKERR assemble_schur(blockMat["precBBOO"], SBB);
if (SOO) {
CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
invBlockMat["OO"]);
CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
"omega", aoSOO, invBlockMat["OO"]);
CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
*dataAtPts->wwMatPtr, eps);
CHKERR assemble_schur(blockMat["precOOww"], SOO);
if (Sww) {
CHKERR invert_symm_schur(blockMat["precOOww"], "w",
invBlockMat["ww"]);
CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
aoSww, invBlockMat["ww"]);
CHKERR assemble_schur(blockMat["ww"], Sww);
}
}
}
}
}
} // namespace EshelbianPlasticity
MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
EntData &data) {
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;
CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
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("Piola2Stress", 9);
Tag th_sigma = create_tag("CauchyStress", 9);
Tag th_u = create_tag("SpatialStreach", 9);
Tag th_h = create_tag("h", 9);
Tag th_x = create_tag("x", 3);
Tag th_X = create_tag("X", 3);
Tag th_detF = create_tag("detF", 1);
Tag th_angular_momentum = create_tag("AngularMomentum", 3);
int nb_gauss_pts = data.getN().size1();
if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Nb. of integration points is not equal to number points on "
"post-processing mesh");
}
auto t_w = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
auto t_x = getFTensor1FromMat<3>(*(dataAtPts->xAtPts));
auto t_omega = getFTensor1FromMat<3>(*(dataAtPts->rotAxisAtPts));
auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
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();
// vectors
FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
t_v(i) = t_d(i);
CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
&*v.data().begin());
};
&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) {
t_m(i, j) = t_d(i, j);
CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
&*m.data().begin());
};
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
// vetors
CHKERR save_vec_tag(th_w, t_w, gg);
CHKERR save_vec_tag(th_x, t_x, gg);
CHKERR save_vec_tag(th_X, t_coords, gg);
CHKERR save_vec_tag(th_omega, t_omega, gg);
// tensors
CHKERR save_mat_tag(th_h, t_h, 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 = dEterminant(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);
CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
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);
CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
&t_leviPRT(0));
++t_w;
++t_x;
++t_h;
++t_u;
++t_omega;
++t_R;
++t_approx_P;
++t_coords;
}
}
OpCalculateStrainEnergy::doWork(int side, EntityType type,
int nb_integration_pts = data.getN().size1();
auto v = getVolume();
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;
}
}
} // namespace EshelbianPlasticity