EshelbianPlasticity::OpTreeSearch Struct Reference

#include <users_modules/eshelbian_plasticit/src/EshelbianContact.hpp>

Inheritance diagram for EshelbianPlasticity::OpTreeSearch:

Collaboration diagram for EshelbianPlasticity::OpTreeSearch:

Public Types
using	UOP = FaceElementForcesAndSourcesCore::UserDataOperator

Public Member Functions
	OpTreeSearch (boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > u_h1_ptr, Range r, moab::Interface *post_proc_mesh_ptr, std::vector< EntityHandle > *map_gauss_pts_ptr)
MoFEMErrorCode	doWork (int side, EntityType type, EntitiesFieldData::EntData &data)

Protected Attributes
boost::shared_ptr< ContactTree >	contactTreePtr
boost::shared_ptr< ContactOps::CommonData >	commonDataPtr
boost::shared_ptr< MatrixDouble >	uH1Ptr
moab::Interface *	postProcMeshPtr = nullptr
std::vector< EntityHandle > *	mapGaussPtsPtr = nullptr
Range	contactRange

Detailed Description

Definition at line 209 of file EshelbianContact.hpp.

Member Typedef Documentation

UOP

using EshelbianPlasticity::OpTreeSearch::UOP = FaceElementForcesAndSourcesCore::UserDataOperator

Definition at line 211 of file EshelbianContact.hpp.

Constructor & Destructor Documentation

OpTreeSearch()

EshelbianPlasticity::OpTreeSearch::OpTreeSearch	(	boost::shared_ptr< ContactTree >	contact_tree_ptr,
		boost::shared_ptr< ContactOps::CommonData >	common_data_ptr,
		boost::shared_ptr< MatrixDouble >	u_h1_ptr,
		Range	r,
		moab::Interface *	post_proc_mesh_ptr,
		std::vector< EntityHandle > *	map_gauss_pts_ptr
	)

Definition at line 1058 of file EshelbianContact.cpp.

    : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
      commonDataPtr(common_data_ptr), uH1Ptr(u_h1_ptr),
      postProcMeshPtr(post_proc_mesh_ptr), mapGaussPtsPtr(map_gauss_pts_ptr),
      contactRange(r) {}

Member Function Documentation

doWork()

MoFEMErrorCode EshelbianPlasticity::OpTreeSearch::doWork	(	int	side,
		EntityType	type,
		EntitiesFieldData::EntData &	data
	)

