gsGenericAssembler.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsAssembler/gsAssembler.h>
#include <gsPde/gsLaplacePde.h>
 
#include <gsAssembler/gsExpressions.h>
#include <gsAssembler/gsExprHelper.h>
#include <gsAssembler/gsExprAssembler.h>
 
namespace gismo
{
 
template <class T>
class gsGenericAssembler : public gsAssembler<T>
{
public:
    typedef gsAssembler<T> Base;
 
public:
 
    static gsOptionList defaultOptions();
 
    gsGenericAssembler( const gsMultiPatch<T>         & patches,
                        const gsMultiBasis<T>         & bases,
                        const gsOptionList            & opt = defaultOptions(),
                        const gsBoundaryConditions<T> * bc = NULL)
    : m_pde(patches)
    {
        if ( bc != NULL)
        {
            m_pde.boundaryConditions() = *bc;
            this->m_ddof.resize(1);
        }
 
        Base::initialize(m_pde, bases, opt);
        gsGenericAssembler::refresh();
    }
 
    /* TODO (id geometry)
    // Constructor with gsMultiBasis
    gsGenericAssembler( gsMultiBasis<T> const         & bases,
                        bool conforming = false,
                        const gsBoundaryConditions<T> * bc = NULL)
    : Base(patches)
    {
        m_bases.push_back(bases);
 
        Base::initialize(m_pde, bases, opt);
        gsGenericAssembler::refresh();
    }
    */
 
    void refresh()
    {
        // Setup sparse system
        gsDofMapper mapper;
        m_bases[0].getMapper(
            (dirichlet::strategy)(m_options.getInt("DirichletStrategy")),
            (iFace::strategy)(m_options.getInt("InterfaceStrategy")),
            this->pde().bc(), mapper, 0);
        m_system = gsSparseSystem<T>(mapper);
        //note: no allocation here
        //        const index_t nz = m_options.numColNz(m_bases[0][0]);
        //        m_system.reserve(nz, 1);
    }
 
    void refresh(gsDofMapper mapper)
    {
         m_system = gsSparseSystem<T>(mapper);
    }
 
    const gsSparseMatrix<T> & assembleMass(const index_t patchIndex = -1, bool refresh = true);
 
    const gsSparseMatrix<T> & assembleMass2()
    {
        typedef typename gsExprAssembler<T>::geometryMap geometryMap;
        typedef typename gsExprAssembler<T>::space       space;
 
        // Elements used for numerical integration
        gsExprAssembler<T> A(1,1);
        A.setIntegrationElements(m_bases.front());
        geometryMap G = A.getMap(m_pde_ptr->patches());
        space u = A.getSpace(m_bases.front());
 
        A.initSystem();
        //m_system.matrix() = A.assemble( u * u.tr() * meas(G) );
        A.assemble( u * u.tr() * meas(G) );
        m_system.matrix() = A.matrix();
 
        return m_system.matrix();
    }
 
    const gsSparseMatrix<T> & assembleStiffness(const index_t patchIndex = -1, const bool refresh = true);
 
    const gsMatrix<T> & assembleMoments(const gsFunction<T> & func, index_t patchIndex = -1, bool refresh = true);
 
    const gsSparseMatrix<T> & assembleDG(const boundaryInterface & iFace, bool refresh = true);
 
    const gsSparseMatrix<T> & assembleNeumann(const boundary_condition<T> & bc, bool refresh = true);
 
    const gsSparseMatrix<T> & assembleNitsche(const boundary_condition<T> & bc, bool refresh = true);
 
    typename gsSparseMatrix<T>::fullView fullMatrix()
    {
        return m_system.matrix().template selfadjointView<Lower>();
    }
 
    const typename gsSparseMatrix<T>::constFullView fullMatrix() const
    {
        return m_system.matrix().template selfadjointView<Lower>();
    }
 
 
private:
 
    // Members from gsAssembler
    using Base::m_pde_ptr;
    using Base::m_bases;
    using Base::m_ddof;
    using Base::m_options;
    using Base::m_system;
 
private:
 
    gsLaplacePde<T> m_pde;
};
 
 
 
} // namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsGenericAssembler.hpp)
#endif