MeshFunction_Level.h

/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_MESH_FUNCTION_LEVEL_H
#define EXTENSIONS_MESH_FUNCTION_LEVEL_H

namespace AMDiS { namespace extensions {


/** \brief
 * Abstract class that can be passed to RefinementLevel* as indicator where
 * to refine the mesh up to which level. It is an AbstractFunction that 
 * overloads the operator() method to return a level or a meshsize depending
 * on the coords/data passed to the operator.
 * You can switch between meshsize and level with the methods hToLevel(double) and 
 * levelToH(int)
 **/
template<typename T, typename T2>
class MeshRefinementFunction : public AbstractFunction<T2, T>
{
public:

  MeshRefinementFunction(Mesh* mesh_) :
    AbstractFunction<T2, T>(0),
    mesh(mesh_), 
    globalSize(0)
  {
      h0 = getMacroMeshSize(mesh);
      reduction = 1.0 / std::sqrt(2.0); // if dim==2
  }

  int getGlobalSize() { return globalSize; }

  double meshSize() { return h0; }
  
  virtual void init(ElInfo* elInfo) {}

  virtual T2 operator()(const T &value) const { return globalSize; }
  virtual T2 getLevel() const { return globalSize; }
  
  virtual double indicator(const T &value) const { return 1.0; }

protected:

  int hToLevel(double h) {
      int level = static_cast<int>(std::floor(std::log(h / h0) / std::log(reduction)));
      return level;
  }

  double levelToH(int level) {
      double h = std::pow(reduction,level)*h0;
      return h;
  }

  double getMacroMeshSize(Mesh* mesh) {
      FixVec<WorldVector<double>, VERTEX> coords = mesh->getMacroElement(0)->getCoord();
      double h = 0.0;
      for (int i = 0; i < coords.getSize(); ++i)
          for (int j = i + 1; j < coords.getSize(); ++j)
              h = std::max(h, norm(coords[i]-coords[j]));
      return h;
  }

protected:

  Mesh* mesh;

  int globalSize;

  double h0;
  double reduction;
};


/** \brief
 * Base-Implementation of a MeshRefinementFunction for refinement of phase-field
 * interfaces. operator() must be overwritten in subclasses. Class provides
 * basic parameters for the refinement like refinement-level values on interface
 * or in the interesting domain and defines the refinement range for the phase-
 * field variable, i.e. interface is defined in [0.05, 0.95]
 **/
template<typename T>
class PhaseFieldRefinementBase : public MeshRefinementFunction<T, int>
{
public:

  PhaseFieldRefinementBase(Mesh *mesh, std::string name = "mesh->refinement") : 
    MeshRefinementFunction<T, int>(mesh),
    lInner(10),
    lOuter(10),
    lInterface(14),
    minPhase(0.05),
    maxPhase(0.95),
    minOuterPhase(0.001),
    maxOuterPhase(0.999)
  {
      Parameters::get(name + "->level in inner domain",lInner);
      Parameters::get(name + "->level in outer domain",lOuter);
      Parameters::get(name + "->level on interface",lInterface);

      lInterface-= mesh->getMacroElementLevel();
      lInner-= mesh->getMacroElementLevel();
      lOuter-= mesh->getMacroElementLevel();

      int local_globalSize = 10;
      Parameters::get(name + "->initial level", local_globalSize);
      MeshRefinementFunction<T, int>::globalSize = local_globalSize;

      Parameters::get(name + "->min interface value",minPhase);
      Parameters::get(name + "->max interface value",maxPhase);
      Parameters::get(name + "->min outer interface value",minOuterPhase);
      Parameters::get(name + "->max outer interface value",maxOuterPhase);
  }

  double meshSize() { return levelToH(MeshRefinementFunction<T, int>::globalSize); }

protected:

  int lInner;
  int lOuter;
  int lInterface;
  
  double minPhase;
  double maxPhase;
  double minOuterPhase;
  double maxOuterPhase;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 **/
class PhaseFieldRefinement : public PhaseFieldRefinementBase<double>
{
public:

    PhaseFieldRefinement(Mesh* mesh_, std::string name = "mesh->refinement") 
      : PhaseFieldRefinementBase<double>(mesh_, name) 
    { };
    
