SpectralLegendreDiscretizer.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2007 Stéphane Del Pino

//  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 2, 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, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: SpectralLegendreDiscretizer.hpp,v 1.11 2008/03/30 20:24:24 delpinux Exp $

#ifndef SPECTRAL_LEGENDRE_DISCRETIZER_HPP
#define SPECTRAL_LEGENDRE_DISCRETIZER_HPP

#include <Problem.hpp>
#include <DegreeOfFreedomSet.hpp>
#include <Structured3DMesh.hpp>
#include <SpectralMesh.hpp>

#include <SpectralFunction.hpp>
#include <ScalarFunctionBase.hpp>

#include <Interval.hpp>
#include <LegendreBasis.hpp>
#include <GaussLobatto.hpp>

#include <SpectralConformTransformation.hpp>

#include <DoubleHashedMatrix.hpp>
#include <UnAssembledMatrix.hpp>

#include <DiscretizationType.hpp>

#include <Discretization.hpp>
#include <ScalarDiscretizationTypeSpectral.hpp>

#include <PDE.hpp>
#include <PDEProblem.hpp>
#include <MassOperator.hpp>
#include <FirstOrderOperator.hpp>
#include <DivMuGrad.hpp>
#include <SecondOrderOperator.hpp>

#include <PDESystem.hpp>

#include <VariationalProblem.hpp>

#include <VariationalOperatorFV.hpp>
#include <VariationalOperatorFdxV.hpp>
#include <VariationalOperatorFdxGV.hpp>
#include <VariationalOperatorFgradGgradV.hpp>

#include <VariationalOperatorAlphaUV.hpp>
#include <VariationalOperatorMuGradUGradV.hpp>
#include <VariationalOperatorAlphaDxUDxV.hpp>
#include <VariationalOperatorNuUdxV.hpp>
#include <VariationalOperatorNuDxUV.hpp>

#include <Mesh.hpp>

#include <ErrorHandler.hpp>

class SpectralLegendreDiscretizer
{
private:
  const Problem& __problem;
  BaseVector& __b;
  // Set of degrees of freedom
  const DegreeOfFreedomSet& __degreeOfFreedomSet;

  const SpectralMesh& __mesh;

  Vector<TinyVector<3,size_t> >  __dofShape;
  


  TinyVector<3, ConstReferenceCounting<Interval> > __interval;
  Vector<TinyVector<3, ConstReferenceCounting<Interval> > > __intervalU;
  TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > __transform;
  Vector<TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > > __transformU;
  TinyVector<3, ConstReferenceCounting<GaussLobatto> > __gaussLobatto;
  Vector<TinyVector<3, ConstReferenceCounting<LegendreBasis>  > > __basisU;
  
