gsExprEvaluator.h

Go to the documentation of this file.

 
#pragma once
 
// #include<gsIO/gsParaviewCollection.h>
#include <fstream>
#include <gsAssembler/gsQuadrature.h>
#include <gsAssembler/gsRemapInterface.h>
#include <gsAssembler/gsCPPInterface.h>
//#include <gsIO/gsWriteParaview.h>
 
namespace gismo
{
 
template<class T>
class gsExprEvaluator
{
private:
    typename gsExprHelper<T>::Ptr m_exprdata;
 
private:
    std::vector<T> m_elWise;
    T              m_value;
 
    gsOptionList m_options;
 
public:
    typedef typename gsBoundaryConditions<T>::bcRefList   bcRefList;
 
    typedef std::vector< boundaryInterface > intContainer;
    typedef std::vector< patchSide > bContainer;
 
    typedef typename gsExprHelper<T>::element     element;
    typedef typename gsExprHelper<T>::geometryMap geometryMap;
    typedef typename gsExprHelper<T>::variable    variable;
 
    typedef typename gsFunction<T>::uPtr ifacemap;
 
public:
 
    gsExprEvaluator() : m_exprdata(gsExprHelper<T>::make()),
    m_options(defaultOptions()) { }
 
    gsExprEvaluator(typename gsExprHelper<T>::Ptr env)
    : m_exprdata(env), m_value(-1), m_options(defaultOptions())
    { }
 
    //gsExprEvaluator(typename gsExprHelper<T> env)
 
    gsExprEvaluator(const gsExprAssembler<T> & o)
    : m_exprdata(o.exprData()), m_options(defaultOptions())
    { }
 
    gsOptionList defaultOptions()
    {
        gsOptionList opt;
        opt.addReal("quA", "Number of quadrature points: quA*deg + quB", 1.0  );
        opt.addInt ("quB", "Number of quadrature points: quA*deg + quB", 1    );
        opt.addInt ("plot.npts", "Number of sampling points for plotting", 3000 );
        opt.addSwitch("plot.elements", "Include the element mesh in plot (when applicable)", false);
        opt.addSwitch("flipSide", "Flip side of interface where evaluation is performed.", false);
        //opt.addSwitch("plot.cnet", "Include the control net in plot (when applicable)", false);
        return opt;
    }
 
    gsOptionList & options() {return m_options;}
 
public:
 
    T value() const { return m_value; }
 
    size_t nValues() const { return m_elWise.size(); }
 
    const std::vector<T> & elementwise() const { return m_elWise; }
 
    gsAsConstVector<T> allValues() const { return gsAsConstVector<T>(m_elWise); }
 
    gsAsConstMatrix<T> allValues(index_t nR, index_t nC) const
    { return gsAsConstMatrix<T>(m_elWise, nR, nC); }
 
    void setIntegrationElements(const gsMultiBasis<T> & mesh)
    { m_exprdata->setMultiBasis(mesh); }
 
    const typename gsExprHelper<T>::Ptr exprData() const { return m_exprdata; }
 
    geometryMap getMap(const gsMultiPatch<T> & mp) //conv->tmp->error
    { return m_exprdata->getMap(mp); }
 
    geometryMap getMap(const gsFunctionSet<T> & gm)
    { return m_exprdata->getMap(gm); }
 
    variable getVariable(const gsFunctionSet<T> & func, index_t dim = 1)
    { return m_exprdata->getVar(func, dim); }
 
    expr::gsComposition<T> getVariable(const gsFunctionSet<T> & func, geometryMap G)
    { return m_exprdata->getVar(func, G); }
 
    element getElement() { return m_exprdata->getElement(); }
 
    void calcSqrt()
    {
        gsAsVector<T> en(m_elWise);
        en.array() = en.array().sqrt();
        m_value = math::sqrt(m_value);
    }
 
    void calcRoot(const index_t p)
    {
        gsAsVector<T> en(m_elWise);
        en.array() = en.array().pow(static_cast<T>(1)/p);
        m_value = math::pow(m_value, static_cast<T>(1)/p );
    }
 
public:
 
    template<class E>
    T integral(const expr::_expr<E> & expr)
    { return compute_impl<E,false,plus_op>(expr); }
 
    template<class E>
    T integralElWise(const expr::_expr<E> & expr)
    { return compute_impl<E,true,plus_op>(expr); }
 
    // overloads for scalars
    T integral(const T & val) { return integral<T>(val); }
    T integralElWise(const T & val) { return integralElWise<T>(val); }
 
    template<class E> // note: integralBdrElWise not offered
    T integralBdr(const expr::_expr<E> & expr)
    { return computeBdr_impl<E,plus_op>(expr,
      m_exprdata->multiBasis().topology().boundaries()); }
 
    template<class E> // note: integralBdrElWise not offered
    T integralBdr(const expr::_expr<E> & expr, const bContainer & bdrlist)
    { return computeBdr_impl<E,plus_op>(expr,bdrlist); }
 
    template<class E> // note: integralBdrElWise not offered
    T integralBdrBc(const bcRefList & BCs, const expr::_expr<E> & expr)
    { return computeBdrBc_impl<E,plus_op>(BCs,expr); }
 
    template<class E> // note: elementwise integral not offered
    T integralInterface(const expr::_expr<E> & expr)
    { return computeInterface_impl<E,plus_op>(expr,
      m_exprdata->multiBasis().topology().interfaces()); }
 
    template<class E> // note: elementwise integral not offered
    T integralInterface(const expr::_expr<E> & expr, const intContainer & iFaces)
    { return computeInterface_impl<E,plus_op>(expr, iFaces); }
 
    template<class E>
    T max(const expr::_expr<E> & expr)
    { return compute_impl<E,false,max_op>(expr); }
 
    template<class E>
    T maxElWise(const expr::_expr<E> & expr)
    { return compute_impl<E,true,max_op>(expr); }
 
    template<class E> // note: elementwise integral not offered
    T maxInterface(const expr::_expr<E> & expr)
    { return computeInterface_impl<E,max_op>(expr, m_exprdata->multiBasis().topology().interfaces()); }
 
    template<class E> // note: elementwise integral not offered
    T maxInterface(const expr::_expr<E> & expr, const intContainer & iFaces)
    { return computeInterface_impl<E,max_op>(expr, iFaces); }
 
    template<class E> // note: elementwise integral not offered
    T minInterface(const expr::_expr<E> & expr)
    { return computeInterface_impl<E,min_op>(expr, m_exprdata->multiBasis().topology().interfaces()); }
 
    template<class E> // note: elementwise integral not offered
    T minInterface(const expr::_expr<E> & expr, const intContainer & iFaces)
    { return computeInterface_impl<E,min_op>(expr, iFaces); }
 
    template<class E>
    T min(const expr::_expr<E> & expr)
    { return compute_impl<E,false,min_op>(expr); }
 
    template<class E>
    T minElWise(const expr::_expr<E> & expr)
    { return compute_impl<E,true,min_op>(expr); }
 
#ifdef __DOXYGEN__
    template<class E> void
#else
    template<class E, int mode, short_t d>
    typename util::enable_if<E::ScalarValued,void>::type
#endif
    eval(const expr::_expr<E> & expr,
         gsGridIterator<T,mode,d> & git,
         const index_t patchInd = 0);
 
    template<class E, int mode, short_t d>
    typename util::enable_if<!E::ScalarValued,void>::type
    eval(const expr::_expr<E> & expr,
              gsGridIterator<T,mode,d> & git,
              const index_t patchInd = 0);
 
    template<class E>
#ifdef __DOXYGEN__
    gsAsConstMatrix<T>
#else
    typename util::enable_if<E::ScalarValued,gsAsConstMatrix<T> >::type
#endif
    eval(const expr::_expr<E> & testExpr, const gsVector<T> & pt,
         const index_t patchInd = 0);
 
    template<class E>
    typename util::enable_if<!E::ScalarValued,gsAsConstMatrix<T> >::type
    eval(const expr::_expr<E> & testExpr, const gsVector<T> & pt,
         const index_t patchInd = 0);
 
    template<class E>
#ifdef __DOXYGEN__
    gsAsConstMatrix<T>
#else
    typename util::enable_if<E::ScalarValued,gsAsConstMatrix<T> >::type
#endif
    evalIfc(const expr::_expr<E> & testExpr, const gsVector<T> & pt,
            const boundaryInterface & ifc);
 
    template<class E>
    typename util::enable_if<!E::ScalarValued,gsAsConstMatrix<T> >::type
    evalIfc(const expr::_expr<E> & testExpr, const gsVector<T> & pt,
            const boundaryInterface & ifc);
 
    template<class E>
#ifdef __DOXYGEN__
    gsAsConstMatrix<T>
#else
    typename util::enable_if<E::ScalarValued,gsAsConstMatrix<T> >::type
#endif
    evalBdr(const expr::_expr<E> & testExpr, const gsVector<T> & pt,
            const patchSide & ps);
 
