src/geometry/geometries1D.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    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 3 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/>.
 */

/*!
  \file geometries1D.cpp
  \authors N. Kielbasiewicz, Y. Lafranche
  \since 18 oct 2012
  \date 29 jul 2015

  \brief This file contains the implementation of methods of classes defined in geometries1D.hpp
*/

#include "geometries1D.hpp"

namespace xlifepp
{

//===========================================================
//
// Curve and child classes member functions
//
//===========================================================

//! default curve is inside the bounding box [0,1]
Curve::Curve()
: Geometry(BoundingBox(0.,1.), 1)
{
  // funh_=0;
}

void Curve::buildParam(const Parameter& p)
{
  trace_p->push("Curve::buildParam");
  ParameterKey key=p.key();

  switch (key)
  {
    case _pk_side_names:
    {
      if (p.type()==_string) { sideNames_.resize(1,p.get_s()); }
      else if (p.type()==_stringVector) { sideNames_=p.get_sv(); }
      else { error("param_badtype",words("value",p.type()),words("param key",key)); }
      break;
    }
    default: Geometry::buildParam(p); break;
  }
  trace_p->pop();
}

void Curve::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Curve::buildDefaultParam");
  switch (key)
  {
    case _pk_side_names: sideNames_.resize(0); break;
    default: Geometry::buildDefaultParam(key); break;
  }
  trace_p->pop();
}

std::set<ParameterKey> Curve::getParamsKeys()
{
  std::set<ParameterKey> params=Geometry::getParamsKeys();
  params.insert(_pk_side_names);
  return params;
}

//! default is segment [0,1]
Segment::Segment() : Curve(), p1_(0.), p2_(1.), n_(2)
{
  shape_=_segment;
  computeMB();
}

void Segment::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Segment::build");
  shape_=_segment;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // managing params
  for (number_t i=0; i < ps.size(); ++i)
  {
    ParameterKey key=ps[i].key();
    buildParam(ps[i]);
    if (params.find(key) != params.end()) { params.erase(key); }
    else
    {
      if (usedParams.find(key) == usedParams.end())
      { error("geom_unexpected_param_key", words("param key",key), words("shape",_segment)); }
      else { warning("param_already_used", words("param key",key)); }
    }
    usedParams.insert(key);
    // user must use nnodes or hsteps, not both
    if (key==_pk_hsteps && usedParams.find(_pk_nnodes) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_nnodes)); }
    if (key==_pk_nnodes && usedParams.find(_pk_hsteps) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_hsteps)); }
    // if user sets v1 or v2, and xmin or xmax, and vice versa it is an error
    if ((key==_pk_xmin || key==_pk_xmax) && usedParams.find(_pk_v1) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_v1)); }
    if ((key==_pk_xmin || key==_pk_xmax) && usedParams.find(_pk_v2) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_v2)); }
    if ((key==_pk_v1 || key==_pk_v2) && usedParams.find(_pk_xmin) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_xmin)); }
    if ((key==_pk_v1 || key==_pk_v2) && usedParams.find(_pk_xmax) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_xmax)); }
  }

  // if hsteps is not used, nnodes is used instead and there is no default value
  if (params.find(_pk_hsteps) != params.end()) { params.erase(_pk_hsteps); }

  // (v1,v2) or (xmin,xmax) has to be set (no default values)
  if (params.find(_pk_v1) != params.end() && params.find(_pk_v2) != params.end() &&
      params.find(_pk_xmin) != params.end() && params.find(_pk_xmax) != params.end())
  { error("param_missing","v1, v2, xmin, xmax"); }
  if (params.find(_pk_v1) != params.end() && params.find(_pk_v2) == params.end()) { error("param_missing","v1"); }
  if (params.find(_pk_v2) != params.end() && params.find(_pk_v1) == params.end()) { error("param_missing","v2"); }
  if (params.find(_pk_xmin) != params.end() && params.find(_pk_xmax) == params.end()) { error("param_missing","xmin"); }
  if (params.find(_pk_xmax) != params.end() && params.find(_pk_xmin) == params.end()) { error("param_missing","xmax"); }

  if (params.find(_pk_xmin) != params.end()) { params.erase(_pk_xmin); params.erase(_pk_xmax); }
  if (params.find(_pk_v1) != params.end()) { params.erase(_pk_v1); params.erase(_pk_v2); }
  std::set<ParameterKey>::const_iterator it_p;
  for (it_p=params.begin(); it_p != params.end(); ++it_p) { buildDefaultParam(*it_p); }
  boundingBox=BoundingBox(p1_,p2_);
  computeMB();
  trace_p->pop();
}

