gsBoundaryConditions.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsBoundary.h>


namespace gismo
{

struct condition_type
{
    enum type
    {
        unknownType = -1,
        dirichlet = 0, 
        weak_dirichlet = 10, 
        neumann   = 1, 
        robin     = 2, 
        clamped   = 3, 
        weak_clamped = 30,
        collapsed = 4,
        laplace = 5 
        //mixed BD means: there are both dirichlet and neumann sides
        //robin: a linear combination of value and derivative
        //cauchy: there are 2 conditions (value+deriv) defined on the same side
    };

};

// Print (as string) a boundary type
inline std::ostream &operator<<(std::ostream &os, const condition_type::type& o)
{
    switch (o)
    {
        case condition_type::dirichlet:
        {
            os<< "Dirichlet";
            break;
        }
        case condition_type::weak_dirichlet:
        {
            os<< "Weak Dirichlet";
            break;
        }
        case condition_type::neumann:
        {
            os<< "Neumann";
            break;
        }
        case condition_type::robin:
        {
            os<< "Robin";
            break;
        }
        case condition_type::clamped:
        {
            os<< "Clamped";
            break;
        }
        case condition_type::weak_clamped:
        {
            os<< "Weak Clamped";
            break;
        }
        case condition_type::collapsed:
        {
            os<< "Collapsed";
            break;
        }
        case condition_type::laplace:
        {
            os<< "Laplace";
            break;
        }
        default:
            os<< "condition type not known.\n";
    };
    return os;
}

template<class T>
struct boundary_condition
{
    typedef typename gsFunctionSet<T>::Ptr function_ptr;

    boundary_condition( int p, boxSide s, const function_ptr & f_shptr,
                        const std::string & label, short_t unknown,
                        short_t unkcomp, bool parametric)
            : ps(p, s),
              m_function(f_shptr),
              m_label(label),
              m_unknown(unknown),
              m_unkcomp(unkcomp),
              m_parametric(parametric)
    {
        if (m_label == "Dirichlet") m_type = condition_type::dirichlet;
        else if (m_label == "Weak Dirichlet") m_type = condition_type::weak_dirichlet;
        else if (m_label == "Neumann")   m_type = condition_type::neumann;
        else if (m_label == "Robin")     m_type = condition_type::robin;
        else if (m_label == "Clamped")   m_type = condition_type::clamped;
        else if (m_label == "Weak Clamped")   m_type = condition_type::weak_clamped;
        else if (m_label == "Collapsed") m_type = condition_type::collapsed;
        else if (m_label == "Laplace") m_type = condition_type::laplace;
        else m_type = condition_type::unknownType;
    }

    boundary_condition( int p, boxSide s, const function_ptr & f_shptr,
                        condition_type::type t, short_t unknown, bool parametric)
            : ps(p, s),
              m_function(f_shptr),
              m_type(t),
              m_unknown(unknown),
              m_unkcomp(-1),
              m_parametric(parametric)
    {
        switch (t)
        {
            case condition_type::dirichlet:
            {
                GISMO_ASSERT(!m_function || m_function->targetDim()==1,"Expecting scalar function");
                m_label = "Dirichlet";
                break;
            }
            case condition_type::weak_dirichlet:
            {
                GISMO_ASSERT(!m_function || m_function->targetDim()==1,"Expecting scalar function");
                m_label = "Weak Dirichlet";
                break;
            }
            case condition_type::neumann:
            {
                m_label = "Neumann";
                break;
            }
            case condition_type::robin:
            {
                m_label = "Robin";
                break;
            }
            case condition_type::clamped:
            {
                m_label = "Clamped";
                break;
            }
            case condition_type::weak_clamped:
            {
                m_label = "weak Clamped";
                break;
            }
            case condition_type::collapsed:
            {
                m_label = "Collapsed";
                break;
            }
            case condition_type::laplace:
            {
                m_label = "Laplace";
                break;
            }
            default:
                m_label = "Unknown";
                break;
        };
    }

