8#ifndef __CONTACTOPS_HPP__
9#define __CONTACTOPS_HPP__
14struct CommonData :
public boost::enable_shared_from_this<CommonData> {
27 static SmartPetscObj<Vec>
31 constexpr int ghosts[] = {0, 1, 2};
33 createGhostVector(m_field.
get_comm(),
56 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
61 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactDisp);
65 return boost::shared_ptr<VectorDouble>(shared_from_this(), &
sdfVals);
69 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
gradsSdf);
73 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
hessSdf);
77 return boost::shared_ptr<VectorDouble>(shared_from_this(), &
constraintVals);
88 SDFPython() =
default;
89 virtual ~SDFPython() =
default;
91 MoFEMErrorCode sdfInit(
const std::string py_file) {
96 auto main_module = bp::import(
"__main__");
97 mainNamespace = main_module.attr(
"__dict__");
98 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
100 sdfFun = mainNamespace[
"sdf"];
101 sdfGradFun = mainNamespace[
"grad_sdf"];
102 sdfHessFun = mainNamespace[
"hess_sdf"];
104 }
catch (bp::error_already_set
const &) {
112 template <
typename T>
113 inline std::vector<T>
114 py_list_to_std_vector(
const boost::python::object &iterable) {
115 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
116 boost::python::stl_input_iterator<T>());
121 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
122 np::ndarray tx, np::ndarray ty, np::ndarray tz,
int block_id,
130 sdf = bp::extract<np::ndarray>(
131 sdfFun(delta_t,
t, x, y, z, tx, ty, tz, block_id));
133 }
catch (bp::error_already_set
const &) {
143 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
144 np::ndarray tx, np::ndarray ty, np::ndarray tz,
int block_id,
145 np::ndarray &grad_sdf
152 grad_sdf = bp::extract<np::ndarray>(
153 sdfGradFun(delta_t,
t, x, y, z, tx, ty, tz, block_id));
155 }
catch (bp::error_already_set
const &) {
165 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
166 np::ndarray tx, np::ndarray ty, np::ndarray tz,
int block_id,
167 np::ndarray &hess_sdf
174 hess_sdf = bp::extract<np::ndarray>(
175 sdfHessFun(delta_t,
t, x, y, z, tx, ty, tz, block_id));
177 }
catch (bp::error_already_set
const &) {
186 bp::object mainNamespace;
188 bp::object sdfGradFun;
189 bp::object sdfHessFun;
192static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
194inline np::ndarray convert_to_numpy(VectorDouble &data,
int nb_gauss_pts,
196 auto dtype = np::dtype::get_builtin<double>();
197 auto size = bp::make_tuple(nb_gauss_pts);
198 auto stride = bp::make_tuple(3 *
sizeof(
double));
199 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
205 double delta_t,
double t,
int nb_gauss_pts, MatrixDouble &spatial_coords,
206 MatrixDouble &normals_at_pts,
int block_id)>;
209 double delta_t,
double t,
int nb_gauss_pts, MatrixDouble &spatial_coords,
210 MatrixDouble &normals_at_pts,
int block_id)>;
213 double delta_t,
double t,
int nb_gauss_pts, MatrixDouble &spatial_coords,
214 MatrixDouble &normals_at_pts,
int block_id)>;
218 MatrixDouble &m_spatial_coords,
219 MatrixDouble &m_normals_at_pts,
223 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
225 VectorDouble v_spatial_coords = m_spatial_coords.data();
226 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
228 bp::list python_coords;
229 bp::list python_normals;
231 for (
int idx = 0; idx < 3; ++idx) {
232 python_coords.append(
233 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
234 python_normals.append(
235 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
238 np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
239 np::dtype::get_builtin<double>());
241 bp::extract<np::ndarray>(python_coords[0]),
242 bp::extract<np::ndarray>(python_coords[1]),
243 bp::extract<np::ndarray>(python_coords[2]),
244 bp::extract<np::ndarray>(python_normals[0]),
245 bp::extract<np::ndarray>(python_normals[1]),
246 bp::extract<np::ndarray>(python_normals[2]),
248 "Failed python call");
250 double *sdf_val_ptr =
reinterpret_cast<double *
>(np_sdf.get_data());
253 v_sdf.resize(nb_gauss_pts,
false);
255 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg)
256 v_sdf[gg] = *(sdf_val_ptr + gg);
262 v_sdf.resize(nb_gauss_pts,
false);
263 m_spatial_coords.resize(3, nb_gauss_pts);
264 auto t_coords = getFTensor1FromMat<3>(m_spatial_coords);
266 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
267 v_sdf[gg] = t_coords(1) + 0.5;
276 MatrixDouble &m_spatial_coords,
277 MatrixDouble &m_normals_at_pts,
int block_id) {
279 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
281 VectorDouble v_spatial_coords = m_spatial_coords.data();
282 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
284 bp::list python_coords;
285 bp::list python_normals;
287 for (
int idx = 0; idx < 3; ++idx) {
288 python_coords.append(
289 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
290 python_normals.append(
291 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
294 np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
295 np::dtype::get_builtin<double>());
297 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
298 bp::extract<np::ndarray>(python_coords[1]),
299 bp::extract<np::ndarray>(python_coords[2]),
300 bp::extract<np::ndarray>(python_normals[0]),
301 bp::extract<np::ndarray>(python_normals[1]),
302 bp::extract<np::ndarray>(python_normals[2]), block_id,
304 "Failed python call");
306 double *grad_ptr =
reinterpret_cast<double *
>(np_grad_sdf.get_data());
308 MatrixDouble m_grad_sdf;
309 m_grad_sdf.resize(3, nb_gauss_pts,
false);
310 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
311 for (
int idx = 0; idx < 3; ++idx)
312 m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
317 MatrixDouble m_grad_sdf;
318 m_grad_sdf.resize(3, nb_gauss_pts,
false);
321 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
323 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
324 t_grad_sdf(
i) = t_grad_sdf_set(
i);
333 MatrixDouble &m_spatial_coords,
334 MatrixDouble &m_normals_at_pts,
int block_id) {
336 if (
auto sdf_ptr = sdfPythonWeakPtr.lock()) {
338 VectorDouble v_spatial_coords = m_spatial_coords.data();
339 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
341 bp::list python_coords;
342 bp::list python_normals;
344 for (
int idx = 0; idx < 3; ++idx) {
345 python_coords.append(
346 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
347 python_normals.append(
348 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
351 np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
352 np::dtype::get_builtin<double>());
354 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
355 bp::extract<np::ndarray>(python_coords[1]),
356 bp::extract<np::ndarray>(python_coords[2]),
357 bp::extract<np::ndarray>(python_normals[0]),
358 bp::extract<np::ndarray>(python_normals[1]),
359 bp::extract<np::ndarray>(python_normals[2]), block_id,
361 "Failed python call");
363 double *hess_ptr =
reinterpret_cast<double *
>(np_hess_sdf.get_data());
365 MatrixDouble m_hess_sdf;
366 m_hess_sdf.resize(6, nb_gauss_pts,
false);
367 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
368 for (
int idx = 0; idx < 6; ++idx)
369 m_hess_sdf(idx, gg) =
370 *(hess_ptr + (6 * gg + idx));
375 MatrixDouble m_hess_sdf;
376 m_hess_sdf.resize(6, nb_gauss_pts,
false);
380 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
382 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
383 t_hess_sdf(
i,
j) = t_hess_sdf_set(
i,
j);
391template <
int DIM, IntegrationType I,
typename BoundaryEleOp>
394template <
int DIM, IntegrationType I,
typename BoundaryEleOp>
397template <
int DIM, IntegrationType I,
typename AssemblyBoundaryEleOp>
400template <
int DIM, IntegrationType I,
typename AssemblyBoundaryEleOp>
403template <
int DIM, IntegrationType I,
typename AssemblyBoundaryEleOp>
406template <
typename T1,
typename T2,
int DIM1,
int DIM2>
409 size_t nb_gauss_pts) {
410 MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
411 m_spatial_coords.clear();
412 auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
414 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
415 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
420 return m_spatial_coords;
423template <
typename T1,
int DIM1>
425 size_t nb_gauss_pts) {
426 MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
427 m_normals_at_pts.clear();
429 auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
430 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
431 t_set_normal(
i) = t_normal_at_pts(
i) / t_normal_at_pts.l2();
435 return m_normals_at_pts;
438template <
int DIM,
typename BoundaryEleOp>
442 boost::shared_ptr<CommonData> common_data_ptr,
double scale = 1,
452template <
int DIM,
typename BoundaryEleOp>
467template <
int DIM,
typename AssemblyBoundaryEleOp>
471 boost::shared_ptr<CommonData> common_data_ptr,
473 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
484template <
int DIM,
typename AssemblyBoundaryEleOp>
488 const std::string col_field_name,
489 boost::shared_ptr<CommonData> common_data_ptr,
491 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
492 EntitiesFieldData::EntData &col_data);
504template <
int DIM,
typename AssemblyBoundaryEleOp>
508 const std::string row_field_name,
const std::string col_field_name,
509 boost::shared_ptr<CommonData> common_data_ptr,
511 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
512 EntitiesFieldData::EntData &col_data);
525 template <
int DIM, IntegrationType I>
529 template <
int DIM, IntegrationType I>
537 template <
int DIM, IntegrationType I>
541 template <
int DIM, IntegrationType I>
545 template <
int DIM, IntegrationType I>
552 constexpr auto eps = std::numeric_limits<float>::epsilon();
553 if (std::abs(x) <
eps)
561inline double w(
const double sdf,
const double tn) {
579template <
int DIM,
typename BoundaryEleOp>
582 boost::shared_ptr<CommonData> common_data_ptr,
double scale,
585 commonDataPtr(common_data_ptr), scaleTraction(
scale),
588template <
int DIM,
typename BoundaryEleOp>
598 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
602 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
603 double jacobian = 1.;
604 if (isAxisymmetric) {
605 jacobian = 2. * M_PI * t_coords(0);
608 t_sum_t(
i) += alpha * t_traction(
i);
614 t_sum_t(
i) *= scaleTraction;
616 constexpr int ind[] = {0, 1, 2};
617 CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
623template <
int DIM,
typename BoundaryEleOp>
625 boost::shared_ptr<CommonData> common_data_ptr)
627 commonDataPtr(common_data_ptr) {}
629template <
int DIM,
typename BoundaryEleOp>
636 auto &sdf_vec = commonDataPtr->sdfVals;
637 auto &grad_mat = commonDataPtr->gradsSdf;
638 auto &hess_mat = commonDataPtr->hessSdf;
639 auto &constraint_vec = commonDataPtr->constraintVals;
640 auto &contactTraction_mat = commonDataPtr->contactTraction;
642 sdf_vec.resize(nb_gauss_pts,
false);
643 grad_mat.resize(DIM, nb_gauss_pts,
false);
644 hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
645 constraint_vec.resize(nb_gauss_pts,
false);
647 auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
649 auto t_sdf = getFTensor0FromVec(sdf_vec);
650 auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
651 auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
652 auto t_constraint = getFTensor0FromVec(constraint_vec);
654 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
666 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
674 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
675 m_spatial_coords, m_normals_at_pts, block_id);
678 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
679 m_spatial_coords, m_normals_at_pts, block_id);
682 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
683 m_spatial_coords, m_normals_at_pts, block_id);
685 auto t_sdf_v = getFTensor0FromVec(v_sdf);
686 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
687 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
701 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
703 auto tn = -t_traction(
i) * t_grad_sdf_v(
i);
707 t_grad_sdf(
i) = t_grad_sdf_v(
i);
708 t_hess_sdf(
i,
j) = t_hess_sdf_v(
i,
j);
717template <
int DIM,
typename AssemblyBoundaryEleOp>
720 boost::shared_ptr<CommonData> common_data_ptr,
726template <
int DIM,
typename AssemblyBoundaryEleOp>
729 EntitiesFieldData::EntData &data) {
737 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
741 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
743 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
744 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
745 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
746 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
748 size_t nb_base_functions = data.getN().size2() / 3;
749 auto t_base = data.getFTensor1N<3>();
753 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
757 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
759 auto ts_time = AssemblyBoundaryEleOp::getTStime();
760 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
766 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
767 m_spatial_coords, m_normals_at_pts, block_id);
770 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
771 m_spatial_coords, m_normals_at_pts, block_id);
773 auto t_sdf = getFTensor0FromVec(v_sdf);
774 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
776 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
778 auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
780 double jacobian = 1.;
781 if (isAxisymmetric) {
782 jacobian = 2. * M_PI * t_coords(0);
784 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
786 auto tn = -t_traction(
i) * t_grad_sdf(
i);
790 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
792 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
801 t_cP(
i,
j) * t_disp(
j) +
802 c * (t_sdf * t_grad_sdf(
i));
806 const double beta = alpha * (t_base(
i) * t_normal(
i));
807 t_nf(
i) -= beta * t_rhs(
i);
812 for (; bb < nb_base_functions; ++bb)
827template <
int DIM,
typename AssemblyBoundaryEleOp>
830 const std::string col_field_name,
831 boost::shared_ptr<CommonData> common_data_ptr,
836 AssemblyBoundaryEleOp::sYmm =
false;
839template <
int DIM,
typename AssemblyBoundaryEleOp>
842 EntitiesFieldData::EntData &row_data,
843 EntitiesFieldData::EntData &col_data) {
850 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
853 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
854 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
855 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
857 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
858 auto t_row_base = row_data.getFTensor1N<3>();
859 size_t nb_face_functions = row_data.getN().size2() / 3;
865 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
869 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
871 auto ts_time = AssemblyBoundaryEleOp::getTStime();
872 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
878 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
879 m_spatial_coords, m_normals_at_pts, block_id);
882 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
883 m_spatial_coords, m_normals_at_pts, block_id);
886 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
887 m_spatial_coords, m_normals_at_pts, block_id);
889 auto t_sdf = getFTensor0FromVec(v_sdf);
890 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
891 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
893 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
895 double jacobian = 1.;
896 if (isAxisymmetric) {
897 jacobian = 2. * M_PI * t_coords(0);
899 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
901 auto tn = -t_traction(
i) * t_grad_sdf(
i);
905 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
907 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
910 t_res_dU(
i,
j) = kronecker_delta(
i,
j) + t_cP(
i,
j);
915 (t_hess_sdf(
i,
j) * (t_grad_sdf(
k) * t_traction(
k)) +
916 t_grad_sdf(
i) * t_hess_sdf(
k,
j) * t_traction(
k)) +
917 c * t_sdf * t_hess_sdf(
i,
j);
923 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
925 const double row_base = t_row_base(
i) * t_normal(
i);
927 auto t_col_base = col_data.getFTensor0N(gg, 0);
929 const double beta = alpha * row_base * t_col_base;
931 t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
939 for (; rr < nb_face_functions; ++rr)
954template <
int DIM,
typename AssemblyBoundaryEleOp>
957 const std::string row_field_name,
const std::string col_field_name,
962 AssemblyBoundaryEleOp::sYmm =
false;
965template <
int DIM,
typename AssemblyBoundaryEleOp>
968 iNtegrate(EntitiesFieldData::EntData &row_data,
969 EntitiesFieldData::EntData &col_data) {
976 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
979 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
980 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
981 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
983 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
984 auto t_row_base = row_data.getFTensor1N<3>();
985 size_t nb_face_functions = row_data.getN().size2() / 3;
989 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
993 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
995 auto ts_time = AssemblyBoundaryEleOp::getTStime();
996 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
1002 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1003 m_spatial_coords, m_normals_at_pts, block_id);
1006 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1007 m_spatial_coords, m_normals_at_pts, block_id);
1009 auto t_sdf = getFTensor0FromVec(v_sdf);
1010 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1012 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1014 double jacobian = 1.;
1015 if (isAxisymmetric) {
1016 jacobian = 2. * M_PI * t_coords(0);
1018 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1020 auto tn = -t_traction(
i) * t_grad_sdf(
i);
1024 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
1026 t_cQ(
i,
j) = kronecker_delta(
i,
j) - t_cP(
i,
j);
1034 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1035 const double row_base = t_row_base(
i) * t_normal(
i);
1037 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1039 const double col_base = t_col_base(
i) * t_normal(
i);
1040 const double beta = alpha * row_base * col_base;
1042 t_mat(
i,
j) -= beta * t_res_dt(
i,
j);
1050 for (; rr < nb_face_functions; ++rr)
1064template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
1066 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1070 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
1072 using OpMixDivURhs =
typename B::template OpMixDivTimesU<3, DIM, DIM>;
1073 using OpMixDivUCylRhs =
1074 typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1076 using OpMixLambdaGradURhs =
typename B::template OpMixTensorTimesGradU<DIM>;
1078 typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1082 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1083 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1084 auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1085 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1090 return 2. * M_PI * r;
1095 pip.push_back(
new OpCalculateVectorFieldValues<DIM>(
1096 u, common_data_ptr->contactDispPtr()));
1098 new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1102 new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1104 pip.push_back(
new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1105 sigma, div_stress_ptr));
1112 new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1114 pip.push_back(
new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1118 pip.push_back(
new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1126 using OpMixLhs::OpMixLhs;
1129 EntitiesFieldData::EntData &row_data,
1130 EntitiesFieldData::EntData &col_data) {
1132 auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1133 auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1135 if (side_fe_entity == side_fe_data) {
1136 CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1143template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEle>
1146 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1147 std::string fe_domain_name, std::string sigma, std::string u,
1148 std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1152 using DomainEleOp =
typename DomainEle::UserDataOperator;
1154 auto op_loop_side =
new OpLoopSide<DomainEle>(
1155 m_field, fe_domain_name, DIM, Sev::noisy,
1156 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1157 pip.push_back(op_loop_side);
1159 CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1162 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
1165 using OpMixDivULhs =
typename B::template OpMixDivTimesVec<DIM>;
1166 using OpMixDivUCylLhs =
1167 typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1168 using OpLambdaGraULhs =
typename B::template OpMixTensorTimesGrad<DIM>;
1174 auto unity = []() {
return 1; };
1178 return 2. * M_PI * r;
1184 op_loop_side->getOpPtrVector().push_back(
1185 new OpMixDivULhsSide(sigma, u, unity, jacobian,
true));
1187 op_loop_side->getOpPtrVector().push_back(
1188 new OpMixDivUCylLhsSide(sigma, u, unity, jacobian,
true));
1190 op_loop_side->getOpPtrVector().push_back(
1191 new OpLambdaGraULhsSide(sigma, u, unity, jacobian,
true));
1193 op_loop_side->getSideFEPtr()->getRuleHook = rule;
1197template <
int DIM, AssemblyType A, IntegrationType I,
typename BoundaryEleOp>
1199 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1205 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1207 pip.push_back(
new OpCalculateVectorFieldValues<DIM>(
1208 u, common_data_ptr->contactDispPtr()));
1209 pip.push_back(
new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1210 sigma, common_data_ptr->contactTractionPtr()));
1212 new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1214 pip.push_back(
new typename C::template Assembly<A>::
1215 template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1221template <
int DIM, AssemblyType A, IntegrationType I,
typename BoundaryEleOp>
1223 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1229 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1231 pip.push_back(
new OpCalculateVectorFieldValues<DIM>(
1232 u, common_data_ptr->contactDispPtr()));
1233 pip.push_back(
new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1234 sigma, common_data_ptr->contactTractionPtr()));
1236 new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1242template <
int DIM, IntegrationType I,
typename BoundaryEleOp>
1244 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1250 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1251 pip.push_back(
new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1252 sigma, common_data_ptr->contactTractionPtr()));
1253 pip.push_back(
new typename C::template OpAssembleTotalContactTraction<DIM, I>(
#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
#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 >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
hess_sdf(t, x, y, z, tx, ty, tz)
grad_sdf(t, x, y, z, tx, ty, tz)
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.
default operator for TRI element
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
double getTStimeStep() const
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column