a_structure.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     <math.h>           // math routines
 
#include "a_structure.h"
#include "a_block_2d4.h"
#include "a_face_2d4.h"
#include <vnl/vnl_vector.h>
#include <vnl/vnl_matrix.h>
#include <vnl/algo/vnl_svd.h>
#include <vnl/vnl_least_squares_function.h>
#include <vnl/algo/vnl_levenberg_marquardt.h>
#include <vnl/vnl_cost_function.h>
#include <vnl/algo/vnl_levenberg_marquardt.h>
#include <vnl/vnl_least_squares_cost_function.h>
#include <vnl/algo/vnl_amoeba.h>
#include <vnl/algo/vnl_powell.h>
#include <vnl/algo/vnl_conjugate_gradient.h>
#include <vnl/algo/vnl_lbfgs.h>
 
#include "config.h"
 
//---------------------------------------------------------------------------
a_structure::a_structure()
{
        objective_ = 0;
        penalty_factor_ = 1.;
}
//---------------------------------------------------------------------------
a_structure::~a_structure()
{
        this->clear();
}
//---------------------------------------------------------------------------
void a_structure::clear()
{
        for(std::vector<a_block *>::iterator it = blocks_.begin(); it != blocks_.end(); ++it) (*it)->clear();
        blocks_.clear();
}
//---------------------------------------------------------------------------
const std::string a_structure::help()
{
        std::ostringstream o;
        o << "************" << std::endl;
        o << "a_structure:" << std::endl;
        o << "************" << std::endl;
        o << "This class is used to model structures made of blocks" << std::endl;
        o << "Commands:" << std::endl;
        o << "--------" << std::endl;
        o << "a_structure S: create a structure named S" << std::endl;
        o << "" << std::endl;
        o << "> building the structure" << std::endl;
        o << "......................" << std::endl;
        o << "clear: clear completely the structure" << std::endl;
        o << "add_block 'block_object': add a block" << std::endl;
        o << "link_blocks 'b1' 'f1' 'b2' 'f2': link blocks 'b1' and 'b2' by their faces 'f1' and 'f2'" << std::endl;
        o << "material 'material': link a material (can also be defined block by block)" << std::endl;
        o << "criteria 'criteria': set a global resistance criteria for the structure (can also be specified block by block)" << std::endl;
        o << "penalty_factor 'value': set penalty factor" << std::endl;
        o << "b 'ref': returns block 'ref' (possibly used to specify parameters for this block)" << std::endl;
        o << "rotate 'b' 'f' 'p' 'a' 'r': rotate 'b' and all connected blocks of an angle 'r' (rad), starting from face 'f', using axis 'a'" << std::endl;
        o << "rotate 'b' 'f' 'e' 'r': rotate 'b' and all connected blocks of an angle 'r' (rad), starting from face 'f' at relative distance (-1 - 1) 'e' from the face center, using ny as an axis" << std::endl;
        o << "slide 'b' 'f' 'u' 'v': slide of ('u','v') all blocks connected to 'b' in the plane of 'f'" << std::endl;
        o << std::endl;
        o << "> export data" << std::endl;
        o << "............." << std::endl;
        o << "convert tri 'filename': write triangle cloud (filename is optional)" << std::endl;
        o << "convert blender 'filename': write structure in  blender format (filename is optional)" << std::endl;
        o << "convert brlcad 'filename': write structure in  brlcad format (filename is optional)" << std::endl;
        o << "trianglecloud 'a_trianglecloud': write structure in a 'a_trianglecloud' object (points&forces)" << std::endl;
        o << "polyline 'a_linecloud': write line of resistance in a 'a_linecloud' object (points&forces)" << std::endl;
        o << std::endl;
        o << "> query the structure" << std::endl;
        o << "....................." << std::endl;
        o << "nb: returns number of blocks in the structure" << std::endl;
        o << "rb 'block_object': returns reference number of block" << std::endl;
        o << "V: returns volume" << std::endl;
        o << "c: returns centre of mass" << std::endl;
        o << "m: returns mass" << std::endl;
        o << "Ws: returns scalar weight" << std::endl;
        o << "W: returns weight (wrench)" << std::endl;
        o << "potential: returns potential energy" << std::endl;
        o << "b 'ref': returns block 'ref'" << std::endl;
        o << "connected: check whether all blocks are connected" << std::endl;
        o << "out: get list of extreme blocks (blocks connected only to 1 other block" << std::endl;
        o << "in: get list of internal blocks (blocks connected at least to 2 other blocks)" << std::endl;
        o << "nodes: get list of blocks connected to more than 2 other blocks (sources of hyperstaticity)" << std::endl;
        o << "clearinternal: clear all the internal face forces" << std::endl;
        o << "clearexternal: clear all the forces on faces" << std::endl;
        o << "clearall: clear all the forces on faces" << std::endl;
        o << "objective 'objective': set objective function (used for optimisation)" << std::endl;
        o << "objective: get value of objective function" << std::endl;
        o << "ok: check whether forces are acceptable" << std::endl;
        o << "penalty: penalty value used for optimization" << std::endl;
        o << "penalty_factor: get penalty factor" << std::endl;
        o << "objectivep: get value of objective function + penalty*penalty_factor" << std::endl;
        o << std::endl;
        o << "> Compute forces" << std::endl;
        o << "................" << std::endl;
        o << "compute 'b' 'f': compute as far as possible the forces on the blocks, starting from block b, entrance face f, return the last block for which it succeeded" << std::endl;
        o << "compute 'b1' 'f1' 'u1' 'v1' 'b2' 'f2' 'u2' 'v2' 'b3' 'f3' 'u3' 'v3': compute inside forces resulting in resistance line passing by 3 points. 'bi': block references, 'fi': faces references, 'ui': excentricity in the direction of u in local reference system of face (def depends of face), 'vj': excentricity in the direction of v in local reference system of face" << std::endl;
        o << "others... work in progress" << std::endl;
        return o.str();
}
//---------------------------------------------------------------------------
int a_structure::add_block(a_block * b) 
{
        blocks_.push_back(b); 
        return this->nb()-1;
}
//---------------------------------------------------------------------------
int a_structure::rb(const a_block* b) const
{
        for (int i=0; i<blocks_.size(); i++)
        {
                a_block * b2 = blocks_[i];
                if (b==b2)
                        return i;
        }
        return -1;
}
//---------------------------------------------------------------------------
double a_structure::V() const
{
        double v;
        for (int i=0; i<blocks_.size(); i++)
                v += blocks_[i]->V();
        return v;
}
//---------------------------------------------------------------------------
a_point a_structure::c() const
{
        double w = this->Ws();
        a_point c;
        for (int i=0; i<blocks_.size(); i++)
        {
                a_point ci = blocks_[i]->c();
                double wi = blocks_[i]->Ws();
                c += (ci*wi);
        }
        return c/w;
}
//---------------------------------------------------------------------------
double a_structure::m() const
{
        double m;
        for (int i=0; i<blocks_.size(); i++)
                m += blocks_[i]->m();
        return m;
}
//---------------------------------------------------------------------------
double a_structure::Ws() const
{
        double w;
        for (int i=0; i<blocks_.size(); i++)
                w += blocks_[i]->Ws();
        return w;
}
//---------------------------------------------------------------------------
a_wrench a_structure::W()
{
        a_wrench w;
        for (int i=0; i<blocks_.size(); i++)
                w += blocks_[i]->W();
        return w;
}
//---------------------------------------------------------------------------
a_wrench a_structure::fe()
{
        a_wrench w;
        for (int i=0; i<blocks_.size(); i++)
                w += blocks_[i]->fe();
        return w;
}
//---------------------------------------------------------------------------
a_wrench a_structure::ft()
{
        a_wrench w;
        for (int i=0; i<blocks_.size(); i++)
                w += blocks_[i]->ft();
        return w;
}
//---------------------------------------------------------------------------
double a_structure::potential() const
{
        double p;
        for (int i=0; i<blocks_.size(); i++)
                p += blocks_[i]->potential();
        return p;
}
//---------------------------------------------------------------------------
void a_structure::link_blocks(int rb1, int rf1, int rb2, int rf2)
{
        a_block * b1 = blocks_[rb1];
        a_block * b2 = blocks_[rb2];
        a_face * f1 = b1->i(rf1);
        a_face * f2 = b2->i(rf2);
        f1->lface(f2);
        f2->lface(f1);
        if (f1->contacttype()==-1)
                std::cerr << "warning! faces of block " << rb1 << " and " << rb2 << " are not touching." << std::endl;
}
//---------------------------------------------------------------------------
a_block * a_structure::b(const int i) const 
{
        if ((i<0)||(i>blocks_.size()-1)) throw range_error();
        return blocks_[i];
}
//---------------------------------------------------------------------------
bool a_structure::connected() const
{
        for (int i=0; i<blocks_.size(); i++)
        {
                if (blocks_[i]->out().size()==4)
                        return false;
        }
        return true;
}
//---------------------------------------------------------------------------
std::vector<int> a_structure::out() const
{
        std::vector<int> list;
        for (int i=0; i<blocks_.size(); i++)
        {
                if (blocks_[i]->in().size()==1)
                        list.push_back(i);
        }
        return list;
}
//---------------------------------------------------------------------------
std::vector<int> a_structure::in() const
{
        std::vector<int> list;
        for (int i=0; i<blocks_.size(); i++)
        {
                if (blocks_[i]->in().size()>1)
                        list.push_back(i);
        }
        return list;
}
//---------------------------------------------------------------------------
std::vector<int> a_structure::nodes() const
{
        std::vector<int> list;
        for (int i=0; i<blocks_.size(); i++)
        {
                if (blocks_[i]->in().size()>2)
                        list.push_back(i);
        }
        return list;
}
//---------------------------------------------------------------------------
void a_structure::clearinternal()
{
        for (int i=0; i<this->blocks_.size(); i++)
                blocks_[i]->clearinternal();
}
//---------------------------------------------------------------------------
void a_structure::clearexternal()
{
        for (int i=0; i<this->blocks_.size(); i++)
                blocks_[i]->clearexternal();
}
//---------------------------------------------------------------------------
void a_structure::clearall()
{
        for (int i=0; i<this->blocks_.size(); i++)
                blocks_[i]->clearall();
}
//---------------------------------------------------------------------------
double a_structure::objective() const
{
        if (objective_==0) return 0;
        return objective_->f(this);
}
//---------------------------------------------------------------------------
double a_structure::objectivep() const 
{
        return this->objective()+penalty_factor_*this->penalty();
}
//---------------------------------------------------------------------------
void a_structure::criteria(a_fcriteria * criteria)
{
        for (int i=0; i<this->blocks_.size(); i++)
                blocks_[i]->criteria(criteria);
}
//---------------------------------------------------------------------------
bool a_structure::ok() const
{
        for (int i=0; i<this->blocks_.size(); i++)
        {
                if (!blocks_[i]->ok())
                        return false;
        }
        return true;
}
//---------------------------------------------------------------------------
double a_structure::penalty() const
{
        double val = 0;
        for (int i=0; i<this->blocks_.size(); i++)
        {
                double v = blocks_[i]->penalty();
//VE(v)
                if (v>val)
                        val = v;
        }
        return val;
}
//---------------------------------------------------------------------------
void a_structure::rotate(int rb, int rf, const a_point& x, const a_point& dir, double r)
{
        if (!this->check_blockref(rb,rf)) return;
        a_block * b = blocks_[rb];
        b->entrance(rf);
        b->pos()->rotate(x,dir,r);
        while (a_block * b2 = b->nextblock())
        {
                b2->pos()->rotate(x,dir,r);
                b = b2;
        }
}
//---------------------------------------------------------------------------
void a_structure::rotate(int rb, int rf, double ex, double r)
{
        if (!this->check_blockref(rb,rf)) return;
        a_block * b = blocks_[rb];
        a_face * f = b->i(rf);
        //a_face_2d4 * f = cf_2d4(b->i(rf));
        a_point x = f->point(ex,0.);
        a_point dir = f->ny();
        this->rotate(rb,rf,x,dir,r);
}
//---------------------------------------------------------------------------
void a_structure::slide(int rb, int rf, double u, double v)
{
        if (!this->check_blockref(rb,rf)) return;
        a_block * b = blocks_[rb];
        b->entrance(rf);
        a_face_2d4 * f = cf_2d4(b->i(rf));
        a_point nx = f->nx();
        a_point ny = f->ny();
        a_point d = u*nx+v*ny;
        b->pos()->translate(d);
        while (a_block * b2 = b->nextblock())
        {
                b2->pos()->translate(d);
                b = b2;
        }
}
//---------------------------------------------------------------------------
void a_structure::material(a_material * material)
{
        for (int i=0; i<blocks_.size(); i++)
                blocks_[i]->material(material);
}
//---------------------------------------------------------------------------
void a_structure::storeexits(std::vector<int>& exits) const
{
        for (int i=0; i<this->nb(); i++)
                exits.push_back(blocks_[i]->exit());
}
//---------------------------------------------------------------------------
void a_structure::restoreexits(std::vector<int>& exits)
{
        for (int i=0; i<this->nb(); i++)
                blocks_[i]->exit(exits[i]);
        exits.clear();
}
//---------------------------------------------------------------------------
bool a_structure::faceonexit(const int rb1, const int rf1, const int rb2, const int rf2)
{
//VE(rb1)
//VE(rf1)
//VE(rb2)
//VE(rf2)
        bool exit;
        std::vector<int> exits;
                this->storeexits(exits);
        a_block * b1 = blocks_[rb1];
//VE(b1->c())
        b1->entrance(rf1);
        a_block * b2 = blocks_[rb2];
        while (true)
        {
                a_block * b = b1->nextblock();
//VE(b1->c())
                if (b==0)
                {
                        throw a_block::no_block_error();
                }
                if (b==b2)
                {
                        if (b->exit()==rf2)
                                exit = true;
                        else
                                exit = false;
                        break;
                }
                b1=b;
        }
        this->restoreexits(exits);
        return exit;
}
//---------------------------------------------------------------------------
a_block * a_structure::compute(int b, int f)
{
        std::vector<int> exits;
        this->storeexits(exits);
        a_block * block = blocks_[b];
        a_block * block0 = block;
        //set the exit face of first block
        int code = block->entrance(f);
        if (code != 0)
        {
                return 0;
        }
        while (a_block * block2 = block->compute())
        {
                block = block2;
        }
        this->restoreexits(exits);
        return block;
}
//---------------------------------------------------------------------------
bool a_structure::check_block(const int rb) const
{
        int nb = this->nb();
        if ((rb<0)||(rb>nb-1))
                return false;
        return true;
}
//---------------------------------------------------------------------------
bool a_structure::check_blockref(const int rb, const int rf) const
{
        int nb = this->nb();
        if ((rb<0)||(rb>nb-1))
                return false;
        if ((rf<0)||(rf>3))
                return false;
        return true;
}
//---------------------------------------------------------------------------
void a_structure::compute(const int rb1, const int rf1, const double u1, const double v1,
             const int rb2, const int rf2, const double u2, const double v2,
             const int rb3, const int rf3, const double u3, const double v3)
{
        if (!this->check_blockref(rb1,rf1)) return;
        if (!this->check_blockref(rb2,rf2)) return;
        if (!this->check_blockref(rb3,rf3)) return;
        std::vector<int> exits;
                this->storeexits(exits);
        a_block * b0 = blocks_.front();
        for (int i =0; i<b0->ni(); i++)
        {
                a_face * f = b0->i(i);
                if ((f->out())&&(f->exit()))
                        f->f(0.,0.,0.);
        }
        a_block * bn = blocks_.back();
        for (int i =0; i<bn->ni(); i++)
        {
                a_face * f = bn->i(i);
                if ((f->out())&&(f->exit()))
                        f->f(0.,0.,0.);
        }
        a_block * b1 = blocks_[rb1];
        b1->entrance(rf1);
        a_block * b2 = blocks_[rb2];
        a_block * bm;
        if (!this->faceonexit(rb1,rf1,rb2,rf2))
        {
                b2->exit(rf2);
                bm = b2->nextblock();
                a_face_2d4 * f2 = cf_2d4(b2->i(rf2));
                int r = f2->lface()->rf();
                bm->exit(r);
        }
        else
        {
                bm = b2;
                bm->exit(rf2);
        }
        a_block * bm2 = bm->nextblock();
        a_block * b3 = blocks_[rb3];
        a_wrench ft1 = b1->ft(*bm); 
        a_wrench ft2 = bm2->ft(*b3);
        a_wrench ftt = ft1+ft2;
        a_face_2d4 * f1 = cf_2d4(b1->i(rf1));
        a_face_2d4 * f2 = cf_2d4(b2->i(rf2));
        a_face_2d4 * f3 = cf_2d4(b3->i(rf3));
        a_point p1 = f1->point(u1,v1);
        a_point p2 = f2->point(u2,v2);
        a_point p3 = f3->point(u3,v3);
        //now use directly a_wrench capacity (not sure it is a good idea, see commit 64503e0a4b243fbc50fced6c738aefec455a8879 in case of problem)
        a_wrench f0(p1,p2,p3,ft1,ft2);
        f1->f(f0);
        this->compute(rb1,rf1);
        a_face * f1b = f1->lface();
        if (f1b!=0)
        {
                a_block * b1b = f1b->block();
                f1b->f(-f0);
                //get reference of face
                int rf1b = f1b->rf();
                //set the exit face block
                b1b->entrance(rf1b);
                int r1b = this->rb(b1b);
                this->compute(r1b,rf1b);
        }
        this->restoreexits(exits);
//VE(this->objective())
//VE(this->b(0)->i(3)->f())
}
//---------------------------------------------------------------------------
void a_structure::midblock(     const int rb1, const int rf1,
                                const int rb2, const int rf2,
                                a_block *& b, int&rf)
{
        if (!this->check_blockref(rb1,rf1)) return;
        if (!this->check_blockref(rb2,rf2)) return;
        a_block * b1 = blocks_[rb1];
        a_block * b0 = b1;
        b1->entrance(rf1);
        a_block * b2 = blocks_[rb2];
        int r = 0;
        while (true)
        {
                r++;
                a_block * b = b1->nextblock(); 
                if (b==0)
                        throw a_block::no_block_error();
                if (b==b2)
                        break;
                b1=b;
        }
        if (r<2)
                return;
        int rb = r/2;
        b1 = b0;
        for (int i=0; i<rb; i++)
        {
                a_block * b = b1->nextblock();
                b1 = b;
        }
        b = b1;
        rf = b1->exit();
}
//---------------------------------------------------------------------------
class vet : public vnl_least_squares_function
{
public:
        vet(a_structure * s,
                int rb1, int rf1,
                int rb2, int rf2,
                int rb3, int rf3):
                vnl_least_squares_function(3, 2, no_gradient), rb1_(rb1), rf1_(rf1), rb2_(rb2), rf2_(rf2), rb3_(rb3), rf3_(rf3),
                s_(s), pfactor_(1.) {}
        void f(const vnl_vector<double>& x, vnl_vector<double>& fx)
        {
                s_->compute(rb1_,rf1_,x[0],0.,rb2_,rf2_,x[1],0.,rb3_,rf3_,x[2],0.);
                double v = s_->objective();
                double p = s_->penalty();
                fx[0] = v;
                fx[1] = p*pfactor_;
//VL(x[0])
//VL(x[1])
//VL(x[2])
//VE(v)
        }
        void pfactor(const double val) {pfactor_=val;}
protected:
        a_structure * s_;
        int rb1_, rf1_, rb2_, rf2_, rb3_, rf3_;
        double pfactor_;
 
};
//---------------------------------------------------------------------------
class a_objective_function : public vnl_cost_function
{
public:
        a_objective_function( a_structure * s,
                                int rb1, int rf1,
                                int rb2, int rf2,
                                int rb3, int rf3) :
                                vnl_cost_function(3),
                                rb1_(rb1), rf1_(rf1), rb2_(rb2), rf2_(rf2), rb3_(rb3), rf3_(rf3), 
                                s_(s) {};
        double f(const vnl_vector<double>& x);
protected:
        a_structure * s_;
        int rb1_, rf1_, rb2_, rf2_, rb3_, rf3_;
};
//---------------------------------------------------------------------------
double a_objective_function::f(const vnl_vector<double>& x)
{
/*      double u0 = sin(x[0]);
        double u1 = sin(x[1]);
        double u2 = sin(x[2]);
*/
        double u0 = x[0];
        double u1 = x[1];
        double u2 = x[2];
        s_->compute(rb1_,rf1_,u0,0.,rb2_,rf2_,u1,0.,rb3_,rf3_,u2,0.);
        return  s_->objectivep();
}
/*
//---------------------------------------------------------------------------
void a_structure::min(  const int rb1, const int rf1,
                        const int rb3, const int rf3)
{
        this->optimise(true,rb1,rf1,rb3,rf3);
}
*/
//---------------------------------------------------------------------------
void a_structure::max(  const int rb1, const int rf1,
                        const int rb3, const int rf3)
{
        this->optimise(false,rb1,rf1,rb3,rf3);
}
//---------------------------------------------------------------------------
double a_structure::grad1(const int rb1, const int rf1,
                                const int rb2, const int rf2,
                                const int rb3, const int rf3,
                                const double * u)
{
        const double small=1.e-4;
        this->compute(rb1,rf1,u[0]+small,0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
        double v2 = this->objectivep();
        this->compute(rb1,rf1,u[0]-small,0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
        double v1 = this->objectivep();
        return (v2-v1)/small/2.;
}
//---------------------------------------------------------------------------
double a_structure::grad2(const int rb1, const int rf1,
                                const int rb2, const int rf2,
                                const int rb3, const int rf3,
                                const double * u)
{
        const double small=1.e-4;
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1]+small,0.,rb3,rf3,u[2],0.);
        double v2 = this->objectivep();
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1]-small,0.,rb3,rf3,u[2],0.);
        double v1 = this->objectivep();
        return (v2-v1)/small/2.;
}
//---------------------------------------------------------------------------
double a_structure::grad3(const int rb1, const int rf1,
                                const int rb2, const int rf2,
                                const int rb3, const int rf3,
                                const double * u)
{
        const double small=1.e-4;
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2]+small,0.);
        double v2 = this->objectivep();
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2]-small,0.);
        double v1 = this->objectivep();
        return (v2-v1)/small/2.;
}
//---------------------------------------------------------------------------
void a_structure::minsteep(     const int rb1, const int rf1,
                                const int rb3, const int rf3)
{
        if (objective_ == 0) return;
        std::vector<int> exits;
        this->storeexits(exits);
        int rf2;
        a_block * b2;
        this->midblock(rb1,rf1,rb3,rf3,b2,rf2);
        int rb2 = this->rb(b2);
        double u[] = {0.,0.,0.};
        double g1 = this->grad1(rb1,rf1,rb2,rf2,rb3,rf3,u);
        double g2 = this->grad2(rb1,rf1,rb2,rf2,rb3,rf3,u);
        double g3 = this->grad3(rb1,rf1,rb2,rf2,rb3,rf3,u);
        a_point g(g1,g2,g3);
        this->mindir(g,rb1,rf1,u[0],rb2,rf2,u[1],rb3,rf3,u[2]);
/*
//      g.max();
        VE(g)
        for (int i=0; i<20; i++)
        {
                this->mindir(g,rb1,rf1,u[0],rb2,rf2,u[1],rb3,rf3,u[2]);
VL(u[0])
VL(u[1])
VE(u[2])
                double g1 = this->grad1(rb1,rf1,rb2,rf2,rb3,rf3,u);
                double g2 = this->grad2(rb1,rf1,rb2,rf2,rb3,rf3,u);
                double g3 = this->grad3(rb1,rf1,rb2,rf2,rb3,rf3,u);
                g.set(g1,g2,g3);
        }
*/
        this->restoreexits(exits);
}
//---------------------------------------------------------------------------
void a_structure::optimise(     const bool is_min,
                                const int rb1, const int rf1,
                                const int rb3, const int rf3)
{
        if (objective_ == 0) return;
        std::vector<int> exits;
        this->storeexits(exits);
        double penalty_factor0 = penalty_factor_;
        const double mult = 10.;
        int rf2;
        a_block * b2;
        this->midblock(rb1,rf1,rb3,rf3,b2,rf2);
        int rb2 = this->rb(b2);
        double u[] = {0.,0.,0.};
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
        penalty_factor_ = penalty_factor_/mult;
        double v0 = this->objectivep();
        for (int i=0;i<5;i++)
        {
VE(this->b(0)->exit())
                penalty_factor_=penalty_factor_*mult;
VE(penalty_factor_)
this->minsteep(rb1,rf1,rb3,rf3);
                while (1)
                {
                        for (int f=0;f<3;f++)
                        {
                                if (is_min)
                                        u[f] = this->min(f,rb1,rf1,u[0],rb2,rf2,u[1],rb3,rf3,u[2]);
                                else
                                        u[f] = this->max(f,rb1,rf1,u[0],rb2,rf2,u[1],rb3,rf3,u[2]);
                        }
                        double v = this->objectivep();
                        if (v==v0)
                                break;
                        v0=v;
                }
VL(u[0])
VL(u[1])
VE(u[2])
        }
        this->restoreexits(exits);
        penalty_factor_ = penalty_factor0;
}
//---------------------------------------------------------------------------
void a_structure::min( const int rb1, const int rf1,
                        const int rb3, const int rf3)
{
        if (objective_ == 0) return;
        std::vector<int> exits;
        this->storeexits(exits);
        double penalty_factor0 = penalty_factor_;
        const double mult = 10.;
        int rf2;
        a_block * b2;
        this->midblock(rb1,rf1,rb3,rf3,b2,rf2);
        int rb2 = this->rb(b2);
 
        vnl_vector<double> x(3);
        x(0) = x(1) = x(2) = 0;
        a_objective_function f(this,rb1,rf1,rb2,rf2,rb3,rf3);
        vnl_lbfgs op(f);
//      op.minimize(x);
 
        this->restoreexits(exits);
        penalty_factor_ = penalty_factor0;
}
//---------------------------------------------------------------------------
double a_structure::min(       const int face,
                        const int rb1, const int rf1, const double u1,
                        const int rb2, const int rf2, const double u2,
                        const int rb3, const int rf3, const double u3)
{
        double u[] = {u1,u2,u3};
        const double gold = (1+sqrt(5.))/2.;
        double x1 = -1;
        double x3 = 1.;
        double x2 = x3-(x3-x1)/gold;
        u[face] = x2;
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
        double f2 = this->objectivep();
        bool grow = true;
        while (x3-x1>0.0001)
        {
                double xn;
                if (grow)
                        xn = x3-(x3-x2)/gold;
                else
                        xn = x1+(x2-x1)/gold;
                u[face] = xn;
                this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
                double fn = this->objectivep();
                if (fn<f2)
                {
                        if (grow)
                                x1 = x2;
                        else
                                x3 = x2;
                        x2 = xn;
                        f2 = fn;
                }
                else
                {
                        if (grow)
                                x3 = xn;
                        else
                                x1 = xn;
                        grow = !grow;
                }
        }
        VE(f2)
        return x2;
}
//---------------------------------------------------------------------------
double a_structure::mindir(const a_point & dir,
                        const int rb1, const int rf1, double & u1,
                        const int rb2, const int rf2, double & u2,
                        const int rb3, const int rf3, double & u3)
{
        a_point u(u1,u2,u3);
        double step = .1;
        double f;
        do 
        {
                VE(step)
                f = 1.e30;
                double f0;
                for (int i=0; i<100; i++)
                {
                        f0 = f;
                        u -= step*dir;
                        this->compute(rb1,rf1,u.x(),0.,rb2,rf2,u.y(),0.,rb3,rf3,u.z(),0.);
                        f = this->objectivep();
                        if (f>f0)
                        {
                                CL(**)
                                break;
                        }
                        VL(f)
                        VL(f0)
                }
                if (f<f0)
                        std::exit(-1);
                u += 2*step*dir;
                step /= 10.;
                CE(-)
        } while (step>.0001);
        u1 = u.x(); u2 = u.y(); u3 = u.z();
        return f;
}
//---------------------------------------------------------------------------
double a_structure::max(const int face,
                        const int rb1, const int rf1, const double u1,
                        const int rb2, const int rf2, const double u2,
                        const int rb3, const int rf3, const double u3)
{
        double u[] = {u1,u2,u3};
        const double gold = (1+sqrt(5.))/2.;
        double x1 = -1;
        double x3 = 1.;
        double x2 = x3-(x3-x1)/gold;
        u[face] = x2;
        this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
        double f2 = this->objectivep();
        bool grow = true;
        while (x3-x1>0.0001)
        {
                double xn;
                if (grow)
                        xn = x3-(x3-x2)/gold;
                else
                        xn = x1+(x2-x1)/gold;
                u[face] = xn;
                this->compute(rb1,rf1,u[0],0.,rb2,rf2,u[1],0.,rb3,rf3,u[2],0.);
                double fn = this->objectivep();
                if (fn>f2)
                {
                        if (grow)
                                x1 = x2;
                        else
                                x3 = x2;
                        x2 = xn;
                        f2 = fn;
                }
                else
                {
                        if (grow)
                                x3 = xn;
                        else
                                x1 = xn;
                        grow = !grow;
                }
        }
        VE(f2)
        return x2;
}
//***************************************************************************
//---------------------------------------------------------------------------
void a_structure::writetri(std::ostream &o) const
{
        a_trianglecloud tri;
        this->trianglecloud(tri);
        tri.write(o);
        tri.clear();
}
//---------------------------------------------------------------------------
void a_structure::writeb(std::ostream &o) const
{
        for (int i=0; i<blocks_.size(); i++)
        {
                blocks_[i]->writeb(o);
                o << std::endl;
        }
}
//---------------------------------------------------------------------------
void a_structure::writeg(std::ostream &o) const
{
        for (int i=0; i<blocks_.size(); i++)
        {
                o << "in b" << i << ".s ";
                blocks_[i]->writeg(o);
                o << std::endl;
                o << "r b" << i << ".r u b" << i << ".s" << std::endl;
        }
        o << "g structure.g ";
        for (int i=0; i<blocks_.size(); i++)
                o << "b" << i << ".r ";
        o << std::endl;
}
//---------------------------------------------------------------------------
void a_structure::trianglecloud(a_trianglecloud& tri) const
{
        tri.clear();
        for (int i=0; i<blocks_.size(); i++)
                blocks_[i]->trianglecloud(tri);
}
//---------------------------------------------------------------------------
void a_structure::polyline(a_linecloud& li) const
{
        li.clear();
        for (int i=0; i<blocks_.size(); i++)
        {
                a_block * block = blocks_[i];
                int k;
                k = 0;
                a_point p[2];
                for (int j=0; j<block->ni(); j++)
                {
                        a_face * fac = blocks_[i]->i(j);
                        //will work only if 2 inside faces -> arch, column...
                        //but polyline make no sense otherwise
                        if ((fac->exit()==1)||(fac->in()==1))
                        {
                                p[k] = fac->x();
                                k += 1;
                        }
                }
                if (k==2)
                        li.addline(p[0],p[1]);
                if (k>2)
                {
                        std::cout << "block " << i << " contains too many inside connections (" << k << ")!" << std::endl;
                        std::cout << "no line was added for this block" << std::endl;
                }
        }
}
//---------------------------------------------------------------------------
std::istream& operator>> (std::istream& in, a_structure& s)
{
        int nb;
        in >> nb;
        for (int i=0; i<nb; i++)
        {
                a_block * b = new a_block_2d4;
                s.blocks_.push_back(b);
                in >> (*b);
        }
        return in;
}
//---------------------------------------------------------------------------
std::ostream& operator<<(std::ostream& o, const a_structure& s)
{
        o << s.nb() << std::endl;
        for (int i=0; i<s.nb(); i++)
                o << s.blocks_[i] << std::endl;
        return o;
}