meshopt_control.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/grid_transfer.hpp>
#include <kernel/assembly/mirror_assembler.hpp>
#include <kernel/assembly/unit_filter_assembler.hpp>
#include <kernel/assembly/slip_filter_assembler.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/sparse_matrix_bwrappedcsr.hpp>
#include <kernel/lafem/transfer.hpp>
#include <kernel/lafem/vector_mirror.hpp>
#include <kernel/geometry/mesh_file_writer.hpp>
#include <kernel/global/filter.hpp>
#include <kernel/global/gate.hpp>
#include <kernel/global/matrix.hpp>
#include <kernel/global/muxer.hpp>
#include <kernel/global/nonlinear_functional.hpp>
#include <kernel/global/transfer.hpp>
#include <kernel/global/vector.hpp>
 
#include <control/domain/domain_control.hpp>
 
namespace FEAT
{
  namespace Control
  {
    namespace Meshopt
    {
      template<typename, typename, template<typename, typename> class>
      class MeshoptSystemLevel;
 
      template<typename, typename, template<typename, typename> class>
      class QuadraticSystemLevel;
 
      template<typename, typename, template<typename, typename> class>
      class NonlinearSystemLevel;
 
 
      template<typename DomainControl_>
      class MeshoptControlBase
      {
        public:
          typedef DomainControl_ DomainControl;
          typedef typename DomainControl_::LevelType DomainLevel;
          typedef typename DomainControl_::MeshType MeshType;
 
          static constexpr int world_dim = DomainControl_::MeshType::world_dim;
          typedef typename DomainControl_::MeshType::CoordType CoordType;
 
          typedef LAFEM::DenseVectorBlocked<CoordType, Index, world_dim> LocalCoordsBuffer;
          typedef LAFEM::VectorMirror<CoordType, Index> CoordsMirror;
          typedef Global::Vector<LocalCoordsBuffer, CoordsMirror> GlobalCoordsBuffer;
 
          typedef Geometry::ExportVTK<MeshType> VTKExporterType;
 
        protected:
          DomainControl& _dom_ctrl;
          std::deque<String> _dirichlet_boundaries;
          std::deque<String> _slip_boundaries;
 
        public:
 
          explicit MeshoptControlBase(DomainControl& dom_ctrl_, const std::deque<String>& dirichlet_boundaries,
          const std::deque<String>& slip_boundaries)
            : _dom_ctrl(dom_ctrl_),
            _dirichlet_boundaries(dirichlet_boundaries),
            _slip_boundaries(slip_boundaries)
            {
            }
 
          virtual ~MeshoptControlBase()
          {
          }
 
          // * \brief Gets a const reference to the mesh at a certain level
          // */
          //const MeshType& get_mesh(Index lvl)
          //{
          //  return _dom_ctrl.at(lvl)->get_mesh();
          //}
 
          virtual CoordType compute_cell_size_defect(CoordType& lambda_min, CoordType& lambda_max,
          CoordType& vol_min, CoordType& vol_max, CoordType& vol) const = 0;
 
          virtual void optimize() = 0;
 
          virtual void prepare(const GlobalCoordsBuffer& DOXY(vec_state)) = 0;
 
          virtual void mesh_to_buffer() = 0;
 
          virtual void buffer_to_mesh() = 0;
 
          virtual GlobalCoordsBuffer& get_coords() = 0;
 
          const std::deque<String>& get_dirichlet_boundaries() const
          {
            return _dirichlet_boundaries;
          }
 
          const std::deque<String>& get_slip_boundaries() const
          {
            return _slip_boundaries;
          }
 
          virtual size_t get_num_levels() const = 0;
 
          virtual String name() const = 0;
 
          virtual void add_to_vtk_exporter
            (Geometry::ExportVTK<MeshType>& DOXY(exporter), const int DOXY(lvl_index)) const
            {
            }
 
