Contact.cc 2.81 KB
Newer Older
1 2 3 4 5
/* 
    This file is a part of eXlibris C++ Library
    under the GNU General Public License:
    See the LICENSE.md files for terms and 
    conditions.
6
*/
7 8


9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
#include "Contact.h"
// Xfem
#include "xSimpleGeometry.h"

using namespace xfem;

OrientedSurface Sphere::operator()(const Trellis_Util::mPoint& p) const
{
  xVector normal(origin,p);
  normal.norm();
  return std::make_pair(Trellis_Util::mPoint(origin(0)+radius*normal(0),origin(1)+radius*normal(1),origin(2)+radius*normal(2)), normal);
}

OrientedSurface Plane::operator()(const Trellis_Util::mPoint& p) const
{
  const double distance = xVector(point_on_plane, p)*normal;
  return std::make_pair(Trellis_Util::mPoint(p(0)-distance*normal(0),p(1)-distance*normal(1),p(2)-distance*normal(2)), normal);
}

OrientedSurface InterVector::operator()(const Trellis_Util::mPoint& p)  const
{
  std::vector<double> ds_val;
  std::vector<PointToOrientedSurface>::const_iterator it;
  double max_val = -std::numeric_limits<double>::max();
  int max_i = -1;
  int i=0;
  for (it=ds.begin(); it!=ds.end(); ++it, ++i)
  {
    OrientedSurface surf = (*it)(p);
    const double distance = xVector(surf.first,p)*surf.second;
    if (distance>max_val) { max_val=distance; max_i=i; }
  }
  return (ds[max_i])(p);
}

OrientedSurface FreePunch::operator()(const Trellis_Util::mPoint& p) const
{
  if (p(0)>=origin(0) && p(1)>=(origin(1)+sqrt(2.)*radius/2.))
  {
    xVector normal(origin,p);
    normal.norm();
    return std::make_pair(Trellis_Util::mPoint(origin(0)+radius*normal(0),origin(1)+radius*normal(1),origin(2)+radius*normal(2)), normal);
  }
  else if (p(0)<origin(0))
  {
    const double distance = xVector(point, p)*normal1;
    return std::make_pair(Trellis_Util::mPoint(p(0)-distance*normal1(0),p(1)-distance*normal1(1),p(2)-distance*normal1(2)), normal1);
  }
  else
  {
    Trellis_Util::mPoint pt(origin(0)+sqrt(2.)*radius/2.,origin(1)+sqrt(2.)*radius/2.,0.);
    const double distance = xVector(pt, p)*normal2;
    return std::make_pair(Trellis_Util::mPoint(p(0)-distance*normal2(0),p(1)-distance*normal2(1),p(2)-distance*normal2(2)), normal2);
  }
}

void EvalContact::operator()(const xGeomElem* geo_appro, const xGeomElem* geo_integ, std::pair<double, xVector>& res) const
{
  Trellis_Util::mPoint pt = geo_integ->getXYZ();
  std::pair<Trellis_Util::mPoint,xVector> oriented_surface = point_to_oriented_surface(pt);
  xVector dir(oriented_surface.first,pt);
  double distance = dir*oriented_surface.second;
  res.first=distance;
  res.second=oriented_surface.second;
}

double Adaptator::operator()(const Trellis_Util::mPoint& p) const
{
  OrientedSurface oriented_surface = point_to_oriented_surface(p);
  xVector vec(oriented_surface.first,p);
  double sign = vec*oriented_surface.second;
  if (sign>0.)
    return xDistance(p,oriented_surface.first);
  else
    return -xDistance(p,oriented_surface.first);
}