Definition at line 1072 of file EshelbianContact.cpp.

                                                                    {
  MoFEMFunctionBegin;
 
  auto &m_field = getPtrFE()->mField;
  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
  auto fe_id = id_from_handle(fe_ent);
 
  if (contactRange.find(fe_ent) == contactRange.end())
    MoFEMFunctionReturnHot(0);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  const auto nb_gauss_pts = getGaussPts().size2();
 
  auto t_disp_h1 = getFTensor1FromMat<3>(*uH1Ptr);
  auto t_coords = getFTensor1CoordsAtGaussPts();
  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
 
  auto next = [&]() {
    ++t_disp_h1;
    ++t_traction;
    ++t_coords;
  };
 
  auto get_ele_centre = [i](auto t_ele_coords) {
    FTensor::Tensor1<double, 3> t_ele_center;
    t_ele_center(i) = 0;
    for (int nn = 0; nn != 3; nn++) {
      t_ele_center(i) += t_ele_coords(i);
      ++t_ele_coords;
    }
    t_ele_center(i) /= 3;
    return t_ele_center;
  };
 
  auto get_ele_radius = [i](auto t_ele_center, auto t_ele_coords) {
    FTensor::Tensor1<double, 3> t_n0;
    t_n0(i) = t_ele_center(i) - t_ele_coords(i);
    return t_n0.l2();
  };
 
  auto get_face_conn = [this](auto face) {
    const EntityHandle *conn;
    int num_nodes;
    CHK_MOAB_THROW(contactTreePtr->getPostProcMesh().get_connectivity(
                       face, conn, num_nodes, true),
                   "get conn");
    if (num_nodes != 3) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "face is not a triangle");
    }
    return conn;
  };
 
  auto get_face_coords = [this](auto conn) {
    std::array<double, 9> coords;
    CHKERR contactTreePtr->getPostProcMesh().get_coords(conn, 3, coords.data());
    return coords;
  };
 
  auto get_closet_face = [this](auto *point_ptr, auto r) {
    FTensor::Tensor1<double, 3> t_point_out;
    std::vector<EntityHandle> faces_out;
    CHK_MOAB_THROW(
        contactTreePtr->getTreeSurfPtr()->sphere_intersect_triangles(
            point_ptr, r / 3, contactTreePtr->getRootSetSurf(), faces_out),
        "get closest faces");
    return faces_out;
  };
 
  auto get_faces_out = [this](auto *point_ptr, auto *unit_ray_ptr, auto radius,
                              auto eps) {
    std::vector<double> distances_out;
    std::vector<EntityHandle> faces_out;
    CHK_MOAB_THROW(
 
        contactTreePtr->getTreeSurfPtr()->ray_intersect_triangles(
            distances_out, faces_out, contactTreePtr->getRootSetSurf(), eps,
            point_ptr, unit_ray_ptr, &radius),
 
        "get closest faces");
    return std::make_pair(faces_out, distances_out);
  };
 
  auto get_normal = [](auto &ele_coords) {
    FTensor::Tensor1<double, 3> t_normal;
    Tools::getTriNormal(ele_coords.data(), &t_normal(0));
    return t_normal;
  };
 
  auto make_map = [&](auto &face_out, auto &face_dist, auto &t_ray_point,
                      auto &t_unit_ray, auto &t_master_coord) {
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    std::map<double, EntityHandle> m;
    for (auto ii = 0; ii != face_out.size(); ++ii) {
      auto face_conn = get_face_conn(face_out[ii]);
      auto face_coords = get_face_coords(face_conn);
      auto t_face_normal = get_normal(face_coords);
      t_face_normal.normalize();
      FTensor::Tensor1<double, 3> t_x;
      t_x(i) = t_ray_point(i) + t_unit_ray(i) * face_dist[ii];
      FTensor::Tensor2<double, 3, 3> t_chi;
      t_chi(i, j) =
          t_x(i) * t_face_normal(j) - t_master_coord(i) * t_unit_ray(j);
      if (t_unit_ray(i) * t_face_normal(i) > std::cos(M_PI / 3)) {
        auto dot = std::sqrt(t_chi(i, j) * t_chi(i, j));
        m[dot] = face_out[ii];
      }
    }
    return m;
  };
 
  auto get_tag_data = [this](auto tag, auto face, auto &vec) {
    MoFEMFunctionBegin;
    int tag_length;
    CHKERR contactTreePtr->getPostProcMesh().tag_get_length(tag, tag_length);
    vec.resize(tag_length);
    CHKERR contactTreePtr->getPostProcMesh().tag_get_data(tag, &face, 1,
                                                          &*vec.begin());
    MoFEMFunctionReturn(0);
  };
 
  auto create_tag = [this](const std::string tag_name, const int size) {
    double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
    Tag th;
    if (postProcMeshPtr) {
      CHKERR postProcMeshPtr->tag_get_handle(
          tag_name.c_str(), size, MB_TYPE_DOUBLE, th,
          MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
      CHKERR postProcMeshPtr->tag_clear_data(th, &*mapGaussPtsPtr->begin(),
                                             mapGaussPtsPtr->size(), def_VAL);
    }
    return th;
  };
 
  auto set_float_precision = [](const double x) {
    if (std::abs(x) < std::numeric_limits<float>::epsilon())
      return 0.;
    else
      return x;
  };
 
  // scalars
  auto save_scal_tag = [&](auto &th, auto v, const int gg) {
    MoFEMFunctionBegin;
    if (postProcMeshPtr) {
      v = set_float_precision(v);
      CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1, &v);
    }
    MoFEMFunctionReturn(0);
  };
 
  // adjacencies
  auto get_fe_adjacencies = [this](auto fe_ent) {
    Range adj_faces;
    CHK_MOAB_THROW(getFaceFE()->mField.get_moab().get_adjacencies(
                       &fe_ent, 1, 2, false, adj_faces, moab::Interface::UNION),
                   "get adj");
    std::set<int> adj_ids;
    for (auto f : adj_faces) {
      adj_ids.insert(id_from_handle(f));
    }
    return adj_ids;
  };
 
  auto get_face_id = [this](auto face) {
    int id;
    if (contactTreePtr->getPostProcMesh().tag_get_data(
            contactTreePtr->thEleId, &face, 1, &id) == MB_SUCCESS) {
      return id;
    }
    return -1;
  };
 
  auto get_body_id = [this](auto face) {
    int id;
    if (contactTreePtr->getPostProcMesh().tag_get_data(
            contactTreePtr->thBodyId, &face, 1, &id) == MB_SUCCESS) {
      return id;
    }
    return -1;
  };
 
  auto get_face_part = [this](auto face) {
    ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
        &(contactTreePtr->getPostProcMesh()), MYPCOMM_INDEX);
    int part;
    if (contactTreePtr->getPostProcMesh().tag_get_data(
            pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
      return part;
    }
    return -1;
  };
 
  auto check_face = [&](auto face, auto fe_id, auto part) {
    auto face_id = get_face_id(face);
    auto face_part = get_face_part(face);
    if (face_id == fe_id && face_part == part)
      return true;
    return false;
  };
 
  // vectors
  VectorDouble3 v(3);
  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
    MoFEMFunctionBegin;
    if (postProcMeshPtr) {
      t_v(i) = t_d(i);
      for (auto &a : v.data())
        a = set_float_precision(a);
      CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1,
                                           &*v.data().begin());
    }
    MoFEMFunctionReturn(0);
  };
 
  Tag th_mark = create_tag("contact_mark", 1);
  Tag th_mark_slave = create_tag("contact_mark_slave", 1);
  Tag th_body_id = create_tag("contact_body_id", 1);
  Tag th_gap = create_tag("contact_gap", 1);
  Tag th_tn_master = create_tag("contact_tn_master", 1);
  Tag th_tn_slave = create_tag("contact_tn_slave", 1);
  Tag th_contact_traction = create_tag("contact_traction", 3);
  Tag th_contact_traction_master = create_tag("contact_traction_master", 3);
  Tag th_contact_traction_slave = create_tag("contact_traction_slave", 3);
  Tag th_c = create_tag("contact_c", 1);
  Tag th_normal = create_tag("contact_normal", 3);
  Tag th_dist = create_tag("contact_dip", 3);
 
  auto t_ele_centre = get_ele_centre(getFTensor1Coords());
  auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
 
  contactTreePtr->shadowDataMap.clear();
  auto &shadow_vec = contactTreePtr->shadowDataMap[fe_ent];
  shadow_vec.clear();
 
  auto adj_fe_ids = get_fe_adjacencies(fe_ent);
 
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
    FTensor::Tensor1<double, 3> t_spatial_coords;
    t_spatial_coords(i) = t_coords(i) + t_disp_h1(i);
 
    if (postProcMeshPtr) {
      CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
    }
 
    auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
    for (auto face_close : faces_close) {
      if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
 
        auto body_id = get_body_id(face_close);
 
        auto master_face_conn = get_face_conn(face_close);
        std::array<double, 9> master_coords;
        CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
            contactTreePtr->thSmallX, master_face_conn, 3,
            master_coords.data());
        std::array<double, 9> master_traction;
        CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
            contactTreePtr->thTraction, master_face_conn, 3,
            master_traction.data());
        auto t_normal_face_close = get_normal(master_coords);
        t_normal_face_close.normalize();
 
        if (postProcMeshPtr) {
          double m = 1;
          CHKERR save_scal_tag(th_mark, m, gg);
          CHKERR save_scal_tag(th_body_id, static_cast<double>(body_id), gg);
          CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
          CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
        }
 
        FTensor::Tensor1<double, 3> t_unit_ray;
        t_unit_ray(i) = -t_normal_face_close(i);
        FTensor::Tensor1<double, 3> t_ray_point;
        t_ray_point(i) = t_spatial_coords(i) - t_unit_ray(i) * ele_radius;
 
        constexpr double eps = 1e-3;
        auto [faces_out, faces_dist] = get_faces_out(
            &t_ray_point(0), &t_unit_ray(0), 2 * ele_radius, eps * ele_radius);
 
        auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
                          t_spatial_coords);
        for (auto m_it = m.begin(); m_it != m.end(); ++m_it) {
          auto face = m_it->second;
          if (face != face_close) {
 
            if (
 
                body_id != get_body_id(face) &&
                (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
                 get_face_part(face) != m_field.get_comm_rank())
 
            ) {
 
              shadow_vec.push_back(ContactTree::FaceData());
              shadow_vec.back().gaussPtNb = gg;
 
              auto slave_face_conn = get_face_conn(face);
              std::array<double, 9> slave_coords;
              CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
                  contactTreePtr->thSmallX, slave_face_conn, 3,
                  slave_coords.data());
              auto t_normal_face = get_normal(slave_coords);
              std::array<double, 9> slave_tractions;
              CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
                  contactTreePtr->thTraction, slave_face_conn, 3,
                  slave_tractions.data());
 
              auto t_master_point =
                  getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
              auto t_slave_point =
                  getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
              auto t_ray_point_data =
                  getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
              auto t_unit_ray_data =
                  getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
 
              t_slave_point(i) = t_ray_point(i) + m_it->first * t_unit_ray(i);
 
              auto eval_position = [&](auto &&t_elem_coords, auto &&t_point) {
                std::array<double, 2> loc_coords;
                CHK_THROW_MESSAGE(
                    Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
                        &t_elem_coords(0, 0), &t_point(0), 1,
                        loc_coords.data()),
                    "get local coords");
                FTensor::Tensor1<double, 3> t_shape_fun;
                CHK_THROW_MESSAGE(Tools::shapeFunMBTRI<0>(&t_shape_fun(0),
                                                          &loc_coords[0],
                                                          &loc_coords[1], 1),
                                  "calc shape fun");
                FTensor::Index<'i', 3> i;
                FTensor::Index<'j', 3> j;
                FTensor::Tensor1<double, 3> t_point_out;
                t_point_out(i) = t_shape_fun(j) * t_elem_coords(j, i);
                return t_point_out;
              };
 
              auto t_master_point_updated = eval_position(
                  getFTensor2FromPtr<3, 3>(master_coords.data()),
                  getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
              t_master_point(i) = t_master_point_updated(i);
 
              auto t_slave_point_updated = eval_position(
                  getFTensor2FromPtr<3, 3>(slave_coords.data()),
                  getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
              t_slave_point(i) = t_slave_point_updated(i);
 
              t_ray_point_data(i) = t_ray_point(i);
              t_unit_ray_data(i) = t_unit_ray(i);
 
              std::copy(master_coords.begin(), master_coords.end(),
                        shadow_vec.back().masterPointNodes.begin());
              std::copy(master_traction.begin(), master_traction.end(),
                        shadow_vec.back().masterTractionNodes.begin());
              std::copy(slave_coords.begin(), slave_coords.end(),
                        shadow_vec.back().slavePointNodes.begin());
              std::copy(slave_tractions.begin(), slave_tractions.end(),
                        shadow_vec.back().slaveTractionNodes.begin());
 
              shadow_vec.back().eleRadius = ele_radius;
 
              // CHKERR get_tag_data(contactTreePtr->thIds, face,
              //                     shadow_vec.back().dofsSlaveIds);
              // CHKERR get_tag_data(contactTreePtr->thCoeff, face,
              //                     shadow_vec.back().dofsSlaveCoeff);
              // CHKERR get_tag_data(contactTreePtr->thBases, face,
              //                     shadow_vec.back().baseSlaveFuncs);
 
              if (postProcMeshPtr) {
                auto [gap, tn_master, tn_slave, c, t_master_traction,
                      t_slave_traction] =
                    multiGetGap(&(shadow_vec.back()), t_spatial_coords);
                FTensor::Tensor1<double, 3> t_gap_vec;
                t_gap_vec(i) = t_slave_point(i) - t_spatial_coords(i);
                CHKERR save_scal_tag(th_gap, gap, gg);
                CHKERR save_scal_tag(th_tn_master, tn_master, gg);
                CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
                CHKERR save_scal_tag(th_c, c, gg);
                double m = 1;
                CHKERR save_scal_tag(th_mark_slave, m, gg);
                CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
                CHKERR save_vec_tag(th_contact_traction_master,
                                    t_master_traction, gg);
                CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,
                                    gg);
              }
 
              break;
            }
          }
        }
        break;
      }
    }
    next();
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

