34 boost::shared_ptr<CommonData> common_data_ptr,
35 moab::Interface *moab_inst);
45 boost::shared_ptr<CommonData> common_data_ptr);
54 const std::string col_field_name,
55 boost::shared_ptr<CommonData> common_data_ptr);
65 const std::string row_field_name,
const std::string col_field_name,
66 boost::shared_ptr<CommonData> common_data_ptr);
77 boost::shared_ptr<CommonData> common_data_ptr,
91 return -
g + t_disp(
i) * t_normal(
i);
94template <
typename T1,
typename T2>
97 return -t_traction(
i) * t_normal(
i);
109inline double w(
const double g,
const double t) {
110 return g - (*cache).cn_cont *
t;
114 return (
w(
g,
t) + std::abs(
w(
g,
t))) / 2 + g0;
118 return -(*cache).cn_cont * (1 +
sign(
w(
g,
t))) / 2;
122 return (
sign(
w(
g,
t)) - 1) / 2;
126 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
128 commonDataPtr(common_data_ptr) {}
133 const size_t nb_gauss_pts = getGaussPts().size2();
134 const size_t nb_dofs = data.
getIndices().size();
138 auto t_normal = getFTensor1Normal();
139 const double ls = sqrt(t_normal(
i) * t_normal(
i));
142 auto t_w = getFTensor0IntegrationWeight();
143 auto t_disp = getFTensor1FromMat<3>(*(
commonDataPtr->mDispPtr));
146 auto t_coords = getFTensor1CoordsAtGaussPts();
148 size_t nb_base_functions = data.
getN().size2() / 3;
151 auto t_contact_normal =
153 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
155 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
157 auto t_nf = getFTensor1FromArray<3, 3>(locF);
159 const double alpha = t_w * getMeasure();
160 if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
161 VectorDouble
n = getNormalsAtGaussPts(gg);
162 auto t_n = getFTensor1FromPtr<3>(&*
n.data().begin());
163 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
165 Tensor2<double, 3, 3> t_P;
166 t_P(
i,
j) = t_contact_normal(
i) * t_contact_normal(
j);
168 Tensor2<double, 3, 3> t_Q;
169 t_Q(
i,
j) = kronecker_delta(
i,
j) - t_P(
i,
j);
171 const double gap_0 =
gap0(t_gap, t_disp, t_contact_normal);
173 t_rhs_constrains(
i) =
174 t_contact_normal(
i) *
179 t_rhs_tangent_disp(
i) = t_Q(
i,
j) * t_disp(
j);
180 t_rhs_tangent_traction(
i) = (*cache).cn_cont * t_Q(
i,
j) * t_traction(
j);
183 for (; bb != nb_dofs / 3; ++bb) {
184 const double beta = alpha * (t_base(
i) * t_normal(
i));
186 t_nf(
i) -= beta * t_rhs_constrains(
i);
187 t_nf(
i) -= beta * t_rhs_tangent_disp(
i);
188 t_nf(
i) += beta * t_rhs_tangent_traction(
i);
193 for (; bb < nb_base_functions; ++bb)
210 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
211 moab::Interface *moab_inst)
213 commonDataPtr(common_data_ptr), moab(moab_inst) {
223 const size_t nb_dofs = data.
getIndices().size();
226 [](
int i) { return i == -1; });
231 auto get_tag = [&](
const std::string name,
size_t size) {
233 double def[3] = {0, 0, 0};
234 CHKERR moab->tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
235 MB_TAG_CREAT | MB_TAG_SPARSE, def);
239 auto th_gap = get_tag(
"GAP", 1);
240 auto th_cons = get_tag(
"CONSTRAINT", 1);
241 auto th_roller = get_tag(
"ROLLER_ID", 1);
242 auto th_traction = get_tag(
"TRACTION", 3);
243 auto th_normal = get_tag(
"NORMAL", 3);
244 auto th_normal_tri = get_tag(
"NORMAL_TRI", 3);
245 auto th_disp = get_tag(
"DISPLACEMENT", 3);
247 auto t_contact_normal =
249 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
250 auto t_traction = getFTensor1FromMat<3>(*(
commonDataPtr->contactTractionPtr));
251 auto t_disp = getFTensor1FromMat<3>(*(
commonDataPtr->mDispPtr));
254 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
256 double coords[] = {0, 0, 0};
258 for (
int dd = 0; dd != 3; dd++) {
261 CHKERR moab->create_vertex(&coords[0], new_vertex);
266 t_contact_normal(2));
267 trac(
i) /= (*cache).young_modulus_inv;
268 double fin_gap = t_gap;
269 const double gap_0 =
gap0(t_gap, t_disp, t_contact_normal);
272 CHKERR moab->tag_set_data(th_gap, &new_vertex, 1, &fin_gap);
273 CHKERR moab->tag_set_data(th_cons, &new_vertex, 1, &constra);
275 CHKERR moab->tag_set_data(th_traction, &new_vertex, 1, &trac(0));
276 CHKERR moab->tag_set_data(th_normal, &new_vertex, 1, &norm(0));
277 CHKERR moab->tag_set_data(th_disp, &new_vertex, 1, &disp(0));
283 CHKERR moab->tag_set_data(th_normal_tri, &new_vertex, 1, &el_normal(0));
285 double d_tag_data = (*cache).closest_roller_id;
286 CHKERR moab->tag_set_data(th_roller, &new_vertex, 1, &d_tag_data);
299 const std::string row_field_name,
const std::string col_field_name,
300 boost::shared_ptr<CommonData> common_data_ptr)
302 commonDataPtr(common_data_ptr) {
309 const size_t nb_gauss_pts = getGaussPts().size2();
310 const size_t row_nb_dofs = row_data.
getIndices().size();
311 const size_t col_nb_dofs = col_data.
getIndices().size();
313 if (row_nb_dofs && col_nb_dofs) {
316 auto t_normal = getFTensor1Normal();
317 const double ls = sqrt(t_normal(
i) * t_normal(
i));
320 auto t_disp = getFTensor1FromMat<3>(*(
commonDataPtr->mDispPtr));
323 auto t_coords = getFTensor1CoordsAtGaussPts();
325 auto t_w = getFTensor0IntegrationWeight();
327 size_t nb_face_functions = row_data.
getN().size2() / 3;
329 auto t_contact_normal =
331 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
335 getFTensor2FromMat<3, 3>(*(
commonDataPtr->contactNormalDiffPtr));
337 auto &
cn = (*cache).cn_cont;
339 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
341 const double alpha = t_w * getMeasure();
343 if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
344 VectorDouble
n = getNormalsAtGaussPts(gg);
345 auto t_n = getFTensor1FromPtr<3>(&*
n.data().begin());
346 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
349 Tensor3<double, 3, 3, 3> t_diff_contact_tangent_tensor;
350 t_diff_contact_tangent_tensor(
i,
j,
k) =
351 -t_diff_normal(
i,
k) * t_contact_normal(
j) -
352 t_contact_normal(
i) * t_diff_normal(
j,
k);
354 Tensor2<double, 3, 3> t_P;
355 t_P(
i,
j) = t_contact_normal(
i) * t_contact_normal(
j);
357 Tensor2<double, 3, 3> t_Q;
358 t_Q(
i,
j) = kronecker_delta(
i,
j) - t_P(
i,
j);
359 const double gap_0 =
gap0(t_gap, t_disp, t_contact_normal);
361 t1(
i) = t_contact_normal(
i) - 0.5 * t_gap_diff(
i) +
362 0.5 *
cn * t_traction(
j) * t_diff_normal(
i,
j) +
363 t_disp(
j) * t_diff_normal(
i,
j) +
364 0.5 *
sign(t_gap +
cn * t_traction(
k) * t_contact_normal(
k)) *
365 (t_gap_diff(
i) +
cn * t_traction(
j) * t_diff_normal(
i,
j));
370 for (; rr != row_nb_dofs / 3; ++rr) {
372 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
373 const double row_base = t_row_base(
i) * t_normal(
i);
377 for (
size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
378 const double beta = alpha * row_base * t_col_base;
381 t_mat(
i,
j) -= beta * t_contact_normal(
i) * t1(
j);
382 t_mat(
i,
j) -= (beta * constrain) * t_diff_normal(
i,
j);
385 t_mat(
i,
j) -= beta * t_Q(
i,
j);
387 beta * t_diff_contact_tangent_tensor(
i,
j,
k) * t_disp(
j);
390 t_mat(
i,
k) += beta * t_diff_contact_tangent_tensor(
i,
j,
k) *
391 (*cache).cn_cont * t_traction(
j);
398 for (; rr < nb_face_functions; ++rr)
417 const std::string row_field_name,
const std::string col_field_name,
418 boost::shared_ptr<CommonData> common_data_ptr)
420 commonDataPtr(common_data_ptr) {
427 const size_t nb_gauss_pts = getGaussPts().size2();
428 const size_t row_nb_dofs = row_data.
getIndices().size();
429 const size_t col_nb_dofs = col_data.
getIndices().size();
431 if (row_nb_dofs && col_nb_dofs) {
433 auto t_normal = getFTensor1Normal();
434 const double ls = sqrt(t_normal(
i) * t_normal(
i));
437 auto t_disp = getFTensor1FromMat<3>(*(
commonDataPtr->mDispPtr));
440 auto t_coords = getFTensor1CoordsAtGaussPts();
442 auto t_w = getFTensor0IntegrationWeight();
444 size_t nb_face_functions = row_data.
getN().size2() / 3;
446 auto t_contact_normal =
448 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
450 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
452 const double alpha = t_w * getMeasure();
454 if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
455 VectorDouble
n = getNormalsAtGaussPts(gg);
456 auto t_n = getFTensor1FromPtr<3>(&*
n.data().begin());
457 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
461 Tensor2<double, 3, 3> t_P;
462 t_P(
i,
j) = t_contact_normal(
i) * t_contact_normal(
j);
464 Tensor2<double, 3, 3> t_Q;
465 t_Q(
i,
j) = kronecker_delta(
i,
j) - t_P(
i,
j);
471 for (; rr != row_nb_dofs / 3; ++rr) {
472 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
473 const double row_base = t_row_base(
i) * t_normal(
i);
477 for (
size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
478 const double col_base = t_col_base(
i) * t_normal(
i);
479 const double beta = alpha * row_base * col_base;
481 t_mat(
i,
j) += (beta * diff_traction) * t_P(
i,
j);
482 t_mat(
i,
j) += beta * (*cache).cn_cont * t_Q(
i,
j);
490 for (; rr < nb_face_functions; ++rr)
508 boost::shared_ptr<CommonData> common_data_ptr,
bool update)
511 commonDataPtr(common_data_ptr), uPdate(update) {
521 const size_t nb_dofs = data.
getIndices().size();
525 auto t_disp = getFTensor1FromMat<3>(*(
commonDataPtr->mDispPtr));
527 commonDataPtr->contactNormalPtr->resize(3, nb_gauss_pts,
false);
529 commonDataPtr->contactNormalDiffPtr->resize(9, nb_gauss_pts,
false);
531 commonDataPtr->contactGapDiffPtr->resize(3, nb_gauss_pts,
false);
536 auto t_contact_normal =
539 getFTensor2FromMat<3, 3>(*(
commonDataPtr->contactNormalDiffPtr));
540 auto t_gap_diff = getFTensor1FromMat<3>(*(
commonDataPtr->contactGapDiffPtr));
541 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
545 int roll_id = (*cache).getClosestRollerID(t_coords, t_disp);
552 (*cache).closest_roller = &(*cache).rollerDataVec.at(tag_data);
553 (*cache).closest_roller_id = tag_data;
555 const bool is_lhs = (int)
getTSCtx() == 3;
556 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
558 auto t_contact = (*cache).closest_roller->getNormal(t_coords, t_disp);
559 t_contact_normal(
i) = t_contact(
i);
560 t_gap = (*cache).closest_roller->getGap(t_coords);
563 auto t_dn_du = ((*cache).closest_roller)
564 ->getDiffNormal(t_coords, t_disp, t_contact_normal);
565 t_diff_normal(
i,
j) = t_dn_du(
i,
j);
567 (*cache).closest_roller->getdGap(t_coords, t_contact_normal);
568 t_gap_diff(
i) = t_d_gap(
i);
586 moab::Interface &post_proc_mesh,
587 std::vector<EntityHandle> &map_gauss_pts,
588 boost::shared_ptr<CommonData> common_data_ptr);
599 const std::string
field_name, moab::Interface &post_proc_mesh,
600 std::vector<EntityHandle> &map_gauss_pts,
601 boost::shared_ptr<CommonData> common_data_ptr)
603 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
604 commonDataPtr(common_data_ptr) {
615 auto get_tag_mat = [&](
const std::string name) {
616 std::array<double, 9> def;
617 std::fill(def.begin(), def.end(), 0);
620 MB_TAG_CREAT | MB_TAG_SPARSE,
625 auto get_tag_vec = [&](
const std::string name) {
626 std::array<double, 3> def;
627 std::fill(def.begin(), def.end(), 0);
630 MB_TAG_CREAT | MB_TAG_SPARSE,
635 MatrixDouble3by3 mat(3, 3);
637 auto set_matrix = [&](
auto &
t) -> MatrixDouble3by3 & {
639 for (
size_t r = 0; r != 3; ++r)
640 for (
size_t c = 0;
c != 3; ++
c)
645 auto set_vector = [&](
auto &
t) -> MatrixDouble3by3 & {
647 for (
size_t r = 0; r != 3; ++r)
652 auto set_tag = [&](
auto th,
auto gg, MatrixDouble3by3 &mat) {
654 &*mat.data().begin());
657 auto th_stress = get_tag_mat(
"SIGMA");
658 auto th_div = get_tag_vec(
"DIV_SIGMA");
661 auto t_stress = getFTensor2FromMat<3, 3>(*(
commonDataPtr->contactStressPtr));
663 getFTensor1FromMat<3>(*(
commonDataPtr->contactStressDivergencePtr));
665 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
667 t_stress(
i,
j) /= (*cache).young_modulus_inv;
668 t_div(
i) /= (*cache).young_modulus_inv;
672 CHKERR set_tag(th_div, gg, set_vector(t_div));
698 boost::shared_ptr<CommonData> common_data_ptr)
707 const size_t nb_dofs = data.
getIndices().size();
711 [](
int i) { return i == -1; });
719 auto t_contact_normal =
721 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
722 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
726 if (t_gap <= (*cache).cn_cont * trac)
745 boost::shared_ptr<CommonData> common_data_ptr)
754 const size_t nb_dofs = data.
getIndices().size();
759 [](
int i) { return i == -1; });
769 auto t_contact_normal =
771 auto t_gap = getFTensor0FromVec(*(
commonDataPtr->contactGapPtr));
778 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
781 const auto sign_w =
sign(
w(t_gap, trac));
783 for (
int dd = 0; dd != 3; ++dd)
784 l_reactions[3 *
id + dd] +=
785 t_traction(dd) * alpha / (*cache).young_modulus_inv;
788 data_vec[0] += alpha;
789 data_vec[1] += alpha;
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr double t
plate stiffness
constexpr auto field_name
virtual moab::Interface & get_moab()=0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
double getMeasure()
get measure of element
VectorDouble & getNormal()
get triangle normal
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
const TSMethod::TSContext getTSCtx() const
ForcesAndSourcesCore * getPtrFE() const
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEM::Interface & mField
OpGetClosestRigidBodyData(MoFEM::Interface &m_field, const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, bool update)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
moab::Interface & postProcMesh
OpPostProcContact(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
[Operators definitions]
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
OpPostProcVertex(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, moab::Interface *moab_inst)