v0.14.0
ContactOps.hpp
/**
* \file ContactOps.hpp
* \example ContactOps.hpp
*/
#ifndef __CONTACTOPS_HPP__
#define __CONTACTOPS_HPP__
namespace ContactOps {
//! [Common data]
struct CommonData : public boost::enable_shared_from_this<CommonData> {
// MatrixDouble contactStress;
VectorDouble sdfVals; ///< size is equal to number of gauss points on element
MatrixDouble gradsSdf; ///< nb of rows is equals to dimension, and nb of cols
///< is equals to number of gauss points on element
MatrixDouble hessSdf; ///< nb of rows is equals to nb of element of symmetric
///< matrix, and nb of cols is equals to number of gauss
///< points on element
static SmartPetscObj<Vec>
totalTraction; // User have to release and create vector when appropiate.
static auto createTotalTraction(MoFEM::Interface &m_field) {
constexpr int ghosts[] = {0, 1, 2, 3, 4};
(m_field.get_comm_rank() == 0) ? 5 : 0, 5,
(m_field.get_comm_rank() == 0) ? 0 : 5, ghosts);
return totalTraction;
}
static auto getFTensor1TotalTraction() {
const double *t_ptr;
CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
"get array");
FTensor::Tensor1<double, 5> t{t_ptr[0], t_ptr[1], t_ptr[2], t_ptr[3],
t_ptr[4]};
CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
"restore array");
return t;
} else {
return FTensor::Tensor1<double, 5>{0., 0., 0., 0., 0.};
}
}
inline auto contactTractionPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto contactDispPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
}
inline auto contactDispGradPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto sdfPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &sdfVals);
}
inline auto gradSdfPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradsSdf);
}
inline auto hessSdfPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &hessSdf);
}
inline auto constraintPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &constraintVals);
}
};
SmartPetscObj<Vec> CommonData::totalTraction;
//! [Common data]
//! [Surface distance function from python]
#ifdef PYTHON_SDF
struct SDFPython {
SDFPython() = default;
virtual ~SDFPython() = default;
MoFEMErrorCode sdfInit(const std::string py_file) {
try {
// create main module
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
// create a reference to python function
sdfFun = mainNamespace["sdf"];
sdfGradFun = mainNamespace["grad_sdf"];
sdfHessFun = mainNamespace["hess_sdf"];
} catch (bp::error_already_set const &) {
// print all other errors to stderr
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 delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &sdf
) {
try {
// call python function
sdf = bp::extract<np::ndarray>(
sdfFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
}
MoFEMErrorCode evalGradSdf(
double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &grad_sdf
) {
try {
// call python function
grad_sdf = bp::extract<np::ndarray>(
sdfGradFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
}
MoFEMErrorCode evalHessSdf(
double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &hess_sdf
) {
try {
// call python function
hess_sdf = bp::extract<np::ndarray>(
sdfHessFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
}
private:
bp::object mainNamespace;
bp::object sdfFun;
bp::object sdfGradFun;
bp::object sdfHessFun;
};
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
int id) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(nb_gauss_pts);
auto stride = bp::make_tuple(3 * sizeof(double));
return (np::from_data(&data[id], dtype, size, stride, bp::object()));
};
#endif
//! [Surface distance function from python]
using SurfaceDistanceFunction = boost::function<VectorDouble(
double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
using GradSurfaceDistanceFunction = boost::function<MatrixDouble(
double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
using HessSurfaceDistanceFunction = boost::function<MatrixDouble(
double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
inline VectorDouble surface_distance_function(double delta_t, double t,
int nb_gauss_pts,
MatrixDouble &m_spatial_coords,
MatrixDouble &m_normals_at_pts,
int block_id) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
}
np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
np::dtype::get_builtin<double>());
CHK_MOAB_THROW(sdf_ptr->evalSdf(delta_t, t,
bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]),
block_id, np_sdf),
"Failed python call");
double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
VectorDouble v_sdf;
v_sdf.resize(nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
v_sdf[gg] = *(sdf_val_ptr + gg);
return v_sdf;
}
#endif
VectorDouble v_sdf;
v_sdf.resize(nb_gauss_pts, false);
auto t_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
v_sdf[gg] = -t_coords(2) - 0.1;
++t_coords;
}
return v_sdf;
}
grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
MatrixDouble &m_spatial_coords,
MatrixDouble &m_normals_at_pts, int block_id) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
}
np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
np::dtype::get_builtin<double>());
CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]), block_id,
np_grad_sdf),
"Failed python call");
double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
MatrixDouble m_grad_sdf;
m_grad_sdf.resize(3, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 3; ++idx)
m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
}
return m_grad_sdf;
}
#endif
MatrixDouble m_grad_sdf;
m_grad_sdf.resize(3, nb_gauss_pts, false);
FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 0.0, -1.0};
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
t_grad_sdf(i) = t_grad_sdf_set(i);
++t_grad_sdf;
}
return m_grad_sdf;
}
hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
MatrixDouble &m_spatial_coords,
MatrixDouble &m_normals_at_pts, int block_id) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
};
np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
np::dtype::get_builtin<double>());
CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]), block_id,
np_hess_sdf),
"Failed python call");
double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
MatrixDouble m_hess_sdf;
m_hess_sdf.resize(6, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 6; ++idx)
m_hess_sdf(idx, gg) =
*(hess_ptr + (6 * gg + idx));
}
return m_hess_sdf;
}
#endif
MatrixDouble m_hess_sdf;
m_hess_sdf.resize(6, nb_gauss_pts, false);
FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
++t_hess_sdf;
}
return m_hess_sdf;
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl;
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactAreaImpl;
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 <typename T1, typename T2, int DIM1, int DIM2>
size_t nb_gauss_pts) {
MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
m_spatial_coords.clear();
auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_spatial_coords(i) = t_coords(i) + t_disp(i);
++t_spatial_coords;
++t_coords;
++t_disp;
}
return m_spatial_coords;
}
template <typename T1, int DIM1>
inline auto get_normalize_normals(FTensor::Tensor1<T1, DIM1> &&t_normal_at_pts,
size_t nb_gauss_pts) {
MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
m_normals_at_pts.clear();
auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
++t_set_normal;
++t_normal_at_pts;
}
return m_normals_at_pts;
}
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>
: public BoundaryEleOp {
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
bool is_axisymmetric = false);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
const double scaleTraction;
bool isAxisymmetric;
};
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>
: public BoundaryEleOp {
OpAssembleTotalContactAreaImpl(
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric = false,
boost::shared_ptr<Range> contact_range_ptr = nullptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
boost::shared_ptr<Range> contactRange;
};
template <int DIM, typename BoundaryEleOp>
struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric = false);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct 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,
bool is_axisymmetric = false);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct 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,
bool is_axisymmetric = false);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <typename BoundaryEleOp> struct ContactIntegrators {
template <int DIM, IntegrationType I>
OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
OpAssembleTotalContactAreaImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
using OpEvaluateSDF = OpEvaluateSDFImpl<DIM, I, BoundaryEleOp>;
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) {
constexpr auto eps = std::numeric_limits<float>::epsilon();
if (std::abs(x) < eps)
return 0;
else if (x > eps)
return 1;
else
return -1;
};
inline double w(const double sdf, const double tn) {
return sdf - cn_contact * tn;
}
/**
* @brief constrain function
*
* return 1 if negative sdf or positive tn
*
* @param sdf signed distance
* @param tn traction
* @return double
*/
inline double constrain(double sdf, double tn) {
const auto s = sign(w(sdf, tn));
return (1 - s) / 2;
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr, double scale,
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
commonDataPtr(common_data_ptr), scaleTraction(scale),
isAxisymmetric(is_axisymmetric) {}
template <int DIM, typename BoundaryEleOp>
int side, EntityType type, EntData &data) {
FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
t_sum_t(i) += alpha * t_traction(i);
++t_w;
++t_traction;
++t_coords;
}
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>
OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactAreaImpl(
boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric,
boost::shared_ptr<Range> contact_range_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric),
contactRange(contact_range_ptr) {}
template <int DIM, typename BoundaryEleOp>
int side, EntityType type, EntData &data) {
auto fe_type = BoundaryEleOp::getFEType();
const auto fe_ent = BoundaryEleOp::getFEEntityHandle();
if (contactRange->find(fe_ent) != contactRange->end()) {
FTensor::Tensor1<double, 2> t_sum_a{0., 0.};
auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->contactDispGrad);
auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = BoundaryEleOp::getTStime();
auto ts_time_step = BoundaryEleOp::getTStimeStep();
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf = gradSurfaceDistanceFunction(
ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0); // Axisymmetric Jacobian
}
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
double alpha = t_w * jacobian;
F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
auto det = determinantTensor(F);
CHKERR invertTensor(F, det, invF);
t_normal_current(i) = det * (invF(j, i) * t_normal_at_pts(j));
alpha *= sqrt(t_normal_current(i) * t_normal_current(i));
if (fe_type == MBTRI) {
alpha /= 2;
}
if (c > 1e-12) {
t_sum_a(0) += alpha; // real area
}
t_sum_a(1) += alpha; // Potential area
++t_w;
++t_traction;
++t_coords;
++t_sdf;
++t_grad_sdf;
++t_grad;
++t_normal_at_pts;
}
constexpr int ind[] = {3, 4};
CHKERR VecSetValues(commonDataPtr->totalTraction, 2, ind, &t_sum_a(0),
ADD_VALUES);
}
}
template <int DIM, typename BoundaryEleOp>
boost::shared_ptr<CommonData> common_data_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
commonDataPtr(common_data_ptr) {}
template <int DIM, typename BoundaryEleOp>
EntData &data) {
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
auto &sdf_vec = commonDataPtr->sdfVals;
auto &grad_mat = commonDataPtr->gradsSdf;
auto &hess_mat = commonDataPtr->hessSdf;
auto &constraint_vec = commonDataPtr->constraintVals;
auto &contactTraction_mat = commonDataPtr->contactTraction;
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);
constraint_vec.resize(nb_gauss_pts, false);
auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
auto t_sdf = getFTensor0FromVec(sdf_vec);
auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
auto t_constraint = getFTensor0FromVec(constraint_vec);
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto ts_time = BoundaryEleOp::getTStime();
auto ts_time_step = BoundaryEleOp::getTStimeStep();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_hess_sdf =
hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf_v = getFTensor0FromVec(v_sdf);
auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
auto next = [&]() {
++t_sdf;
++t_sdf_v;
++t_grad_sdf;
++t_grad_sdf_v;
++t_hess_sdf;
++t_hess_sdf_v;
++t_disp;
++t_traction;
++t_constraint;
};
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
auto tn = -t_traction(i) * t_grad_sdf_v(i);
auto c = constrain(t_sdf_v, tn);
t_sdf = t_sdf_v;
t_grad_sdf(i) = t_grad_sdf_v(i);
t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
t_constraint = c;
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), isAxisymmetric(is_axisymmetric) {}
template <int DIM, typename AssemblyBoundaryEleOp>
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>();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
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_rhs(i) =
t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
+
t_cP(i, j) * t_disp(j) +
c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
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;
++t_sdf;
++t_grad_sdf;
}
}
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(row_field_name, col_field_name,
AssemblyBoundaryEleOp::OPROWCOL),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
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;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_hess_sdf =
hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
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);
if (c > 0) {
t_res_dU(i, j) +=
(c * cn_contact) *
(t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
c * t_sdf * t_hess_sdf(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.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_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
++t_hess_sdf;
}
}
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, bool is_axisymmetric)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
AssemblyBoundaryEleOp::OPROWCOL),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
iNtegrate(EntitiesFieldData::EntData &row_data,
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
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;
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
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_dt(i, j) = -cn_contact * t_cQ(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_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
}
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u, bool is_axisymmetric = false) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
using OpMixDivUCylRhs =
typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
using OpMixUTimesDivLambdaRhs =
typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
using OpMixUTimesLambdaRhs =
typename B::template OpGradTimesTensor<1, DIM, 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>();
auto jacobian = [is_axisymmetric](const double r, const double,
const double) {
return 2. * M_PI * r;
else
return 1.;
};
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));
} else {
pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
sigma, div_stress_ptr));
}
pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
pip.push_back(
new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
} else {
pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
jacobian));
}
pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
}
template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
using OpMixLhs::OpMixLhs;
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
// Only assemble side which correspond to edge entity on boundary
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>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string fe_domain_name, std::string sigma, std::string u,
std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
bool is_axisymmetric = false) {
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<
A>::template BiLinearForm<I>;
using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
using OpMixDivUCylLhs =
typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
auto unity = []() { return 1; };
auto jacobian = [is_axisymmetric](const double r, const double,
const double) {
return 2. * M_PI * r;
else
return 1.;
};
op_loop_side->getOpPtrVector().push_back(
new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
} else {
op_loop_side->getOpPtrVector().push_back(
new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
}
op_loop_side->getOpPtrVector().push_back(
new OpLambdaGraULhsSide(sigma, u, unity, jacobian, 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, bool is_axisymmetric = false) {
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, is_axisymmetric));
pip.push_back(new typename C::template Assembly<A>::
template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
sigma, sigma, common_data_ptr, is_axisymmetric));
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u, bool is_axisymmetric = false) {
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, is_axisymmetric));
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, bool is_axisymmetric = false) {
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, is_axisymmetric));
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
OpLoopSide<SideEle> *op_loop_side, std::string sigma, std::string u,
bool is_axisymmetric = false,
boost::shared_ptr<Range> contact_range_ptr = nullptr) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
op_loop_side->getOpPtrVector().push_back(
"U", common_data_ptr->contactDispGradPtr()));
if (contact_range_ptr) {
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(op_loop_side);
pip.push_back(new typename C::template OpAssembleTotalContactArea<DIM, I>(
common_data_ptr, is_axisymmetric, contact_range_ptr));
}
}
}; // namespace ContactOps
#endif // __CONTACTOPS_HPP__
ContactOps::get_spatial_coords
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
Definition: ContactOps.hpp:414
NOSPACE
@ NOSPACE
Definition: definitions.h:83
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
ContactOps::CommonData::contactTractionPtr
auto contactTractionPtr()
Definition: ContactOps.hpp:57
OpMixLhs
ContactOps::sign
double sign(double x)
Definition: ContactOps.hpp:580
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
ContactOps::CommonData::constraintVals
VectorDouble constraintVals
Definition: ContactOps.hpp:26
ContactOps::CommonData::contactDisp
MatrixDouble contactDisp
Definition: ContactOps.hpp:17
sdf.hess_sdf
def hess_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:19
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
ContactOps::hess_surface_distance_function
MatrixDouble hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:338
ContactOps::w
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:590
ContactOps
Definition: contact.cpp:99
ContactOps::constrain
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:603
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:93
ContactOps::OpEvaluateSDFImpl< DIM, GAUSS, BoundaryEleOp >::OpEvaluateSDFImpl
OpEvaluateSDFImpl(boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:745
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
ContactOps::CommonData::contactDispGrad
MatrixDouble contactDispGrad
Definition: ContactOps.hpp:18
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
ContactOps::OpConstrainBoundaryRhsImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:849
ContactOps::ContactIntegrators::Assembly::OpConstrainBoundaryLhs_dU
OpConstrainBoundaryLhs_dUImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryLhs_dU
Definition: ContactOps.hpp:572
sdf
Definition: sdf.py:1
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ContactOps::CommonData::contactDispGradPtr
auto contactDispGradPtr()
Definition: ContactOps.hpp:66
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
ContactOps::SurfaceDistanceFunction
boost::function< VectorDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> SurfaceDistanceFunction
[Common data]
Definition: ContactOps.hpp:213
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ContactOps.hpp:619
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
ContactOps::HessSurfaceDistanceFunction
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> HessSurfaceDistanceFunction
Definition: ContactOps.hpp:221
ContactOps::ContactIntegrators::OpAssembleTotalContactArea
OpAssembleTotalContactAreaImpl< DIM, I, BoundaryEleOp > OpAssembleTotalContactArea
Definition: ContactOps.hpp:556
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
ContactOps::opFactoryBoundaryToDomainLhs
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1265
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
ContactOps::ContactIntegrators::OpAssembleTotalContactTraction
OpAssembleTotalContactTractionImpl< DIM, I, BoundaryEleOp > OpAssembleTotalContactTraction
Definition: ContactOps.hpp:552
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:23
MoFEM::createGhostVector
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Definition: PetscSmartObj.hpp:179
BoundaryEleOp
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:29
ContactOps::ContactIntegrators::Assembly::OpConstrainBoundaryLhs_dTraction
OpConstrainBoundaryLhs_dTractionImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryLhs_dTraction
Definition: ContactOps.hpp:576
ContactOps::get_normalize_normals
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
Definition: ContactOps.hpp:431
double
ContactOps::opFactoryBoundaryRhs
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1343
convert.type
type
Definition: convert.py:64
ContactOps::opFactoryCalculateArea
MoFEMErrorCode opFactoryCalculateArea(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, OpLoopSide< SideEle > *op_loop_side, std::string sigma, std::string u, bool is_axisymmetric=false, boost::shared_ptr< Range > contact_range_ptr=nullptr)
Definition: ContactOps.hpp:1381
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::invertTensor
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition: Templates.hpp:1771
MoFEM::determinantTensor
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition: Templates.hpp:1571
ContactOps::CommonData::hessSdf
MatrixDouble hessSdf
Definition: ContactOps.hpp:23
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
ContactOps::OpAssembleTotalContactAreaImpl< DIM, GAUSS, BoundaryEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ContactOps.hpp:662
ContactOps::CommonData::gradSdfPtr
auto gradSdfPtr()
Definition: ContactOps.hpp:75
ContactOps::AssemblyBoundaryEleOp
FormsIntegrators< BoundaryEleOp >::Assembly< A >::OpBase AssemblyBoundaryEleOp
Definition: EshelbianContact.hpp:18
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
ContactOps::CommonData::sdfVals
VectorDouble sdfVals
size is equal to number of gauss points on element
Definition: ContactOps.hpp:20
ContactOps::CommonData::gradsSdf
MatrixDouble gradsSdf
Definition: ContactOps.hpp:21
ContactOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1186
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
ContactOps::CommonData::sdfPtr
auto sdfPtr()
Definition: ContactOps.hpp:71
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
ContactOps::CommonData::contactDispPtr
auto contactDispPtr()
Definition: ContactOps.hpp:62
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
ContactOps::CommonData::createTotalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:31
ContactOps::cn_contact
double cn_contact
Definition: contact.cpp:100
ContactOps::GradSurfaceDistanceFunction
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> GradSurfaceDistanceFunction
Definition: ContactOps.hpp:217
ContactOps::BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: EshelbianContact.hpp:14
DomainEleOp
ContactOps::OpEvaluateSDFImpl< DIM, GAUSS, BoundaryEleOp >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ContactOps.hpp:752
ContactOps::opFactoryBoundaryLhs
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1319
ContactOps::CommonData::constraintPtr
auto constraintPtr()
Definition: ContactOps.hpp:83
ContactOps::CommonData::contactTraction
MatrixDouble contactTraction
Definition: ContactOps.hpp:16
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition: Kronecker_Delta.hpp:81
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
ContactOps::CommonData::getFTensor1TotalTraction
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:42
ContactOps::ContactIntegrators::Assembly::OpConstrainBoundaryRhs
OpConstrainBoundaryRhsImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryRhs
Definition: ContactOps.hpp:568
ContactOps::surface_distance_function
VectorDouble surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:223
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
AssemblyBoundaryEleOp
sdf_wavy_2d.ind
float ind
Definition: sdf_wavy_2d.py:7
ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:962
ContactOps::OpMixLhsSide::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:1248
ContactOps::grad_surface_distance_function
MatrixDouble grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:281
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
ContactOps::opFactoryCalculateTraction
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1364
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
ContactOps::ContactIntegrators::OpEvaluateSDF
OpEvaluateSDFImpl< DIM, I, BoundaryEleOp > OpEvaluateSDF
Definition: ContactOps.hpp:559
ContactOps::CommonData::hessSdfPtr
auto hessSdfPtr()
Definition: ContactOps.hpp:79
sdf.grad_sdf
def grad_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:15
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
F
@ F
Definition: free_surface.cpp:394