a_rib.cxx

Go to the documentation of this file.

/*
         Copyright 2002-2019 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_rib.h"
#include <sstream>
#include <iomanip>
 
// .PORTABILITY : ansi C++
 
//***************************************************************************
//---------------------------------------------------------------------------
a_rib::a_rib()
{
}
//---------------------------------------------------------------------------
a_rib::~a_rib()
{
        for (auto prof:profiles_)
                delete prof;
        profiles_.clear();
}
//---------------------------------------------------------------------------
void a_rib::init()
{
        this->construct_profile_pos();//find the position of the profiles along the path
        a_point z = path_.z_axis();
        a_point z1 = profiles_[0]->z_axis();//z axis of the first profile
        double a = profile_pos_[0];//position of the first profile on the path
        a_point z2 = path_.dz_axis(a);//z axis at the origin
        if (z1*z2<0)//not the same direction
                this->invert();
        //some new profiles are possibly made to have a complete set of data for processing
        if (profiles_.size()==1)// just 1 profile -> make 1 or 2 new profiles at the extremities
                this->extrapolate_1profile_setup();
        else //if no profile at extremities, extrapolate new ones
                this->extrapolate_profiles_setup();
        if (z1*z2<0)//not the same direction
                this->invert();//necessary to invert the new profiles
        //std::cout << *this;
        //std::exit(0);
}
//---------------------------------------------------------------------------
void a_rib::invert()
{
        int n_pat = profiles_.size();
        for (int i = 0; i < n_pat; i++)//for each profile
                profiles_[i]->invert();
}
//---------------------------------------------------------------------------
int a_rib::ref(double a) const//find the previous profile for a given relative position along the path
{
        int i = -1;
        double a2;
        do
                a2 = profile_pos_[++i];//profiles' relative positions
        while ((a2 <= a)&&(i<profile_pos_.size()));
        return i-1;
}//---------------------------------------------------------------------------
void a_rib::construct_profile_pos()
{
        int n_pat = profiles_.size();
        int ref0 = -1;
        for (int i = 0; i < n_pat; i++)//for each profile
        {
                double m1;
                a_curve_lin * pat = profiles_[i];
                int ref1;
                //find the reference of the point on the path closest to the profile
                //& find the reference of the point on the profile closest to the path
                double dist = path_.dist2(*pat,ref1,path_pos_,m1,m2_);
                if (ref1<ref0)
                {
                        int test = path_.size();
                        if (ref1==0)
                                profile_pos_.push_back(1);
                        else if (ref1==path_.size()-2)
                                profile_pos_.push_back(0);
                        else //for closed paths
                        {
                                std::cerr << "profiles should be ordered along the path" << std::endl;
                                std::cerr << "ref0 " << ref0 << std::endl;
                                std::cerr << "ref1 " << ref1 << std::endl;
                                throw profile_order_error();
                        }
                }
                else
                {
                        double a1 = path_.a()[ref1];
                        double a2 = path_.a()[ref1+1];
                        profile_pos_.push_back(a1+m1*(a2-a1));
                }
                ref0 = ref1;
        }
}
//---------------------------------------------------------------------------
double gia_f(double a, double a1, double a2, bool path_is_closed)
{
        if (path_is_closed)
        {
                if (a1<a2)
                        return (a-a1)/(a2-a1);
                else
                {
                        double base = (1-a1)+a2;
                        if (a<=a2)
                                return 1.+(a-a2)/base;
                        else
                                return (a-a1)/base;
                }
        }
        else
                return (a-a1)/(a2-a1);
}
//---------------------------------------------------------------------------
a_curve_lin a_rib::gia_profile(double a, int ref1, int ref2) const
{
        a_curve_lin pat1 = *profiles_[ref1];//the 2 profiles to interpolate
        a_curve_lin pat2 = *profiles_[ref2];
 
        double a1 = profile_pos_[ref1];//position of the profile on the path
        double a2 = profile_pos_[ref2];
 
        double m1,m2;
        int r1,r2;
        int rz1,rz2;
        double mz1,mz2;
 
        a_point gap1_3D = path_.shortest(pat1,r1,rz1,m1,mz1).vec();//shortest segment between the path & profile
        a_point gap2_3D = path_.shortest(pat2,r1,rz2,m1,mz2).vec();
 
        a_point x1 = pat1.x_axis();//reference system of the 1st profile
        a_point y1 = pat1.y_axis();
        a_point z1 = pat1.z_axis();
 
        a_point gap1_2D(x1*gap1_3D,y1*gap1_3D,0.);//find the gap coordinates in the plane of the profile
 
        a_point dx1 = path_.dx_axis(a1).rotate(x1,y1,z1);//angles with the reference system on the path at a1
        a_point dy1 = path_.dy_axis(a1).rotate(x1,y1,z1);
        a_point dz1 = path_.dz_axis(a1).rotate(x1,y1,z1);
 
        a_point x2 = pat2.x_axis();//reference system of the 2d profile
        a_point y2 = pat2.y_axis();
        a_point z2 = pat2.z_axis();
 
        a_point gap2_2D(x2*gap2_3D,y2*gap2_3D,0.);//find the gap coordinates in the plane of the profile
 
        a_point dx2 = path_.dx_axis(a2).rotate(x2,y2,z2);//angles with the reference system on the path at a2
        a_point dy2 = path_.dy_axis(a2).rotate(x2,y2,z2);
        a_point dz2 = path_.dz_axis(a2).rotate(x2,y2,z2);
 
        double f = gia_f(a,a1,a2,path_.is_closed());//relative position of the 2 profiles
 
        a_point dx3 = average_rot(dx1,dx2,f);//interpolated angle with the reference system on the path at a
        a_point dy3 = average_rot(dy1,dy2,f);
        a_point dz3 = average_rot(dz1,dz2,f);
 
        pat1.flatten2D();//the 2 profiles are flatten in 2D
        pat2.flatten2D();
 
        a_curve_lin pat = interpolate(pat1,pat2,f);//an interpolated flatten path is found
 
        a_point ref1_2D = *pat1[rz1]+mz1*(*pat1[rz1+1]-*pat1[rz1]);//vector between the origin of the profile and the closest point to the path
        a_point ref2_2D = *pat2[rz2]+mz2*(*pat2[rz2+1]-*pat2[rz2]);
 
        a_point gap_2D = average(gap1_2D,gap2_2D,f);
        a_point ref_2D =  average(ref1_2D,ref2_2D,f);
 
        a_point p = path_(a);//position of the new profile on the path
 
        a_point x0 = path_.dx_axis(a);
        a_point y0 = path_.dy_axis(a);
        a_point z0 = path_.dz_axis(a);
 
        a_point x = (x0*dx3.x()+y0*dy3.x()+z0*dz3.x()).normalise();
        a_point y = (x0*dx3.y()+y0*dy3.y()+z0*dz3.y()).normalise();
 
        a_point ref_3D =  ref_2D.x()*x+ref_2D.y()*y;//place the interpolated reference shift in the plane of the new profile
        a_point gap_3D = gap_2D.x()*x+gap_2D.y()*y;//place the interpolated gap in the plane of the new profile
 
        a_point prov = p+gap_3D-ref_3D;
        pat.place3D(prov,x,y);
 
        pat.set_horiz(path_.get_horiz());
        pat.construct_axis();
 
        return pat;
}
//---------------------------------------------------------------------------
void a_rib::extrapolate_1profile_setup()
{
        a_curve_lin pat = *profiles_[0];//the profile to extrapolate
 
        double a = profile_pos_[0];//position of the profile on the path
 
        double m1,m2;
        int r1,r2;
        a_point gap_3D = path_.shortest(pat,r1,r2,m1,m2).vec();
 
        a_point x = pat.x_axis();//reference system of the profile
        a_point y = pat.y_axis();
        a_point z = pat.z_axis();
 
        a_point gap_2D(x*gap_3D,y*gap_3D,0.);//find the gap coordinates in the plane of the profile
 
        a_point dx = path_.dx_axis(a).rotate(x,y,z);//angles with the reference system on the path at a
        a_point dy = path_.dy_axis(a).rotate(x,y,z);
        a_point dz = path_.dz_axis(a).rotate(x,y,z);
        if (a > 0)//there is no profile at the path's beginning
        {
                auto patn = new a_curve_lin(pat);
                patn->flatten2D();
                a_point ref_2D = *(*patn)[r2]+m2*(*(*patn)[r2+1]-*(*patn)[r2]);//vector between the origin of the profile and the closest point to the path
                a_point pa = path_(0.);//position of the profile on the path
 
                a_point x0 = path_.dx_axis(0.);
                a_point y0 = path_.dy_axis(0.);
                a_point z0 = path_.dz_axis(0.);
 
                a_point dxa = (x0*dx.x()+y0*dy.x()+z0*dz.x()).normalise();
                a_point dya = (x0*dx.y()+y0*dy.y()+z0*dz.y()).normalise();
 
                a_point ref_3D =  ref_2D.x()*dxa+ref_2D.y()*dya;//place the interpolated reference shift in the plane of the new profile
                a_point gapn_3D = gap_2D.x()*dxa+gap_2D.y()*dya;//place the interpolated gap in the plane of the new profile
                a_point prov = pa+gap_3D-ref_3D;
                patn->place3D(prov,dxa,dya);
                //patn->set_horiz(path_.get_horiz());
                //patn->construct_axis();
                profiles_.insert(profiles_.begin(),patn);
                profile_pos_.insert(profile_pos_.begin(),0.);
        }
 
        if ((a < 1)&&(!path_.is_closed()))//there is no profile at the path's end
        {
                auto patn = new a_curve_lin(pat);
                patn->flatten2D();
                a_point ref_2D = *(*patn)[r2]+m2*(*(*patn)[r2+1]-*(*patn)[r2]);//vector between the origin of the profile and the closest point to the path
                a_point pb = path_(1.);//position of the profile on the path
 
                a_point x0 = path_.dx_axis(1.);
                a_point y0 = path_.dy_axis(1.);
                a_point z0 = path_.dz_axis(1.);
 
                a_point dxb = (x0*dx.x()+y0*dy.x()+z0*dz.x()).normalise();
                a_point dyb = (x0*dx.y()+y0*dy.y()+z0*dz.y()).normalise();
 
                //              a_point dxb = path_.dx_axis(1.).rotate(dx,dy,dz).normalise();
                //              a_point dyb = path_.dy_axis(1.).rotate(dx,dy,dz).normalise();
 
                a_point ref_3D =  ref_2D.x()*dxb+ref_2D.y()*dyb;//place the interpolated reference shift in the plane of the new profile
                a_point gapn_3D = gap_2D.x()*dxb+gap_2D.y()*dyb;//place the interpolated gap in the plane of the new profile
                a_point prov = pb+gap_3D-ref_3D;
                patn->place3D(prov,dxb,dyb);
                //patn->set_horiz(path_.get_horiz());
                //patn->construct_axis();
                profiles_.push_back(patn);
                profile_pos_.push_back(1.);
        }
}
//---------------------------------------------------------------------------
void a_rib::extrapolate_profiles_setup()
{
        if (!path_.is_closed())
        {
                if (profile_pos_[0] != 0)//no profile at the beginning
                {
                        auto pat = new a_curve_lin(this->gia_profile(0,0,1));
                        profiles_.insert(profiles_.begin(),pat);
                        profile_pos_.insert(profile_pos_.begin(),0.);
                }
                int n_pat = profiles_.size();
                if (profile_pos_[n_pat-1] < 1)//no profile at the end
                {
                        auto pat = new a_curve_lin(this->gia_profile(1,n_pat-2,n_pat-1));
                        profiles_.push_back(pat);
                        profile_pos_.push_back(1.);
                }
        }
}
//---------------------------------------------------------------------------
a_curve_lin a_rib::profile(double a) const
{
        if (path_.is_closed())
        {
                if ((a>=profile_pos_[0])&&(a<=*(profile_pos_.end()-1)))
                {
                        int ref = this->ref(a);
                        if (ref==profiles_.size()-1)
                                ref--;
                        return this->gia_profile(a,ref,ref+1);
                }
                else 
                        return this->gia_profile(a,profiles_.size()-1,0);
        }
        else
        {
                int ref = this->ref(a);
                if (ref==profiles_.size()-1)
                        ref--;
                return this->gia_profile(a,ref,ref+1);
        }
}
//---------------------------------------------------------------------------
void a_rib::output_pts(std::ostream& o, int& imax, int& jmax, bool is_closed_volume) const
{
        imax = path_.size();
        jmax = profiles_[0]->size();
        if (path_.is_closed())
        {
                imax--;
                o << imax*jmax << std::endl;//number of points + 2 in the middle of the end faces
        }
        else
        {
                if (is_closed_volume) //closed triangle volume will be produced
                        o << imax*jmax+2 << std::endl;//number of points + 2 in the middle of the end faces
                else //open triangle volume will be produced
                        o << imax*jmax << std::endl;//number of points
        }
 
        //points 
        std::ostringstream out;
        for (int i = 0; i < imax; i++)
        {
                a_point * p0 = path_[i];//a point along the path
                double a = path_.a()[i];//its relative position (0-1)
                a_curve_lin pat = this->profile(a);//the profile at that place
                if (((i==0)||(i==imax-1))&&(is_closed_volume))//extreme section and closed set
                {
                        // a 'centre' point is computed for the extremity's sections
                        a_point mid;
                        for (int j = 0; j < jmax; j++)//write the points of the profile
                        {
                                a_point * x = pat[j];
                                mid += *x;
                                o << *x << std::endl;
                        }
                        mid /= jmax;
                        out << mid << std::endl;//mid point is put in a provisory place
                }
                else
                {
                        for (int j = 0; j < jmax; j++)//write the points of the profile
                                o << *pat[j] << std::endl;
                }
        }
        if ((is_closed_volume)&&(!path_.is_closed()))
                o << out.str();
}
//---------------------------------------------------------------------------
void a_rib::output_pts_vol(std::ostream& o, int& imax, int& jmax) const
{
        imax = path_.size();
        jmax = profiles_[0]->size();
 
        if (path_.is_closed())
                imax--;
        o << imax*(jmax+1) << std::endl;
        //points 
        for (int i = 0; i < imax; i++)
        {
                a_point * p0 = path_[i];//a point along the path
                double a = path_.a()[i];//its relative position (0-1)
                a_curve_lin pat = this->profile(a);//the profile at that place
                // a 'centre' point is computed
                a_point mid;
                for (int j = 0; j < jmax; j++)//write the points of the profile
                {
                        a_point * x = pat[j];
                        mid += *x;
                        o << *x << std::endl;
                }
                mid /= jmax;
                o << mid << std::endl;//mid point
        }
}
//---------------------------------------------------------------------------
void a_rib::output_tri_open(std::ostream& o) const
{
        int imax, jmax;
        bool is_closed_volume = false;//open structure
        this->output_pts(o,imax,jmax,is_closed_volume);
        //triangles
        if (path_.is_closed())
                o << imax*(jmax-1)*2 << std::endl;//number of triangles
        else
                o << (imax-1)*(jmax-1)*2 << std::endl;//number of triangles
        for (int i = 0; i < imax-1; i++)//write the topology
        {
                for (int j = 0; j < jmax-1; j++)
                {
                        int i1 = i*jmax+j;                      int i2 = i1+1;
                        int j1 = i1+jmax;                       int j2 = j1+1;
                        o << i1 << "\t" << j2 << "\t" << j1 << std::endl;
                        o << i1 << "\t" << i2 << "\t" << j2 << std::endl;
                }
        }
        if (path_.is_closed())//path is closed
        {
                for (int j = 0; j < jmax-1; j++)
                {
                        int i1 = (imax-1)*jmax+j;       int i2 = i1+1;
                        int j1 = j;                                     int j2 = j1+1;
                        o << i1 << "\t" << j2 << "\t" << j1 << std::endl;
                        o << i1 << "\t" << i2 << "\t" << j2 << std::endl;
                }
        }
}
//---------------------------------------------------------------------------
void a_rib::output_tri_closed(std::ostream& o) const
{
        int imax, jmax;
        bool is_closed_volume = true;
        this->output_pts(o,imax,jmax,is_closed_volume);
        //triangles
        o << imax*jmax*2 << std::endl;//number of triangles
        for (int i = 0; i < imax-1; i++)//write the topology
        {
                for (int j = 0; j < jmax-1; j++)
                {
                        int i1 = i*jmax+j;                      int i2 = i1+1;
                        int j1 = i1+jmax;                       int j2 = j1+1;
                        o << i1 << "\t" << j2 << "\t" << j1 << std::endl;
                        o << i1 << "\t" << i2 << "\t" << j2 << std::endl;
                }
                int i1 = (i+1)*jmax-1;                  int i2 = i*jmax;
                int j1 = i1+jmax;                               int j2 = i2+jmax;
                o << i1 << "\t" << j2 << "\t" << j1 << std::endl;
                o << i1 << "\t" << i2 << "\t" << j2 << std::endl;
        }
        if (path_.is_closed())
        {
                //triangles connecting the sections of the extremities
                for (int j = 0; j < jmax-1; j++)
                {
                        int i1 = (imax-1)*jmax+j;       int i2 = i1+1;
                        int j1 = j;                                     int j2 = j1+1;
                        o << i1 << "\t" << j2 << "\t" << j1 << std::endl;
                        o << i1 << "\t" << i2 << "\t" << j2 << std::endl;
                }
                int i1 = imax*jmax-1;                   int i2 = i1-jmax+1;
                int j1 = jmax-1;                                int j2 = 0;
                o << i1 << "\t" << j2 << "\t" << j1 << std::endl;
                o << i1 << "\t" << i2 << "\t" << j2 << std::endl;
        }
        else//path is open
        {
                //triangles closing the sections of the extremities
                for (int j = 0; j < jmax-1; j++)
                {
                        o << imax*jmax << "\t" << j+1 << "\t" << j << std::endl;
                        o << imax*jmax+1 << "\t" << (imax-1)*jmax+j << "\t" << (imax-1)*jmax+j+1 << std::endl;
                }
                o << imax*jmax << "\t" << 0 << "\t" << jmax-1 << std::endl;
                o << imax*jmax+1 << "\t" << (imax-1)*jmax+jmax-1 << "\t" << (imax-1)*jmax << std::endl;
        }
}
//---------------------------------------------------------------------------
void a_rib::output_vol(std::ostream& o) const
{
        int imax, jmax;
        bool is_closed_volume = true;
        this->output_pts_vol(o,imax,jmax);
        //prisms
        o << (imax-1)*jmax << std::endl;//number of prisms
        for (int i = 0; i < imax-1; i++)//write the topology
        {
                int k1 = i*(jmax+1)+jmax;
                int k2 = (i+1)*(jmax+1)+jmax;
                for (int j = 0; j < jmax; j++)
                {
                        int i1 = i*(jmax+1)+j;
                        int i2;
                        if (j==jmax-1)
                                i2 = i*(jmax+1);
                        else
                                i2 = i1 + 1;
                        int j1 = i1 + jmax + 1;
                        int j2 = i2 + jmax + 1;
                        o << k1 << "\t" << i1 << "\t" << i2 << "\t" << k2 << "\t" << j1 << "\t" << j2 << std::endl;
                }
        }
        if (path_.is_closed())
        {
                o << "to complete" << std::endl;
        }
}
//***************************************************************************
//---------------------------------------------------------------------------
std::istream& operator>> (std::istream& in, a_rib & r)
{
        in >> r.path_;//read the path
        a_point h = r.path_.z_axis();
        h.z(0);
        r.path_.set_horiz(h);
 
        int n_pats;
        in >> n_pats;//read all the profiles
 
        for (int i = 0; i < n_pats; i++)
        {
                auto pat = new a_curve_lin;
                in >> *pat;
                pat->set_horiz(h);
                pat->construct_axis();
                r.profiles_.push_back(pat);
        }
        r.init();
        return in;
}
//---------------------------------------------------------------------------
std::ostream& operator<< (std::ostream& o, const a_rib & r)
{
        o << r.path_;
        int n_pats = r.profiles_.size();
        o << n_pats << std::endl;
        for (int i = 0; i < n_pats; i++)
                o << *r.profiles_[i];
        return o;
}