11struct CommonData :
public boost::enable_shared_from_this<CommonData> {
21 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
mD);
25 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
mGrad);
29 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactStress);
33 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
38 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
43 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactDisp);
51 SDFPython() =
default;
52 virtual ~SDFPython() =
default;
54 MoFEMErrorCode sdfInit(
const std::string py_file) {
59 auto main_module = bp::import(
"__main__");
60 mainNamespace = main_module.attr(
"__dict__");
61 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
63 sdfFun = mainNamespace[
"sdf"];
64 sdfGradFun = mainNamespace[
"grad_sdf"];
65 sdfHessFun = mainNamespace[
"hess_sdf"];
67 }
catch (bp::error_already_set
const &) {
77 py_list_to_std_vector(
const boost::python::object &iterable) {
78 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
79 boost::python::stl_input_iterator<T>());
84 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
92 sdf = bp::extract<double>(sdfFun(
t, x, y, z, tx, ty, tz));
94 }
catch (bp::error_already_set
const &) {
104 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
105 std::vector<double> &grad_sdf
113 py_list_to_std_vector<double>(sdfGradFun(
t, x, y, z, tx, ty, tz));
115 }
catch (bp::error_already_set
const &) {
130 double t,
double x,
double y,
double z,
double tx,
double ty,
double tz,
131 std::vector<double> &hess_sdf
139 py_list_to_std_vector<double>(sdfHessFun(
t, x, y, z, tx, ty, tz));
141 }
catch (bp::error_already_set
const &) {
155 bp::object mainNamespace;
157 bp::object sdfGradFun;
158 bp::object sdfHessFun;
161static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
173 boost::shared_ptr<CommonData> common_data_ptr);
174 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &data);
182 const std::string col_field_name,
183 boost::shared_ptr<CommonData> common_data_ptr);
184 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &row_data,
185 EntitiesFieldData::EntData &col_data);
193 const std::string row_field_name,
const std::string col_field_name,
194 boost::shared_ptr<CommonData> common_data_ptr);
195 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &row_data,
196 EntitiesFieldData::EntData &col_data);
202template <
typename T1,
typename T2>
207 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
209 CHK_MOAB_THROW(sdf_ptr->evalSdf(
t, t_coords(0), t_coords(1), t_coords(2),
210 t_traction(0), t_traction(1),
212 "Failed python call");
216 return t_coords(1) + 0.5;
219template <
typename T1,
typename T2>
224 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
225 std::vector<double> grad_sdf;
227 sdf_ptr->evalGradSdf(
t, t_coords(0), t_coords(1), t_coords(2),
228 t_traction(0), t_traction(1),
229 (
SPACE_DIM == 3) ? t_traction(2) : 0, grad_sdf),
230 "Failed python call");
237template <
typename T1,
typename T2>
242 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
243 std::vector<double> hess_sdf;
245 sdf_ptr->evalHessSdf(
t, t_coords(0), t_coords(1), t_coords(2),
246 t_traction(0), t_traction(1),
247 (
SPACE_DIM == 3) ? t_traction(2) : 0, hess_sdf),
248 "Failed python call");
250 hess_sdf[2], hess_sdf[3],
251 hess_sdf[4], hess_sdf[5]};
266inline double w(
const double sdf,
const double tn) {
return sdf -
cn * tn; }
274 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
276 commonDataPtr(common_data_ptr) {}
282 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
284 auto &nf = AssemblyBoundaryEleOp::locF;
286 auto t_normal = getFTensor1Normal();
287 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
289 auto t_w = getFTensor0IntegrationWeight();
290 auto t_disp = getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactDisp);
292 getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactTraction);
293 auto t_coords = getFTensor1CoordsAtGaussPts();
295 size_t nb_base_functions = data.getN().size2() / 3;
296 auto t_base = data.getFTensor1N<3>();
297 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
299 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
300 const double alpha = t_w * getMeasure();
303 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
308 getTStime(), t_spatial_coords, t_traction);
309 auto tn = -t_traction(
i) * t_grad_sdf(
i);
313 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
315 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
320 t_cQ(
i,
j) * (t_disp(
j) -
cn * t_traction(
j))
324 t_cP(
i,
j) * t_disp(
j) +
325 c * (
sdf * t_grad_sdf(
i));
328 for (; bb != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++bb) {
329 const double beta = alpha * (t_base(
i) * t_normal(
i));
330 t_nf(
i) -= beta * t_rhs(
i);
335 for (; bb < nb_base_functions; ++bb)
348 const std::string row_field_name,
const std::string col_field_name,
349 boost::shared_ptr<CommonData> common_data_ptr)
352 commonDataPtr(common_data_ptr) {
358 EntitiesFieldData::EntData &col_data) {
361 const size_t nb_gauss_pts = getGaussPts().size2();
362 auto &locMat = AssemblyBoundaryEleOp::locMat;
364 auto t_normal = getFTensor1Normal();
365 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
367 auto t_disp = getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactDisp);
369 getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactTraction);
370 auto t_coords = getFTensor1CoordsAtGaussPts();
372 auto t_w = getFTensor0IntegrationWeight();
373 auto t_row_base = row_data.getFTensor1N<3>();
374 size_t nb_face_functions = row_data.getN().size2() / 3;
378 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
380 const double alpha = t_w * getMeasure();
383 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
388 getTStime(), t_spatial_coords, t_traction);
390 auto tn = -t_traction(
i) * t_grad_sdf(
i);
394 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
396 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
400 kronecker_delta(
i,
j) + t_cP(
i,
j);
404 getTStime(), t_spatial_coords, t_traction);
406 (
c *
cn) * (t_hess_sdf(
i,
j) * (t_grad_sdf(
k) * t_traction(
k)) +
407 t_grad_sdf(
i) * t_hess_sdf(
k,
j) * t_traction(
k)) +
408 c *
sdf * t_hess_sdf(
i,
j);
412 for (; rr != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++rr) {
414 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
417 const double row_base = t_row_base(
i) * t_normal(
i);
419 auto t_col_base = col_data.getFTensor0N(gg, 0);
420 for (
size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols /
SPACE_DIM;
422 const double beta = alpha * row_base * t_col_base;
424 t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
432 for (; rr < nb_face_functions; ++rr)
445 const std::string row_field_name,
const std::string col_field_name,
446 boost::shared_ptr<CommonData> common_data_ptr)
449 commonDataPtr(common_data_ptr) {
454 EntitiesFieldData::EntData &row_data,
455 EntitiesFieldData::EntData &col_data) {
458 const size_t nb_gauss_pts = getGaussPts().size2();
459 auto &locMat = AssemblyBoundaryEleOp::locMat;
461 auto t_normal = getFTensor1Normal();
462 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
464 auto t_disp = getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactDisp);
466 getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactTraction);
467 auto t_coords = getFTensor1CoordsAtGaussPts();
469 auto t_w = getFTensor0IntegrationWeight();
470 auto t_row_base = row_data.getFTensor1N<3>();
471 size_t nb_face_functions = row_data.getN().size2() / 3;
473 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
475 const double alpha = t_w * getMeasure();
478 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
483 getTStime(), t_spatial_coords, t_traction);
484 auto tn = -t_traction(
i) * t_grad_sdf(
i);
488 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
490 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
493 t_res_dt(
i,
j) = -
cn * t_cQ(
i,
j);
496 for (; rr != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++rr) {
498 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
500 const double row_base = t_row_base(
i) * t_normal(
i);
502 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
503 for (
size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols /
SPACE_DIM;
505 const double col_base = t_col_base(
i) * t_normal(
i);
506 const double beta = alpha * row_base * col_base;
508 t_mat(
i,
j) -= beta * t_res_dt(
i,
j);
516 for (; rr < nb_face_functions; ++rr)
#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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
def grad_sdf(t, x, y, z, tx, ty, tz)
def hess_sdf(t, x, y, z, tx, ty, tz)
constexpr double t
plate stiffness
constexpr auto field_name