extractor.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-GPL-2.0-only-with-dune-grid-glue-exception
/*
 *  Filename:    extractor.hh
 *  Version:     1.0
 *  Created on:  Oct 05, 2009
 *  Author:      Christian Engwer
 *  ---------------------------------
 *  Project:     dune-grid-glue
 *  Description: base class for all grid extractors
 *
 */
#ifndef DUNE_GRIDGLUE_EXTRACTORS_EXTRACTOR_HH
#define DUNE_GRIDGLUE_EXTRACTORS_EXTRACTOR_HH
 
#include <vector>
#include <map>
#include <algorithm>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/version.hh>
#include <dune/grid/common/geometry.hh>
#include <dune/grid/common/grid.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/geometry/multilineargeometry.hh>
 
namespace Dune {
 
  namespace GridGlue {
 
template<typename GV, int cd>
class Extractor
{
 
public:
 
  static constexpr auto dimworld = GV::dimensionworld;
  static constexpr auto dim      = GV::dimension;
  static constexpr auto codim    = cd;
 
  static constexpr int cube_corners = 1 << (dim-codim);
 
  typedef GV GridView;
  typedef typename GridView::Grid Grid;
 
  typedef typename GV::Grid::ctype ctype;
  typedef Dune::FieldVector<ctype, dimworld>                       Coords;
  typedef Dune::FieldVector<ctype, dim>                            LocalCoords;
 
  typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
  typedef typename Vertex::EntitySeed VertexSeed;
 
  typedef typename GV::Traits::template Codim<0>::Entity Element;
  typedef typename Element::EntitySeed ElementSeed;
 
  typedef std::vector<unsigned int>                                VertexVector;
  using CellMapper = MultipleCodimMultipleGeomTypeMapper<GridView>;
 
  // typedef typename CellMapper::IndexType                               IndexType;
  typedef int IndexType;
public:
 
  // transformations
  typedef Dune::MultiLinearGeometry<ctype, dim-codim, dimworld> Geometry;
  typedef Dune::MultiLinearGeometry<ctype, dim-codim, dim>      LocalGeometry;
 
protected:
  /************************** PRIVATE SUBCLASSES **********************/
 
  struct CornerInfo
  {
    unsigned int idx : 28;   
    unsigned int num : 4;   
  };
 
  struct CoordinateInfo
  {
    CoordinateInfo()
    {}
 
    CoordinateInfo(unsigned int index_, IndexType vtxindex_)
      : vtxindex(vtxindex_), index(index_)
    {}
 
    IndexType vtxindex;
 
    Coords coord;
 
    unsigned int index;
  };
 
  struct VertexInfo
  {
    VertexInfo(unsigned int idx_, const Vertex& p_) : idx(idx_), p(p_.seed())
    {}
    unsigned int idx;
    VertexSeed p;
  };
 
 
  struct ElementInfo
  {
    ElementInfo(unsigned int idx_, const Element& p_, unsigned int f_) : idx(idx_), faces(f_), p(p_.seed())
    {}
 
    unsigned int idx : 28;
 
    unsigned int faces : 4;
 
    ElementSeed p;
  };
 
 
  struct SubEntityInfo
  {
    SubEntityInfo()
    /*
     * TODO: move default value of `geometryType_` to member declaration
     * when removing support for older dune-geometry
     */
      : geometryType_(GeometryTypes::simplex(dim-codim))
      {}
 
    SubEntityInfo(IndexType parent_, unsigned int num_in_parent_,
                  const Dune::GeometryType& geometryType)
      : parent(parent_), num_in_parent(num_in_parent_), geometryType_(geometryType)
    {}
 
    unsigned int nCorners() const
    {
      return Dune::ReferenceElements<ctype, dim-codim>::general(geometryType_).size(dim-codim);
    }
 
    IndexType parent;
 
    unsigned int num_in_parent : 3;
 
    Dune::GeometryType geometryType_;
 
    CornerInfo corners[cube_corners];     // sim = number of vertices in a simplex
  };
 
 
  typedef std::map<IndexType, ElementInfo>  ElementInfoMap;
  typedef std::map<IndexType, VertexInfo>   VertexInfoMap;
 
  /************************** MEMBER VARIABLES ************************/
 
  const GridView                gv_;
 
  /*        Geometrical and Topological Information                */
 
  std::vector<CoordinateInfo>   coords_;
 
  std::vector<SubEntityInfo>    subEntities_;
 
  VertexInfoMap vtxInfo_;
 
  ElementInfoMap elmtInfo_;
 
  CellMapper cellMapper_;
 
public:
 
  /*  C O N S T R U C T O R S   A N D   D E S T R U C T O R S  */
 