    boundary_condition( int p, boxSide s, const function_ptr & f_shptr,
                        condition_type::type t, int unknown, int unkcomp, bool parametric)
            : ps(p, s),
              m_function(f_shptr),
              m_type(t),
              m_unknown(unknown),
              m_unkcomp(unkcomp),
              m_parametric(parametric)
    {
        switch (t)
        {
            case condition_type::dirichlet:
            {
                GISMO_ASSERT(!m_function || m_function->targetDim()==1,"Expecting scalar function");
                m_label = "Dirichlet";
                break;
            }
            case condition_type::weak_dirichlet:
            {
                // GISMO_ASSERT(!m_function || m_function->targetDim()==1,"Expecting scalar function");
                m_label = "Weak Dirichlet";
                break;
            }
            case condition_type::neumann:
            {
                m_label = "Neumann";
                break;
            }
            case condition_type::robin:
            {
                m_label = "Robin";
                break;
            }
            case condition_type::clamped:
            {
                m_label = "Clamped";
                break;
            }
            case condition_type::weak_clamped:
            {
                m_label = "Weak Clamped";
                break;
            }
            case condition_type::collapsed:
            {
                m_label = "Collapsed";
                break;
            }
            case condition_type::laplace:
            {
                m_label = "Laplace";
                break;
            }
            default:
                gsWarn<<"gsBoundaryConditions: Unknown boundary condition.\n";
                m_label = "Unknown";
                break;
        };
    }

    bool isHomogeneous() const { return m_function.get() == NULL; }

    function_ptr function() const { return m_function; }

    // Returns a reference to the function data
    //const gsFunctionSet<T> & function() const { return *m_function; }

    condition_type::type  type() const { return m_type; }

    const std::string & ctype() const { return m_label; }

    index_t patch()    const { return ps.patch; }

    boxSide side()     const { return ps.side(); }

    short_t     unknown()  const { return m_unknown; }

    short_t     unkComponent()  const { return m_unkcomp; }

    bool    parametric()  const { return m_parametric; }

    bool isSame(const boundary_condition & other) {return ((other.ps == ps) && (other.unknown() == unknown()) && (other.unkComponent() == unkComponent())) ;}

    patchSide ps;                

    function_ptr m_function;     

    // TO DO : robin coefficients?

    condition_type::type m_type;// todo: remove

    std::string m_label;         

    short_t m_unknown;               

    short_t m_unkcomp;               

    bool m_parametric;
};

    template<class T>
    struct corner_value
    {
        corner_value(index_t p, boxCorner c, T v, short_t unk = 0, int comp = -1)
                : patch(p), corner(c), value(v), unknown(unk), component(comp) { }

        index_t patch;     
        boxCorner corner; 
        T value;          
        short_t   unknown;    
        int   component;    
    };

    template<class T>
    struct coupled_boundary
    {
        coupled_boundary(index_t p1, boxSide s1, index_t p2, boxSide s2, short_t dim, short_t unk = 0, int comp = -1)
                : ifc(patchSide(p1,s1), patchSide(p2,s2), dim), unknown(unk), component(comp)
                {}

        boundaryInterface ifc;
        short_t   unknown;    
        int   component;    
    };

template<class T>
class GISMO_EXPORT gsBoundaryConditions
{

public:

    typedef typename std::deque<boundary_condition<T> > bcContainer;
    typedef typename bcContainer::iterator iterator;
    typedef typename bcContainer::const_iterator const_iterator;

    typedef typename std::deque<corner_value<T> >       cornerContainer;
    typedef typename cornerContainer::iterator citerator;
    typedef typename cornerContainer::const_iterator const_citerator;

    typedef typename std::deque<coupled_boundary<T>> cplContainer;
    typedef typename cplContainer::iterator cpliterator;
    typedef typename cplContainer::const_iterator const_cpliterator;

    typedef typename std::deque<boundaryInterface> ppContainer;
    typedef typename ppContainer::iterator ppiterator;
    typedef typename ppContainer::const_iterator const_ppiterator;

    // Format: std::pair<type,bcContainer>
    typedef std::map<std::string,bcContainer> bcData;
    typedef typename bcData::iterator bciterator;
    typedef typename bcData::const_iterator const_bciterator;