          virtual String info() const = 0;
 
 
        void compute_mesh_quality(CoordType& edge_angle, CoordType& qi_min, CoordType& qi_mean,
          CoordType* edge_angle_cellwise = nullptr, CoordType* qi_cellwise = nullptr,
          int lvl_index = -1) const
        {
          // max_level_index cannot be called for the default argument, so we do it here
          if(lvl_index == -1)
          {
            lvl_index = _dom_ctrl.max_level_index();
          }
          XASSERT(lvl_index >= _dom_ctrl.min_level_index());
          XASSERT(lvl_index <= _dom_ctrl.max_level_index());
 
          const Dist::Comm& comm = _dom_ctrl.comm();
 
          CoordType qi_sum(0);
 
          for(std::size_t i(0u); i < _dom_ctrl.size_physical(); ++i)
          {
            const auto& dom_lvl = *_dom_ctrl.at(i);
            if(dom_lvl.get_level_index() == lvl_index)
            {
              const auto& my_mesh = dom_lvl.get_mesh();
 
              Index ncells(my_mesh.get_num_entities(MeshType::shape_dim));
              comm.allreduce(&ncells, &ncells, std::size_t(1), Dist::op_sum);
 
              Geometry::MeshQualityHeuristic<typename MeshType::ShapeType>::compute(qi_min, qi_sum,
              my_mesh.template get_index_set<MeshType::shape_dim, 0>(), my_mesh.get_vertex_set(), qi_cellwise);
 
              comm.allreduce(&qi_min, &qi_min, std::size_t(1), Dist::op_min);
              comm.allreduce(&qi_sum, &qi_sum, std::size_t(1), Dist::op_sum);
 
              edge_angle = Geometry::MeshQualityHeuristic<typename MeshType::ShapeType>::angle(
              my_mesh.template get_index_set<MeshType::shape_dim, 0>(), my_mesh.get_vertex_set(), edge_angle_cellwise);
              comm.allreduce(&edge_angle, &edge_angle, std::size_t(1), Dist::op_min);
 
              qi_mean = qi_sum/CoordType(ncells);
 
              return;
            }
          }
 
          // We should never get to this point
          XABORTM("Could not find level with index "+stringify(lvl_index)+"!\n");
 
        }
      }; // class MeshoptControlBase
 
 
      template
      <
        typename DT_, typename IT_,
        template<typename, typename> class Functional_
      >
      class MeshoptSystemLevel
      {
        public:
          typedef DT_ DataType;
          typedef IT_ IndexType;
 
          typedef Functional_<DT_, IT_> LocalFunctional;
 
          typedef typename LocalFunctional::SpaceType SpaceType;
          typedef typename LocalFunctional::TrafoType TrafoType;
          typedef typename TrafoType::MeshType MeshType;
 
          typedef typename LocalFunctional::VectorTypeL LocalSystemVectorL;
          typedef typename LocalFunctional::VectorTypeR LocalSystemVectorR;
          typedef typename LocalFunctional::ScalarVectorType LocalScalarVector;
          typedef typename LocalFunctional::CoordsBufferType LocalCoordsBuffer;
          typedef LAFEM::SparseMatrixBWrappedCSR<DT_, IT_, LocalFunctional::MeshType::world_dim>
            LocalSystemTransferMatrix;
          typedef LAFEM::Transfer<LocalSystemTransferMatrix> LocalSystemTransfer;
 
          typedef typename LocalFunctional::FilterType LocalSystemFilter;
          typedef typename LocalFunctional::SlipFilterSequence LocalSlipFilterSequence;
          typedef typename LocalFunctional::DirichletFilterSequence LocalDirichletFilterSequence;
 
          //typedef LAFEM::VectorMirrorBlocked<DT_, IT_, LocalFunctional::BlockHeight> SystemMirror;
          typedef LAFEM::VectorMirror<DT_, IT_> SystemMirror;
          typedef Global::Gate<LocalSystemVectorR, SystemMirror> SystemGate;
          typedef LAFEM::VectorMirror<DT_, IT_> ScalarMirror;
          typedef Global::Gate<LocalScalarVector, ScalarMirror> ScalarGate;
 
