a_face_2d4.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_2d4.h"
#include "a_block_2d3.h"
#include "a_block_2d4.h"
#include "a_block_2d5.h"
 
#include "a_debug.h"
 
const double verysmall = 1.e-6;
 
//---------------------------------------------------------------------------
a_face_2d4 * cf_2d4(a_face * f) {return dynamic_cast<a_face_2d4*>(f); }
const a_face_2d4 * cf_2d4(const a_face * f) {return dynamic_cast<const a_face_2d4*>(f); }
//---------------------------------------------------------------------------
double a_face_2d4::thickness() const 
{
        if (this->block_->name()=="a_block_2d3")
                return cb_2d3(this->block_)->thickness();
        if (this->block_->name()=="a_block_2d4")
                return cb_2d4(this->block_)->thickness();
        if (this->block_->name()=="a_block_2d5")
                return cb_2d5(this->block_)->thickness();
        return 0.; //problem!
}
//---------------------------------------------------------------------------
double a_face_2d4::dyM() const
{
        return this->thickness()/2.;
}
//---------------------------------------------------------------------------
double a_face_2d4::dym() const
{
        return this->thickness()/2.;
}
//---------------------------------------------------------------------------
double a_face_2d4::length() const 
{
        a_point * p1 = this->vl(0);
        a_point * p2 = this->vl(1);
        return (*p2-*p1).norm();
}
//---------------------------------------------------------------------------
void a_face_2d4::slide(const double& u, const double& v)
{
        if (v!=0.)
                std::cout << "do not expect much with v different from 0. with this type of face" << std::endl;
        a_face::slide(u,v);
}
//---------------------------------------------------------------------------
double a_face_2d4::S() const 
{
        return (this->length())*(this->thickness());
}
//---------------------------------------------------------------------------
a_point a_face_2d4::cl() const
{
        a_point * p1 = this->vl(0);
        a_point * p2 = this->vl(1);
        return ((*p1)+(*p2))/2.;
}
//---------------------------------------------------------------------------
a_point a_face_2d4::v3d(const int i) const
{
        a_point n = this->ny();//normal to the block
        n *= this->thickness()/2.;
        if (i==0)
                return v(0)-n;
        if (i==1)
                return v(0)+n;
        if (i==2)
                return v(1)+n;
        return v(1)-n;
}
//---------------------------------------------------------------------------
a_point a_face_2d4::nx() const
{
        a_point p1,p2;
        if (block_->name()=="a_block_2d4")
        {
                a_block_2d4 * b = dynamic_cast<a_block_2d4*>(block_);
                p1 = (*(b->pos())) * (*(this->vl(0)));
                p2 = (*(b->pos())) * (*(this->vl(1)));
        }
        if (block_->name()=="a_block_2d3")
        {
                a_block_2d3 * b = dynamic_cast<a_block_2d3*>(block_);
                p1 = (*(b->pos())) * (*(this->vl(0)));
                p2 = (*(b->pos())) * (*(this->vl(1)));
        }
        if (block_->name()=="a_block_2d5")
        {
                a_block_2d5 * b = dynamic_cast<a_block_2d5*>(block_);
                p1 = (*(b->pos())) * (*(this->vl(0)));
                p2 = (*(b->pos())) * (*(this->vl(1)));
        }
        return (p2-p1).normalise();
}
//---------------------------------------------------------------------------
a_point a_face_2d4::ny() const
{
        if (block_->name()=="a_block_2d4")
        {
                a_block_2d4 * b = dynamic_cast<a_block_2d4*>(block_);
                return b->normal();
        }
        if (block_->name()=="a_block_2d3")
        {
                a_block_2d3 * b = dynamic_cast<a_block_2d3*>(block_);
                return b->normal();
        }
        if (block_->name()=="a_block_2d5")
        {
                a_block_2d5 * b = dynamic_cast<a_block_2d5*>(block_);
                return b->normal();
        }
        return a_point(); // problem
}
//---------------------------------------------------------------------------
a_point a_face_2d4::normal() const
{
        a_point px = this->nx();
        a_point py = this->ny();
        return cross(px,py).normalise();
}
//---------------------------------------------------------------------------
int a_face_2d4::ptsinplane(const a_plane& p) const
{
        int n = 0;
        for (int i=0; i<4; i++)
        {
                if (p.contains(this->v3d(i)))
                        n++;
        }
        return n;
}
//---------------------------------------------------------------------------
double a_face_2d4::exrx() const
{
        return 2*(this->exx()/this->length());
}
//---------------------------------------------------------------------------
double a_face_2d4::exry() const
{
        return 2*(this->exy()/this->length());
}
//---------------------------------------------------------------------------
std::vector<a_triangle> a_face_2d4::triangles() const
{
        //find the 4 corners of face
        a_point p01 = v(0);
        a_point p02 = v(1);
        a_point n = this->ny();//normal to the block
        double th = this->thickness()/2.;
        a_point p1 = p01-n*th; 
        a_point p2 = p01+n*th; 
        a_point p3 = p02+n*th; 
        a_point p4 = p02-n*th; 
        a_triangle t1(p1,p2,p3);
        a_triangle t2(p1,p3,p4);
        std::vector<a_triangle> list;
        list.push_back(t1);
        list.push_back(t2);
        return list;
}
//---------------------------------------------------------------------------
int a_face_2d4::contacttype() const
{
        if (lface_==0)
                return -2;
        a_face_2d4 * f2 = cf_2d4(lface_);
        int n1 = this->ptsinplane(f2->plane());
        int n2 = f2->ptsinplane(this->plane());
        int n=n2;
        if (n1>n2)
                n=n1;
        if (n==0)
                return -1;
        if (n>2)
                return 2;
        return n-1;
}
//---------------------------------------------------------------------------
bool a_face_2d4::contains(const a_point& p) const
{
        std::vector<a_triangle> t = this->triangles();
        bool i1 = t[0].contains(p);
        bool i2 = t[1].contains(p);
        return (i1||i2);
}
//---------------------------------------------------------------------------
a_segment a_face_2d4::segment(const int i) const
{
        return a_segment(this->v3d(i%4),this->v3d((i+1)%4));
}
//---------------------------------------------------------------------------
void a_face_2d4::grow(const double& val)
{
        a_point c = this->c();
        a_point d = (val/2.)*this->nx();
        int r0 = r_vertices_[0];
        int r1 = r_vertices_[1];
        a_point p0n = this->v(0)-d;
        a_point p1n = this->v(1)+d;
        this->block()->v(r0)->set(p0n.x(),p0n.y(),p0n.z());
        this->block()->v(r1)->set(p1n.x(),p1n.y(),p1n.z());
}
//---------------------------------------------------------------------------
a_segment a_face_2d4::intersect(const a_segment& s) const
{
        a_segment s1;
        a_segment s2;
        bool t1 = true;
        bool t2 = true;
        std::vector<a_triangle> t = this->triangles();
        try {
                s1 = t[0].intersect(s);
        }
        catch(a_segment::null_segment_error) {
                t1 =false;
        }
        try {
                s2 = t[1].intersect(s);
        }
        catch(a_segment::null_segment_error) {
                t2 =false;
        }
        if (!t1&&!t2)
                throw a_face_2d4::no_contact_error();
        if (t1&&t2)
        {
                a_point v1 = s1.p1();
                a_point v2 = s1.p2();
                a_point v3 = s2.p1();
                a_point v4 = s2.p2();
                if ((v2-v3).norm()<verysmall)
                        return a_segment(v1,v4);
                if ((v1-v3).norm()<verysmall)
                        return a_segment(v2,v4);
                if ((v1-v4).norm()<verysmall)
                        return a_segment(v2,v3);
                return a_segment(v1,v3);
        }
        if (t1) return s1;
        //if (t2) return s2; // test not necessary
        return s2;
}
//---------------------------------------------------------------------------
a_point a_face_2d4::contactpoint() const
{
        int code = this->contacttype();
        if (code != 0)
                throw a_face::bad_contact_type_error();
        a_face_2d4 * f2 = cf_2d4(lface_);
        std::vector<a_triangle> t1 = this->triangles();
        std::vector<a_triangle> t2 = f2->triangles();
        a_point p0;
        for (int i=0; i<2; i++)
        {
                for (int j=0; j<3; j++)
                {
                        a_point p = t1[i][j];
                        double d = t2[0].distp(p);
                        if (d<verysmall)
                        {
                                p0 = p;
                                break;
                        }
                }
        }
        //check if p0 is inside f2
        if (f2->contains(p0))
                return p0;
        else
                throw a_face_2d4::no_contact_error();
}
//---------------------------------------------------------------------------
a_segment a_face_2d4::contactsegment() const
{
        int code = this->contacttype();
        if (code != 1)
                throw a_face::bad_contact_type_error();
        const a_face_2d4 * f2 = cf_2d4(lface_);
        const a_face_2d4 * f;
        int n1 = this->ptsinplane(f2->plane());
        if (n1==2)
                f = this;
        else
        {
                f = f2;
                f2 = this;
        }
        //find edge on f
        std::vector<a_point> pts;
        a_plane p = f2->plane();
        for (int i=0; i<4; i++)
        {
                a_point x = f->v3d(i);
                if (p.contains(x))
                        pts.push_back(x);
        }
        a_segment s1(pts[0],pts[1]);
        a_segment s2 = f2->intersect(s1);
        return s2;
}
//---------------------------------------------------------------------------
a_contact a_face_2d4::contactface() const
{
        int code = this->contacttype();
        if (code != 2)
                throw a_face::bad_contact_type_error();
        a_face_2d4 * f2 = cf_2d4(lface_);
        std::vector<a_triangle> t1 = this->triangles();
        std::vector<a_triangle> t2 = f2->triangles();
        std::vector<a_segment> segs;
        for (int i=0; i<4; i++)
        {
                a_segment s = this->segment(i);
                for (int k=0; k<2; k++)
                {
                        try {
                                a_segment s2 = t2[k].intersect(s);
                                segs.push_back(s2);     
                        }
                        catch(a_segment::null_segment_error) { }
                }
        }
        for (int i=0; i<4; i++)
        {
                a_segment s = f2->segment(i);
                for (int k=0; k<2; k++)
                {
                        try {
                                a_segment s2 = t1[k].intersect(s);
                                segs.push_back(s2);     
                        }
                        catch(a_segment::null_segment_error) { }
                }
        }
        int nsegs = segs.size();
        if (nsegs==0)
                throw a_face_2d4::no_contact_error();
        a_contact faces;
        faces.f(this->f());
        faces.av(segs[0].p1());
        int i = 0;
        int npts = 1;
        while (npts!=nsegs)
        {
                a_point p = segs[i].p2();
                bool stop;
                for (int j=0; j<nsegs; j++)
                {
                        if (j != i)
                        {
                                a_point p2 = segs[j].p1();
                                double d = (p-p2).norm();
                                if (d<verysmall)
                                {
                                        i=j;
                                        faces.av(p2);
                                        npts++;
                                        break;
                                }
                                p2 = segs[j].p2();
                                d = (p-p2).norm();
                                if (d<verysmall)
                                {
                                        i=j;
                                        faces.av(p2);
                                        npts++;
                                        a_point p3 = segs[j].p2();
                                        segs[j].p2(segs[j].p1());
                                        segs[j].p1(p3);
                                        break;
                                }
 
                        }
                }
        }
        faces.clean();
        return faces;
}
/*
// ***************************************************************************
//-operator>>----------------------------------------------------------------
std::istream& operator>> (std::istream& in, a_face& f)
{
        int nv;
        in >> nv;
        f.r_vertices_.clear();
        for (int i=0; i<nv; i++)
        {
                int ref;
                in >> ref;
                f.r_vertices_.push_back(ref);
        }
        return in;
}
//-operator<<----------------------------------------------------------------
std::ostream& operator<< (std::ostream& o, const a_face& f)
{
        std::cout << f.nv() << " ";
        for (int i=0; i<f.nv(); i++)
                std::cout << f.rv(i) << " ";
        return o;
}
*/