    typedef memory::shared_ptr<gsBoundaryConditions> Ptr;
    typedef memory::unique_ptr<gsBoundaryConditions> uPtr;

    typedef typename boundary_condition<T>::function_ptr function_ptr;

    typedef std::list<util::reference_wrapper<const boundary_condition<T> > > bcRefList;
public:

    gsBoundaryConditions() : m_patches(nullptr) { }

    /*
    gsBoundaryConditions & operator= (uPtr other)
    {
        if ( other.get() != NULL )
        {
            this->swap(*other);
            other.reset();
        }
        return *this;
    }
    */

    void swap(gsBoundaryConditions & other)
    {
        m_bc.swap(other.m_bc);
        corner_values.swap(other.corner_values);
        m_periodicPairs.swap(other.m_periodicPairs);
        coupled_boundaries.swap(other.coupled_boundaries);
    }

public:

    void clear()
    {
        m_bc.clear();
        corner_values.clear();
        m_periodicPairs.clear();
        coupled_boundaries.clear();
    }

    size_t size() const
    {
        size_t sz = 0;
        for (typename bcData::const_iterator it = m_bc.begin(); it != m_bc.end(); ++it)
            sz += it->second.size();
        return sz + corner_values.size() + coupled_boundaries.size();
    }

    const bcContainer & container(const std::string & label) const {return m_bc[label]; }

    bcRefList get(const std::string & label, const short_t unk = 0, int comp = -1) const
    {
        bcRefList result;
        const const_bciterator it = m_bc.find(label);
        if ( it != m_bc.end() )
            for (const_iterator c = it->second.begin(); c!= it->second.end(); ++c)
                if ( c->m_unknown == unk )
                    if ( c->m_unkcomp == comp || comp == -1 )
                        result.push_back(*c);
        return result;
    }

    const bcContainer & dirichletSides() const {return m_bc["Dirichlet"]; }

    const bcContainer & weakDirichletSides() const {return m_bc["Weak Dirichlet"]; }

    const bcContainer & neumannSides()   const {return m_bc["Neumann"]; }

    const bcContainer & robinSides()     const {return m_bc["Robin"]; }

    const cornerContainer & cornerValues() const  {return corner_values;  }

    const cplContainer & coupledBoundaries() const {return coupled_boundaries; }

    bcContainer reducedContainer(const bcContainer & container, short_t unknown) const
    {
        bcContainer red;
        //red.reserve(container.size());
        for(typename bcContainer::const_iterator iter=container.begin(); iter!=container.end();++iter)
        {
            if(iter->unknown()==unknown)
                red.push_back(*iter);
        }
        return red;
    }

    bcContainer allConditions() const
    {
        bcContainer all;
        //all.reserve( size() - corner_values.size() );
        for (typename bcData::const_iterator it = m_bc.begin(); it != m_bc.end(); ++it)
            all.insert( all.end(), it->second.begin(), it->second.end() );
        return all;
    }

    const_iterator begin(const std::string & label) const {return m_bc[label].begin(); }

    iterator begin(const std::string & label) { return m_bc[label].begin(); }

    const_iterator end(const std::string & label) const {return m_bc[label].end(); }

    iterator end(const std::string & label) { return m_bc[label].end(); }

    const_bciterator beginAll() const {return m_bc.begin(); }
    bciterator beginAll() {return m_bc.begin(); }

    const_bciterator endAll() const {return m_bc.end(); }
    bciterator endAll() {return m_bc.end(); }

    const_iterator dirichletBegin() const
    { return m_bc["Dirichlet"].begin(); }

    const_iterator dirichletEnd() const
    { return m_bc["Dirichlet"].end(); }

    iterator dirichletBegin()
    { return m_bc["Dirichlet"].begin(); }

    iterator dirichletEnd()
    { return m_bc["Dirichlet"].end(); }

    const_iterator weakDirichletBegin() const
    { return m_bc["Weak Dirichlet"].begin(); }

    const_iterator weakDirichletEnd() const
    { return m_bc["Weak Dirichlet"].end(); }

    iterator weakDirichletBegin()
    { return m_bc["Weak Dirichlet"].begin(); }

