trace_assembler.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
 
// includes, FEAT
#include <kernel/assembly/asm_traits.hpp>
#include <kernel/geometry/mesh_part.hpp>
#include <kernel/geometry/intern/face_ref_trafo.hpp>
#include <kernel/geometry/intern/congruency_sampler.hpp>
#include <kernel/geometry/intern/congruency_trafo.hpp>
#include <kernel/geometry/intern/index_representative.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
 
#include <vector>
 
namespace FEAT
{
  namespace Assembly
  {
    namespace Intern
    {
      template<typename Outer_, typename Shape_, int cell_dim_>
      class CompIndexMap
      {
      protected:
        const Outer_& _outer;
        int _idx;
 
      public:
        explicit CompIndexMap(const Outer_& outer, int idx) : _outer(outer), _idx(idx) {}
 
        Index operator[](int i) const
        {
          return _outer[Geometry::Intern::FaceIndexMapping<Shape_, cell_dim_, 0>::map(_idx,i)];
        }
      }; // class CompIndexMap
    } // namespace Intern
 
#ifdef DOXYGEN
    class TraceAssemblyJob
    {
    public:
      class Task
      {
      public:
        static constexpr bool assemble_pairwise = true or false;
 
        static constexpr bool need_scatter = true or false;
 
        static constexpr bool need_combine = true or false;
 
        explicit Task(TraceAssemblyJob& job);
 
        void prepare(Index facet, Index cell, int local_facet, int facet_ori);
 
        void prepare(Index facet, Index cell_a, Index cell_b, int local_facet_a, int local_facet_b, int facet_ori_a, int facet_ori_b);
 
        void assemble();
 
        void scatter();
 
        void finish();
 
        void combine();
      }; // class Task
    }; // class TraceAssemblyJob
#endif // DOXYGEN
 
    template<typename Trafo_>
    class TraceAssembler
    {
    public:
      typedef Trafo_ TrafoType;
      typedef typename TrafoType::MeshType MeshType;
      typedef typename TrafoType::ShapeType ShapeType;
 
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int facet_dim = shape_dim-1;
 
    protected:
      const TrafoType& _trafo;
      std::vector<int> _facet_mask, _cell_facet, _facet_ori;
      std::vector<Index> _cells, _facets;
 
    public:
      explicit TraceAssembler(const TrafoType& trafo) :
        _trafo(trafo),
        _facet_mask(trafo.get_mesh().get_num_entities(facet_dim), 0)
      {
      }
 
      const TrafoType& get_trafo() const
      {
        return _trafo;
      }
 
      void clear()
      {
        _cell_facet.clear();
        _facet_ori.clear();
        _cells.clear();
        _facets.clear();
        for(auto& f : _facets)
          f = Index(0);
      }
 
      void add_facet(Index ifacet)
      {
        ASSERTM(ifacet < Index(_facet_mask.size()), "invalid facet index");
        _facet_mask.at(ifacet) = 1;
      }
 
      void add_mesh_part(const Geometry::MeshPart<MeshType>& mesh_part)
      {
        const auto& trg = mesh_part.template get_target_set<facet_dim>();
        for(Index i(0); i < trg.get_num_entities(); ++i)
          _facet_mask.at(trg[i]) = 1;
      }
 
      void compile()
      {
        _cell_facet.clear();
        _facet_ori.clear();
        _cells.clear();
        _facets.clear();
 
        // build elements-at-facet graph
        Adjacency::Graph elem_at_facet(Adjacency::RenderType::injectify_transpose,
          _trafo.get_mesh().template get_index_set<shape_dim, facet_dim>());
 
        // loop over all facets
        for(Index iface(0); iface < Index(_facet_mask.size()); ++iface)
        {
          // use this facet?
          if(_facet_mask[iface] == 0)
            continue;
 
          // ensure that this is a boundary facet if required
          //if(only_boundary && (elem_at_facet.degree(iface) != Index(1)))
            //XABORTM("facet is adjacent to more than 1 element");
 
          // add all elements
          for(auto it = elem_at_facet.image_begin(iface); it != elem_at_facet.image_end(iface); ++it)
          {
            Index icell = *it;
 
            // try to compute local facet index and orientation
            int loc_face(0), face_ori(0);
            if(!_find_local_facet(iface, icell, loc_face, face_ori))
              XABORTM("failed to find local facet");
 
            // alright, add this facet to our list
            _facets.push_back(iface);
            _cells.push_back(icell);
            _cell_facet.push_back(loc_face);
            _facet_ori.push_back(face_ori);
          }
        }
      }
 
      void compile_all_facets(bool inner, bool outer)
      {
        _cell_facet.clear();
        _facet_ori.clear();
        _cells.clear();
        _facets.clear();
 
        // build elements-at-facet graph
        Adjacency::Graph elem_at_facet(Adjacency::RenderType::injectify_transpose,
          _trafo.get_mesh().template get_index_set<shape_dim, facet_dim>());
 
        // loop over all facets
        for(Index iface(0); iface < Index(_facet_mask.size()); ++iface)
        {
          // how many elements are adjacent to this facet?
          const int degree = int(elem_at_facet.degree(iface));
          XASSERT(degree > 0);
          XASSERT(degree < 3);
 
          // outer facet with only 1 adjacent element?
          if((degree == 1) && !outer)
            continue;
 
          // inner facet with 2 adjacent elements?
          if((degree == 2) && !inner)
            continue;
 
          // add all elements
          for(auto it = elem_at_facet.image_begin(iface); it != elem_at_facet.image_end(iface); ++it)
          {
            Index icell = *it;
 
            // try to compute local facet index and orientation
            int loc_face(0), face_ori(0);
            if(!_find_local_facet(iface, icell, loc_face, face_ori))
              XABORTM("failed to find local facet");
 
            // alright, add this facet to our list
            _facets.push_back(iface);
            _cells.push_back(icell);
            _cell_facet.push_back(loc_face);
            _facet_ori.push_back(face_ori);
          }
        }
      }
 
      template<typename Job_>
      void assemble(Job_& job)
      {
        typedef typename Job_::Task TaskType;
 
        //XASSERTM(_compiled, "assembler has not been compiled yet");
 
        // no facets to assemble on?
        if(this->_facets.empty())
          return;
 
        // create a task for this thread
        std::unique_ptr<TaskType> task(new TaskType(job));
 
        // note: in the case of inner facets, a single inner facet is stored twice in the _facets
        // vector in two successive elements, i.e. one has _facets[i] == _facets[i+1], and the
        // corresponding two adjacent cells are given by _cells[i] and _cells[i+1]
 
        // loop over all facets that have been added to the assembler
        for(Index fi(0); fi < Index(_facets.size()); ++fi)
        {
          // prepare task -- either pairwise or not
          if constexpr(TaskType::assemble_pairwise)
          {
            // task supports pairwise assembly, so let's check if the next facet is an inner facet
            Index fj = fi+1u;
            if((fj < Index(_facets.size())) && (_facets[fi] == _facets[fj]))
            {
              // that's an inner facet, so assemble pairwise here
              task->prepare(_facets[fi], _cells[fi], _cells[fj], _cell_facet[fi], _cell_facet[fj], _facet_ori[fi], _facet_ori[fj]);
 
              // make sure we don't try to assemble on fi+1 again
              ++fi;
            }
            else
            {
              // not an inner facet, so call single-cell prepare overload
              task->prepare(_facets[fi], _cells[fi], _cell_facet[fi], _facet_ori[fi]);
            }
          }
          else
          {
            // task does not support pairwise assembly, so always call the single-cell prepare overload
            task->prepare(_facets[fi], _cells[fi], _cell_facet[fi], _facet_ori[fi]);
          }
 
          // assemble task
          task->assemble();
 
          // scatter
          if(task->need_scatter)
          {
            task->scatter();
          }
 
          // finish
          task->finish();
        }
 
        // do we have to combine the assembly?
        if(task->need_combine)
        {
          // combine the assembly
          task->combine();
        }
 
        // delete task object
        task.reset();
      }
 
    protected:
      bool _find_local_facet(Index face, Index cell, int& facet, int& ori)
      {
        typedef typename Shape::FaceTraits<ShapeType, facet_dim>::ShapeType FacetType;
        static constexpr int num_facets = Shape::FaceTraits<ShapeType, shape_dim-1>::count;
        static constexpr int num_vaf = Shape::FaceTraits<FacetType, 0>::count;
 
        typedef Geometry::Intern::IndexRepresentative<FacetType> FacetRepr;
 
        const auto& vert_at_elem = _trafo.get_mesh().template get_index_set<shape_dim, 0>();
        const auto& vert_at_face = _trafo.get_mesh().template get_index_set<facet_dim, 0>();
 
        Index face_verts[num_vaf];
        FacetRepr::compute(face_verts, vert_at_face[face]);
 
        for(int li(0); li < num_facets; ++li)
        {
          Intern::CompIndexMap<decltype(vert_at_elem[0]), ShapeType, shape_dim-1> cim(vert_at_elem[cell], li);
          Index lf_verts[num_vaf];
          FacetRepr::compute(lf_verts, cim);
          bool match = true;
          for(int k(0); k < num_vaf; ++k)
          {
            if(lf_verts[k] != face_verts[k])
            {
              match = false;
              break;
            }
          }
          if(match)
          {
            ori = Geometry::Intern::CongruencySampler<FacetType>::compare(vert_at_face[face], cim);
            facet = li;
            return true;
          }
        }
        return false;
      }
    }; // class TraceAssembler<...>
  } // namespace Assembly
} // namespace FEAT