gsThinShellDWRHelper.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsKLShell/src/gsThinShellAssemblerDWR.h>
 
namespace gismo
{
 
template <class T>
class gsThinShellDWRHelper
 
{
public:
    typedef typename gsThinShellAssemblerDWRBase<T>::bContainer bContainer;
 
 
    // empty constructor
    gsThinShellDWRHelper() {};
 
 
    gsThinShellDWRHelper(gsThinShellAssemblerDWRBase<T> * assembler)
    :
    m_assembler(assembler)
    {
 
    }
 
    void computeError(const gsVector<T> & U)
    {
        gsMultiPatch<T> primalL, deformed;
 
        m_assembler->constructSolutionL(U,deformed);
        m_assembler->constructMultiPatchL(U,primalL);
 
        this->computeError(deformed,primalL);
    }
 
    void computeError(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> primalL, bool withLoads = false)
    {
        gsMatrix<T> points;
        bContainer bnds;
        this->computeError(deformed,primalL,bnds,points,true,withLoads);
    }
 
    void computeError(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> primalL, bool withLoads = false,
                        std::string filename = std::string(), unsigned np=1000, bool parametric=false, bool mesh=false)
    {
        gsMatrix<T> points;
        bContainer bnds;
        this->computeError(deformed,primalL,bnds,points,true,withLoads,filename,np,parametric,mesh);
    }
 
    void computeError(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> primalL,
                        const bContainer & bnds, const gsMatrix<T> & points, bool interior = true,bool withLoads=false,
                        std::string filename = std::string(), unsigned np=1000, bool parametric=false, bool mesh=false)
    {
        gsMultiPatch<T> dualL, dualH;
        gsVector<T> solVector;
        gsVector<T> rhsL(m_assembler->numDofsL());
        rhsL.setZero();
        gsVector<T> rhsH(m_assembler->numDofsH());
        rhsH.setZero();
 
        gsSparseSolver<real_t>::uPtr solver;
        #ifdef GISMO_WITH_PARDISO
            solver = gsSparseSolver<real_t>::get( "PardisoLU");
        #else
            solver = gsSparseSolver<real_t>::get( "LU");
        #endif
 
        gsInfo << "Assembling dual matrix (L)... "<< std::flush;
        m_assembler->assembleMatrixL(deformed);
        gsInfo << "done.\n";
        gsInfo << "Assembling dual vector (L)... "<< std::flush;
        if (interior)
        {
            m_assembler->assembleDualL(primalL,deformed);
            rhsL+=m_assembler->dualL();
        }
        if (bnds.size()!=0)
        {
            m_assembler->assembleDualL(bnds,primalL,deformed);
            rhsL+=m_assembler->dualL();
        }
        if (points.cols()!=0)
        {
            m_assembler->assembleDualL(points,primalL,deformed);
            rhsL+=m_assembler->dualL();
        }
        gsInfo << "done.\n";
        gsInfo << "Solving dual (L), size = "<<m_assembler->matrixL().rows()<<","<<m_assembler->matrixL().cols()<<"... "<< std::flush;
        solver->compute(m_assembler->matrixL());
        solVector = solver->solve(rhsL);
        m_assembler->constructMultiPatchL(solVector,dualL);
        gsInfo << "done.\n";
 
        gsInfo << "Assembling dual matrix (H)... "<< std::flush;
        m_assembler->assembleMatrixH(deformed);
        gsInfo << "done.\n";
        gsInfo << "Assembling dual vector (H)... "<< std::flush;
        if (interior)
        {
            m_assembler->assembleDualH(primalL,deformed);
            rhsH+=m_assembler->dualH();
        }
        if (bnds.size()!=0)
        {
            m_assembler->assembleDualH(bnds,primalL,deformed);
            rhsH+=m_assembler->dualH();
        }
        if (points.cols()!=0)
        {
            m_assembler->assembleDualH(points,primalL,deformed);
            rhsH+=m_assembler->dualH();
        }
 
        gsInfo << "done.\n";
        gsInfo << "Solving dual (H), size = "<<m_assembler->matrixH().rows()<<","<<m_assembler->matrixH().cols()<<"... "<< std::flush;
        solver->compute(m_assembler->matrixH());
        solVector = solver->solve(rhsH);
        m_assembler->constructMultiPatchH(solVector,dualH);
        gsInfo << "done.\n";
 
        m_error = m_assembler->computeError(dualL,dualH,deformed,withLoads,filename,np,parametric,mesh);
        m_errors = m_assembler->computeErrorElements(dualL,dualH,deformed,withLoads);
        m_sqerrors = m_assembler->computeSquaredErrorElements(dualL,dualH,deformed,withLoads);
    }
 
    T error() const { return  m_error; }
    std::vector<T> errors(bool normalize=false) const
    {
        if (normalize)
        {
            std::vector<T> result = m_errors;
            for (typename std::vector<T>::iterator it = result.begin(); it!=result.end(); it++)
                *it = std::abs(*it) / std::abs(m_error);
            return result;
        }
        else
            return m_errors;
 
    }
    std::vector<T> absErrors(bool normalize=false) const
    { 
        std::vector<T> result = m_errors;
        if (normalize)
        {
            for (typename std::vector<T>::iterator it = result.begin(); it!=result.end(); it++)
                *it = std::abs(*it)  / std::abs(m_error);
        }
        else
        {
            for (typename std::vector<T>::iterator it = result.begin(); it!=result.end(); it++)
                *it = std::abs(*it);
        }
        return result;
    }
    std::vector<T> sqErrors(bool normalize=false) const
    {
        if (normalize)
        {
            std::vector<T> result = m_sqerrors;
            for (typename std::vector<T>::iterator it = result.begin(); it!=result.end(); it++)
                *it = std::abs(*it) / (math::pow(m_error,2));
            return result;
        }
        else
            return m_sqerrors;
    }
 
protected:
    // bool m_verbose;
    gsThinShellAssemblerDWRBase<T> * m_assembler;
    T m_error;
    std::vector<T> m_errors, m_sqerrors;
};
 
} // namespace gismo