    iterator weakDirichletEnd()
    { return m_bc["Weak Dirichlet"].end(); }


    const_iterator neumannBegin() const
    { return m_bc["Neumann"].begin(); }

    const_iterator neumannEnd() const
    { return m_bc["Neumann"].end(); }

    iterator neumannBegin()
    { return m_bc["Neumann"].begin(); }

    iterator neumannEnd()
    { return m_bc["Neumann"].end(); }

    const_iterator robinBegin() const
    { return m_bc["Robin"].begin(); }

    const_iterator robinEnd() const
    { return m_bc["Robin"].end(); }

    const_citerator cornerBegin() const
    { return corner_values.begin(); }

    const_citerator cornerEnd() const
    { return corner_values.end(); }

    const_cpliterator coupledBegin() const
    { return coupled_boundaries.begin(); }

    const_cpliterator coupledEnd() const
    { return coupled_boundaries.end(); }

    iterator robinBegin()
    { return m_bc["Robin"].begin(); }

    iterator robinEnd()
    { return m_bc["Robin"].end(); }

    citerator cornerBegin()
    { return corner_values.begin(); }

    citerator cornerEnd()
    { return corner_values.end(); }

    cpliterator coupledBegin()
    { return coupled_boundaries.begin(); }

    cpliterator coupledEnd()
    { return coupled_boundaries.end(); }


    void add(int p, boxSide s, const std::string & label,
             const function_ptr & f_ptr, short_t unknown = 0,
             int comp = -1, bool parametric = false)
    {
        m_bc[label].push_back(
                boundary_condition<T>(p, s, f_ptr, label, unknown, comp, parametric) );
    }

    void add(int p, boxSide s, const std::string & label,
             gsFunctionSet<T> * f, short_t unknown = 0,
             int comp = -1, bool parametric = false)
    {
        function_ptr f_ptr = memory::make_shared_not_owned(f);
        m_bc[label].push_back(
                boundary_condition<T>(p, s, f_ptr, label, unknown, comp, parametric) );
    }

    void add(int p, boxSide s, const std::string & label,
             const gsFunctionSet<T> & f, short_t unknown = 0,
             int comp = -1, bool parametric = false)
    {
        function_ptr fun = memory::make_shared(f.clone().release());
        add(p,s,label,fun,unknown,comp,parametric);
    }

    void addCondition(int p, boxSide s, condition_type::type t,
                      gsFunctionSet<T> * f, short_t unknown = 0, bool parametric = false, int comp = -1)
    {
        function_ptr fun = memory::make_shared_not_owned(f);
        addCondition(p,s,t,fun,unknown,parametric,comp);
    }

    void addConditions(const bcRefList & bcrf)
    {
        for (auto bc : bcrf)
            m_bc[bc.get().ctype()].push_back(bc);
    }

    void addCondition(int p, boxSide s, condition_type::type t,
                      const function_ptr & f_shptr, short_t unknown = 0,
                      bool parametric = false, int comp = -1)
    {
        auto bc = boundary_condition<T>(p,s,f_shptr,t,unknown,comp,parametric);
        // Every condition type can only be applied once
        const auto & c = m_bc[bc.ctype()];
        auto it = std::find_if(c.begin(), c.end(), [&bc](const boundary_condition<T> & b){ return bc.isSame(b); });
        if (it==c.end())
            m_bc[bc.ctype()].push_back(bc);
        else // bc of this type is already defined on patch p and side s
            gsWarn<<"Condition of type "<<bc.ctype()<<" on patch "<<bc.patch()<<" side "<<bc.side()<<" of unknown "<<bc.unknown()<<" with component "<<bc.unkComponent()<<" ignored, because it has already been defined\n";
    }

    void addCondition(int p, boxSide s, condition_type::type t,
                      const gsFunctionSet<T> & func, short_t unknown = 0,
                      bool parametric = false, int comp = -1)
    {
        function_ptr fun(func.clone().release());
        addCondition(p,s,t,fun,unknown,parametric,comp);
    }

