gsApproxC1GluingData.h

 
#pragma once
 
#include <gsUnstructuredSplines/src/gsApproxC1Utils.h>
 
namespace gismo
{
 
template<short_t d, class T>
class gsApproxC1GluingData
{
private:
    typedef typename std::vector<gsPatchReparameterized<d,T>> C1AuxPatchContainer;
 
    typedef memory::shared_ptr<gsApproxC1GluingData> Ptr;
 
    typedef memory::unique_ptr<gsApproxC1GluingData> uPtr;
 
public:
    gsApproxC1GluingData()
    { }
 
 
    gsApproxC1GluingData(C1AuxPatchContainer const & auxPatchContainer,
                       gsOptionList const & optionList,
                       std::vector<patchSide> sidesContainer,
                       std::vector<bool> isInterface = std::vector<bool>{},
                       gsTensorBSplineBasis<d, T> basis = gsTensorBSplineBasis<d, T>())
        : m_auxPatches(auxPatchContainer), m_optionList(optionList)
    {
        alphaSContainer.resize(2);
        betaSContainer.resize(2);
        if (m_auxPatches.size() == 2) // Interface
        {
            setGlobalGluingData(1,0); // u
            setGlobalGluingData(0,1); // v
        }
        else if (m_auxPatches.size() == 1 && sidesContainer.size() == 2) // Vertex
        {
            for (size_t dir = 0; dir < sidesContainer.size(); dir++)
            {
                // Map global side to local side
                index_t localSide = auxPatchContainer[0].getMapIndex(sidesContainer[dir].index());
                //gsInfo << "Global: " << sidesContainer[dir] << " : " << localSide << "\n";
                index_t localDir = localSide < 3 ? 1 : 0;
 
                createGluingDataSpace(m_auxPatches[0].getPatchRotated(), basis,
                                      localDir, bsp_gD, m_optionList.getInt("gluingDataDegree"), m_optionList.getInt("gluingDataSmoothness"));
 
                if(isInterface[dir]) // West
                {
                    setGlobalGluingData(0, localDir);
                }
 
                else
                {
                    // empty
                }
            }
        }
        //else
        //    gsInfo << "I am here \n";
 
    }
 
    // Computed the gluing data globally
    void setGlobalGluingData(index_t patchID = 0,  index_t dir = 1);
 
    gsBSpline<T> & alphaS(index_t patchID) { return alphaSContainer[patchID]; }
    gsBSpline<T> & betaS(index_t patchID) { return betaSContainer[patchID]; }
 
protected:
 
    // Spline space for the gluing data + multiPatch
    C1AuxPatchContainer m_auxPatches;
 
    gsBSplineBasis<T> bsp_gD;
 
    const gsOptionList m_optionList;
 
    // Result
    std::vector<gsBSpline<T>> alphaSContainer, betaSContainer;
 
}; // class gsGluingData
 
 
template<short_t d, class T>
void gsApproxC1GluingData<d, T>::setGlobalGluingData(index_t patchID, index_t dir)
{
    // Interpolate boundary yes or no //
    bool interpolate_boundary = false;
    // Interpolate boundary yes or no //
 
    gsTensorBSplineBasis<d, T> basis = dynamic_cast<const gsTensorBSplineBasis<d, T> &>(m_auxPatches[patchID].getBasisRotated().piece(0));;
    if (m_auxPatches.size() == 2) // Interface
    {
        gsTensorBSplineBasis<d, T> basis2 = dynamic_cast<const gsTensorBSplineBasis<d, T> &>(m_auxPatches[1-patchID].getBasisRotated().piece(0));
        if (basis.component(dir).numElements() > basis2.component(1-dir).numElements())
            basis.component(dir) = basis2.component(1-dir);
 
        // ======== Space for gluing data : S^(p_tilde, r_tilde) _k ========
        createGluingDataSpace(m_auxPatches[patchID].getPatchRotated(), basis,
                              dir, bsp_gD, m_optionList.getInt("gluingDataDegree"), m_optionList.getInt("gluingDataSmoothness"));
    }
 
 
 
    typename gsSparseSolver<T>::SimplicialLDLT solver;
    gsExprAssembler<T> A(1,1);
 
    // Elements used for numerical integration
    gsMultiBasis<T> BsplineSpace(bsp_gD);
    A.setIntegrationElements(BsplineSpace);
    gsExprEvaluator<T> ev(A);
 
    gsAlpha<T> alpha(m_auxPatches[patchID].getPatchRotated(), dir);
    auto aa = A.getCoeff(alpha);
 
    // Set the discretization space
    auto u = A.getSpace(BsplineSpace);
 
    // Create Mapper
    gsDofMapper map(BsplineSpace);
    gsMatrix<index_t> act(2,1);
    act(0,0) = 0;
    act(1,0) = BsplineSpace[0].size()-1; // First and last
    if (interpolate_boundary)
        map.markBoundary(0, act); // Patch 0
    map.finalize();
 
    u.setupMapper(map);
 
    gsMatrix<T> & fixedDofs = const_cast<expr::gsFeSpace<T>&>(u).fixedPart();
    fixedDofs.setZero( u.mapper().boundarySize(), 1 );
 
    // For the boundary
    gsMatrix<T> points_bdy(1,2);
    points_bdy << 0.0, 1.0;
    if (interpolate_boundary)
        fixedDofs = alpha.eval(points_bdy).transpose();
 
    A.initSystem();
    A.assemble(u * u.tr(), u * aa);
 
    solver.compute( A.matrix() );
    gsMatrix<T> solVector = solver.solve(A.rhs());
 
    auto u_sol = A.getSolution(u, solVector);
    gsMatrix<T> sol;
    u_sol.extract(sol);
 
    typename gsGeometry<T>::uPtr tilde_temp;
    tilde_temp = bsp_gD.makeGeometry(sol);
    alphaSContainer[dir] = dynamic_cast<gsBSpline<T> &> (*tilde_temp);
 
    gsBeta<T> beta(m_auxPatches[patchID].getPatchRotated(), dir);
    auto bb = A.getCoeff(beta);
 
    // For the boundary
    if (interpolate_boundary)
        fixedDofs = beta.eval(points_bdy).transpose();
 
    A.initSystem();
    A.assemble(u * u.tr(), u * bb);
 
    solver.compute( A.matrix() );
    solVector = solver.solve(A.rhs());
 
    auto u_sol2 = A.getSolution(u, solVector);
    u_sol2.extract(sol);
 
    tilde_temp = bsp_gD.makeGeometry(sol);
    betaSContainer[dir] = dynamic_cast<gsBSpline<T> &> (*tilde_temp);
 
} // setGlobalGluingData
 
 
} // namespace gismo