term/computation/FEMatrixComputation.hpp

/*
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 FEMatrixComputation.hpp
  \author E. Lunéville
  \since 9 jan 2013
  \date 25 jun 2013

  \brief Implementation of template FE TermMatrix computation functions

  this file is included by TermMatrix.hpp
  do not include it elsewhere
*/

/* ================================================================================================
                                      tools
   ================================================================================================ */

namespace xlifepp
{

//data structure to store FE computation parameters
class FEcomputationParameters
{
  public:
    const IntgBilinearForm*  ibf;     // basic bilinear form as simple integral
    const OperatorOnUnknown* op_u;    // operator on unknown
    const OperatorOnUnknown* op_v;    // operator on test function
    bool id_u;                        // opu is Id
    bool id_v;                        // opv is Id
    dimen_t ord_opu;                  // derivative order of op_u
    dimen_t ord_opv;                  // derivative order of op_v
    bool has_fun;                     // true if one of operators has a data function
    bool same_val;                    // true if opu and opv have same values
    bool mapsh_u;                     // true if mapping of u-shapevalues is required
    bool mapsh_v;                     // true if mapping of v-shapevalues is required
    AlgebraicOperator aop;            // algebraic operator between differential operator
    const QuadratureIM* qim;          // quadrature pointer related to basic bilinear form

    FEcomputationParameters()
      : ibf(0), op_u(0), op_v(0), id_u(false), id_v(false),ord_opu(0), ord_opv(0), has_fun(false), same_val(true),
        mapsh_u(false), mapsh_v(false), aop(_product), qim(0) {}

    FEcomputationParameters(const IntgBilinearForm*  ib, const OperatorOnUnknown* opu, const OperatorOnUnknown* opv,
                            bool idu, bool idv, dimen_t ordopu, dimen_t ordopv, bool hasfun, bool sameval,
                            bool mapshu, bool mapshv, AlgebraicOperator ao, const QuadratureIM* q)
      : ibf(ib), op_u(opu), op_v(opv), id_u(idu), id_v(idv),ord_opu(ordopu), ord_opv(ordopv), has_fun(hasfun), same_val(sameval),
        mapsh_u(mapshu), mapsh_v(mapshv), aop(ao), qim(q) {}
};

//map shapevalues in physical space regarding some rules
// refElt  : FE ref elment
// melt    : geometric element
// mapdata : structure handling geometric stuff (jacobian, ...)
// mapsh   : true if mapping
// femt    : type of mapping
// rotsh   : if true shapevalues are rotated along a rotation given by the ref element
// ord     : if 0 apply mapping only on shape values, if 1 apply mapping on shape values and derivatives
// changesign : if true, sign of some shape values have to be reversed according to normal or tangent orientation
// sign    : vector of signs of dofs
// dimfun  : dimension of shape functions in reference element
// dimfunp : dimension of shape functions in physical element, set by this function
// sh      : shapeValues in reference element
// shmap   : mapped shapevalues in physical element
// NOTE : dimfunp may be modified by this function
inline void mapShapeValues(RefElement& refElt, MeshElement& melt, GeomMapData& mapdata, bool mapsh, FEMapType femt,
                           bool rotsh, number_t ord, bool changeSign, const Vector<real_t>* sign, dimen_t dimfun,
                           dimen_t& dimfunp, const ShapeValues& sh, ShapeValues& shmap)
{
  dimfunp=dimfun;
  if(&sh!=&shmap) shmap.assign(sh);
  if(!mapsh && !changeSign && !rotsh) return;  //nothing to do
  if(rotsh) refElt.rotateDofs(melt.verticesNumbers(), shmap, ord>0);    //rotate shape values
  if(mapsh)   //map the shapevalues according to FE map type
    {
      number_t nbfun = sh.w.size()/dimfun;  //nb of shape functions
      switch(femt)         //shape values in physical space
        {
          case _contravariantPiolaMap :
            shmap.contravariantPiolaMap(shmap, mapdata, ord);
            break;
          case _covariantPiolaMap :
            shmap.covariantPiolaMap(shmap, mapdata, ord);
            break;
          default :
            shmap.map(shmap, mapdata, ord);
        }
      dimfunp = shmap.w.size()/nbfun;  //new dimfun when 2d-3d mapping or 1d-2d mapping
    }
  if(changeSign) shmap.changeSign(*sign, dimfunp, ord);  //change sign of shape functions according to sign vector
}

// precompute geometry data and shape values on any element for any quadrature involved
// take place but computation is then faster in case of one order geometric element (not recomputed)
//
//   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
//   mapquad : map quadrature points to physical space if true, take place in memory !
//   useAux  : use auxiliary shapevalues map (relt->qshvs_aux) instead of standard (relt->qshvs_)
//             required when same quadrature on same reference element but number of component differs
inline void preComputationFE(Space* subsp, const MeshDomain* dom, std::set<Quadrature*>& quads, dimen_t nbc,
                             bool der, bool invertJacobian, bool nor, bool mapquad, bool useAux)
{

  //shapevalues for all  quadratures
  //=====  warning : computation of shapevalue derivatives is enforced ====//
  der = true;

  const std::set<RefElement*>& refelts=subsp->refElements();
  std::set<RefElement*>::const_iterator itref=refelts.begin();

  for(; itref!=refelts.end(); ++itref)
    {
      RefElement* relt = *itref;                              //reference element
      if(relt->hasShapeValues) // Ref Element has to propose shape functions!
        {
          std::set<Quadrature*>::iterator itq=quads.begin();
          for(; itq!=quads.end(); ++itq)
            {
              bool compute = false;
              std::map<Quadrature*,std::vector<ShapeValues> >::iterator itmq;
              if(!useAux)
                {
                  itmq =relt->qshvs_.find(*itq);
                  if(itmq==relt->qshvs_.end())  compute = true; //shapevalues not computed
                }
              else
                {
                  itmq =relt->qshvs_aux.find(*itq);
                  if(itmq==relt->qshvs_aux.end())  compute = true; //shapevalues not computed
                }
              if(!compute) //may be some attributes may have changed
                {
                  std::vector<ShapeValues>& shs = itmq->second;
                  if(shs.size()==0) compute = true;                            // hazardous empty shapevalues
                  else if(shs[0].w.size()!= relt->nbDofs()*nbc) compute=true;  // not the same size
                }
              if(compute) // shapevalues have to be (re)computed
                {
                  dimen_t dq = (*itq)->dim();                                 //dimension of quadrature points
                  number_t nq = (*itq)->numberOfPoints();                     //number of quadrature points
                  std::vector<ShapeValues>::iterator itshv;
                  if(!useAux)
                    {
                      relt->qshvs_[*itq]=std::vector<ShapeValues>(nq);            //allocate shapevalues at quadrature points
                      itshv = relt->qshvs_[*itq].begin();
                    }
                  else
                    {
                      relt->qshvs_aux[*itq]=std::vector<ShapeValues>(nq);            //allocate shapevalues at quadrature points
                      itshv = relt->qshvs_aux[*itq].begin();
                    }
                  std::vector<real_t>::const_iterator itp = (*itq)->point();
                  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
                    }
                }
            }
        }
    }

//  if(!mapquad &&  dom->jacobianComputed && dom->diffEltComputed
//      && (!nor || (nor && dom->jacobianComputed))
//      && (!invertJacobian || (invertJacobian && dom->inverseJacobianComputed))) return;

