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