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

  \brief Implementation of xlifepp::GeomMapData functions
*/

#include "GeomMapData.hpp"
#include "GeomElement.hpp"

namespace xlifepp
{

//constructors
GeomMapData::GeomMapData(const MeshElement* melt, const Point& x)
  : geomElement_p(melt), currentPoint(x)
{
  elementDim = geomElement_p->elementDim();
  spaceDim = geomElement_p->spaceDim();
  jacobianDeterminant = 0;
  differentialElement=0;
  metricTensorDeterminant=0;
  jacobianMatrix.clear();           // to reset the size of vector to 0 (default 1x1)
  inverseJacobianMatrix.clear();    // to reset the size of vector to 0 (default 1x1)
  normalVector.clear();             // to reset the size of vector to 0 (default 1)
  extdata=0;
}

GeomMapData::GeomMapData(const MeshElement* melt, std::vector<real_t>::const_iterator itp)
  : geomElement_p(melt)
{
  elementDim = geomElement_p->elementDim();
  spaceDim = geomElement_p->spaceDim();
  currentPoint = Point(itp, elementDim);
  jacobianDeterminant=0;
  differentialElement=0;
  metricTensorDeterminant=0;
  jacobianMatrix.clear();           // to reset the size of vector to 0 (default 1x1)
  inverseJacobianMatrix.clear();    // to reset the size of vector to 0 (default 1x1)
  normalVector.clear();             // to reset the size of vector to 0 (default 1)
  extdata=0;
}

GeomMapData::GeomMapData(const MeshElement* melt)    //origin as point
  : geomElement_p(melt)
{
  elementDim = geomElement_p->elementDim();
  spaceDim = geomElement_p->spaceDim();
  currentPoint = Point(std::vector<real_t>(elementDim, 0.));
  jacobianDeterminant = 0;
  differentialElement=0;
  metricTensorDeterminant=0;
  jacobianMatrix.clear();           // to reset the size of vector to 0 (default 1x1)
  inverseJacobianMatrix.clear();    // to reset the size of vector to 0 (default 1x1)
  normalVector.clear();             // to reset the size of vector to 0 (default 1)
  extdata=0;
}

/*
--------------------------------------------------------------------------------
  mapping a point from reference element onto geometric element point
--------------------------------------------------------------------------------
*/
//mapping from ref. elt onto geom. elt. Element at a given point x
Point GeomMapData::geomMap(const std::vector<real_t>& x, number_t side)
{
  currentPoint = x;
  //geomElement_p->refElement()->computeShapeValues(x.begin(), false);
  ShapeValues shv = geomElement_p->refElement()->shapeValues;
  geomElement_p->refElement()->computeShapeValues(x.begin(), shv, false);
  return geomMap(shv, side);
}

Point GeomMapData::geomMap(std::vector<real_t>::const_iterator itx, number_t side)
{
  Point x(itx, geomElement_p->refElement()->dim());
  return geomMap(x,side);
}

Point GeomMapData::geomMap(number_t side)
{
  return geomMap(geomElement_p->refElement()->shapeValues, side);
}

//mapping from ref. elt onto geom. elt. Element (shape functions are given, thread safe)
// side = 0, transformation on whole element
// side > 0, transformation on side side of element of dimension reference_dim_
Point GeomMapData::geomMap(const ShapeValues& shv, number_t side)
{
  // spaceDim   : dimension of mapped point
  // elementDim : dimension of reference point
  Point pt(std::vector<real_t>(spaceDim, 0.));
  Point::iterator it_pt = pt.begin(); //first coordinate of pt
  std::vector<Point*>::const_iterator it_pb = geomElement_p->nodes.begin(), it_p; // iterator on geometric element point coordinates
  std::vector<real_t>::const_iterator
  it_wb = shv.w.begin(), // iterators on reference element shape functions at current reference point
  it_we = shv.w.end(), it_w;

  if (side == 0)     //element map
    {
      for(dimen_t i = 0; i < spaceDim; i++, it_pt++)  // loop on coordinates
        {
          it_p = it_pb;  //first node (dof)
          for(it_w = it_wb; it_w != it_we; it_w++, it_p++)  // loop on reference d.o.f (shape values)
            { *it_pt += (**it_p)[i] **it_w; }
        }
    }
  else  //side element map
    {
      std::vector<number_t>::const_iterator it_dofb = geomElement_p->refElement()->sideDofNumbers_[side - 1].begin(),
                                            it_dofe = geomElement_p->refElement()->sideDofNumbers_[side - 1].end(), it_dof;
      for (dimen_t i = 0; i < spaceDim; i++)   // loop on coordinates
        {
          it_w = it_wb;
          for (it_dof = it_dofb; it_dof != it_dofe; it_dof++) // loop on reference side d.o.f
            { *it_pt += (**(it_pb + *it_dof - 1))[i] **(it_w + *it_dof - 1); }
        } // end for(dimen_t
    } // end of else
  return pt;
}

/*!mapping from elt onto ref. elt. (inverse of GeomMap)
  mapping may be non linear : order >1 or not simplicial element (quadrangle, hexahedron)
  so to inverse it we have to use an iterative method (Newton for instance)
    first step : initialize Newton's method with first order map say P_0=(JF1)^-1(P-F1(0))
             if element is a 1 order simplicial element return P_0
    second step: Newton's iteration P_k+1=P_k - JF(P_k)^-1(P-F(Pk))
             stop if ||F(P_k+1)-X||< eps or ||P_k+1-P_k+1||< eps
             stop also if the number of iteration k>50 (error)
*/
Point GeomMapData::geomMapInverse(const std::vector<real_t>& p, real_t eps, number_t maxIter)
{
  //if (geomElement_p->isSimplex() && geomElement_p->order()==1)  //first order
  if (geomElement_p->linearMap)  // linear map
    {
      if (inverseJacobianMatrix.size()==0) //compute it
        {
          if(jacobianMatrix.size()==0)
            {
              currentPoint = Point(std::vector<real_t>(elementDim,0.));
              geomElement_p->refElement()->computeShapeValues(currentPoint.begin(), true);
              computeJacobianMatrix();
            }
          invertJacobianMatrix();
        }
      Point q=geomMap(Point(std::vector<real_t>(elementDim,0.)));
      Vector<real_t> qp = p - q;
      return (inverseJacobianMatrix * qp);
    }
  else  //non linear map to inverse, to be improved
    {
      //create first order map using first order element
      std::vector<MeshElement*> eltsP1=geomElement_p->splitP1();
      GeomMapData gmd1(*eltsP1.begin(),p);
      gmd1.computeJacobianMatrix();
      gmd1.invertJacobianMatrix();
      Point q=gmd1.geomMap(Point(std::vector<real_t>(p.size(),0.)));
      Vector<real_t> qp = p - q;
      if (geomElement_p->spaceDim()== geomElement_p->elementDim())
        {
          invertJacobianMatrix();
          q=inverseJacobianMatrix * qp;
        }
      else //non square system, use least square method
        {
          if (inverseJacobianMatrix.size()==0) //compute it
            {
              currentPoint = p;
              computeJacobianMatrix();
            }
          Matrix<real_t> Jt=transpose(jacobianMatrix);
          q = inverse(Jt*jacobianMatrix)*(Jt*qp);
        }

      //Newton algorithm
      bool cont=true;
      for (number_t k=0; k<maxIter && cont; k++)
        {
          qp = p - geomMap(q);
          if (norminfty(qp)>eps)
            {
              currentPoint = q;
              computeJacobianMatrix();
              Vector<real_t> t;
              if (geomElement_p->spaceDim()== geomElement_p->elementDim())
                {
                  invertJacobianMatrix();
                  t=inverseJacobianMatrix * qp;
                }
              else //non square system, use least square method
                {
                  Matrix<real_t> JtJ=transpose(jacobianMatrix) * jacobianMatrix;
                  t = inverse(JtJ)*transpose(jacobianMatrix)*qp;
                }
              if (norminfty(t)>eps) q+=Point(t);
              else cont = false;
            }
          else cont=false;
        }
      std::vector<MeshElement*>::iterator ite=eltsP1.begin();
      for (; ite!=eltsP1.end(); ite++) delete *ite;
      return q;
    }
  return p;
}

//Piola mapping from ref. elt onto geom. elt. Element (shape functions knowns)
Point GeomMapData::piolaMap(number_t side)
{
  error("not_yet_implemented", "GeomMapData::piolaMap");
  return Point();
}

/*! return the covariant Piola map matrix at current point : J^{-T}
    when mapping from 1d->2d or 2d->3d the matrix J^{-T} is replace by J*(J{T}*J)^{-1}{T}
    note : in that case invertJacobian() compute J*(J{T}*J)^{-1}
    */
Matrix<real_t> GeomMapData::covariantPiolaMap(const Point& P)
{
  if(P.size()!=0)
    {
      computeJacobianMatrix(P.begin());
      invertJacobianMatrix();
    }
  else //use current point
    {
      if(jacobianMatrix.size()==0) computeJacobianMatrix();
      if(inverseJacobianMatrix.size()==0) invertJacobianMatrix();
    }
  return transpose(inverseJacobianMatrix);
}

/*! return the contravariant Piola map matrix at current point : J/|J|
    when mapping from 1d->2d or 2d->3d the matrix J/|J| is replace by J/sqrt|Jt*J|
    note : in that case computeJacobianDeterminant() compute sqrt|Jt*J|
*/
Matrix<real_t> GeomMapData::contravariantPiolaMap(const Point& P)
{
  if(P.size()!=0) computeJacobianMatrix(P.begin());
  else if(jacobianMatrix.size()==0) computeJacobianMatrix();//use current point
  real_t j=computeJacobianDeterminant();
  return jacobianMatrix/j;
}

/*
--------------------------------------------------------------------------------
  computes normal vector at a (boundary) point of element
  assuming that jacobian matrix is up to date
  only when space dimension = element dimension + 1
  or        space dimension = element dimension + 2 (special process)
--------------------------------------------------------------------------------
*/
void GeomMapData::computeNormalVector()
{
  dimen_t sd=geomElement_p->spaceDim();   //space dimension
  normalVector.resize(sd);
  if (spaceDim == elementDim+1)
    {
      switch (spaceDim)
        {
        case 1:
          normalVector[0] = 1.;
          break;
        case 2:
          normalVector[0] = jacobianMatrix[1];
          normalVector[1] = -jacobianMatrix[0];
          break;
        case 3:
          normalVector[0] = jacobianMatrix[2] * jacobianMatrix[5] - jacobianMatrix[4] * jacobianMatrix[3];
          normalVector[1] = jacobianMatrix[4] * jacobianMatrix[1] - jacobianMatrix[0] * jacobianMatrix[5];
          normalVector[2] = jacobianMatrix[0] * jacobianMatrix[3] - jacobianMatrix[2] * jacobianMatrix[1];
          break;
        }
      return;
    }
  where("GeomMapData::computeNormalVector()");
  error("geoelt_space_mismatch_dims", elementDim+1, spaceDim);
}

/*
--------------------------------------------------------------------------------
  normalize normal vector into a outward unit normal vector
--------------------------------------------------------------------------------
*/
void GeomMapData::normalize()
{
  // normalize vector
  real_t norm = 0.;
  norm = inner_product(normalVector.begin(), normalVector.end(), normalVector.begin(), norm);
  if (norm < theZeroThreshold) { error("is_null", words("norm")); }
  norm = 1. / std::sqrt(norm);
  if (geomElement_p->orientation < 0) { norm *= -1.; }
  for (Vector<real_t>::iterator it = normalVector.begin(); it != normalVector.end(); it++)
    { *it *= norm; }
}

/*
--------------------------------------------------------------------------------
  compute oriented unit normal vector, assuming jacobian matrix is update
--------------------------------------------------------------------------------
*/
void GeomMapData::computeOutwardNormal()
{
  computeNormalVector();
  normalize();
}

/*
--------------------------------------------------------------------------------
  jacobianMatrix computes Jacobian Matrix at Reference Element point x
    side = 0, jacobian matrix on whole element
    side > 0, jacobian matrix on side of element of dimension reference_dim
    spDim  : number of rows of jacobian matrix    (dim of geometric element points)
    elDim :  number of columns of jacobian matrix (dim of reference element points)
--------------------------------------------------------------------------------
*/
void GeomMapData::computeJacobianMatrix(const std::vector<real_t>& x, number_t side)
{
  currentPoint = x;
  geomElement_p->refElement(side)->computeShapeValues(x.begin(), true);
  computeJacobianMatrix(side);
}

void GeomMapData::computeJacobianMatrix(std::vector<real_t>::const_iterator itx, number_t side)
{
  currentPoint = Point(itx, geomElement_p->elementDim());
  geomElement_p->refElement(side)->computeShapeValues(itx, true);
  computeJacobianMatrix(side);
}

// assuming that derivatives of shape functions at x are known and stored in RefElement.shapeValues
void GeomMapData::computeJacobianMatrix(number_t side)
{
  computeJacobianMatrix(geomElement_p->refElement(side)->shapeValues, side);
}

// assuming that derivatives of shape functions at x are given by shv
void GeomMapData::computeJacobianMatrix(const ShapeValues& shv, number_t side)
{
  if (geomElement_p->refElement()->shapeType() == _point)
    {
      jacobianMatrix.changesize(1,1,1.);
      return;
    }
  //initialization
  if (side > 0) { elementDim--; }
  if (jacobianMatrix.size() != dimen_t(spaceDim * elementDim)) jacobianMatrix.changesize(spaceDim,elementDim);
  //if(jacobianMatrix.size() != dimen_t(spaceDim * elementDim)) { jacobianMatrix = Matrix<real_t>(spaceDim, elementDim, 0.); }

  //iterator definitions
  Matrix<real_t>::iterator it_jmb = jacobianMatrix.begin(), it_jm = it_jmb;        // iterators on matrix entries
  std::vector<Point*>::const_iterator it_pb = geomElement_p->nodes.begin(), it_p;  // iterator on geometric element point coordinates
  std::vector< std::vector<real_t> >::const_iterator
  //it_svb = geomElement_p->refElement(side)->shapeValues.dw.begin(), it_sv; // iterators on reference element shape functions derivatives
  it_svb = shv.dw.begin(), it_sv;                                            // iterators on reference element shape functions derivatives
  std::vector<real_t>::const_iterator it_dw;

  if (side == 0)    //jacobian matrix on whole element
    {
      for (dimen_t i = 0; i < spaceDim; i++)    // loop on rows of jacobian matrix
        {
          it_sv=it_svb;
          for (dimen_t j=0; j<elementDim; j++, ++it_sv, ++it_jm ) // loop on columns of jacobian matrix
            {
              *it_jm = 0.;             // iterator on geometric element point coordinates
              it_p = it_pb;            // loop on reference d.o.f (associated shape function derivative)
              for (it_dw = it_sv->begin(); it_dw != it_sv->end(); ++it_dw, ++it_p)
                *it_jm += (**it_p)[i] **it_dw;
            }
        } // end for(dimen_t
    } // end of if(Side

  else  //jacobian matrix on side
    {
      //std::vector< std::vector<real_t> >::const_iterator it_sve = geomElement_p->refElement(side)->shapeValues.dw.end();
      std::vector< std::vector<real_t> >::const_iterator it_sve = shv.dw.end();
      //travel dofs on side
      std::vector<number_t>::const_iterator it_dofb = geomElement_p->refElement()->sideDofNumbers_[side - 1].begin(), it_dof;
      for (dimen_t i = 0; i < spaceDim; i++)   // loop on rows of jacobian matrix
        {
          for (it_sv = it_svb; it_sv != it_sve; it_sv++, it_jm++)    // loop on columns of jacobian matrix
            {
              *it_jm = 0.;
              it_dof = it_dofb;      //begin of dof numbers
              for (it_dw = it_sv->begin(); it_dw != it_sv->end(); it_dw++, it_dof++)
                {
                  it_p = it_pb + *it_dof - 1;     // node on side
                  *it_jm += (**it_p)[i] **it_dw;  // update jacobian coefficient
                }
            }
        } // end for(dimen_t
    } // end of else
}

/*
--------------------------------------------------------------------------------------
  jacobianDeterminant computes Jacobian determinant, jacobianMatrix has to be computed
--------------------------------------------------------------------------------------
*/
real_t GeomMapData::computeJacobianDeterminant()
{
  if (spaceDim==elementDim) //standard jacobian
    {
      switch (spaceDim)
        {
        case 1:
          jacobianDeterminant = jacobianMatrix[0];
          break;
        case 2:
          jacobianDeterminant = jacobianMatrix[0] * jacobianMatrix[3] - jacobianMatrix[1] * jacobianMatrix[2];
          break;
        case 3:
          jacobianDeterminant = jacobianMatrix[0] * (jacobianMatrix[4] * jacobianMatrix[8] - jacobianMatrix[5] * jacobianMatrix[7])
                                + jacobianMatrix[1] * (jacobianMatrix[5] * jacobianMatrix[6] - jacobianMatrix[3] * jacobianMatrix[8])
                                + jacobianMatrix[2] * (jacobianMatrix[3] * jacobianMatrix[7] - jacobianMatrix[4] * jacobianMatrix[6]);
          break;
        default:
          where("GeomMapData::computeJacobianDeterminant()");
          error("dim_not_in_range", 1, 3);
        }
      differentialElement = std::abs(jacobianDeterminant);
      return jacobianDeterminant;
    }

  //case eltDim < spaceDim : choose as determinant the sqrt(abs(det(Jt*J))), determinant of the metric tensor
  Matrix<real_t> JtJ = transpose(jacobianMatrix)*jacobianMatrix;
  switch (elementDim)
    {
    case 1:
      jacobianDeterminant = std::sqrt(std::abs(JtJ[0]));
      break;
    case 2:
      jacobianDeterminant = std::sqrt(std::abs(JtJ[0] * JtJ[3] - JtJ[1] * JtJ[2]));
      break;
    default:
      where("GeomMapData::computeJacobianDeterminant()");
      error("dim_not_in_range", 1, 2);
    }
  computeDifferentialElement();
  return jacobianDeterminant;
}

/*
--------------------------------------------------------------------------------
  invertJacobianMatrix inverts Jacobian Matrix !
  note that this function also sets the jacobian determinant and the differential element
  Jacobian may not be a square matrix (spdim x eldim matrix with spdim > eldim)
  in that case the pseudo inverse (Jt * J)^-1 * Jt (eldim x spdim matrix) is returned
  and the differential element is computed using computeDifferentialElement function
--------------------------------------------------------------------------------
*/
void GeomMapData::invertJacobianMatrix()
{
  dimen_t elDim = jacobianMatrix.numberOfColumns();
  if (jacobianMatrix.size() == 0) { error("geoelt_noJacobian"); } //compute jacobian matrix before
  real_t det = 0.;
  Matrix<real_t>& jm = jacobianMatrix, &invjm = inverseJacobianMatrix;
  if (spaceDim == elDim)
    {
      if (invjm.size() != dimen_t(spaceDim * spaceDim)) { invjm = Matrix<real_t>(spaceDim, spaceDim, 0.); }

      switch (spaceDim)
        {
        case 1:
          det = jm[0];
          if (std::abs(det) < theZeroThreshold) { error("geoelt_nulldet", det, theZeroThreshold); }
          invjm[0] = 1. / det;
          break;
        case 2:
          det = jm[0] * jm[3] - jm[1] * jm[2];
          if (std::abs(det) < theZeroThreshold) { error("geoelt_nulldet", det, theZeroThreshold); }
          invjm[0] = jm[3] / det;
          invjm[3] = jm[0] / det;
          invjm[1] = -jm[1] / det;
          invjm[2] = -jm[2] / det;
          break;
        case 3:
          // set 1st column to signed cofactors of 1st row
          invjm[0] = jm[4] * jm[8] - jm[5] * jm[7];
          invjm[3] = jm[5] * jm[6] - jm[3] * jm[8];
          invjm[6] = jm[3] * jm[7] - jm[4] * jm[6];
          // compute determinant
          det = jm[0] * invjm[0] + jm[1] * invjm[3] + jm[2] * invjm[6];
          // singular matrix ?
          if (std::abs(det) < theZeroThreshold) { error("geoelt_nulldet", det, theZeroThreshold); }
          // divide 1st column by determinant
          invjm[0] = invjm[0] / det;
          invjm[3] = invjm[3] / det;
          invjm[6] = invjm[6] / det;
          // set 2nd column to signed cofactors of 2nd row divided by determinant
          invjm[1] = (jm[2] * jm[7] - jm[1] * jm[8]) / det;
          invjm[4] = (jm[0] * jm[8] - jm[2] * jm[6]) / det;
          invjm[7] = (jm[1] * jm[6] - jm[0] * jm[7]) / det;
          // set 3rd column to signed cofactors of 3rd row  divided by determinant
          invjm[2] = (jm[1] * jm[5] - jm[2] * jm[4]) / det;
          invjm[5] = (jm[2] * jm[3] - jm[0] * jm[5]) / det;
          invjm[8] = (jm[0] * jm[4] - jm[1] * jm[3]) / det;
          break;
        }
      jacobianDeterminant = det;
      differentialElement = std::abs(det);
    }
  else  //non square jacobian, use pseudoinverse (Jt*J)^-1 * Jt (matrix eltdim x spdim)
    {
      Matrix<real_t> jmt=transpose(jm);
      Matrix<real_t> jmtjm=jmt*jm;
      inverseJacobianMatrix = inverse(jmtjm)*jmt;
      switch (elementDim)
        {
        case 1:
          jacobianDeterminant = std::sqrt(std::abs(jmtjm[0]));
          break;
        case 2:
          jacobianDeterminant = std::sqrt(std::abs(jmtjm[0] * jmtjm[3] - jmtjm[1] * jmtjm[2]));
          break;
        default:
          where("GeomMapData::invertJacobianDeterminant()");
          error("dim_not_in_range", 1, 2);
        }
      computeDifferentialElement();  //compute the differential element
    }
}

/*
--------------------------------------------------------------------------------
  differentialElement computation at the current point
  used in computation of integrals over element
--------------------------------------------------------------------------------
*/
//compute differential element assuming jacobian matrix and jacobian determinant are update
void GeomMapData::computeDifferentialElement()
{
  dimen_t elDim = jacobianMatrix.numberOfColumns(); //in case of side of element, number of columns of jacobian is the side element dimension

  if (spaceDim == elDim)            // 3D volume integrals, 2D surface integrals, 1D line integrals
    {
      computeJacobianDeterminant();
    }
  else if (spaceDim == elDim + 1)  // 3d surface integrals (3D element face), 2D line integrals (2D element edge)
    {
      computeNormalVector();    // compute outward normal vector
      differentialElement = std::sqrt(std::inner_product(normalVector.begin(), normalVector.end(),
                                      normalVector.begin(), 0.));
      normalize();              //normalize normal vector
    }
  else if (spaceDim == elDim + 2)   // 3d line integrals, compute length of tangent vector
    {
      differentialElement = std::sqrt(std::inner_product(jacobianMatrix.begin(), jacobianMatrix.end(),
                                      jacobianMatrix.begin(), 0.));
    }
}

//compute differential element updating jacobian matrix and jacobian determinant
real_t GeomMapData::diffElement()
{
  computeJacobianMatrix();
  if (spaceDim == elementDim)  //  Reference Element dimension is space dimension ('Volume Domain computation")
    { invertJacobianMatrix(); }
  else   //  Reference Element dimension is less than space dimension ('Side Domain computation")
    {
      computeDifferentialElement();
      if(elementDim == spaceDim - 1) { normalize(); }
    }
  return differentialElement;
}

//compute differential element on side, updating jacobian matrix and jacobian determinant
//and normalize the normal to side
real_t GeomMapData::diffElement(number_t side)
{
  computeJacobianMatrix(side);
  computeDifferentialElement();
  if(jacobianMatrix.numberOfColumns() == geomElement_p->spaceDim() - 1) { normalize(); }
  return differentialElement;
}

/*
--------------------------------------------------------------------------------
 compute metric tensor and its determinant for surface differential geometry
 metric tensor is only defined for spDim=3 and elDim=2 : mt=J^tJ (matrix 2x2)
--------------------------------------------------------------------------------
*/
void GeomMapData::computeMetricTensor()
{
  if (spaceDim != 3)  { error("3D_only","GeomMapData::computeMetricTensor"); }
  if (elementDim != 2)  { error("geoelt_2D_only","GeomMapData::computeMetricTensor"); }
  if (metricTensor.size() != 4) { metricTensor = Matrix<real_t>(2,2,0.); }
  std::vector<real_t>::const_iterator it_jac = jacobianMatrix.begin();
  std::vector<real_t>::iterator it_mt;
  for (it_mt = metricTensor.begin(); it_mt != metricTensor.end(); it_mt++) { *it_mt = 0.; }
  it_mt = metricTensor.begin();
  for (dimen_t i = 0; i < spaceDim; i++)
    {
      *it_mt += *it_jac **it_jac;
      *(it_mt + 1) += *(it_jac + 1) **it_jac++;
      *(it_mt + 3) += *it_jac **it_jac++;
    }
  *(it_mt + 2) = *(it_mt + 1);
  metricTensorDeterminant = *it_mt **(it_mt + 3) - *(it_mt + 1) **(it_mt + 2);
}
/*
--------------------------------------------------------------------------------
 compute surface gradient gradS from 2D reference gradient (grad0, grad1)
 (grad0, grad1)^t is reference gradient, gradS is output surface gradient
    Compute z0 = mt^{-1} t0 and z1 = mt^{-1} t1
     where t0, t1 are reference tangent vectors (columns of rectangular jacobian matrix)
     and mt is the metric tensor : {mt}_ij = ti.tj
    division by tensor matrix determinant is made further down in final linear combination
 only 3D
--------------------------------------------------------------------------------
*/
void GeomMapData::computeSurfaceGradient(real_t grad0, real_t grad1, std::vector<real_t>& gradS)
{
  std::vector<real_t> z0(spaceDim, 0.), z1(spaceDim, 0.);
  std::vector<real_t>::iterator it_gs = gradS.begin(), it_mt = metricTensor.begin(),
                                it_jac = jacobianMatrix.begin(), it_z0 = z0.begin(), it_z1 = z1.begin();

  for (dimen_t i = 0; i < spaceDim; i++, it_jac += elementDim)
    {
      *it_z0++ = -*(it_mt + 1) **it_jac + *it_mt **(it_jac + 1);
      *it_z1++ =  *(it_mt + 3) **it_jac - *(it_mt + 1) **(it_jac + 1);
    }

  //gradS = grad0 z1 + grad1 z0
  it_z0 = z0.begin();
  it_z1 = z1.begin();
  real_t a0 = grad0 / metricTensorDeterminant, a1 = grad1 / metricTensorDeterminant;
  for (dimen_t i = 0; i < spaceDim; i++, it_z1++, it_z0++, it_gs++)
    { *it_gs = a0 **it_z1 + a1 **it_z0; }
}

/*
--------------------------------------------------------------------------------
 print GeomMapadata
--------------------------------------------------------------------------------
*/
void GeomMapData::print(std::ostream& os) const
{
  os<<"MeshElement pointer="<<geomElement_p;
  if(theVerboseLevel>10) os<<" MeshElement : "<<(*geomElement_p)<<eol;
  os<<"current point="<<currentPoint<<eol;
  os<<"jacobianMatrix="<<jacobianMatrix<<eol;
  os<<"jacobian Determinant="<<jacobianDeterminant<<eol;
  os<<"differential element="<<differentialElement<<eol;
  os<<"inverse jacobian Matrix="<< inverseJacobianMatrix<<eol;
  os<<"normal vector="<<normalVector<<eol;
  os<<"metric tensor="<<metricTensor<<eol;
  os<<"metric tensor determinant="<<metricTensorDeterminant<<eol;
}

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

std::vector< Vector<real_t> >& GeomMapData::sideNV() const
{
  if (sideNormalVectors.size()==0)
    {
      number_t ns=geomElement_p->numberOfSides();
      sideNormalVectors.resize(ns);
      for (number_t s=0; s<ns; s++)
        {
          sideNormalVectors[s]=geomElement_p->normalVector(s+1);
          sideNormalVectors[s]/=norm2(sideNormalVectors[s]);
        }
    }
  return sideNormalVectors;
}

} // end of namespace xlifepp