    void addCondition( boxSide s, condition_type::type t,
                       gsFunctionSet<T> * f, short_t unknown = 0, bool parametric = false, int comp = -1)
    {
        // for single-patch only
        addCondition(0,s,t,f,unknown,parametric,comp);
    }

    void addCondition(const patchSide& ps, condition_type::type t,
                      gsFunctionSet<T> * f, short_t unknown = 0, bool parametric = false, int comp = -1)
    {
        addCondition(ps.patch, ps.side(), t, f, unknown,parametric,comp);
    }

    void addCondition(const patchSide& ps, condition_type::type t,
                      const function_ptr & f_shptr, short_t unknown = 0, bool parametric = false, int comp = -1)
    {
        addCondition(ps.patch, ps.side(), t, f_shptr, unknown,parametric,comp);
    }

    void addCondition(const patchSide& ps, condition_type::type t,
                      const gsFunctionSet<T> & func, short_t unknown = 0, bool parametric = false, int comp = -1)
    {
        addCondition(ps.patch, ps.side(), t, func, unknown,parametric,comp);
    }

    void addCondition(int p, boundary::side s, condition_type::type t,
                      const gsFunctionSet<T> & func, short_t unknown = 0,
                      bool parametric = false, int comp = -1)
    {
        function_ptr fun(func.clone().release());
        addCondition(p,boxSide(s),t,fun,unknown,parametric,comp);
    }

    void addCondition(int p, boundary::side s, condition_type::type t,
                      gsFunctionSet<T> * func, short_t unknown = 0,
                      bool parametric = false, int comp = -1)
    {
        addCondition(p,boxSide(s),t,func,unknown,parametric,comp);
    }

    void addCornerValue(boxCorner c, T value, int p = 0, short_t unknown = 0, int component = -1)
    {
        corner_values.push_back( corner_value<T>(p,c,value,unknown,component) );
    }

    void addCornerValue(boundary::corner c, T value, int p = 0, short_t unknown = 0, int component = -1)
    {
        corner_values.push_back( corner_value<T>(p,boxCorner(c),value,unknown,component) );
    }

    void addCoupled(int p1, boundary::side s1, int p2, boundary::side s2, index_t dim,
                      short_t unknown = 0, int comp = -1)
    {
        coupled_boundaries.push_back( coupled_boundary<T>(p1,boxSide(s1),p2,boxSide(s2),dim,unknown,comp));
    }

    void addCoupled(int p1, boxSide s1, int p2, boxSide s2, index_t dim,
                      short_t unknown = 0, int comp = -1)
    {
        coupled_boundaries.push_back( coupled_boundary<T>(p1,s1,p2,s2,dim,unknown,comp));
    }

    std::ostream &print(std::ostream &os, const bool verbose = false) const {
      // os << "gsBoundaryConditions :\n";
      for (typename bcData::const_iterator it = m_bc.begin(); it != m_bc.end();
           ++it)
        os << "* " << std::setw(13) << std::left << it->first << " : "
           << it->second.size() << "\n";

      if (!corner_values.empty())
        os << "* Corner values : " << corner_values.size() << "\n";

      // This block prints out all boundary conditions with more information
      if (verbose) {
        os << "*\n* Summary\n*\n* " << std::right << std::setw(15) << "Type"
           << std::setw(8) << "Patch" << std::setw(7) << "Side" << std::setw(9)
           << "Unknown" << std::setw(13) << "Components" << std::setw(11)
           << "Function" << std::endl;
        for (auto a = beginAll(); a != endAll(); a++) {
          for (const auto &element : a->second) {
            os << "* " << std::right << std::setw(15) << element.ctype()
               << std::setw(8) << element.patch() << std::setw(7)
               << element.side().index() << std::setw(9) << element.unknown();

            const auto &component = element.unkComponent();
            if (component == -1) {
              os << std::setw(13) << "all";
            } else {
              os << std::setw(13) << (element.unkComponent());
            }
            if (element.function()) {
              os << "   " << *element.function();
            } else {
              os << "   Homogeneous";
            }
            os << "\n";
          }
        }
        os << "*" << std::endl;
      }

      return os;
    }