    template<class E>
    typename util::enable_if<!E::ScalarValued,gsAsConstMatrix<T> >::type
    evalBdr(const expr::_expr<E> & testExpr, const gsVector<T> & pt,
            const patchSide & ps);
 
    template<class E> void
    testEval(const expr::_expr<E> & expr,
             const gsVector<T> & pt, const index_t patchInd = 0)
    {
        expr.print(gsInfo);
        gsInfo << "Result:\n"<< eval(expr,pt,patchInd) <<"\n";
    }
 
    void info() const { m_exprdata->print(); }
 
    // Interpolates the expression \a expr over the isogeometric domain \a G
    template<class E> void interpolate(const expr::_expr<E> &)
    {
        GISMO_NO_IMPLEMENTATION
        // for all patches
        //   get anchors of patch
        //   evaluate expr
        //   solve system
        //   add to multipatch
        // return multipatch
    }
 
    //( expression \a expr over the isogeometric domain \a G.
    // template<class E>
    // void writeParaview(const expr::_expr<E> & expr, const gsMatrix<T> & uv,
    //                    geometryMap G, std::string const & fn)
    // { writeParaview_impl<E,true>(expr,uv,G,fn); }
 
    template<class E>
    void writeParaview(const expr::_expr<E> & expr,
                       geometryMap G, std::string const & fn)
    { writeParaview_impl<E,true>(expr,G,fn); }
 
    //( expression \a expr over the parametric domain.
    template<class E>
    void writeParaview(const expr::_expr<E> & expr,
                       std::string const & fn)
    { writeParaview_impl<E,false>(expr,m_exprdata->getMap(),fn); }
 
private:
 
    template<class E, bool gmap>
    void writeParaview_impl(const expr::_expr<E> & expr,
                            geometryMap G, std::string const & fn);
 
 
 
    template<class E, bool storeElWise, class _op>
    T compute_impl(const expr::_expr<E> & expr);
 
    template<class E, class _op>
    T computeBdr_impl(const expr::_expr<E> & expr, const bContainer & bdrlist);
 
    template<class E, class _op>
    T computeBdrBc_impl(const bcRefList & BCs, const expr::_expr<E> & expr);
 
    template<class E, class _op>
    T computeInterface_impl(const expr::_expr<E> & expr, const intContainer & iFaces);
 
    template<class E>
    void computeGrid_impl(const expr::_expr<E> & expr, const index_t patchInd);
 
    struct plus_op
    {
        static inline T init() { return 0; }
        static inline void acc(const T contrib, const T w, T & res)
        { res += w * contrib; }
 
        static inline void acc_global(const T contrib, T & res)
        {
#           pragma omp atomic update
            res += contrib;
        }
    };
    struct min_op
    {
        static inline T init() { return math::limits::max(); }
        static inline void acc (const T contrib, const T, T & res)
        {res = math::min(contrib, res); }
        static inline void acc_global(const T contrib, T & res)
        {
#           pragma omp atomic write
            res = math::min(contrib, res);
        }
 
    };
    struct max_op
    {
        static inline T init() { return math::limits::min(); }
        static inline void acc (const T contrib, const T, T & res)
        { res = math::max(contrib, res); }
        static inline void acc_global(const T contrib, T & res)
        {
#           pragma omp atomic write
            res = math::max(contrib, res);
        }
    };
 
};
 
template<class T>
template<class E, bool storeElWise, class _op>
T gsExprEvaluator<T>::compute_impl(const expr::_expr<E> & expr)
{
    m_value = _op::init();
    m_elWise.clear();
    if ( storeElWise )
        m_elWise.resize(m_exprdata->multiBasis().totalElements());
 
#pragma omp parallel
{
#ifdef _OPENMP
        const int tid = omp_get_thread_num();
        const int nt  = omp_get_num_threads();
        T thValue = _op::init();
#endif
        gsQuadRule<T> QuRule;  // Quadrature rule
        auto _arg = expr.val();
        m_exprdata->parse(_arg);
        m_exprdata->activateFlags(SAME_ELEMENT);
 
        // Computed value on element
        T elVal;
        index_t poffset = 0;
        typename gsBasis<T>::domainIter domIt;
        for (unsigned patchInd=0; patchInd < m_exprdata->multiBasis().nBases(); ++patchInd)
        {
            // Quadrature rule
            QuRule =  gsQuadrature::get(m_exprdata->multiBasis().basis(patchInd), m_options);
 
            // Initialize domain element iterator
            domIt = m_exprdata->multiBasis().piece(patchInd).makeDomainIterator();
#ifdef _OPENMP
            for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
#else
                for (; domIt->good(); domIt->next() )
#endif
                {
                    // Map the Quadrature rule to the element
                    QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(),
                                  m_exprdata->points(), m_exprdata->weights());
 
                    // Perform required pre-computations on the quadrature nodes
                    m_exprdata->precompute(patchInd);
 
                    // Compute on element
                    elVal = _op::init();
                    for (index_t k = 0; k != m_exprdata->weights().rows(); ++k) // loop over quad. nodes
                        _op::acc(_arg.eval(k), m_exprdata->weights()[k], elVal);
                    _op::acc(elVal, (T)1,
#ifdef _OPENMP
                    thValue);
#else
                    m_value);
#endif
                    if ( storeElWise )
                    {
                        m_elWise[poffset+domIt->id()] = elVal;
                    }
               }
               poffset += m_exprdata->multiBasis().basis(patchInd).numElements();
        }
#ifdef _OPENMP
    _op::acc_global(thValue, m_value);
#endif
}//omp parallel
    return m_value;
}
 