  void _interpolateLagrange(const ScalarFunctionBase& alpha,
                            Structured3DVector<real_t>& alpha_rst) const
  {
    TinyVector<3,Vector<real_t> > nodes;
    for( size_t i=0; i<3; ++i){
      nodes[i].resize((*__gaussLobatto[i]).numberOfPoints());
      for (size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints(); ++j) {
        nodes[i][j] =(*__transform[i])((*__gaussLobatto[i])(j)); 
      }
    }
    for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
      const real_t& x = nodes[0][r];
      for (size_t s=0; s < __gaussLobatto[1]->numberOfPoints(); ++s) {
        const real_t& y = nodes[1][s];
        for (size_t t=0; t < __gaussLobatto[2]->numberOfPoints(); ++t) {
          const real_t& z = nodes[2][t];
          alpha_rst(r,s,t) = alpha(x,y,z);
        }
      }
    }
  }


  void _getAu(const Structured3DVector<real_t>& coef,
              const real_t& detXU, const real_t& detYU, const real_t& detZU,
              const real_t& detXV, const real_t& detYV, const real_t& detZV,
              const Vector<Vector<real_t> >& basisOrDerivativeValueU0,
              const Vector<Vector<real_t> >& basisOrDerivativeValueU1,
              const Vector<Vector<real_t> >& basisOrDerivativeValueU2,
              const Vector<Vector<real_t> >& basisOrDerivativeValueV0,
              const Vector<Vector<real_t> >& basisOrDerivativeValueV1,
              const Vector<Vector<real_t> >& basisOrDerivativeValueV2,
              const size_t& unknownNumber, const size_t& testFunctionNumber,
              const Vector<real_t>& u,
              Vector<real_t>& v) const
  {
    const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
      = __basisU[unknownNumber];
    
    const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& testBasis
      = __basisU[testFunctionNumber];
    
    Structured3DVector<real_t> u_ijt(Array3DShape(unknownBasis[0]->dimension(),
                                                  unknownBasis[1]->dimension(),
                                                  __gaussLobatto[2]->numberOfPoints()));
    
    u_ijt=0;
    
    
    const TinyVector<3,size_t>& dofShapeU = __dofShape[unknownNumber];
    
    {
      for (size_t i=0; i< unknownBasis[0]->dimension(); ++i) {
        for (size_t j=0; j< unknownBasis[1]->dimension(); ++j) {
          for (size_t k=0; k< unknownBasis[2]->dimension(); ++k) {
            size_t l  = i*dofShapeU[1]*dofShapeU[2] + j*dofShapeU[2] + k;
            const real_t& u_ijk = u[__degreeOfFreedomSet(unknownNumber,l)];
            for (size_t t=0; t< __gaussLobatto[2]->numberOfPoints(); ++t) {
              u_ijt(i,j,t) += u_ijk * basisOrDerivativeValueU2[t][k];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_ist(Array3DShape(unknownBasis[0]->dimension(),
                                                  __gaussLobatto[1]->numberOfPoints(),
                                                  __gaussLobatto[2]->numberOfPoints()));
    
    u_ist=0;

    for (size_t i=0; i< unknownBasis[0]->dimension(); ++i) {
      for (size_t j=0; j< unknownBasis[1]->dimension(); ++j) {
        for (size_t t=0; t<__gaussLobatto[2]->numberOfPoints(); ++t) {
          const real_t& u_ijt_value = u_ijt(i,j,t);
          for (size_t s=0; s<__gaussLobatto[1]->numberOfPoints(); ++s) {
            u_ist(i,s,t) += u_ijt_value*basisOrDerivativeValueU1[s][j];
          }
        }
      }
    }
    

    Structured3DVector<real_t> u_rst(Array3DShape(__gaussLobatto[0]->numberOfPoints(),
                                                  __gaussLobatto[1]->numberOfPoints(),
                                                  __gaussLobatto[2]->numberOfPoints()));
    u_rst=0;
    
    for (size_t i=0; i< unknownBasis[0]->dimension(); ++i) {
      for (size_t s=0; s<__gaussLobatto[1]->numberOfPoints(); ++s) {
        for (size_t t=0; t<__gaussLobatto[2]->numberOfPoints(); ++t) {
          const real_t& u_ist_value = u_ist(i,s,t);
          for (size_t r=0; r<__gaussLobatto[0]->numberOfPoints(); ++r) {
            u_rst(r,s,t) += u_ist_value*basisOrDerivativeValueU0[r][i];
          }
        }
      }
    }

    Structured3DVector<real_t> v_rsk(Array3DShape(__gaussLobatto[0]->numberOfPoints(),
                                                  __gaussLobatto[1]->numberOfPoints(),
                                                  testBasis[2]->dimension()));
    v_rsk=0;    


    
    {
      for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
        for (size_t s=0; s < __gaussLobatto[1]->numberOfPoints(); ++s) {
          for (size_t t=0; t < __gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& wt = __gaussLobatto[2]->weight(t);
            const real_t value_rst = coef(r,s,t)*u_rst(r,s,t)
              *(__transform[2]->determinant())*wt;
            for (size_t k=0; k< testBasis[2]->dimension(); ++k) {
              v_rsk(r,s,k) += value_rst * basisOrDerivativeValueV2[t][k];
            }
          }
        }
      }
      
      Structured3DVector<real_t> v_rjk(Array3DShape(__gaussLobatto[0]->numberOfPoints(),
                                                    testBasis[1]->dimension(),
                                                    testBasis[2]->dimension()));
      v_rjk=0;    

      {
        for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
          for (size_t s=0; s < __gaussLobatto[1]->numberOfPoints(); ++s) {
            const real_t& ws = __gaussLobatto[1]->weight(s);
            for (size_t k=0; k < testBasis[2]->dimension(); ++k) {
              const real_t value_rsk = v_rsk(r,s,k)
                *(__transform[1]->determinant())*ws;
              for (size_t j=0; j< testBasis[1]->dimension(); ++j) {
                v_rjk(r,j,k) += value_rsk * basisOrDerivativeValueV1[s][j];
              }
            }
          }
        }
      }
      
      real_t detUV = detXU *  detYU *detZU *detXV *detYV *detZV;  
      
      {
        const TinyVector<3,size_t>& dofShape = __dofShape[testFunctionNumber];
        
        for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
          const real_t& wr = __gaussLobatto[0]->weight(r);
          for (size_t j=0; j < testBasis[1]->dimension(); ++j) {
            for (size_t k=0; k < testBasis[2]->dimension(); ++k) {
              const real_t value_rjk = v_rjk(r,j,k)*(__transform[0]->determinant())*wr*detUV;
              for (size_t i=0; i< testBasis[0]->dimension(); ++i) {
                size_t dofNumber  = i*dofShape[1]*dofShape[2] + j*dofShape[2] + k;
                v[__degreeOfFreedomSet(testFunctionNumber, dofNumber)] +=
                  value_rjk * basisOrDerivativeValueV0[r][i];
              }
            }
          }
        }
      }
    }
  }

  void _getA(const Structured3DVector<real_t>& coef,
             const real_t& detXU, const real_t& detYU, const real_t& detZU,
             const real_t& detXV, const real_t& detYV, const real_t& detZV,
             const Vector<Vector<real_t> >& basisOrDerivativeValueU0,
             const Vector<Vector<real_t> >& basisOrDerivativeValueU1,
             const Vector<Vector<real_t> >& basisOrDerivativeValueU2,
             const Vector<Vector<real_t> >& basisOrDerivativeValueV0,
             const Vector<Vector<real_t> >& basisOrDerivativeValueV1,
             const Vector<Vector<real_t> >& basisOrDerivativeValueV2,
             const size_t& unknownNumber,
             DoubleHashedMatrix& A) const
  {
  }

  void _getD(const Structured3DVector<real_t>& coef,
             const real_t& detXU, const real_t& detYU, const real_t& detZU,
             const real_t& detXV, const real_t& detYV, const real_t& detZV,
             const Vector<Vector<real_t> >& basisOrDerivativeValueU0,
             const Vector<Vector<real_t> >& basisOrDerivativeValueU1,
             const Vector<Vector<real_t> >& basisOrDerivativeValueU2,
             const Vector<Vector<real_t> >& basisOrDerivativeValueV0,
             const Vector<Vector<real_t> >& basisOrDerivativeValueV1,
             const Vector<Vector<real_t> >& basisOrDerivativeValueV2,
             const size_t& unknownNumber,
             Vector<real_t>& v) const
  {
    const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
      = __basisU[unknownNumber];
    

    TinyVector<3, Vector<real_t> > diag;
    for (size_t i=0; i< 3; ++i) {
      diag[i].resize(__gaussLobatto[i]->numberOfPoints());
      diag[i] =0;
    }
    
    for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
      const real_t& wr = __gaussLobatto[0]->weight(r);
      for (size_t i=0; i < unknownBasis[0]->dimension(); ++i) {
        diag[0][i] += basisOrDerivativeValueU0[r][i]* basisOrDerivativeValueV0[r][i]* wr;
      }
    }

    for (size_t s=0; s < __gaussLobatto[1]->numberOfPoints(); ++s) {
      const real_t& ws = __gaussLobatto[1]->weight(s);
      for (size_t j=0; j < unknownBasis[1]->dimension(); ++j) {
        diag[1][j] += basisOrDerivativeValueU1[s][j]* basisOrDerivativeValueV1[s][j]* ws;
      }
    }
    
    for (size_t t=0; t < __gaussLobatto[2]->numberOfPoints(); ++t) {
      const real_t& wt = __gaussLobatto[2]->weight(t);
      for (size_t k=0; k < unknownBasis[2]->dimension(); ++k) {
        diag[2][k] += basisOrDerivativeValueU2[t][k]* basisOrDerivativeValueV2[t][k]* wt;
      }
    }
    const TinyVector<3,size_t>& dofShapeU = __dofShape[unknownNumber];
    const real_t det = detXU *  detYU *detZU *detXV *detYV *detZV 
      * __transform[0]->determinant() 
      * __transform[1]->determinant()
      * __transform[2]->determinant() ;

    for (size_t i=0; i < unknownBasis[0]->dimension(); ++i){
      const real_t  diag1= diag[0][i]*det;
      for (size_t j=0; j < unknownBasis[1]->dimension(); ++j){
        const real_t diag2 = diag1 *diag[1][j];
        for (size_t k=0; k < unknownBasis[2]->dimension(); ++k){
          size_t l  = i*dofShapeU[1]*dofShapeU[2] + j*dofShapeU[2] + k;
          v[__degreeOfFreedomSet(unknownNumber,l)] += diag2*diag[2][k];
        }
      }
    }
  }

public:  
  void timesX(const BaseVector& U, BaseVector& V) const
  {
    const Vector<real_t>& u = dynamic_cast<const Vector<real_t>&>(U);
    Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(V);

    const size_t& numberOfUnknown = __problem.numberOfUnknown();
    
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseValuesU(numberOfUnknown);
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseDerivativeValuesU(numberOfUnknown);
    
    for (size_t unknownNumber =0; unknownNumber < numberOfUnknown;  ++unknownNumber){
      const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
        = __basisU[unknownNumber];
      
      for(size_t i=0; i<3; ++i){ 
        quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        quadratureBaseDerivativeValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        for(size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints();++j){
          real_t xi=(*__transformU[unknownNumber][i]).inverse((*__transform[i])((*__gaussLobatto[i])(j)));
          quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          quadratureBaseDerivativeValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          unknownBasis[i]->getValues(xi,quadratureBaseValuesU[unknownNumber][i][j]);
          unknownBasis[i]->getDerivativeValues(xi,quadratureBaseDerivativeValuesU[unknownNumber][i][j]);
        }
      }
    }  
    
    switch (__problem.type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem = dynamic_cast<const PDESystem&>(__problem);
      
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented yet!",
                         ErrorHandler::unexpected);
      break;
    }
    case Problem::variational: {
      const VariationalProblem& P
        = dynamic_cast<const VariationalProblem&>(__problem);
      
      for (VariationalProblem::bilinearOperatorConst_iterator
             iOperator = P.beginBilinearOperator();
           iOperator != P.endBilinearOperator(); ++iOperator) {
        
        switch ((**iOperator).type()) {  
        case VariationalBilinearOperator::alphaUV: {
          const VariationalAlphaUVOperator& alphaUV
            = dynamic_cast<const VariationalAlphaUVOperator&>(**iOperator);
          
          const ScalarFunctionBase& alpha = *alphaUV.alpha();
          
          Structured3DVector<real_t> alpha_rst(Array3DShape(__gaussLobatto[0]->numberOfPoints(),
                                                            __gaussLobatto[1]->numberOfPoints(),
                                                            __gaussLobatto[2]->numberOfPoints())); 


          _interpolateLagrange(alpha,alpha_rst);
          _getAu(alpha_rst,
                 1.,1.,1.,
                 1.,1.,1.,
                 quadratureBaseValuesU[alphaUV.unknownNumber()][0],
                 quadratureBaseValuesU[alphaUV.unknownNumber()][1],
                 quadratureBaseValuesU[alphaUV.unknownNumber()][2],
                 quadratureBaseValuesU[alphaUV.testFunctionNumber()][0],
                 quadratureBaseValuesU[alphaUV.testFunctionNumber()][1],
                 quadratureBaseValuesU[alphaUV.testFunctionNumber()][2],
                 alphaUV.unknownNumber(),alphaUV.testFunctionNumber(),
                 u,v);
          break;         
        }
          
        case VariationalBilinearOperator::muGradUGradV: {
          const VariationalMuGradUGradVOperator& muGradUgradV
            = dynamic_cast<const VariationalMuGradUGradVOperator&>(**iOperator);
          const ScalarFunctionBase& mu = *muGradUgradV.mu();
          
          Structured3DVector<real_t> mu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());
          _interpolateLagrange(mu,mu_rst);
          
          _getAu(mu_rst,
                 __transformU[muGradUgradV.unknownNumber()][0]->inverseDeterminant(),1.,1.,
                 __transformU[muGradUgradV.testFunctionNumber()][0]->inverseDeterminant(),1.,1.,
                 quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                 quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                 muGradUgradV.unknownNumber(), muGradUgradV.testFunctionNumber(),
                 u,v);
          
          _getAu(mu_rst,
                 1.,__transformU[muGradUgradV.unknownNumber()][1]->inverseDeterminant(),1.,
                 1.,__transformU[muGradUgradV.testFunctionNumber()][1]->inverseDeterminant(),1.,
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                 quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                 quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                 muGradUgradV.unknownNumber(), muGradUgradV.testFunctionNumber(),
                 u,v);
          
          _getAu(mu_rst,
                 1., 1.,__transformU[muGradUgradV.unknownNumber()][2]->inverseDeterminant(),
                 1., 1.,__transformU[muGradUgradV.testFunctionNumber()][2]->inverseDeterminant(),
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                 quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][2],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                 quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][2],
                 muGradUgradV.unknownNumber(), muGradUgradV.testFunctionNumber(),
                 u,v);
          break;
        }
        case VariationalBilinearOperator::alphaDxUDxV: {
          const VariationalAlphaDxUDxVOperator& alphaDxUDxV
            = dynamic_cast<const VariationalAlphaDxUDxVOperator&>(**iOperator);
          const ScalarFunctionBase& alpha = *alphaDxUDxV.alpha();
          
        
        
        
          // sets default determinant values
          
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          // Set default legendre base quadrature values
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[alphaDxUDxV.unknownNumber()][i];
            pQuadratureVValues[i] = &quadratureBaseValuesU[alphaDxUDxV.testFunctionNumber()][i];
          }
          
          detsU[alphaDxUDxV.j()] *=
            __transformU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()]->inverseDeterminant();
          
          detsV[alphaDxUDxV.i()] *=
            __transformU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()]->inverseDeterminant();
          pQuadratureUValues[alphaDxUDxV.j()] =
            &quadratureBaseDerivativeValuesU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()];
          
          pQuadratureVValues[alphaDxUDxV.i()] =
            &quadratureBaseDerivativeValuesU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()];
        
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());
          _interpolateLagrange(alpha,alpha_rst);
          
          _getAu(alpha_rst,
                 detsU[0],detsU[1], detsU[2],
                 detsV[0],detsV[1], detsV[2],
                 *pQuadratureUValues[0],*pQuadratureUValues[1],*pQuadratureUValues[2],
                 *pQuadratureVValues[0],*pQuadratureVValues[1],*pQuadratureVValues[2],
                 alphaDxUDxV.unknownNumber(), alphaDxUDxV.testFunctionNumber(),
                 u,v);
          
          break;
        }
          
        case VariationalBilinearOperator::nuUdxV: {
          const VariationalNuUdxVOperator& nuUdxV
            = dynamic_cast<const VariationalNuUdxVOperator&>(**iOperator);
          
          const ScalarFunctionBase& nu = *nuUdxV.nu();
          
          
          
          // sets default determinant values
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            pQuadratureVValues[i] = &quadratureBaseValuesU[nuUdxV.testFunctionNumber()][i];
          }
          
          detsV[nuUdxV.i()] *=
            __transformU[nuUdxV.testFunctionNumber()][nuUdxV.i()]->inverseDeterminant();
          
          pQuadratureVValues[nuUdxV.i()] =
            &quadratureBaseDerivativeValuesU[nuUdxV.testFunctionNumber()][nuUdxV.i()];
          
          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());         

          _interpolateLagrange(nu,nu_rst);
          
          _getAu(nu_rst,
                 1.,1.,1.,
                 detsV[0], detsV[1], detsV[2],
                 quadratureBaseValuesU[nuUdxV.unknownNumber()][0], 
                 quadratureBaseValuesU[nuUdxV.unknownNumber()][1], 
                 quadratureBaseValuesU[nuUdxV.unknownNumber()][2],
                 *pQuadratureVValues[0],
                 *pQuadratureVValues[1],
                 *pQuadratureVValues[2],
                 nuUdxV.unknownNumber(), nuUdxV.testFunctionNumber(),
                 u,v);
          
          break;
        }
          
          
        case VariationalBilinearOperator::nuDxUV: {
          const VariationalNuDxUVOperator& nuDxUV
            = dynamic_cast<const VariationalNuDxUVOperator&>(**iOperator);
          
          const ScalarFunctionBase& nu = *nuDxUV.nu();
          
          
          
          // sets default determinant values

          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          
          // Set default legendre base quadrature values
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[nuDxUV.unknownNumber()][i];
          }
          
          detsU[nuDxUV.i()] *=
            __transformU[nuDxUV.unknownNumber()][nuDxUV.i()]->inverseDeterminant();
          
          pQuadratureUValues[nuDxUV.i()] = 
            &quadratureBaseDerivativeValuesU[nuDxUV.unknownNumber()][nuDxUV.i()];
          
          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());
         
          _interpolateLagrange(nu,nu_rst);
          
          _getAu(nu_rst,
                 detsU[0], detsU[1], detsU[2],
                 1.,1.,1.,
                 *pQuadratureUValues[0],
                 *pQuadratureUValues[1],
                 *pQuadratureUValues[2],
                 quadratureBaseValuesU[nuDxUV.testFunctionNumber()][0], 
                 quadratureBaseValuesU[nuDxUV.testFunctionNumber()][1], 
                 quadratureBaseValuesU[nuDxUV.testFunctionNumber()][2],
                 nuDxUV.unknownNumber(), nuDxUV.testFunctionNumber(),
                 u,v);
          
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unexpected operator type",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected problem type",
                         ErrorHandler::unexpected);
    }
    }
  }
  
  void transposedTimesX(const BaseVector& X, BaseVector& Z) const
  {
    const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
    
    const size_t& numberOfUnknown = __problem.numberOfUnknown();

    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseValuesU(numberOfUnknown);
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseDerivativeValuesU(numberOfUnknown);
    
    for (size_t unknownNumber =0; unknownNumber < numberOfUnknown; ++ unknownNumber){
      const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
        = __basisU[unknownNumber];
      
      for(size_t i=0; i<3; ++i){ 
        quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        quadratureBaseDerivativeValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        for(size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints();++j){
          real_t xi=(*__transformU[unknownNumber][i]).inverse((*__transform[i])((*__gaussLobatto[i])(j)));
                  
          quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          quadratureBaseDerivativeValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          unknownBasis[i]->getValues(xi,quadratureBaseValuesU[unknownNumber][i][j]);
          unknownBasis[i]->getDerivativeValues(xi,quadratureBaseDerivativeValuesU[unknownNumber][i][j]);
        }
      }
    }
    
    switch (__problem.type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem = dynamic_cast<const PDESystem&>(__problem);
      
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented yet!",
                         ErrorHandler::unexpected);
      break;
    }
    case Problem::variational: {
      const VariationalProblem& P
        = dynamic_cast<const VariationalProblem&>(__problem);
      
      for (VariationalProblem::bilinearOperatorConst_iterator
             iOperator = P.beginBilinearOperator();
           iOperator != P.endBilinearOperator(); ++iOperator) {
        
        switch ((**iOperator).type()) {  
        case VariationalBilinearOperator::alphaUV: {
          const VariationalAlphaUVOperator& alphaUV
            = dynamic_cast<const VariationalAlphaUVOperator&>(**iOperator);
          
          const ScalarFunctionBase& alpha = *alphaUV.alpha();
          
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());
         
          _interpolateLagrange(alpha,alpha_rst);
          
          _getAu(alpha_rst,
                 1.,1.,1.,
                 1.,1.,1.,
                 quadratureBaseValuesU[alphaUV.testFunctionNumber()][0],
                 quadratureBaseValuesU[alphaUV.testFunctionNumber()][1],
                 quadratureBaseValuesU[alphaUV.testFunctionNumber()][2],
                 quadratureBaseValuesU[alphaUV.unknownNumber()][0],
                 quadratureBaseValuesU[alphaUV.unknownNumber()][1],
                 quadratureBaseValuesU[alphaUV.unknownNumber()][2],
                 alphaUV.testFunctionNumber(),alphaUV.unknownNumber(),
                 x,z);
          break;
        }
          
        case VariationalBilinearOperator::muGradUGradV: {
          const VariationalMuGradUGradVOperator& muGradUgradV
            = dynamic_cast<const VariationalMuGradUGradVOperator&>(**iOperator);
          const ScalarFunctionBase& mu = *muGradUgradV.mu();
          
          Structured3DVector<real_t> mu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());        
         
          _interpolateLagrange(mu,mu_rst);
          
          _getAu(mu_rst,
                 __transformU[muGradUgradV.testFunctionNumber()][0]->inverseDeterminant(),1.,1.,
                 __transformU[muGradUgradV.unknownNumber()][0]->inverseDeterminant(),1.,1.,
                 quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                 quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                 muGradUgradV.testFunctionNumber(), muGradUgradV.unknownNumber(),
                 x,z);
                
          _getAu(mu_rst,
                 1.,__transformU[muGradUgradV.testFunctionNumber()][1]->inverseDeterminant(),1.,
                 1.,__transformU[muGradUgradV.unknownNumber()][1]->inverseDeterminant(),1.,
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                 quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                 quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][1],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                 muGradUgradV.testFunctionNumber(), muGradUgradV.unknownNumber(),
                 x,z);
          
          _getAu(mu_rst,
                 1., 1.,__transformU[muGradUgradV.testFunctionNumber()][2]->inverseDeterminant(),
                 1., 1.,__transformU[muGradUgradV.unknownNumber()][2]->inverseDeterminant(),
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                 quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][2],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                 quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                 quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][2],
                 muGradUgradV.testFunctionNumber(), muGradUgradV.unknownNumber(),
                 x,z);
          break;
        }
          
        case VariationalBilinearOperator::alphaDxUDxV: {
          const VariationalAlphaDxUDxVOperator& alphaDxUDxV
            = dynamic_cast<const VariationalAlphaDxUDxVOperator&>(**iOperator);
          const ScalarFunctionBase& alpha = *alphaDxUDxV.alpha();
          
         
         
          // sets default determinant values
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          // Set default legendre base quadrature values
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[alphaDxUDxV.unknownNumber()][i];
            pQuadratureVValues[i] = &quadratureBaseValuesU[alphaDxUDxV.testFunctionNumber()][i];
          }
          
          detsU[alphaDxUDxV.j()] *= __transformU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()]->inverseDeterminant();
          
          detsV[alphaDxUDxV.i()] *= __transformU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()]->inverseDeterminant();
          
          pQuadratureUValues[alphaDxUDxV.j()] = &quadratureBaseDerivativeValuesU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()];
          
          pQuadratureVValues[alphaDxUDxV.i()] = &quadratureBaseDerivativeValuesU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()];
         
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());     
         
          _interpolateLagrange(alpha,alpha_rst);
          
          _getAu(alpha_rst,
                 detsV[0], detsV[1], detsV[2],
                 detsU[0], detsU[1], detsU[2],
                 *pQuadratureVValues[0],*pQuadratureVValues[1],*pQuadratureVValues[2],
                 *pQuadratureUValues[0],*pQuadratureUValues[1],*pQuadratureUValues[2],  
                 alphaDxUDxV.testFunctionNumber(), alphaDxUDxV.unknownNumber(),
                 x,z);
          break;
        }
          
        case VariationalBilinearOperator::nuUdxV: {
          const VariationalNuUdxVOperator& nuUdxV
            = dynamic_cast<const VariationalNuUdxVOperator&>(**iOperator);
          
          const ScalarFunctionBase& nu = *nuUdxV.nu();
                          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            
            pQuadratureVValues[i] = &quadratureBaseValuesU[nuUdxV.testFunctionNumber()][i];
          }
          
          detsV[nuUdxV.i()] *= __transformU[nuUdxV.testFunctionNumber()][nuUdxV.i()]->inverseDeterminant();
          pQuadratureVValues[nuUdxV.i()] = &quadratureBaseDerivativeValuesU[nuUdxV.testFunctionNumber()][nuUdxV.i()];
         
          
          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());
          
          _interpolateLagrange(nu,nu_rst);
           
          _getAu(nu_rst,
                 detsV[0], detsV[1], detsV[2],
                 1.,1.,1.,
                 *pQuadratureVValues[0],
                 *pQuadratureVValues[1], 
                 *pQuadratureVValues[2],
                 quadratureBaseValuesU[ nuUdxV.unknownNumber()][0],
                 quadratureBaseValuesU[ nuUdxV.unknownNumber()][1],
                 quadratureBaseValuesU[ nuUdxV.unknownNumber()][2],
                 nuUdxV.testFunctionNumber(),nuUdxV.unknownNumber(),
                 x,z);
          break;
        }
        case VariationalBilinearOperator::nuDxUV: {
          const VariationalNuDxUVOperator& nuDxUV
            = dynamic_cast<const VariationalNuDxUVOperator&>(**iOperator);
          
          const ScalarFunctionBase& nu = *nuDxUV.nu();
                   
          // sets default determinant values
          
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[nuDxUV.unknownNumber()][i];
          }
          
          detsU[nuDxUV.i()] *= __transformU[nuDxUV.unknownNumber()][nuDxUV.i()]->inverseDeterminant();
          pQuadratureUValues[nuDxUV.i()] = &quadratureBaseDerivativeValuesU[nuDxUV.unknownNumber()][nuDxUV.i()];
          
          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());      
         
          _interpolateLagrange(nu,nu_rst);
          
          _getAu(nu_rst,
                 detsU[0], detsU[1], detsU[2],
                 1.,1.,1.,
                 quadratureBaseValuesU[nuDxUV.testFunctionNumber()][0], 
                 quadratureBaseValuesU[nuDxUV.testFunctionNumber()][1], 
                 quadratureBaseValuesU[nuDxUV.testFunctionNumber()][2],
                 *pQuadratureUValues[0],
                 *pQuadratureUValues[1],
                 *pQuadratureUValues[2],
                 nuDxUV.testFunctionNumber(),nuDxUV.unknownNumber(),
                 x,z);
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unexpected operator type",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected problem type",
                         ErrorHandler::unexpected);
    }
    }
  }

  void getMultiDiagonal(BaseMatrix& ABase) const
  {
    DoubleHashedMatrix& A = dynamic_cast <DoubleHashedMatrix&> (ABase);
 
    const size_t& numberOfUnknown = __problem.numberOfUnknown();
    
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseValuesU(numberOfUnknown);
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseDerivativeValuesU(numberOfUnknown);
    
    for (size_t unknownNumber =0; unknownNumber < numberOfUnknown;  ++unknownNumber){
      const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
        = __basisU[unknownNumber];
    
      for(size_t i=0; i<3; ++i){ 
        quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        quadratureBaseDerivativeValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        for(size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints();++j){
          real_t xi=(*__transformU[unknownNumber][i]).inverse((*__transform[i])((*__gaussLobatto[i])(j)));
          quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          quadratureBaseDerivativeValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          unknownBasis[i]->getValues(xi,quadratureBaseValuesU[unknownNumber][i][j]);
          unknownBasis[i]->getDerivativeValues(xi,quadratureBaseDerivativeValuesU[unknownNumber][i][j]);
        }
      }
    }
    switch (__problem.type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem = dynamic_cast<const PDESystem&>(__problem);
  
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented yet!",
                         ErrorHandler::unexpected);
      break;
    }
    case Problem::variational: {
      const VariationalProblem& P
        = dynamic_cast<const VariationalProblem&>(__problem);
      
      for (VariationalProblem::bilinearOperatorConst_iterator
             iOperator = P.beginBilinearOperator();
           iOperator != P.endBilinearOperator(); ++iOperator) {
        
        switch ((**iOperator).type()) {  
        case VariationalBilinearOperator::alphaUV: {
          const VariationalAlphaUVOperator& alphaUV
            = dynamic_cast<const VariationalAlphaUVOperator&>(**iOperator);
          
          const ScalarFunctionBase& alpha = *alphaUV.alpha();
          
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());      
         
          _interpolateLagrange(alpha,alpha_rst);
          
          _getA(alpha_rst,
                1.,1.,1.,
                1.,1.,1.,
                quadratureBaseValuesU[alphaUV.unknownNumber()][0],
                quadratureBaseValuesU[alphaUV.unknownNumber()][1],
                quadratureBaseValuesU[alphaUV.unknownNumber()][2],
                quadratureBaseValuesU[alphaUV.testFunctionNumber()][0],
                quadratureBaseValuesU[alphaUV.testFunctionNumber()][1],
                quadratureBaseValuesU[alphaUV.testFunctionNumber()][2],
                alphaUV.unknownNumber(),
                A);
          break;
        }
        
        case VariationalBilinearOperator::muGradUGradV: {
          const VariationalMuGradUGradVOperator& muGradUgradV
            = dynamic_cast<const VariationalMuGradUGradVOperator&>(**iOperator);
          const ScalarFunctionBase& mu = *muGradUgradV.mu();
          
          Structured3DVector<real_t> mu_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());
          
          
          _interpolateLagrange(mu,mu_rst);
          
          _getA(mu_rst,
                __transformU[muGradUgradV.unknownNumber()][0]->inverseDeterminant(),1.,1.,
                __transformU[muGradUgradV.testFunctionNumber()][0]->inverseDeterminant(),1.,1.,
                quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][0],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][0],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                muGradUgradV.unknownNumber(),
                A);
          
          _getA(mu_rst,
                1.,__transformU[muGradUgradV.unknownNumber()][1]->inverseDeterminant(),1.,
                1.,__transformU[muGradUgradV.testFunctionNumber()][1]->inverseDeterminant(),1.,
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][1],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][1],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                muGradUgradV.unknownNumber(),
                A);
          
          _getA(mu_rst,
                1., 1.,__transformU[muGradUgradV.unknownNumber()][2]->inverseDeterminant(),
                1., 1.,__transformU[muGradUgradV.testFunctionNumber()][2]->inverseDeterminant(),
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][2],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][2],
                muGradUgradV.unknownNumber(),
                A);
          break;
        }
        case VariationalBilinearOperator::alphaDxUDxV: {
          const VariationalAlphaDxUDxVOperator& alphaDxUDxV
            = dynamic_cast<const VariationalAlphaDxUDxVOperator&>(**iOperator);
          const ScalarFunctionBase& alpha = *alphaDxUDxV.alpha();
          
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          // Set default legendre base quadrature values
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[alphaDxUDxV.unknownNumber()][i];
            pQuadratureVValues[i] = &quadratureBaseValuesU[alphaDxUDxV.unknownNumber()][i];
          }
          
          detsU[alphaDxUDxV.j()] *=
            __transformU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()]->inverseDeterminant();
          
          detsV[alphaDxUDxV.i()] *=
            __transformU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()]->inverseDeterminant();
          pQuadratureUValues[alphaDxUDxV.j()] =
            &quadratureBaseDerivativeValuesU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()];
          
          pQuadratureVValues[alphaDxUDxV.i()] =
            &quadratureBaseDerivativeValuesU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()];
          
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());      
         
          _interpolateLagrange(alpha,alpha_rst);
          
          _getA(alpha_rst,
                detsU[0],detsU[1], detsU[2],
                detsV[0],detsV[1], detsV[2],
                *pQuadratureUValues[0],*pQuadratureUValues[1],*pQuadratureUValues[2],
                *pQuadratureVValues[0],*pQuadratureVValues[1],*pQuadratureVValues[2],
                alphaDxUDxV.unknownNumber(),
                A);
          break;
        }
        case VariationalBilinearOperator::nuUdxV: {
          const VariationalNuUdxVOperator& nuUdxV
            = dynamic_cast<const VariationalNuUdxVOperator&>(**iOperator);
          
          const ScalarFunctionBase& nu = *nuUdxV.nu();
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            pQuadratureVValues[i] = &quadratureBaseValuesU[nuUdxV.unknownNumber()][i];
          }
          
          detsV[nuUdxV.i()] *=
            __transformU[nuUdxV.testFunctionNumber()][nuUdxV.i()]->inverseDeterminant();
          
          pQuadratureVValues[nuUdxV.i()] =
            &quadratureBaseDerivativeValuesU[nuUdxV.testFunctionNumber()][nuUdxV.i()];

          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());     
        
          _interpolateLagrange(nu,nu_rst);
          
          _getA(nu_rst,
                1.,1.,1.,
                detsV[0], detsV[1], detsV[2],
                quadratureBaseValuesU[nuUdxV.unknownNumber()][0], 
                quadratureBaseValuesU[nuUdxV.unknownNumber()][1], 
                quadratureBaseValuesU[nuUdxV.unknownNumber()][2],
                *pQuadratureVValues[0],  *pQuadratureVValues[1],  *pQuadratureVValues[2],
                nuUdxV.unknownNumber(),
                A);
          break;
        }
        case VariationalBilinearOperator::nuDxUV: {
          const VariationalNuDxUVOperator& nuDxUV
            = dynamic_cast<const VariationalNuDxUVOperator&>(**iOperator);

          const ScalarFunctionBase& nu = *nuDxUV.nu();

          // sets default determinant values
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[nuDxUV.unknownNumber()][i];
          }
          detsU[nuDxUV.i()] *=
            __transformU[nuDxUV.unknownNumber()][nuDxUV.i()]->inverseDeterminant();
          
          pQuadratureUValues[nuDxUV.i()] = 
            &quadratureBaseDerivativeValuesU[nuDxUV.unknownNumber()][nuDxUV.i()];
          
          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());        

                  
          _interpolateLagrange(nu,nu_rst);
          
          _getA(nu_rst,
                detsU[0], detsU[1], detsU[2],
                1.,1.,1.,
                *pQuadratureUValues[0],
                *pQuadratureUValues[1],
                *pQuadratureUValues[2],
                quadratureBaseValuesU[nuDxUV.testFunctionNumber()][0], 
                quadratureBaseValuesU[nuDxUV.testFunctionNumber()][1], 
                quadratureBaseValuesU[nuDxUV.testFunctionNumber()][2],
                nuDxUV.unknownNumber(),
                A);
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unexpected operator type",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected problem type",
                         ErrorHandler::unexpected);
    }
    }
  }
 
  void getDiagonal(BaseVector& Z) const
  {
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
 
    const size_t& numberOfUnknown = __problem.numberOfUnknown();
    
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseValuesU(numberOfUnknown);
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseDerivativeValuesU(numberOfUnknown);
    
    for (size_t unknownNumber =0; unknownNumber < numberOfUnknown;  ++unknownNumber){
      const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
        = __basisU[unknownNumber];
    
      for(size_t i=0; i<3; ++i){ 
        quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        quadratureBaseDerivativeValuesU[unknownNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        for(size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints();++j){
          real_t xi=(*__transformU[unknownNumber][i]).inverse((*__transform[i])((*__gaussLobatto[i])(j)));
          quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          quadratureBaseDerivativeValuesU[unknownNumber][i][j].resize(unknownBasis[i]->dimension());
          unknownBasis[i]->getValues(xi,quadratureBaseValuesU[unknownNumber][i][j]);
          unknownBasis[i]->getDerivativeValues(xi,quadratureBaseDerivativeValuesU[unknownNumber][i][j]);
        }
      }
    }  
    switch (__problem.type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem = dynamic_cast<const PDESystem&>(__problem);
  
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented yet!",
                         ErrorHandler::unexpected);
      break;
    }
    case Problem::variational: {
      const VariationalProblem& P
        = dynamic_cast<const VariationalProblem&>(__problem);
      
      for (VariationalProblem::bilinearOperatorConst_iterator
             iOperator = P.beginBilinearOperator();
           iOperator != P.endBilinearOperator(); ++iOperator) {
        
        switch ((**iOperator).type()) {  
        case VariationalBilinearOperator::alphaUV: {
          const VariationalAlphaUVOperator& alphaUV
            = dynamic_cast<const VariationalAlphaUVOperator&>(**iOperator);
          
          const ScalarFunctionBase& alpha = *alphaUV.alpha();
          
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());      
         
          _interpolateLagrange(alpha,alpha_rst);
          
          _getD(alpha_rst,
                1.,1.,1.,
                1.,1.,1.,
                quadratureBaseValuesU[alphaUV.unknownNumber()][0],
                quadratureBaseValuesU[alphaUV.unknownNumber()][1],
                quadratureBaseValuesU[alphaUV.unknownNumber()][2],
                quadratureBaseValuesU[alphaUV.testFunctionNumber()][0],
                quadratureBaseValuesU[alphaUV.testFunctionNumber()][1],
                quadratureBaseValuesU[alphaUV.testFunctionNumber()][2],
                alphaUV.unknownNumber(),
                z);
          break;
        }
        
        case VariationalBilinearOperator::muGradUGradV: {
          const VariationalMuGradUGradVOperator& muGradUgradV
            = dynamic_cast<const VariationalMuGradUGradVOperator&>(**iOperator);
          const ScalarFunctionBase& mu = *muGradUgradV.mu();
          
          Structured3DVector<real_t> mu_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());
          
          
          _interpolateLagrange(mu,mu_rst);
          
          _getD(mu_rst,
                __transformU[muGradUgradV.unknownNumber()][0]->inverseDeterminant(),1.,1.,
                __transformU[muGradUgradV.testFunctionNumber()][0]->inverseDeterminant(),1.,1.,
                quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][0],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][0],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                muGradUgradV.unknownNumber(),
                z);
          
          _getD(mu_rst,
                
                1.,__transformU[muGradUgradV.unknownNumber()][1]->inverseDeterminant(),1.,
                1.,__transformU[muGradUgradV.testFunctionNumber()][1]->inverseDeterminant(),1.,
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][1],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][2],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][1],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][2],
                muGradUgradV.unknownNumber(),
                z);
          
          _getD(mu_rst,
                1., 1.,__transformU[muGradUgradV.unknownNumber()][2]->inverseDeterminant(),
                1., 1.,__transformU[muGradUgradV.testFunctionNumber()][2]->inverseDeterminant(),
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][0],
                quadratureBaseValuesU[muGradUgradV.unknownNumber()][1],
                quadratureBaseDerivativeValuesU[muGradUgradV.unknownNumber()][2],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][0],
                quadratureBaseValuesU[muGradUgradV.testFunctionNumber()][1],
                quadratureBaseDerivativeValuesU[muGradUgradV.testFunctionNumber()][2],
                muGradUgradV.unknownNumber(),
                z);
          break;
        }
        case VariationalBilinearOperator::alphaDxUDxV: {
          const VariationalAlphaDxUDxVOperator& alphaDxUDxV
            = dynamic_cast<const VariationalAlphaDxUDxVOperator&>(**iOperator);
          const ScalarFunctionBase& alpha = *alphaDxUDxV.alpha();
          
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          // Set default legendre base quadrature values
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[alphaDxUDxV.unknownNumber()][i];
            pQuadratureVValues[i] = &quadratureBaseValuesU[alphaDxUDxV.unknownNumber()][i];
          }
          
          detsU[alphaDxUDxV.j()] *=
            __transformU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()]->inverseDeterminant();
          
          detsV[alphaDxUDxV.i()] *=
            __transformU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()]->inverseDeterminant();
          pQuadratureUValues[alphaDxUDxV.j()] =
            &quadratureBaseDerivativeValuesU[alphaDxUDxV.unknownNumber()][alphaDxUDxV.j()];
          
          pQuadratureVValues[alphaDxUDxV.i()] =
            &quadratureBaseDerivativeValuesU[alphaDxUDxV.testFunctionNumber()][alphaDxUDxV.i()];
          
          Structured3DVector<real_t> alpha_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());      
         
          _interpolateLagrange(alpha,alpha_rst);
          
          _getD(alpha_rst,
                detsU[0],detsU[1], detsU[2],
                detsV[0],detsV[1], detsV[2],
                *pQuadratureUValues[0],*pQuadratureUValues[1],*pQuadratureUValues[2],
                *pQuadratureVValues[0],*pQuadratureVValues[1],*pQuadratureVValues[2],
                alphaDxUDxV.unknownNumber(),
                z);
          break;
        }

        case VariationalBilinearOperator::nuUdxV: {
          const VariationalNuUdxVOperator& nuUdxV
            = dynamic_cast<const VariationalNuUdxVOperator&>(**iOperator);
          
          const ScalarFunctionBase& nu = *nuUdxV.nu();
          
          TinyVector<3, real_t> detsV;
          detsV[0] = 1.;
          detsV[1] = 1.;
          detsV[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureVValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            pQuadratureVValues[i] = &quadratureBaseValuesU[nuUdxV.unknownNumber()][i];
          }
          
          detsV[nuUdxV.i()] *=
            __transformU[nuUdxV.testFunctionNumber()][nuUdxV.i()]->inverseDeterminant();
          
          pQuadratureVValues[nuUdxV.i()] =
            &quadratureBaseDerivativeValuesU[nuUdxV.testFunctionNumber()][nuUdxV.i()];

          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                               __gaussLobatto[1]->numberOfPoints(),
                                               __gaussLobatto[2]->numberOfPoints());     
        
          _interpolateLagrange(nu,nu_rst);
          
          _getD(nu_rst,
                1.,1.,1.,
                detsV[0], detsV[1], detsV[2],
                quadratureBaseValuesU[nuUdxV.unknownNumber()][0], 
                quadratureBaseValuesU[nuUdxV.unknownNumber()][1], 
                quadratureBaseValuesU[nuUdxV.unknownNumber()][2],
                *pQuadratureVValues[0],  *pQuadratureVValues[1],  *pQuadratureVValues[2],
                nuUdxV.unknownNumber(),
                z);
          break;
        }
        case VariationalBilinearOperator::nuDxUV: {
          const VariationalNuDxUVOperator& nuDxUV
            = dynamic_cast<const VariationalNuDxUVOperator&>(**iOperator);

          const ScalarFunctionBase& nu = *nuDxUV.nu();

          // sets default determinant values
          TinyVector<3, real_t> detsU;
          detsU[0] = 1.;
          detsU[1] = 1.;
          detsU[2] = 1.;
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureUValues(0,0,0);
          
          for (size_t i=0; i<3;++i) {
            pQuadratureUValues[i] = &quadratureBaseValuesU[nuDxUV.unknownNumber()][i];
          }
          detsU[nuDxUV.i()] *=
            __transformU[nuDxUV.unknownNumber()][nuDxUV.i()]->inverseDeterminant();
          
          pQuadratureUValues[nuDxUV.i()] = 
            &quadratureBaseDerivativeValuesU[nuDxUV.unknownNumber()][nuDxUV.i()];
          
          Structured3DVector<real_t> nu_rst(__gaussLobatto[0]->numberOfPoints(),
                                            __gaussLobatto[1]->numberOfPoints(),
                                            __gaussLobatto[2]->numberOfPoints());        

                  
          _interpolateLagrange(nu,nu_rst);
          
          _getD(nu_rst,
                detsU[0], detsU[1], detsU[2],
                1.,1.,1.,
                *pQuadratureUValues[0],
                *pQuadratureUValues[1],
                *pQuadratureUValues[2],
                quadratureBaseValuesU[nuDxUV.testFunctionNumber()][0], 
                quadratureBaseValuesU[nuDxUV.testFunctionNumber()][1], 
                quadratureBaseValuesU[nuDxUV.testFunctionNumber()][2],
                nuDxUV.unknownNumber(),
                z);
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unexpected operator type",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected problem type",
                         ErrorHandler::unexpected);
    }
    }
  }

  void assembleSecondMember()
  {
    Vector<real_t>& b = (static_cast<Vector<real_t>&>(this->__b));

    //    const ScalarDegreeOfFreedomPositionsSet& dofPositions
    //= this->__degreeOfFreedomSet.positionsSet(0);
    
    const size_t& numberOfUnknown = __problem.numberOfUnknown();
    
    Vector<TinyVector<3,Vector<Vector<real_t> > > > quadratureBaseValuesU(numberOfUnknown); 
    Vector<TinyVector<3,Vector<Vector<real_t> > > >  quadratureBaseDerivativeValuesU(numberOfUnknown); 
    
    for (size_t testFunctionNumber =0; testFunctionNumber < numberOfUnknown; ++ testFunctionNumber){
      
      const  TinyVector<3, ConstReferenceCounting<LegendreBasis> >& testBasis
        = __basisU[testFunctionNumber];    
      for(size_t i=0; i<3; ++i){ 
        quadratureBaseValuesU[testFunctionNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        quadratureBaseDerivativeValuesU[testFunctionNumber][i].resize((*__gaussLobatto[i]).numberOfPoints());
        for(size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints();++j){
          real_t xi= (*__transformU[testFunctionNumber][i]).inverse((*__transform[i])((*__gaussLobatto[i])(j)));
          quadratureBaseValuesU[testFunctionNumber][i][j].resize(testBasis[i]->dimension());
          quadratureBaseDerivativeValuesU[testFunctionNumber][i][j].resize(testBasis[i]->dimension());
          testBasis[i]->getValues(xi,quadratureBaseValuesU[testFunctionNumber][i][j]);
          testBasis[i]->getDerivativeValues(xi,quadratureBaseDerivativeValuesU[testFunctionNumber][i][j]);
        }
      }
    }  
    
    // Transformations of the nodes 
    TinyVector<3,Vector<real_t> > nodes;
    for( size_t i=0; i<3; ++i){
      nodes[i].resize((*__gaussLobatto[i]).numberOfPoints());   
      for (size_t j=0; j<(*__gaussLobatto[i]).numberOfPoints(); ++j) {
        nodes[i][j] =(*__transform[i])((*__gaussLobatto[i])(j)); 
      }
    }
    
    switch(__problem.type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem
        = static_cast<const PDESystem&>(__problem);
      
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented yet!",
                         ErrorHandler::unexpected);
      break;
    }
    case Problem::variational: {
      const VariationalProblem& variationalProblem
        =  dynamic_cast<const VariationalProblem&>(__problem);
      
      for (VariationalProblem::linearOperatorConst_iterator i
             = variationalProblem.beginLinearOperator();
           i != variationalProblem.endLinearOperator(); ++i) {
        switch ((*(*i)).type()) {
        case VariationalLinearOperator::FV: {
          const VariationalOperatorFV& fv
            = dynamic_cast<const VariationalOperatorFV&>(*(*i));
          
          const ScalarFunctionBase& f = fv.f();
          
          for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
            const real_t& x = nodes[0][r];
            const real_t& wr = __gaussLobatto[0]->weight(r);
            
            for (size_t s=0; s < __gaussLobatto[1]->numberOfPoints(); ++s) {
              const real_t& y = nodes[1][s];
              const real_t& wrs = __gaussLobatto[1]->weight(s)*wr;
              
              for (size_t t=0; t < __gaussLobatto[2]->numberOfPoints(); ++t) {    
                const real_t& z = nodes[2][t];
                const real_t& wrst = __gaussLobatto[2]->weight(t)*wrs;
                
                const real_t fweight = f(x,y,z)* wrst *
                  (*__transform[0]).determinant()*  
                  (*__transform[1]).determinant()*
                  (*__transform[2]).determinant();
                
                const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& testBasis
                  = __basisU[fv.testFunctionNumber()];          
                
                const TinyVector<3,size_t>& dofShape = __dofShape[fv.testFunctionNumber()];
                for (size_t m=0; m< testBasis[0]->dimension(); ++m) {
                  const real_t p1 = fweight*quadratureBaseValuesU[fv.testFunctionNumber()][0][r][m];
                  for (size_t p=0; p< testBasis[1]->dimension(); ++p) {
                    const real_t p2 = quadratureBaseValuesU[fv.testFunctionNumber()][1][s][p] *p1;
                    for (size_t q=0; q< testBasis[2]->dimension(); ++q) {
                      const real_t p3 = quadratureBaseValuesU[fv.testFunctionNumber()][2][t][q]*p2;
                      size_t l  = m*dofShape[1]*dofShape[2] + p*dofShape[2] + q;
                      b[__degreeOfFreedomSet(fv.testFunctionNumber(),l)] +=p3;
                    }
                  }
                }
              }
            }
          }
              
          break;
        }
          
        case VariationalLinearOperator::FdxV: {
          const VariationalOperatorFdxV& fdxV
            = dynamic_cast<const VariationalOperatorFdxV&>(*(*i));
          
          const ScalarFunctionBase& f = fdxV.f();
          
          const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& testBasis
            = __basisU[fdxV.testFunctionNumber()];
          
          TinyVector<3,real_t> detsV;
          detsV[0] =1.;
          detsV[1] =1.;
          detsV[2] =1.;
          
          TinyVector<3,Vector<Vector<real_t> >*>pQuadratureVValues(0,0,0);
          for (size_t i=0;i<3;i++){
            pQuadratureVValues[i] = &quadratureBaseValuesU[fdxV.testFunctionNumber()][i];
          }
          
          detsV[fdxV.number()] *=
            __transformU[fdxV.testFunctionNumber()][fdxV.number()]->inverseDeterminant();
          
          pQuadratureVValues[fdxV.number()] = 
            &quadratureBaseDerivativeValuesU[fdxV.testFunctionNumber()][fdxV.number()];
          
          for (size_t r=0; r < __gaussLobatto[0]->numberOfPoints(); ++r) {
            const real_t& x = nodes[0][r];
            const real_t& wr = __gaussLobatto[0]->weight(r);
            
            for (size_t s=0; s < __gaussLobatto[1]->numberOfPoints(); ++s) {
              const real_t& y = nodes[1][s];
              const real_t& wrs = __gaussLobatto[1]->weight(s)*wr;
              
              for (size_t t=0; t < __gaussLobatto[2]->numberOfPoints(); ++t) {    
                const real_t& z = nodes[2][t];
                const real_t& wrst = __gaussLobatto[2]->weight(t)*wrs;
                
                const real_t fweight = f(x,y,z)* wrst * 
                  (*__transform[0]).determinant()*  
                  (*__transform[1]).determinant()*
                  (*__transform[2]).determinant()*
                  detsV[0]*detsV[1]*detsV[2];
                
                const TinyVector<3,size_t>& dofShape = __dofShape[fdxV.testFunctionNumber()];
                
                for (size_t m=0; m< testBasis[0]->dimension(); ++m) {
                  const real_t p1 = fweight*(*pQuadratureVValues[0])[r][m];
                  for (size_t p=0; p< testBasis[1]->dimension(); ++p) {
                    const real_t p2=(*pQuadratureVValues[1])[s][p] *p1;
                    for (size_t q=0; q< testBasis[2]->dimension(); ++q) {
                      const real_t p3= (*pQuadratureVValues[2])[t][q]*p2;
                      size_t l  = m*dofShape[1]*dofShape[2] + p*dofShape[2] + q;                                      b[__degreeOfFreedomSet(fdxV.testFunctionNumber(),l)] +=p3;
                    }
                  }
                }
              }
            }
          }
          
          break;        
        }
        case VariationalLinearOperator::FdxGV: {
          const VariationalOperatorFdxGV& fdxGV
            = dynamic_cast<const VariationalOperatorFdxGV&>(*(*i));
          
          const ScalarFunctionBase& f = fdxGV.f();
          const ScalarFunctionBase& g = fdxGV.g();

          TinyVector<3, real_t> A;
          TinyVector<3, real_t> B;
          for (size_t direction=0; direction<3; ++direction) {
            A[direction] = __intervalU[fdxGV.testFunctionNumber()][direction]->a();
            B[direction] = __intervalU[fdxGV.testFunctionNumber()][direction]->b();
          }

          Structured3DMeshShape s3dM(TinyVector<3,size_t>(__basisU[fdxGV.testFunctionNumber()][0]->dimension()-1,
                                                          __basisU[fdxGV.testFunctionNumber()][1]->dimension()-1,
                                                          __basisU[fdxGV.testFunctionNumber()][2]->dimension()-1),
                                     A, B);
          ReferenceCounting<VerticesCorrespondance> correspondance
            = new VerticesCorrespondance((__basisU[fdxGV.testFunctionNumber()][0]->dimension()+1)*
                                         (__basisU[fdxGV.testFunctionNumber()][1]->dimension()+1)*
                                         (__basisU[fdxGV.testFunctionNumber()][2]->dimension()+1));

          ReferenceCounting<SpectralMesh> uMesh = new SpectralMesh(s3dM, correspondance);
          SpectralFunction gProjection(uMesh, g);
          
          TinyVector<3,real_t> detsV;
          detsV[0] =1.;
          detsV[1] =1.;
          detsV[2] =1.;
          
          detsV[fdxGV.number()] *=
            __transformU[fdxGV.testFunctionNumber()][fdxGV.number()]->inverseDeterminant();
          
          TinyVector<3,Vector<Vector<real_t> >*> pQuadratureGValues(0,0,0);
          for (size_t i=0; i<3; ++i) {
            pQuadratureGValues[i] = &quadratureBaseValuesU[fdxGV.testFunctionNumber()][i];
          }
          
          pQuadratureGValues[fdxGV.number()] = 
            &quadratureBaseDerivativeValuesU[fdxGV.testFunctionNumber()][fdxGV.number()];
          
          Structured3DVector<real_t> f_rst(__gaussLobatto[0]->numberOfPoints(),
                                           __gaussLobatto[1]->numberOfPoints(),
                                           __gaussLobatto[2]->numberOfPoints());                  

          _interpolateLagrange(f,f_rst);
          
          _getAu(f_rst,
                 1.,1.,1.,
                 detsV[0], detsV[1], detsV[2],
                 *pQuadratureGValues[0],
                 *pQuadratureGValues[1],
                 *pQuadratureGValues[2],
                 quadratureBaseValuesU[fdxGV.testFunctionNumber()][0],
                 quadratureBaseValuesU[fdxGV.testFunctionNumber()][1], 
                 quadratureBaseValuesU[fdxGV.testFunctionNumber()][2],
                 fdxGV.testFunctionNumber(),fdxGV.testFunctionNumber(),
                 gProjection.values(),b);
          
          break;
        }
          
        case VariationalLinearOperator::FgradGgradV: {
          const VariationalOperatorFgradGgradV& fgradGgradV
            = dynamic_cast<const VariationalOperatorFgradGgradV&>(*(*i));
          
          const ScalarFunctionBase& f = fgradGgradV.f();
          const ScalarFunctionBase& g = fgradGgradV.g();
          
          
           Structured3DVector<real_t> f_rst(Array3DShape(__gaussLobatto[0]->numberOfPoints(),
                                                         __gaussLobatto[1]->numberOfPoints(),
                                                         __gaussLobatto[2]->numberOfPoints()));




          _interpolateLagrange(f,f_rst);
          
          TinyVector<3, real_t> A;
          TinyVector<3, real_t> B;
          for (size_t direction=0; direction<3; ++direction) {
            A[direction] = __intervalU[fgradGgradV.testFunctionNumber()][direction]->a();
            B[direction] = __intervalU[fgradGgradV.testFunctionNumber()][direction]->b();
          }

          Structured3DMeshShape s3dM(TinyVector<3,size_t>(__basisU[fgradGgradV.testFunctionNumber()][0]->dimension()-1,
                                                          __basisU[fgradGgradV.testFunctionNumber()][1]->dimension()-1,
                                                          __basisU[fgradGgradV.testFunctionNumber()][2]->dimension()-1),
                                     A, B);
          ReferenceCounting<VerticesCorrespondance> correspondance
            = new VerticesCorrespondance((__basisU[fgradGgradV.testFunctionNumber()][0]->dimension()+1)*
                                         (__basisU[fgradGgradV.testFunctionNumber()][1]->dimension()+1)*
                                         (__basisU[fgradGgradV.testFunctionNumber()][2]->dimension()+1));

          ReferenceCounting<SpectralMesh> uMesh = new SpectralMesh(s3dM, correspondance);
          SpectralFunction gProjection(uMesh, g);
          
          _getAu(f_rst,
                 __transformU[fgradGgradV.testFunctionNumber()][0]->inverseDeterminant(),1.,1.,
                 __transformU[fgradGgradV.testFunctionNumber()][0]->inverseDeterminant(),1.,1.,
                 quadratureBaseDerivativeValuesU[fgradGgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][2],
                 quadratureBaseDerivativeValuesU[fgradGgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][2],
                 fgradGgradV.testFunctionNumber(), fgradGgradV.testFunctionNumber(),
                 gProjection.values(),b);
          
          _getAu(f_rst,
                 1., __transformU[fgradGgradV.testFunctionNumber()][1]->inverseDeterminant(),1.,
                 1., __transformU[fgradGgradV.testFunctionNumber()][1]->inverseDeterminant(),1.,
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][0],
                 quadratureBaseDerivativeValuesU[fgradGgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][2],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][0],
                 quadratureBaseDerivativeValuesU[fgradGgradV.testFunctionNumber()][1],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][2],
                 fgradGgradV.testFunctionNumber(), fgradGgradV.testFunctionNumber(),
                 gProjection.values(),b);
          
          _getAu(f_rst,
                 1.,1.,__transformU[fgradGgradV.testFunctionNumber()][2]->inverseDeterminant(),
                 1.,1.,__transformU[fgradGgradV.testFunctionNumber()][2]->inverseDeterminant(),
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][1],
                 quadratureBaseDerivativeValuesU[fgradGgradV.testFunctionNumber()][2],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][0],
                 quadratureBaseValuesU[fgradGgradV.testFunctionNumber()][1],
                 quadratureBaseDerivativeValuesU[fgradGgradV.testFunctionNumber()][2],
                 fgradGgradV.testFunctionNumber(), fgradGgradV.testFunctionNumber(),
                 gProjection.values(),b);
          
          break;
        }
          
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unknown variational operator",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unknown problem type",
                         ErrorHandler::unexpected);
    }
    }
  }
  
