v0.13.2
Loading...
Searching...
No Matches
ContactOps.hpp
/**
* \file ContactOps.hpp
* \example ContactOps.hpp
*/
namespace ContactOps {
//! [Common data]
struct CommonData : public boost::enable_shared_from_this<CommonData> {
MatrixDouble mD;
MatrixDouble mGrad;
MatrixDouble contactStress;
MatrixDouble contactTraction;
MatrixDouble contactDisp;
inline auto mDPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
}
inline auto mGradPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mGrad);
}
inline auto contactStressPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactStress);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto contactTractionPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto contactDispPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
}
};
//! [Common data]
//! [Surface distance function from python]
#ifdef PYTHON_SFD
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 t, double x, double y, double z, double &sdf
) {
try {
// call python function
sdf = bp::extract<double>(sdfFun(t, x, y, z));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
}
MoFEMErrorCode evalGradSdf(
double t, double x, double y, double z, std::vector<double> &grad_sdf
) {
try {
// call python function
grad_sdf = py_list_to_std_vector<double>(sdfGradFun(t, x, y, z));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
if (grad_sdf.size() != 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 3");
}
}
MoFEMErrorCode evalHessSdf(
double t, double x, double y, double z, std::vector<double> &hess_sdf
) {
try {
// call python function
hess_sdf = py_list_to_std_vector<double>(sdfHessFun(t, x, y, z));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
if (hess_sdf.size() != 6) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 6");
}
}
private:
bp::object mainNamespace;
bp::object sdfFun;
bp::object sdfGradFun;
bp::object sdfHessFun;
};
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
#endif
//! [Surface distance function from python]
struct OpConstrainBoundaryRhs : public AssemblyBoundaryEleOp {
OpConstrainBoundaryRhs(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpConstrainBoundaryLhs_dU : public AssemblyBoundaryEleOp {
OpConstrainBoundaryLhs_dU(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;
};
struct OpConstrainBoundaryLhs_dTraction : public AssemblyBoundaryEleOp {
OpConstrainBoundaryLhs_dTraction(
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 T>
inline double surface_distance_function(double t,
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
double sdf;
sdf_ptr->evalSdf(t, t_coords(0), t_coords(1), t_coords(2), sdf),
"Failed python call");
return sdf;
}
#endif
return t_coords(1) + 0.5;
};
template <typename T>
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> grad_sdf;
CHK_MOAB_THROW(sdf_ptr->evalGradSdf(t, t_coords(0), t_coords(1),
t_coords(2), grad_sdf),
"Failed python call");
return FTensor::Tensor1<double, 3>{grad_sdf[0], grad_sdf[1], grad_sdf[2]};
}
#endif
return FTensor::Tensor1<double, 3>{0., 1., 0.};
};
template <typename T>
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> hess_sdf;
CHK_MOAB_THROW(sdf_ptr->evalHessSdf(t, t_coords(0), t_coords(1),
t_coords(2), hess_sdf),
"Failed python call");
return FTensor::Tensor2_symmetric<double, 3>{hess_sdf[0], hess_sdf[1],
hess_sdf[2], hess_sdf[3],
hess_sdf[4], hess_sdf[5]};
}
#endif
return FTensor::Tensor2_symmetric<double, 3>{0., 0., 0., 0., 0., 0.};
};
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 sdf - cn * tn; }
inline double constrain(double sdf, double tn) {
const auto s = sign(w(sdf, tn));
return (1 - s) / 2;
}
const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
MoFEMErrorCode
OpConstrainBoundaryRhs::iNtegrate(EntitiesFieldData::EntData &data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_w = getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactDisp);
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
auto t_coords = 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) {
auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
const double alpha = t_w * getMeasure();
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto sdf = surface_distance_function(getTStime(), t_spatial_coords);
auto t_grad_sdf =
grad_surface_distance_function(getTStime(), t_spatial_coords);
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(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 * t_traction(j))
+
t_cP(i, j) * t_disp(j) +
c * (sdf * t_grad_sdf(i)); // add gap0 displacements
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows / SPACE_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;
}
}
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,
commonDataPtr(common_data_ptr) {
sYmm = false;
}
MoFEMErrorCode
OpConstrainBoundaryLhs_dU::iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactDisp);
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = 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>();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto sdf = surface_distance_function(getTStime(), t_spatial_coords);
auto t_grad_sdf =
grad_surface_distance_function(getTStime(), t_spatial_coords);
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(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) {
auto t_hess_sdf =
hess_surface_distance_function(getTStime(), t_spatial_coords);
t_res_dU(i, j) +=
(c * cn) * (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 * sdf * t_hess_sdf(i, j);
}
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
locMat, SPACE_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 / SPACE_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;
}
}
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,
commonDataPtr(common_data_ptr) {
sYmm = false;
}
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactDisp);
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = 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) {
const double alpha = t_w * getMeasure();
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto sdf = surface_distance_function(getTStime(), t_spatial_coords);
auto t_grad_sdf =
grad_surface_distance_function(getTStime(), t_spatial_coords);
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(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 * t_cQ(i, j);
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
locMat, SPACE_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 / SPACE_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;
}
}
}; // namespace ContactOps
constexpr int SPACE_DIM
Kronecker Delta class.
FormsIntegrators< BoundaryEleOp >::Assembly< A >::OpBase AssemblyBoundaryEleOp
Definition: contact.cpp:53
double cn
Definition: contact.cpp:124
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
constexpr auto t_kd
const double c
speed of light (cm/ns)
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:163
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:162
FTensor::Tensor2_symmetric< double, 3 > hess_surface_distance_function(double t, FTensor::Tensor1< T, 3 > &t_coords)
Definition: ContactOps.hpp:229
double surface_distance_function(double t, FTensor::Tensor1< T, 3 > &t_coords)
Definition: ContactOps.hpp:198
FTensor::Tensor1< double, 3 > grad_surface_distance_function(double t, FTensor::Tensor1< T, 3 > &t_coords)
Definition: ContactOps.hpp:214
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:161
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:164
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:253
double constrain(double sdf, double tn)
Definition: ContactOps.hpp:255
double sign(double x)
Definition: ContactOps.hpp:244
Definition: sdf.py:1
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
MatrixDouble contactDisp
Definition: ContactOps.hpp:18
MatrixDouble contactTraction
Definition: ContactOps.hpp:17
auto contactStressDivergencePtr()
Definition: ContactOps.hpp:32
MatrixDouble contactStressDivergence
Definition: ContactOps.hpp:16
MatrixDouble contactStress
Definition: ContactOps.hpp:15
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:429
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:194
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:438
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:183
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:333
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:343
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:260
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:172
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:266
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns