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