gsVisitorNitscheBiharmonic.h

Go to the documentation of this file.

#pragma once

namespace gismo
{

template <class T>
class gsVisitorNitscheBiharmonic
{
public:

    gsVisitorNitscheBiharmonic(const gsFunction<T> & dirdata, T penalty, boxSide s)
    : dirdata_ptr(&dirdata), m_penalty(penalty), side(s)
    { }

    void initialize(const gsBasis<T> & basis,
                    gsQuadRule<T>    & rule)
    {
        const int dir = side.direction();
        gsVector<int> numQuadNodes(basis.dim());
        for (int i = 0; i < basis.dim(); ++i)
            numQuadNodes[i] = basis.degree(i) + 1;
        numQuadNodes[dir] = 1;

        // Setup Quadrature
        rule = gsGaussRule<T>(numQuadNodes);// harmless slicing occurs here

        // 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,
                         // todo: add element here for efficiency
                         gsMatrix<T>         & quNodes)
    {
        md.points = quNodes;
        // Compute the active basis functions
        // Assumes actives are the same for all quadrature points on the current element
        basis.active_into(md.points.col(0), actives);
        const index_t numActive = actives.rows();

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

        // Compute geometry related values
        geo.computeMap(md);

        // Evaluate the Dirichlet data
        dirdata_ptr->eval_into(md.values[0], dirData);

        // Initialize local matrix/rhs
        localMat.setZero(numActive, numActive);
        localRhs.setZero(numActive, dirdata_ptr->targetDim());
    }

    inline void assemble(gsDomainIterator<T> & element,
                         const gsVector<T>   & quWeights)
    {
        const unsigned d = element.dim();

        const index_t numActive = actives.rows();

        for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
        {
            //const typename gsMatrix<T>::Block basisVals  = basisData.topRows(numActive);
            const typename gsMatrix<T>::Block basisGrads =
                basisData.middleRows(numActive, numActive * d);
            const typename gsMatrix<T>::Block basis2ndDerivs =
                basisData.bottomRows(numActive * (d * (d + 1)) / 2);

            // Compute the outer normal vector on the side
            outerNormal(md, k, side, unormal);

            // Multiply quadrature weight by the geometry measure
            const T weight = quWeights[k] * unormal.norm();

            // Compute the unit normal vector : Dim x 1
            unormal.normalize();

            // Compute physical gradients at k as a Dim x NumActive matrix
            transformGradients(md, k, basisGrads, physBasisGrads);
            // Compute physical laplacian at k as a 1 x numActive matrix
            transformLaplaceHgrad(md, k, basisGrads, basis2ndDerivs, physBasisLaplace);

            // Get penalty parameter
            const T h = element.getCellSize();
            const T mu = m_penalty / (0 != h ? h : 1);

            // Sum up quadrature point evaluations
            localRhs.noalias() += weight * ((physBasisLaplace.transpose() + mu * physBasisGrads.transpose() * unormal)
                * dirData.col(k).transpose());

            localMat.noalias() += weight * (physBasisGrads.transpose() * unormal * physBasisLaplace
                + (physBasisGrads.transpose() * unormal * physBasisLaplace).transpose()
                - mu * physBasisGrads.transpose() * physBasisGrads);
        }
    }

    void localToGlobal(const gsDofMapper & mapper,
                       const gsMatrix<T> & eliminatedDofs,
                       const index_t       patchIndex,
                       gsSparseMatrix<T> & sysMatrix,
                       gsMatrix<T>       & rhsMatrix )
    {
        // Local DoFs to global DoFs
        mapper.localToGlobal(actives, patchIndex, actives);
        const index_t numActive = actives.rows();

        // Push element contribution to the global matrix and load vector
        for (index_t j=0; j!=numActive; ++j)
        {
            const unsigned jj = actives(j);
            rhsMatrix.row(jj) -= localRhs.row(j);
            for (index_t i=0; i!=numActive; ++i)
            {
                const unsigned ii = actives(i);
                if ( jj <= ii ) // assuming symmetric problem (!) probably we must not.
                    sysMatrix( ii, jj ) -= localMat(i,j);
            }
        }

    }

private:
    // Dirichlet function
    const gsFunction<T> * dirdata_ptr;

    // Penalty constant
    T m_penalty;

    // Side
    boxSide side;

private:
    // Basis values
    gsMatrix<T>       basisData;
    gsMatrix<T>       physBasisGrads;
    gsMatrix<T>       physBasisLaplace;
    gsMatrix<index_t> actives;

    // Normal and Neumann values
    gsVector<T> unormal;
    gsMatrix<T> dirData;

    // Local  matrix and rhs
    gsMatrix<T> localMat;
    gsMatrix<T> localRhs;

    gsMapData<T> md;
};


} // namespace gismo