mllgeometry.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef MLLGEOMETRY_HH
#define MLLGEOMETRY_HH
 
#include <limits>
#include <stack>
 
#include "dune/common/fmatrix.hh"
#include <dune/geometry/referenceelements.hh>
 
#include "fem/fixdune.hh"
 
namespace Kaskade
{
  template <class Grid>
  class MultiLevelLocalGeometry
  {
    typedef typename Grid::GlobalIdSet::IdType Id;
    typedef typename Grid::template Codim<0>::Entity       Cell;
 
  protected:
    Cell child, ancestor;
 
  private:
 
    // Changed in DUNE 2.1
    // typename Grid::template Codim<0>::GlobalIdSet const& idSet;
    typename Grid::GlobalIdSet const& idSet;
    int spanwidth;
    bool geometricallyNested, transformViaGlobal;
 
  public:
 
    enum Direction{
      ChildIsGlobal,
      FatherIsGlobal
    };
 
    typedef typename Grid::template Codim<0>::Entity Entity;
 
 
  private:
 
    Direction dir;
 
    typedef typename Grid::ctype ctype;
    static int const dim = Grid::dimension;
    static int const dimw = Grid::dimensionworld;
 
  public:
 
    MultiLevelLocalGeometry(Grid const &grid,
                            Cell const& child_,
                            Cell const& ancestor_,
                            MultiLevelLocalGeometry<Grid>::Direction direction_,
                            bool geometricallyNested_= (dim==dimw),
                            ctype tol_=-100):
      child(child_),
      ancestor(ancestor_),
      idSet(grid.globalIdSet()),
      geometricallyNested(geometricallyNested_),
      dir(direction_),
      tol(tol_<=-100? 100*std::numeric_limits<ctype>::epsilon(): tol_)
    {
      Cell ci(child);
      spanwidth = 0;
      while (idSet.id(ci) != idSet.id(ancestor))
      {
        ci = ci.father();
        spanwidth++;
      }
      transformViaGlobal = geometricallyNested && (spanwidth >= 3);
    }
 
    Entity const& localEntity() const { return dir==ChildIsGlobal? ancestor: *child; }
 
 
    void global(std::vector<Dune::FieldVector<ctype,dim> > & x,
                std::vector<int> &componentsInside) const
    {
      if (dir==ChildIsGlobal) {
        if(transformViaGlobal) childToFatherViaGlobal(x);
        else                   childToFather(x);
      }
      else {
        if(transformViaGlobal) fatherToChildViaGlobal(x,componentsInside);
        else                   fatherToChild(x,componentsInside);
      }
    }
 
    void local(std::vector<Dune::FieldVector<ctype,dim>>& x,
               std::vector<int>& componentsInside) const
    {
      if (dir==FatherIsGlobal) {
        if(transformViaGlobal) fatherToChildViaGlobal(x,componentsInside);
        else                   fatherToChild(x,componentsInside);
      }
      else {
        if(transformViaGlobal) childToFatherViaGlobal(x);
        else                   childToFather(x);
      }
    }
 
    Dune::FieldMatrix<ctype,dim,dim> jacobianInverseTransposed (Dune::FieldVector<ctype,dim> const& xLocal) const
    {
      if (transformViaGlobal)
        if (dir==FatherIsGlobal) {
          Dune::FieldVector<ctype,dim> xGlobal = ancestor.geometry().local(child.geometry().global(xLocal));
          return ancestor.geometry().jacobianTransposed(xGlobal) * child.geometry().jacobianInverseTransposed(xLocal);
        } else {
          Dune::FieldVector<ctype,dim> xGlobal = child.geometry().local(ancestor.geometry().global(xLocal));
          Dune::FieldMatrix<ctype,dim,dim> Jct = child.geometry().jacobianInverseTransposed(xGlobal);
          Jct.invert();
          return Jct * ancestor.geometry().jacobianInverseTransposed(xLocal);
        }
      else {
        Dune::FieldMatrix<ctype,dim,dim> invt = unitMatrix<ctype,dim>();
        Dune::FieldVector<ctype,dim> x = xLocal;
        Cell ci(child);
        if (dir==FatherIsGlobal) {
          while(idSet.id(ci) != idSet.id(ancestor)) {
            invt.leftmultiply(ci.geometryInFather().jacobianInverseTransposed(x));
            x = ci.geometryInFather().global(x);
          }
        } else {
          std::stack<Cell> stack;
          while(idSet.id(ci) != idSet.id(ancestor)) {
            stack.push(ci);
            ci = ci.father();
          }
          while (!stack.empty()) {
            x = stack.top()->geometryInFather().local(x);
            // TODO: throw exception if x not in ci
            invt.rightmultiply(stack.top()->geometryInFather().jacobianInverseTransposed(x));
            stack.pop();
          }
          invt.invert();
        }
        return invt;
      }
    }
 