    int operator()(const double& phase) const {
      int result= lOuter;
      if (minPhase < phase && phase < maxPhase)
        result= lInterface;         // auf dem Interface
      else if (phase > minOuterPhase && phase <= minPhase)
        result= std::max(lOuter, static_cast<int>(std::floor((lOuter+lInterface)/2.0)));
      else if (phase < maxOuterPhase && phase >= maxPhase)
        result= std::max(lInner, static_cast<int>(std::floor((lInner+lInterface)/2.0)));
      else if (phase > (minPhase+maxPhase)/2.0)
        result= lInner;
      return result;
    }
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 * Additionally to the phase-field value a list of refinement points can be given
 * where the mesh should be refined up to the a given level, e.g. in corners of 
 * the macro-mesh.
 **/
class PhaseFieldCoordsRefinement : public PhaseFieldRefinementBase< std::pair<WorldVector<double>, double> >
{
public:
    PhaseFieldCoordsRefinement(Mesh* mesh_, std::vector<WorldVector<double> > points_, double radius_, std::string name = "mesh->refinement") 
      : PhaseFieldRefinementBase< std::pair<WorldVector<double>, double> >(mesh_, name),
	lPoints(14),
	points(points_),
	radius(radius_)
    {
        Parameters::get(name + "->level on points", lPoints);
        lPoints-= mesh->getMacroElementLevel();
    }
    
    int operator()(const std::pair<WorldVector<double>, double>& data) const {
      double phase = data.second;
      int result= lOuter;
      if (minPhase < phase && phase < maxPhase)
        result = lInterface;         // auf dem Interface
      else if ((phase > minOuterPhase && phase <= minPhase) || (phase < maxOuterPhase && phase >= maxPhase))
        result = std::max(lInner, static_cast<int>(std::floor((lOuter+lInterface)/2.0)));
      else if (phase > 0.5)
        result = lInner;

      WorldVector<double> x = data.first;
      double minDist = 1.e15;
      for (unsigned i = 0; i < points.size(); ++i) {
        minDist = std::min(minDist, norm(points[i]-x));
      }
      double lambda = std::max(0.0, 1.0 - minDist/radius);
      result = std::max(result, static_cast<int>(std::floor(lOuter+lambda*(lPoints-lOuter))));

      return result;
    }
    
private:

    int lPoints;

    std::vector<WorldVector<double> > points;

    double radius;
};


/** \brief
 * Implementation of MeshRefinementFunction. Provides the method operator()
 * that defines a list of refinement points where the mesh should be refined 
 * up to the a given level, e.g. in corners of the macro-mesh.
 **/
class CoordsRefinement : public MeshRefinementFunction< WorldVector<double>, int >
{
public:
    CoordsRefinement(Mesh* mesh_, std::vector<WorldVector<double> > points_, double radius_, std::string name = "mesh->refinement") 
      : MeshRefinementFunction< WorldVector<double>, int >(mesh_),
	points(points_),
	radius(radius_)
    {
        lInner= 15; lPoints = 20;
        Parameters::get(name + "->level on points",lPoints);
        Parameters::get(name + "->level in inner domain",lInner);
        lPoints-= mesh->getMacroElementLevel();
        lInner-= mesh->getMacroElementLevel();

        globalSize= lInner;
    }
    
    int operator()(const WorldVector<double>& x) const {
      double minDist = 1.e15;
      for (unsigned i = 0; i < points.size(); ++i) {
        minDist = std::min(minDist, norm(points[i]-x));
      }
      int result = lInner;
      double lambda = std::max(0.0, 1.0 - minDist/radius);
      result = std::max(result, static_cast<int>(std::floor(lInner+lambda*(lPoints-lInner))));

      return result;
    }
    double meshSize() { return levelToH(lInner); }
    
private:
    int lInner, lPoints;
    std::vector<WorldVector<double> > points;
    double radius;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the level of refinement depending on a list of phase-fields, e.g.
 * a diffuse-domain representation and a Cahn-Hilliard field.
 **/
class PhaseFieldChRefinement : public PhaseFieldRefinementBase< std::vector<double> >
{
public:
  PhaseFieldChRefinement(Mesh* mesh_, std::string name = "mesh->refinement") 
    : PhaseFieldRefinementBase< std::vector<double> >(mesh_, name),
      lInterfaceDomain(14)
  {
    Parameters::get(name + "->level on domain interface",lInterfaceDomain);
    lInterfaceDomain-= mesh->getMacroElementLevel();
    
    Parameters::get(name + "->initial level", globalSize);
  }
  