commonDataPtr

boost::shared_ptr<ContactOps::CommonData> EshelbianPlasticity::OpTreeSearch::commonDataPtr

protected

Definition at line 226 of file EshelbianContact.hpp.

contactRange

Range EshelbianPlasticity::OpTreeSearch::contactRange

protected

Definition at line 232 of file EshelbianContact.hpp.

contactTreePtr

boost::shared_ptr<ContactTree> EshelbianPlasticity::OpTreeSearch::contactTreePtr

protected

Definition at line 225 of file EshelbianContact.hpp.

mapGaussPtsPtr

std::vector<EntityHandle>* EshelbianPlasticity::OpTreeSearch::mapGaussPtsPtr = nullptr

protected

Definition at line 230 of file EshelbianContact.hpp.

postProcMeshPtr

moab::Interface* EshelbianPlasticity::OpTreeSearch::postProcMeshPtr = nullptr

protected

Definition at line 229 of file EshelbianContact.hpp.

uH1Ptr

boost::shared_ptr<MatrixDouble> EshelbianPlasticity::OpTreeSearch::uH1Ptr

protected

Definition at line 227 of file EshelbianContact.hpp.

---

The documentation for this struct was generated from the following files:

• users_modules/eshelbian_plasticit/src/EshelbianContact.hpp
• users_modules/eshelbian_plasticit/src/impl/EshelbianContact.cpp