#ifndef __CONTACTOPS_HPP__
#define __CONTACTOPS_HPP__
struct CommonData :
public boost::enable_shared_from_this<CommonData> {
static SmartPetscObj<Vec>
constexpr int ghosts[] = {0, 1, 2};
}
const double *t_ptr;
"get array");
"restore array");
} else {
}
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactDisp);
}
return boost::shared_ptr<VectorDouble>(shared_from_this(), &
sdfVals);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
gradsSdf);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
hessSdf);
}
};
#ifdef PYTHON_SFD
struct SDFPython {
SDFPython() = default;
virtual ~SDFPython() = default;
MoFEMErrorCode sdfInit(const std::string py_file) {
try {
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
sdfFun = mainNamespace["sdf"];
sdfGradFun = mainNamespace["grad_sdf"];
sdfHessFun = mainNamespace["hess_sdf"];
} catch (bp::error_already_set const &) {
PyErr_Print();
}
};
template <typename T>
inline std::vector<T>
py_list_to_std_vector(const boost::python::object &iterable) {
return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
boost::python::stl_input_iterator<T>());
}
MoFEMErrorCode evalSdf(
double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
) {
try {
sdf = bp::extract<double>(sdfFun(
t, x, y, z, tx, ty, tz));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
MoFEMErrorCode evalGradSdf(
double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
std::vector<double> &grad_sdf
) {
try {
grad_sdf =
py_list_to_std_vector<double>(sdfGradFun(
t, x, y, z, tx, ty, tz));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
if (grad_sdf.size() != 3) {
}
}
MoFEMErrorCode evalHessSdf(
double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
std::vector<double> &hess_sdf
) {
try {
hess_sdf =
py_list_to_std_vector<double>(sdfHessFun(
t, x, y, z, tx, ty, tz));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
if (hess_sdf.size() != 6) {
}
}
private:
bp::object mainNamespace;
bp::object sdfFun;
bp::object sdfGradFun;
bp::object sdfHessFun;
};
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
#endif
double t,
double x,
double y,
double z,
double tx,
double ty,
double tz)>;
double t,
double x,
double y,
double z,
double tx,
double ty,
double tz)>;
boost::function<FTensor::Tensor2_symmetric<double, 3>(
double t,
double x,
double y,
double z,
double tx,
double ty,
double tz)>;
double tx, double ty, double tz) {
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
"Failed python call");
}
#endif
return y + 0.5;
}
double tx, double ty, double tz) {
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> grad_sdf;
"Failed python call");
}
#endif
}
double tx, double ty, double tz) {
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> hess_sdf;
"Failed python call");
hess_sdf[2], hess_sdf[3],
hess_sdf[4], hess_sdf[5]};
}
#endif
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl;
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpEvaluateSDFImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl;
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl<DIM, GAUSS,
BoundaryEleOp>
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr,
double scale = 1);
private:
boost::shared_ptr<CommonData> commonDataPtr;
const double scaleTraction;
};
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <typename BoundaryEleOp> struct ContactIntegrators {
template <int DIM, IntegrationType I>
OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
template <AssemblyType A> struct Assembly {
typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
};
};
inline double sign(
double x) {
if (x == 0)
return 0;
else if (x > 0)
return 1;
else
return -1;
};
inline double w(
const double sdf,
const double tn) {
}
return (1 - s) / 2;
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr,
double scale)
commonDataPtr(common_data_ptr), scaleTraction(
scale) {}
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_sum_t(
i) += alpha * t_traction(
i);
++t_w;
++t_traction;
}
t_sum_t(
i) *= scaleTraction;
constexpr int ind[] = {0, 1, 2};
CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
ADD_VALUES);
}
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
auto &sdf_vec = commonDataPtr->sdfVals;
auto &grad_mat = commonDataPtr->gradsSdf;
auto &hess_mat = commonDataPtr->hessSdf;
sdf_vec.resize(nb_gauss_pts, false);
grad_mat.resize(DIM, nb_gauss_pts, false);
hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
auto t_sdf = getFTensor0FromVec(sdf_vec);
auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto next = [&]() {
++t_sdf;
++t_grad_sdf;
++t_hess_sdf;
++t_disp;
++t_coords;
};
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto sdf_v = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf_v = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_hess_sdf_v = hessSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
t_sdf = sdf_v;
t_grad_sdf(
i) = t_grad_sdf_v(
i);
t_hess_sdf(
i,
j) = t_hess_sdf_v(
i,
j);
next();
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr)
AssemblyBoundaryEleOp::OPROW),
commonDataPtr(common_data_ptr) {}
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
size_t nb_base_functions = data.getN().size2() / 3;
auto t_base = data.getFTensor1N<3>();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_normal_at_pts(
i) / std::sqrt(t_normal_at_pts(
i) * t_normal_at_pts(
i));
auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto sdf = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto tn = -t_traction(
i) * t_grad_sdf(
i);
t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
+
c * (
sdf * t_grad_sdf(
i));
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
const double beta = alpha * (t_base(
i) * t_normal(
i));
t_nf(
i) -= beta * t_rhs(
i);
++t_nf;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
++t_normal_at_pts;
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
AssemblyBoundaryEleOp::OPROWCOL),
commonDataPtr(common_data_ptr) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_normal_at_pts(
i) / std::sqrt(t_normal_at_pts(
i) * t_normal_at_pts(
i));
const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto sdf = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto tn = -t_traction(
i) * t_grad_sdf(
i);
t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
t_res_dU(
i,
j) = kronecker_delta(
i,
j) + t_cP(
i,
j);
auto t_hess_sdf = hessSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1),
t_spatial_coords(2), t_total_traction(0), t_total_traction(1),
t_total_traction(2));
(t_hess_sdf(
i,
j) * (t_grad_sdf(
k) * t_traction(
k)) +
t_grad_sdf(
i) * t_hess_sdf(
k,
j) * t_traction(
k)) +
}
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
const double row_base = t_row_base(
i) * t_normal(
i);
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
const double beta = alpha * row_base * t_col_base;
t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
++t_normal_at_pts;
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryLhs_dTractionImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
AssemblyBoundaryEleOp::OPROWCOL),
commonDataPtr(common_data_ptr) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_normal_at_pts(
i) / std::sqrt(t_normal_at_pts(
i) * t_normal_at_pts(
i));
const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto sdf = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto tn = -t_traction(
i) * t_grad_sdf(
i);
t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
const double row_base = t_row_base(
i) * t_normal(
i);
auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
const double col_base = t_col_base(
i) * t_normal(
i);
const double beta = alpha * row_base * col_base;
t_mat(
i,
j) -= beta * t_res_dt(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
++t_normal_at_pts;
}
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u) {
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
using OpMixDivURhs =
typename B::template OpMixDivTimesU<3, DIM, DIM>;
using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto div_stress_ptr = boost::make_shared<MatrixDouble>();
auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(
new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
pip.push_back(
new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
pip.push_back(
[](double, double, double) constexpr { return 1; }));
pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr));
}
template <
typename OpMixLhs>
struct OpMixLhsSide :
public OpMixLhs {
using OpMixLhs::OpMixLhs;
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
if (side_fe_entity == side_fe_data) {
CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
col_data);
}
}
};
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string fe_domain_name, std::string sigma, std::string u,
std::string geom, ForcesAndSourcesCore::RuleHookFun rule) {
using DomainEleOp =
typename DomainEle::UserDataOperator;
auto op_loop_side = new OpLoopSide<DomainEle>(
m_field, fe_domain_name, DIM, Sev::noisy,
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
pip.push_back(op_loop_side);
CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
{H1, HDIV}, geom);
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
using OpMixDivULhs =
typename B::template OpMixDivTimesVec<DIM>;
using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
auto unity = []() { return 1; };
op_loop_side->getOpPtrVector().push_back(
new OpMixDivULhsSide(sigma, u, unity, true));
op_loop_side->getOpPtrVector().push_back(
new OpLambdaGraULhsSide(sigma, u, unity, true));
op_loop_side->getSideFEPtr()->getRuleHook = rule;
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(
new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
DIM, GAUSS>(sigma, u, common_data_ptr));
pip.push_back(new typename C::template Assembly<A>::
template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
sigma, sigma, common_data_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(
new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
DIM, GAUSS>(sigma, common_data_ptr));
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
common_data_ptr, 1. /
scale));
}
};
#endif
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ 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.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr double t
plate stiffness
constexpr auto field_name
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element