public:
  
  SpectralLegendreDiscretizer(const Problem& p,
                              const SpectralMesh& m,
                              BaseMatrix& a,
                              BaseVector& bb,
                              const DegreeOfFreedomSet& dof,
                              const DiscretizationType& discretizationType)
    : __problem(p),
      __b(bb),
      __degreeOfFreedomSet(dof),
      __mesh(m),
      __dofShape(__problem.numberOfUnknown()),
      __interval(new Interval(__mesh.shape().a()[0],__mesh.shape().b()[0]),
                 new Interval(__mesh.shape().a()[1],__mesh.shape().b()[1]),
                 new Interval(__mesh.shape().a()[2],__mesh.shape().b()[2])),
      
      __intervalU(__problem.numberOfUnknown()),
      
      __transform(new SpectralConformTransformation(*__interval[0]),
                  new SpectralConformTransformation(*__interval[1]),
                  new SpectralConformTransformation(*__interval[2])),
      
      __transformU(__problem.numberOfUnknown()),
      __gaussLobatto(GaussLobattoManager::instance().getReference(__mesh.degree(0)+1),
                     GaussLobattoManager::instance().getReference(__mesh.degree(1)+1),
                     GaussLobattoManager::instance().getReference(__mesh.degree(2)+1)),
      __basisU(__problem.numberOfUnknown())
  {
    for (size_t i=0; i < __problem.numberOfUnknown(); i++) {
      if (discretizationType[i].type() != ScalarDiscretizationTypeBase::spectralLegendre) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "discretization '"
                           +ScalarDiscretizationTypeBase::name(discretizationType[i])
                           +"' is incompatible with spectral method",
                           ErrorHandler::unexpected);
      }
      const ScalarDiscretizationTypeSpectral& discretization
        = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(discretizationType[i]);

      __dofShape[i] = discretization.degrees()+TinyVector<3,size_t>(1,1,1);

      __intervalU[i][0]= new Interval(discretization.a()[0],discretization.b()[0]);
      __intervalU[i][1]= new Interval(discretization.a()[1],discretization.b()[1]);
      __intervalU[i][2]= new Interval(discretization.a()[2],discretization.b()[2]);
      __transformU[i][0] = new SpectralConformTransformation(*__intervalU[i][0]);
      __transformU[i][1] = new SpectralConformTransformation(*__intervalU[i][1]);
      __transformU[i][2] = new SpectralConformTransformation(*__intervalU[i][2]); 
      __basisU[i][0]= new LegendreBasis(discretization.degrees()[0]);
      __basisU[i][1]= new LegendreBasis(discretization.degrees()[1]);
      __basisU[i][2]= new LegendreBasis(discretization.degrees()[2]);
    }
    
  }

  ~SpectralLegendreDiscretizer()
  {
    ;
  }
};

#endif // SPECTRAL_LEGENDRE_DISCRETIZER_HPP

---