gsC1SurfBasisVertex.h

#pragma once


#include <gsUnstructuredSplines/src/gsC1SurfGluingData.h>
#include <gsUnstructuredSplines/src/gsC1SurfVisitorBasisVertex.h>



namespace gismo
{
    template<class T, class bhVisitor = gsG1ASVisitorBasisVertex<T>>
    class gsC1SurfBasisVertex : public gsAssembler<T>
    {
    public:
        typedef gsAssembler<T> Base;

    public:
        gsC1SurfBasisVertex(gsMultiPatch<T> mp, // Single Patch
                          gsMultiBasis<T> basis, // Single Basis
                          std::vector<bool> isBoundary,
                          gsMatrix<T> &Phi,
                          gsMatrix<T> gluingD)
                : m_mp(mp), m_basis(basis), m_isBoundary(isBoundary), m_Phi(Phi), m_gD(gluingD)
        {

            for (index_t dir = 0; dir < m_mp.parDim(); dir++) // For the TWO directions
            {
                // Computing the G1 - basis function at the edge
                // Spaces for computing the g1 basis
                gsBSplineBasis<T> basis_edge = dynamic_cast<gsBSplineBasis<T> &>(m_basis.basis(0).component(dir)); // 0 -> u, 1 -> v

                gsBSplineBasis<T> basis_plus(basis_edge);
                basis_plus.elevateContinuity(1);
                m_basis_plus.push_back(basis_plus);

                gsBSplineBasis<T> basis_minus(basis_edge);
                basis_minus.degreeReduce(1);
                m_basis_minus.push_back(basis_minus);
            }

            // Basis for the G1 basis
            m_basis_g1 = m_basis.basis(0);


        }


        void refresh();
        void assemble();
        inline void apply(bhVisitor & visitor, index_t patchIndex);
        void solve();

        void constructSolution(gsMultiPatch<T> & result);

        void setG1BasisVertex(gsMultiPatch<T> & result)
        {
            m_geo = m_basis_g1;
            refresh();
            assemble();
            solve();

            constructSolution(result);
        }


    private:
        // Avoid hidden overloads w.r.t. gsAssembler
        void constructSolution(const gsMatrix<T>& solVector,
                               gsMultiPatch<T>& result, short_t unk) const
        { GISMO_NO_IMPLEMENTATION; }

        void constructSolution(const gsMatrix<T>& solVector,
                               gsMultiPatch<T>& result,
                               const gsVector<index_t>  & unknowns) const
        { GISMO_NO_IMPLEMENTATION; }

    protected:

        // Input
        gsMultiPatch<T> m_mp;
        gsMultiBasis<T> m_basis;
        std::vector<bool> m_isBoundary;
        gsMatrix<T> m_Phi;

        // Gluing data
        gsMatrix<T> m_gD;

        // Basis for getting the G1 Basis
        std::vector<gsBSplineBasis<T>> m_basis_plus;
        std::vector<gsBSplineBasis<T>> m_basis_minus;

        // Basis for the G1 Basis
        gsMultiBasis<T> m_basis_g1;

        // Basis for Integration
        gsMultiBasis<T> m_geo;

        // System
        std::vector<gsSparseSystem<T> > m_f;

        // For Dirichlet boundary
        using Base::m_ddof;

        std::vector<gsMatrix<T>> solVec;

    }; // class gsG1BasisEdge


    template <class T, class bhVisitor>
    void gsC1SurfBasisVertex<T,bhVisitor>::constructSolution(gsMultiPatch<T> & result)
    {

        result.clear();

        // Dim is the same for all basis functions
        const index_t dim = ( 0!=solVec.at(0).cols() ? solVec.at(0).cols() :  m_ddof[0].cols() );

        gsMatrix<T> coeffs;
        for (index_t p = 0; p < 6; ++p)
        {
            const gsDofMapper & mapper = m_f.at(p).colMapper(0); // unknown = 0

            // Reconstruct solution coefficients on patch p
            index_t sz;
            sz = m_basis.basis(0).size();

            coeffs.resize(sz, dim);

            for (index_t i = 0; i < sz; ++i)
            {
                if (mapper.is_free(i, 0)) // DoF value is in the solVector // 0 = unitPatch
                {
                    coeffs.row(i) = solVec.at(p).row(mapper.index(i, 0));
                }
                else // eliminated DoF: fill with Dirichlet data
                {
                    //gsInfo << "mapper index dirichlet: " << m_ddof[unk].row( mapper.bindex(i, p) ).head(dim) << "\n";
                    coeffs.row(i) = m_ddof[0].row( mapper.bindex(i, 0) ).head(dim); // = 0
                }
            }
            result.addPatch(m_basis_g1.basis(0).makeGeometry(give(coeffs)));
        }
    }