template<class T>
template<class E, class _op>
T gsExprEvaluator<T>::computeBdr_impl(const expr::_expr<E> & expr,
                                      const bContainer & bdrlist)
{
    // GISMO_ASSERT( expr.isScalar(),
    //               "Expecting scalar expression instead of "
    //               <<expr.cols()<<" x "<<expr.rows() );
 
    //expr.print(gsInfo);
 
    gsQuadRule<T> QuRule;  // Quadrature rule
    auto _arg = expr.val();
    m_exprdata->parse(_arg);
    m_exprdata->activateFlags(SAME_ELEMENT);
 
    // Computed value
    T elVal;
    m_value = _op::init();
    m_elWise.clear();
    for (typename gsBoxTopology::const_biterator bit = 
             bdrlist.begin(); bit != bdrlist.end(); ++bit)
    {
        // Quadrature rule
        QuRule = gsQuadrature::get(m_exprdata->multiBasis().basis(bit->patch), m_options,bit->direction());
 
        // Initialize domain element iterator
        typename gsBasis<T>::domainIter domIt =
            m_exprdata->multiBasis().piece(bit->patch).makeDomainIterator(bit->side());
 
        // Start iteration over elements
        for (; domIt->good(); domIt->next() )
        {
            // Map the Quadrature rule to the element
            QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(),
                          m_exprdata->points(), m_exprdata->weights());
 
            // Perform required pre-computations on the quadrature nodes
            m_exprdata->precompute(bit->patch, bit->side() );
 
            // Compute on element
            elVal = _op::init();
            for (index_t k = 0; k != m_exprdata->weights().rows(); ++k) // loop over quadrature nodes
                _op::acc(_arg.eval(k), m_exprdata->weights()[k], elVal);
 
            _op::acc(elVal, 1, m_value);
            //if ( storeElWise ) m_elWise.push_back( elVal );
        }
    }
 
    return m_value;
}
 
template<class T>
template<class E, class _op>
T gsExprEvaluator<T>::computeBdrBc_impl(const bcRefList & BCs,
                                      const expr::_expr<E> & expr)
{
    // GISMO_ASSERT( expr.isScalar(),
    //               "Expecting scalar expression instead of "
    //               <<expr.cols()<<" x "<<expr.rows() );
 
    //expr.print(gsInfo);
 
    if ( BCs.empty() ) return 0;
    m_exprdata->setMutSource(*BCs.front().get().function()); //initialize once
 
    typename gsQuadRule<T>::uPtr QuRule; // Quadrature rule  ---->OUT
    auto _arg = expr.val();
    m_exprdata->parse(_arg);
    m_exprdata->activateFlags(SAME_ELEMENT);
 
    // Computed value
    T elVal;
    m_value = _op::init();
    m_elWise.clear();
 
    for (typename bcRefList::const_iterator iit = BCs.begin(); iit!= BCs.end(); ++iit)
    {
        const boundary_condition<T> * it = &iit->get();
 
        // Quadrature rule
        QuRule = gsQuadrature::getPtr(m_exprdata->multiBasis().basis(it->patch()), m_options, it->side().direction());
 
        // Update boundary function source
        m_exprdata->setMutSource(*it->function());
 
        // Initialize domain element iterator
        typename gsBasis<T>::domainIter domIt =
            m_exprdata->multiBasis().basis(it->patch()).makeDomainIterator(it->side());
 
        // Start iteration over elements
        for (; domIt->good(); domIt->next() )
        {
            // Map the Quadrature rule to the element
            QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(),
                          m_exprdata->points(), m_exprdata->weights());
 
            if (m_exprdata->points().cols()==0)
                continue;
 
            // Perform required pre-computations on the quadrature nodes
            m_exprdata->precompute(it->patch(), it->side() );
 
            // Compute on element
            elVal = _op::init();
            for (index_t k = 0; k != m_exprdata->weights().rows(); ++k) // loop over quadrature nodes
                _op::acc(_arg.eval(k), m_exprdata->weights()[k], elVal);
 
            _op::acc(elVal, 1, m_value);
            //if ( storeElWise ) m_elWise.push_back( elVal );
        }
    }
 
    return m_value;
}
 