          typedef Global::Filter<LocalSystemFilter, SystemMirror> GlobalSystemFilter;
          typedef Global::Vector<LocalScalarVector, ScalarMirror> GlobalScalarVector;
          typedef Global::Vector<LocalSystemVectorL, SystemMirror> GlobalSystemVectorL;
          typedef Global::Vector<LocalSystemVectorR, SystemMirror> GlobalSystemVectorR;
          typedef Global::Vector<LocalCoordsBuffer, SystemMirror> GlobalCoordsBuffer;
          typedef Global::Transfer<LocalSystemTransfer, SystemMirror> GlobalSystemTransfer;
          typedef Global::Muxer<LocalSystemVectorR, SystemMirror> GlobalSystemMuxer;
 
        public:
          LocalFunctional local_functional;
          ScalarGate gate_scalar;
          SystemGate gate_sys;
          GlobalSystemFilter filter_sys;
          GlobalCoordsBuffer coords_buffer;
          GlobalSystemMuxer coarse_muxer_sys;
          GlobalSystemTransfer transfer_sys;
 
        private:
          const int _level_index;
 
        public:
          template<typename... Args_>
          explicit MeshoptSystemLevel(const int level_index,
          Geometry::RootMeshNode<MeshType>* rmn,
          TrafoType& trafo,
          const std::deque<String>& dirichlet_list,
          const std::deque<String>& slip_list,
          Args_&&... args) :
            local_functional(rmn, trafo, dirichlet_list, slip_list, std::forward<Args_>(args)...),
            gate_sys(),
            filter_sys(),
            coords_buffer(&gate_sys, local_functional.get_coords().clone(LAFEM::CloneMode::Shallow)),
            coarse_muxer_sys(),
            transfer_sys(&coarse_muxer_sys),
            _level_index(level_index)
            {
 
              LocalSlipFilterSequence slip_sequence(slip_list);
              LocalDirichletFilterSequence dirichlet_sequence(dirichlet_list);
 
              LocalSystemFilter local_filter (std::move(slip_sequence), std::move(dirichlet_sequence));
 
              filter_sys.local() = std::move(local_filter);
 
            }
 
          MeshoptSystemLevel() = delete;
          MeshoptSystemLevel(const MeshoptSystemLevel&) = delete;
          MeshoptSystemLevel(MeshoptSystemLevel&&) = delete;
 
          virtual ~MeshoptSystemLevel()
          {
          }
 
          int get_level_index() const
          {
            return _level_index;
          }
 
          bool empty() const
          {
            return local_functional.empty();
          }
 
          GlobalSystemVectorL create_vector_l() const
          {
            return GlobalSystemVectorL(&gate_sys, local_functional.create_vector_l());
          }
 
          GlobalSystemVectorR create_vector_r() const
          {
            return GlobalSystemVectorR(&gate_sys, local_functional.create_vector_r());
          }
 
 
          void sync_system_filter()
          {
            // Sync the filter vectors in the SlipFilters
            auto& slip_filters = filter_sys.local().template at<0>();
            {
              for(auto& it : slip_filters)
              {
                // Get the filter vector
                auto& slip_filter_vector = it.second.get_filter_vector();
 
                if(slip_filter_vector.used_elements() > 0)
                {
                  // Temporary DenseVector for syncing
                  LocalSystemVectorL tmp(slip_filter_vector.size());
                  auto* tmp_elements = tmp.template elements<LAFEM::Perspective::native>();
                  auto* sfv_elements = slip_filter_vector.template elements<LAFEM::Perspective::native>();
 
                  // Copy sparse filter vector contents to DenseVector
                  for(Index isparse(0); isparse < slip_filter_vector.used_elements(); ++isparse)
                  {
                    Index idense(slip_filter_vector.indices()[isparse]);
                    tmp_elements[idense] = sfv_elements[isparse];
                  }
 
                  gate_sys.sync_0(tmp);
 
                  // Copy sparse filter vector contents to DenseVector
                  for(Index isparse(0); isparse < slip_filter_vector.used_elements(); ++isparse)
                  {
                    Index idense(slip_filter_vector.indices()[isparse]);
                    sfv_elements[isparse] = tmp_elements[idense];
                  }
                }
                else
                {
                  // Temporary DenseVector for syncing
                  LocalSystemVectorL tmp(slip_filter_vector.size());
                  gate_sys.sync_0(tmp);
                }
              }
 
            } // evil slip filter sync
 
          }
 
          GlobalSystemVectorL assemble_rhs_vector() const
          {
            XASSERTM(!empty(), "Assemble_rhs_vector for empty functional");
 
            // create new vector
            GlobalSystemVectorL vec_rhs(create_vector_l());
 
            vec_rhs.format();
            filter_sys.filter_rhs(vec_rhs);
 
            return vec_rhs;
          }
 
          GlobalSystemVectorR assemble_sol_vector() const
          {
            XASSERTM(!empty(), "Assemble_sol_vector for empty functional");
            GlobalSystemVectorR vec_sol(create_vector_r());
 
            vec_sol.local().copy(coords_buffer.local());
            filter_sys.local().template at<1>().filter_sol(vec_sol.local());
 
            return vec_sol;
          }
 
          template<typename DomainLayer_>
          void assemble_gates(const DomainLayer_& dom_layer)
          {
            // set the gate comm
            gate_sys.set_comm(dom_layer.comm_ptr());
            gate_scalar.set_comm(dom_layer.comm_ptr());
 
            // Loop over all ranks
            for(Index i(0); i < dom_layer.neighbor_count(); ++i)
            {
              int rank(dom_layer.neighbor_rank(i));
 
              // try to find our halo
              auto* halo = local_functional.get_mesh_node()->get_halo(rank);
              XASSERTM(halo != nullptr, "Halo not found.");
 
              // assemble the system mirror
              SystemMirror sys_mirror;
              Assembly::MirrorAssembler::assemble_mirror(sys_mirror, local_functional.trafo_space, *halo);
 
              // push mirror into gates
              gate_sys.push(rank, std::move(sys_mirror));
 
              // assemble the scalar mirror
              ScalarMirror scalar_mirror;
              Assembly::MirrorAssembler::assemble_mirror(scalar_mirror, local_functional.trafo_space, *halo);
 
              // push mirror into gates
              gate_scalar.push(rank, std::move(scalar_mirror));
            }
 
            // create local template vectors
            LocalSystemVectorR tmp_sys(local_functional.trafo_space.get_num_dofs());
            LocalScalarVector tmp_scalar(local_functional.trafo_space.get_num_dofs());
 
            // compile gates
            gate_sys.compile(std::move(tmp_sys));
            gate_scalar.compile(std::move(tmp_scalar));
          }
 
          void assemble_system_transfer(const MeshoptSystemLevel& level_coarse)
          {
            // get global transfer operator
            GlobalSystemTransfer& glob_trans = this->transfer_sys;
 
            // get local transfer operator
            LocalSystemTransfer& loc_trans = glob_trans.local();
 
            // get local transfer matrices
            typename LocalSystemTransferMatrix::BaseClass& loc_prol = loc_trans.get_mat_prol();
            typename LocalSystemTransferMatrix::BaseClass& loc_rest = loc_trans.get_mat_rest();
 
            // assemble structure?
            if (loc_prol.empty())
            {
              Assembly::SymbolicAssembler::assemble_matrix_2lvl(
                loc_prol, local_functional.trafo_space, level_coarse.local_functional.trafo_space);
            }
 
            // Create a global scalar weight vector
            GlobalScalarVector glob_vec_weight( &this->gate_scalar, loc_prol.create_vector_l());
            // Get local scalar weight vector
            auto& loc_vec_weight = glob_vec_weight.local();
 
            // Assemble the underlying scalar prolongation matrix
            {
              Cubature::DynamicFactory cubature_factory("auto-degree:3");
 
              loc_prol.format();
              loc_vec_weight.format();
 
              // Assemble prolongation matrix
              Assembly::GridTransfer::assemble_prolongation(loc_prol, loc_vec_weight,
              local_functional.trafo_space, level_coarse.local_functional.trafo_space, cubature_factory);
 
              // Synchronize weight vector
              glob_vec_weight.sync_0();
 
              // Invert components
              loc_vec_weight.component_invert(loc_vec_weight);
 
              // Scale prolongation matrix
              loc_prol.scale_rows(loc_prol, loc_vec_weight);
 
              // Copy and transpose
              loc_rest = loc_prol.transpose();
            }
          }
 
      }; // class MeshoptSystemLevel
 
      template
      <
        typename DT_, typename IT_,
        template<typename, typename> class Functional_
      >
      class QuadraticSystemLevel : public MeshoptSystemLevel<DT_, IT_, Functional_>
      {
        public:
          typedef DT_ DataType;
          typedef IT_ IndexType;
 
          typedef MeshoptSystemLevel<DT_, IT_, Functional_> BaseClass;
 
          typedef Functional_<DT_, IT_> LocalFunctional;
 
          typedef typename LocalFunctional::SpaceType SpaceType;
          typedef typename LocalFunctional::TrafoType TrafoType;
          typedef typename TrafoType::MeshType MeshType;
 
          typedef typename LocalFunctional::MatrixType LocalMatrix;
          typedef typename BaseClass::LocalSystemVectorL LocalSystemVectorL;
          typedef typename BaseClass::LocalSystemVectorR LocalSystemVectorR;
          typedef typename BaseClass::LocalScalarVector LocalScalarVector;
          typedef typename BaseClass::LocalCoordsBuffer LocalCoordsBuffer;
          typedef LAFEM::SparseMatrixBWrappedCSR<DT_, IT_, LocalFunctional::MeshType::world_dim>
            LocalSystemTransferMatrix;
          typedef LAFEM::Transfer<LocalSystemTransferMatrix> LocalSystemTransfer;
 
          typedef typename LocalFunctional::FilterType LocalSystemFilter;
          typedef typename LocalFunctional::SlipFilterSequence LocalSlipFilterSequence;
          typedef typename LocalFunctional::DirichletFilterSequence LocalDirichletFilterSequence;
 
          typedef typename BaseClass::SystemMirror SystemMirror;
          typedef typename BaseClass::SystemGate SystemGate;
          typedef LAFEM::VectorMirror<DT_, IT_> ScalarMirror;
          typedef Global::Gate<LocalScalarVector, ScalarMirror> ScalarGate;
 
          typedef Global::Matrix<LocalMatrix, SystemMirror, SystemMirror> GlobalSystemMatrix;
          typedef Global::Filter<LocalSystemFilter, SystemMirror> GlobalSystemFilter;
          typedef Global::Vector<LocalScalarVector, ScalarMirror> GlobalScalarVector;
          typedef Global::Vector<LocalSystemVectorL, SystemMirror> GlobalSystemVectorL;
          typedef Global::Vector<LocalSystemVectorR, SystemMirror> GlobalSystemVectorR;
          typedef Global::Vector<LocalCoordsBuffer, SystemMirror> GlobalCoordsBuffer;
          typedef Global::Transfer<LocalSystemTransfer, SystemMirror> GlobalSystemTransfer;
          typedef Global::Muxer<LocalSystemVectorR, SystemMirror> GlobalSystemMuxer;
 
        public:
          GlobalSystemMatrix matrix_sys;
 
        public:
          template<typename... Args_>
          explicit QuadraticSystemLevel(const int level_index,
          Geometry::RootMeshNode<MeshType>* rmn,
          TrafoType& trafo,
          const std::deque<String>& dirichlet_list,
          const std::deque<String>& slip_list,
          Args_&&... args) :
            BaseClass(level_index, rmn, trafo, dirichlet_list, slip_list, std::forward<Args_>(args)...),
            matrix_sys(&(BaseClass::gate_sys), &(BaseClass::gate_sys),
            BaseClass::local_functional.matrix_sys.clone(LAFEM::CloneMode::Shallow))
            {
            }
 
