struct CommonData :
public boost::enable_shared_from_this<CommonData> {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
mD);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
mGrad);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactStress);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactDisp);
}
};
#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 &
sdf
) {
try {
sdf = bp::extract<double>(sdfFun(
t, x, y, z));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
MoFEMErrorCode evalGradSdf(
double t,
double x,
double y,
double z, std::vector<double> &grad_sdf
) {
try {
grad_sdf = py_list_to_std_vector<double>(sdfGradFun(
t, x, y, z));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
if (grad_sdf.size() != 3) {
}
}
MoFEMErrorCode evalHessSdf(
double t,
double x,
double y,
double z, std::vector<double> &hess_sdf
) {
try {
hess_sdf = py_list_to_std_vector<double>(sdfHessFun(
t, x, y, z));
} 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
OpConstrainBoundaryRhs(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &data);
private:
};
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:
};
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:
};
template <typename T>
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
sdf_ptr->evalSdf(
t, t_coords(0), t_coords(1), t_coords(2),
sdf),
"Failed python call");
}
#endif
return t_coords(1) + 0.5;
};
template <typename T>
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> grad_sdf;
t_coords(2), grad_sdf),
"Failed python call");
}
#endif
};
template <typename T>
#ifdef PYTHON_SFD
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> hess_sdf;
t_coords(2), hess_sdf),
"Failed python call");
hess_sdf[2], hess_sdf[3],
hess_sdf[4], hess_sdf[5]};
}
#endif
};
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; }
return (1 - s) / 2;
}
const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
MoFEMErrorCode
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 =
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();
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto t_grad_sdf =
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_cQ(
i,
j) * (t_disp(
j) -
cn * t_traction(
j))
+
c * (
sdf * t_grad_sdf(
i));
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)
commonDataPtr(common_data_ptr) {
sYmm = false;
}
MoFEMErrorCode
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 =
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();
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto t_grad_sdf =
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);
kronecker_delta(
i,
j) + t_cP(
i,
j);
auto t_hess_sdf =
(
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)) +
}
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
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)
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 =
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();
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
auto t_grad_sdf =
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_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>(
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;
}
}
};
#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()
const double c
speed of light (cm/ns)
constexpr double t
plate stiffness
constexpr auto field_name
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns