dudv_functional.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/base_header.hpp>
#include <kernel/assembly/bilinear_operator_assembler.hpp> // for BilinearOperatorAssembler
#include <kernel/assembly/common_operators.hpp>            // for DuDvOperator
#include <kernel/assembly/symbolic_assembler.hpp>          // for SymbolicMatrixAssembler
#include <kernel/assembly/slip_filter_assembler.hpp>
#include <kernel/assembly/unit_filter_assembler.hpp>
#include <kernel/cubature/dynamic_factory.hpp>             // for DynamicFactory
#include <kernel/geometry/conformal_mesh.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/lafem/filter_chain.hpp>
#include <kernel/lafem/filter_sequence.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/slip_filter.hpp>
#include <kernel/lafem/unit_filter_blocked.hpp>
#include <kernel/meshopt/mesh_quality_functional.hpp>
#include <kernel/space/lagrange1/element.hpp>
 
namespace FEAT
{
  namespace Meshopt
  {
    template
    <
      typename DT_, typename IT_, typename TrafoType_,
      template<typename, typename, int, int> class MatrixType_ = LAFEM::SparseMatrixBCSR
    >
    class DuDvFunctional:
      public MeshQualityFunctional<typename TrafoType_::MeshType>
    {
      public:
        typedef DT_ DataType;
        typedef IT_ IndexType;
        typedef TrafoType_ TrafoType;
        typedef typename TrafoType::MeshType MeshType;
        typedef typename MeshType::ShapeType ShapeType;
 
        typedef MatrixType_<DT_, IT_, MeshType::world_dim, MeshType::world_dim> MatrixType;
        static constexpr int BlockHeight = MatrixType::BlockHeight;
        static constexpr int BlockWidth = MatrixType::BlockWidth;
 
        template <typename DT2_ = DT_, typename IT2_ = IT_>
        using ContainerType = DuDvFunctional<DT2_, IT2_, TrafoType_, MatrixType_>;
 
        template <typename DataType2_, typename IndexType2_>
        using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
        typedef MeshQualityFunctional<MeshType> BaseClass;
        typedef typename BaseClass::CoordsBufferType CoordsBufferType;
 
        typedef typename MeshType::CoordType CoordType;
 
        //template<typename A, typename B, typename C>
        //using MatrixTemplate = MatrixType_<A, B, C, MeshType::world_dim, MeshType::world_dim>;
 
        typedef typename MatrixType::VectorTypeL VectorTypeL;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef LAFEM::DenseVector<DT_, IT_> ScalarVectorType;
 
        typedef LAFEM::UnitFilterBlocked<DT_, IT_, MeshType::world_dim> DirichletFilterType;
        typedef LAFEM::FilterSequence<DirichletFilterType> DirichletFilterSequence;
        typedef LAFEM::SlipFilter<DT_, IT_, MeshType::world_dim> SlipFilterType;
        typedef LAFEM::FilterSequence<SlipFilterType> SlipFilterSequence;
        typedef LAFEM::FilterChain<SlipFilterSequence, DirichletFilterSequence> FilterType;
 
        typedef typename Intern::TrafoFE<TrafoType>::Space SpaceType;
        // 2 * (degree of trafo) for both trial and test spaces, +1 for safety reasons
        // This could be decreased by the degree of the operator, i.e. 2 for Du:Dv
        static constexpr int _local_degree = 4*SpaceType::local_degree + 1;
 
        MatrixType matrix_sys;
 
      protected:
        TrafoType& _trafo;
        ScalarVectorType _lambda;
        std::map<String, std::shared_ptr<Assembly::UnitFilterAssembler<MeshType>>> _dirichlet_asm;
        std::map<String, std::shared_ptr<Assembly::SlipFilterAssembler<TrafoType>>> _slip_asm;
        Cubature::DynamicFactory _cubature_factory;
 
      public:
        SpaceType trafo_space;
 
      public:
        explicit DuDvFunctional(
          Geometry::RootMeshNode<MeshType>* rmn_,
          TrafoType& trafo,
          const std::deque<String>& dirichlet_list, const std::deque<String>& slip_list):
          BaseClass(rmn_),
          matrix_sys(),
          _trafo(trafo),
          _lambda(rmn_->get_mesh()->get_num_entities(MeshType::shape_dim)),
          _dirichlet_asm(),
          _slip_asm(),
          _cubature_factory("auto-degree:"+stringify(int(_local_degree))),
          trafo_space(trafo)
          {
            // The symbolic assembly has to happen in the constructor because the matrix is shallow copied immediately
            // after construction in the Control::QuadraticSystemLevel
            Assembly::SymbolicAssembler::assemble_matrix_std1(matrix_sys, trafo_space);
 
            // For every MeshPart specified in the list of Dirichlet boundaries, create a UnitFilterAssembler and
            // insert it into the std::map
            for(auto& it : dirichlet_list)
            {
              // Create empty assembler
              auto new_asm = std::make_shared<Assembly::UnitFilterAssembler<MeshType>>();
 
              // Add the MeshPart to the assembler if it is there. There are legimate reasons for it NOT to be
              // there, i.e. we are in parallel and our patch is not adjacent to that MeshPart
              auto* mpp = rmn_->find_mesh_part(it);
              if(mpp != nullptr)
              {
                new_asm->add_mesh_part(*mpp);
              }
 
              // Insert into the map
              String identifier(it);
              _dirichlet_asm.emplace(identifier, new_asm);
            }
 
            // For every MeshPart specified in the list of slip boundaries, create a SlipFilterAssembler and
            // insert it into the std::map
            for(auto& it : slip_list)
            {
              // Create empty assembler
              auto new_asm = std::make_shared<Assembly::SlipFilterAssembler<TrafoType>>(this->_trafo);
 
              // Add the MeshPart to the assembler if it is there. There are legimate reasons for it NOT to be
              // there, i.e. we are in parallel and our patch is not adjacent to that MeshPart
              auto* mpp = rmn_->find_mesh_part(it);
              if(mpp != nullptr)
              {
                new_asm->add_mesh_part(*mpp);
              }
 
              // Insert into the map
              String identifier(it);
              _slip_asm.emplace(identifier, new_asm);
            }
 
          }
 
        DuDvFunctional(const DuDvFunctional&) = delete;
 
        virtual ~DuDvFunctional()
        {
        }
 
        virtual void init() override
        {
          //Assembly::SymbolicAssembler::assemble_matrix_std1(matrix_sys, trafo_space);
          assemble_system_matrix();
        }
 
        static String name()
        {
          return "DuDvFunctional<"+MeshType::name()+">";
        }
 
        virtual void add_to_vtk_exporter(Geometry::ExportVTK<MeshType>& exporter) const override
        {
          exporter.add_cell_scalar("lambda", this->_lambda.elements());
        }
 
        virtual void assemble_system_matrix()
        {
          matrix_sys.format();
 
          Assembly::Common::DuDvOperatorBlocked<MeshType::world_dim> my_operator;
 
          Assembly::BilinearOperatorAssembler::assemble_matrix1(
            matrix_sys,           // the matrix that receives the assembled operator
            my_operator, // the operator that is to be assembled
            trafo_space,            // the finite element space in use
            _cubature_factory  // the cubature factory to be used for integration
            );
 
        }
 
        virtual void prepare(VectorTypeR& vec_state, FilterType& filter)
        {
          for(Index cell(0); cell < this->get_mesh()->get_num_entities(MeshType::shape_dim); ++cell)
          {
            _lambda(cell, DataType(_trafo.template compute_vol<typename MeshType::ShapeType>(cell)));
          }
 
          auto& dirichlet_filters = filter.template at<1>();
 
          for(auto& it : dirichlet_filters)
          {
            const auto& assembler = _dirichlet_asm.find(it.first);
            if(assembler == _dirichlet_asm.end())
            {
              XABORTM("Could not find dirichlet assembler for filter with key "+it.first);
            }
 
            assembler->second->assemble(it.second, trafo_space, vec_state);
          }
 
          // The slip filter contains the outer unit normal, so reassemble it
          auto& slip_filters = filter.template at<0>();
 
          for(auto& it : slip_filters)
          {
            const auto& assembler = _slip_asm.find(it.first);
            if(assembler == _slip_asm.end())
            {
              XABORTM("Could not find slip filter assembler for filter with key "+it.first);
            }
 
            assembler->second->assemble(it.second, trafo_space);
          }
 
          for(const auto& it:slip_filters)
          {
            this->_mesh_node->adapt_by_name(it.first);
          }
 
        }
 
        virtual CoordType compute_cell_size_defect(CoordType& lambda_min, CoordType& lambda_max, CoordType& vol_min, CoordType& vol_max, CoordType& vol) const override
        {
 
          compute_cell_size_defect_pre_sync(vol_min, vol_max, vol);
          return compute_cell_size_defect_post_sync(lambda_min, lambda_max, vol_min, vol_max, vol);
 
        }
 
        void compute_cell_size_defect_pre_sync(CoordType& vol_min, CoordType& vol_max, CoordType& vol) const
        {
          vol = CoordType(0);
 
          vol_min = Math::huge<CoordType>();
          vol_max = CoordType(0);
 
          for(Index cell(0); cell < this->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
          {
            CoordType my_vol = this->_trafo.template compute_vol<ShapeType, CoordType>(cell);
            vol_min = Math::min(vol_min, my_vol);
            vol_max = Math::max(vol_min, my_vol);
            vol += my_vol;
          }
        }
 
        virtual CoordType compute_cell_size_defect_post_sync(CoordType& lambda_min, CoordType& lambda_max, CoordType& vol_min, CoordType& vol_max, const CoordType& vol) const
        {
          CoordType size_defect(0);
          lambda_min = Math::huge<CoordType>();
          lambda_max = CoordType(0);
 
          for(Index cell(0); cell < this->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
          {
            size_defect += Math::abs(this->_trafo.template compute_vol<ShapeType, CoordType>(cell)/vol - this->_lambda(cell));
            lambda_min = Math::min(lambda_min, CoordType(this->_lambda(cell)));
            lambda_max = Math::max(lambda_max, CoordType(this->_lambda(cell)));
          }
 
          vol_min /= vol;
          vol_max/= vol;
 
          return size_defect;
        }
 
        bool empty() const
        {
          return matrix_sys.empty();
        }
 
        VectorTypeL create_vector_l() const
        {
          return VectorTypeL(trafo_space.get_num_dofs());
        }
 
        VectorTypeR create_vector_r() const
        {
          return VectorTypeR(trafo_space.get_num_dofs());
        }
 
 
    }; // class DuDvFunctional
 
#ifdef FEAT_EICKT
    extern template class DuDvFunctional
    <
      double, Index,
      Trafo::Standard::Mapping<Geometry::ConformalMesh<Shape::Simplex<2>, 2, double>>,
      LAFEM::SparseMatrixBCSR
    >;
 
    extern template class DuDvFunctional
    <
      double, Index,
      Trafo::Standard::Mapping<Geometry::ConformalMesh<Shape::Hypercube<2>, 2, double>>,
      LAFEM::SparseMatrixBCSR
    >;
#endif // FEAT_EICKT
  } // namespace Meshopt
} // namespace FEAT