xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/geometry/GeomElement.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 GeomElement.cpp
  \authors D.Martin, E. Lunéville, N. Kielbasiewicz
  \since 14 apr 2012
  \date 5 jul 2013

  \brief Implementation of xlifepp::GeomElement classes and functionnalities
*/

#include "GeomElement.hpp"
#include "Mesh.hpp"

namespace xlifepp
{

//===========================================================================
// GeomElement member functions and related external functions
//===========================================================================

//construct a plain element from ReferenceElement & space dimension
GeomElement::GeomElement(const Mesh* msh, const RefElement* ref, dimen_t d, number_t n)
  : mesh_p(msh), number_(n), twin_p(0)
{
  meshElement_p = new MeshElement(ref, d, n);
  materialId=0; domainId=0;
  theta=0; phi=0;
  color=0; flag=0;
}

//construct a side element from GeomElement and a side number
GeomElement::GeomElement(GeomElement* parent, number_t sideNum, number_t index)
  : number_(index), twin_p(0)
{
  mesh_p = parent->mesh_p;
  parentSides_.push_back(GeoNumPair(parent, sideNum));
  meshElement_p = 0;
  materialId=0; domainId=0;
  theta=0; phi=0.;
  color=0;flag=0;
}

//copy constructor
GeomElement::GeomElement(const GeomElement& ge)
{
  meshElement_p = 0;
  if(ge.meshElement_p != 0) { meshElement_p = new MeshElement(*ge.meshElement_p);}
  if(ge.parentSides_.size() > 0) { parentSides_ = ge.parentSides_; }
  number_ = ge.number_;
  mesh_p = ge.mesh_p;
  materialId = ge.materialId;
  domainId=ge.domainId;
  theta=ge.theta; phi=ge.phi;
  color=ge.color;flag=ge.flag;
  twin_p = ge.twin_p;
}

//assignment operator
GeomElement& GeomElement::operator=(const GeomElement& ge)
{
  if(this == &ge) { return *this; }
  if(meshElement_p != 0) { delete meshElement_p; }
  meshElement_p = 0;
  if(ge.meshElement_p != 0) { meshElement_p = new MeshElement(*ge.meshElement_p); }
  parentSides_.clear();
  if(ge.parentSides_.size() > 0) { parentSides_ = ge.parentSides_; }
  number_ = ge.number_;
  materialId = ge.materialId;
  domainId=ge.domainId;
  theta=ge.theta; phi=ge.phi;
  color=ge.color; flag=ge.flag;
  twin_p = ge.twin_p;
  return *this;
}

//delete meshElement information
void GeomElement::deleteMeshElement()
{
  if(meshElement_p != 0) { delete meshElement_p; }
  meshElement_p = 0;
}

//---------------------------------------------------------------------------
//build MeshElement structure for a side element if required
//---------------------------------------------------------------------------
MeshElement* GeomElement::buildSideMeshElement() const
{
  if(meshElement_p != 0) { return meshElement_p;}   //already built
  trace_p->push("GeomElement::buildSideMeshElement");

  //new MeshElement
  RefElement* ref = findRefElement(shapeType(), &refElement()->interpolation());
  //meshElement_p = new MeshElement(ref, geomRefElement()->dim(), 1);
  dimen_t dim=geomRefElement()->dim(); //case of mesh_p=0 (normally does not occur)
  if(mesh_p!=0) { dim= mesh_p->spaceDim();}
  meshElement_p = new MeshElement(ref,dim,1);
  //set nodenumbers from sideDofNumbers of RefElement parent (assuming a Lagrange element)
  GeomElement* parent = parentSides_[0].first;
  number_t side = parentSides_[0].second;
  if(parent->meshElement_p != 0) {
      // set nodeNumbers and nodes
      std::vector<number_t> locnum = parent->refElement()->sideDofNumbers_[side - 1]; //local node numbering on side
      size_t lns = locnum.size();
      meshElement_p->nodeNumbers.resize(lns);
      meshElement_p->nodes.resize(lns);
      for(number_t i = 0; i < lns; i++) {
        meshElement_p->nodeNumbers[i] = parent->meshElement_p->nodeNumbers[locnum[i] - 1];
        meshElement_p->nodes[i] = parent->meshElement_p->nodes[locnum[i] - 1];
      }
      // set vertexNumbers
      const std::vector<number_t>& svn = parent->geomRefElement()->sideVertexNumbers()[side - 1];
      for(number_t i = 0; i < svn.size(); i++)
        { meshElement_p->vertexNumbers[i] = parent->meshElement_p->vertexNumber(svn[i]); }
      //update measures
      meshElement_p->computeMeasures();
      //update linearMap flag
      if(!meshElement_p->linearMap) meshElement_p->linearMap=meshElement_p->checkLinearMap();
      trace_p->pop();
      return meshElement_p;
    }
  // parent has no MeshElement structure, go to grand parent (case of a side of side, edge in 3D)
  number_t sideOfSide = parent->parentSides_[0].second;
  parent = parent->parentSides_[0].first;
  if(parent->meshElement_p == 0) { error("geoelt_sosos", "GeomElement::buildSideMeshElement"); }   //stop recursion of side of side
  // set nodeNumbers and nodes
  number_t sos = std::abs(parent->geomRefElement()->sideOfSideNumber(sideOfSide, side));
  std::vector<number_t> locnum = parent->refElement()->sideOfSideDofNumbers_[sos - 1]; //local node numbering on side of side
  size_t lns = locnum.size();
  meshElement_p->nodeNumbers.resize(lns);
  meshElement_p->nodes.resize(lns);
  for(number_t i = 0; i < lns; i++) {
    meshElement_p->nodeNumbers[i] = parent->meshElement_p->nodeNumbers[locnum[i] - 1];
    meshElement_p->nodes[i] = parent->meshElement_p->nodes[locnum[i] - 1];
  }
  //set vertexNumbers
  const std::vector<number_t>& svn = parent->geomRefElement()->sideOfSideVertexNumbers()[sos - 1];
  for(number_t i = 0; i < svn.size(); i++) { meshElement_p->vertexNumbers[i] = parent->meshElement_p->vertexNumber(svn[i]); }
  //update measures
  meshElement_p->computeMeasures();
  //update linearMap flag
  if(!meshElement_p->linearMap) meshElement_p->linearMap=meshElement_p->checkLinearMap();

  trace_p->pop();
  return meshElement_p;
}

//---------------------------------------------------------------------------
//  complex accessors
//---------------------------------------------------------------------------
//return i-th parent, return 0 if not a side element
const GeomElement* GeomElement::parent(number_t i) const
{
  if(i < parentSides_.size()) { return parentSides_[i].first; }
  else { return 0; }
}

//return i-th parent and its side as a GeoNumPair, (0,0) if not a side element
const GeoNumPair GeomElement::parentSide(number_t i) const
{
  if(i < parentSides_.size()) { return parentSides_[i]; }
  else { return GeoNumPair(static_cast<GeomElement*>(0), 0); }
}

//return first parent in a domain - not optimized
const GeomElement* GeomElement::parentInDomain(const MeshDomain* mdom) const
{
  if(mdom==0) return parent();   //no domain return first parent
  std::vector<GeoNumPair>::const_iterator itp=parentSides_.begin();
  for(; itp!=parentSides_.end(); ++itp)
    {
      const GeomElement* gelt=itp->first;
      std::vector<GeomElement*>::const_iterator itg=mdom->geomElements.begin(), itge=mdom->geomElements.end();
      while(itg!=itge)  //loop on domain elements - may be expansive
        {
          if(*itg==gelt) return gelt;
          itg++;
        }
    }
  return 0;
}

//return meshElement_p if not a null pointer
const MeshElement* GeomElement::meshElement() const
{
  //if (meshElement_p == 0) { error("geoelt_not_meshelement"); }
  return meshElement_p;
}

//return meshElement_p if not a null pointer (writable)
MeshElement* GeomElement::meshElement()
{
  //if (meshElement_p == 0) { error("geoelt_not_meshelement"); }
  return meshElement_p;
}

//return space dimension of the element
dimen_t GeomElement::elementDim() const
{
  if(parentSides_.size() > 0 && parentSides_[0].first->elementDim()==1) return 0;     //case of a 0 dim element !!!
  return geomRefElement()->dim();
}

//return space dimension, i.e nodes dimension
dimen_t GeomElement::spaceDim() const
{
  return mesh_p->spaceDim();
}

//return reference element if i=0 else refelement of side i
const RefElement* GeomElement::refElement(number_t i) const
{
  if(meshElement_p != 0) { return meshElement_p->refElement(i); }
  // side element with no MeshElement structure,
  return parentSides_[0].first->refElement(parentSides_[0].second)->refElement(i);
}

// returns number of vertices of element
number_t GeomElement::numberOfVertices() const
{
  if(meshElement_p != 0) { return meshElement_p->vertexNumbers.size(); }
  if(elementDim() < 1) { return 1; }   // one vertex for a point element
  return geomRefElement()->nbVertices();
}

// returns number of nodes of element
number_t GeomElement::numberOfNodes() const
{
  if(meshElement_p != 0) { return meshElement_p->nodes.size(); }
  if(elementDim() < 1) { return 1; }   //1 node for a point element
  return refElement()->nbDofs();     //nb of dofs is the number of nodes for a Lagrange element
}

// returns number of sides (faces in 3D, edges in 2D)
number_t GeomElement::numberOfSides() const
{
  if(meshElement_p != 0) {return meshElement_p->numberOfSides(); }
  if(elementDim() < 1) { return 0; }   //0 sides for a point element
  return geomRefElement()->nbSides();
}

// returns number of side of sides (edges in 3D)
number_t GeomElement::numberOfSideOfSides() const
{
  if(meshElement_p != 0) { return meshElement_p->numberOfSideOfSides(); }
  if(elementDim() < 1) { return 0; }   //0 sides of sides for a point element
  return geomRefElement()->nbSideOfSides();
}

// return the shape type of element (i=0) and its sides (_segment, _triangle, ...)
ShapeType GeomElement::shapeType(number_t i) const
{
  if(meshElement_p != 0) { return meshElement_p->shapeType(i); }
  return geomRefElement()->shapeType(i);
}

// return the global number of vertex i (i=1,...)
number_t GeomElement::vertexNumber(number_t i) const
{
  if(meshElement_p != 0) { return meshElement_p->vertexNumber(i); }
  // is a side element, go to parent
  GeomElement* parent = parentSides_[0].first;
  number_t sideNum = parentSides_[0].second;
  if(parent->meshElement_p != 0)
    {
      return parent->meshElement_p->vertexNumber(parent->meshElement_p->geomRefElement()->sideVertexNumbers()[sideNum - 1][i - 1]);
    }

  // is a side of side element, go to grandparent
  number_t sideNumgp = parent->parentSides_[0].second;
  parent = parent->parentSides_[0].first;
  if(parent->meshElement_p == 0) { error("geoelt_sosos", "GeomElement::vertexNumber"); }
  number_t s = std::abs(parent->meshElement_p->geomRefElement()->sideOfSideNumber(sideNum, sideNumgp)); //sideofside number
  return parent->meshElement_p->vertexNumber(parent->meshElement_p->geomRefElement()->sideOfSideVertexNumbers()[s - 1][i - 1]);
}

// return the global number of node i (i=1,...)
number_t GeomElement::nodeNumber(number_t i) const
{
  if(meshElement_p != 0) { return meshElement_p->nodeNumbers[i-1]; }
  // is a side element, go to parent
  GeomElement* parent = parentSides_[0].first;
  number_t sideNum = parentSides_[0].second;
  if(parent->meshElement_p != 0)
    {
      return parent->meshElement_p->nodeNumber(parent->meshElement_p->refElement()->sideDofNumbers_[sideNum - 1][i - 1]);
    }

  // is a side of side element, go to grandparent
  number_t sideNumgp = parent->parentSides_[0].second;
  parent = parent->parentSides_[0].first;
  if(parent->meshElement_p == 0) { error("geoelt_sosos", "GeomElement::nodeNumber"); }
  number_t s = std::abs(parent->meshElement_p->geomRefElement()->sideOfSideNumber(sideNum, sideNumgp)); //sideofside number
  return parent->meshElement_p->nodeNumber(parent->meshElement_p->refElement()->sideOfSideDofNumbers_[s - 1][i - 1]);
}

// return the global numbers of vertices (s=0) or the global numbers of vertices on side s
std::vector<number_t> GeomElement::vertexNumbers(number_t s) const
{
  if(s == 0)   //global numbers of vertices
    {
      if(meshElement_p != 0) { return meshElement_p->vertexNumbers; }
      std::vector<number_t> vnum(numberOfVertices());
      for(number_t i = 0; i < numberOfVertices(); i++) { vnum[i] = vertexNumber(i + 1); }
      return vnum;
    }
  // global numbers of vertices on side s
  std::vector<number_t> vnum = geomRefElement()->sideVertexNumbers()[s - 1];
  for(number_t i = 0; i < vnum.size(); i++) { vnum[i] = vertexNumber(vnum[i]); }
  return vnum;
}

// return the global numbers of nodes (s=0) or the global numbers of nodes on side s
std::vector<number_t> GeomElement::nodeNumbers(number_t s) const
{
  if(s == 0)   //global numbers of nodes
    {
      if(meshElement_p != 0) { return meshElement_p->nodeNumbers;}
      std::vector<number_t> vnum(numberOfNodes());
      for(number_t i = 0; i < numberOfNodes(); i++) {vnum[i] = nodeNumber(i + 1);}
      return vnum;
    }
  // global numbers of vertices on side s
  std::vector<number_t> vnum = refElement()->sideDofNumbers_[s - 1];
  for(number_t i = 0; i < vnum.size(); i++) { vnum[i] = nodeNumber(vnum[i]); }
  return vnum;
}

//utility to encode element or side element : ascending order vertices global numbers as string
string_t GeomElement::encodeElement(number_t s) const
{
  string_t key = "";
  std::vector<number_t> vnum = vertexNumbers(s);
  std::sort(vnum.begin(), vnum.end());
  for(number_t i = 0; i < vnum.size(); i++) { key += tostring(vnum[i]) + " "; }
  return key;
}

//utility to encode element or side element : ascending order nodes global numbers as string
string_t GeomElement::encodeElementNodes(number_t s) const
{
  string_t key = "";
  std::vector<number_t> vnum = nodeNumbers(s);
  std::sort(vnum.begin(), vnum.end());
  for(number_t i = 0; i < vnum.size(); i++) { key += tostring(vnum[i]) + " "; }
  return key;
}

//! return measure of element if side_no = 0, measure of side number side_no > 0
real_t GeomElement::measure(const number_t sideNo) const
{
  if(meshElement_p!=0) return meshElement_p->measure(sideNo);
  if(sideNo==0) return parentSides_[0].first->measure(parentSides_[0].second);
  buildSideMeshElement();
  return meshElement_p->measure(sideNo);
}

//! characteristic size of element (computed if not available)
real_t GeomElement::size() const
{
  if(meshElement_p!=0) return meshElement_p->size;
  buildSideMeshElement();
  return meshElement_p->size;
}

//! centroid of element (computed if not available)
Point GeomElement::center() const
{
    if(meshElement_p==0) buildSideMeshElement();
    return meshElement_p->center();
}


//test if a point belongs to current geom element (use MeshElement::contains)
bool GeomElement::contains(const std::vector<real_t>& p)
{
  if(meshElement_p==0) buildSideMeshElement();
  return meshElement_p->contains(p);
}

//project a point to current geom element (use MeshElement::projection), h is the distance to element
Point GeomElement::projection(const std::vector<real_t>& p, real_t& h)
{
  if(meshElement_p==0) buildSideMeshElement();
  return meshElement_p->projection(p,h);
}

//adjacent element by the side s, do not work for side elements!
//Mesh::sides_ must have been constructed before
//return element and its side number if exists, else return (0,0)
GeoNumPair GeomElement::elementSharingSide(number_t s) const
{
  if(isSideElement()) { error("geoelt_isside", "elementSharingSide"); }
  if(mesh_p->sides().size() == 0) { error("geoelt_sidesnotbuilt", "elementOnSide", "sides"); }
  const std::vector<GeoNumPair>& parSide = mesh_p->sides()[meshElement()->sideNumbers[s - 1] - 1]->parentSides();
  if(parSide.size() > 0 && parSide[0].first != this) { return parSide[0]; }
  if(parSide.size() > 1 && parSide[1].first != this) { return parSide[1]; }
  return GeoNumPair(static_cast<GeomElement*>(0), 0);
}

//adjacent element by the side of side s, only in 3D, do not work for side elements!
//when obs=true, list of elements sharing only the side of side s
//Mesh::sideOfSides_ must have been constructed before
//return element and its side number if exists, else return void
std::vector<GeoNumPair> GeomElement::elementsSharingSideOfSide(number_t s, bool obs) const
{
  if(isSideElement()) { error("geoelt_isside", "elementSharingSideOfSide"); }
  if(mesh_p->sideOfSides().size() == 0) { error("geoelt_sidesnotbuilt", "elementSharingSideOfSide", "sideOfSides"); }
  std::map<GeomElement*, number_t> elts;

  const std::vector<GeoNumPair>& par = mesh_p->sideOfSides()[meshElement()->sideOfSideNumbers[s - 1] - 1]->parentSides();
  std::vector<GeoNumPair>::const_iterator itpar, itgpar;
  for(itpar = par.begin(); itpar != par.end(); itpar++)
    {
      GeomElement* pelt = itpar->first;
      if(!pelt->isSideElement())
        {
          where("GeomElement::elementsSharingSideOfSide");
          error("geoelt_notside");         // parent element is not a side element
        }
      for(itgpar = pelt->parentSides().begin(); itgpar != pelt->parentSides().end(); itgpar++)   // add grandparents
        {
          if(itgpar->first != this && elts.find(itgpar->first) == elts.end()) { elts[itgpar->first] = itgpar->second; }   // add new grand parent element
        }
    }

  // if obs=true remove adjacent elements by side
  if(obs)
    {
      for(number_t s = 1; s <= numberOfSides(); s++) { elts.erase(elementSharingSide(s).first); }
    }

  // copy map in result vector
  std::map<GeomElement*, number_t>::iterator itm;
  std::vector<GeoNumPair> res(elts.size());
  std::vector<GeoNumPair>::iterator itres = res.begin();
  for(itm = elts.begin(); itm != elts.end(); itm++, itres++) { *itres = GeoNumPair(itm->first, itm->second); }
  return res;
}

//list of adjacent elements by the vertex v, do not work for side elements!
//when obv=true, list of elements sharing only the vertex v
//Mesh::vertexElements_ must have been constructed before
//return a vector of GeoNumPair, may be a void list if there is no adjacent element
std::vector<GeoNumPair> GeomElement::elementsSharingVertex(number_t v, bool obv) const
{
  if(isSideElement()) { error("geoelt_isside", "elementsOnVertex"); }
  if(mesh_p->vertexElements().size() == 0) { error("geoelt_sidesnotbuilt", "elementsOnVertex", "vertexElement"); }
  const std::vector<GeoNumPair>& elts = mesh_p->vertexElements()[vertexNumber(v) - 1]; //list of all elements sharing vertex v
  std::vector<GeoNumPair> res(elts.size() - 1);
  std::vector<GeoNumPair>::const_iterator itelts;
  std::vector<GeoNumPair>::iterator itres = res.begin();
  for(itelts = elts.begin(); itelts != elts.end(); itelts++)
    {
      if((*itelts).first != this) {*itres = *itelts; itres++;}
    }
  if(!obv) { return res; }

  // if obv=true remove adjacent elements by side or side by side
  std::map<GeomElement*, number_t> selts;
  std::vector<GeoNumPair>::iterator itg;
  for(number_t sos = 1; sos <= numberOfSideOfSides(); sos++)
    {
      std::vector<GeoNumPair> soselts = elementsSharingSideOfSide(sos); //elements sharing side of side sos
      for(itg = soselts.begin(); itg != soselts.end(); itg++) { selts[itg->first] = itg->second; }
    }
  std::vector<GeoNumPair> res2;
  for(itres = res.begin(); itres != res.end(); itres++)
    {
      if(selts.find(itres->first) == selts.end()) { res2.push_back(*itres); }
    }
  return res2;
}


// split mesh element in first order elements (P1) and return split elements in a list
std::vector<GeomElement*> GeomElement::splitP1() const
{
  std::vector<GeomElement*> listel;
  if(meshElement_p == 0)
    {
      warning("free_warning", "try to split a side element in GeomElement::split, nothing is done");
      return listel;
    }

  if(meshElement_p->order() == 1 && meshElement_p->isSimplex())    //no new element, copy of pointer
    {
      GeomElement* geoelt= new GeomElement(*this);
      listel.push_back(geoelt);
      return listel;
    }
  // split element
  std::vector<std::vector<number_t> > subeltnum=meshElement_p->refElement()->splitP1();
  std::vector<std::vector<number_t> >::iterator ite;
  Interpolation* intp=findInterpolation(_Lagrange, _standard, 1,_H1);
  RefElement* refelt=0;
  if(meshElement_p->elementDim()==1) refelt=findRefElement(_segment,intp);
  if(meshElement_p->elementDim()==2) refelt=findRefElement(_triangle,intp);
  if(meshElement_p->elementDim()==3) refelt=findRefElement(_tetrahedron,intp);
  for(ite=subeltnum.begin(); ite!=subeltnum.end(); ite++)
    {
      GeomElement* geoelt= new GeomElement(0, refelt,mesh_p->spaceDim(), 0);  //update later mesh pointer and number
      listel.push_back(geoelt);
      //update meshElement
      MeshElement* melt=geoelt->meshElement_p;
      melt->firstOrderParent_= meshElement_p;
      number_t nbnodes=ite->size();
      melt->nodeNumbers.resize(nbnodes);
      melt->nodes.resize(nbnodes);
      std::vector<number_t>::iterator itn, itk=melt->nodeNumbers.begin();
      std::vector<Point*>::iterator itp=melt->nodes.begin();
      for(std::vector<number_t>::iterator itn=ite->begin(); itn!=ite->end(); itn++,itk++,itp++)
        {
          *itk=meshElement_p->nodeNumbers[*itn-1];
          *itp=meshElement_p->nodes[*itn-1];
        }
      melt->vertexNumbers=melt->nodeNumbers;
      melt->computeMeasures();
      melt->computeOrientation();
    }
  return listel;
}

/*
// split mesh element in first order elements of same shape and return split elements in a list
std::vector<GeomElement*> GeomElement::splitO1() const
{
  std::vector<GeomElement*> listel;
  if(meshElement_p == 0)
    {
      warning("free_warning", "try to split a side element in GeomElement::split, nothing is done");
      return listel;
    }
  if(meshElement_p->order() == 1)    //no new element, copy of pointer
    {
      GeomElement* geoelt= new GeomElement(*this);
      listel.push_back(geoelt);
      return listel;
    }
  // split element
  std::vector<std::pair<ShapeType,std::vector<number_t> > > subeltnum = meshElement_p->refElement()->splitO1();
  std::vector<std::pair<ShapeType,std::vector<number_t> > >::iterator ite;
  Interpolation* intp=findInterpolation(_Lagrange, _standard, 1,_H1);
  RefElement* refelt=0;
  for(ite=subeltnum.begin(); ite!=subeltnum.end(); ite++)
    {
      refelt=findRefElement(ite->first,intp);
      GeomElement* geoelt= new GeomElement(0, refelt,mesh_p->spaceDim(), 0);  //update later mesh pointer and number
      listel.push_back(geoelt);
      //update meshElement
      MeshElement* melt=geoelt->meshElement_p;
      melt->firstOrderParent_= meshElement_p;
      number_t nbnodes=(ite->second).size();
      melt->nodeNumbers.resize(nbnodes);
      melt->nodes.resize(nbnodes);
      std::vector<number_t>::iterator itn, itk=melt->nodeNumbers.begin();
      std::vector<Point*>::iterator itp=melt->nodes.begin();
      for(std::vector<number_t>::iterator itn=(ite->second).begin(); itn!=(ite->second).end(); itn++,itk++,itp++)
        {
          *itk=meshElement_p->nodeNumbers[*itn-1];
          *itp=meshElement_p->nodes[*itn-1];
        }
      melt->vertexNumbers=melt->nodeNumbers;
      melt->computeMeasures();
      melt->computeOrientation();
    }
  return listel;
}
*/

//! return outward normal vector on side
std::vector<real_t> GeomElement::normalVector(number_t s) const
{
  if(meshElement_p==0)  buildSideMeshElement();
  return  meshElement_p->normalVector(s);
}

//! return outward normal vector on side
std::vector<real_t> GeomElement::tangentVector(number_t s) const
{
  if(meshElement_p==0)  buildSideMeshElement();
  return  meshElement_p->tangentVector(s);
}

//---------------------------------------------------------------------------
//print utilities
//---------------------------------------------------------------------------
void GeomElement::print(std::ostream& os) const
{
  if(theVerboseLevel == 0) { return; }
  if(meshElement_p != 0)
    {
      os << words("geometric element") << " " << number_ <<" (material "<<materialId<<", color "<<color<<") : " << (*meshElement_p);
      if(theVerboseLevel < 2) return;

      // print adjacent elements if possible
      if(mesh_p != 0 && mesh_p->sides().size() > 0 && !isSideElement())
        {
          os << "\n   " << words("adjacent elements") << " " << words("by side") << " :";
          bool first=true;
          for(number_t s = 1; s <= numberOfSides(); s++)
            {
              GeoNumPair g = elementSharingSide(s);
              if(g.first != 0)
                {
                  if(!first)  os << ",";
                  else  first = false;
                  os << " " << s << " -> " << g.first->number();
                }
            }
        }
      if(mesh_p != 0 && mesh_p->sideOfSides().size() > 0 && !isSideElement())
        {
          os << "\n   " << words("adjacent elements") << " " << words("by side of side") << " :";
          bool first = true;
          for(number_t s = 1; s <= numberOfSideOfSides(); s++)
            {
              const std::vector<GeoNumPair> elts = elementsSharingSideOfSide(s);
              if(elts.size() > 0)
                {
                  if(!first) os << ",";
                  else first = false;
                  os << "  " << s << "->";
                  for(number_t e = 0; e < elts.size(); e++) { os << " " << elts[e].first->number(); }
                }
            }
        }
      if(mesh_p != 0 && mesh_p->vertexElements().size() > 0 && !isSideElement())
        {
          os << "\n   " << words("adjacent elements") << " " << words("by vertex") << " :";
          bool first = true;
          for(number_t v = 1; v <= numberOfVertices(); v++)
            {
              const std::vector<GeoNumPair> elts = elementsSharingVertex(v);
              if(elts.size() > 0)
                {
                  if(!first) os << ",";
                  else first = false;
                  os << " " << v << "->";
                  for(number_t e = 0; e < elts.size(); e++) { os << " " << elts[e].first->number(); }
                }
            }
        }
    }
  if(isSideElement())
    {
      if (meshElement_p != 0) os << std::endl;
      os << words("geometric side element") << " " << number_ << " : ";
      for(number_t i = 0; i < parentSides_.size(); i++)
        {
          if (i > 0) { os << ", ";}
          os << words("side") << " " << parentSides_[i].second
             << " " << words("of element") << " " << parentSides_[i].first->number_;
          if(theVerboseLevel > 2) { os << "\n  " << (*parentSides_[i].first); }
        }
    }
}

//---------------------------------------------------------------------------
// external functions involving GeomElement
//---------------------------------------------------------------------------
//overload operators < and == to sort elements using STL vector sort method
bool operator== (const GeomElement& x, const GeomElement& y)
{ return (x.number() == y.number()); }

bool operator< (const GeomElement& x, const GeomElement& y)
{ return x.number() < y.number(); }

//output stream operator <<
std::ostream& operator<< (std::ostream& os, const GeomElement& geo)
{
  geo.print(os);
  return os;
}

/*! create an index of all sides of a list of elements
   very similar to buildSides function except the fact that no new side element are created
   do not clear current sideIndex !
*/
void createSideIndex(const std::vector<GeomElement*>& elements,
                     std::map<string_t,std::vector<GeoNumPair> >& sideIndex)
{
  trace_p->push("createSideIndex of mesh");
  std::vector<GeomElement*>::const_iterator itel=elements.begin();
  for(; itel != elements.end(); itel++)
    {
      MeshElement* melt = (*itel)->meshElement();
      std::map<string_t,std::vector<GeoNumPair> >::iterator itsn;
      if(melt!=0 && melt->elementDim() == melt->spaceDim())  //plain element
        if(!(*itel)->isSideElement())  //not side element
          {
            for(number_t s = 0; s < (*itel)->numberOfSides(); s++)  // loop on element sides
              {
                string_t key = (*itel)->encodeElement(s + 1); // encode vertex numbers of side s
                itsn = sideIndex.find(key);                   // search side in sideIndex map
                if(itsn == sideIndex.end()) sideIndex[key]=std::vector<GeoNumPair>(1,GeoNumPair(*itel,s+1));
                else sideIndex[key].push_back(GeoNumPair(*itel,s+1));
              } // end for (number_t s
          }// end if(!itel
    } //end for (itel
  trace_p->pop();
}

/*! create an index of all side elements in given side element list, do not clear current sideEltIndex !
*/
void createSideEltIndex(const std::vector<GeomElement*>& elements, std::map<string_t, GeomElement* >& sideEltIndex)
{
  trace_p->push("Mesh::createSideEltIndex");
  std::map<string_t, GeomElement* >::iterator its;
  std::vector<GeomElement*>::const_iterator itel=elements.begin();
  for(; itel != elements.end(); itel++)
    {
      if((*itel)->isSideElement())  //side element
        {
          string_t key = (*itel)->encodeElement(); // encode vertex numbers of side s
          its = sideEltIndex.find(key);            // search side in sideEltIndex map
          if(its == sideEltIndex.end()) sideEltIndex[key] = *itel; //insert sidelt in sideEltIndex map
        }// end if(!itel
    } //end for (itel
  trace_p->pop();
}

//===========================================================================
// MeshElement member functions and related external functions
//===========================================================================
//default constructor
MeshElement::MeshElement()
  : refElt_p(0), index_(0), spaceDim_(0) {orientation = 0; geomMapData_p = 0; firstOrderParent_=0; linearMap=true;}

//basic constructor
MeshElement::MeshElement(const RefElement* re, dimen_t sdim, number_t n)
  :  refElt_p(re), index_(n), spaceDim_(sdim)
{
  orientation = 0;
  geomMapData_p = 0;
  firstOrderParent_=0;
  // Create arrays for global numbering ...
  const GeomRefElement* gre_p = refElt_p->geomRefElem_p;
  nodes.resize(refElt_p->nbPts());                          // ... of node pointers
  nodeNumbers.resize(refElt_p->nbPts());                    // ... of node numbers
  vertexNumbers.resize(gre_p->nbVertices());                // ... of vertices
  sideNumbers.resize(gre_p->nbSides());                     // ... of sides (faces in 3D, edges in 2D)
  sideOfSideNumbers.resize(gre_p->nbSideOfSides());         // ... of faces
  measures.resize(1 + gre_p->nbSides(), 0.);                // ... of measures
  linearMap = (refElt_p->isSimplex()) && (refElt_p->order()<2);
}

//check if linear map is available, return true if it is
bool MeshElement::checkLinearMap() const
{
   if(nodes.size()==0) return false; //prevent undefined nodes definition
   std::vector<Point*>::const_iterator itn=nodes.begin();
   dimen_t dn=(*itn)->dim();
   if(dn==0) return false;   //prevent unfinished nodes definition, nodes vector allocated but not filled
   itn++;
   for(;itn!=nodes.end();++itn)
       if((*itn)->dim()!=dn) return false;  //prevent partial nodes definition, nodes vector allocated but not completely filled

   //build linear map f(x,y,z)=0+x*A+y*B+z*C
   std::vector<number_t> sns=refElt_p->geomRefElem_p->simplexNodes(); //node numbers (>=1) defining first simplex of ref element
   number_t d=sns.size()-1;
   Point O = *nodes[sns[0]-1];
   Point A = *nodes[sns[1]-1]-O;
   Point B,C;
   if(d>1) B=  *nodes[sns[2]-1]-O;
   if(d>2) C=  *nodes[sns[3]-1]-O;
   //test f(ref_node_i)==node_i
   real_t eps=characteristicSize()/1000;
   number_t i=0;
   std::vector<RefDof*>::const_iterator itd=refElt_p->refDofs.begin();
   for(;itd!=refElt_p->refDofs.end();++itd,++i)
     {
         std::vector<real_t> Q=(*itd)->point();
         Point P = O + Q[0]*A;
         if(d>1)  P += Q[1]*B;
         if(d>2)  P += Q[2]*C;
         if(dist(P,*nodes[i])>eps) return false;
     }
  return true;
}

// update node pointers vector from Mesh::nodes and nodeNumbers
// do not use if nodeNumbers has not to be built
void MeshElement::setNodes(std::vector<Point>& points)
{
  for(number_t i = 0; i < nodeNumbers.size(); i++) { nodes[i] = &points[nodeNumbers[i] - 1]; }
  //update linearMap flag
  if(!linearMap) linearMap=checkLinearMap();
}

// return the global numbers of vertices of element (s=0), of side s>0
std::vector<number_t> MeshElement::verticesNumbers(number_t s) const
{
  if(s == 0) { return vertexNumbers; }
  //global numbers of vertices on side s
  std::vector<number_t> vnum = geomRefElement()->sideVertexNumbers()[s - 1];
  for(number_t i = 0; i < vnum.size(); i++) { vnum[i] = vertexNumber(vnum[i]); }
  return vnum;
}

/*
--------------------------------------------------------------------------------
  set measure of elements and its side and set orientation (sign of the jacobioan)
--------------------------------------------------------------------------------
*/
void MeshElement::computeMeasures()
{
  computeMeasure();
  computeMeasureOfSides();
  centroid = center();
  size = characteristicSize();
}

//compute a characteristic size assuming measures of element and sides are up to date
//    for simplex : diameter of the circumscribed ball
//    other cases : max of length of segments joining two vertices
real_t MeshElement::characteristicSize() const
{
  number_t nbv=0;
  switch(shapeType())
    {
      case _point      :
      case _segment    : return measures[0];
      case _triangle   : return measures[1]* measures[2]* measures[3]/(2*measures[0]);  //circumscribed circle diameter
      case _tetrahedron: //circumscribed sphere diameter
        {
          real_t a=pointDistance(*nodes[0], *nodes[1]), a1=pointDistance(*nodes[2], *nodes[3]);
          real_t b=pointDistance(*nodes[1], *nodes[2]), b1=pointDistance(*nodes[0], *nodes[3]);
          real_t c=pointDistance(*nodes[0], *nodes[2]), c1=pointDistance(*nodes[1], *nodes[3]);
          real_t aa1=a*a1, bb1=b*b1, cc1=c*c1, p=(aa1+bb1+cc1)/2;
          return std::sqrt(p*(p-aa1)*(p-bb1)*(p-cc1))/(3*measures[0]);
        }
      case _quadrangle : nbv=4; break;
      case _prism      : nbv=6; break;
      case _hexahedron : nbv=8; break;
      case _pyramid    : nbv=5; break;
      default:
        error("geoelt_noshapetype", shapeType(), "MeshElement::characteristicSize");
        break;
    }

  //max of all segments joining two vertices
  real_t h=0;
  for(number_t i=0; i<nbv; i++)
    for(number_t j=i+1; j<nbv; j++)
      h=std::max(h, pointDistance(*nodes[i], *nodes[j]));
  return h;
}

/*compute measure of elements (length, area or volume of implied assembly of straight simplices)
  Should be called only for Lagrange elements, assuming vertices are first points of nodes vector
*/
void MeshElement::computeMeasure()
{
  switch(shapeType())
    {
      case _point:
        measures[0] = 1.; // Dirac
        break;
      case _segment:
        measures[0] = pointDistance(*nodes[0], *nodes[1]);
        break;
      case _triangle:
        measures[0] = triangleArea(*nodes[0], *nodes[1], *nodes[2]);
        break;
      case _quadrangle:
        measures[0] = triangleArea(*nodes[0], *nodes[1], *nodes[2])
                      + triangleArea(*nodes[0], *nodes[2], *nodes[3]);
        break;
      case _tetrahedron:
        measures[0] = tetrahedronVolume(*nodes[0], *nodes[1], *nodes[2], *nodes[3]);
        break;
      case _prism:
        measures[0] = tetrahedronVolume(*nodes[0], *nodes[1], *nodes[5], *nodes[2])
                      + tetrahedronVolume(*nodes[2], *nodes[3], *nodes[5], *nodes[4])
                      + tetrahedronVolume(*nodes[0], *nodes[2], *nodes[4], *nodes[3]);
        break;
      case _hexahedron:
        measures[0] = tetrahedronVolume(*nodes[0], *nodes[2], *nodes[7], *nodes[3])
                      + tetrahedronVolume(*nodes[3], *nodes[4], *nodes[7], *nodes[6])
                      + tetrahedronVolume(*nodes[0], *nodes[3], *nodes[6], *nodes[4])
                      + tetrahedronVolume(*nodes[2], *nodes[0], *nodes[5], *nodes[1])
                      + tetrahedronVolume(*nodes[1], *nodes[6], *nodes[5], *nodes[4])
                      + tetrahedronVolume(*nodes[2], *nodes[1], *nodes[4], *nodes[6]);
        break;
      case _pyramid:
        measures[0] = tetrahedronVolume(*nodes[0], *nodes[1], *nodes[2], *nodes[4])
                      + tetrahedronVolume(*nodes[1], *nodes[2], *nodes[3], *nodes[4]);
        break;
      default:
        error("geoelt_noshapetype", shapeType(), "MeshElement::computeMeasure");
        break;
    }
}

/*compute measure of side of elements (length, area or volume of implied assembly of straight simplices)
  Should be called only for Lagrange elements, assuming vertices are first points of nodes vector
*/
void MeshElement::computeMeasureOfSides()
{
  for(number_t s = 1; s <= numberOfSides(); s++)
    {
      //vertex local numbers on side s
      const std::vector<number_t>& vn = geomRefElement()->sideVertexNumbers()[s - 1];

      switch(shapeType(s))
        {
          case _point :
            measures[s] = 0.;
            break;
          case _segment:
            measures[s] = pointDistance(*nodes[vn[0] - 1], *nodes[vn[1] - 1]);
            break;
          case _triangle:
            measures[s] = triangleArea(*nodes[vn[0] - 1], *nodes[vn[1] - 1], *nodes[vn[2] - 1]);
            break;
          case _quadrangle:
            measures[s] = triangleArea(*nodes[vn[0] - 1], *nodes[vn[1] - 1], *nodes[vn[2] - 1])
                          + triangleArea(*nodes[vn[2] - 1], *nodes[vn[3] - 1], *nodes[vn[0] - 1]);
            break;
          default: error("geoelt_noshapetype", shapeType(s), "MeshElement::computeMeasureOfSides");
            break;
        }
    }
}

/*
--------------------------------------------------------------------------------
 compute element orientation
 if necessary, a GeomMapData object is created but it is deleted at the end
 if a GeomMapData object already exists, it is not deleted
 if jacobianMatrix does not exist it is created from element centroid and then cleared
 NOTE : when domain is a manifold (elementDim()< spaceDim_)
      the orientation is computed using a general process, see setOrientationManifold, setOrientationBoundary
--------------------------------------------------------------------------------
*/
void MeshElement::computeOrientation()
{
  if(geomMapData_p == 0)  {geomMapData_p = new GeomMapData(this);}
  std::vector<real_t>::const_iterator itc = geomRefElement()->centroid();
  std::vector<real_t> center(itc, itc + geomRefElement()->dim());
  geomMapData_p->computeJacobianMatrix(center);
  orientation=0;
  if(elementDim()==spaceDim_)
    orientation = geomMapData_p->computeJacobianDeterminant() > 0 ? 1 : -1;
}

/*! return the center of element based on vertices : sum_i vertex(i) / nb_vertex */
Point MeshElement::center() const
{
  Point G=*nodes[0];
  number_t n=verticesNumbers().size();
  for(number_t i=1; i<n; i++) G+= *nodes[i];
  return G/=n;
}

/*! return outward normal vector on side s (edge in 2D, face in 3D)
   normal is computed from vertices of side, exact only for 1 order element
   Is not normalized !
   if s=0 return the normal to element using the jacobian and the orientation
*/
std::vector<real_t> MeshElement::normalVector(number_t s) const
{
  dimen_t edim = elementDim();
  if(s>0)  //normal vector to a side
  {
      switch(edim)
    {
      case 1 : //for consistancy
        {
          return std::vector<real_t>(1,1.);
        }

      case 2 : //general formula that not assumes that element "lives" in R2
        {
          number_t ss=s+1;
          if(s==numberOfSides()) ss=1;
          Point S12=vertexOnSide(2,s)-vertexOnSide(1,s), S13=vertexOnSide(2,ss)-vertexOnSide(1,s);
          return (dot(S12,S13)/dot(S12,S12)*S12-S13);
        }

      case 3 : //use cross product
        {
          Point S12=vertexOnSide(2,s)-vertexOnSide(1,s), S23=vertexOnSide(3,s)-vertexOnSide(2,s);
          Point T=vertexOppositeToSide(s)-vertexOnSide(1,s);
          Point n=cross3D(S12,S23);
          if(dot(T,n) > 0) return -n;
          return n;
        }

      default :
        where("MeshElement::normalVector(Number)");
        error("index_out_of_range","dimension",1,3);
    }
  }

  //normal vector to element (edim should be <=2)
  if(edim==spaceDim()) // it is not a manifold, return as normal vector (0,1) for a 1D element or (0,0,1) for a 2D element
  {
      std::vector<real_t> n(edim+1,0.);
      n[edim]=1.;
      return n;
  }
  //case of a manifold, use jacobian and orientation
  if(geomMapData_p==0) geomMapData_p = new GeomMapData(this);
  if(geomMapData_p->normalVector.size()==0)
  {
      if(geomMapData_p->jacobianMatrix.size()==0) geomMapData_p->computeJacobianMatrix(center());
      geomMapData_p->computeNormalVector();                       // compute unnormalized normal vector
      if(orientation==0)
      {
        where("MeshElement::normalVector(Number)");
        error("is_null",words("orientation"));
      }
      geomMapData_p->normalize();                                 // normalize normal vector
  }
  return geomMapData_p->normalVector;
}

/*! return tangent vector on edge e (side in 2D, side on side in 3D)
   tangent vector  is computed from vertices of edge, exact only for 1 order element
*/
std::vector<real_t> MeshElement::tangentVector(number_t e) const
{
  switch(elementDim())
    {
      case 2 : return vertexOnSide(2,e)-vertexOnSide(1,e);
      case 3 : return vertexOnSideOfSide(2,e)-vertexOnSideOfSide(1,e);
      default :
        where("MeshElement::tangentVector(Number)");
        error("dim_not_in_range",2,3);
    }
  return std::vector<real_t>();  //fake return
}

//--------------------------------------------------------------------------------
/*! test if a point belongs to current mesh element
   if element is a simplex and is one order , we use the map from element to ref element
   else the element is split in P1 elements and test is performed on each P1 sub element
   The tolerance is chosen as size(element)/1000
*/
//--------------------------------------------------------------------------------
bool MeshElement::contains(const std::vector<real_t>& p)
{
  bool inelt = false;
  real_t tol = measure();
  if(tol==0.) {computeMeasure(); tol = measure();}
  tol=std::max(measure()/1000, theTolerance);
  dimen_t edim=elementDim();
  if(order()==1 && isSimplex())
    {
      if(p.size()==edim)   //use inverse map
        {
          GeomMapData gd(this,p);
          Point q=gd.geomMapInverse(p);
          return refElement()->geomRefElement()->contains(q, tol);
        }
      else //p.size > elementDim
        {
          real_t s=0;
          if(edim==1)  //use cross product S1S2 x S1P to check if P belongs to the line containing segment
            {
              Point p01=*nodes[1]-*nodes[0], p0=p-*nodes[0];
              s=norm2(crossProduct(p01,p0));
              if(s>tol) return false;
              real_t a=dot(p01,p0)/norm2(p01);
              return (a>-tol && a<1+tol);
            }
          if(edim==2)  //use mix product (S1S2 x S1S3).S1P to check if P belongs to the plane containing triangle
            {
              Point n10=*nodes[1]-*nodes[0], n20=*nodes[2]-*nodes[0], p0= p-*nodes[0];
              Point u=crossProduct(n10,n20);
              s=std::abs(dot(u,p0));
              if(s>tol) return false;
              //use barycentric coordinates
              real_t nu=dot(u,u);
              real_t a= dot(crossProduct(p0,n20),u)/ nu;
              if(a<-tol || a>1+tol) return false;
              real_t b= dot(crossProduct(n10,p0),u)/ nu;
              if(b<-tol || b>1+tol) return false;
              real_t c=1-a-b;
              if(c<-tol || c>1+tol) return false;
              return true;
            }
          //old method
//       GeomMapData gd(this,p);
//       Point q=gd.geomMapInverse(p);
//       return= refElement()->geomRefElement()->contains(q, tol);
        }
    }
  //split in P1 elements
  std::vector<MeshElement*> eltsP1=splitP1();
  std::vector<MeshElement*>::iterator ite=eltsP1.begin();
  for(; ite!=eltsP1.end(); ite++)
    {
      if(!inelt && (*ite)->contains(p)) inelt=true;
      delete *ite;
    }
  return inelt;
}

// projection of a point onto a geometric element
Point MeshElement::projection(const std::vector<real_t>& p, real_t& h) const
{
  dimen_t edim=elementDim();
  if(p.size()==edim)  //use reference element
    {
      GeomMapData gd(this,p);
      Point q=gd.geomMapInverse(p);
      Point r=refElement()->geomRefElement()->projection(q,h);
      return gd.geomMap(r);
    }

  // case p.size()>elementDim(), use projection on one order element
  // !!! has to be extended in future for element of order greater than 1 !!!
  switch(edim)
    {
      case 1 : //projection on a segment [S1 S2]
        return projectionOnSegment(p,*nodes[0],*nodes[1],h);
      case 2 :
        {
          if(isSimplex())  return projectionOnTriangle(p,*nodes[0],*nodes[1],*nodes[2],h);  //projection on a triangle [S1 S2 S3]
          else //split element in triangles
            {
              std::vector<MeshElement*> eltsP1=splitP1();
              std::vector<MeshElement*>::iterator ite=eltsP1.begin();
              h=theRealMax;
              real_t he;
              Point q, qe;
              for(; ite!=eltsP1.end(); ++ite)
                {
                  std::vector<Point*>& ns = (*ite)->nodes;
                  qe = projectionOnTriangle(p,*ns[0],*ns[1],*ns[2],he);
                  if(he<h) {h=he; q=qe;}
                  delete *ite;
                }
              return q;
            }
        }
      default :
        where("MeshElement::projection(Vector<Real>, Real)");
        error("index_out_of_range","dimension", 1, 2);
    }
  return Point();// dummy return
}

// split mesh element in first order elements (P1) and return split elements in a list
// it is a "clone" of  std::vector<GeomElement*> GeomElement::splitP1() const
std::vector<MeshElement*> MeshElement::splitP1() const
{
  std::vector<MeshElement*> listel;

  if(order() == 1 && isSimplex())    //no new element, copy of MeshElement
    {
      listel.push_back(new MeshElement(*this));
      return listel;
    }
  // split element
  std::vector<std::vector<number_t> > subeltnum=refElement()->splitP1();
  std::vector<std::vector<number_t> >::iterator ite;
  Interpolation* intp=findInterpolation(_Lagrange, _standard, 1,_H1);
  RefElement* refelt=0;
  if(elementDim()==1) refelt=findRefElement(_segment,intp);
  if(elementDim()==2) refelt=findRefElement(_triangle,intp);
  if(elementDim()==3) refelt=findRefElement(_tetrahedron,intp);
  for(ite=subeltnum.begin(); ite!=subeltnum.end(); ite++)
    {
      MeshElement* melt= new MeshElement(refelt,spaceDim_, 0);  //update later mesh pointer and number
      listel.push_back(melt);
      melt->firstOrderParent_= this;
      number_t nbnodes=ite->size();
      melt->nodeNumbers.resize(nbnodes);
      melt->nodes.resize(nbnodes);
      std::vector<number_t>::iterator itn, itk=melt->nodeNumbers.begin();
      std::vector<Point*>::iterator itp=melt->nodes.begin();
      for(std::vector<number_t>::iterator itn=ite->begin(); itn!=ite->end(); itn++,itk++,itp++)
        {
          *itk=nodeNumbers[*itn-1];
          *itp=nodes[*itn-1];
        }
      melt->vertexNumbers=melt->nodeNumbers;
      melt->computeMeasures();
      melt->computeOrientation();
    }
  return listel;
}

//--------------------------------------------------------------------------------
// prints Geom Element characteristics and point numbers to opened ostream or default print file
//--------------------------------------------------------------------------------

void MeshElement::print(std::ostream& os) const
{
  if(theVerboseLevel == 0) { return; }
  os << refElt_p->name() << ", orientation " << std::showpos << orientation << std::noshowpos<<", ";
  if(!linearMap) os<<" non";
  os<<" linear map";
  if(measures.size() > 0) { os << ", measure = " << measures[0]; }

  if(theVerboseLevel < 2) { return; }

  std::vector<number_t>::const_iterator it;
  os << "\n   " << words("nodes") << " : ";
  for(number_t n = 0; n < nodeNumbers.size(); n++)
    {
      os << nodeNumbers[n] << " ";
      if(theVerboseLevel > 2)
        {
          os << "-> " << (*nodes[n]);
          if(n != nodeNumbers.size() - 1) { os << ", "; }
        }
    }

  os << "\n   " << words("vertices") << " : ";
  if(vertexNumbers.size() > 0 && vertexNumbers[0] > 0)
    {
      for(number_t n = 0; n < vertexNumbers.size(); n++)
        {
          os << vertexNumbers[n] << " ";
          if(theVerboseLevel > 2)
            {
              os << "-> " << (*nodes[n]);
              if(n != vertexNumbers.size() - 1) { os << ", "; }
            }
        }
    }
  else { os << words("unset"); }

  if(theVerboseLevel < 5) { return; }
  if(measures.size() > 1)
    {
      os << "\n   " << words("measure of sides") << " =";
      for(dimen_t i = 0; i < numberOfSides(); i++) { os << " " << measures[i + 1]; }
    }

  os << "\n   " << words("sides") << " : ";
  if(sideNumbers.size() > 0 && sideNumbers[0] > 0)
    {
      for(number_t n = 0; n < sideNumbers.size(); n++) { os << sideNumbers[n] << " "; }
    }
  else { os << words("unset"); }

  os << "\n   " << words("sideOfSides") << " : ";
  if(sideOfSideNumbers.size() > 0 && sideOfSideNumbers[0] > 0)
    {
      for(it = sideOfSideNumbers.begin(); it != sideOfSideNumbers.end(); it++) { os << (*it) << " "; }
    }
  else { os << words("unset"); }

  if(theVerboseLevel < 10000)  return;
  os << "\n" << *(refElt_p);
}

std::ostream& operator<< (std::ostream& os, const MeshElement& me)
{
  me.print(os);
  return os;
}

/*! distance from two MeshElement, assuming not intersection and using first order, i.e vertices
    under these assumptions
        dist(e1,e2)= min{dist(p1,p2), p1 in e1, p2 in e2}
                   = min{dist(p1,p2), p1 any vertex of e1, p2 any vertex of e2}
*/
real_t distance(const MeshElement& elt1, const MeshElement& elt2)
{
  real_t d = theRealMax;
  for(number_t i=1; i<=elt1.numberOfVertices(); i++)
    {
      Point & vi = elt1.vertex(i);
      for(number_t j=1; j<=elt2.numberOfVertices(); j++)
        {
          d=std::min(d,vi.distance(elt2.vertex(j)));
          if(d<=theTolerance) return 0.;
        }
    }
  return d;
}


//--------------------------------------------------------------------------------
// area and volume of simplices
//--------------------------------------------------------------------------------
  real_t triangleArea(const Point& p1, const Point& p2, const Point& p3)
  { return 0.5 * norm(crossProduct(p2 - p1, p3 - p1)); }