  //update  geometric data 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 u
    {
      const Element* elt = subsp->element_p(k);
      GeomElement* gelt = elt->geomElt_p;

      //geometry stuff
      if(gelt->meshElement()==0) gelt->buildSideMeshElement();
      MeshElement* melt = gelt->meshElement();
      GeomMapData* mapdata=melt->geomMapData_p;
      if(mapdata==0)
        {
          mapdata=new GeomMapData(melt);
          melt->geomMapData_p = mapdata;
          mapdata->computeJacobianMatrix(std::vector<real_t>(elt->dimElt(),0.));
          mapdata->computeJacobianDeterminant();// induced  mapdata->computeDifferentialElement();
          if(invertJacobian) mapdata->invertJacobianMatrix();
          if(nor) mapdata->computeOutwardNormal();
        }
      else
        {
          if(mapdata->jacobianMatrix.size()==0) mapdata->computeJacobianMatrix(std::vector<real_t>(elt->dimElt(),0.));
          if(mapdata->jacobianDeterminant==0) mapdata->computeJacobianDeterminant();
          //if(mapdata->differentialElement==0) mapdata->computeDifferentialElement();
          if(invertJacobian && mapdata->inverseJacobianMatrix.size()==0) mapdata->invertJacobianMatrix();
          if(nor && mapdata->normalVector.size()==0) mapdata->computeOutwardNormal();
        }

      if(mapquad)  //map quadrature points onto physical element
        {
          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);                      //quadrature point in physical space
                }
            }
        }
    }
}

//template matrix assembly
// mat : matrix to fill (standard matrix)
// itM : iterator set to the beginning of the block matrix to add to mat (it is not a matrix but a part of matrix)
// nbu : increment to apply when moving to the next line of block to add
//
//   assemblyMat functions support both omp and no omp computations, but when omp is available omp critical command is used to prevent data races
//   assemblyMatNoCritical do the same with no critical omp command, computation may not be thread safe in omp

template <typename K, typename IteratorM >
inline void assemblyMat(Matrix<K>& mat, IteratorM itM, number_t nbu)
{
#ifdef XLIFEPP_WITH_OMP
  where("assemblyMat(Matrix<K>)");
  error("not_thread_safe");
#endif // XLIFEPP_WITH_OMP
  typename Matrix<K>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      IteratorM it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++)  *itmat +=  *it;
    }
}

//degenerating general scalar case
template <typename K, typename IteratorM >
inline void assemblyMat(K& mat, IteratorM itM, number_t nbu)  //note that nbu=1 in that case
{
#ifdef XLIFEPP_WITH_OMP
  where("assemblyMat(K)");
  error("not_thread_safe");
#endif // XLIFEPP_WITH_OMP
  mat += *itM;
}

//degenerating real scalar case
inline void assemblyMat(real_t& mat, Matrix<real_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
#ifdef XLIFEPP_WITH_OMP
  #pragma omp atomic
#endif // XLIFEPP_WITH_OMP
  mat += *itM;
}

//degenerating real/complex scalar case - forbidden
inline void assemblyMat(real_t& mat, Matrix<complex_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
  error("forbidden","assemblyMat(Real, Matrix<Complex>, Number)");
}

//complex scalar case
inline void assemblyMat(complex_t& mat, Matrix<complex_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
#ifdef XLIFEPP_WITH_OMP
  real_t* mr = reinterpret_cast<real_t*>(&mat), *mc = mr+1;  //trick to move to real in order atomic works
  #pragma omp atomic
  *mr+=itM->real();
  #pragma omp atomic
  *mc+=itM->imag();
#else
  mat += *itM;
#endif // XLIFEPP_WITH_OMP
}

