term/computation/termUtils.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 termUtils.cpp
  \author N. Kielbasiewicz
  \since 27 feb 2013
  \date 1 mar 2013

  \brief Implementation of term utility functions
*/

#include "termUtils.hpp"
#include "../TermMatrix.hpp"
#include "../KernelOperatorOnTermVector.hpp"
#include "../../form/BilinearForm.hpp"
#include "SauterSchwabIM.hpp"
#include "DuffyIM.hpp"
#include <vector>
#include <sstream>

namespace xlifepp
{

//! special cast function used by FEComputation's
real_t toRealComplex(const Vector<real_t>& v,const complex_t& c) {return c.real();}
complex_t toRealComplex(const Vector<complex_t>&,const complex_t& c) {return c;}
real_t toRealComplex(const Vector<Vector<real_t> >& v,const complex_t& c) {return c.real();}
complex_t toRealComplex(const Vector<Vector<complex_t> >& v,const complex_t& c) {return c;}

// build a storage from a pair of Space
// for FE space it is assumed that only connected dofs in contiguous elements give non zero values
// FE spaces must live on same domain
MatrixStorage* buildStorage(const Space& rs, const Space& cs, StorageType st, AccessType at)
{
  trace_p->push("buildStorage");

  // the storage id is built from the adresses of Space objects
  std::stringstream ss;
  ss << &rs << "-" << &cs;
  MatrixStorage* res = findMatrixStorage(ss.str(), st, at);

  if(res != 0)
    {
      // the storage already exists, then we return a pointer to
      trace_p->pop();
      return res;
    }

  // we have to build a new storage
  //detect special case : one domain is the extension of the other
  const MeshDomain* cdom = cs.domain()->meshDomain(), *rdom = rs.domain()->meshDomain();
  string_t extension="";
  if(cdom!=0 && rdom!=0 && cdom!=rdom)
    {
      if(cdom->extensionof_p == rdom) extension="c=ext_r";
      if(rdom->extensionof_p == cdom) extension="r=ext_c";
    }
  std::vector<number_t>::const_iterator cit_r, cit_c;

  //build col index
  number_t dimr=rs.nbDofs(), dimc=cs.nbDofs();
  std::vector<std::set<number_t> > colset;      // temporary vector of col index
  if(st!=_dense)
    {
      colset.resize(dimr);
      if(rs.typeOfSubSpace() == _spSpace || cs.typeOfSubSpace() == _spSpace || (cs.domain()!=rs.domain() && extension==""))
        {
          std::set<number_t> colInd;
          std::set<number_t>::iterator jt=colInd.begin();
          for(number_t j = 0; j < dimc; j++) jt=colInd.insert(jt, j + 1);
          std::vector<std::set<number_t> >::iterator it=colset.begin();
          for(number_t i = 0; i < dimr; i++, it++) { *it = colInd; }
        }

      else if(extension=="c=ext_r") // case of extension of row domain
        {

          std::set<GeomElement*> elts_ext(cdom->geomElements.begin(),cdom->geomElements.end()); //set of column geom elements
          number_t n = rs.nbOfElements();
          for(number_t k = 0; k < n; k++)  //loop on side domain elements
            {
              std::vector<number_t> red = rs.elementDofs(k);            // row dof ranks for element k
              GeomElement* gelt = rs.element_p(k)->geomElt_p;           // geometric element supporting side element
              std::vector<GeoNumPair>& parents=gelt->parentSides();     // parent of side element
              std::vector<GeoNumPair>::iterator itp=parents.begin();
              for(; itp!=parents.end(); ++itp)  //loop on parents
                if(elts_ext.find(itp->first)!=elts_ext.end())     //parent elts in extended domain
                  {
                    std::vector<number_t> ced = cs.elementDofs(cs.numElement(itp->first)); // col dof ranks for element k ### A VERIFIER
                    for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
                      for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
                        colset[*cit_r - 1].insert(*cit_c);
                  }
            }
        }
      else if(extension=="r=ext_c") // case of a extension of col domain
        {
          std::set<GeomElement*> elts_ext(rdom->geomElements.begin(),rdom->geomElements.end());//set of row geom elements
          number_t n = cs.nbOfElements();
          for(number_t k = 0; k < n; k++)  //loop on side domain elements
            {
              std::vector<number_t> ced = cs.elementDofs(k);            // col dof ranks for element k
              GeomElement* gelt = cs.element_p(k)->geomElt_p;           // geometric element supporting side element
              std::vector<GeoNumPair>& parents=gelt->parentSides();     // parent of side element
              std::vector<GeoNumPair>::iterator itp=parents.begin();
              for(; itp!=parents.end(); ++itp)  //loop on parents
                if(elts_ext.find(itp->first)!=elts_ext.end())     //parent elts in extended domain
                  {
                    std::vector<number_t> red = rs.elementDofs(rs.numElement(itp->first)); // row dof ranks for element k ### A VERIFIER
                    for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
                      for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
                        colset[*cit_r - 1].insert(*cit_c);
                  }
            }
        }
      else     // both spaces are standard finite element spaces
        {
          number_t n = rs.nbOfElements();
          // build list of col index for each rows
          if(at!=_sym)
            {
              for(number_t k = 0; k < n; k++)
                {
                  // dof ranks for element k
                  std::vector<number_t> red = rs.elementDofs(k);               // row dof ranks for element k
                  GeomElement* gelt = rs.element_p(k)->geomElt_p;              // geometric element supporting element
                  number_t kc=k;
                  if(cs.element_p(k)->geomElt_p!=gelt) kc=cs.numElement(gelt); //not the same support, find the right one
                  std::vector<number_t> ced = cs.elementDofs(kc);              // col dof ranks for element k
                  for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
                    for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
                      colset[*cit_r - 1].insert(*cit_c);
                }
            }
          else //symmetric storage : keep only  strict lower part indices
            {
              for(number_t k = 0; k < n; k++)
                {
                  std::vector<number_t> red = rs.elementDofs(k);               // row dof ranks for element k
                  GeomElement* gelt = rs.element_p(k)->geomElt_p;              // geometric element supporting element
                  number_t kc=k;
                  if(cs.element_p(k)->geomElt_p!=gelt) kc=cs.numElement(gelt); //not the same support, find the right one
                  std::vector<number_t> ced = cs.elementDofs(kc);              // col dof ranks for element k
                  for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
                    for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
                      if(*cit_r > *cit_c) colset[*cit_r - 1].insert(*cit_c);
                }
            }
        }
    }
  // construct the storage with the column indices by rows
  switch(st)
    {
      case _cs :
        switch(at)
          {
            case _sym :  trace_p->pop(); return new SymCsStorage(dimr, colset, _lower, ss.str());
            case _row :  trace_p->pop(); return new RowCsStorage(dimr, dimc, colset, ss.str());
            case _col :  trace_p->pop(); return new ColCsStorage(dimr, dimc, colset, ss.str());
            case _dual :  trace_p->pop(); return new DualCsStorage(dimr, dimc, colset, ss.str());
            default : error("storage_bad_access", words("access type",at), words("storage type",_cs));
          }
        break;
      case _skyline :
        switch(at)
          {
            case _sym :  trace_p->pop(); return new SymSkylineStorage(dimr, colset, ss.str());
            case _dual :  trace_p->pop(); return new DualSkylineStorage(dimr, dimc, colset, ss.str());
            default : error("storage_bad_access", words("access type",at), words("storage type",_skyline));
          }
        break;
      case _dense :
        switch(at)
          {
            case _sym :  trace_p->pop(); return new SymDenseStorage(dimr, ss.str());
            case _row :  trace_p->pop(); return new RowDenseStorage(dimr, dimc, ss.str());
            case _col :  trace_p->pop(); return new ColDenseStorage(dimr, dimc, ss.str());
            case _dual :  trace_p->pop(); return new DualDenseStorage(dimr, dimc, ss.str());
            default : error("storage_bad_access", words("access type",at), words("storage type",_dense));
          }
        break;
      default : error("storage_not_implemented","buildStorage",words("storage type",st));
    }
  trace_p->pop();
  return 0;
}

/*! modify MatrixStorage stm to include dense submatrix given by its row and column indices (rows and cols)
    the new storage is of type (st,at)
    rows (resp. cols) has to be a subset of row indices (resp. column indices) of the stm storage
    if overwrite=true the current storage is overwritten and returned
    else a new MatrixStorage is returned, the stm storage being not cleared!
*/

MatrixStorage* updateStorage(MatrixStorage& stm, const std::vector<number_t> rows, const std::vector<number_t> cols,
                             StorageType st, AccessType at, bool overwrite)
{
  trace_p->push("updateStorage(MatrixStorage, rows, cols,...)");
  MatrixStorage* newstm=&stm;
  if(!overwrite) newstm= stm.clone();
  if(stm.storageType()==st && stm.accessType()==at)  //target storage is same as stm storage
    {
      switch(st)
        {
          case _cs :
          case _skyline : newstm->addSubMatrixIndices(rows,cols);
          default :  //dense in dense, nothing to do
            trace_p->pop();
            return newstm;
        }
    }

  error("storage_not_handled", words("storage type",stm.storageType()),words("access type",stm.accessType()));
  trace_p->pop();
  return newstm;
}

// build a storage from  type and dimension, no allocation of pointers (void matrix)
MatrixStorage* buildStorage(StorageType st, AccessType at, number_t dimr, number_t dimc)
{
  switch(st)
    {
      case _cs :
        switch(at)
          {
            case _sym :  return new SymCsStorage(dimr);
            case _row :  return new RowCsStorage(dimr, dimc);
            case _col :  return new ColCsStorage(dimr, dimc);
            case _dual :  return new DualCsStorage(dimr, dimc);
            default : error("storage_bad_access", words("access type",at), words("storage type",_cs));
          }
        break;
      case _skyline :
        switch(at)
          {
            case _sym :  return new SymSkylineStorage(dimr);
            case _dual : return new DualSkylineStorage(dimr, dimc);
            default : error("storage_bad_access", words("access type",at), words("storage type",_skyline));
          }
        break;
      case _dense :
        switch(at)
          {
            case _sym :  return new SymDenseStorage(dimr);
            case _row :  return new RowDenseStorage(dimr, dimc);
            case _col :  return new ColDenseStorage(dimr, dimc);
            case _dual : return new DualDenseStorage(dimr, dimc);
            default : error("storage_bad_access", words("access type",at), words("storage type",_dense));
          }
        break;
      default : error("storage_not_implemented","buildStorage",words("storage type",st));
    }
  return 0;
}

// build a storage from  type, dimension and column indices (vector of vectors)
MatrixStorage* buildStorage(StorageType st, AccessType at, number_t dimr, number_t dimc, const std::vector<std::vector<number_t> >& indices)
{
  trace_p->push("buildStorage");
  string_t str="";
  switch(st)
    {
      case _cs :
        switch(at)
          {
            case _sym :  trace_p->pop(); return new SymCsStorage(dimr, indices, _all, str);
            case _row :  trace_p->pop(); return new RowCsStorage(dimr, dimc, indices, str);
            case _col :
              {
                trace_p->pop();
                std::vector<std::vector<number_t> > rowindices(dimc);
                std::vector<std::vector<number_t> >::const_iterator itr=indices.begin();
                std::vector<number_t>::const_iterator itc;
                number_t r=1;
                for(; itr!=indices.end(); itr++, r++)
                  for(itc=itr->begin(); itc!= itr->end(); itc++) rowindices[*itc - 1].push_back(r);
                return new ColCsStorage(dimr, dimc, rowindices, str);
              }
            case _dual :  trace_p->pop(); return new DualCsStorage(dimr, dimc, indices, str);
            default : error("storage_bad_access", words("access type",at), words("storage type",_cs));
          }
        break;
      case _skyline :
        switch(at)
          {
            case _sym :  trace_p->pop(); return new SymSkylineStorage(dimr, indices, str);
            case _dual :  trace_p->pop(); return new DualSkylineStorage(dimr, dimc, indices, str);
            default : error("storage_bad_access", words("access type",at), words("storage type",_skyline));
          }
        break;
      case _dense :
        switch(at)
          {
            case _sym :  trace_p->pop(); return new SymDenseStorage(dimr, str);
            case _row :  trace_p->pop(); return new RowDenseStorage(dimr, dimc, str);
            case _col :  trace_p->pop(); return new ColDenseStorage(dimr, dimc, str);
            case _dual :  trace_p->pop(); return new DualDenseStorage(dimr, dimc, str);
            default : error("storage_bad_access", words("access type",at), words("storage type",_dense));
          }
        break;
      default : error("storage_not_implemented","buildStorage",words("storage type",st));
    }

  trace_p->pop();
  return 0;
}

// build a storage from  type, dimension and column indices (vector of sets)
MatrixStorage* buildStorage(StorageType st, AccessType at, number_t dimr, number_t dimc, const std::vector<std::set<number_t> >& colIndices)
{
  // convert set to vector
  std::vector<std::vector<number_t> > indices(colIndices.size());
  std::vector<std::set<number_t> >::const_iterator itc;
  std::vector<std::vector<number_t> >::iterator iti = indices.begin();
  for(itc = colIndices.begin(); itc != colIndices.end(); ++itc, ++iti)
    *iti = std::vector<number_t>(itc->begin(), itc->end());
  return buildStorage(st,at,dimr,dimc,indices);
}

/*! compute interpolated normals on Lagrange Dofs of space sp
    solve the projection problem Mk*nk = M0*n0 where
       Mk is the mass matrix associated to sp interpolation
       M0 is the "hybrid" mass matrix between sp interpolation and P0 interpolation
    the normals are returned as a vector of Vector

    NOTE : this method is well adapted for planar element (same normal every where in element)
           in case of curved geometric element, the method may be improved by using Gauss Lobatto interpolation of higher degree
*/
void interpolatedNormals(Space& sp, std::vector<Vector<real_t> >& ns)
{
  if(!sp.isFE())
    {
      where("interpolatedNormals(Space, vector<Vector<Real> >");
      error("not_fe_space_type", sp.name());
    }
  const Interpolation* intp=sp.rootSpace()->interpolation();
  if(intp->type!=_Lagrange) error("lagrange_fe_space_only",sp.name());

  trace_p->push("interpolatedNormals");

  if(intp->numtype > 0)  //non P0 case
    {

      const GeomDomain& dom=*sp.domain();
      dimen_t d=dom.mesh()->spaceDim();
      Unknown n(*sp.rootSpace(),"n",d);
      TestFunction t(n,"t");

      //define matrix
      TermMatrix Mk(intg(dom,n|t),"normalPkM");
      Space V0(dom,P0,"V0",false);
      Unknown n0(V0,"n0",d);
      TermMatrix M0(intg(dom,n0|t),"normalPkP0M");

      //fill normal P0 in a TermVector
      TermVector N0(n0, dom, Vector<real_t>(d,0.),"normalsP0");
      Vector<Vector<real_t> >::iterator it=N0.begin()->second->entries()->beginrv();
      for(number_t k=0; k<sp.nbOfElements(); k++, it++)
        {
          const Element* elt=sp.element_p(k);
          MeshElement* melt = elt->geomElt_p->meshElement();
          if(melt==0) melt=elt->geomElt_p->buildSideMeshElement();
          GeomMapData gmap(melt);
          gmap.computeJacobianMatrix(elt->refElt_p->geomRefElem_p->centroid());
          gmap.computeOutwardNormal();
          *it = gmap.normalVector;
        }

      //solve projection problem
      TermVector Nk=directSolve(Mk,M0*N0);

      //return normals
      ns= *Nk.begin()->second->entries()->rvEntries_p;
      std::vector<Vector<real_t> >::iterator itn=ns.begin();
      for(; itn!=ns.end(); itn++)  *itn/=itn->norm2();
    }

  else //P0 interpolation
    {
      verboseLevel(10);
      ns.resize(sp.nbOfElements());
      std::vector<Vector<real_t> >::iterator it=ns.begin();
      for(number_t k=0; k<sp.nbOfElements(); k++, it++)
        {
          const Element* elt=sp.element_p(k);
          MeshElement* melt = elt->geomElt_p->meshElement();
          if(melt==0) melt=elt->geomElt_p->buildSideMeshElement();
          GeomMapData gmap(melt);
          gmap.computeJacobianMatrix(elt->refElt_p->geomRefElem_p->centroid());
          gmap.computeOutwardNormal();
          *it = gmap.normalVector;
        }
    }
  trace_p->pop();
  return;
}

//compute normals of a side domain using interpolation given by an unknown
TermVector normalsOn(GeomDomain& dom, const Unknown& u)
{
  if (!dom.isSideDomain() && ((dom.dim()+1)!=dom.spaceDim()))
  {
    where("normalsOn(Domain, Unknown)");
    error("domain_not_side", dom.name());
  }
  const Space* spu = u.space();
  const Interpolation* itp0 = findInterpolation(Lagrange, _standard, 0, L2);
  const Space* sp0=0; bool sub = false;
  if (spu->interpolation()==itp0)
  {
    if (spu->domain() == &dom) sp0 = spu;  // u space is a P0 space
    else if (spu->domain()->include(dom)) {sp0 = spu; sub=true;}
  }
  if (sp0==0) // find or create space P0
  {
    std::vector<Space*>::iterator it;
    for (it = Space::theSpaces.begin(); it != Space::theSpaces.end() && (sp0==0 || sub==false); ++it)
    {
      const Space* spit = (*it)->space();
      if (spit != 0)
      {
        const Interpolation* itp = spit->interpolation();
        if (itp == itp0)
        {
          if (spit->domain()==&dom) sp0=spit;
          else if (spit->domain()->include(dom)) {sp0=spit; sub=true;}
        }
      }
    }
  }
  if (sp0==0) sp0 = new Space(dom, P0, "P0");           // no P0 space found
  else if (sub) sp0 = Space::findSubSpace(&dom, sp0);    // P0 parent space found
  dimen_t d = dom.spaceDim();
  Unknown n0(*const_cast<Space*>(sp0), "n_"+dom.name(),d);

  // fill normal P0 in a TermVector
  TermVector N0(n0, dom, Vector<real_t>(d,0.),"normalsP0");
  Vector<Vector<real_t> >::iterator it=N0.begin()->second->entries()->beginrv();
  for (number_t k=0; k<sp0->nbOfElements(); k++, it++)
  {
    if (!dom.meshDomain()->normalComputed)
    {
      const Element* elt=sp0->element_p(k);
      MeshElement* melt = elt->geomElt_p->meshElement();
      if (melt==0) melt=elt->geomElt_p->buildSideMeshElement();
      GeomMapData gmap(melt);
      gmap.computeJacobianMatrix(elt->refElt_p->geomRefElem_p->centroid());
      gmap.computeOutwardNormal();
      *it = gmap.normalVector;
    }
    else *it = sp0->element_p(k)->geomElt_p->meshElement()->geomMapData_p->normalVector;

  }

  if (spu->interpolation()==itp0) return N0;

  //project normal on space spu
  Unknown* uc = const_cast<Unknown*>(&u);
  //if(spu->dimFun()==1 && u.nbOfComponents()==1)  uc=new Unknown(*const_cast<Space*>(spu),"uc_"+tostring(spu),d); // permissive, user pass a scalar unknown
  if (uc->nbOfComponents()!=d && spu->dimFun()!=d) error("unknown_inconsistent_size", uc->name(), d, uc->nbOfComponents());
  TestFunction t(*uc,"t_u");
  TermMatrix Mk(intg(dom,(*uc)|t),"u_t_mass");
  TermMatrix M0(intg(dom,n0|t),"n0_t_mass");
  TermVector Nk=directSolve(Mk,M0*N0);
  return Nk;
}

/*! precompute geometry data and shape values on any element for any quadrature involved
    restricted to one order geometric element, take large place but computation is then faster in case of IE!!!

    dom   : integration domain
    subsp : subspace providing the list of elements
    quads : list of quadratures
    nbc   : number of components of the unknown
    der   : compute shape derivatives if true
    nor   : compute normal if true
    extendedDomain : true when dealing with an extended domain
*/
void preComputationIE(const GeomDomain& dom, const Space* subsp, std::set<Quadrature*>& quads, dimen_t nbc, bool der, bool nor,
                      bool extendedDomain)
{
  der=true;  //force computation of derivatives to avoid misfit when recall

  //shapevalues for all single quadrature (non parallel)
  for(number_t k = 0; k < subsp->nbOfElements(); k++) // data related to element
    {
      const Element* elt = subsp->element_p(k);
      RefElement* relt = elt->refElt_p;                              //reference element related to element
      GeomElement* gelt = elt->geomElt_p;
      const GeomRefElement* grelt = gelt->geomRefElement();
      number_t sw = relt->nbDofs()*nbc;
      if(extendedDomain) //build geometric data for extended element
        {
          if(gelt->meshElement()==0) gelt->buildSideMeshElement();
          MeshElement* melt = gelt->meshElement();
          GeomMapData* mapdata = melt->geomMapData_p;
          if(mapdata==0) mapdata=new GeomMapData(melt);
          mapdata->computeJacobianMatrix(std::vector<real_t>(gelt->elementDim(),0.));
          mapdata->invertJacobianMatrix();
          melt->geomMapData_p = mapdata;
        }
      //precompute reference shape values
      std::set<Quadrature*>::iterator itq=quads.begin();
      for(; itq!=quads.end(); ++itq)
        {
          std::map<Quadrature*,std::vector<ShapeValues> >::iterator itmq=relt->qshvs_.find(*itq);
          if(itmq==relt->qshvs_.end() || itmq->second.begin()->w.size()!=sw) //shapevalues re-computed
            {
              dimen_t dq = (*itq)->dim();                               //dimension of quadrature points
              number_t nq = (*itq)->numberOfPoints();                   //number of quadrature points
              if(!extendedDomain)  //no extension
                {
                  std::vector<real_t>::const_iterator itp = (*itq)->point();
                  relt->qshvs_[*itq]=std::vector<ShapeValues>(nq);          //to store all shapevalues at quadrature points
                  std::vector<ShapeValues>::iterator itshv = relt->qshvs_[*itq].begin();
                  for(number_t i = 0; i < nq; i++, itp += dq, itshv++)
                    {
                      relt->computeShapeValues(itp, der);
                      *itshv = relt->shapeValues;
                      if(nbc>1) itshv->extendToVector(nbc);    //extend scalar shape functions to nbc vector shape functions
                    }
                }
              else // in case of extension build shapevalues for each side of ref element stored sequentially
                {
                  number_t nbs = relt->geomRefElem_p->nbSides();
                  relt->qshvs_[*itq]=std::vector<ShapeValues>(nbs*nq);       //to store all shapevalues at quadrature points
                  std::vector<ShapeValues>::iterator itshv = relt->qshvs_[*itq].begin();
                  for(number_t s=1; s<=nbs; ++s) //loop on side
                    {
                      std::vector<real_t>::const_iterator itp = (*itq)->point();
                      for(number_t i = 0; i < nq; i++, itp += dq, itshv++)
                        {
                          relt->computeShapeValues(grelt->sideToElt(s,itp).begin(), der);
                          *itshv = relt->shapeValues;
                          if(nbc>1) itshv->extendToVector(nbc);    //extend scalar shape functions to nbc vector shape functions
                        }
                    }
                }
            }
        }
    }
    //precompute signDofs of element if required
    #ifdef XLIFEPP_WITH_OMP
    #pragma omp parallel for
    #endif // XLIFEPP_WITH_OMP
    for(number_t k = 0; k < subsp->nbOfElements(); k++) // data related to element
    {
      const Element* elt = subsp->element_p(k);
      if(elt->refElt_p->dofCompatibility == _signDofCompatibility)  elt->setDofSigns();
    }

  MeshDomain* mdom=const_cast<MeshDomain*>(dom.meshDomain());
  //compute map of quadrature points (parallel)
#ifdef XLIFEPP_WITH_OMP
  #pragma omp parallel for
#endif // XLIFEPP_WITH_OMP
  for(number_t k = 1; k <= mdom->numberOfElements(); k++)
    {
      GeomElement* gelt =mdom->element(k);
      //geometry information
      if(gelt->meshElement()==0) gelt->buildSideMeshElement();
      MeshElement* melt = gelt->meshElement();
      GeomMapData* mapdata=melt->geomMapData_p;
      if(mapdata==0) mapdata=new GeomMapData(melt);
      mapdata->computeJacobianMatrix(std::vector<real_t>(gelt->elementDim(),0.));
      mapdata->invertJacobianMatrix();
      mapdata->computeDifferentialElement();
      melt->geomMapData_p = mapdata;
      if(nor) mapdata->computeOutwardNormal();
      //physical points
      std::set<Quadrature*>::iterator itq=quads.begin();
      for(; itq!=quads.end(); ++itq)
        {
          dimen_t dq = (*itq)->dim();                                 //dimension of quadrature points
          number_t nq = (*itq)->numberOfPoints();                     //number of quadrature points
          std::map<Quadrature*,std::vector<Point> >::iterator itmqphy= mapdata->phyPoints.find(*itq);
          if(itmqphy==mapdata->phyPoints.end())  //create physical points
            {
              mapdata->phyPoints[*itq] = std::vector<Point>(nq);
              std::vector<Point>::iterator itphy=mapdata->phyPoints[*itq].begin();
              std::vector<real_t>::const_iterator itp = (*itq)->point();
              for(number_t i = 0; i < nq; i++, itp += dq, ++itphy) *itphy = mapdata->geomMap(itp);
            }
        }
    }
}

/*! precompute geometry data and shape values on any element for any quadrature involved
    restricted to one order geometric element, take large place but computation is then faster in case of IE!!!

    subsp : subspace providing the list of elements
    quads : list of quadratures
    nbc   : number of components of the unknown
    der   : compute shape derivatives if true
    nor   : compute normal if true
*/
void preComputationIE_old(const Space* subsp, std::set<Quadrature*>& quads, dimen_t nbc, bool der, bool nor)
{
  der=true;  //force computation of derivatives to avoid misfit when recall

  for(number_t k = 0; k < subsp->nbOfElements(); k++) // data related to u
    {
      const Element* elt = subsp->element_p(k);
      GeomElement* gelt = elt->geomElt_p;

      //geometry information
      if(gelt->meshElement()==0) gelt->buildSideMeshElement();
      MeshElement* melt = gelt->meshElement();
      GeomMapData* mapdata=melt->geomMapData_p;
      if(mapdata==0) mapdata=new GeomMapData(melt);
      mapdata->computeJacobianMatrix(std::vector<real_t>(elt->dim(),0.));
      mapdata->invertJacobianMatrix();
      mapdata->computeDifferentialElement();
      melt->geomMapData_p = mapdata;
      if(nor) mapdata->computeOutwardNormal();

      //shapevalues and physical points for all single quadrature
      RefElement* relt = elt->refElt_p;                              //reference element related to element
      std::set<Quadrature*>::iterator itq=quads.begin();
      for(; itq!=quads.end(); ++itq)
        {
          std::map<Quadrature*,std::vector<ShapeValues> >::iterator itmq=relt->qshvs_.find(*itq);
          dimen_t dq = (*itq)->dim();                                 //dimension of quadrature points
          number_t nq = (*itq)->numberOfPoints();                     //number of quadrature points
          if(itmq==relt->qshvs_.end()) //shapevalues not computed
            {
              relt->qshvs_[*itq]=std::vector<ShapeValues>(nq);          //to store all shapevalues at quadrature points
              std::vector<real_t>::const_iterator itp = (*itq)->point();
              std::vector<ShapeValues>::iterator itshv = relt->qshvs_[*itq].begin();
              for(number_t i = 0; i < nq; i++, itp += dq, itshv++)
                {
                  relt->computeShapeValues(itp, der);
                  *itshv = relt->shapeValues;
                  if(nbc>1) itshv->extendToVector(nbc);    //extend scalar shape functions to nbc vector shape functions
                }
            }
          //compute map of quadrature points
          std::map<Quadrature*,std::vector<Point> >::iterator itmqphy= mapdata->phyPoints.find(*itq);
          if(itmqphy==mapdata->phyPoints.end())  //create physical points
            {
              mapdata->phyPoints[*itq] = std::vector<Point>(nq);
              std::vector<Point>::iterator itphy=mapdata->phyPoints[*itq].begin();
              std::vector<real_t>::const_iterator itp = (*itq)->point();
              for(number_t i = 0; i < nq; i++, itp += dq, ++itphy)
                *itphy = mapdata->geomMap(itp);                      //quadrature point in physical space
            }
        }
    }
}

//===================================================================================================================================
/*! check and set quadrature pointers of IntegrationsMethods if they are not set

  if intgMethods is empty
      if kernel is singular (Kernel::singularType!=_notSingular) add if possible
          - a well suited singular quadrature method applied to adjacent elements (dist=0)
          - some regular quadrature products applied to non adjacent elements (one for 0<dist<=1, one for 1<dist<=2 and one for dist>2)
      if kernel is regular (Kernel::singularType==_notSingular) add
          - a regular quadrature products
  if intgMethods contains only a singular quadrature
      if kernel is regular stop on error
      if kernel is singular, check if singular quadrature is consistent
         if not stop on error
         else add some regular quadrature products applied to non adjacent elements (one for 0<dist<=1, one for 1<dist<=2 and one for dist>2)
  if intgMethods contains only a regular quadrature
      if kernel is regular check if it is consistent
      if kernel is singular add a well suited singular quadrature method applied to adjacent elements (dist=0)
  if intgMethods are more than one integration method
     if kernel is singular,
        check if there is only one consistent singular quadrature for dist=0, if not stop on error
        if there is no regular quadrature, add some
     if kernel is regular, check that all are regular quadratures

 At the end of the process, produce an intgMethods vector consistent with the Kernel
 In a second step, update pointers of quadrature if required
*/
// Implemented here due to a cyclic dependance library
//===================================================================================================================================
void DoubleIntgBilinearForm::setIntegrationMethods()
{
  const Kernel* kerp = kopus_p->kernel();
  string_t mes="in DoubleIntgBilinearForm::setIntegrationMethods(), ";
  const std::set<ShapeType>& shapes_u = domainu_p->meshDomain()->shapeTypes;
  const std::set<ShapeType>& shapes_v = domainv_p->meshDomain()->shapeTypes;
  if (shapes_u.size()!=1 || shapes_v.size()!=1)
  {
    where("DoubleIntgBilinearForm::setIntegrationMethods()");
    error("domain_multiple_shapes");
  }
  ShapeType shu = *shapes_u.begin(), shv = *shapes_v.begin();

  //analyze given integration methods
  SingularityType st=_notsingular;
  real_t so = 0;
  bool hasSingPart = false, hasRegPart = false;
  string_t sn ="Id", ssn="Id";
  number_t nbsq=0, nbrq=0, i=0;
  bool waitRegularPart = false;
  if (kerp!=0) //if kerp not allocated, assume ker=Id
  {
    st=kerp->singularType;
    so=kerp->singularOrder;
    hasSingPart = kerp->singPart!=0;
    hasRegPart  = kerp->regPart!=0;
    sn=kerp->shortname;
    if (hasSingPart) ssn = kerp->singPart->shortname;
  }
  std::multimap<real_t, int_t> ims;  //intg methods number indexed by distance (negative number for singular quadrature method)
  std::multimap<real_t, int_t>::iterator itm, itme;
  std::vector<IntgMeth>::iterator iti=intgMethods.begin();
  for (; iti!=intgMethods.end(); ++iti, ++i)  //loop on integration methods
  {
    const IntegrationMethod* im =iti->intgMeth;
    if (im->imType==_quadratureIM) //single quadrature, build the product of quadrature
    {
      const QuadratureIM* qim = static_cast<const QuadratureIM*>(im);
      iti->intgMeth = new ProductIM(qim->quadRule,qim->degree);
      im = iti->intgMeth;
    }
    if (!im->isDoubleIM())
    {
      where("DoubleIntgBilinearForm::setIntegrationMethods()");
      error("im_not_single");
    }
    FunctionPart pai=iti->functionPart;
    if (pai==_regularPart && !hasRegPart)
    {
      where("DoubleIntgBilinearForm::setIntegrationMethods()");
      error("null_pointer","kerp->regPart");
    }
    if (pai==_singularPart && !hasSingPart)
    {
      where("DoubleIntgBilinearForm::setIntegrationMethods()");
      error("null_pointer","kerp->singPart");
    }
    SingularityType sti=im->singularType;
    real_t soi=im->singularOrder;
    real_t di = iti->bound;
    if (sti ==_notsingular)   //regular quadrature
    {
      if (st!=_notsingular && pai!=_regularPart && di==0)
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_not_handled",im->name);
      }
      ims.insert(std::pair<real_t,number_t>(di,i+1));
      nbrq++;  //regular quadrature is valid
    }
    else //singular quadrature
    {
      if (sti!=st || soi!=so)  //check singularity consitency
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_not_handled",im->name);
      }
      string_t nai=im->kerName;
      string_t na= sn;
      if (pai==_singularPart) {na=ssn; waitRegularPart = true;}
      if (nai!="" && nai!=na)  // check if intg method is dedicated to particular kernel
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_not_handled",im->name);
      }
      if (di>0)
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_not_handled",im->name);
      }
      IntegrationMethodType imt=im->type();
      string_t mes2="";
      if (imt==_SauterSchwabIM && (shu!=_triangle || shv!=_triangle))
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_shape_only", "Sauter-Schwab", words("shape", _triangle), words("shape",shu), words("shape",shv));
      }
      if (imt==_DuffyIM && (shu!=_segment || shv!=_segment))
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_shape_only", "Duffy", words("shape", _segment), words("shape",shu), words("shape",shv));
      }
      if (imt==_LenoirSalles3dIM && (shu!=_triangle || shv!=_triangle))
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_shape_only", "Lenoir-Salles3D", words("shape", _triangle), words("shape",shu), words("shape",shv));
      }
      if (imt==_LenoirSalles2dIM && (shu!=_segment || shv!=_segment))
      {
        where("DoubleIntgBilinearForm::setIntegrationMethods()");
        error("im_shape_only", "Lenoir-Salles2D", words("shape", _segment), words("shape",shu), words("shape",shv));
      }
      ims.insert(std::pair<real_t,number_t>(di,-(i+1)));
      nbsq++;  //singular quadrature is valid
    }
  }

  // finalize analyze and add omitted integration methods by using default ones
  if (st!=_notsingular && so<0 && domainu_p==domainv_p)  //kernel is singular, nbsq should be 1
  {
    if (nbsq>1) //too much singular quadrature
    {
      where("DoubleIntgBilinearForm::setIntegrationMethods()");
      string_t mes="";
      itm=ims.lower_bound(0.); itme = ims.upper_bound(0.);
      for (; itm!=itme; ++itm)
        if (itme->second<0) mes+=" "+intgMethods[std::abs(itme->second)-1].intgMeth->name;
      error("multiple_im",mes);
    }
    if (nbsq==0) //add default singular method  (NOT WORK FOR MOMENT BECAUSE DuffyIMNew and SauterSchwabNew are not included, lib dependance !)
    {
      if (shu==_segment &&  shv==_segment)  //use Duffy order 6
        intgMethods.push_back(IntgMeth(DuffyIM(),_allFunction, 0.));
      else
      {
        if (shu==_triangle && shv==_triangle)  //use Sauter-Schwab order 3
          intgMethods.push_back(IntgMeth(SauterSchwabIM(),_allFunction, 0.));
        else
        {
          where("DoubleIntgBilinearForm::setIntegrationMethods()");
          error("im_not_found");
        }
      }
      nbsq=1;
    }
  }

  //deal with regular quadratures
  if (st!=_notsingular && so<0)  //kernel is singular
  {
    if (waitRegularPart)  //find if there is a regular quadrature for regularPart and dist=0
    {
      itm=ims.lower_bound(0.); itme = ims.upper_bound(0.);
      for(; itm!=itme && waitRegularPart ; ++itm)
        if(itme->second>0 && intgMethods[itme->second-1].functionPart==_regularPart)  waitRegularPart=false;
    }
    if (waitRegularPart) // create a quadrature for regular part
      intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,3),_regularPart, 0.));
    if (nbrq==0) //add additional regular quadrature
    {
      intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,5),_allFunction, 1.));
      intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,4),_allFunction, 2.));
      intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,3),_allFunction, theRealMax));
    }
  }
  else  //kernel is regular
  {
    if(nbrq==0) intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,3),_allFunction, theRealMax));
  }

  // update QuadratureIM pointers of regular quadratures
  for (iti=intgMethods.begin(); iti!=intgMethods.end(); ++iti)
  {
    const IntegrationMethod* im =iti->intgMeth;
    if (im->type()==_productIM) //update single quadrature pointers
    {
      const ProductIM* pim=static_cast<const ProductIM*>(im);
      if (pim->getxIM()->type()==_quadratureIM) static_cast<QuadratureIM*>(pim->getxIM())->setQuadratures(shapes_v);
      if (pim->getyIM()->type()==_quadratureIM) static_cast<QuadratureIM*>(pim->getyIM())->setQuadratures(shapes_u);
      pim->name = pim->getxIM()->name+" x "+pim->getyIM()->name;
    }
  }
}

