gsVisitorCDR.h

Go to the documentation of this file.

#pragma once

namespace gismo
{

template <class T>
class gsVisitorCDR
{
public:

    gsVisitorCDR(const gsPde<T> & pde)
    {
        const gsConvDiffRePde<T>* cdr =
            static_cast<const gsConvDiffRePde<T>*>(&pde);

        coeff_A_ptr = cdr->diffusion ();
        coeff_b_ptr = cdr->convection();
        coeff_c_ptr = cdr->reaction  ();
        rhs_ptr     = cdr->rhs       ();

        flagStabType = stabilizerCDR::none;

        GISMO_ASSERT( rhs_ptr->targetDim() == 1 ,
                      "Not yet tested for multiple right-hand-sides");
    }

    gsVisitorCDR(const gsFunction<T> & rhs,
                 const gsFunction<T> & coeff_A,
                 const gsFunction<T> & coeff_b,
                 const gsFunction<T> & coeff_c,
                 stabilizerCDR::method flagStabilization = stabilizerCDR::SUPG) :
        rhs_ptr(&rhs),
        coeff_A_ptr( & coeff_A),coeff_b_ptr( & coeff_b),coeff_c_ptr( & coeff_c),
        flagStabType( flagStabilization )
    {
        GISMO_ASSERT( rhs.targetDim() == 1 ,"Not yet tested for multiple right-hand-sides");
        GISMO_ASSERT( flagStabilization == stabilizerCDR::none || flagStabilization == stabilizerCDR::SUPG, "flagStabilization not known");
    }

    void initialize(const gsBasis<T> & basis,
                    const index_t ,
                    const gsOptionList & options,
                    gsQuadRule<T>    & rule)
    {
        // Setup Quadrature
        rule = gsQuadrature::get(basis, options); // harmless slicing occurs here

        //flagStabType = static_cast<unsigned>(options.askSwitch("SUPG", false));
        flagStabType = static_cast<stabilizerCDR::method>(options.askInt("Stabilization", stabilizerCDR::none));

        // Set Geometry evaluation flags
        md.flags = NEED_VALUE | NEED_MEASURE | NEED_GRAD_TRANSFORM | NEED_2ND_DER;
    }

    inline void evaluate(const gsBasis<T>       & basis, // to do: more unknowns
                         const gsGeometry<T>    & geo,
                         const gsMatrix<T>      & quNodes)
    {
        base = &geo;
        md.points = quNodes;

        // Compute the active basis functions
        // Assumes actives are the same for all quadrature points on the elements
        basis.active_into(md.points.col(0), actives);
        numActive = actives.rows();

        // Evaluate basis functions on element
        basis.evalAllDers_into(md.points, 2, basisData, true);

        // Compute image of Gauss nodes under geometry mapping as well as Jacobians
        geo.computeMap(md);

        // Evaluate the coefficients
        coeff_A_ptr->eval_into(md.values[0], coeff_A_vals);
        coeff_b_ptr->eval_into(md.values[0], coeff_b_vals);
        coeff_c_ptr->eval_into(md.values[0], coeff_c_vals);

        // Evaluate right-hand side at the geometry points
        rhs_ptr->eval_into(md.values[0], rhsVals); // to do: parametric rhs ?

        // Initialize local matrix/rhs
        localMat.setZero(numActive, numActive);
        localRhs.setZero(numActive, rhsVals.rows());//multiple right-hand sides
    }

    inline void assemble(gsDomainIterator<T>    & element,
                         const gsVector<T>      & quWeights)
    {

        const index_t N = numActive;

        gsMatrix<T> & basisVals  = basisData[0];
        gsMatrix<T> & basisGrads = basisData[1];
        gsMatrix<T> & basis2ndDerivs = basisData[2];

        const unsigned d = element.dim();

        // supgMat will contain the contributions to the assembled matrix
        // that come from the SUPG stabilization. It is initialized whether
        // SUPG-stabilization is used or not, because the SUPG-parameter
        // has to be computed AFTER the loop over the quadrature points
        // (see below).
        gsMatrix<T> supgMat( localMat.rows(), localMat.cols() );
        supgMat.setZero();

        for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
        {
            // Multiply weight by the geometry measure
            const T weight = quWeights[k] * md.measure(k);

            // Compute physical gradients at k as a Dim x numActive matrix
            transformGradients   (md, k, basisGrads, physBasisGrad);
            transformDeriv2Hgrad (md, k, basisGrads, basis2ndDerivs, physBasisd2);

            // d ... dim
            // N ... numActive

            // physBasisGrad : d x N
            // A.col(k)      : d^2 x 1
            // tmp_A         : d x d
            gsMatrix<T> tmp_A = coeff_A_vals.col(k);
            tmp_A.resize(d,d);

            // b.col(k)      : d x 1
            // physBasisGrad : d x N
            // tmp_b         : 1 x N
            gsMatrix<T> b_basisGrads = coeff_b_vals.col(k).transpose() * physBasisGrad;

            // basisVals.col(k): N x 1
            // rhsVals.col(k)  : 1 x 1
            // result:         : N x 1
            localRhs.noalias() += weight * ( basisVals.col(k) * rhsVals.col(k).transpose() ) ;

            // ( N x d ) * ( d x d ) * ( d x N ) = N x N
            localMat.noalias() += weight * (physBasisGrad.transpose() * ( tmp_A * physBasisGrad) );
            // ( N x 1 ) * ( 1 x N) = N x N
            localMat.noalias() += weight * (basisVals.col(k) * b_basisGrads);
            // ( scalar ) * ( N x 1 ) * ( 1 x N ) = N x N
            localMat.noalias() += weight * coeff_c_vals(0,k) * (basisVals.col(k) * basisVals.col(k).transpose());


            if( flagStabType == stabilizerCDR::SUPG ) // 1: SUPG
            {
                //const typename gsMatrix<T>::constColumns J = geoEval.jacobian(k); //todo: correct?
                const typename gsFuncData<T>::matrixTransposeView J = md.jacobian(k);
                gsMatrix<T> Jinv = J.inverse();

                gsMatrix<T> grad_b_basisGradsT(N,d);
                grad_b_basisGradsT.setZero();

                // loop over all basis functions
                for( index_t fct_i = 0; fct_i < N; ++fct_i  )
                {
                    gsMatrix<T> tmp_basis2ndDerivs(d,d);
                    tmp_basis2ndDerivs.setZero();
                    if( d == 2 )
                    {
                        tmp_basis2ndDerivs(0,0) = physBasisd2(fct_i, 0);
                        tmp_basis2ndDerivs(1,1) = physBasisd2(fct_i, 1);
                        tmp_basis2ndDerivs(0,1) =
                        tmp_basis2ndDerivs(1,0) = physBasisd2(fct_i, 2);
                    }
                    else if (d == 3 )
                    {
                        tmp_basis2ndDerivs(0,0) = physBasisd2(fct_i, 0);
                        tmp_basis2ndDerivs(1,1) = physBasisd2(fct_i, 1);
                        tmp_basis2ndDerivs(2,2) = physBasisd2(fct_i, 2);
                        tmp_basis2ndDerivs(0,1) =
                        tmp_basis2ndDerivs(1,0) = physBasisd2(fct_i, 3);
                        tmp_basis2ndDerivs(0,2) =
                        tmp_basis2ndDerivs(2,0) = physBasisd2(fct_i, 4);
                        tmp_basis2ndDerivs(1,2) =
                        tmp_basis2ndDerivs(2,1) = physBasisd2(fct_i, 5);
                    }
                    else
                    {
                        GISMO_ERROR("What kind of dimension are you using? Should be 2 or 3.");
                    }

                    //grad_b_basisGradsT.row(fct_i) = (tmp_basis2ndDerivs * coeff_b_vals.col(k) ).transpose();
                    for( unsigned i = 0; i < d; ++i )
                        for( unsigned j = 0; j < d; ++j )
                            grad_b_basisGradsT(fct_i, i) += coeff_b_vals(j,k) * tmp_basis2ndDerivs(j,i);
                }

                supgMat.noalias() += weight * grad_b_basisGradsT * ( tmp_A * physBasisGrad );
                supgMat.noalias() += weight * (b_basisGrads.transpose() * b_basisGrads);
                supgMat.noalias() += weight * coeff_c_vals(0,k) * ( b_basisGrads.transpose() * basisVals.col(k).transpose());

            }
        }

        if( flagStabType == stabilizerCDR::SUPG ) // 1: SUPG
        {
            // Calling getSUPGParameter re-evaluates the (*base) geometry. // todo: is that correct so?
            // Thus, it has to be called AFTER geo (*base) has been used.
            T supgParam = getSUPGParameter( element.lowerCorner(),
                                            element.upperCorner());
            // Add the contributions from the SUPG-stabilization.
            localMat.noalias() += supgParam * supgMat;
        }
    }

    inline void localToGlobal(const index_t                     patchIndex,
                              const std::vector<gsMatrix<T> > & eliminatedDofs,
                              gsSparseSystem<T>               & system)
    {
        // Map patch-local DoFs to global DoFs
        system.mapColIndices(actives, patchIndex, actives);

        // Add contributions to the system matrix and right-hand side
        system.push(localMat, localRhs, actives, eliminatedDofs.front(), 0, 0);
    }

