xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/geometry/geometries3D.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 geometries3D.cpp
  \authors N. Kielbasiewicz, Y. Lafranche
  \since 18 oct 2012
  \date 30 jul 2015

  \brief This file contains the implementation of methods of Volume class and Polyhedron and Ellipsoid classes
  defined in geometries3D.hpp
*/

#include "geometries2D.hpp"
#include "geometries3D.hpp"
#include "geometries_utils.hpp"

namespace xlifepp
{

//===========================================================
//
// Volume and child classes member functions
//
//===========================================================

//! default volume is inside bounding box [0,1]^3
Volume::Volume()
  : Geometry(BoundingBox(0.,1.,0.,1.,0.,1.), 3) {}

void Volume::buildParam(const Parameter& p)
{
  trace_p->push("Volume::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 Volume::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Volume::buildDefaultParam");
  switch(key)
    {
      case _pk_side_names: sideNames_.resize(0); break;
      default: Geometry::buildDefaultParam(key); break;
    }
  trace_p->pop();
}

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

//! default polyhedron is tetrahedron (0,0,0) (1,0,0) (0,1,0) (0,0,1)
Polyhedron::Polyhedron() : Volume()
{
  p_.resize(4);
  p_[0] = Point(0.,0.,0.);
  p_[1] = Point(1.,0.,0.);
  p_[2] = Point(0.,1.,0.);
  p_[3] = Point(0.,0.,1.);
  faces_.resize(4);
  faces_[0] = new Triangle(p_[0], p_[1], p_[2]);
  faces_[1] = new Triangle(p_[0], p_[1], p_[3]);
  faces_[2] = new Triangle(p_[1], p_[2], p_[3]);
  faces_[3] = new Triangle(p_[2], p_[0], p_[3]);
  shape_=_polyhedron;
  computeMB();
}

void Polyhedron::buildP()
{
  for(number_t i=0; i < faces_.size(); ++i)
    {
      boundingBox+=faces_[i]->boundingBox;
      for(number_t j=0; j < faces_[i]->p().size(); ++j)
        {
          int_t foundId=-1;
          for(number_t k=0; k < p_.size(); ++k)
            {
              if(p_[k] == faces_[i]->p(j+1)) { foundId=j; }
            }
          if(foundId == -1) { p_.push_back(faces_[i]->p(j+1)); }
        }
    }
}

void Polyhedron::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Polyhedron::build");
  shape_=_polyhedron;
  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",_polyhedron)); }
          else { warning("param_already_used", words("param key",key)); }
        }
      usedParams.insert(key);
    }

  // faces has to be set (no default values)
  if(params.find(_pk_faces) != params.end()) { error("param_missing","faces"); }

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

  buildP();

  computeMB();
  trace_p->pop();
}

void Polyhedron::buildParam(const Parameter& p)
{
  trace_p->push("Polyhedron::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_faces:
        {
          switch(p.type())
            {
              case _pointer:
                {
                  const std::vector<Polygon>& faces = *reinterpret_cast<const std::vector<Polygon>*>(p.get_p());
                  faces_.resize(faces.size());
                  for(number_t i=0; i < faces_.size(); ++i) { faces_[i]=faces[i].clonePG(); }
                  break;
                }
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Volume::buildParam(p); break;
    }
  trace_p->pop();
}

void Polyhedron::buildDefaultParam(ParameterKey key)
{
  Volume::buildDefaultParam(key);
}

std::set<ParameterKey> Polyhedron::getParamsKeys()
{
  std::set<ParameterKey> params=Volume::getParamsKeys();
  params.insert(_pk_faces);
  return params;
}

Polyhedron::Polyhedron(const Parameter& p1) : Volume()
{
  std::vector<Parameter> ps(1,p1);
  build(ps);
}

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

Polyhedron::Polyhedron(const Polyhedron& ph) : Volume(ph)
{
  p_=ph.p_;
  n_=ph.n_;
  h_=ph.h_;
  faces_.resize(ph.faces_.size());
  for(number_t i=0; i < faces_.size(); ++i)
    {
      faces_[i]=ph.faces_[i]->clonePG();
    }
}

string_t Polyhedron::asString() const
{
  string_t s("Polyhedron (");
  s += tostring(faces_.size()) + " faces )";
  return s;
}

std::vector<const Point*> Polyhedron::boundNodes() const
{
  std::vector<const Point*> nodes(p_.size());
  for(number_t i=0; i < p_.size(); ++i) { nodes[i]=&p_[i]; }
  return nodes;
}

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

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

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Polyhedron::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves, curves2;

  for(number_t i=0; i < faces_.size(); ++i)
    {
      for(number_t j=0; j < faces_[i]->curves().size(); ++j)
        {
          curves2=faces_[i]->curves();
          if(findBorder(curves2[j], curves) == -1)
            {
              curves.push_back(curves2[j]);
            }
        }
    }
  return curves;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Polyhedron::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(faces_.size());
  for(number_t i=0; i < faces_.size(); ++i)
    {
      surfs[i]=faces_[i]->surfs()[0];
    }
  return surfs;
}

//! collect in a list all canonical geometry's with name n
void Polyhedron::collect(const string_t& n, std::list<Geometry*>& geoms) const
{
  if(domName_==n) geoms.push_back(const_cast<Geometry*>(static_cast<const Geometry*>(this)));
  //loop on sides
  std::vector<Polygon*>::const_iterator its=faces_.begin();
  for(; its!=faces_.end(); ++its)
    {
      if((*its)->domName()==n) geoms.push_back((*its)->clone()); //create geometry related to side by cloning
    }
}

//! default tetrahedron is tetrahedron (0,0,0) (1,0,0) (0,1,0) (0,0,1)
Tetrahedron::Tetrahedron() : Polyhedron()
{
  n_.resize(6,2);
  shape_=_tetrahedron;
  computeMB();
  setFaces();
}

void Tetrahedron::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Tetrahedron::build");
  shape_=_tetrahedron;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // _faces key is replaced by _v1, _v2 and _v3
  params.erase(_pk_faces);

  p_.resize(4);
  // 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",_tetrahedron)); }
          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); }

  // (v1,v2,v3,v4) has to be set (no default values)
  if(params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v3) != params.end()) { error("param_missing","v3"); }
  if(params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }

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

  // checking nnodes and hsteps size
  if(n_.size() != 0)
    {
      if(n_.size() == 1) { number_t n0=n_[0]; n_.resize(6,n0); }
      else if(n_.size() != 6) {error("bad_size","nnodes",6,n_.size()); }
    }
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(4,h0); }
      else if(h_.size() != 4) {error("bad_size","hsteps",4,h_.size()); }
    }

  boundingBox=BoundingBox(p_[0], p_[1], p_[2], p_[3]);
  computeMB();
  setFaces();
  trace_p->pop();
}

void Tetrahedron::setFaces()
{
  faces_.resize(4);
  std::vector<number_t> ns(4,2);  // ns is set to 2
  if(sideNames_.size()>=4)
    {
      faces_[0] = new Triangle(p_[0], p_[1], p_[2], ns, sideNames_[0]);
      faces_[1] = new Triangle(p_[0], p_[1], p_[3], ns, sideNames_[1]);
      faces_[2] = new Triangle(p_[1], p_[2], p_[3], ns, sideNames_[2]);
      faces_[3] = new Triangle(p_[2], p_[0], p_[3], ns, sideNames_[3]);
    }
  else
    {
      string_t sn="";
      if(sideNames_.size()>=1) sn=sideNames_[0];
      faces_[0] = new Triangle(p_[0], p_[1], p_[2], ns, sn);
      faces_[1] = new Triangle(p_[0], p_[1], p_[3], ns, sn);
      faces_[2] = new Triangle(p_[1], p_[2], p_[3], ns, sn);
      faces_[3] = new Triangle(p_[2], p_[0], p_[3], ns, sn);
    }
}

real_t Tetrahedron::measure() const
{
  real_t h;
  Point P=projectionOnStraightLine(p_[2], p_[0], p_[1], h);
  real_t area=0.5*(p_[0].distance(p_[1])) * h;
  P=projectionOnTriangle(p_[3], p_[0], p_[1], p_[2], h);
  return area*h/3.;
}