    const boundary_condition<T>* getConditionFromSide (patchSide ps) const
    {
        const_iterator beg, end, cur;
        patchSideComparison psRef(ps);
        for (auto & bcc : m_bc)
        {
            beg = bcc.second.begin();
            end = bcc.second.end();
            cur = std::find_if(beg,end,psRef);
            if (cur != end)
                return &(*cur);
        }
        return nullptr;
    }

    void getConditionFromSide (patchSide ps, bcContainer& result) const
    {
        result.clear();
        const_iterator beg, end, cur;
        for (auto & bcc : m_bc)
        {
            beg = bcc.second.begin();
            end = bcc.second.end();
            for(cur=beg; cur!=end; ++cur)
            if (cur->ps == ps)
                result.push_back(*cur);
        }
    }

    void getConditionsForPatch(const index_t np, gsBoundaryConditions& result) const
    {
        result.clear();
        bcContainer bc_all = allConditions(); //inefficient, but fewer code
        for(const_iterator it = bc_all.begin(); it!= bc_all.end();it++)
        {
            if((*it).patch()==np)
                result.add(0,(*it).side(),it->ctype(),(*it).function(),(*it).unknown());
        }

        for(const_citerator it = cornerBegin(); it!= cornerEnd();it++)
        {
            if((*it).patch==np)
                result.addCornerValue( (*it).corner, (*it).value,  0, (*it).unknown);
        }
    }

    // Periodic conditions

    size_t numPeriodic() const { return m_periodicPairs.size(); }

    const ppContainer & periodicPairs() const {return m_periodicPairs; }

    const_ppiterator periodicBegin() const
    { return m_periodicPairs.begin(); }

    const_ppiterator periodicEnd() const
    { return m_periodicPairs.end(); }

    ppiterator periodicBegin()
    { return m_periodicPairs.begin(); }

    ppiterator periodicEnd()
    { return m_periodicPairs.end(); }

    void addPeriodic(int p1, boxSide s1, int p2, boxSide s2, short_t dim)
    { m_periodicPairs.push_back( boundaryInterface(patchSide(p1,s1), patchSide(p2,s2), dim) ); }

    void clearPeriodicPairs() { m_periodicPairs.clear(); }

    void setTransformMatrix(gsMatrix<T> trMatrix)
    { m_trMatrix = trMatrix; }

    void setIdentityMatrix(short_t dim)
    { m_trMatrix = gsMatrix<T>::Identity(dim, dim); }

    gsMatrix<T> getTransformMatrix() const
    {
        GISMO_ASSERT(m_trMatrix.rows() > 0, "Transformation matrix for periodic conditions not set!");
        return m_trMatrix;
    }

    void setGeoMap(const gsFunctionSet<T> & gm)
    {
        m_patches = &gm;
    }

    bool hasGeoMap() const
    {
        return nullptr!=m_patches;
    }

    const gsFunctionSet<T> & geoMap() const
    {
        GISMO_ASSERT(nullptr!=m_patches, "Geometry map was not provided in BC.");
        return *m_patches;
    }

private: // Data members
    struct patchSideComparison
    {
        const patchSide m_ps;
        patchSideComparison(patchSide ps)
                : m_ps(ps)
        {}

        bool operator() (const boundary_condition<T> &bc) const
        {
            return bc.ps==m_ps;
        }
    };

    cornerContainer corner_values; 

    cplContainer coupled_boundaries; 

    mutable bcData m_bc;  

    ppContainer m_periodicPairs; // TODO: add read from xml
    gsMatrix<T> m_trMatrix;

    // Pointer to associated multipatch domain
    const gsFunctionSet<T> * m_patches;

}; // class gsBoundaryConditions

template<class T>
std::ostream &operator<<(std::ostream &os, const gsBoundaryConditions<T>& bvp)
{return bvp.print(os); }

#ifdef GISMO_WITH_PYBIND11

    void pybind11_init_gsBoundaryConditions(pybind11::module &m);
    void pybind11_enum_gsBoundaryConditions(pybind11::module &m);

#endif

} // namespace gismo

#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsBoundaryConditions.hpp)
#endif