  Extractor(const GV& gv)
    :  gv_(gv)
    , cellMapper_(gv, mcmgElementLayout())
  {}
 
  /*  F U N C T I O N A L I T Y  */
 
  void clear()
  {
    coords_.clear();
    coords_.shrink_to_fit();
 
    subEntities_.clear();
    subEntities_.shrink_to_fit();
 
    // ...then clear the maps themselves, too
    vtxInfo_.clear();
    elmtInfo_.clear();
  }
 
 
  /*  G E T T E R S  */
 
  void getCoords(std::vector<Dune::FieldVector<ctype, dimworld> >& coords) const
  {
    coords.resize(coords_.size());
    for (unsigned int i = 0; i < coords_.size(); ++i)
      coords[i] = coords_[i].coord;
  }
 
 
  unsigned int nCoords() const
  {
    return coords_.size();
  }
 
  void getGeometryTypes(std::vector<Dune::GeometryType>& geometryTypes) const
  {
    geometryTypes.resize(subEntities_.size());
    for (size_t i=0; i<subEntities_.size(); i++)
      geometryTypes[i] = subEntities_[i].geometryType_;
  }
 
 
  void getFaces(std::vector<VertexVector>& faces) const
  {
    faces.resize(subEntities_.size());
    for (unsigned int i = 0; i < subEntities_.size(); ++i) {
      faces[i].resize(subEntities_[i].nCorners());
      for (unsigned int j = 0; j < subEntities_[i].nCorners(); ++j)
        faces[i][j] = subEntities_[i].corners[j].idx;
    }
  }
 
 
  bool faceIndices(const Element& e, int& first, int& count) const
  {
    typename ElementInfoMap::const_iterator it =
      elmtInfo_.find(cellMapper_.map(e));
    if (it == elmtInfo_.end())
    {
      first = -1;
      count = 0;
      return false;
    }
    // the iterator is valid, fill the out params
    first = it->second.idx;
    count = it->second.faces;
    return true;
  }
 
 
  int indexInInside(unsigned int index) const
  {
    return index < subEntities_.size() ? subEntities_[index].num_in_parent : -1;
  }
 
  // /**
  //  * @brief tests that a given entry in the extraction set does have local couplings
  //  * @todo parallel interface
  //  */
  // bool contains (unsigned int global, unsigned int & local) const
  // {
  //     local = global;
  //     return true;
  // }
 
  const GridView & gridView() const
  {
    return gv_;
  }
 
  const Grid& grid() const
  {
    return gv_.grid();
  }
 
  Element
  element(unsigned int index) const
  {
    if (index >= subEntities_.size())
      DUNE_THROW(Dune::GridError, "invalid face index");
    const ElementSeed seed = elmtInfo_.at(subEntities_[index].parent).p;
    return grid().entity(seed);
  }
 
#if 1
  Vertex
  vertex(unsigned int index) const
  {
    if (index >= coords_.size())
      DUNE_THROW(Dune::GridError, "invalid coordinate index");
    const VertexSeed seed = vtxInfo_.at(coords_[index].vtxindex).p;
    return grid().entity(seed);
  }
#endif
 
  Geometry geometry(unsigned int index) const;
 
  LocalGeometry geometryLocal(unsigned int index) const;
 
};
 
 
template<typename GV, int cd>
typename Extractor<GV,cd>::Geometry Extractor<GV,cd>::geometry(unsigned int index) const
{
  std::vector<Coords> corners(subEntities_[index].nCorners());
  for (unsigned int i = 0; i < subEntities_[index].nCorners(); ++i)
    corners[i] = coords_[subEntities_[index].corners[i].idx].coord;
 
  return Geometry(subEntities_[index].geometryType_, corners);
}
 
 
template<typename GV, int cd>
typename Extractor<GV,cd>::LocalGeometry Extractor<GV,cd>::geometryLocal(unsigned int index) const
{
  std::vector<LocalCoords> corners(subEntities_[index].nCorners());
 
  // get face info
  const SubEntityInfo & face = subEntities_[index];
  Dune::GeometryType facetype = subEntities_[index].geometryType_;
 
  // get reference element
  const auto elmtseed = elmtInfo_.at(face.parent).p;
  const auto elmt = grid().entity(elmtseed);
  const Dune::GeometryType celltype = elmt.type();
  const auto& re = Dune::ReferenceElements<ctype, dim>::general(celltype);
  for (unsigned int i = 0; i < subEntities_[index].nCorners(); ++i)
    corners[i] = re.position(face.corners[i].num,dim);
 
  return LocalGeometry(facetype, corners);
}
 
}  // namespace GridGlue
 
}  // namespace Dune
 
#endif // DUNE_GRIDGLUE_EXTRACTORS_EXTRACTOR_HH