void Tetrahedron::buildParam(const Parameter& p)
{
  trace_p->push("Tetrahedron::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_v1:
        {
          switch(p.type())
            {
              case _pt: p_[0]=p.get_pt(); break;
              case _integer: p_[0]=Point(real_t(p.get_i())); break;
              case _real: p_[0]=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: p_[1]=p.get_pt(); break;
              case _integer: p_[1]=Point(real_t(p.get_i())); break;
              case _real: p_[1]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v3:
        {
          switch(p.type())
            {
              case _pt: p_[2]=p.get_pt(); break;
              case _integer: p_[2]=Point(real_t(p.get_i())); break;
              case _real: p_[2]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v4:
        {
          switch(p.type())
            {
              case _pt: p_[3]=p.get_pt(); break;
              case _integer: p_[3]=Point(real_t(p.get_i())); break;
              case _real: p_[3]=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 _integerVector:
                {
                  std::vector<number_t> n=p.get_nv();
                  n_.resize(n.size());
                  for(number_t i=0; i < n.size(); ++i) { n_[i] = n[i] > 2 ? n[i] : 2; }
                  break;
                }
              case _integer: n_=std::vector<number_t>(1,p.get_n() > 2 ? p.get_n() : 2); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_hsteps:
        {
          switch(p.type())
            {
              case _realVector: h_=p.get_rv(); break;
              case _integer: h_=std::vector<real_t>(1,real_t(p.get_i())); break;
              case _real: h_=std::vector<real_t>(1,p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Polyhedron::buildParam(p); break;
    }
  trace_p->pop();
}

void Tetrahedron::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Tetrahedron::buildDefaultParam");
  switch(key)
    {
      case _pk_nnodes: n_=std::vector<number_t>(6,2); break;
      default: Polyhedron::buildDefaultParam(key); break;
    }
  trace_p->pop();
}

std::set<ParameterKey> Tetrahedron::getParamsKeys()
{
  std::set<ParameterKey> params=Polyhedron::getParamsKeys();
  params.insert(_pk_v1);
  params.insert(_pk_v2);
  params.insert(_pk_v3);
  params.insert(_pk_v4);
  params.insert(_pk_nnodes);
  params.insert(_pk_hsteps);
  return params;
}

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

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

Tetrahedron::Tetrahedron(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                         const Parameter& p6) : Polyhedron()
{
  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);
}

Tetrahedron::Tetrahedron(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                         const Parameter& p6, const Parameter& p7) : Polyhedron()
{
  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 Tetrahedron::asString() const
{
  string_t s("Tetrahedron (");
  s += p_[0].toString() +", " + p_[1].toString() + ", " + p_[2].toString() + ", " + p_[3].toString() + ")";
  return s;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Tetrahedron::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(6);
  std::vector<const Point*> vertices(2);
  vertices[0]=&p_[0]; vertices[1]=&p_[1];
  curves[0]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2];
  curves[1]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[0];
  curves[2]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[3];
  curves[3]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[3];
  curves[4]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[3];
  curves[5]=std::make_pair(_segment, vertices);

  return curves;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Tetrahedron::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(4);
  std::vector<const Point*> vertices(3);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[2];
  surfs[0]=std::make_pair(_triangle, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[3];
  surfs[1]=std::make_pair(_triangle, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2]; vertices[2]=&p_[3];
  surfs[2]=std::make_pair(_triangle, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[0]; vertices[2]=&p_[3];
  surfs[3]=std::make_pair(_triangle, vertices);

  return surfs;
}

Hexahedron::Hexahedron() : Polyhedron()
{
  p_.resize(8);
  p_[0]=Point(0.,0.,0.);
  p_[1]=Point(1.,0.,0.);
  p_[2]=Point(1.,1.,0.);
  p_[3]=Point(0.,1.,0.);
  p_[4]=Point(0.,0.,1.);
  p_[5]=Point(1.,0.,1.);
  p_[6]=Point(1.,1.,1.);
  p_[7]=Point(0.,1.,1.);
  n_.resize(12,2);
  shape_=_hexahedron;
  computeMB();
  setFaces();
}

void Hexahedron::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Hexahedron::build");
  shape_=_hexahedron;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // _faces key is replaced by _v1, _v2, _v3, _v4, _v5, _v6, _v7 and _v8
  params.erase(_pk_faces);

  p_.resize(8);
  // 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",_hexahedron)); }
          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); }

  // (v1,v2,v3,v4) has to be set (no default values)
  if(params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v3) != params.end()) { error("param_missing","v3"); }
  if(params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_v5) != params.end()) { error("param_missing","v5"); }
  if(params.find(_pk_v6) != params.end()) { error("param_missing","v6"); }
  if(params.find(_pk_v7) != params.end()) { error("param_missing","v7"); }
  if(params.find(_pk_v8) != params.end()) { error("param_missing","v8"); }

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

  // checking nnodes and hsteps size
  if(n_.size() != 0)
    {
      if(n_.size() == 1) { number_t n0=n_[0]; n_.resize(12,n0); }
      else if(n_.size() != 12) {error("bad_size","nnodes",12,n_.size()); }
    }
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(8,h0); }
      else if(h_.size() != 8) {error("bad_size","hsteps",8,h_.size()); }
    }

  boundingBox=BoundingBox(p_[0], p_[1], p_[4], p_[4]);
  computeMB();
  setFaces();
  trace_p->pop();
}

void Hexahedron::buildParam(const Parameter& p)
{
  trace_p->push("Hexahedron::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_v1:
        {
          switch(p.type())
            {
              case _pt: p_[0]=p.get_pt(); break;
              case _integer: p_[0]=Point(real_t(p.get_i())); break;
              case _real: p_[0]=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: p_[1]=p.get_pt(); break;
              case _integer: p_[1]=Point(real_t(p.get_i())); break;
              case _real: p_[1]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v3:
        {
          switch(p.type())
            {
              case _pt: p_[2]=p.get_pt(); break;
              case _integer: p_[2]=Point(real_t(p.get_i())); break;
              case _real: p_[2]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v4:
        {
          switch(p.type())
            {
              case _pt: p_[3]=p.get_pt(); break;
              case _integer: p_[3]=Point(real_t(p.get_i())); break;
              case _real: p_[3]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v5:
        {
          switch(p.type())
            {
              case _pt: p_[4]=p.get_pt(); break;
              case _integer: p_[4]=Point(real_t(p.get_i())); break;
              case _real: p_[4]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v6:
        {
          switch(p.type())
            {
              case _pt: p_[5]=p.get_pt(); break;
              case _integer: p_[5]=Point(real_t(p.get_i())); break;
              case _real: p_[5]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v7:
        {
          switch(p.type())
            {
              case _pt: p_[6]=p.get_pt(); break;
              case _integer: p_[6]=Point(real_t(p.get_i())); break;
              case _real: p_[6]=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_v8:
        {
          switch(p.type())
            {
              case _pt: p_[7]=p.get_pt(); break;
              case _integer: p_[7]=Point(real_t(p.get_i())); break;
              case _real: p_[7]=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 _integerVector:
                {
                  std::vector<number_t> n=p.get_nv();
                  n_.resize(n.size());
                  for(number_t i=0; i < n.size(); ++i) { n_[i] = n[i] > 2 ? n[i] : 2; }
                  break;
                }
              case _integer: n_=std::vector<number_t>(1,p.get_n()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_hsteps:
        {
          switch(p.type())
            {
              case _realVector: h_=p.get_rv(); break;
              case _integer: h_=std::vector<real_t>(1,real_t(p.get_i())); break;
              case _real: h_=std::vector<real_t>(1,p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Polyhedron::buildParam(p); break;
    }
  trace_p->pop();
}

void Hexahedron::setFaces()
{
  faces_.resize(6);
  std::vector<number_t> ns(4,2);  // ns is set to 2
  if(sideNames_.size()>=6)
    {
      faces_[0] = new Quadrangle(p_[0], p_[1], p_[2], p_[4], ns, sideNames_[0]);
      faces_[1] = new Quadrangle(p_[4], p_[5], p_[6], p_[7], ns, sideNames_[1]);
      faces_[2] = new Quadrangle(p_[0], p_[1], p_[5], p_[4], ns, sideNames_[2]);
      faces_[3] = new Quadrangle(p_[3], p_[2], p_[6], p_[7], ns, sideNames_[3]);
      faces_[4] = new Quadrangle(p_[0], p_[3], p_[7], p_[4], ns, sideNames_[4]);
      faces_[5] = new Quadrangle(p_[1], p_[2], p_[6], p_[5], ns, sideNames_[5]);
    }
  else
    {
      string_t sn="";
      if(sideNames_.size()>0) sn=sideNames_[0];
      faces_[0] = new Quadrangle(p_[0], p_[1], p_[2], p_[4], ns, sn);
      faces_[1] = new Quadrangle(p_[4], p_[5], p_[6], p_[7], ns, sn);
      faces_[2] = new Quadrangle(p_[0], p_[1], p_[5], p_[4], ns, sn);
      faces_[3] = new Quadrangle(p_[3], p_[2], p_[6], p_[7], ns, sn);
      faces_[4] = new Quadrangle(p_[0], p_[3], p_[7], p_[4], ns, sn);
      faces_[5] = new Quadrangle(p_[1], p_[2], p_[6], p_[5], ns, sn);
    }
}

void Hexahedron::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Hexahedron::buildDefaultParam");
  switch(key)
    {
      case _pk_nnodes: n_=std::vector<number_t>(12,2); break;
      default: Polyhedron::buildDefaultParam(key); break;
    }
  trace_p->pop();
}

std::set<ParameterKey> Hexahedron::getParamsKeys()
{
  std::set<ParameterKey> params=Polyhedron::getParamsKeys();
  params.insert(_pk_v1);
  params.insert(_pk_v2);
  params.insert(_pk_v3);
  params.insert(_pk_v4);
  params.insert(_pk_v5);
  params.insert(_pk_v6);
  params.insert(_pk_v7);
  params.insert(_pk_v8);
  params.insert(_pk_nnodes);
  params.insert(_pk_hsteps);
  return params;
}

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

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

Hexahedron::Hexahedron(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                       const Parameter& p6, const Parameter& p7, const Parameter& p8, const Parameter& p9, const Parameter& p10)
  : Polyhedron()
{
  std::vector<Parameter> ps(10);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5; ps[5]=p6; ps[6]=p7; ps[7]=p8; ps[8]=p9; ps[9]=p10;
  build(ps);
}

Hexahedron::Hexahedron(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                       const Parameter& p6, const Parameter& p7, const Parameter& p8, const Parameter& p9, const Parameter& p10,
                       const Parameter& p11) : Polyhedron()
{
  std::vector<Parameter> ps(11);
  ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4; ps[4]=p5; ps[5]=p6; ps[6]=p7; ps[7]=p8; ps[8]=p9; ps[9]=p10; ps[10]=p11;
  build(ps);
}

string_t Hexahedron::asString() const
{
  string_t s("Hexahedron (");
  s += p_[0].toString() +", " + p_[1].toString() + ", " + p_[2].toString() + ", " + p_[3].toString();
  s += p_[4].toString() +", " + p_[5].toString() + ", " + p_[6].toString() + ", " + p_[7].toString() + ")";
  return s;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Hexahedron::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(12);
  std::vector<const Point*> vertices(2);
  vertices[0]=&p_[0]; vertices[1]=&p_[1];
  curves[0]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2];
  curves[1]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[3];
  curves[2]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[3]; vertices[1]=&p_[0];
  curves[3]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[4]; vertices[1]=&p_[5];
  curves[4]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[5]; vertices[1]=&p_[6];
  curves[5]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[6]; vertices[1]=&p_[7];
  curves[6]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[7]; vertices[1]=&p_[4];
  curves[7]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[4];
  curves[8]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[5];
  curves[9]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[6];
  curves[10]=std::make_pair(_segment, vertices);
  vertices[0]=&p_[3]; vertices[1]=&p_[7];
  curves[11]=std::make_pair(_segment, vertices);

  return curves;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Hexahedron::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(6);
  std::vector<const Point*> vertices(4);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[2]; vertices[3]=&p_[3];
  surfs[0]=std::make_pair(_quadrangle, vertices);
  vertices[0]=&p_[4]; vertices[1]=&p_[5]; vertices[2]=&p_[6]; vertices[3]=&p_[7];
  surfs[1]=std::make_pair(_quadrangle, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[5]; vertices[3]=&p_[4];
  surfs[2]=std::make_pair(_quadrangle, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[3]; vertices[2]=&p_[7]; vertices[3]=&p_[6];
  surfs[3]=std::make_pair(_quadrangle, vertices);
  vertices[0]=&p_[3]; vertices[1]=&p_[0]; vertices[2]=&p_[4]; vertices[3]=&p_[7];
  surfs[4]=std::make_pair(_quadrangle, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2]; vertices[2]=&p_[6]; vertices[3]=&p_[5];
  surfs[5]=std::make_pair(_quadrangle, vertices);

  return surfs;
}

Parallelepiped::Parallelepiped() : Hexahedron()
{
  shape_=_parallelepiped;
  computeMB();
}

void Parallelepiped::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Parallelepiped::build");
  shape_=_parallelepiped;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // _faces key is replaced by _v1, _v2, _v4 and _v5
  params.erase(_pk_faces);
  // _v3, _v6, _v7, _v8 are disabled
  params.erase(_pk_v3);
  params.erase(_pk_v6);
  params.erase(_pk_v7);
  params.erase(_pk_v8);

  p_.resize(8);
  // 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",_parallelepiped)); }
          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); }

  // (v1,v2,v4,v5) has to be set (no default values)
  if(params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_v5) != params.end()) { error("param_missing","v5"); }

  // p_ is not fully built
  p_[2]=p_[1]+p_[3]-p_[0];
  p_[5]=p_[1]+p_[4]-p_[0];
  p_[6]=p_[2]+p_[4]-p_[0];
  p_[7]=p_[3]+p_[4]-p_[0];

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

  // checking nnodes and hsteps size
  if(n_.size() != 0)
    {
      if(n_.size() == 1) { number_t n0=n_[0]; n_.resize(12,n0); }
      else if(n_.size() == 3)
        {
          number_t n0=n_[0], n1=n_[1], n2=n_[2];
          n_.clear();
          n_.resize(12, n0);
          n_[1]=n1;
          n_[3]=n1;
          n_[5]=n1;
          n_[7]=n1;
          for(number_t i=8; i < 12; ++i) { n_[i]=n2; }
        }
      else if(n_.size() != 12) {error("bad_size","nnodes",12,n_.size()); }
    }
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(8,h0); }
      else if(h_.size() != 8) {error("bad_size","hsteps",8,h_.size()); }
    }

  boundingBox=BoundingBox(p_[0], p_[1], p_[3], p_[4]);
  computeMB();
  setFaces();
  trace_p->pop();
}

void Parallelepiped::setFaces()
{
  faces_.resize(6);
  if(sideNames_.size()>=6)
    {
      faces_[0] = new Parallelogram(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[0]);
      faces_[1] = new Parallelogram(_v1=p_[4], _v2=p_[5], _v4=p_[7], _nnodes=2, _domain_name=sideNames_[1]);
      faces_[2] = new Parallelogram(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[2]);
      faces_[3] = new Parallelogram(_v1=p_[3], _v2=p_[2], _v4=p_[7], _nnodes=2, _domain_name=sideNames_[3]);
      faces_[4] = new Parallelogram(_v1=p_[0], _v2=p_[3], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[4]);
      faces_[5] = new Parallelogram(_v1=p_[1], _v2=p_[2], _v4=p_[5], _nnodes=2, _domain_name=sideNames_[5]);
    }
  else
    {
      string_t sn="";
      if(sideNames_.size()>=1) sn=sideNames_[0];
      faces_[0] = new Parallelogram(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[1] = new Parallelogram(_v1=p_[4], _v2=p_[5], _v4=p_[7], _nnodes=2, _domain_name=sn);
      faces_[2] = new Parallelogram(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[3] = new Parallelogram(_v1=p_[3], _v2=p_[2], _v4=p_[7], _nnodes=2, _domain_name=sn);
      faces_[4] = new Parallelogram(_v1=p_[0], _v2=p_[3], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[5] = new Parallelogram(_v1=p_[1], _v2=p_[2], _v4=p_[5], _nnodes=2, _domain_name=sn);
    }
}

void Parallelepiped::buildParam(const Parameter& p)
{
  Hexahedron::buildParam(p);
}

void Parallelepiped::buildDefaultParam(ParameterKey key)
{
  Hexahedron::buildDefaultParam(key);
}

std::set<ParameterKey> Parallelepiped::getParamsKeys()
{
  std::set<ParameterKey> params=Hexahedron::getParamsKeys();
  return params;
}

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

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

Parallelepiped::Parallelepiped(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                               const Parameter& p6) : Hexahedron()
{
  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);
}

Parallelepiped::Parallelepiped(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                               const Parameter& p6, const Parameter& p7) : Hexahedron()
{
  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 Parallelepiped::asString() const
{
  string_t s("Parallelepiped (");
  s += p_[0].toString() +", " + p_[1].toString() + ", " + p_[3].toString() + ", " + p_[4].toString() + ")";
  return s;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Parallelepiped::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(6);
  std::vector<const Point*> vertices(4);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[2]; vertices[3]=&p_[3];
  surfs[0]=std::make_pair(_parallelogram, vertices);
  vertices[0]=&p_[4]; vertices[1]=&p_[5]; vertices[2]=&p_[6]; vertices[3]=&p_[7];
  surfs[1]=std::make_pair(_parallelogram, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[5]; vertices[3]=&p_[4];
  surfs[2]=std::make_pair(_parallelogram, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[3]; vertices[2]=&p_[7]; vertices[3]=&p_[6];
  surfs[3]=std::make_pair(_parallelogram, vertices);
  vertices[0]=&p_[3]; vertices[1]=&p_[0]; vertices[2]=&p_[4]; vertices[3]=&p_[7];
  surfs[4]=std::make_pair(_parallelogram, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2]; vertices[2]=&p_[6]; vertices[3]=&p_[5];
  surfs[5]=std::make_pair(_parallelogram, vertices);

  return surfs;
}

real_t Parallelepiped::measure() const
{
  real_t h;
  Point P=projectionOnStraightLine(p_[2], p_[0], p_[1], h);
  real_t area=(p_[0].distance(p_[1])) * h;
  P=projectionOnTriangle(p_[4], p_[0], p_[1], p_[2], h);
  return area*h;
}

Cuboid::Cuboid() : Parallelepiped(), center_(Point(0.5,0.5,0.5)), origin_(Point(0.,0.,0.)), isCenter_(false), isOrigin_(false),
  xlength_(1.), ylength_(1.), zlength_(1.), xmin_(0.), xmax_(1.), ymin_(0.), ymax_(1.), zmin_(0.), zmax_(1.),
  isBounds_(false)
{
  shape_=_cuboid;
  computeMB();
}

void Cuboid::buildP()
{
  if(isBounds_)
    {
      p_[0]=Point(xmin_,ymin_,zmin_);
      p_[1]=Point(xmax_,ymin_,zmin_);
      p_[2]=Point(xmax_,ymax_,zmin_);
      p_[3]=Point(xmin_,ymax_,zmin_);
      p_[4]=Point(xmin_,ymin_,zmax_);
      p_[5]=Point(xmax_,ymin_,zmax_);
      p_[6]=Point(xmax_,ymax_,zmax_);
      p_[7]=Point(xmin_,ymax_,zmax_);
      origin_=p_[0];
      center_=(p_[0]+p_[6])/2.;
      xlength_=xmax_-xmin_;
      ylength_=ymax_-ymin_;
      zlength_=zmax_-zmin_;
    }
  else if(isCenter_)
    {
      Point shift0(-xlength_/2., -ylength_/2., -zlength_/2.);
      p_[0]=center_+shift0;
      Point shift1(xlength_/2., -ylength_/2., -zlength_/2.);
      p_[1]=center_+shift1;
      Point shift2(xlength_/2., ylength_/2., -zlength_/2.);
      p_[2]=center_+shift2;
      Point shift3(-xlength_/2., ylength_/2., -zlength_/2.);
      p_[3]=center_+shift3;
      Point shift4(-xlength_/2., -ylength_/2., zlength_/2.);
      p_[4]=center_+shift4;
      Point shift5(xlength_/2., -ylength_/2., zlength_/2.);
      p_[5]=center_+shift5;
      Point shift6(xlength_/2., ylength_/2., zlength_/2.);
      p_[6]=center_+shift6;
      Point shift7(-xlength_/2., ylength_/2., zlength_/2.);
      p_[7]=center_+shift7;
      origin_=p_[0];
      xmin_=xmax_=ymin_=ymax_=zmin_=zmax_=0.;
    }
  else if(isOrigin_)
    {
      p_[0]=origin_;
      Point shift1(xlength_,0.,0.);
      p_[1]=origin_+shift1;
      Point shift2(xlength_, ylength_, 0.);
      p_[2]=origin_+shift2;
      Point shift3(0., ylength_, 0.);
      p_[3]=origin_+shift3;
      Point shift4(0., 0., zlength_);
      p_[4]=origin_+shift4;
      Point shift5(xlength_, 0., zlength_);
      p_[5]=origin_+shift5;
      Point shift6(xlength_, ylength_, zlength_);
      p_[6]=origin_+shift6;
      Point shift7(0., ylength_, zlength_);
      p_[7]=origin_+shift7;
      center_=(p_[0]+p_[6])/2.;
      xmin_=xmax_=ymin_=ymax_=zmin_=zmax_=0.;
    }
  else // vertices keys are used so p_[2] is not defined
    {
      p_[2]=p_[1]+p_[3]-p_[0];
      p_[5]=p_[1]+p_[4]-p_[0];
      p_[6]=p_[2]+p_[4]-p_[0];
      p_[7]=p_[3]+p_[4]-p_[0];
      origin_=p_[0];
      center_=(p_[0]+p_[6])/2.;
      xlength_=p_[0].distance(p_[1]);
      ylength_=p_[0].distance(p_[3]);
      zlength_=p_[0].distance(p_[4]);
      xmin_=xmax_=ymin_=ymax_=zmin_=zmax_=0.;
      if(dot(p_[3]-p_[0],p_[1]-p_[0]) > theTolerance ||
          dot(p_[4]-p_[0],p_[1]-p_[0]) > theTolerance ||
          dot(p_[4]-p_[0],p_[3]-p_[0]) > theTolerance) { error("geometry_incoherent_points", words("shape",_cuboid)); }
    }
}

void Cuboid::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Cuboid::build");
  shape_=_cuboid;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // _faces key is replaced by (_v1,_v2,_v4,_v5), (_center|_origin, _xlength, _ylength, _zlength)
  // or (_xmin, _xmax, _ymin, _ymax, _zmin, _zmax)
  params.erase(_pk_faces);
  // _v3, _v6, _v7, _v8 are disabled
  params.erase(_pk_v3);
  params.erase(_pk_v6);
  params.erase(_pk_v7);
  params.erase(_pk_v8);

  p_.resize(8);
  // 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",_cuboid)); }
          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)); }
      // user must use only one of the following combination: (v1,v2,v4,v5), (xmin,xmax,ymin,ymax,zmin,zmax),
      // (center,xlength,ylength,zlength), (origin,xlength,ylength,zlength)
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && 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 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_xmax) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_xmax)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_ymin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ymin)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_ymax) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ymax)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_zmin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zmin)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_zmax) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zmax)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_center) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_center)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_origin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_origin)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_xlength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_xlength)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_ylength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ylength)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_zlength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zlength)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          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 || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_v2) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v2)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_v4) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v4)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_v5) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v5)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_center) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_center)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_origin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_origin)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_xlength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_xlength)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_ylength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ylength)); }
      if((key==_pk_xmin || key==_pk_xmax || key==_pk_ymin || key==_pk_ymax || key==_pk_zmin || key==_pk_zmax) &&
          usedParams.find(_pk_zlength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zlength)); }
      if(key==_pk_center && usedParams.find(_pk_origin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_origin)); }
      if(key==_pk_origin && usedParams.find(_pk_center) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_center)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_v1) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v1)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_v2) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v2)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_v4) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v4)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_v5) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v5)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_xmin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_xmin)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_xmax) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_xmax)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_ymin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ymin)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_ymax) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ymax)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_zmin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zmin)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) &&
          usedParams.find(_pk_zmax) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zmax)); }
    }

  // 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,v4) or (xmin,xmax,ymin,ymax) or (center,xlength,ylength) has to be set (no default values)
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v5) != params.end()) { error("param_missing","v5"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v5) != params.end()) { error("param_missing","v5"); }
  if(params.find(_pk_v4) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v4) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v4) == params.end() && params.find(_pk_v5) != params.end()) { error("param_missing","v5"); }
  if(params.find(_pk_v5) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v5) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v5) == params.end() && params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_xmin) == params.end() && params.find(_pk_xmax) != params.end()) { error("param_missing","xmax"); }
  if(params.find(_pk_xmin) == params.end() && params.find(_pk_ymin) != params.end()) { error("param_missing","ymin"); }
  if(params.find(_pk_xmin) == params.end() && params.find(_pk_ymax) != params.end()) { error("param_missing","ymax"); }
  if(params.find(_pk_xmin) == params.end() && params.find(_pk_zmin) != params.end()) { error("param_missing","zmin"); }
  if(params.find(_pk_xmin) == params.end() && params.find(_pk_zmax) != params.end()) { error("param_missing","zmax"); }
  if(params.find(_pk_xmax) == params.end() && params.find(_pk_xmin) != params.end()) { error("param_missing","xmin"); }
  if(params.find(_pk_xmax) == params.end() && params.find(_pk_ymin) != params.end()) { error("param_missing","ymin"); }
  if(params.find(_pk_xmax) == params.end() && params.find(_pk_ymax) != params.end()) { error("param_missing","ymax"); }
  if(params.find(_pk_xmax) == params.end() && params.find(_pk_zmin) != params.end()) { error("param_missing","zmin"); }
  if(params.find(_pk_xmax) == params.end() && params.find(_pk_zmax) != params.end()) { error("param_missing","zmax"); }
  if(params.find(_pk_ymin) == params.end() && params.find(_pk_xmin) != params.end()) { error("param_missing","xmin"); }
  if(params.find(_pk_ymin) == params.end() && params.find(_pk_xmax) != params.end()) { error("param_missing","xmax"); }
  if(params.find(_pk_ymin) == params.end() && params.find(_pk_ymax) != params.end()) { error("param_missing","ymax"); }
  if(params.find(_pk_ymin) == params.end() && params.find(_pk_zmin) != params.end()) { error("param_missing","zmin"); }
  if(params.find(_pk_ymin) == params.end() && params.find(_pk_zmax) != params.end()) { error("param_missing","zmax"); }
  if(params.find(_pk_ymax) == params.end() && params.find(_pk_xmin) != params.end()) { error("param_missing","xmin"); }
  if(params.find(_pk_ymax) == params.end() && params.find(_pk_xmax) != params.end()) { error("param_missing","xmax"); }
  if(params.find(_pk_ymax) == params.end() && params.find(_pk_ymin) != params.end()) { error("param_missing","ymin"); }
  if(params.find(_pk_ymax) == params.end() && params.find(_pk_zmin) != params.end()) { error("param_missing","zmin"); }
  if(params.find(_pk_ymax) == params.end() && params.find(_pk_zmax) != params.end()) { error("param_missing","zmax"); }
  if(params.find(_pk_zmin) == params.end() && params.find(_pk_xmin) != params.end()) { error("param_missing","xmin"); }
  if(params.find(_pk_zmin) == params.end() && params.find(_pk_xmax) != params.end()) { error("param_missing","xmax"); }
  if(params.find(_pk_zmin) == params.end() && params.find(_pk_ymin) != params.end()) { error("param_missing","ymin"); }
  if(params.find(_pk_zmin) == params.end() && params.find(_pk_ymax) != params.end()) { error("param_missing","ymax"); }
  if(params.find(_pk_zmin) == params.end() && params.find(_pk_zmax) != params.end()) { error("param_missing","zmax"); }
  if(params.find(_pk_zmax) == params.end() && params.find(_pk_xmin) != params.end()) { error("param_missing","xmin"); }
  if(params.find(_pk_zmax) == params.end() && params.find(_pk_xmax) != params.end()) { error("param_missing","xmax"); }
  if(params.find(_pk_zmax) == params.end() && params.find(_pk_ymin) != params.end()) { error("param_missing","ymin"); }
  if(params.find(_pk_zmax) == params.end() && params.find(_pk_ymax) != params.end()) { error("param_missing","ymax"); }
  if(params.find(_pk_zmax) == params.end() && params.find(_pk_zmin) != params.end()) { error("param_missing","zmin"); }
  if((params.find(_pk_center) == params.end() || params.find(_pk_origin) == params.end()) && params.find(_pk_xlength) != params.end())
    { error("param_missing","xlength"); }
  if((params.find(_pk_center) == params.end() || params.find(_pk_origin) == params.end()) && params.find(_pk_ylength) != params.end())
    { error("param_missing","ylength"); }
  if((params.find(_pk_center) == params.end() || params.find(_pk_origin) == params.end()) && params.find(_pk_zlength) != params.end())
    { error("param_missing","zlength"); }
  if(params.find(_pk_xlength) == params.end() && params.find(_pk_ylength) != params.end()) { error("param_missing","ylength"); }
  if(params.find(_pk_xlength) == params.end() && params.find(_pk_zlength) != params.end()) { error("param_missing","zlength"); }
  if(params.find(_pk_ylength) == params.end() && params.find(_pk_xlength) != params.end()) { error("param_missing","xlength"); }
  if(params.find(_pk_ylength) == params.end() && params.find(_pk_zlength) != params.end()) { error("param_missing","zlength"); }
  if(params.find(_pk_zlength) == params.end() && params.find(_pk_xlength) != params.end()) { error("param_missing","xlength"); }
  if(params.find(_pk_zlength) == params.end() && params.find(_pk_ylength) != params.end()) { error("param_missing","ylength"); }
  if(params.find(_pk_xlength) == params.end() && params.find(_pk_center) != params.end() && params.find(_pk_origin) != params.end())
    { error("param_missing","center or origin"); }
  isCenter_=true;
  isOrigin_=true;
  isBounds_=true;
  // now, we clean unwanted keys
  if(params.find(_pk_xlength) != params.end())
    { params.erase(_pk_xlength); params.erase(_pk_ylength);  params.erase(_pk_zlength); }
  if(params.find(_pk_center) != params.end()) { params.erase(_pk_center); isCenter_=false; }
  if(params.find(_pk_origin) != params.end()) { params.erase(_pk_origin); isOrigin_=false; }
  if(params.find(_pk_xmin) != params.end())
    {
      params.erase(_pk_xmin); params.erase(_pk_xmax);
      params.erase(_pk_ymin); params.erase(_pk_ymax);
      params.erase(_pk_zmin); params.erase(_pk_zmax);
      isBounds_=false;
    }
  if(params.find(_pk_v1) != params.end()) { params.erase(_pk_v1); params.erase(_pk_v2); params.erase(_pk_v4); params.erase(_pk_v5); }

  // p_ is not (fully) built so we do
  buildP();

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

  // checking nnodes and hsteps size
  if(n_.size() != 0)
    {
      if(n_.size() == 1) { number_t n0=n_[0]; n_.resize(12,n0); }
      else if(n_.size() == 3)
        {
          number_t n0=n_[0], n1=n_[1], n2=n_[2];
          n_.clear();
          n_.resize(12, n0);
          n_[1]=n1;
          n_[3]=n1;
          n_[5]=n1;
          n_[7]=n1;
          for(number_t i=8; i < 12; ++i) { n_[i]=n2; }
        }
      else if(n_.size() != 12) {error("bad_size","nnodes",12,n_.size()); }
    }
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(8,h0); }
      else if(h_.size() != 8) {error("bad_size","hsteps",8,h_.size()); }
    }
  boundingBox=BoundingBox(p_[0], p_[1], p_[3], p_[4]);
  computeMB();
  setFaces();
  trace_p->pop();
}

