term/computation/FESPMatrixComputation.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 FESPMatrixComputation.hpp
  \author E. Lunéville
  \since 9 jan 2013
  \date 25 juin 2013

  \brief Implementation of template FE-SP TermMatrix computation functions

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


/* ================================================================================================
                                      FE-SP computation algorithm
   ================================================================================================ */

namespace xlifepp
{

/*! FE-SP hybrid computation of a scalar SuBiinearForm on a == unique domain ==
   ---- u space is a FE/SP space and v space is a SP/FE space ---------
   sublf : a single unknown bilinear form defined on a unique domain (assumed), single integral
   v     : 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<_FESPComputation>::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::computeFESP");

   if(theVerboseLevel>0) std::cout<<"computing FE-SP 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
   const Space* fe_sp=0, *sp_sp=0;
   Space * space_fe=0;
   const Unknown* u_fe=0, *u_sp=0;
   bool uisfe=true;
   if( u_p->isFE())
   {
      fe_sp=sp_u;
      space_fe=space_u_p;
      u_fe=u_p;
      if(v_p->isSpectral())
      {
         sp_sp=sp_v;
         u_sp=v_p;
      }
   }
   else if( u_p->isSpectral())
   {
      sp_sp=sp_u;
      u_sp=u_p;
      uisfe=false;
      if(v_p->isFE())
      {
         fe_sp=sp_v;
         space_fe=space_v_p;
         u_fe=v_p;
      }
   }
   if( fe_sp==0 || sp_sp ==0 ) error("not_fe_sp_pair");

   const GeomDomain* dom = it->first->asIntgForm()->domain();  //common domain of integrals
   Space* subsp_fe = fe_sp->findSubSpace(dom,space_fe);        //find subspace (linked to domain dom) of space_fe
   if(subsp_fe == 0) subsp_fe = space_fe;                      //is a FESpace, not a FESubspace

   const Function* mapTo=findMap(*dom,*sp_sp->domain());

   bool doflag = (subsp_fe == space_fe);
   dimen_t dimfun_fe = fe_sp->dimFun(), dimfun_sp = sp_sp->dimFun();
   dimen_t nbc_fe = u_fe->nbOfComponents(), nbc_sp = u_sp->nbOfComponents(); //number of components of unknowns
   if(nbc_fe > 1) dimfun_fe = nbc_fe;                                      //case of vector dofs (scalar shape functions will be extended to vector shape functions)
   if(nbc_sp > 1) dimfun_sp = nbc_sp;                                      //case of vector dofs (scalar shape functions will be extended to vector shape functions)
   number_t nb_sp = sp_sp->dimSpace();
   number_t nb_fe = fe_sp->dimSpace();
   number_t nb_fet= nbc_fe*nb_fe;
   number_t nb_spt= nbc_sp*nb_sp;

   //copy of sublf to be thread safe
   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
   dimen_t ordu = 0, ordv = 0;
   for(cit_vblfp itf = sublf.begin(); itf != sublf.end(); itf++) {
      const IntgBilinearForm* ibf = itf->first->asIntgForm(); //basic bilinear form as simple integral
      ordu = std::max(ordu, dimen_t(ibf->opu().diffOrder()));   //order of u differential involved in operators
      ordv = std::max(ordv, dimen_t(ibf->opv().diffOrder()));   //order of v differential involved in operators
      sublf_copy.push_back(std::pair<IntgBilinearForm, complex_t>(*ibf,itf->second));
   }

   //global information for spectral unknown and trivial numbering
   const SpectralBasis* basis = sp_sp->spSpace()->spectralBasis();    //spectral basis
   std::vector<number_t> dofNum_sp(nb_sp);
   std::vector<number_t>::iterator itd = dofNum_sp.begin();
   for(number_t k = 1; k <= nb_sp; k++, itd++) *itd = k;
   std::vector<number_t> dofNum_fe;

   //temporary data structures to store shape values
   std::map<const Quadrature*, std::vector<ShapeValues> > shapevalues;      //shape values for fe unknown at quadrature points
   bool sym = false;

   //link by reference (u,v) -> (fe,sp)
   dimen_t ord_fe=ordu/*, ord_sp=ordv*/;
   number_t   nb_u=nb_fe,   nbc_u=nbc_fe, nb_ut=nb_fet;
   number_t /*nb_v=nb_sp,*/ nbc_v=nbc_sp, nb_vt=nb_spt;
   Vector<K> val_fe, val_sp;
   if(!uisfe) {
      nb_u=nb_sp,   nbc_u=nbc_sp, nb_ut=nb_spt;
    /*nb_v=nb_fe,*/ nbc_v=nbc_fe, nb_vt=nb_fet;
      ord_fe=ordv; /*ord_sp=ordu;*/
   }

   bool analytic=(basis->funcFormType()==_analytical);
   number_t nbelt = subsp_fe->nbOfElements();

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