//complex scalar/real case
inline void assemblyMat(complex_t& mat, Matrix<real_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
#ifdef XLIFEPP_WITH_OMP
  real_t* mr = reinterpret_cast<real_t*>(&mat);   //trick to move to real in order atomic works
  #pragma omp atomic
  *mr += *itM;
#else
  mat += *itM;
#endif // XLIFEPP_WITH_OMP
}

//real matrix case
inline void assemblyMat(Matrix<real_t>& mat, Matrix<real_t>::iterator itM, number_t nbu)
{
  Matrix<real_t>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      Matrix<real_t>::iterator it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++)
        {
#ifdef XLIFEPP_WITH_OMP
          #pragma omp atomic
#endif // XLIFEPP_WITH_OMP
          *itmat +=  *it;
        }
    }
}

//complex matrix case
inline void assemblyMat(Matrix<complex_t>& mat, Matrix<complex_t>::iterator itM, number_t nbu)
{
  Matrix<complex_t>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      Matrix<complex_t>::iterator it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++)
        {
#ifdef XLIFEPP_WITH_OMP
          real_t* mr = reinterpret_cast<real_t*>(&*itmat), *mc = mr+1;  //trick to move to real in order atomic works
          #pragma omp atomic
          *mr+=it->real();
          #pragma omp atomic
          *mc+=it->imag();
#else
          *itmat +=  *it;
#endif // XLIFEPP_WITH_OMP
        }
    }
}

//complex/real matrix case
inline void assemblyMat(Matrix<complex_t>& mat, Matrix<real_t>::iterator itM, number_t nbu)
{
  Matrix<complex_t>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      Matrix<real_t>::iterator it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++)
        {
#ifdef XLIFEPP_WITH_OMP
          real_t* mr = reinterpret_cast<real_t*>(&*itmat);  //trick to move to real in order atomic works
          #pragma omp atomic
          *mr+=*it;
#else
          *itmat +=  *it;
#endif // XLIFEPP_WITH_OMP
        }
    }
}

//real/complex matrix case
inline void assemblyMat(Matrix<real_t>& mat, Matrix<complex_t>::iterator itM, number_t nbu)
{
  error("forbidden","assemblyMat(Matrix<Real>,Matrix<Complex>,Number)");
}

// same assemblyMat with no critical omp, should be used when no data races
// ========================================================================

template <typename K, typename IteratorM >
inline void assemblyMatNoCritical(Matrix<K>& mat, IteratorM itM, number_t nbu)
{
  typename Matrix<K>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      IteratorM it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++)  *itmat +=  *it;
    }
}

//degenerating general scalar case
template <typename K, typename IteratorM >
inline void assemblyMatNoCritical(K& mat, IteratorM itM, number_t nbu)  //note that nbu=1 in that case
{
  mat += *itM;
}

//degenerating real scalar case
inline void assemblyMatNoCritical(real_t& mat, Matrix<real_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
  mat += *itM;
}

//degenerating real/complex scalar case - forbidden
inline void assemblyMatNoCritical(real_t& mat, Matrix<complex_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
  error("forbidden","assemblyMatNoCritical(Real, Matrix<Complex>, Number)");
}

//complex scalar case
inline void assemblyMatNoCritical(complex_t& mat, Matrix<complex_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
  mat += *itM;
}

//complex scalar/real case
inline void assemblyMatNoCritical(complex_t& mat, Matrix<real_t>::iterator itM, number_t nbu)  //note that nbu=1 in that case
{
  mat += *itM;
}

//real matrix case
inline void assemblyMatNoCritical(Matrix<real_t>& mat, Matrix<real_t>::iterator itM, number_t nbu)
{
  Matrix<real_t>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      Matrix<real_t>::iterator it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++)
        {
          *itmat +=  *it;
        }
    }
}

//complex matrix case
inline void assemblyMatNoCritical(Matrix<complex_t>& mat, Matrix<complex_t>::iterator itM, number_t nbu)
{
  Matrix<complex_t>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      Matrix<complex_t>::iterator it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++) *itmat +=  *it;
    }
}

//complex/real matrix case
inline void assemblyMatNoCritical(Matrix<complex_t>& mat, Matrix<real_t>::iterator itM, number_t nbu)
{
  Matrix<complex_t>::iterator itmat = mat.begin();
  for(dimen_t i = 0; i < mat.numberOfRows(); i++, itM += nbu)
    {
      Matrix<real_t>::iterator it = itM;
      for(dimen_t j = 0; j < mat.numberOfColumns(); j++, itmat++, it++) *itmat +=  *it;
    }
}

//real/complex matrix case
inline void assemblyMatNoCritical(Matrix<real_t>& mat, Matrix<complex_t>::iterator itM, number_t nbu)
{
  error("forbidden","assemblyMat(Matrix<Real>,Matrix<Complex>,Number)");
}

/* ================================================================================================
                                    FE computation algorithm
   ================================================================================================ */