void Cuboid::setFaces()
{
  faces_.resize(6);
  if(sideNames_.size()>=6)
    {
      faces_[0] = new Rectangle(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[0]);
      faces_[1] = new Rectangle(_v1=p_[4], _v2=p_[5], _v4=p_[7], _nnodes=2, _domain_name=sideNames_[1]);
      faces_[2] = new Rectangle(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[2]);
      faces_[3] = new Rectangle(_v1=p_[3], _v2=p_[2], _v4=p_[7], _nnodes=2, _domain_name=sideNames_[3]);
      faces_[4] = new Rectangle(_v1=p_[0], _v2=p_[3], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[4]);
      faces_[5] = new Rectangle(_v1=p_[1], _v2=p_[2], _v4=p_[5], _nnodes=2, _domain_name=sideNames_[5]);
    }
  else
    {
      string_t sn="";
      if(sideNames_.size()>=1) sn=sideNames_[0];
      faces_[0] = new Rectangle(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[1] = new Rectangle(_v1=p_[4], _v2=p_[5], _v4=p_[7], _nnodes=2, _domain_name=sn);
      faces_[2] = new Rectangle(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[3] = new Rectangle(_v1=p_[3], _v2=p_[2], _v4=p_[7], _nnodes=2, _domain_name=sn);
      faces_[4] = new Rectangle(_v1=p_[0], _v2=p_[3], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[5] = new Rectangle(_v1=p_[1], _v2=p_[2], _v4=p_[5], _nnodes=2, _domain_name=sn);
    }
}


void Cuboid::buildParam(const Parameter& p)
{
  trace_p->push("Cuboid::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_center:
        {
          switch(p.type())
            {
              case _pt: center_=p.get_pt(); break;
              case _integer: center_=Point(real_t(p.get_i())); break;
              case _real: center_=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_origin:
        {
          switch(p.type())
            {
              case _pt: origin_=p.get_pt(); break;
              case _integer: origin_=Point(real_t(p.get_i())); break;
              case _real: origin_=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_xlength:
        {
          switch(p.type())
            {
              case _integer: xlength_=real_t(p.get_n()); break;
              case _real: xlength_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_ylength:
        {
          switch(p.type())
            {
              case _integer: ylength_=real_t(p.get_n()); break;
              case _real: ylength_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_zlength:
        {
          switch(p.type())
            {
              case _integer: zlength_=real_t(p.get_n()); break;
              case _real: zlength_=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: xmin_=real_t(p.get_i()); break;
              case _real: xmin_=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: xmax_=real_t(p.get_i()); break;
              case _real: xmax_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_ymin:
        {
          switch(p.type())
            {
              case _integer: ymin_=real_t(p.get_i()); break;
              case _real: ymin_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_ymax:
        {
          switch(p.type())
            {
              case _integer: ymax_=real_t(p.get_i()); break;
              case _real: ymax_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_zmin:
        {
          switch(p.type())
            {
              case _integer: zmin_=real_t(p.get_i()); break;
              case _real: zmin_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_zmax:
        {
          switch(p.type())
            {
              case _integer: zmax_=real_t(p.get_i()); break;
              case _real: zmax_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Parallelepiped::buildParam(p); break;
    }
  trace_p->pop();
}

void Cuboid::buildDefaultParam(ParameterKey key)
{
  Parallelepiped::buildDefaultParam(key);
}

std::set<ParameterKey> Cuboid::getParamsKeys()
{
  std::set<ParameterKey> params=Parallelepiped::getParamsKeys();
  params.insert(_pk_center);
  params.insert(_pk_origin);
  params.insert(_pk_xlength);
  params.insert(_pk_ylength);
  params.insert(_pk_zlength);
  params.insert(_pk_xmin);
  params.insert(_pk_xmax);
  params.insert(_pk_ymin);
  params.insert(_pk_ymax);
  params.insert(_pk_zmin);
  params.insert(_pk_zmax);
  return params;
}

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

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

Cuboid::Cuboid(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
               const Parameter& p6) : Parallelepiped()
{
  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);
}

Cuboid::Cuboid(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
               const Parameter& p6, const Parameter& p7) : Parallelepiped()
{
  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);
}

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

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

Cuboid::Cuboid(const Cuboid& c) : Parallelepiped(c), center_(c.center_), origin_(c.origin_), isCenter_(c.isCenter_),
  isOrigin_(c.isOrigin_), xlength_(c.xlength_), ylength_(c.ylength_), zlength_(c.zlength_),
  xmin_(c.xmin_), xmax_(c.xmax_), ymin_(c.ymin_), ymax_(c.ymax_), zmin_(c.zmin_), zmax_(c.zmax_),
  isBounds_(c.isBounds_)
{}

string_t Cuboid::asString() const
{
  string_t s("Cuboid (");
  s += p_[0].toString() +", " + p_[1].toString() + ", " + p_[3].toString() + ", " + p_[4].toString() + ")";
  return s;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Cuboid::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(6);
  std::vector<const Point*> vertices(4);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[2]; vertices[3]=&p_[3];
  surfs[0]=std::make_pair(_rectangle, vertices);
  vertices[0]=&p_[4]; vertices[1]=&p_[5]; vertices[2]=&p_[6]; vertices[3]=&p_[7];
  surfs[1]=std::make_pair(_rectangle, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[5]; vertices[3]=&p_[4];
  surfs[2]=std::make_pair(_rectangle, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[3]; vertices[2]=&p_[7]; vertices[3]=&p_[6];
  surfs[3]=std::make_pair(_rectangle, vertices);
  vertices[0]=&p_[3]; vertices[1]=&p_[0]; vertices[2]=&p_[4]; vertices[3]=&p_[7];
  surfs[4]=std::make_pair(_rectangle, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2]; vertices[2]=&p_[6]; vertices[3]=&p_[5];
  surfs[5]=std::make_pair(_rectangle, vertices);

  return surfs;
}

Cube::Cube() : Cuboid(), nboctants_(8)
{ shape_=_cube; }

void Cube::buildP()
{
  if(isCenter_)
    {
      Point shift0(-xlength_/2., -ylength_/2., -zlength_/2.);
      p_[0]=center_+shift0;
      Point shift1(xlength_/2., -ylength_/2., -zlength_/2.);
      p_[1]=center_+shift1;
      Point shift2(xlength_/2., ylength_/2., -zlength_/2.);
      p_[2]=center_+shift2;
      Point shift3(-xlength_/2., ylength_/2., -zlength_/2.);
      p_[3]=center_+shift3;
      Point shift4(-xlength_/2., -ylength_/2., zlength_/2.);
      p_[4]=center_+shift4;
      Point shift5(xlength_/2., -ylength_/2., zlength_/2.);
      p_[5]=center_+shift5;
      Point shift6(xlength_/2., ylength_/2., zlength_/2.);
      p_[6]=center_+shift6;
      Point shift7(-xlength_/2., ylength_/2., zlength_/2.);
      p_[7]=center_+shift7;
      origin_=p_[0];
    }
  else if(isOrigin_)
    {
      p_[0]=origin_;
      Point shift1(xlength_, 0., 0.);
      p_[1]=origin_+shift1;
      Point shift2(xlength_, ylength_, 0.);
      p_[2]=origin_+shift2;
      Point shift3(0., ylength_, 0.);
      p_[3]=origin_+shift3;
      Point shift4(0., 0., zlength_);
      p_[4]=origin_+shift4;
      Point shift5(xlength_, 0., zlength_);
      p_[5]=origin_+shift5;
      Point shift6(xlength_, ylength_, zlength_);
      p_[6]=origin_+shift6;
      Point shift7(0., ylength_, zlength_);
      p_[7]=origin_+shift7;
      center_=(p_[0]+p_[6])/2.;
    }
  else // vertices keys are used so p_[2] is not defined
    {
      p_[2]=p_[1]+p_[3]-p_[0];
      p_[5]=p_[1]+p_[4]-p_[0];
      p_[6]=p_[2]+p_[4]-p_[0];
      p_[7]=p_[3]+p_[4]-p_[0];
      origin_=p_[0];
      center_=(p_[0]+p_[6])/2.;
      xlength_=p_[0].distance(p_[1]);
      ylength_=p_[0].distance(p_[3]);
      zlength_=p_[0].distance(p_[4]);
      if(dot(p_[3]-p_[0],p_[1]-p_[0]) > theTolerance ||
          dot(p_[4]-p_[0],p_[1]-p_[0]) > theTolerance ||
          dot(p_[4]-p_[0],p_[3]-p_[0]) > theTolerance) { error("geometry_incoherent_points", words("shape",_cube)); }
    }
}

void Cube::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Cube::build");
  shape_=_cube;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // _facces key is replaced by (_v1,_v2,_v4), (_center|_origin, _xlength, _ylength) or (_xmin, _xmax, _ymin, _ymax)
  params.erase(_pk_faces);
  // _v3, _v6, _v7, _v8, _xlength, _ylength, _zlength, _xmin, _xmax, _ymin, _ymax, _zmin and _zmax are disabled
  params.erase(_pk_v3);
  params.erase(_pk_v6);
  params.erase(_pk_v7);
  params.erase(_pk_v8);
  params.erase(_pk_xlength);
  params.erase(_pk_ylength);
  params.erase(_pk_zlength);
  params.erase(_pk_xmin);
  params.erase(_pk_xmax);
  params.erase(_pk_ymin);
  params.erase(_pk_ymax);
  params.erase(_pk_zmin);
  params.erase(_pk_zmax);

  p_.resize(8);
  // 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",_cube)); }
          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)); }
      // user must use only one of the following combination: (v1,v2,v4), (center,length) or (origin,length)
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_center) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_center)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_origin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_origin)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v4 || key==_pk_v5) && usedParams.find(_pk_length) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_length)); }
      if(key==_pk_center && usedParams.find(_pk_origin) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_origin)); }
      if(key==_pk_origin && usedParams.find(_pk_center) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_center)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_length) && usedParams.find(_pk_v1) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v1)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_length) && usedParams.find(_pk_v2) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v2)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_length) && usedParams.find(_pk_v4) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v4)); }
      if((key==_pk_center || key==_pk_origin || key==_pk_length) && usedParams.find(_pk_v5) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v5)); }
    }

  // 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,v4) or (xmin,xmax,ymin,ymax) or (center,xlength,ylength) has to be set (no default values)
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v4) != params.end()) { error("param_missing","v4"); }
  if(params.find(_pk_v4) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v4) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if((params.find(_pk_center) == params.end() || params.find(_pk_origin) == params.end()) && params.find(_pk_length) != params.end())
    { error("param_missing","length"); }
  if(params.find(_pk_length) == params.end() && params.find(_pk_center) != params.end() && params.find(_pk_origin) != params.end())
    { error("param_missing","center or origin"); }
  isCenter_=true;
  isOrigin_=true;
  // now, we clean unwanted keys
  if(params.find(_pk_length) != params.end()) { params.erase(_pk_length); }
  if(params.find(_pk_center) != params.end()) { params.erase(_pk_center); isCenter_=false; }
  if(params.find(_pk_origin) != params.end()) { params.erase(_pk_origin); isOrigin_=false; }
  if(params.find(_pk_v1) != params.end()) { params.erase(_pk_v1); params.erase(_pk_v2); params.erase(_pk_v4); params.erase(_pk_v5); }

  // p_ is not (fully) built so we do
  buildP();
  if(std::abs(p_[0].distance(p_[1])-p_[0].distance(p_[3])) > theTolerance ||
      std::abs(p_[0].distance(p_[1])-p_[0].distance(p_[4])) > theTolerance ||
      std::abs(p_[0].distance(p_[3])-p_[0].distance(p_[4])) > theTolerance)
    { error("geometry_incoherent_points", words("shape",_cube)); }

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

  // checking nnodes and hsteps size
  if(n_.size() != 0)
    {
      if(n_.size() == 1) { number_t n0=n_[0]; n_.resize(12,n0); }
      else if(n_.size() == 3)
        {
          number_t n0=n_[0], n1=n_[1], n2=n_[2];
          n_.clear();
          n_.resize(12, n0);
          n_[1]=n1;
          n_[3]=n1;
          n_[5]=n1;
          n_[7]=n1;
          for(number_t i=8; i < 12; ++i) { n_[i]=n2; }
        }
      else if(n_.size() != 12) {error("bad_size","nnodes",12,n_.size()); }
    }
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(8,h0); }
      else if(h_.size() != 8) {error("bad_size","hsteps",8,h_.size()); }
    }

  boundingBox=BoundingBox(p_[0], p_[1], p_[3], p_[4]);
  computeMB();
  setFaces();
  trace_p->pop();
}