    T getSUPGParameter( const gsVector<T> & lo,
                        const gsVector<T> & up)
    {
        const index_t N = 2;

        const index_t d = lo.size();

        gsMatrix<T> b_at_phys_pts;

        // compute the center point of the cell...

        // ...get the points map it to the physical space...
        gsMatrix<T> phys_pts = md.values[0];
        // ...evaluate the convection coefficient there, ...
        coeff_b_ptr->eval_into( phys_pts, b_at_phys_pts );
        // ...and get it's norm.
        T b_norm = 0;
        for( index_t i=0; i < d; i++)
            b_norm += b_at_phys_pts(i,0) * b_at_phys_pts(i,0);
        b_norm = math::sqrt( b_norm );

        T SUPG_param = (T)(0.0);
        if( b_norm > 0 )
        {
            gsMatrix<T> aMat;

            if( d == 2 )
            {
                index_t N1 = N+1;
                md.points.resize( 2, 4*N1 );
                aMat.resize( 2, 4*N1 );

                for( index_t i = 0; i <= N; ++i )
                {
                    T a = (T)(i)/(T)(N);
                    aMat(0,i) = a;
                    aMat(1,i) = (T)(0.0);
                    aMat(0,i+N1) = a;
                    aMat(1,i+N1) = (T)(1.0);

                    aMat(0,i+2*N1) = (T)(0.0);
                    aMat(1,i+2*N1) = a;
                    aMat(0,i+3*N1) = (T)(1.0);
                    aMat(1,i+3*N1) = a;
                }
            }
            else if( d == 3 )
            {

                GISMO_ASSERT(false,"NOT IMLEMENTED YET, Mark m271");

                /*
                md.points.resize( 3, 6*(N+1)*(N+1) );
                aMat.resize( 3, 6*(N+1)*(N+1) );

                index_t N1 = N+1;
                md.points.resize( 2, 4*N1 );
                aMat.resize( 2, 4*N1 );

                index_t ij = 0;
                for( index_t i = 0; i <= N; ++i )
                    for( index_t j = 0; j <= N; ++j )
                    {
                        T ai = (T)(i)/(T)(N);
                        T aj = (T)(j)/(T)(N);
                        aMat(0,ij) = (T)(0.0);
                        aMat(1,ij) = (T)(0.0);
                        aMat(2,ij) = (T)(0.0);
                        aMat(0,ij+N1) = (T)(1.0);
                        aMat(1,ij+N1) = (T)(1.0);
                        aMat(1,ij+N1) = (T)(1.0);

                    }
                */

            }
            else
            {
                GISMO_ASSERT(false,"WRONG DIMENSION. Mark m243");
            }


            for( index_t di = 0; di < d; ++di )
                for( index_t i = 0; i < aMat.cols(); ++i)
                {
                    md.points(di,i) = ( (T)(1) - aMat(di,i) )*lo[di] + aMat(di,i) * up[di];
                }

            base->computeMap(md);
            gsMatrix<T> b_proj = md.values[0].transpose() * b_at_phys_pts;

            T b_proj_min = b_proj(0,0);
            T b_proj_max = b_proj(0,0);
            for( index_t i = 0; i < b_proj.size(); i++)
            {
                if( b_proj_min > b_proj(i) )
                    b_proj_min = b_proj(i);
                if( b_proj_max < b_proj(i) )
                    b_proj_max = b_proj(i);
            }
            SUPG_param = ( b_proj_max - b_proj_min ) / ( (T)(2) * b_norm );
        }

        return SUPG_param;
    }



protected:
    // Right hand side
    const gsFunction<T> * rhs_ptr;
    // PDE Coefficient
    const gsFunction<T> * coeff_A_ptr;
    const gsFunction<T> * coeff_b_ptr;
    const gsFunction<T> * coeff_c_ptr;
    // flag for stabilization method
    stabilizerCDR::method flagStabType;

protected:
    // Basis values
    std::vector<gsMatrix<T> > basisData;
    gsMatrix<T>        physBasisGrad, physBasisd2;
    gsMatrix<index_t> actives;
    index_t numActive;

    gsMatrix<T> coeff_A_vals;
    gsMatrix<T> coeff_b_vals;
    gsMatrix<T> coeff_c_vals;

protected:
    // Local values of the right hand side
    gsMatrix<T> rhsVals;

protected:
    // Local matrices
    gsMatrix<T> localMat;
    gsMatrix<T> localRhs;

    const gsGeometry<T> * base;
    gsMapData<T> md;
};


} // namespace gismo