template<class T>
template<class E, class _op>
T gsExprEvaluator<T>::computeInterface_impl(const expr::_expr<E> & expr, const intContainer & iFaces)
{
    auto arg_tpl = expr.val();
    m_exprdata->parse(arg_tpl);
    // m_exprdata->activateFlags(SAME_ELEMENT);
 
    typename gsQuadRule<T>::uPtr QuRule;
    // Computed value
    T elVal;
    m_value = _op::init();
    //if ( storeElWise )
    m_elWise.reserve(m_exprdata->multiBasis().topology().nInterfaces());
    m_elWise.clear();
 
    ifacemap interfaceMap;
    for (typename gsBoxTopology::const_iiterator iit =
             iFaces.begin(); iit != iFaces.end(); ++iit)
    {
        const boundaryInterface & iFace = *iit;
        const index_t patch1 = iFace.first().patch;
        const index_t patch2 = iFace.second().patch;
 
        if (iFace.type() == interaction::conforming)
            interfaceMap = gsAffineFunction<T>::make( iFace.dirMap(), iFace.dirOrientation(),
                                                      m_exprdata->multiBasis().basis(patch1).support(),
                                                      m_exprdata->multiBasis().basis(patch2).support() );
        else
            interfaceMap = gsCPPInterface<T>::make(m_exprdata->multiPatch(), m_exprdata->multiBasis(), iFace);
 
        //gsRemapInterface<T> interfaceMap(m_exprdata->multiPatch(),
        //                                 m_exprdata->multiBasis(),
        //                                 *iit);//,opt
 
        // Quadrature rule
        QuRule = gsQuadrature::getPtr(m_exprdata->multiBasis().basis(patch1),
                                      m_options, iFace.first().side().direction());
 
        // Initialize domain element iterator
        typename gsBasis<T>::domainIter domIt =
            //interfaceMap.makeDomainIterator();
            m_exprdata->multiBasis().piece(patch1).makeDomainIterator(iFace.first().side());
 
        // Start iteration over elements
        elVal = _op::init();
        for (; domIt->good(); domIt->next() )
        {
            // Map the Quadrature rule to the element
            QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(),
                           m_exprdata->points(), m_exprdata->weights());
            interfaceMap->eval_into(m_exprdata->points(), m_exprdata->pointsIfc());
 
            // Perform required pre-computations on the quadrature nodes
            m_exprdata->precompute(iFace);
 
            // Compute on element
            for (index_t k = 0; k != m_exprdata->weights().rows(); ++k) // loop over qu-nodes
            {
                _op::acc(arg_tpl.eval(k), m_exprdata->weights()[k], elVal);
            }
        }
        _op::acc(elVal, 1, m_value);
        //if ( storeElWise )
            m_elWise.push_back( elVal );
    }
 
    return m_value;
}
 
 
template<class T>
template<class E, int mode, short_t d>
typename util::enable_if<E::ScalarValued,void>::type
gsExprEvaluator<T>::eval(const expr::_expr<E> & expr,
                         gsGridIterator<T,mode,d> & git,
                         const index_t patchInd)
{ // to remove
    // bug: fails due to gsFeVariable::rows() before evaluation
    // GISMO_ASSERT( expr.isScalar(), "Expecting scalar");
 
    auto _arg = expr.val();
    m_exprdata->parse(_arg);
    m_elWise.clear();
    m_elWise.reserve(git.numPoints());
 
    for( git.reset(); git; ++git )
    {
        m_exprdata->points() = *git;
        m_exprdata->precompute(patchInd);
        m_elWise.push_back( _arg.eval(0) );
 
        // equivalent:
        //m_elWise.push_back( m_exprdata->eval(expr).value() );
 
    }
    m_value = m_elWise.back(); // not used
}
 