  private:
 
    void childToFather(std::vector< Dune::FieldVector<ctype,dim> >& x) const
    {
      Cell ci(child);
      while(idSet.id(ci) != idSet.id(ancestor))
      {
        auto const& geo = ci.geometryInFather();
        for(auto& xi: x)
          xi = geo.global(xi);
        ci = ci.father();
      }
    }
 
    void childToFatherViaGlobal(std::vector<Dune::FieldVector<ctype,dim>>& x) const
    {
      auto const& geoChild = child.geometry();
      auto const& geoAnces = ancestor.geometry();
      for(auto& xi: x) 
        xi = geoAnces.local(geoChild.global(xi));
    }
 
    void fatherToChild(std::vector<Dune::FieldVector<ctype,dim>>& x,
                       std::vector<int>& componentsInside) const
    {
      assert(x.size()==componentsInside.size());
 
      // We need to go from father cells up to child cells. However, we can only obtain the father
      // directly. Thus we have to traverse from child to ancestor and store the whole line
      // of fathers to be traced back in reverse direction.
      Cell ci(child);
      std::stack<Cell> ciStack;
      while(idSet.id(ci) != idSet.id(ancestor))
      {
        ciStack.push(ci);
        ci = ci.father();
      }
      // Now ciStack.top()->father is ancestor.
 
      // Trace back the line from ancestor to child and transform the coordinates given in x
      // to the childrens' local coordinates on the way.
      while(!ciStack.empty())
      {
        ci = ciStack.top();
        auto vi = x.begin();
        auto compi = componentsInside.begin();
        while (vi!=x.end())
        {
          *vi = ci.geometryInFather().local(*vi);
          if (checkInside(ci.type(),*vi)<tol)   // (probably) inside:
          {                                     //     go to next location
            ++vi;
            ++compi;
          }
          else                                  // (almost sure) outside:
          {                                     //     remove the location
            vi = x.erase(vi);
            compi = componentsInside.erase(compi);
          }
        }
        ciStack.pop();
      }
    }
 
    void fatherToChildViaGlobal(std::vector<Dune::FieldVector<ctype,dim>>& x,
                                std::vector<int>& componentsInside) const
    {
      for(auto& xi: x)
        xi = child.geometry().local(ancestor.geometry().global(xi));
 
      auto vi = x.begin();
      auto compi = componentsInside.begin(); 
      while (vi!=x.end()) {
        if (checkInside(child.type(),*vi)<tol) { // (probably) inside
          ++vi;
          ++compi;
        } else { // (almost sure) outside
          vi = x.erase(vi);
          compi = componentsInside.erase(compi);
        }
      }
    }
 
  private:
    ctype tol;
  };
 
  template<class Grid>
  class LocalGeometryInCoarseGridAncestor : public MultiLevelLocalGeometry<Grid>
  {
    typedef typename Grid::template Codim<0>::Entity Cell;
    typedef typename Grid::ctype ctype;
    static int const dim = Grid::dimension;
  public:
 
    LocalGeometryInCoarseGridAncestor(Grid const &grid,
                                      Cell const& child,
                                      bool geometricallyNested=false) :
    MultiLevelLocalGeometry<Grid>(grid,child,coarseGridAncestor(child),MultiLevelLocalGeometry<Grid>::ChildIsGlobal,geometricallyNested)
    {
    }
 
    Dune::FieldVector<ctype,dim> global(Dune::FieldVector<ctype,dim>const& x) const
    {
      std::vector<Dune::FieldVector<ctype,dim> > xi(1,x);
      std::vector<int> isi(1,0);
      MultiLevelLocalGeometry<Grid>::global(xi,isi);
      return xi[0];
    }
 
    Cell const& getFather() const  { return this->ancestor; }
 
  private:
 
    static Cell coarseGridAncestor(Cell const& child)
    {
      Cell ancestor(child);
 
      while (ancestor.level()>0) 
        ancestor = ancestor.father();
      return ancestor;
    }
  };
} // end of namespace Kaskade
 
#endif