          virtual ~QuadraticSystemLevel()
          {
          }
 
      }; // class QuadraticSystemLevel<...>
 
      template
      <
        typename DT_, typename IT_,
        template<typename, typename> class Functional_
      >
      class NonlinearSystemLevel : public MeshoptSystemLevel<DT_, IT_, Functional_>
      {
        public:
          typedef DT_ DataType;
          typedef IT_ IndexType;
 
          typedef Functional_<DT_, IT_> LocalFunctional;
 
          typedef MeshoptSystemLevel<DT_, IT_, Functional_> BaseClass;
 
          typedef typename LocalFunctional::SpaceType SpaceType;
          typedef typename LocalFunctional::TrafoType TrafoType;
          typedef typename TrafoType::MeshType MeshType;
 
          typedef typename LocalFunctional::VectorTypeL LocalSystemVectorL;
          typedef typename LocalFunctional::VectorTypeR LocalSystemVectorR;
          typedef typename LocalFunctional::ScalarVectorType LocalScalarVector;
          typedef typename LocalFunctional::CoordsBufferType LocalCoordsBuffer;
          typedef LAFEM::SparseMatrixBWrappedCSR<DT_, IT_, LocalFunctional::MeshType::world_dim>
            LocalSystemTransferMatrix;
          typedef LAFEM::Transfer<LocalSystemTransferMatrix> LocalSystemTransfer;
 
          typedef typename LocalFunctional::FilterType LocalSystemFilter;
          typedef typename LocalFunctional::SlipFilterSequence LocalSlipFilterSequence;
          typedef typename LocalFunctional::DirichletFilterSequence LocalDirichletFilterSequence;
 
          typedef LAFEM::VectorMirror<DT_, IT_> SystemMirror;
          typedef Global::Gate<LocalSystemVectorR, SystemMirror> SystemGate;
          typedef LAFEM::VectorMirror<DT_, IT_> ScalarMirror;
          typedef Global::Gate<LocalScalarVector, ScalarMirror> ScalarGate;
 
          typedef Global::NonlinearFunctional<LocalFunctional, SystemMirror, SystemMirror> GlobalFunctional;
          typedef Global::Filter<LocalSystemFilter, SystemMirror> GlobalSystemFilter;
          typedef Global::Vector<LocalScalarVector, ScalarMirror> GlobalScalarVector;
          typedef Global::Vector<LocalSystemVectorL, SystemMirror> GlobalSystemVectorL;
          typedef Global::Vector<LocalSystemVectorR, SystemMirror> GlobalSystemVectorR;
          typedef Global::Vector<LocalCoordsBuffer, SystemMirror> GlobalCoordsBuffer;
          typedef Global::Transfer<LocalSystemTransfer, SystemMirror> GlobalSystemTransfer;
          typedef Global::Muxer<LocalSystemVectorR, SystemMirror> GlobalSystemMuxer;
 
        public:
          GlobalFunctional global_functional;
 
        public:
          template<typename... Args_>
          explicit NonlinearSystemLevel(const int level_index,
          Geometry::RootMeshNode<MeshType>* rmn,
          TrafoType& trafo,
          const std::deque<String>& dirichlet_list,
          const std::deque<String>& slip_list,
          Args_&&... args) :
            BaseClass(level_index, rmn, trafo, dirichlet_list, slip_list, std::forward<Args_>(args)...),
            global_functional(&(BaseClass::gate_sys), &(BaseClass::gate_sys), BaseClass::local_functional)
            {
            }
 
          virtual ~NonlinearSystemLevel()
          {
          }
 
      }; // class NonlinearSystemLevel<...>
 
    } // namespace Meshopt
 
  } // namespace Control
} //namespace FEAT