template<class T>
template<class E, int mode, short_t d>
typename util::enable_if<!E::ScalarValued,void>::type
gsExprEvaluator<T>::eval(const expr::_expr<E> & expr,
                         gsGridIterator<T,mode,d> & git,
                         const index_t patchInd)
{
    // bug: fails due to gsFeVariable::rows() before evaluation
    // GISMO_ASSERT( expr.isScalar(), "Expecting scalar");
 
    m_exprdata->parse(expr);
    m_elWise.clear();
    m_elWise.reserve(git.numPoints());
    gsMatrix<T> tmp;
    for( git.reset(); git; ++git )
    {
        m_exprdata->points() = *git; // ..
        m_exprdata->precompute(patchInd);
        // m_exprdata->eval_into(expr, tmp);
        tmp = expr.eval(0);
        m_elWise.insert(m_elWise.end(), tmp.data(), tmp.data()+tmp.size());
    }
    m_value = 0; // not used
}
 
 
template<class T>
template<class E>
typename util::enable_if<E::ScalarValued,gsAsConstMatrix<T> >::type
gsExprEvaluator<T>::eval(const expr::_expr<E> & expr, const gsVector<T> & pt,
                         const index_t patchInd)
{
    auto _arg = expr.val();
    m_exprdata->parse(_arg);
    m_elWise.clear();
    m_exprdata->points() = pt;
    m_exprdata->precompute(patchInd);
 
    // expr.printDetail(gsInfo); //
 
    m_value = _arg.eval(0);
    return gsAsConstMatrix<T>(&m_value,1,1);
}
 
template<class T>
template<class E>
typename util::enable_if<!E::ScalarValued,gsAsConstMatrix<T> >::type
gsExprEvaluator<T>::eval(const expr::_expr<E> & expr, const gsVector<T> & pt,
                         const index_t patchInd)
{
    m_exprdata->parse(expr);
    m_exprdata->points() = pt;
    m_exprdata->precompute(patchInd);
 
    // expr.printDetail(gsInfo); //after precompute
 
    gsMatrix<T> tmp = expr.eval(0);
    // const index_t r = expr.rows();
    // const index_t c = expr.cols();
    // gsInfo <<"tmp - "<< tmp.dim() <<" rc - "<< r <<", "<<c<<"\n";
    const index_t r = tmp.rows();
    const index_t c = tmp.cols();
    m_elWise.resize(r*c);
    gsAsMatrix<T>(m_elWise, r, c) = tmp; //expr.eval(0);
    return gsAsConstMatrix<T>(m_elWise, r, c);
}
 
// template<class T>
// template<class E, bool gmap>
// void gsExprEvaluator<T>::writeParaview_impl(const expr::_expr<E> & expr,
//                                             geometryMap G,
//                                             std::string const & fn)
//     {
//         m_exprdata->parse(expr);
 
//         //if false, embed topology ?
//         const index_t n = m_exprdata->multiBasis().nBases();
//         gsParaviewCollection collection(fn);
//         std::string fileName;
 
//         gsMatrix<T> pts, vals, ab;
 
//         const bool mesh = m_options.askSwitch("plot.elements");
 
//         for ( index_t i=0; i != n; ++i )
//         {
//             fileName = fn + util::to_string(i);
//             unsigned nPts = m_options.askInt("plot.npts", 1000);
//             ab = m_exprdata->multiBasis().piece(i).support();
//             gsGridIterator<T,CUBE> pt(ab, nPts);
//             eval(expr, pt, i);
//             nPts = pt.numPoints();
//             vals = allValues(m_elWise.size()/nPts, nPts);
 
//             if (gmap) // Forward the points ?
//             {
//                 eval(G, pt, i);
//                 pts = allValues(m_elWise.size()/nPts, nPts);
//             }
 
//             gsWriteParaviewTPgrid( gmap ? pts : pt.toMatrix(), // parameters
//                                   vals,
//                                   pt.numPointsCwise(), fileName );
//             collection.addPart(fileName + ".vts");
 
//             if ( mesh )
//             {
//                 fileName+= "_mesh";
//                 gsMesh<T> msh(m_exprdata->multiBasis().basis(i), 2);
//                 static_cast<const gsGeometry<T>&>(G.source().piece(i)).evaluateMesh(msh);
//                 gsWriteParaview(msh, fileName, false);
//                 collection.addPart(fileName+ ".vtp");
//             }
//         }
//         collection.save();
//     }
 
 
 