void DoubleIntgBilinearForm::setHMIntegrationMethods()
{
  if (intgMethod_p==0) return;  //no integration method
  if (intgMethod_p->type()!=_HMatrixIM) return;  //not a HMatrixIM method
  const HMatrixIM* him= static_cast<const HMatrixIM*>(intgMethod_p);
  intgMethods=him->intgMethods;
  setIntegrationMethods();
}

// IntegrationMethods constructors
// implemented here due to cyclic dependance library
IntegrationMethods::IntegrationMethods(IntegrationMethodType imt)
{
  switch (imt)
  {
    case _SauterSchwabIM :  add(SauterSchwabIM(),_allFunction, theRealMax); break;
    case _DuffyIM : add(DuffyIM(),_allFunction,0.); break;
    case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,theRealMax); break;
    case _LenoirSalles3dIM : add(LenoirSalles3dIM(),_allFunction,theRealMax); break;
    case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,theRealMax); break;
    case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,theRealMax); break;
    default:
      where("IntegrationMethods::IntegrationMethods(IntegrationMethodType)");
      error("im_not_handled", words("imtype", imt));
  }
}

IntegrationMethods::IntegrationMethods(QuadRule qr, number_t ord)
{
  add(QuadratureIM(qr,ord),_allFunction,theRealMax);
}

IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, QuadRule qr2, number_t ord2)
{
  switch (imt)
  {
    case _SauterSchwabIM :  add(SauterSchwabIM(),_allFunction, 0.); break;
    case _DuffyIM : add(DuffyIM(),_allFunction,0.); break;
    case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,0.); break;
    case _LenoirSalles3dIM : add(LenoirSalles3dIM(),_allFunction,0.); break;
    case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,0.); break;
    case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,0.); break;
    default:
      where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, QuadRule, Number)");
      error("im_not_handled", words("imtype", imt));
  }
  add(QuadratureIM(qr2,ord2),_allFunction,theRealMax);
}

IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, number_t ord1, real_t b1, QuadRule qr2, number_t ord2)
{
  switch (imt)
  {
    case _SauterSchwabIM :  add(SauterSchwabIM(ord1),_allFunction,b1); break;
    case _DuffyIM : add(DuffyIM(ord1),_allFunction,b1); break;
    case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
    case _LenoirSalles3dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
    case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,b1); break;
    case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,b1); break;
    default:
      where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, Number, Real, QuadRule, Number)");
      error("im_not_handled", words("imtype", imt));
  }
  add(QuadratureIM(qr2,ord2),_allFunction, theRealMax);
}

IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, number_t ord1, real_t b1,
                                       QuadRule qr2, number_t ord2, real_t b2, QuadRule qr3, number_t ord3)
{
  switch (imt)
  {
    case _SauterSchwabIM :  add(SauterSchwabIM(ord1),_allFunction,b1); break;
    case _DuffyIM : add(DuffyIM(ord1),_allFunction,b1); break;
    case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
    case _LenoirSalles3dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
    case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,b1); break;
    case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,b1); break;
    default:
      where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, Number, Real, QuadRule, Number, Real, QuadRule, Number)");
      error("im_not_handled", words("imtype", imt));
  }
  add(QuadratureIM(qr2,ord2),_allFunction,b2);
  add(QuadratureIM(qr3,ord3),_allFunction,theRealMax);
}

IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, number_t ord1, real_t b1,
                                       QuadRule qr2, number_t ord2, real_t b2, QuadRule qr3, number_t ord3, real_t b3, QuadRule qr4, number_t ord4)
{
  switch (imt)
  {
    case _SauterSchwabIM :  add(SauterSchwabIM(ord1),_allFunction,b1); break;
    case _DuffyIM : add(DuffyIM(ord1),_allFunction,b1); break;
    case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
    case _LenoirSalles3dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
    case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,b1); break;
    case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,b1); break;
    default:
      where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, Number, Real, QuadRule, Number, Real, QuadRule, Number, Real, QuadRule, Number)");
      error("im_not_handled", words("imtype", imt));
  }
  add(QuadratureIM(qr2,ord2),_allFunction,b2);
  add(QuadratureIM(qr3,ord3),_allFunction,b3);
  add(QuadratureIM(qr4,ord4),_allFunction,theRealMax);
}

IntegrationMethods::IntegrationMethods(QuadRule qr1, number_t ord1, real_t b1, QuadRule qr2, number_t ord2)
{
  add(QuadratureIM(qr1,ord1),_allFunction,b1);
  add(QuadratureIM(qr2,ord2),_allFunction,theRealMax);
}

IntegrationMethods::IntegrationMethods(QuadRule qr1, number_t ord1, real_t b1, QuadRule qr2, number_t ord2, real_t b2, QuadRule qr3, number_t ord3)
{
  add(QuadratureIM(qr1,ord1),_allFunction,b1);
  add(QuadratureIM(qr2,ord2),_allFunction,b2);
  add(QuadratureIM(qr3,ord3),_allFunction,theRealMax);
}

} // end of namespace xlifepp