  real_t tetrahedronVolume(const Point& p1, const Point& p2, const Point& p3, const Point& p4)
  {
    std::vector<real_t> v = crossProduct(p1 - p4, p1 - p3);
    std::vector<real_t> w = p2 - p3;
    return std::abs(std::inner_product(v.begin(), v.end(), w.begin(), 0.)) / 6.;
  }

//--------------------------------------------------------------------------------
// default GeomElement coloring rule
//--------------------------------------------------------------------------------
  /*!
    the default GeomElement coloring rule is the following
      let np the number of strictly positive values (v_i) on the n vertices of the geomelement
      if(np>n/2) then color =1 else color =0
      for a triangle : color is 1 if at least 2 vertex values are >0
      for a quadrangle  : color is 1 if at least 3 vertex values are >0
      for a tetrahedron  : color is 1 if at least 3 vertex values are >0
      for a hexahedron  : color is 1 if at least 5 vertex values are >0

      gelt : geom element
      val : real values on vertices

  */
  real_t defaultColoringRule(const GeomElement& gelt, const std::vector<real_t>& val)
  {
    number_t np=0, n=std::min(val.size(),gelt.numberOfVertices());
    std::vector<real_t>::const_iterator itv=val.begin();
    for(number_t i=0; i<n; ++i, ++itv)
      if(*itv>0) np++;
    if(np>n/2) return 1;
    return 0;
  }


} // end of namespace xlifepp