void Cube::setFaces()
{
  faces_.resize(6);
  if(sideNames_.size()>=6)
    {
      faces_[0] = new Square(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[0]);
      faces_[1] = new Square(_v1=p_[4], _v2=p_[5], _v4=p_[7], _nnodes=2, _domain_name=sideNames_[1]);
      faces_[2] = new Square(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[2]);
      faces_[3] = new Square(_v1=p_[3], _v2=p_[2], _v4=p_[7], _nnodes=2, _domain_name=sideNames_[3]);
      faces_[4] = new Square(_v1=p_[0], _v2=p_[3], _v4=p_[4], _nnodes=2, _domain_name=sideNames_[4]);
      faces_[5] = new Square(_v1=p_[1], _v2=p_[2], _v4=p_[5], _nnodes=2, _domain_name=sideNames_[5]);
    }
  else
    {
      string_t sn="";
      if(sideNames_.size()>=1) sn=sideNames_[0];
      faces_[0] = new Square(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[1] = new Square(_v1=p_[4], _v2=p_[5], _v4=p_[7], _nnodes=2, _domain_name=sn);
      faces_[2] = new Square(_v1=p_[0], _v2=p_[1], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[3] = new Square(_v1=p_[3], _v2=p_[2], _v4=p_[7], _nnodes=2, _domain_name=sn);
      faces_[4] = new Square(_v1=p_[0], _v2=p_[3], _v4=p_[4], _nnodes=2, _domain_name=sn);
      faces_[5] = new Square(_v1=p_[1], _v2=p_[2], _v4=p_[5], _nnodes=2, _domain_name=sn);
    }
}

void Cube::buildParam(const Parameter& p)
{
  trace_p->push("Cube::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_length:
        {
          switch(p.type())
            {
              case _integer: xlength_=ylength_=zlength_=real_t(p.get_n()); break;
              case _real: xlength_=ylength_=zlength_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_nboctants:
        {
          switch(p.type())
            {
              case _integer: nboctants_ = p.get_n(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Cuboid::buildParam(p); break;
    }
  trace_p->pop();
}

void Cube::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Cube::buildDefaultParam");
  switch(key)
    {
      case _pk_nboctants: nboctants_=8; break;
      default: Cuboid::buildDefaultParam(key); break;
    }
  trace_p->pop();
}

std::set<ParameterKey> Cube::getParamsKeys()
{
  std::set<ParameterKey> params=Cuboid::getParamsKeys();
  params.insert(_pk_length);
  params.insert(_pk_nboctants);
  return params;
}

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

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

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

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

Cube::Cube(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
           const Parameter& p6) : Cuboid()
{
  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);
}

Cube::Cube(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
           const Parameter& p6, const Parameter& p7) : Cuboid()
{
  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);
}

Cube::Cube(const Cube& c) : Cuboid(c), nboctants_(c.nboctants_)
{}

number_t Cube::nbSubdiv() const
{
  number_t n0=n(1);
  for(number_t i=1; i < 12; ++i) { if(n(i+1) > n0) { n0=n(i+1); } }
  return static_cast<number_t>(theTolerance+std::log(n0-1)/std::log(2));
}

//! Format Cube as string : "center = (.,.,.), edge length = L"
string_t Cube::asString() const
{
  string_t s("Cube (");
  s += p_[0].toString() +", " + p_[1].toString() + ", " + p_[3].toString() + ", " + p_[4].toString() + ")";
  return s;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Cube::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(6);
  std::vector<const Point*> vertices(4);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[2]; vertices[3]=&p_[3];
  surfs[0]=std::make_pair(_square, vertices);
  vertices[0]=&p_[4]; vertices[1]=&p_[5]; vertices[2]=&p_[6]; vertices[3]=&p_[7];
  surfs[1]=std::make_pair(_square, vertices);
  vertices[0]=&p_[0]; vertices[1]=&p_[1]; vertices[2]=&p_[5]; vertices[3]=&p_[4];
  surfs[2]=std::make_pair(_square, vertices);
  vertices[0]=&p_[2]; vertices[1]=&p_[3]; vertices[2]=&p_[7]; vertices[3]=&p_[6];
  surfs[3]=std::make_pair(_square, vertices);
  vertices[0]=&p_[3]; vertices[1]=&p_[0]; vertices[2]=&p_[4]; vertices[3]=&p_[7];
  surfs[4]=std::make_pair(_square, vertices);
  vertices[0]=&p_[1]; vertices[1]=&p_[2]; vertices[2]=&p_[6]; vertices[3]=&p_[5];
  surfs[5]=std::make_pair(_square, vertices);

  return surfs;
}

//! default ellipsoid is sphere of center (0,0,0) and radius 1
Ellipsoid::Ellipsoid() : Volume(), center_(Point(0.,0.,0.)), p1_(Point(1.,0.,0.)), p2_(Point(0.,1.,0.)), p3_(Point(-1.,0.,0.)),
  p4_(Point(0.,-1.,0.)), p5_(Point(0.,0.,-1.)), p6_(Point(0.,0.,1.)), xlength_(1.), ylength_(1.), zlength_(1.),
  isAxis_(false), n1_(2), n2_(2), n3_(2), n4_(2), n5_(2), n6_(2), n7_(2), n8_(2), n9_(2), n10_(2), n11_(2),
  n12_(2), nboctants_(8), type_(1)
{
  shape_=_ellipsoid;
  computeMB();
}

void Ellipsoid::buildP()
{
  if(isAxis_)
    {
      Point shift1(xlength_/2., 0., 0.);
      p1_=center_+shift1;
      Point shift2(0., ylength_/2., 0.);
      p2_=center_+shift2;
      Point shift3(-xlength_/2., 0., 0.);
      p3_=center_+shift3;
      Point shift4(0., -ylength_/2., 0.);
      p4_=center_+shift4;
      Point shift5(0., 0., -zlength_/2.);
      p5_=center_+shift5;
      Point shift6(0., 0., zlength_/2.);
      p6_=center_+shift6;
    }
  else // vertices keys are used so p3_ and p4_ are not defined
    {
      p3_=2.*center_-p1_;
      p4_=2.*center_-p2_;
      p5_=2.*center_-p6_;
      xlength_=2.*center_.distance(p1_);
      ylength_=2.*center_.distance(p2_);
      zlength_=2.*center_.distance(p6_);
      if(dot(p2_-center_,p1_-center_) > theTolerance ||
          dot(p6_-center_,p1_-center_) > theTolerance ||
          dot(p6_-center_,p2_-center_) > theTolerance) { error("geometry_incoherent_points", words("shape",shape_)); }
    }
}

void Ellipsoid::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Ellipsoid::build");
  shape_=_ellipsoid;
  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",_ellipsoid)); }
          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)); }
      // user must use only one of the following combination: (center,v1,v2) or (center,xlength,ylength)
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v6) && usedParams.find(_pk_xlength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_xlength)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v6) && usedParams.find(_pk_ylength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_ylength)); }
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v6) && usedParams.find(_pk_zlength) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_zlength)); }
      if((key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) && usedParams.find(_pk_v1) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v1)); }
      if((key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) && usedParams.find(_pk_v2) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v2)); }
      if((key==_pk_xlength || key==_pk_ylength || key==_pk_zlength) && usedParams.find(_pk_v6) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v6)); }
    }

  // 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,v2,v6) or (center,xlength,ylength,zlength) has to be set (no default values)
  if(params.find(_pk_center) != params.end()) { error("param_missing","center"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v6) != params.end()) { error("param_missing","v6"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v6) != params.end()) { error("param_missing","v6"); }
  if(params.find(_pk_v6) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v6) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_xlength) == params.end() && params.find(_pk_ylength) != params.end()) { error("param_missing","ylength"); }
  if(params.find(_pk_xlength) == params.end() && params.find(_pk_zlength) != params.end()) { error("param_missing","zlength"); }
  if(params.find(_pk_ylength) == params.end() && params.find(_pk_xlength) != params.end()) { error("param_missing","xlength"); }
  if(params.find(_pk_ylength) == params.end() && params.find(_pk_zlength) != params.end()) { error("param_missing","zlength"); }
  if(params.find(_pk_zlength) == params.end() && params.find(_pk_xlength) != params.end()) { error("param_missing","xlength"); }
  if(params.find(_pk_zlength) == params.end() && params.find(_pk_ylength) != params.end()) { error("param_missing","ylength"); }

  isAxis_=true;
  // now, we clean unwanted keys
  if(params.find(_pk_xlength) != params.end())
    { params.erase(_pk_xlength); params.erase(_pk_ylength); params.erase(_pk_zlength); isAxis_=false; }
  if(params.find(_pk_v1) != params.end()) { params.erase(_pk_v1); params.erase(_pk_v2); params.erase(_pk_v6); }
  // p_ is not (fully) built so we do
  buildP();

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

  // checking hsteps size
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(6,h0); }
      else if(h_.size() != 6) {error("bad_size","hsteps",6,h_.size()); }
    }

  boundingBox=BoundingBox(4.*center_-p1_-p2_-p6_, 2.*center_+p1_-p2_-p6_, 2.*center_+p2_-p1_-p6_, 2.*center_+p6_-p1_-p2_);
  computeMB();
  trace_p->pop();
}