// This is a copy of the above (commented out ) function, which has parts
// of gsParaviewCollection pasted in it. In this way the inclusion of
// gsParaviewCollection.h is prevented. This is a temporary modification.
template<class T>
template<class E>
typename util::enable_if<E::ScalarValued,gsAsConstMatrix<T> >::type
gsExprEvaluator<T>::evalIfc(const expr::_expr<E> & expr, const gsVector<T> & pt,
                            const boundaryInterface & ifc)
{
    auto _arg = expr.val();
    m_exprdata->parse(_arg);
    m_elWise.clear();
 
    const bool flipSide = m_options.askSwitch("flipSide", false);
    const boundaryInterface & iFace =  flipSide ? ifc.getInverse() : ifc;
 
    const index_t patch1 = iFace.first().patch;
    const index_t patch2 = iFace.second().patch;
 
    ifacemap interfaceMap;
    if (iFace.type() == interaction::conforming)
        interfaceMap = gsAffineFunction<T>::make( iFace.dirMap(), iFace.dirOrientation(),
                                                  m_exprdata->multiBasis().basis(patch1).support(),
                                                  m_exprdata->multiBasis().basis(patch2).support() );
    else
        interfaceMap = gsCPPInterface<T>::make(m_exprdata->multiPatch(), m_exprdata->multiBasis(), iFace);
 
    m_exprdata->points() = pt;
    interfaceMap->eval_into(m_exprdata->points(), m_exprdata->pointsIfc());
    m_exprdata->precompute(iFace);
 
    // expr.printDetail(gsInfo); //
 
    m_value = _arg.eval(0);
    return gsAsConstMatrix<T>(&m_value,1,1);
}
 
template<class T>
template<class E>
typename util::enable_if<!E::ScalarValued,gsAsConstMatrix<T> >::type
gsExprEvaluator<T>::evalIfc(const expr::_expr<E> & expr, const gsVector<T> & pt,
                            const boundaryInterface & ifc)
{
    m_exprdata->parse(expr);
 
    const bool flipSide = m_options.askSwitch("flipSide", false);
    const boundaryInterface & iFace =  flipSide ? ifc.getInverse() : ifc;
 
    const index_t patch1 = iFace.first().patch;
    const index_t patch2 = iFace.second().patch;
 
    ifacemap interfaceMap;
    if (iFace.type() == interaction::conforming)
        interfaceMap = gsAffineFunction<T>::make( iFace.dirMap(), iFace.dirOrientation(),
                                                  m_exprdata->multiBasis().basis(patch1).support(),
                                                  m_exprdata->multiBasis().basis(patch2).support() );
    else
        interfaceMap = gsCPPInterface<T>::make(m_exprdata->multiPatch(), m_exprdata->multiBasis(), iFace);
 
    m_exprdata->points() = pt;
    interfaceMap->eval_into(m_exprdata->points(), m_exprdata->pointsIfc());
 
    m_exprdata->precompute(iFace);
 
    // expr.printDetail(gsInfo); //after precompute
 
    gsMatrix<T> tmp = expr.eval(0);
    // const index_t r = expr.rows();
    // const index_t c = expr.cols();
    // gsInfo <<"tmp - "<< tmp.dim() <<" rc - "<< r <<", "<<c<<"\n";
    const index_t r = tmp.rows();
    const index_t c = tmp.cols();
    m_elWise.resize(r*c);
    gsAsMatrix<T>(m_elWise, r, c) = tmp; //expr.eval(0);
    return gsAsConstMatrix<T>(m_elWise, r, c);
}
 
// This is a copy of the above (commented out ) function, which has parts
// of gsParaviewCollection pasted in it. In this way the inclusion of
// gsParaviewCollection.h is prevented. This is a temporary modification.
template<class T>
template<class E>
typename util::enable_if<E::ScalarValued,gsAsConstMatrix<T> >::type
gsExprEvaluator<T>::evalBdr(const expr::_expr<E> & expr, const gsVector<T> & pt,
                            const patchSide & ps)
{
    GISMO_ASSERT(pt(ps.side().direction())==ps.side().parameter(),"Point "<<pt.transpose()<<" is not on boundary "<<ps.side());
    auto _arg = expr.val();
    m_exprdata->parse(_arg);
    m_elWise.clear();
 
    m_exprdata->points() = pt;
    m_exprdata->precompute(ps.patch, ps.side());
 
    m_value = _arg.eval(0);
    return gsAsConstMatrix<T>(&m_value,1,1);
}
 