   //main loop on finite elements of fe subspace
   #ifdef XLIFEPP_WITH_OMP
    //to be thread safe, because the parameter of analytic basis is updated for each index basis, use copies of analytic basis
    //no copy used when interpolated basis !
    SpectralBasisFun sbf(Function(),0);   //fake spectral basis functions
    if(analytic) sbf=SpectralBasisFun(reinterpret_cast<const SpectralBasisFun&>(*basis));
    #pragma omp parallel firstprivate(sublf_copy, sbf, dofNum_fe, dofNum_sp, shapevalues, val_fe, val_sp) shared(mat)
    #pragma omp for schedule(dynamic) nowait
   #endif
   for(number_t k = 0; k < nbelt; k++)
   {
      // data related to fe unknown
      const Element* elt = subsp_fe->element_p(k);
      RefElement* relt = elt->refElt_p;       //reference element
      GeomElement* gelt = elt->geomElt_p;     //geometric element
      ShapeType sh = gelt->shapeType();       //shape type of element
      if(doflag) dofNum_fe = std::vector<number_t>(subsp_fe->elementDofs(k)); //dof numbers (local numbering) of element
      else dofNum_fe = std::vector<number_t>(subsp_fe->elementParentDofs(k)); //dof numbers (parent numbering) of element
      nb_fe = dofNum_fe.size();
      nb_fet= nbc_fe*nb_fe;
      std::vector<number_t>* dofNum_u, *dofNum_v;
      Vector<K>* val_u, *val_v;
      if(uisfe)
      {
          nb_u=nb_fe; nb_ut=nb_fet;
          dofNum_u=&dofNum_fe;dofNum_v=&dofNum_sp;
          val_u=&val_fe; val_v=&val_sp;
       }
      else
      {/*nb_v=nb_fe;*/ nb_vt=nb_fet;
       dofNum_u=&dofNum_sp;dofNum_v=&dofNum_fe;
        val_u=&val_sp; val_v=&val_fe;}
      //link geomElement to basis functions
      #ifdef XLIFEPP_WITH_OMP
          if(analytic) sbf.setGeomElementPointer(gelt);
      #else
          if(analytic) basis->setGeomElementPointer(gelt);
      #endif

      // local to global numbering
      std::vector<number_t> adrs(nb_fe * nb_sp, 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;
      std::map<const Quadrature*, std::vector<ShapeValues> >::iterator itsh;

      //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();             //operator on unknown
         const OperatorOnUnknown& op_v = ibf->opv();             //operator on test function
         if(op_u.elementRequired() || op_v.elementRequired()) setElement(gelt);  // transmit gelt to thread data
         AlgebraicOperator aop = ibf->algop();                   //algebraic operator between differential operator
         K coef = complexToT<K>(itf->second);                    //coefficient of linear form
         const OperatorOnUnknown* op_fe=&op_u, *op_sp=&op_v;
         if(!uisfe) {op_fe=&op_v; op_sp=&op_u;}
         const QuadratureIM* qim = dynamic_cast<const QuadratureIM*>(ibf->intgMethod());
         const Quadrature* quad = qim->getQuadrature(sh);        //quadrature rule pointer attached to integral

         //compute shape values on quadrature points or get them if already computed
         number_t nbquad = quad->numberOfPoints();
         itsh= shapevalues.find(quad);
         if(itsh == shapevalues.end()) //compute new shape values on quadrature points
         {
            std::vector<ShapeValues> shvs(nbquad);
            for(number_t q = 0; q < nbquad; q++)
            {
//               relt->computeShapeValues(quad->point(q), ord_fe > 0);
               shvs[q] = relt->shapeValues;  //to size shvs[q]
               relt->computeShapeValues(quad->point(q), shvs[q], ord_fe > 0);
               if(nbc_fe > 1) shvs[q].extendToVector(nbc_fe);  //extend scalar shape functions to nbc vector shape functions
            }
            itsh = shapevalues.insert(std::pair<const Quadrature*, std::vector<ShapeValues> >(quad, shvs)).first;
         }

         //compute integrand on quadrature points
         MeshElement* melt = gelt->meshElement();
         if(melt == 0) melt = gelt->buildSideMeshElement(); //it is a side element with no meshElement structure, build one
         GeomMapData mapdata(melt);                         //structure to store jacobian
         bool invertJacobian = (op_u.diffOrder() > 0 || op_v.diffOrder() > 0);

         //loop on quadrature points
         Matrix<K> matel(nb_vt, nb_ut, K(0));
         dimen_t d,m;  //block size
         for(number_t q = 0; q < nbquad; q++)
         {
            mapdata.computeJacobianMatrix(quad->point(q));         //compute jacobian at q th quadrature point
            mapdata.computeDifferentialElement();                  //compute differential element
            ShapeValues sv(itsh->second[q]);                       //transportation of shapevalues
            if(invertJacobian) mapdata.invertJacobianMatrix();     //compute inverse of jacobian
            sv.map(itsh->second[q], mapdata, op_fe->diffOrder());  //shape values in physical space
            K alpha = coef * mapdata.differentialElement **(quad->weight(q));
            Point xq(quad->point(q), quad->dim()), x = mapdata.geomMap(xq);
            Point xs=x;
            if(mapTo != 0) xs = (*mapTo)(x, xs);                                                                         //map physical point to point in the reference space of basis function
            if(op_fe->hasFunction()) op_fe->eval(x, sv.w, sv.dw, dimfun_fe, val_fe, d, m);     //evaluate differential operator with function
            else                     op_fe->eval(sv.w, sv.dw, dimfun_fe, val_fe, d, m);                      //evaluate differential operator
            if(analytic)
                #ifdef XLIFEPP_WITH_OMP
                computeSPOperator(*op_sp, &sbf, xs, sp_sp->valueType(), nb_sp, dimfun_sp, val_sp);//evaluate differential operator for spectral basis
                #else
                computeSPOperator(*op_sp, basis, xs, sp_sp->valueType(), nb_sp, dimfun_sp, val_sp);//evaluate differential operator for spectral basis
                #endif
            else computeSPOperator(*op_sp, basis, xs, sp_sp->valueType(), nb_sp, dimfun_sp, val_sp);//evaluate differential operator for spectral basis
            tensorOpAdd(aop, *val_v, nb_vt, *val_u, nb_ut, matel, alpha);                      //elementary matrix, note the transposition of u and v
         }//end of quadrature points loop

         //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(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(show_status) std::cout<<" done"<<eol<<std::flush;
   trace_p->pop();
}

} // end of namespace xlifepp