    template <class T, class bhVisitor>
    void gsC1SurfBasisVertex<T,bhVisitor>::refresh()
    {
        // 1. Obtain a map from basis functions to matrix columns and rows
        gsDofMapper map(m_basis.basis(0));

        // SET THE DOFS

        map.finalize();


        // 2. Create the sparse system
        gsSparseSystem<T> m_system = gsSparseSystem<T>(map);
        for (index_t i = 0; i < 6; i++)
            m_f.push_back(m_system);

    } // refresh()

    template <class T, class bhVisitor>
    void gsC1SurfBasisVertex<T,bhVisitor>::assemble()
    {
        // Reserve sparse system
        const index_t nz = gsAssemblerOptions::numColNz(m_basis[0],2,1,0.333333);
        for (size_t i = 0; i < m_f.size(); i++)
            m_f.at(i).reserve(nz, 1);


        if(m_ddof.size()==0)
            m_ddof.resize(1); // 0,1

        const gsDofMapper & map = m_f.at(0).colMapper(0); // Map same for every

        m_ddof[0].setZero(map.boundarySize(), 1 ); // plus

        // Assemble volume integrals
        bhVisitor visitor;
        apply(visitor,0); // patch 0

        for (size_t i = 0; i < m_f.size(); i++)
            m_f.at(i).matrix().makeCompressed();

    } // assemble()

    template <class T, class bhVisitor>
    void gsC1SurfBasisVertex<T,bhVisitor>::apply(bhVisitor & visitor, index_t patchIndex)
    {
#pragma omp parallel
        {

            gsQuadRule<T> quRule ; // Quadrature rule
            gsMatrix<T> quNodes  ; // Temp variable for mapped nodes
            gsVector<T> quWeights; // Temp variable for mapped weights

            bhVisitor
#ifdef _OPENMP
            // Create thread-private visitor
    visitor_(visitor);
    const int tid = omp_get_thread_num();
    const int nt  = omp_get_num_threads();
#else
                    &visitor_ = visitor;
#endif

            gsBasis<T> & basis_g1 = m_basis_g1.basis(0); // basis for construction

            // Same for all patches
            gsBasis<T> & basis_geo = m_basis.basis(0); // 2D

            // Initialize reference quadrature rule and visitor data
            visitor_.initialize(basis_g1, quRule);

            const gsGeometry<T> & patch = m_mp.patch(0);

            // Initialize domain element iterator
            typename gsBasis<T>::domainIter domIt = m_geo.basis(0).makeDomainIterator(boundary::none);

#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(), quNodes, quWeights );
#pragma omp critical(evaluate)
                // Perform required evaluations on the quadrature nodes
                visitor_.evaluate(basis_g1, basis_geo, m_basis_plus, m_basis_minus, patch, quNodes, m_gD, m_isBoundary, m_Phi);

                // Assemble on element
                visitor_.assemble(*domIt, quWeights);

                // Push to global matrix and right-hand side vector
#pragma omp critical(localToGlobal)
                visitor_.localToGlobal(patchIndex, m_ddof, m_f); // omp_locks inside // patchIndex == 0
            }
        }//omp parallel
    } // apply

    template <class T, class bhVisitor>
    void gsC1SurfBasisVertex<T,bhVisitor>::solve()
    {
        typename gsSparseSolver<T>::SimplicialLDLT solver;
//    typename gsSparseSolver<T>::LU solver;


//    gsInfo << "rhs: " << m_f.at(4).rhs() << "\n";

        for (index_t i = 0; i < 6; i++) // Tilde
        {
            solver.compute(m_f.at(i).matrix());
            solVec.push_back(solver.solve(m_f.at(i).rhs()));
        }
    } // solve

} // namespace gismo