void Ellipsoid::buildParam(const Parameter& p)
{
  trace_p->push("Ellipsoid::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_center:
        {
          switch(p.type())
            {
              case _pt: center_=p.get_pt(); break;
              case _integer: center_=Point(real_t(p.get_i())); break;
              case _real: center_=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_v6:
        {
          switch(p.type())
            {
              case _pt: p6_=p.get_pt(); break;
              case _integer: p6_=Point(real_t(p.get_i())); break;
              case _real: p6_=Point(p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_xlength:
        {
          switch(p.type())
            {
              case _integer: xlength_=real_t(p.get_n()); break;
              case _real: xlength_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_ylength:
        {
          switch(p.type())
            {
              case _integer: ylength_=real_t(p.get_n()); break;
              case _real: ylength_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_zlength:
        {
          switch(p.type())
            {
              case _integer: zlength_=real_t(p.get_n()); break;
              case _real: zlength_=p.get_r(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_type:
        {
          switch(p.type())
            {
              case _integer: type_=dimen_t(p.get_n()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_nnodes:
        {
          switch(p.type())
            {
              case _integerVector:
                {
                  std::vector<number_t> n=p.get_nv();
                  if(n.size() == 1) { n1_=n2_=n3_=n4_=n5_=n6_=n7_=n8_=n9_=n10_=n11_=n12_=n[0] > 2 ? n[0] : 2; }
                  else if(n.size() == 3)
                    {
                      n1_=n2_=n3_=n4_=n[0] > 2 ? n[0] : 2;
                      n5_=n6_=n7_=n8_=n[1] > 2 ? n[1] : 2;
                      n9_=n10_=n11_=n12_=n[2] > 2 ? n[2] : 2;
                    }
                  else if(n.size() != 12) { error("bad_size","nnodes",12,n.size()); }
                  else
                    {
                      n1_=n[0] > 2 ? n[0] : 2;
                      n2_=n[1] > 2 ? n[1] : 2;
                      n3_=n[2] > 2 ? n[2] : 2;
                      n4_=n[3] > 2 ? n[3] : 2;
                      n5_=n[4] > 2 ? n[4] : 2;
                      n6_=n[5] > 2 ? n[5] : 2;
                      n7_=n[6] > 2 ? n[6] : 2;
                      n8_=n[7] > 2 ? n[7] : 2;
                      n9_=n[8] > 2 ? n[8] : 2;
                      n10_=n[9] > 2 ? n[9] : 2;
                      n11_=n[10] > 2 ? n[10] : 2;
                      n12_=n[11] > 2 ? n[11] : 2;
                    }
                  break;
                }
              case _integer: n1_=n2_=n3_=n4_=n5_=n6_=n7_=n8_=n9_=n10_=n11_=n12_=p.get_n() > 2 ? p.get_n() : 2; break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_hsteps:
        {
          switch(p.type())
            {
              case _realVector: h_=p.get_rv(); break;
              case _integer: h_=std::vector<real_t>(1,real_t(p.get_n())); break;
              case _real: h_=std::vector<real_t>(1,p.get_r()); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      case _pk_nboctants:
        {
          switch(p.type())
            {
              //case _integer: nboctants_=p.get_d(); break;
              case _integer: nboctants_ = p.get_n(); break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Volume::buildParam(p); break;
    }
  trace_p->pop();
}

void Ellipsoid::buildDefaultParam(ParameterKey key)
{
  trace_p->push("Ellipsoid::buildDefaultParam");
  switch(key)
    {
      case _pk_nnodes: n1_=n2_=n3_=n4_=n5_=n6_=n7_=n8_=n9_=n10_=n11_=n12_=2; break;
      case _pk_type: type_=1; break;
      case _pk_nboctants: nboctants_=0; break;
      default: Volume::buildDefaultParam(key); break;
    }
  trace_p->pop();
}

std::set<ParameterKey> Ellipsoid::getParamsKeys()
{
  std::set<ParameterKey> params=Volume::getParamsKeys();
  params.insert(_pk_center);
  params.insert(_pk_v1);
  params.insert(_pk_v2);
  params.insert(_pk_v6);
  params.insert(_pk_xlength);
  params.insert(_pk_ylength);
  params.insert(_pk_zlength);
  params.insert(_pk_nnodes);
  params.insert(_pk_hsteps);
  params.insert(_pk_nboctants);
  params.insert(_pk_type);
  return params;
}

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

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

Ellipsoid::Ellipsoid(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                     const Parameter& p6) : Volume()
{
  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);
}

Ellipsoid::Ellipsoid(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
                     const Parameter& p6, const Parameter& p7) : Volume()
{
  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);
}

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

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

Ellipsoid::Ellipsoid(const Ellipsoid& e) : Volume(e), center_(e.center_), p1_(e.p1_), p2_(e.p2_), p3_(e.p3_), p4_(e.p4_), p5_(e.p5_),
  p6_(e.p6_), xlength_(e.xlength_), ylength_(e.ylength_), zlength_(e.zlength_),
  isAxis_(e.isAxis_), n1_(e.n1_), n2_(e.n2_), n3_(e.n3_), n4_(e.n4_), n5_(e.n5_), n6_(e.n6_),
  n7_(e.n7_), n8_(e.n8_), n9_(e.n9_), n10_(e.n10_), n11_(e.n11_), n12_(e.n12_), h_(e.h_),
  nboctants_(e.nboctants_), type_(e.type_)
{}

number_t Ellipsoid::nbSubdiv() const
{
  number_t n=n1_;
  if(n2_ > n) { n=n2_; }
  if(n3_ > n) { n=n3_; }
  if(n4_ > n) { n=n4_; }
  if(n5_ > n) { n=n5_; }
  if(n6_ > n) { n=n6_; }
  if(n7_ > n) { n=n7_; }
  if(n8_ > n) { n=n8_; }
  if(n9_ > n) { n=n9_; }
  if(n10_ > n) { n=n10_; }
  if(n11_ > n) { n=n11_; }
  if(n12_ > n) { n=n12_; }
  return static_cast<number_t>(theTolerance+std::log(n-1)/std::log(2));
}

string_t Ellipsoid::asString() const
{
  string_t s("Ellipsoid (center = ");
  s += center_.toString() +", ";
  s += "1st apogee = " + p1_.toString() + ", 2nd apogee = " + p2_.toString() + ", 3rd apogee = " + p6_.toString() + ")";
  return s;
}

std::vector<const Point*> Ellipsoid::boundNodes() const
{
  std::vector<const Point*> nodes(6);
  nodes[0]=&p1_; nodes[1]=&p2_; nodes[2]=&p3_; nodes[3]=&p4_; nodes[4]=&p5_; nodes[5]=&p6_;
  return nodes;
}

std::vector<Point*> Ellipsoid::nodes()
{
  std::vector<Point*> nodes(7);
  nodes[0]=&center_; nodes[1]=&p1_; nodes[2]=&p2_; nodes[3]=&p3_; nodes[4]=&p4_; nodes[5]=&p5_; nodes[6]=&p6_;
  return nodes;
}

std::vector<const Point*> Ellipsoid::nodes() const
{
  std::vector<const Point*> nodes(7);
  nodes[0]=&center_; nodes[1]=&p1_; nodes[2]=&p2_; nodes[3]=&p3_; nodes[4]=&p4_; nodes[5]=&p5_; nodes[6]=&p6_;
  return nodes;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Ellipsoid::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(12);
  std::vector<const Point*> vertices(2);
  vertices[0]=&p1_; vertices[1]=&p2_;
  curves[0]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p2_; vertices[1]=&p3_;
  curves[1]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p3_; vertices[1]=&p4_;
  curves[2]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p4_; vertices[1]=&p1_;
  curves[3]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p1_; vertices[1]=&p6_;
  curves[4]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p6_; vertices[1]=&p3_;
  curves[5]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p3_; vertices[1]=&p5_;
  curves[6]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p5_; vertices[1]=&p1_;
  curves[7]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p2_; vertices[1]=&p6_;
  curves[8]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p6_; vertices[1]=&p4_;
  curves[9]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p4_; vertices[1]=&p5_;
  curves[10]=std::make_pair(_ellArc,vertices);
  vertices[0]=&p5_; vertices[1]=&p2_;
  curves[11]=std::make_pair(_ellArc,vertices);
  return curves;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Ellipsoid::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > surfs(8);
  std::vector<const Point*> vertices(3);
  vertices[0]=&p1_; vertices[1]=&p2_; vertices[2]=&p6_;
  surfs[0]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p1_; vertices[1]=&p2_; vertices[2]=&p5_;
  surfs[1]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p3_; vertices[1]=&p2_; vertices[2]=&p6_;
  surfs[2]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p3_; vertices[1]=&p2_; vertices[2]=&p5_;
  surfs[3]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p1_; vertices[1]=&p4_; vertices[2]=&p6_;
  surfs[4]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p1_; vertices[1]=&p4_; vertices[2]=&p5_;
  surfs[5]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p3_; vertices[1]=&p4_; vertices[2]=&p6_;
  surfs[6]=std::make_pair(_ellipsoidPart,vertices);
  vertices[0]=&p3_; vertices[1]=&p4_; vertices[2]=&p5_;
  surfs[7]=std::make_pair(_ellipsoidPart,vertices);
  return surfs;
}

real_t Ellipsoid::measure() const
{
  number_t nboct=nboctants_;
  if (nboctants_ == 0) { nboct=8; }
  return pi_*xlength_*ylength_*zlength_/46.*nboct;
}

//! default is ball of center (0,0,0) and radius 1
Ball::Ball() : Ellipsoid()
{ shape_ = _ball; }

void Ball::build(const std::vector<Parameter>& ps)
{
  trace_p->push("Ball::build");
  shape_=_ball;
  std::set<ParameterKey> params=getParamsKeys(), usedParams;
  // xlength, ylength and zlength are replaced by radius
  params.erase(_pk_xlength);
  params.erase(_pk_ylength);
  params.erase(_pk_zlength);

  // 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",_ball)); }
          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)); }
      // user must use only one of the following combination: (center,v1,v2) or (center,xlength,ylength)
      if((key==_pk_v1 || key==_pk_v2 || key==_pk_v6) && usedParams.find(_pk_radius) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_radius)); }
      if((key==_pk_radius) && usedParams.find(_pk_v1) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v1)); }
      if((key==_pk_radius) && usedParams.find(_pk_v2) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v2)); }
      if((key==_pk_radius) && usedParams.find(_pk_v6) != usedParams.end())
        { error("param_conflict", words("param key",key), words("param key",_pk_v6)); }
    }

  // 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,v2) or (center,xlength,ylength) has to be set (no default values)
  if(params.find(_pk_center) != params.end()) { error("param_missing","center"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }
  if(params.find(_pk_v1) == params.end() && params.find(_pk_v6) != params.end()) { error("param_missing","v6"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v2) == params.end() && params.find(_pk_v6) != params.end()) { error("param_missing","v6"); }
  if(params.find(_pk_v6) == params.end() && params.find(_pk_v1) != params.end()) { error("param_missing","v1"); }
  if(params.find(_pk_v6) == params.end() && params.find(_pk_v2) != params.end()) { error("param_missing","v2"); }

  isAxis_=true;
  // now, we clean unwanted keys
  if(params.find(_pk_radius) != params.end()) { params.erase(_pk_radius); isAxis_=false; }
  if(params.find(_pk_v1) != params.end()) { params.erase(_pk_v1); params.erase(_pk_v2); params.erase(_pk_v6); }

  // p_ is not (fully) built so we do
  buildP();

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

  if((std::abs(center_.distance(p1_) - center_.distance(p2_)) > theTolerance) ||
      (std::abs(center_.distance(p1_) - center_.distance(p6_)) > theTolerance) ||
      (std::abs(center_.distance(p2_) - center_.distance(p6_)) > theTolerance))
    { error("geometry_incoherent_points", words("shape",_ball)); }

  // checking hsteps size
  if(h_.size() != 0)
    {
      if(h_.size() == 1) { real_t h0=h_[0]; h_.resize(6,h0); }
      else if(h_.size() != 6) {error("bad_size","hsteps",6,h_.size()); }
    }

  boundingBox=BoundingBox(4.*center_-p1_-p2_-p6_, 2.*center_+p1_-p2_-p6_, 2.*center_+p2_-p1_-p6_, 2.*center_+p6_-p1_-p2_);
  computeMB();
  trace_p->pop();
}

void Ball::buildParam(const Parameter& p)
{
  trace_p->push("Ball::buildParam");
  ParameterKey key=p.key();
  switch(key)
    {
      case _pk_radius:
        {
          switch(p.type())
            {
              case _integer: xlength_=ylength_=zlength_=real_t(p.get_i())*2.; break;
              case _real: xlength_=ylength_=zlength_=p.get_r()*2.; break;
              default: error("param_badtype",words("value",p.type()),words("param key",key));
            }
          break;
        }
      default: Ellipsoid::buildParam(p); break;
    }
  trace_p->pop();
}

void Ball::buildDefaultParam(ParameterKey key)
{
  Ellipsoid::buildDefaultParam(key);
}

std::set<ParameterKey> Ball::getParamsKeys()
{
  std::set<ParameterKey> params=Ellipsoid::getParamsKeys();
  params.insert(_pk_radius);
  return params;
}

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

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

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

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

Ball::Ball(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
           const Parameter& p6) : Ellipsoid()
{
  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);
}

Ball::Ball(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4, const Parameter& p5,
           const Parameter& p6, const Parameter& p7) : Ellipsoid()
{
  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);
}

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

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

Ball::Ball(const Ball& b) : Ellipsoid(b) {}

//! Format Ball as string : "Ball (center = (.,.,.), radius = R)"
string_t Ball::asString() const
{
  string_t s("Ball (center = ");
  s += center_.toString() +", ";
  s += "1st point = " + p1_.toString() + ", 2nd point = " + p2_.toString() + ", 3rd point = " + p6_.toString() + ")";
  return s;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Ball::curves() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(12);
  std::vector<const Point*> vertices(2);
  vertices[0]=&p1_; vertices[1]=&p2_;
  curves[0]=std::make_pair(_circArc,vertices);
  vertices[0]=&p2_; vertices[1]=&p3_;
  curves[1]=std::make_pair(_circArc,vertices);
  vertices[0]=&p3_; vertices[1]=&p4_;
  curves[2]=std::make_pair(_circArc,vertices);
  vertices[0]=&p4_; vertices[1]=&p1_;
  curves[3]=std::make_pair(_circArc,vertices);
  vertices[0]=&p1_; vertices[1]=&p6_;
  curves[4]=std::make_pair(_circArc,vertices);
  vertices[0]=&p6_; vertices[1]=&p3_;
  curves[5]=std::make_pair(_circArc,vertices);
  vertices[0]=&p3_; vertices[1]=&p5_;
  curves[6]=std::make_pair(_circArc,vertices);
  vertices[0]=&p5_; vertices[1]=&p1_;
  curves[7]=std::make_pair(_circArc,vertices);
  vertices[0]=&p2_; vertices[1]=&p6_;
  curves[8]=std::make_pair(_circArc,vertices);
  vertices[0]=&p6_; vertices[1]=&p4_;
  curves[9]=std::make_pair(_circArc,vertices);
  vertices[0]=&p4_; vertices[1]=&p5_;
  curves[10]=std::make_pair(_circArc,vertices);
  vertices[0]=&p5_; vertices[1]=&p2_;
  curves[11]=std::make_pair(_circArc,vertices);
  return curves;
}

std::vector<std::pair<ShapeType,std::vector<const Point*> > > Ball::surfs() const
{
  std::vector<std::pair<ShapeType,std::vector<const Point*> > > curves(8);
  std::vector<const Point*> vertices(3);
  vertices[0]=&p1_; vertices[1]=&p2_; vertices[2]=&p6_;
  curves[0]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p1_; vertices[1]=&p2_; vertices[2]=&p5_;
  curves[1]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p3_; vertices[1]=&p2_; vertices[2]=&p6_;
  curves[2]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p3_; vertices[1]=&p2_; vertices[2]=&p5_;
  curves[3]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p1_; vertices[1]=&p4_; vertices[2]=&p6_;
  curves[4]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p1_; vertices[1]=&p4_; vertices[2]=&p5_;
  curves[5]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p3_; vertices[1]=&p4_; vertices[2]=&p6_;
  curves[6]=std::make_pair(_spherePart,vertices);
  vertices[0]=&p3_; vertices[1]=&p4_; vertices[2]=&p5_;
  curves[7]=std::make_pair(_spherePart,vertices);
  return curves;
}

} // end of namespace xlifepp