template<class T>
template<class E>
typename util::enable_if<!E::ScalarValued,gsAsConstMatrix<T> >::type
gsExprEvaluator<T>::evalBdr(const expr::_expr<E> & expr, const gsVector<T> & pt,
                            const patchSide & ps)
{
    GISMO_ASSERT(pt(ps.side().direction())==ps.side().parameter(),"Point "<<pt.transpose()<<" is not on boundary "<<ps.side());
    m_exprdata->parse(expr);
    m_exprdata->points() = pt;
    m_exprdata->precompute(ps.patch, ps.side());
 
    gsMatrix<T> tmp = expr.eval(0);
    const index_t r = tmp.rows();
    const index_t c = tmp.cols();
    m_elWise.resize(r*c);
    gsAsMatrix<T>(m_elWise, r, c) = tmp; //expr.eval(0);
    return gsAsConstMatrix<T>(m_elWise, r, c);
}
 
template<class T>
template<class E, bool gmap>
void gsExprEvaluator<T>::writeParaview_impl(const expr::_expr<E> & expr,
                                            geometryMap G,
                                            std::string const & fn)
    {
        //if gmap is false, embed topology ?
        m_exprdata->parse(expr);
 
        //if false, embed topology ?
        const index_t n = m_exprdata->multiBasis().nBases();
 
        // Snippet from gsParaviewCollection
        // gsParaviewCollection collection(fn);
        std::stringstream file;
        int counter = 0;
        file <<"<?xml version=\"1.0\"?>\n";
        file <<"<VTKFile type=\"Collection\" version=\"0.1\">";
        file <<"<Collection>\n";
        // End snippet from gsParaviewCollection
 
        //const index_t n = G.source().nPieces();
        //gsParaviewCollection collection(fn);
 
        std::string fileName;
 
        gsMatrix<T> pts, vals, ab;
 
        const bool mesh = m_options.askSwitch("plot.elements");
 
        for ( index_t i=0; i != n; ++i )
        {
            fileName = fn + util::to_string(i);
            unsigned nPts = m_options.askInt("plot.npts", 1000);
            //ab = m_exprdata->multiBasis().piece(i).support();
            ab = G.source().piece(i).support();
            gsGridIterator<T,CUBE> pt(ab, nPts);
            eval(expr, pt, i);
            nPts = pt.numPoints();
            vals = allValues(m_elWise.size()/nPts, nPts);
 
            if (gmap) // Forward the points ?
            {
                eval(G, pt, i);
                pts = allValues(m_elWise.size()/nPts, nPts);
            }
 
            gsWriteParaviewTPgrid( gmap ? pts : pt.toMatrix(), // parameters
                                  vals,
                                  pt.numPointsCwise(), fileName );
 
            // Snippet from gsParaviewCollection
            // collection.addPart(fileName+ ".vts");
            GISMO_ASSERT(counter!=-1, "Error: collection has been already saved." );
            file << "<DataSet part=\""<< counter++ <<"\" file=\""<<fileName<<".vts"<<"\"/>\n";
            // End snippet from gsParaviewCollection
 
            if ( mesh )
            {
                gsMesh<T> msh(m_exprdata->multiBasis().basis(i), 2);
                static_cast<const gsGeometry<T>&>(G.source().piece(i)).evaluateMesh(msh);
                gsWriteParaview(msh, fileName + "_mesh", false);
                // Snippet from gsParaviewCollection
                // collection.addPart(fileName+ ".vtp");
                GISMO_ASSERT(counter!=-1, "Error: collection has been already saved." );
                file << "<DataSet part=\""<< counter++ <<"\" file=\""<<fileName<<"_mesh.vtp"<<"\"/>\n";
                // End snippet from gsParaviewCollection
            }
        }
 
        // Snippet from gsParaviewCollection
        // collection.save();
        GISMO_ASSERT(counter!=-1, "Error: gsParaviewCollection::save() already called." );
        file <<"</Collection>\n";
        file <<"</VTKFile>\n";
 
        std::string mfn = fn + ".pvd";
        // gsInfo << mfn << "\n";
        std::ofstream f( mfn.c_str() );
        GISMO_ASSERT(f.is_open(), "Error creating "<< mfn );
        f << file.rdbuf();
        f.close();
        file.str("");
        counter = -1;
        // End snippet from gsParaviewCollection
    }
 
} //namespace gismo