a_face.cxx

Go to the documentation of this file.

/*
Copyright 2010-2016 Pierre SMARS (smars@yuntech.edu.tw)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
 
You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "a_face.h"
#include "a_block.h"
 
#include "config.h"
 
const double verysmall = 1.e-8;
 
//---------------------------------------------------------------------------
a_face::a_face() : exit_face_(false), lface_(0) 
{
        r_vertices_.clear();
        criteria_ = 0; 
        f_ = new a_wrench;
        f_original_ = true; 
}
//---------------------------------------------------------------------------
a_face::a_face(a_block * b, int ref_face) : exit_face_(false), lface_(0), ref_(ref_face) 
{
        r_vertices_.clear(); 
        block_=b;
        criteria_ = 0; 
        f_ = new a_wrench;
        f_original_ = true; 
}
//---------------------------------------------------------------------------
a_face::~a_face()
{
        this->clear();
        if (f_original_)
                delete f_;
}
//---------------------------------------------------------------------------
void a_face::clear()
{
        //nothing special to do, just a vector of int
        r_vertices_.clear();
}
//---------------------------------------------------------------------------
const std::string a_face::help()
{
        std::ostringstream o;
        o << "********" << std::endl;
        o << "a_face:" << std::endl;
        o << "********" << std::endl;
        o << "This class is used to model the face of a 'a_block' (see a_block_help for more information)" << std::endl;
        o << "This is a virtual class, it cannot be used directly, see a_face_... for usable classes" << std::endl;
        o << "Commands:" << std::endl;
        o << "--------" << std::endl;
        o << "c: centre of gravity" << std::endl;
        o << "cl: centre of gravity (in local coordinate system)" << std::endl;
        o << "dx: x dimension of face" << std::endl;
        o << "dy: y dimension of face" << std::endl;
        o << "dxM: maximum x-distance between centre of gravity and vertices" << std::endl;
        o << "dxm: minimum x-distance between centre of gravity and vertices" << std::endl;
        o << "dyM: maximum y-distance between centre of gravity and vertices" << std::endl;
        o << "dym: minimum y-distance between centre of gravity and vertices" << std::endl;
        o << "point 'u' 'v': get point given excentricities u and v" << std::endl;
        o << "exit 'bool': set face as an exit face" << std::endl;
        o << "exit: check whether face is an exit face" << std::endl;
        o << "S: surface" << std::endl;
        o << "Sx: projection of surface perpendicular to x" << std::endl;
        o << "Sy: projection of surface perpendicular to y" << std::endl;
        o << "Sz: projection of surface perpendicular to z" << std::endl;
        o << "nx: direction x" << std::endl;
        o << "ny: direction y" << std::endl;
        o << "normal: normal to the face" << std::endl;
        o << "normale 'i': normal to edge i" << std::endl;
        o << "plane: gets a a_plane" << std::endl;
        o << "N: normal force" << std::endl;
        o << "Qx: shear force" << std::endl;
        o << "Qy: shear force" << std::endl;
        o << "Mx: torque (bending)" << std::endl;
        o << "My: torque (bending)" << std::endl;
        o << "Mt: torque (torsion)" << std::endl;
        o << "ax: angle of incidence" << std::endl;
        o << "ay: angle of incidence" << std::endl;
        o << "x: application point of f on face" << std::endl;
        o << "exx: excentricity in x direction" << std::endl;
        o << "exy: excentricity in y direction" << std::endl;
        o << "exn: excentricity in normal direction (should be 0)" << std::endl;
        o << "nv: get number of vertices" << std::endl;
        o << "av 'ref': add reference to a vertex" << std::endl;
        o << "rv 'i': gets vertex reference" << std::endl;
        o << "block: gets block to which the face belongs" << std::endl;
        o << "f 'a_wrench': apply wrench to face" << std::endl;
        o << "f 'a_point': apply force to centre of gravity of face" << std::endl;
        o << "f 'fx' 'fy' 'fz': apply force to centre of gravity of face" << std::endl;
        o << "f: gets force applied to face" << std::endl;
        o << "lface: gets connected face" << std::endl;
        o << "lblock: gets connected block" << std::endl;
        o << "in: face is connected to other faces, inside the structure" << std::endl;
        o << "out: face is not connected to other faces, on the boudary" << std::endl;
        o << "hinge: hinge opening: rotation vector" << std::endl;
        o << "contacttype: contact type: (-2: no sister face, -1: not touching, 2: plane contact, 1: edge contact, 0: point contact)" << std::endl;
        o << "contactpoint: return point of contact if faces just touch by one point" << std::endl;
        o << "contactsegment: return edge of contact if faces touch by one edge: hinge" << std::endl;
        o << "contactface: return common surface of contact" << std::endl;
        o << "grow 'val': incread the size of the face" << std::endl;
        o << "slide 'u' 'v': slide face in its own plane (values are relative excentricities)" << std::endl;
        o << "deform 'i' 'val': move edge 'i' perpendicularly to itself" << std::endl;
        return o.str();
}
//---------------------------------------------------------------------------
a_point a_face::vc(const int i) const
{
        return this->v(i)-this->c();
}
//---------------------------------------------------------------------------
double a_face::dxM() const
{
        double dmax = 0.;
        a_point c = this->c();
        for (unsigned short i=0; i<this->nv(); i++)
        {
                double d;
                d = (this->vc(i))*(this->nx());
                if (d>dmax)
                        dmax = d;
        }
        return dmax;
}
//---------------------------------------------------------------------------
double a_face::dxm() const
{
        double dmin = 0.;
        for (unsigned short i=0; i<this->nv(); i++)
        {
                double d;
                d = (this->vc(i))*(this->nx());
                if (d<dmin)
                        dmin = d;
        }
        return -dmin;
}
//---------------------------------------------------------------------------
double a_face::dyM() const
{
        double dmax = 0.;
        for (unsigned short i=0; i<this->nv(); i++)
        {
                double d;
                d = (this->vc(i))*(this->ny());
                if (d>dmax)
                        dmax = d;
        }
        return dmax;
}
//---------------------------------------------------------------------------
double a_face::dym() const
{
        double dmin = 0.;
        for (unsigned short i=0; i<this->nv(); i++)
        {
                double d;
                d = (this->vc(i))*(this->ny());
                if (d<dmin)
                        dmin = d;
        }
        return -dmin;
}
//---------------------------------------------------------------------------
double a_face::dx() const
{
        return this->dxM()+this->dxm();
}
//---------------------------------------------------------------------------
double a_face::dy() const
{
        return this->dyM()+this->dym();
}
//---------------------------------------------------------------------------
a_point a_face::normale(const unsigned short ref)
{
        if (ref>this->nv()) throw a_face::bad_indexing_error();
  a_point c = this->c();
  a_point p0 = this->v(ref);
        if (this->nv()==2)
        {
                //we are then in 2D, edge is a point
                return (p0-c).normalise();
        }
        else
        {
                //we are then in 3D, find perp to edge
                int ref1 = (ref+1)%this->nv();
        a_point p1 = this->v(ref1);
                //direction of the edge
                a_point l = p1-p0;
                //normal to the face
                a_point n = this->normal();
                a_point ne = cross(n,l);
                //direction has to be checked...
                return ne.normalise();
        }
}
//---------------------------------------------------------------------------
a_point a_face::point(const double& u, const double& v) const
{
  return this->c()+(u*this->dx()/2.)*this->nx()+(v*this->dy()/2.)*this->ny();
}
//---------------------------------------------------------------------------
void a_face::slide(const double& u, const double& v)
{
        a_point d = (u*this->dx()/2.)*this->nx()+(v*this->dy()/2.)*this->ny();
        for (unsigned short i=0; i<this->nv(); i++)
        {
                a_point pn = this->v(i)+d;
                int r = r_vertices_[i];
        this->block()->v(r)->set(pn.x(),pn.y(),pn.z());
        }
}
//---------------------------------------------------------------------------
void a_face::deform(const unsigned short& ref, const double& val)
{
        if (ref>this->nv()) throw a_face::bad_indexing_error();
        int r0 = r_vertices_[ref];
        a_point d = val*this->normale(ref);
        a_point p0n = this->v(ref)+d;
        this->block()->v(r0)->set(p0n.x(),p0n.y(),p0n.z());
        if (this->nv()>2)
  {
                //we are then in 3D, next point also need to be moved
                int ref1 = (ref+1)%this->nv();
                int r1 = r_vertices_[ref1];
        a_point p1n = this->v(ref1)+d;
        this->block()->v(r1)->set(p1n.x(),p1n.y(),p1n.z());
        }
}
//---------------------------------------------------------------------------
bool a_face::ok() const
{
        //if no criteria was linked, should be always good!
        if (criteria_==0)
                return true;
        else
                return criteria_->ok(this);
}
//---------------------------------------------------------------------------
double a_face::penalty() const
{
        //if no criteria was linked, should be always good!
        if (criteria_==0)
                return 0.; 
        else
                return criteria_->penalty(this);
}
//---------------------------------------------------------------------------
void a_face::f(a_wrench * f)
{
//      if (f_ != 0)
//              delete f_;
        f_original_ = false;
        f_ = f;
}
//---------------------------------------------------------------------------
double a_face::Sx() const 
{
        double c = fabs(this->normal().x());
        return (this->S())*c;
}
//---------------------------------------------------------------------------
double a_face::Sy() const 
{
        double c = fabs(this->normal().y());
        return (this->S())*c;
}
//---------------------------------------------------------------------------
double a_face::Sz() const 
{
        double c = fabs(this->normal().z());
        return (this->S())*c;
}
//---------------------------------------------------------------------------
a_point a_face::v(const int i) const 
{
        return (*(block_->pos())) * (*(this->vl(i)));
}
//---------------------------------------------------------------------------
int a_face::rv(const int i) const 
{
        if ((i<0)||(i>r_vertices_.size()-1))
                return -1;
        return r_vertices_[i];
}
//---------------------------------------------------------------------------
a_point * a_face::vl(const int i) 
{
        return block_->v(this->rv(i));
}
//---------------------------------------------------------------------------
a_point * a_face::vl(const int i) const 
{
        return block_->v(this->rv(i));
}
//---------------------------------------------------------------------------
a_point a_face::c() const 
{
        return (*block_->pos())*this->cl();
}
//---------------------------------------------------------------------------
void a_face::f(const a_point& f)
{
        a_point c = this->c();
        f_->set(c,f);
}
//---------------------------------------------------------------------------
void a_face::f(const double fx, const double fy, const double fz)
{
        a_point f(fx,fy,fz);
        this->f(f);
}
//---------------------------------------------------------------------------
double a_face::N() const
{
        a_point f = f_->p1();
        return f*this->normal();
}
//---------------------------------------------------------------------------
double a_face::N(const a_wrench& fe) const
{
        a_point f = fe.p1();
        return f*this->normal();
}
//---------------------------------------------------------------------------
double a_face::Qx() const
{
        a_point f = f_->p1();
        return -f*this->nx();
}
//---------------------------------------------------------------------------
double a_face::Qy() const
{
        a_point f = f_->p1();
        return f*this->ny();
}
//---------------------------------------------------------------------------
double a_face::Mt() const
{
        a_point m = f_->M();
        return m*this->normal();
}
//---------------------------------------------------------------------------
double a_face::Mx() const
{
        a_point m = f_->M();
        m -= cross(this->c(),f_->F());
        return m*this->nx();
}
//---------------------------------------------------------------------------
double a_face::Mx(const a_wrench& f) const
{
        a_point m = f.M();
        m -= cross(this->c(),f.F());
        return m*this->nx();
}
//---------------------------------------------------------------------------
double a_face::My() const
{
        a_point m = f_->M();
        m -= cross(this->c(),f_->F());
        return m*this->ny();
}
//---------------------------------------------------------------------------
double a_face::My(const a_wrench& f) const
{
        a_point m = f.M();
        m -= cross(this->c(),f.F());
        return m*this->ny();
}
//---------------------------------------------------------------------------
double a_face::ax() const
{
        double n = this->N();   
        double qx = this->Qx();
        return atan2(-qx,-n);   
}
//---------------------------------------------------------------------------
double a_face::ay() const
{
        double n = this->N();   
        double qy = this->Qx();
        return atan2(-qy,-n);   
}
//---------------------------------------------------------------------------
a_point a_face::x() const
{
        a_point c = this->c();
        double N = this->N();
        double Mx = this->Mx();
        double My = this->My();
        double x = 0.;
        double y = 0.;
        if (fabs(N)>verysmall)
        {
                x = -My/N;
                y = Mx/N;
        }
        return c+x*this->nx()+y*this->ny();
}
//---------------------------------------------------------------------------
a_point a_face::x(const a_wrench& f) const
{
        a_point c = this->c();
        double N = this->N(f);
        double Mx = this->Mx(f);
        double My = this->My(f);
        double x = 0.;
        double y = 0.;
        if (fabs(N)>verysmall)
        {
                x = -My/N;
                y = Mx/N;
        }
        return c+x*this->nx()+y*this->ny();
}
//---------------------------------------------------------------------------
double a_face::exx() const
{
        a_point x = this->x()-this->c();
        a_point d = this->nx();
        return d*x;
}
//---------------------------------------------------------------------------
double a_face::exx(const a_wrench& f) const
{
        a_point x = this->x(f)-this->c();
        a_point d = this->nx();
        return d*x;
}
//---------------------------------------------------------------------------
double a_face::exy() const
{
        a_point x = this->x()-this->c();
        a_point d = this->ny();
        return d*x;
}
//---------------------------------------------------------------------------
double a_face::exy(const a_wrench& f) const
{
        a_point x = this->x(f)-this->c();
        a_point d = this->ny();
        return d*x;
}
//---------------------------------------------------------------------------
double a_face::exn() const
{
        a_point x = this->x()-this->c();
        a_point d = this->normal();
        return d*x;
}
//---------------------------------------------------------------------------
a_block * a_face::lblock() const
{
        if (lface_ == 0)
                return 0;
        else
                return lface_->block();
}
//---------------------------------------------------------------------------
a_point a_face::hinge() const
{
        a_point n1 = this->normal();
        a_face * f2 = this->lface();
        if (f2==0)
                return a_point(0,0,0);
        a_point n2 = f2->normal();
        //takes care of numerical errors
        double v = -n1*n2;
        if (v>1.) v=1.;
        if (v<-1.) v=-1.;
        double angle = acos(v);
        a_point n = cross(n1,-n2);
        return n*angle;
}
//***************************************************************************
//---------------------------------------------------------------------------
void a_face::read(std::istream &in)
{
        int nv;
        in >> nv;
        r_vertices_.clear();
        for (int i=0; i<nv; i++)
        {
                int ref;
                in >> ref;
                r_vertices_.push_back(ref);
        }
}
//---------------------------------------------------------------------------
void a_face::write(std::ostream &o) const
{
        o << this->nv() << " ";
        for (int i=0; i<this->nv(); i++)
                o << this->rv(i) << " ";
        o << std::endl;
}
//-operator>>----------------------------------------------------------------
std::istream& operator>> (std::istream& in, a_face& f)
{
        f.read(in);
        return in;
}
//-operator<<----------------------------------------------------------------
std::ostream& operator<< (std::ostream& o, const a_face& f)
{
        f.write(o);
        return o;
}