  int operator()(const std::vector<double>& phases) const 
  {
    int result= lOuter;
    for (unsigned i = 0; i < phases.size(); ++i) {
      int localResult= lOuter;
      if (minPhase < phases[i] && phases[i] < maxPhase
        && (i < 1 || phases[i-1] > 0.5))
        localResult= (i==0 ? lInterfaceDomain : lInterface);         // auf dem Interface
      else if (phases[i] > 0.5
        && (i < 1 || phases[i-1] > 0.5))
        localResult= lInner;            // im Innern des Gebietes
      if (((phases[i] > minOuterPhase && phases[i] <= minPhase) || (phases[i] < maxOuterPhase && phases[i] >= maxPhase))
        && (i < 1 || phases[i-1] > 0.5))
        localResult= static_cast<int>(std::floor((lInner + lInterface) / 2.0));
      result= std::max(result, localResult);
    }
    return result;
  }
    
private:

  int lInterfaceDomain;
};


/** \brief
 * Base-Implementation of a MeshRefinementFunction for refinement of phase-field
 * interfaces. operator() must be overwritten in subclasses. Class provides
 * basic parameters for the refinement like refinement-level values on interface
 * or in the interesting domain and defines the refinement range for the phase-
 * field variable, i.e. interface is defined in [0.05, 0.95]
 **/
template<typename T>
class SignedDistRefinementBase : public MeshRefinementFunction<T, int>
{
public:

  SignedDistRefinementBase(Mesh *mesh, std::string name = "mesh->refinement") :
    MeshRefinementFunction<T, int>(mesh),
    name(name),
    lInner(10),
    lOuter(10),
    lInterface(14),
    interfaceWidth(0.2),
    fadeOutWidth(0.0),
    signInInnerDomain(-1.0)
  {
      Parameters::get(name + "->level in inner domain",lInner);
      Parameters::get(name + "->level in outer domain",lOuter);
      Parameters::get(name + "->level on interface",lInterface);

      lInterface-= mesh->getMacroElementLevel();
      lInner-= mesh->getMacroElementLevel();
      lOuter-= mesh->getMacroElementLevel();

      int local_globalSize = 10;
      Parameters::get(name + "->initial level", local_globalSize);
      MeshRefinementFunction<T, int>::globalSize = local_globalSize;

      Parameters::get(name + "->interface width",interfaceWidth);
      Parameters::get(name + "->fade out width",fadeOutWidth);
      Parameters::get(name + "->sign in inner domain",signInInnerDomain);
  }

  double meshSize() { return levelToH(MeshRefinementFunction<T, int>::globalSize); }

protected:
  std::string name;

  int lInner;
  int lOuter;
  int lInterface;

  double interfaceWidth;
  double fadeOutWidth;
  double signInInnerDomain;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 **/
class SignedDistRefinement : public SignedDistRefinementBase<double>
{
public:

    SignedDistRefinement(Mesh* mesh_, std::string name = "mesh->refinement") : SignedDistRefinementBase<double>(mesh_, name) {};

    int operator()(const double& dist) const {
      int result= lOuter;
      if (std::abs(dist) < interfaceWidth)
        result= lInterface;         // auf dem Interface
      else if (std::abs(dist) < fadeOutWidth+interfaceWidth && signInInnerDomain*dist < 0.0) {
	double lambda = std::abs(dist)/fadeOutWidth - interfaceWidth/fadeOutWidth;
        result= static_cast<int>(std::floor(lambda*lInner+(1.0-lambda)*lInterface));
      } else if (std::abs(dist) < fadeOutWidth+interfaceWidth && signInInnerDomain*dist > 0.0) {
	double lambda = std::abs(dist)/fadeOutWidth - interfaceWidth/fadeOutWidth;
        result= static_cast<int>(std::floor(lambda*lOuter+(1.0-lambda)*lInterface));
      } else if (signInInnerDomain*dist < 0.0)
        result= lInner;
      return result;
    }
};


class ESIndicator : public MeshRefinementFunction< std::vector<double>, int >
{
public:
  ESIndicator(Mesh* mesh_, std::string name = "mesh->refinement") 
    : MeshRefinementFunction< std::vector<double>, int >(mesh_)
  {
    Parameters::get(name + "->initial level", globalSize);
  }
  
  int operator()(const std::vector<double>& means) const 
  {
    double tol = means[0];
    return (std::abs(means[1]-means[2]) > tol ? 0 : -1);
  }
    
  double indicator(const std::vector<double>& means) const 
  {
    return std::abs(means[1]-means[2]);
  }
};

} }

#endif // EXTENSIONS_MESH_FUNCTION_LEVEL_H