//FE computation of a SuBilinearForm on a == unique domain ==
// u space and v space have to be FE spaces
// sublf : a single unknown bilinear form defined on a unique domain (assumed)
// mat    : reference to matrix entries of type T
// vt     : to pass as template the scalar type (real or complex)
template <>
template <typename T, typename K>
void ComputationAlgorithm<_FEComputation>::compute(const SuBilinearForm& sublf, LargeMatrix<T>& mat, K& vt,
    Space* space_u_p, Space* space_v_p, const Unknown* u_p, const TestFct* v_p)
{
  if(sublf.size() == 0) return;  // void bilinear form, nothing to compute !!!
  trace_p->push("SuTermMatrix::computeFE");

  if(theVerboseLevel>0) std::cout<<"computing FE term "<< sublf.asString()<<", using "<<numberOfThreads()<<" threads : "<<std::flush;

  //retry spaces
  cit_vblfp it = sublf.begin();
  const Space* sp_u = sublf.up()->space();  //space of unknown
  const Space* sp_v = sublf.vp()->space();  //space of test function
  if(sp_u->typeOfSpace() != _feSpace) error("not_fe_space_type", sp_u->name());
  if(sp_v->typeOfSpace() != _feSpace) error("not_fe_space_type", sp_v->name());
  bool same_interpolation = (sp_u == sp_v);
  const GeomDomain* dom = it->first->asIntgForm()->domain(); //common domain of integrals
  const MeshDomain* mdom = dom->meshDomain();                //common domain of integrals as mesh domain

  //domain info
  bool same_eltshape = (dom->meshDomain()->shapeTypes.size()==1);

  //retry numbering stuff
  //space_u_p and space_v_p are the largest subspaces involved in the computation, may be different from the whole spaces
  Space* subsp_u = Space::findSubSpace(dom, space_u_p);     //find subspace (linked to domain dom) of space_u_p
  if(subsp_u == 0) subsp_u = space_u_p;                     //is a FESpace, not a FESubspace
  Space* subsp_v = subsp_u;
  if(space_u_p!=space_v_p)  //not same interpolation or different parent spaces
    {
      subsp_v = Space::findSubSpace(dom, space_v_p);        //find subspace (linked to domain dom) of space_v_p
      if(subsp_v == 0) subsp_v = space_v_p;                 //is a FESpace, not a FESubspace
      subsp_v->buildgelt2elt();                            //gelt to elt map is useful when element numberings are different
    }
  if(!same_interpolation && subsp_u->nbOfElements() != subsp_v->nbOfElements())  //consistency
    error("term_mismatch_nb_components", subsp_u->nbOfElements(), subsp_v->nbOfElements());
  bool doflag_u = (subsp_u == space_u_p), doflag_v = (subsp_v == space_v_p);
  if(!doflag_u) space_u_p->builddofid2rank();               //build map dofid_u->rank
  if(!doflag_v) space_v_p->builddofid2rank();               //build map dofid_v->rank

  //global information
  bool invertJacobian = false;
  RefElement* relt_u = subsp_u->element_p(number_t(0))->refElt_p;
  RefElement* relt_v = subsp_v->element_p(number_t(0))->refElt_p;
  FEMapType femt_u=relt_u->mapType, femt_v=relt_v->mapType;
  if(femt_u==_covariantPiolaMap || femt_v==_covariantPiolaMap) invertJacobian = true;
  //dimen_t dimfun_u = sp_u->dimFun(), dimfun_v = sp_v->dimFun();  //not correct when trace of vector space
  dimen_t dimfun_u=relt_u->dimShapeFunction;
  dimen_t dimfun_v=relt_v->dimShapeFunction;
  dimen_t nbc_u = u_p->nbOfComponents(), nbc_v = v_p->nbOfComponents(); //number of components of unknown
  if(nbc_u > 1) dimfun_u = nbc_u;                                       //case of vector dofs (scalar shape functions will be extended to vector shape functions)
  if(nbc_v > 1) dimfun_v = nbc_v;                                       //case of vector dofs (scalar shape functions will be extended to vector shape functions)
  bool same_shv = (same_interpolation && nbc_v==nbc_u);
  bool sym = same_shv && !(mat.sym == _noSymmetry);


  //set computation info regarding each bilinear forms
  dimen_t ord_u = 0, ord_v = 0;
  bool nor = false;//, extension = false;
  std::set<Quadrature*> quads; //list of all quadratures required
  std::vector<std::pair<IntgBilinearForm, complex_t> > sublf_copy;    //full copy of subilinear form in order to be thread safe when updating parameters in function
  for(cit_vblfp itf = sublf.begin(); itf != sublf.end(); itf++)
    {
      const IntgBilinearForm* ibf = itf->first->asIntgForm();      //basic bilinear form as simple integral
      const OperatorOnUnknown& op_u = ibf->opu();                  //operator on unknown
      const OperatorOnUnknown& op_v = ibf->opv();                  //operator on test function
      ord_u = std::max(ord_u, dimen_t(op_u.diffOrder()));          //order of u differential involved in operators
      ord_v = std::max(ord_v, dimen_t(op_v.diffOrder()));          //order of v differential involved in operators
      if(op_u.normalRequired() || op_v.normalRequired()) nor=true;  // normal vector required
      if(ord_u >0 || ord_v >0) invertJacobian = true;
      setDomain(const_cast<GeomDomain*>(dom));    //transmit domain pointer to thread data
      std::list<Quadrature*> quadl=ibf->intgMethod()->quadratures();
      quads.insert(quadl.begin(),quadl.end());
      sublf_copy.push_back(std::pair<IntgBilinearForm, complex_t>(*ibf,itf->second));
      //if(op_u.hasExtension() || op_v.hasExtension()) extension = true;
    }

  //precompute , geometry stuff on each element, shapevalues on ref element, no mapping of quadrature points
  preComputationFE(subsp_u, mdom, quads, nbc_u, ord_u>0, invertJacobian, nor, false, false);       //
  if(!same_interpolation) preComputationFE(subsp_v, mdom, quads, nbc_v, ord_v>0, invertJacobian, nor, false, false);    //use structure refElt->qshvs
  else if(nbc_v!=nbc_u)   preComputationFE(subsp_v, mdom, quads, nbc_v, ord_v>0, invertJacobian, nor, false, true);     //use auxiliary structure refElt->qshvs_aux

  //update domain attribute
  mdom->jacobianComputed=true;
  mdom->diffEltComputed=true;
  if(invertJacobian) mdom->inverseJacobianComputed=true;
  if(nor) mdom->normalComputed=true;

  const QuadratureIM* qim = dynamic_cast<const QuadratureIM*>(sublf.begin()->first->asIntgForm()->intgMethod());
  const Element* elt_u = subsp_u->element_p(number_t(0)), *elt_v = subsp_v->element_p(number_t(0));
  Quadrature* quad = qim->getQuadrature(elt_u->geomElt_p->shapeType());
  number_t nbquad = quad->numberOfPoints();
  std::vector<ShapeValues>* shv_u = &elt_u->refElt_p->qshvs_[quad];
  std::vector<ShapeValues>* shv_v = shv_u;
  if(!same_interpolation) shv_v = &elt_v->refElt_p->qshvs_[quad];
  else if(nbc_u!=nbc_v)   shv_v = &elt_v->refElt_p->qshvs_aux[quad];
  bool oneIntg = (sublf.size()==1);

  //temporary structures
  ShapeValues svt_u, svt_v, svt_g;
  Vector<K> val_u, val_v;
  number_t nbelt = subsp_u->nbOfElements();
  ExtensionData extdata;

  //elapsedTime("to prepare FE computation");

  //print status
  number_t nbeltdiv10 = nbelt/10;
  bool show_status = (theVerboseLevel > 0 && mat.nbRows > 1000 &&  nbeltdiv10 > 1);

  //------------------------------------------
  //main loop on finite elements of fesubspace
  //------------------------------------------
#ifdef XLIFEPP_WITH_OMP
  #pragma omp parallel for firstprivate(val_u, val_v, svt_u, svt_v, svt_g, shv_u, shv_v, elt_u, elt_v, sublf_copy, qim, quad, nbquad, extdata)
#endif // XLIFEPP_WITH_OMP
  for(number_t k = 0; k < nbelt; k++)
    {
//      thePrintStream<<"=============== FE compute on element "<<k+1<<" ======================="<<eol<<std::flush;

      // data related to u
      elt_u = subsp_u->element_p(k);
      RefElement* relt_u = elt_u->refElt_p;       //reference element
      GeomElement* gelt = elt_u->geomElt_p;       //geometric element
      ShapeType sh = gelt->shapeType();           //shape type of element
      bool linearMap = gelt->meshElement()->linearMap;
      std::vector<number_t>* dofNum_u = 0;
      if(doflag_u) dofNum_u = const_cast<std::vector<number_t>* >(&subsp_u->elementDofs(k));  //dof numbers (local numbering) of element for u
      else
        {
          dofNum_u = new std::vector<number_t>(elt_u->dofNumbers.size());
          ranks(space_u_p->dofid2rank(), elt_u->dofNumbers, *dofNum_u);            //use ranks function with map dofid2rank (faster)
        }
      number_t nb_u = dofNum_u->size();
//      thePrintStream<<*elt_u<<eol<<std::flush;

      // data related to v
      elt_v = elt_u;
      RefElement* relt_v = relt_u ;                //reference element
      std::vector<number_t>* dofNum_v = dofNum_u;  //dof numbers of element for v
      number_t nb_v = nb_u;
      GeomElement* gelt_v = gelt;
      if(space_u_p!=space_v_p)  //find element of space_v having gelt as geometrical support
        {
          number_t kv = subsp_v->numElement(gelt);
          elt_v = subsp_v->element_p(kv); //try with gelt
          if(elt_v->geomElt_p!=gelt)
            {
              gelt_v = gelt->twin_p;
              if(gelt_v!=0)//try with gelt->twin_p
                {
                  kv = subsp_v->numElement(gelt_v);
                  elt_v = subsp_v->element_p(kv);
                  if(elt_v->geomElt_p!=gelt_v) error("geoelt_not_found");
                }
            }
          relt_v = elt_v->refElt_p;       //reference element
          if(doflag_v) dofNum_v = const_cast<std::vector<number_t>* >(&subsp_v->elementDofs(kv));  //dof numbers (local numbering) of element for v
          else
            {
              dofNum_v = new std::vector<number_t>(elt_v->dofNumbers.size());
              ranks(space_v_p->dofid2rank(), elt_v->dofNumbers, *dofNum_v);            //use ranks function with map dofid2rank (faster)
            }
          nb_v = dofNum_v->size(); //number of DoFs for v
        }
//      thePrintStream<<*elt_v<<eol<<std::flush;

      // local to global numbering
      std::vector<number_t> adrs(nb_u * nb_v, 0);
      mat.positions(*dofNum_v, *dofNum_u, adrs, true); //adresses in mat of elementary matrix stored as a row vector
      std::vector<number_t>::iterator itadrs = adrs.begin(), itk;
      number_t nb_ut = nb_u * nbc_u, nb_vt = nb_v * nbc_v;

      //dof sign correction for div or rot element
      Vector<real_t>* sign_u=0, *sign_v=0;
      bool changeSign_u=false, changeSign_v=false;
      if(relt_u->dofCompatibility == _signDofCompatibility)
        {
          sign_u=&elt_u->getDofSigns();
          changeSign_u = sign_u->size()!=0;
        }
      if(same_shv)
        {
          changeSign_v=changeSign_u;
          sign_v=sign_u;
        }
      else if(relt_v->dofCompatibility == _signDofCompatibility)
        {
          sign_v=&elt_v->getDofSigns();
          changeSign_v = sign_v->size()!=0;
        }

      // dof rotation flag
      bool rotsh_u=elt_u->refElt_p->rotateDof;
      bool rotsh_v=elt_v->refElt_p->rotateDof;

      //geometric stuff
      MeshElement* melt = gelt->meshElement();
      GeomMapData* mapdata= melt->geomMapData_p;
      const RefElement* relt_g = melt->refElement();    //reference element of geom element, may differ from relt_u
      Vector<real_t>* np =0;
      if(nor) np = &mapdata->normalVector;

      // if geom elements supporting finite elements are not the same, build the the additional map
      //      Gv = (Fv)^-1 F  with F : relt_g -> gelt and Fv : relt_g -> gelt_v
      //      Gv(M) = Av * M + Bv  with  Av = Jv^{-1} * J and  Bv = Jv^{-1}*(F(0)-Fv(0))
      // Gv maps relt_g to relt_g. If gelt=gelt_v, then Fv=F and Gv=Id
      // this map is required only for side elements (1d or 2d),  regarding the order of vertices there are only few configurations
      // (2 for a segment, 6 for a triangle, 24 for a quadrangle), so these maps should be precomputed
      // Note : linear map can be used even the gelt/gelt_v are curvilinear!
      Matrix<real_t> Av;
      Vector<real_t> Bv;
      MeshElement* melt_v=0;
      GeomMapData* mapdata_v=0;
      if(gelt_v!=gelt)
        {
          melt_v = gelt_v->meshElement();
          if(melt_v==0) melt_v=gelt_v->buildSideMeshElement();
          mapdata_v= melt_v->geomMapData_p;
          if(mapdata_v==0) {mapdata_v=new GeomMapData(melt_v); melt_v->geomMapData_p = mapdata_v;}
          std::vector<real_t> O(relt_g->dim(),0.);
          Point R=mapdata->geomMap(O), Rv=mapdata_v->geomMap(O);
          mapdata_v->computeJacobianMatrix();
          mapdata_v->computeJacobianDeterminant();  //compute differential element at quadrature point
          mapdata_v->invertJacobianMatrix();        //compute inverse of jacobian
          Bv= mapdata_v->inverseJacobianMatrix*(R-Rv);
          Av= mapdata_v->inverseJacobianMatrix*mapdata->jacobianMatrix;
          //compute v-shape values at quadrature points AV*xq+Bv
          std::set<Quadrature*>::iterator itq=quads.begin();
          std::map<Quadrature*,std::vector<ShapeValues> >::iterator itmq;
          for(; itq!=quads.end(); ++itq)
            {
              itmq =relt_v->qshvs_.find(*itq);
              std::vector<ShapeValues>& shs = itmq->second;
              dimen_t dq = (*itq)->dim();                             //dimension of quadrature points
              number_t nq = (*itq)->numberOfPoints();                 //number of quadrature points
              relt_v->qshvs_[*itq]=std::vector<ShapeValues>(nq);      //allocate shapevalues at quadrature points
              std::vector<ShapeValues>::iterator itshv = relt_v->qshvs_[*itq].begin();
              std::vector<real_t>::const_iterator itp = (*itq)->point();
              for(number_t i = 0; i < nq; i++, itp += dq, itshv++) //compute the shape values at the mapped quadrature point
                {
                  Vector<real_t> x(itp,itp+dq);
                  x=Av*x+Bv;  // map quadrature point
                  relt_v->computeShapeValues(x.begin(), ord_v>0);
                  *itshv = relt_v->shapeValues;
                  if(nbc_v>1) itshv->extendToVector(nbc_v);    //extend scalar shape functions to nbc vector shape functions
                }
            }
        }

      //loop on basic bilinear forms
      //----------------------------
      std::vector<std::pair<IntgBilinearForm, complex_t> > ::iterator itf;
      for(itf = sublf_copy.begin(); itf != sublf_copy.end(); itf++)
        {
          const IntgBilinearForm* ibf = &itf->first;           // basic bilinear form as simple integral
          const OperatorOnUnknown* op_u = ibf->opu_p();        // operator on unknown
          const OperatorOnUnknown* op_v = ibf->opv_p();        // operator on test function
          const Extension* ext_u = op_u->extension(), *ext_v = op_v->extension();
          if((ext_u==0 || (ext_u!=0 && ext_u->domToSide.find(gelt)!=ext_u->domToSide.end())) &&
              (ext_v==0 || (ext_v!=0 && ext_v->domToSide.find(gelt)!=ext_v->domToSide.end())))
            //do computation if no extension or element in extension
            {
//            thePrintStream<<"   compute "; ibf->print(thePrintStream);
              bool id_u = op_u->isId();                            // opu is Id
              bool id_v = op_v->isId();                            // opv is Id
              number_t ord_opu = op_u->diffOrder();                // derivative order of op_u
              number_t ord_opv = op_v->diffOrder();                // derivative order of op_v
              bool hasf = op_u->hasFunction() || op_v->hasFunction();   // true if data function
              bool same_val = (op_u->difOpType() == op_v->difOpType()
                               && !op_u->hasOperand() && !op_v->hasOperand()
                               && same_shv);
              bool mapsh_u = (femt_u!=_standardMap || ord_opu>0);
              bool mapsh_v = (femt_v!=_standardMap || ord_opv>0);
              AlgebraicOperator aop = ibf->algop();                  // algebraic operator between differential operator
              K coef = complexToT<K>(itf->second);                   // coefficient of linear form
              qim = dynamic_cast<const QuadratureIM*>(ibf->intgMethod());
              if(!oneIntg || !same_eltshape)   //update quadrature
                {
                  quad = qim->getQuadrature(sh);
                  nbquad = quad->numberOfPoints();
                  shv_u = &relt_u->qshvs_[quad];
                  shv_v=shv_u;
                  if(!same_interpolation) shv_v = &elt_v->refElt_p->qshvs_[quad];
                  else if(nbc_u!=nbc_v)   shv_v = &elt_v->refElt_p->qshvs_aux[quad];
                }
              if(op_u->elementRequired() || op_v->elementRequired()) setElement(gelt);  // transmit gelt to thread data
              std::vector<ShapeValues>::iterator itsh_u=shv_u->begin(), itsh_v=shv_v->begin();

              //compute integrand on quadrature points
              Matrix<K> matel(nb_vt, nb_ut, K(0));
              for(number_t q = 0; q < nbquad; q++, ++itsh_u, ++itsh_v)    //loop on quadrature points
                {
//                  thePrintStream<<"-------------- q="<<q<<" ----------------"<<eol<<std::flush;
                  if(!linearMap)   //update geometric stuff on each quadrature points
                    {
                      svt_g.resize(relt_g->nbDofs(),relt_g->dim());
                      relt_g->computeShapeValues(quad->point(q), svt_g);   //compute geometric shape values at quadrature point
                      mapdata->computeJacobianMatrix(svt_g);
                      mapdata->computeJacobianDeterminant();               //compute differential element at quadrature point
                      if(invertJacobian) mapdata->invertJacobianMatrix();   //compute inverse of jacobian
                      if(nor)
                        {
                          mapdata->computeOutwardNormal();                 //compute normal vector
                          np = &mapdata->normalVector;
                        }
                    }
                  //map u-shapevalues in physical space if required
                  ShapeValues* sv_u = &*itsh_u;
                  dimen_t dimfunp_u = dimfun_u;
                  if(!relt_u->hasShapeValues) //compute shape values in physical space using dof
                    {
                      Point x = mapdata->geomMap(Point(quad->point(q), quad->dim()));
                      svt_u=elt_u->computeShapeValues(x,ord_opu>0);
                      sv_u=&svt_u;
                    }
                  else
                    {
                      if(mapsh_u || changeSign_u || rotsh_u)
                        {
                          number_t ord=ord_opu;
                          if(same_shv) ord=std::max(ord_opu,ord_opv);
                          mapShapeValues(*relt_u, *melt, *mapdata, mapsh_u, femt_u, rotsh_u,
                                         ord, changeSign_u, sign_u, dimfun_u, dimfunp_u, *itsh_u, svt_u);
                          sv_u=&svt_u;
                        }
                    }
//                  thePrintStream<<"mapped shv_u="<<*sv_u<<" sign_u="<<*sign_u<<eol<<std::flush;
                  //map v-shapevalues in physical space if required
                  ShapeValues* sv_v = sv_u;
                  dimen_t dimfunp_v = dimfun_v;
                  if(!relt_v->hasShapeValues) //compute shape values in physical space using dof
                    {
                      if(!same_shv)
                        {
                          Point x = mapdata->geomMap(Point(quad->point(q), quad->dim()));
                          svt_v=elt_v->computeShapeValues(x,ord_opv>0);
                          sv_v=&svt_v;
                        }
                    }
                  else
                    {
                      if(!same_shv) sv_v = &*itsh_v;
                      if(mapsh_v || changeSign_v || rotsh_v) // v-shapevalues are not the same
                        {
                          if(gelt==gelt_v)
                            mapShapeValues(*relt_v, *melt, *mapdata, mapsh_v, femt_v, rotsh_v,
                                           ord_opv, changeSign_v, sign_v, dimfun_v, dimfunp_v, *itsh_v, svt_v);
                          else
                            mapShapeValues(*relt_v, *melt_v, *mapdata_v, mapsh_v, femt_v, rotsh_v,
                                           ord_opv, changeSign_v, sign_v, dimfun_v, dimfunp_v, *itsh_v, svt_v);
                          sv_v=&svt_v;
                        }
                    }
//                  thePrintStream<<"mapped shv_v="<<*sv_v<<" sign_v="<<*sign_v<<eol<<std::flush;
                  K alpha = coef * mapdata->differentialElement** (quad->weight(q));

                  //compute integrand (OperatorOnUnknown)
                  dimen_t du,mu,dv,mv; //to store block sizes
                  if(same_val)  //same values in u and v (and no value or function) - faster computation
                    {
                      if(id_u)   //shortcut, bypass operator evaluate
                        {
                          tensorOpAdd(aop, sv_u->w, nb_ut, sv_u->w, nb_ut, matel, alpha);  //elementary matrix, note the transposition of u and v
                        }
                      else
                        {
                          op_u->eval(sv_u->w, sv_u->dw, dimfunp_u, val_u, du, mu, np);      //evaluate differential operator
                          tensorOpAdd(aop, val_u, nb_ut, val_u, nb_ut, matel, alpha);      //elementary matrix, note the transposition of u and v
                        }
                    }
                  else //general case
                    {
                      if(!hasf)   //no function to evaluate
                        {
                          if(id_u)
                            {
                              if(id_v) tensorOpAdd(aop, sv_v->w, nb_vt, sv_u->w, nb_ut, matel, alpha);   //elementary matrix, note the transposition of u and v
                              else
                                {
                                  op_v->eval(sv_v->w, sv_v->dw, dimfunp_v, val_v, dv, mv, np);          //evaluate differential operator
                                  tensorOpAdd(aop, val_v, nb_vt, sv_u->w, nb_ut, matel, alpha);         //elementary matrix, note the transposition of u and v
                                }
                            }
                          else
                            {
                              op_u->eval(sv_u->w, sv_u->dw, dimfunp_u, val_u, du, mu, np);              //evaluate differential operator
                              if(id_v) tensorOpAdd(aop, sv_v->w, nb_vt, val_u, nb_ut, matel, alpha);     //elementary matrix, note the transposition of u and v
                              else
                                {
                                  op_v->eval(sv_v->w, sv_v->dw, dimfunp_v, val_v, dv, mv, np);          //evaluate differential operator
                                  tensorOpAdd(aop, val_v, nb_vt, val_u, nb_ut, matel, alpha);           //elementary matrix, note the transposition of u and v
                                }
                            }
                        }
                      else // general case
                        {
                          Point xq(quad->point(q), quad->dim()), x;
                          x = mapdata->geomMap(xq);
                          if(ext_u)
                            {
                              extdata.compute(op_u->extension(), gelt, xq.begin());  //update the extension data
                              op_u->eval(x, sv_u->w, sv_u->dw, dimfunp_u, val_u, du, mu, np, &extdata);
                            }
                          else
                            {
                              if(op_u->hasFunction()) op_u->eval(x, sv_u->w, sv_u->dw, dimfunp_u, val_u, du, mu, np);   //evaluate differential operator with function
                              else op_u->eval(sv_u->w, sv_u->dw, dimfunp_u, val_u, du, mu, np);     //evaluate differential operator
                            }
                          if(ext_v)
                            {
                              extdata.compute(op_v->extension(), gelt, xq.begin());  //update the extension data
                              op_v->eval(x, sv_v->w, sv_v->dw, dimfunp_v, val_v, dv, mv, np, &extdata);   //evaluate differential operator with function
                            }
                          else
                            {
                              if(op_v->hasFunction()) op_v->eval(x, sv_v->w, sv_v->dw, dimfunp_v, val_v, dv, mv, np);   //evaluate differential operator with function
                              else op_v->eval(sv_v->w, sv_v->dw, dimfunp_v, val_v, dv, mv, np);      //evaluate differential operator
                            }
                          tensorOpAdd(aop, val_v, nb_vt, val_u, nb_ut, matel, alpha);           //elementary matrix, note the transposition of u and v
//                          thePrintStream<<" val_u="<<val_u<<"val_v="<<val_v<<eol;
                        }
                    }
                }//end of quadrature points loop

//              thePrintStream<<"dof u="<<*dofNum_u<<" dof v="<<*dofNum_v<<eol;
//              thePrintStream<<"adrs="<<adrs<<eol<<"matel="<<matel<<eol<<std::flush;

              //assembling matel in global matrix
              //matel is never a block matrix as mat is a block matrix in case of vector unknowns
              std::vector<number_t>::iterator itd_u, itd_v, it_ue;   //dof numbering iterators
              typename Matrix<K>::iterator itm = matel.begin(), itm_u;
              number_t incr = nbc_u * nbc_v * nb_u; //block row increment in matel
              number_t i = 0;
              //loop on (vector) dofs in u and v
              for(itd_v = dofNum_v->begin(); itd_v != dofNum_v->end() ; itd_v++, itm += incr, i++)
                {
                  itm_u = itm;
                  itk = itadrs + i * nb_u;
                  for(itd_u = dofNum_u->begin(); itd_u != dofNum_u->end() ; itd_u++, itm_u += nbc_u, itk++)
                    {
                      if(!sym || *itd_v >= *itd_u)
                        {
                          assemblyMat(mat(*itk), itm_u, nb_ut);
                        }
                    }
                }
            }
        } //end of bilinear forms loop

      if(!doflag_u) delete dofNum_u; //clean temporary dofNum structure
      if(!same_interpolation && !doflag_v) delete dofNum_v;

      if(show_status)  //progress status
        {
          number_t numthread=0;
#ifdef XLIFEPP_WITH_OMP
          numthread = omp_get_thread_num();
#endif // XLIFEPP_WITH_OMP
          if(numthread==0 && k!=0 && k % nbeltdiv10 == 0) { std::cout<< k/nbeltdiv10 <<"0% "<<std::flush; }
        }
    } //end of elements loop

  if(theVerboseLevel > 0) std::cout<<" done"<<eol<<std::flush;
  trace_p->pop();
}

} // end of namespace xlifepp