void Segment::buildParam(const Parameter& p)
{
  trace_p->push("Segment::buildParam");
  ParameterKey key=p.key();
  switch (key)
  {
    case _pk_v1:
    {
      switch (p.type())
      {
        case _pt: p1_=p.get_pt(); break;
        case _integer: p1_=Point(real_t(p.get_i())); break;
        case _real: p1_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_v2:
    {
      switch (p.type())
      {
        case _pt: p2_=p.get_pt(); break;
        case _integer: p2_=Point(real_t(p.get_i())); break;
        case _real: p2_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_xmin:
    {
      switch (p.type())
      {
        case _integer: p1_=Point(real_t(p.get_i())); break;
        case _real: p1_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_xmax:
    {
      switch (p.type())
      {
        case _integer: p2_=Point(real_t(p.get_i())); break;
        case _real: p2_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_nnodes:
    {
      switch (p.type())
      {
        case _integer:
        {
          number_t n=p.get_n();
          n_=n > 2 ? n : 2;
          break;
        }
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_hsteps:
    {
      switch (p.type())
      {
        case _integer: h_=std::vector<real_t>(2,real_t(p.get_n())); break;
        case _real: h_=std::vector<real_t>(2,p.get_r()); break;
        case _realVector:
        {
          h_=p.get_rv();
          if (h_.size() != 2) { error("bad_size"); }
          break;
        }
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    default: Curve::buildParam(p); break;
  }
  trace_p->pop();
}

void Segment::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Segment::buildDefaultParam");
  switch (key)
  {
    case _pk_nnodes: n_=2; break;
    default: Curve::buildDefaultParam(key); break;
  }
  trace_p->pop();
}

std::set<ParameterKey> Segment::getParamsKeys()
{
  std::set<ParameterKey> params=Curve::getParamsKeys();
  params.insert(_pk_v1);
  params.insert(_pk_v2);
  params.insert(_pk_xmin);
  params.insert(_pk_xmax);
  params.insert(_pk_nnodes);
  params.insert(_pk_hsteps);
  return params;
}

Segment::Segment(const Parameter& p1, const Parameter& p2) : Curve()
{
  std::vector<Parameter> ps(2);
  ps[0]=p1; ps[1]=p2;
  build(ps);
}

Segment::Segment(const Parameter& p1, const Parameter& p2, const Parameter& p3) : Curve()
{
  std::vector<Parameter> ps(3);
  ps[0]=p1; ps[1]=p2; ps[2]=p3;
  build(ps);
}

Segment::Segment(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4) : Curve()
{
  std::vector<Parameter> ps(4);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4;
  build(ps);
}

Segment::Segment(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4,
                 const Parameter& p5) : Curve()
{
  std::vector<Parameter> ps(5);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5;
  build(ps);
}

string_t Segment::asString() const
{
  string_t s("Segment (");
  s += p1_.toString() +", " + p2_.toString() + ")";
  return s;
}

Segment& Segment::reverse()
{
  std::swap(p1_,p2_);
  if (h_.size()==2) { std::swap(h_[0],h_[1]); }
  return *this;
}

Segment operator~(const Segment& s)
{
  Segment s2=s;
  s2.reverse();
  return s2;
}

std::vector<const Point*> Segment::boundNodes() const
{
  std::vector<const Point*> nodes(2);
  nodes[0]=&p1_; nodes[1]=&p2_;
  return nodes;
}

std::vector<Point*> Segment::nodes()
{
  std::vector<Point*> nodes(2);
  nodes[0]=&p1_; nodes[1]=&p2_;
  return nodes;
}

std::vector<const Point*> Segment::nodes() const
{
  std::vector<const Point*> nodes(2);
  nodes[0]=&p1_; nodes[1]=&p2_;
  return nodes;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Segment::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(1);
  curves[0]=std::make_pair(_segment,boundNodes());
  return curves;
}

//! default ellipse arc is quarter of circle of center (0,0,0) and radius 1
EllArc::EllArc() : Curve(), c_(0.,0.), a_(1.,0.), p1_(1.,0.), p2_(0.,1.), n_(2)
{
  shape_=_ellArc;
  computeMB();
}

void EllArc::build(const std::vector<Parameter>& ps)
{
  trace_p->push("EllArc::build");
  shape_=_ellArc;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // managing params
  for (number_t i=0; i < ps.size(); ++i)
  {
    ParameterKey key=ps[i].key();
    buildParam(ps[i]);
    if (params.find(key) != params.end()) { params.erase(key); }
    else
    {
      if (usedParams.find(key) == usedParams.end())
      { error("geom_unexpected_param_key", words("param key",key), words("shape",_ellArc)); }
      else { warning("param_already_used", words("param key",key)); }
    }
    usedParams.insert(key);
    // user must use nnodes or hsteps, not both
    if (key==_pk_hsteps && usedParams.find(_pk_nnodes) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_nnodes)); }
    if (key==_pk_nnodes && usedParams.find(_pk_hsteps) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_hsteps)); }
  }

  // if hsteps is not used, nnodes is used instead and there is no default value
  if (params.find(_pk_hsteps) != params.end()) { params.erase(_pk_hsteps); }

  // center, v1 and v2 have to be set
  // apogee is optional (default value is v1)
  if (params.find(_pk_center) != params.end()) { error("param_missing","center"); }
  if (params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if (params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }

  std::set<ParameterKey>::const_iterator it_p;
  for (it_p=params.begin(); it_p != params.end(); ++it_p) { buildDefaultParam(*it_p); }
  computeBB();
  computeMB();
  trace_p->pop();
}

void EllArc::buildParam(const Parameter& p)
{
  trace_p->push("EllArc::buildParam");
  ParameterKey key=p.key();
  switch (key)
  {
    case _pk_center:
    {
      switch (p.type())
      {
        case _pt: c_=p.get_pt(); break;
        case _integer: c_=Point(real_t(p.get_i())); break;
        case _real: c_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_apogee:
    {
      switch (p.type())
      {
        case _pt: a_=p.get_pt(); break;
        case _integer: a_=Point(real_t(p.get_i())); break;
        case _real: a_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_v1:
    {
      switch (p.type())
      {
        case _pt: p1_=p.get_pt(); break;
        case _integer: p1_=Point(real_t(p.get_i())); break;
        case _real: p1_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_v2:
    {
      switch (p.type())
      {
        case _pt: p2_=p.get_pt(); break;
        case _integer: p2_=Point(real_t(p.get_i())); break;
        case _real: p2_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_nnodes:
    {
      switch (p.type())
      {
        case _integer:
        {
          number_t n=p.get_n();
          n_=n > 2 ? n : 2;
          break;
        }
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_hsteps:
    {
      switch (p.type())
      {
        case _integer: h_=std::vector<real_t>(2,real_t(p.get_n())); break;
        case _real: h_=std::vector<real_t>(2,p.get_r()); break;
        case _realVector:
        {
          h_=p.get_rv();
          if (h_.size() != 2) { error("bad_size"); }
          break;
        }
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    default: Curve::buildParam(p); break;
  }
  trace_p->pop();
}

void EllArc::buildDefaultParam(ParameterKey key)
{
  trace_p->push("EllArc::buildDefaultParam");
  switch (key)
  {
    case _pk_nnodes: n_=2; break;
    case _pk_apogee: a_=p1_; break;
    default: Curve::buildDefaultParam(key); break;
  }
  trace_p->pop();
}

std::set<ParameterKey> EllArc::getParamsKeys()
{
  std::set<ParameterKey> params=Curve::getParamsKeys();
  params.insert(_pk_v1);
  params.insert(_pk_v2);
  params.insert(_pk_center);
  params.insert(_pk_apogee);
  params.insert(_pk_nnodes);
  params.insert(_pk_hsteps);
  return params;
}

EllArc::EllArc(const Parameter& p1, const Parameter& p2, const Parameter& p3) : Curve()
{
  std::vector<Parameter> ps(3);
  ps[0]=p1; ps[1]=p2; ps[2]=p3;
  build(ps);
}

EllArc::EllArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4) : Curve()
{
  std::vector<Parameter> ps(4);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4;
  build(ps);
}

EllArc::EllArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4,
               const Parameter& p5) : Curve()
{
  std::vector<Parameter> ps(5);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5;
  build(ps);
}

EllArc::EllArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4,
               const Parameter& p5, const Parameter& p6) : Curve()
{
  std::vector<Parameter> ps(6);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5; ps[5]=p6;
  build(ps);
}

EllArc::EllArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4,
               const Parameter& p5, const Parameter& p6, const Parameter& p7) : Curve()
{
  std::vector<Parameter> ps(7);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5; ps[5]=p6; ps[6]=p7;
  build(ps);
}

string_t EllArc::asString() const
{
  string_t s("Elliptic arc (bounds = {");
  s += p1_.toString() +", " + p2_.toString() + "}, center = " + c_.toString() + ", apogee = " + a_.toString() + ")";
  return s;
}

EllArc& EllArc::reverse()
{
  std::swap(p1_,p2_);
  if (h_.size()==2) { std::swap(h_[0],h_[1]); }
  return *this;
}

EllArc operator~(const EllArc& e)
{
  EllArc e2=e;
  e2.reverse();
  return e2;
}

void EllArc::computeMB()
{
  Point i=(p1_+p2_)/2.;
  real_t r=c_.distance(a_);
  Point p=c_+(r/c_.distance(i))*(i-c_);
  minimalBox= MinimalBox(p1_,p2_,p1_+p-i);
}

void EllArc::computeBB()
{
  Point cp1=p1_-c_;
  Point cp2=p2_-c_;
  real_t r=c_.distance(a_);
  real_t h;

  Point ph=projectionOnStraightLine(p2_, c_, p1_, h);
  Point pp1=ph;
  Point pp2=c_+(r/c_.distance(p1_))*p1_;
  Point pp3=p2_;
  if (dot(ph-c_,cp1) < 0.)
  {
    Point pt=c_+p2_-ph;
    pp3=(r/c_.distance(pt))*pt;
  }
  boundingBox=BoundingBox(pp1, pp2, pp3);
}

std::vector<const Point*> EllArc::boundNodes() const
{
  std::vector<const Point*> nodes(2);
  nodes[0]=&p1_; nodes[1]=&p2_;
  return nodes;
}

std::vector<Point*> EllArc::nodes()
{
  std::vector<Point*> nodes(4);
  nodes[0]=&c_; nodes[1]=&a_; nodes[2]=&p1_; nodes[3]=&p2_;
  return nodes;
}

std::vector<const Point*> EllArc::nodes() const
{
  std::vector<const Point*> nodes(4);
  nodes[0]=&c_; nodes[1]=&a_; nodes[2]=&p1_; nodes[3]=&p2_;
  return nodes;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > EllArc::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(1);
  curves[0]=std::make_pair(_ellArc,boundNodes());
  return curves;
}

//! default circular arc is quarter of circle of center (0,0,0) and radius 1
CircArc::CircArc() : Curve(), c_(0.,0.), p1_(1.,0.), p2_(0.,1.), n_(2)
{
  shape_=_circArc;
  computeMB();
}

void CircArc::build(const std::vector<Parameter>& ps)
{
  trace_p->push("CircArc::build");
  shape_=_circArc;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // managing params
  for (number_t i=0; i < ps.size(); ++i)
  {
    ParameterKey key=ps[i].key();
    buildParam(ps[i]);
    if (params.find(key) != params.end()) { params.erase(key); }
    else
    {
      if (usedParams.find(key) == usedParams.end())
      { error("geom_unexpected_param_key", words("param key",key), words("shape",_circArc)); }
      else { warning("param_already_used", words("param key",key)); }
    }
    usedParams.insert(key);
    // user must use nnodes or hsteps, not both
    if (key==_pk_hsteps && usedParams.find(_pk_nnodes) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_nnodes)); }
    if (key==_pk_nnodes && usedParams.find(_pk_hsteps) != usedParams.end())
    { error("param_conflict", words("param key",key), words("param key",_pk_hsteps)); }
  }

  // if hsteps is not used, nnodes is used instead and there is no default value
  if (params.find(_pk_hsteps) != params.end()) { params.erase(_pk_hsteps); }

  // center, v1 and v2 have to be set
  if (params.find(_pk_center) != params.end()) { error("param_missing","center"); }
  if (params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if (params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }

  std::set<ParameterKey>::const_iterator it_p;
  for (it_p=params.begin(); it_p != params.end(); ++it_p) { buildDefaultParam(*it_p); }
  computeBB();
  computeMB();
  trace_p->pop();
}

void CircArc::buildParam(const Parameter& p)
{
  trace_p->push("CircArc::buildParam");
  ParameterKey key=p.key();
  switch (key)
  {
    case _pk_center:
    {
      switch (p.type())
      {
        case _pt: c_=p.get_pt(); break;
        case _integer: c_=Point(real_t(p.get_i())); break;
        case _real: c_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_v1:
    {
      switch (p.type())
      {
        case _pt: p1_=p.get_pt(); break;
        case _integer: p1_=Point(real_t(p.get_i())); break;
        case _real: p1_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_v2:
    {
      switch (p.type())
      {
        case _pt: p2_=p.get_pt(); break;
        case _integer: p2_=Point(real_t(p.get_i())); break;
        case _real: p2_=Point(p.get_r()); break;
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_nnodes:
    {
      switch (p.type())
      {
        case _integer:
        {
          number_t n=p.get_n();
          n_=n > 2 ? n : 2;
          break;
        }
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    case _pk_hsteps:
    {
      switch (p.type())
      {
        case _integer: h_=std::vector<real_t>(2,real_t(p.get_n())); break;
        case _real: h_=std::vector<real_t>(2,p.get_r()); break;
        case _realVector:
        {
          h_=p.get_rv();
          if (h_.size() != 2) { error("bad_size"); }
          break;
        }
        default: error("param_badtype",words("value",p.type()),words("param key",key));
      }
      break;
    }
    default: Curve::buildParam(p); break;
  }
  trace_p->pop();
}

void CircArc::buildDefaultParam(ParameterKey key)
{
  trace_p->push("CircArc::buildDefaultParam");
  switch (key)
  {
    case _pk_nnodes: n_=2; break;
    default: Curve::buildDefaultParam(key); break;
  }
  trace_p->pop();
}

std::set<ParameterKey> CircArc::getParamsKeys()
{
  std::set<ParameterKey> params=Curve::getParamsKeys();
  params.insert(_pk_v1);
  params.insert(_pk_v2);
  params.insert(_pk_center);
  params.insert(_pk_nnodes);
  params.insert(_pk_hsteps);
  return params;
}

CircArc::CircArc(const Parameter& p1, const Parameter& p2, const Parameter& p3) : Curve()
{
  std::vector<Parameter> ps(3);
  ps[0]=p1; ps[1]=p2; ps[2]=p3;
  build(ps);
}

CircArc::CircArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4) : Curve()
{
  std::vector<Parameter> ps(4);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4;
  build(ps);
}

CircArc::CircArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4,
                 const Parameter& p5) : Curve()
{
  std::vector<Parameter> ps(5);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5;
  build(ps);
}

CircArc::CircArc(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4,
                 const Parameter& p5, const Parameter& p6) : Curve()
{
  std::vector<Parameter> ps(6);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5; ps[5]=p6;
  build(ps);
}

void CircArc::computeMB()
{
  Point i=(p1_+p2_)/2.;
  real_t r=c_.distance(p1_);
  Point p=c_+(r/c_.distance(i))*(i-c_);
  minimalBox= MinimalBox(p1_,p2_,p1_+p-i);
}

void CircArc::computeBB()
{
  Point cp1=p1_-c_;
  Point cp2=p2_-c_;
  real_t r=c_.distance(p1_);
  real_t h;

  Point ph=projectionOnStraightLine(p2_, c_, p1_, h);
  Point pp1=ph;
  Point pp2=c_+(r/c_.distance(p1_))*p1_;
  Point pp3=p2_;
  if (dot(ph-c_,cp1) < 0.)
  {
    Point pt=c_+p2_-ph;
    pp3=(r/c_.distance(pt))*pt;
  }
  boundingBox=BoundingBox(pp1, pp2, pp3);
}

string_t CircArc::asString() const
{
  string_t s("Circular arc (bounds = {");
  s += p1_.toString() +", " + p2_.toString() + "}, center = " + c_.toString() + ")";
  return s;
}

CircArc& CircArc::reverse()
{
  std::swap(p1_,p2_);
  if (h_.size()==2) { std::swap(h_[0],h_[1]); }
  return *this;
}

CircArc operator~(const CircArc& c)
{
  CircArc c2=c;
  c2.reverse();
  return c2;
}

std::vector<const Point*> CircArc::boundNodes() const
{
  std::vector<const Point*> nodes(2);
  nodes[0]=&p1_; nodes[1]=&p2_;
  return nodes;
}

std::vector<Point*> CircArc::nodes()
{
  std::vector<Point*> nodes(3);
  nodes[0]=&c_; nodes[1]=&p1_; nodes[2]=&p2_;
  return nodes;
}

std::vector<const Point*> CircArc::nodes() const
{
  std::vector<const Point*> nodes(3);
  nodes[0]=&c_; nodes[1]=&p1_; nodes[2]=&p2_;
  return nodes;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > CircArc::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(1);
  curves[0]=std::make_pair(_circArc,boundNodes());
  return curves;
}

real_t CircArc::measure() const
{
  real_t theta;
  real_t radius=c_.distance(p1_);
  real_t scal=dot(p1_-c_,p2_-c_);
  if (std::abs(scal) < theEpsilon) { return 0.5*pi_*radius; }
  theta=std::atan(norm2(crossProduct(p1_-c_,p2_-c_))/scal);
  if (theta < 0) { theta+=pi_; }
  return theta;
}

string_t SetOfPoints::asString() const
{
  string_t s("SetOfPoints (");
  s += tostring(pts_.size()) + " points)";
  return s;
}

std::vector<Point*> SetOfPoints::nodes()
{
  std::vector<Point*> nodes(pts_.size());
  for (size_t i=0; i < pts_.size(); i++) { nodes[i]=&pts_[i]; }
  return nodes;
}

std::vector<const Point*> SetOfPoints::nodes() const
{
  std::vector<const Point*> nodes(pts_.size());
  for (size_t i=0; i < pts_.size(); i++) { nodes[i]=&pts_[i]; }
  return nodes;
}

real_t SetOfPoints::measure() const
{
  real_t length=0;
  for (number_t i=1; i < pts_.size(); ++i)
  {
    length += pts_[i-1].distance(pts_[i]);
  }
  return length;
}

} // end of namespace xlifepp