level_set_user_data.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2015:  G-CSC, Goethe University Frankfurt
 * Author: Christian Wehner
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __LEVEL_SET__LEVEL_SET_USER_DATA__
#define __LEVEL_SET__LEVEL_SET_USER_DATA__
 
#include "common/common.h"
 
#include "lib_grid/tools/subset_group.h"
 
#include "lib_disc/lib_disc.h"
#include "lib_disc/common/function_group.h"
#include "lib_disc/common/groups_util.h"
#include "lib_disc/quadrature/quadrature.h"
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/spatial_disc/user_data/user_data.h"
#include "lib_disc/reference_element/reference_mapping_provider.h"
#include "lib_disc/operator/non_linear_operator/newton_solver/newton_update_interface.h"
 
 
namespace ug{
namespace LevelSet{
 
template <typename TData, int dim, typename TImpl>
class LevelSetUserDataBase
:   public StdUserData<LevelSetUserDataBase<TData,dim,TImpl>, TData, dim>
{
  public:
  virtual void operator() (TData& value,
                           const MathVector<dim>& globIP,
                           number time, int si) const
  {
    UG_THROW("LevelSetUserData: Need element.");
  }
 
  virtual void operator() (TData vValue[],
                           const MathVector<dim> vGlobIP[],
                           number time, int si, const size_t nip) const
  {
    UG_THROW("LevelSetUserData: Need element.");
  }
 
  template <int refDim>
  void evaluate(TData vValue[],
                  const MathVector<dim> vGlobIP[],
                  number time, int si,
                  GridObject* elem,
                  const MathVector<dim> vCornerCoords[],
                  const MathVector<refDim> vLocIP[],
                  const size_t nip,
                  LocalVector* u,
                  const MathMatrix<refDim, dim>* vJT = NULL) const
  {
    getImpl().template evaluate<refDim>(vValue,vGlobIP,time,si,elem,
                                        vCornerCoords,vLocIP,nip,u,vJT);
  }
 
  virtual void compute(LocalVector* u, GridObject* elem,
                       const MathVector<dim> vCornerCoords[], bool bDeriv = false)
  {
    const number t = this->time();
    const int si = this->subset();
    for(size_t s = 0; s < this->num_series(); ++s)
      getImpl().template evaluate<dim>(this->values(s), this->ips(s), t, si,
                                       elem, vCornerCoords, this->template local_ips<dim>(s),
                                       this->num_ip(s), u);
  }
 
  virtual void compute(LocalVectorTimeSeries* u, GridObject* elem,
                       const MathVector<dim> vCornerCoords[], bool bDeriv = false)
  {
    const int si = this->subset();
    for(size_t s = 0; s < this->num_series(); ++s)
      getImpl().template evaluate<dim>(this->values(s), this->ips(s),  this->time(s), si,
                                       elem, vCornerCoords, this->template local_ips<dim>(s),
                                       this->num_ip(s), &(u->solution(this->time_point(s))));
  }
 
  virtual bool continuous() const {return false;}
 
  virtual bool requires_grid_fct() const {return true;}
 
  protected:
  TImpl& getImpl() {return static_cast<TImpl&>(*this);}
 
  const TImpl& getImpl() const {return static_cast<const TImpl&>(*this);}
};
 
 
 
template <typename TGridFunction>
class LevelSetUserData
: public LevelSetUserDataBase<number, TGridFunction::dim,
  LevelSetUserData<TGridFunction> >, virtual public INewtonUpdate
  {
  typedef typename TGridFunction::domain_type domain_type;
 
  typedef typename TGridFunction::algebra_type algebra_type;
 
  typedef typename domain_type::position_accessor_type position_accessor_type;
 
  static const int dim = domain_type::dim;
 
  typedef typename domain_type::grid_type grid_type;
 
  typedef typename TGridFunction::template dim_traits<dim>::grid_base_object elem_type;
 
  enum eval_type {sharp,cr_ip_average};
 
    public:
  // set evaluation type, implemented so far:
  // 0 sharp (compute lsf in ip and give out inside value if inside or outside value if outside)
  // 1 cr_ip_average (compute averaged value in ip, if in CR-FV-Geometry from-value in ip and to-value in ip have different signs,
  //                  harmonic average is computed using intersection position from from and to node which is computed from level set function,
  //                  the ips given in evaluate must be the CR-FV ips)
  void set_eval_type(int type){
    if (type==0) m_eval_type=sharp;
    if (type==1) m_eval_type=cr_ip_average;
  }
 
  eval_type m_eval_type;
 
    private:
  // level set grid function
  SmartPtr<TGridFunction> m_phi;
 
  //  approximation space for level and surface grid
  SmartPtr<ApproximationSpace<domain_type> > m_spApproxSpace;
 
  //  grid
  grid_type* m_grid;
 
  /*  //  component of function
    size_t m_fct;
 
  //  local finite element id
    LFEID m_lfeID;               */
 
    public:
  void set_inside_data(SmartPtr<CplUserData<number, dim> > user){
    m_imInsideData = user;
  }
  void set_inside_data(number val){
    set_inside_data(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_inside_data(const char* fctName){
    set_inside_data(LuaUserDataFactory<number, dim>::create(fctName));
  }
#endif
  void set_outside_data(SmartPtr<CplUserData<number, dim> > user){
    m_imOutsideData = user;
  }
  void set_outside_data(number val){
    set_outside_data(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_outside_data(const char* fctName){
    set_outside_data(LuaUserDataFactory<number, dim>::create(fctName));
  }
#endif
 
    private:
  SmartPtr<CplUserData<number,dim> > m_imInsideData;
  SmartPtr<CplUserData<number,dim> > m_imOutsideData;
 
    public:
  LevelSetUserData(SmartPtr<ApproximationSpace<domain_type> > approxSpace,SmartPtr<TGridFunction> spGridFct){
    m_phi = spGridFct;
    domain_type& domain = *m_phi->domain().get();
    grid_type& grid = *domain.grid();
    m_grid = &grid;
    m_spApproxSpace = approxSpace;
    m_eval_type = cr_ip_average;
  }
 
  virtual ~LevelSetUserData(){};
 
  template <int refDim>
  inline void evaluate(number vValue[],
                       const MathVector<dim> vGlobIP[],
                       number time, int si,
                       GridObject* elem,
                       const MathVector<dim> vCornerCoords[],
                       const MathVector<refDim> vLocIP[],
                       const size_t nip,
                       LocalVector* u,
                       const MathMatrix<refDim, dim>* vJT = NULL) const
  {
    UG_ASSERT(dynamic_cast<elem_type*>(elem) != NULL, "Unsupported element type");
    elem_type* element = static_cast<elem_type*>(elem);
 
    const size_t numVertices = element->num_vertices();
    std::vector<Vertex*> childVertex(numVertices);
    std::vector<number> phi(numVertices);
 
  //  for (size_t i=0;i<nip;i++){
  //    UG_LOG("co(" << i << ",:)=" << vGlobIP[i] << "\n");
  //  }
 
    // find child vertices by injection
    for(size_t i = 0; i < numVertices; ++i){
      childVertex[i] = element->vertex(i);
    };
 
    // find out level, then get child vertices on finest level
    size_t numChildren = m_grid->template num_children<elem_type>(element);
    if (numChildren!=0){
      size_t lowerLevel=0;
      elem_type* childElem = m_grid->template get_child<elem_type>(element,0);
      numChildren = m_grid->template num_children<elem_type>(childElem);
      lowerLevel++;
      while (numChildren>0){
        childElem = m_grid->template get_child<elem_type>(childElem,0);
        numChildren = m_grid->template num_children<elem_type>(childElem);
        lowerLevel++;
      };
      for (size_t i=0;i<nip;i++){
        for (size_t j=0;j<lowerLevel;j++){
          childVertex[i]=m_grid->template get_child<Vertex>(childVertex[i],0);
        }
      }
    }
    //    create Multiindex
    std::vector<DoFIndex> ind;
    for (size_t i=0;i<numVertices;i++){
      m_phi->dof_indices(childVertex[i], 0, ind);
      phi[i]=DoFRef(*m_phi, ind[0]);
    };
    bool onls=false;
    bool inside=false;
    for (size_t i=0;i<numVertices;i++){
      if (phi[i]==0){
        continue;
      };
      if (phi[i]<0) inside=true;
      for (size_t j=i+1;j<numVertices;j++){
        if (phi[i]*phi[j]<0){
          onls = true;
          break;
        };
      }
    }
    if (onls==false){
      if (inside==false){
        (*m_imInsideData)(vValue,
            vGlobIP,
            time, si,
            elem,
            vCornerCoords,
            vLocIP,
            nip,
            u,
            vJT);
      } else {
        (*m_imOutsideData)(vValue,
            vGlobIP,
            time, si,
            elem,
            vCornerCoords,
            vLocIP,
            nip,
            u,
            vJT);
      };
      return;
    };
    number vValueInside[max_number_of_ips];
    number vValueOutside[max_number_of_ips];
    (*m_imInsideData)(vValueInside,
        vGlobIP,
        time, si,
        elem,
        vCornerCoords,
        vLocIP,
        nip,
        u,
        vJT);
    (*m_imOutsideData)(vValueOutside,
        vGlobIP,
        time, si,
        elem,
        vCornerCoords,
        vLocIP,
        nip,
        u,
        vJT);
    if (m_eval_type == sharp){
      for(size_t ip = 0; ip < nip; ++ip)
      {
        //  reference object id
        ReferenceObjectID roid = elem->reference_object_id();
 
        //  memory for shapes
        std::vector<number> vShape;
 
        // compute lsf value in ip (from Lagrange-1 shape function) and set value according to position (inside or outside)
        const LocalShapeFunctionSet<refDim>& rTrialSpace =
            LocalFiniteElementProvider::get<refDim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
 
        //  evaluate shapes at ip
        rTrialSpace.shapes(vShape, vLocIP[ip]);
 
        //    get multiindices of element
        std::vector<DoFIndex> ind;
        m_phi->dof_indices(elem, 0, ind);
 
        //     compute lsf at integration point
        number phiValue = 0.0;
        for(size_t sh = 0; sh < vShape.size(); ++sh)
        {
          phiValue += phi[sh] * vShape[sh];
        }
        if (phiValue<0) vValue[ip] = vValueInside[ip];
        else  vValue[ip] = vValueOutside[ip];
        //                        UG_LOG("vValue(" << ip << ")=" << vValueOutside[ip] << "\n");
      };
    }
 
    if (m_eval_type == cr_ip_average){
      //  reference object id
      ReferenceObjectID roid = elem->reference_object_id();
 
      //  memory for shapes
      std::vector<number> vShape;
 
      //    get domain of grid function
      const domain_type& domain = *m_phi->domain().get();
 
      //    get position accessor
      typedef typename domain_type::position_accessor_type position_accessor_type;
      const position_accessor_type& posAcc = domain.position_accessor();
 
      //position_accessor_type aaPos = m_phi->domain()->position_accessor();
 
      //    coord and vertex array
      MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
      Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
      DimCRFVGeometry<dim> geo;
 
      for(size_t i = 0; i < numVertices; ++i){
        vVrt[i] = element->vertex(i);
        coCoord[i] = posAcc[vVrt[i]];
        // UG_LOG("co_coord(" << i<< "+1,:)=" << coCoord[i] << "\n");
      };
      //    evaluate finite volume geometry
      geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
      std::vector<number> phiSideValue(geo.num_sh());
      // compute interpolated lsf values in cr dofs
      for (size_t i=0;i<geo.num_scv();i++){
        // compute lsf value in ip (from Lagrange-1 shape function) and set value according to position (inside or outside)
        const LocalShapeFunctionSet<dim>& rTrialSpace =
            LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
 
        //  evaluate shapes at ip
        rTrialSpace.shapes(vShape, geo.scv(i).local_ip());
 
        //    get multiindices of element
        std::vector<DoFIndex> ind;
        m_phi->dof_indices(elem, 0, ind);
 
        //     compute lsf at integration point
        phiSideValue[i] = 0.0;
        for(size_t sh = 0; sh < vShape.size(); ++sh)
        {
          phiSideValue[i] += phi[sh] * vShape[sh];
        }
      }
      for (size_t ip=0;ip<nip;ip++){
        const typename DimCRFVGeometry<dim>::SCVF& scvf = geo.scvf(ip);
        number phiFrom = phiSideValue[scvf.from()];
        number phiTo = phiSideValue[scvf.to()];
        //                        UG_LOG("ip=" << ip << " inside=" << vValueInside[ip] << " outside=" << vValueOutside[ip] << "\n");
        if (phiFrom==0){
          if (phiTo<=0){
            vValue[ip] = vValueInside[ip];
          } else {
            vValue[ip] = vValueOutside[ip];
          }
        } else {
          number theta;
          if (phiFrom<0){
            if (phiTo<=0){
              vValue[ip] = vValueInside[ip];
            } else {
              theta = phiFrom/(phiFrom-phiTo);
              vValue[ip] = vValueInside[ip]*vValueOutside[ip]/(theta*vValueInside[ip]+(1-theta)*vValueOutside[ip]);
              //                                    UG_LOG("theta=" << theta << " vValue=" << vValue[ip] << "\n");
            }
          } else {
            if (phiTo>=0){
              vValue[ip] = vValueOutside[ip];
            } else {
              theta = phiFrom/(phiFrom-phiTo);
              vValue[ip] = vValueOutside[ip]*vValueInside[ip]/(theta*vValueOutside[ip]+(1-theta)*vValueInside[ip]);
              //                                    UG_LOG("theta=" << theta << " vValue=" << vValue[ip] << "\n");
            }
          }
        }
      }
    }; // ip-cr-average
  }; // evaluate
 
  void update(){}
 
    private:
  static const size_t max_number_of_ips = 20;
  };
 
template <typename TGridFunction>
class LevelSetUserVectorData
: public LevelSetUserDataBase<MathVector<TGridFunction::dim>, TGridFunction::dim,
  LevelSetUserVectorData<TGridFunction> >, virtual public INewtonUpdate
  {
  typedef typename TGridFunction::domain_type domain_type;
 
  typedef typename TGridFunction::algebra_type algebra_type;
 
  typedef typename domain_type::position_accessor_type position_accessor_type;
 
  static const int dim = domain_type::dim;
 
  typedef typename domain_type::grid_type grid_type;
 
  typedef typename TGridFunction::template dim_traits<dim>::grid_base_object elem_type;
 
  enum eval_type {sharp,cr_ip_average};
 
    public:
  // set evaluation type, implemented so far:
  // 0 sharp (compute lsf in ip and give out inside value if inside or outside value if outside)
  // 1 cr_ip_average (compute averaged value in ip, if in CR-FV-Geometry from-value in ip and to-value in ip have different signs,
  //                  harmonic average is computed using intersection position from from and to node which is computed from level set function,
  //                  the ips given in evaluate must be the CR-FV ips)
  void set_eval_type(int type){
    if (type==0) m_eval_type=sharp;
    if (type==1) m_eval_type=cr_ip_average;
  }
 
  eval_type m_eval_type;
 
    private:
  // level set grid function
  SmartPtr<TGridFunction> m_phi;
 
  //  approximation space for level and surface grid
  SmartPtr<ApproximationSpace<domain_type> > m_spApproxSpace;
 
  //  grid
  grid_type* m_grid;
 
  /*  //  component of function
      size_t m_fct;
 
    //  local finite element id
      LFEID m_lfeID;               */
 
    public:
  void set_inside_data(SmartPtr<CplUserData<MathVector<dim>, dim> > data)
  {
    m_imInsideData = data;
  }
 
  void set_inside_data(number f_x)
  {
    SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
    for (int i=0;i<dim;i++){
      f->set_entry(i, f_x);
    }
    set_inside_data(f);
  }
 
  void set_inside_data(number f_x, number f_y)
  {
    if (dim!=2){
      UG_THROW("NavierStokes: Setting source vector of dimension 2"
          " to a Discretization for world dim " << dim);
    } else {
      SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
      f->set_entry(0, f_x);
      f->set_entry(1, f_y);
      set_inside_data(f);
    }
  }
 
  void set_inside_data(number f_x, number f_y, number f_z)
  {
    if (dim<3){
      UG_THROW("NavierStokes: Setting source vector of dimension 3"
          " to a Discretization for world dim " << dim);
    }
    else
    {
      SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
      f->set_entry(0, f_x);
      f->set_entry(1, f_y);
      f->set_entry(2, f_z);
      set_inside_data(f);
    }
  }
 
#ifdef UG_FOR_LUA
  void set_inside_data(const char* fctName)
  {
    set_inside_data(LuaUserDataFactory<MathVector<dim>, dim>::create(fctName));
  }
#endif
 
  void set_outside_data(SmartPtr<CplUserData<MathVector<dim>, dim> > data)
  {
    m_imOutsideData = data;
  }
 
  void set_outside_data(number f_x)
  {
    SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
    for (int i=0;i<dim;i++){
      f->set_entry(i, f_x);
    }
    set_outside_data(f);
  }
 
  void set_outside_data(number f_x, number f_y)
  {
    if (dim!=2){
      UG_THROW("NavierStokes: Setting source vector of dimension 2"
          " to a Discretization for world dim " << dim);
    } else {
      SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
      f->set_entry(0, f_x);
      f->set_entry(1, f_y);
      set_outside_data(f);
    }
  }
 
  void set_outside_data(number f_x, number f_y, number f_z)
  {
    if (dim<3){
      UG_THROW("NavierStokes: Setting source vector of dimension 3"
          " to a Discretization for world dim " << dim);
    }
    else
    {
      SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
      f->set_entry(0, f_x);
      f->set_entry(1, f_y);
      f->set_entry(2, f_z);
      set_outside_data(f);
    }
  }
 
#ifdef UG_FOR_LUA
  void set_outside_data(const char* fctName)
  {
    set_outside_data(LuaUserDataFactory<MathVector<dim>, dim>::create(fctName));
  }
#endif
 
    private:
  SmartPtr<CplUserData<MathVector<dim>,dim> > m_imInsideData;
  SmartPtr<CplUserData<MathVector<dim>,dim> > m_imOutsideData;
 
    public:
  LevelSetUserVectorData(SmartPtr<ApproximationSpace<domain_type> > approxSpace,SmartPtr<TGridFunction> spGridFct){
    m_phi = spGridFct;
    domain_type& domain = *m_phi->domain().get();
    grid_type& grid = *domain.grid();
    m_grid = &grid;
    m_spApproxSpace = approxSpace;
    m_eval_type = cr_ip_average;
  }
 
  virtual ~LevelSetUserVectorData(){};
 
  template <int refDim>
  inline void evaluate(MathVector<dim> vValue[],
                       const MathVector<dim> vGlobIP[],
                       number time, int si,
                       GridObject* elem,
                       const MathVector<dim> vCornerCoords[],
                       const MathVector<refDim> vLocIP[],
                       const size_t nip,
                       LocalVector* u,
                       const MathMatrix<refDim, dim>* vJT = NULL) const
  {
    UG_ASSERT(dynamic_cast<elem_type*>(elem) != NULL, "Unsupported element type");
    elem_type* element = static_cast<elem_type*>(elem);
 
    const size_t numVertices = element->num_vertices();
    std::vector<Vertex*> childVertex(numVertices);
    std::vector<number> phi(numVertices);
 
    for (size_t i=0;i<nip;i++){
      //                  UG_LOG("co(" << i << ",:)=" << vGlobIP[i] << "\n");
    }
 
    // find child vertices by injection
    for(size_t i = 0; i < numVertices; ++i){
      childVertex[i] = element->vertex(i);
    };
 
    // find out level, then get child vertices on finest level
    size_t numChildren = m_grid->template num_children<elem_type>(element);
    if (numChildren!=0){
      size_t lowerLevel=0;
      elem_type* childElem = m_grid->template get_child<elem_type>(element,0);
      numChildren = m_grid->template num_children<elem_type>(childElem);
      lowerLevel++;
      while (numChildren>0){
        childElem = m_grid->template get_child<elem_type>(childElem,0);
        numChildren = m_grid->template num_children<elem_type>(childElem);
        lowerLevel++;
      };
      for (size_t i=0;i<nip;i++){
        for (size_t j=0;j<lowerLevel;j++){
          childVertex[i]=m_grid->template get_child<Vertex>(childVertex[i],0);
        }
      }
    }
    //    create Multiindex
    std::vector<DoFIndex> ind;
    for (size_t i=0;i<numVertices;i++){
      m_phi->dof_indices(childVertex[i], 0, ind);
      phi[i]=DoFRef(*m_phi, ind[0]);
    };
    bool onls=false;
    bool inside=false;
    for (size_t i=0;i<numVertices;i++){
      if (phi[i]==0){
        continue;
      };
      if (phi[i]<0) inside=true;
      for (size_t j=i+1;j<numVertices;j++){
        if (phi[i]*phi[j]<0){
          onls = true;
          break;
        };
      }
    }
    if (onls==false){
      if (inside==false){
        (*m_imInsideData)(vValue,
            vGlobIP,
            time, si,
            elem,
            vCornerCoords,
            vLocIP,
            nip,
            u,
            vJT);
      } else {
        (*m_imOutsideData)(vValue,
            vGlobIP,
            time, si,
            elem,
            vCornerCoords,
            vLocIP,
            nip,
            u,
            vJT);
      };
      return;
    };
    MathVector<dim>  vValueInside[max_number_of_ips];
    MathVector<dim>  vValueOutside[max_number_of_ips];
    (*m_imInsideData)(vValueInside,
        vGlobIP,
        time, si,
        elem,
        vCornerCoords,
        vLocIP,
        nip,
        u,
        vJT);
    (*m_imOutsideData)(vValueOutside,
        vGlobIP,
        time, si,
        elem,
        vCornerCoords,
        vLocIP,
        nip,
        u,
        vJT);
    if (m_eval_type == sharp){
      for(size_t ip = 0; ip < nip; ++ip)
      {
        //  reference object id
        ReferenceObjectID roid = elem->reference_object_id();
 
        //  memory for shapes
        std::vector<number> vShape;
 
        // compute lsf value in ip (from Lagrange-1 shape function) and set value according to position (inside or outside)
        const LocalShapeFunctionSet<refDim>& rTrialSpace =
            LocalFiniteElementProvider::get<refDim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
 
        //  evaluate shapes at ip
        rTrialSpace.shapes(vShape, vLocIP[ip]);
 
        //    get multiindices of element
        std::vector<DoFIndex> ind;
        m_phi->dof_indices(elem, 0, ind);
 
        //     compute lsf at integration point
        number phiValue = 0.0;
        for(size_t sh = 0; sh < vShape.size(); ++sh)
        {
          phiValue += phi[sh] * vShape[sh];
        }
        if (phiValue<0) vValue[ip] = vValueInside[ip];
        else  vValue[ip] = vValueOutside[ip];
        //                        UG_LOG("vValue(" << ip << ")=" << vValueOutside[ip] << "\n");
      };
    }
 
    if (m_eval_type == cr_ip_average){
      //  reference object id
      ReferenceObjectID roid = elem->reference_object_id();
 
      //  memory for shapes
      std::vector<number> vShape;
 
      //    get domain of grid function
      const domain_type& domain = *m_phi->domain().get();
 
      //    get position accessor
      typedef typename domain_type::position_accessor_type position_accessor_type;
      const position_accessor_type& posAcc = domain.position_accessor();
 
      //position_accessor_type aaPos = m_phi->domain()->position_accessor();
 
      //    coord and vertex array
      MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
      Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
      DimCRFVGeometry<dim> geo;
 
      for(size_t i = 0; i < numVertices; ++i){
        vVrt[i] = element->vertex(i);
        coCoord[i] = posAcc[vVrt[i]];
        // UG_LOG("co_coord(" << i<< "+1,:)=" << coCoord[i] << "\n");
      };
      //    evaluate finite volume geometry
      geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
      std::vector<number> phiSideValue(geo.num_sh());
      // compute interpolated lsf values in cr dofs
      for (size_t i=0;i<geo.num_scv();i++){
        // compute lsf value in ip (from Lagrange-1 shape function) and set value according to position (inside or outside)
        const LocalShapeFunctionSet<dim>& rTrialSpace =
            LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
 
        //  evaluate shapes at ip
        rTrialSpace.shapes(vShape, geo.scv(i).local_ip());
 
        //    get multiindices of element
        std::vector<DoFIndex> ind;
        m_phi->dof_indices(elem, 0, ind);
 
        //     compute lsf at integration point
        phiSideValue[i] = 0.0;
        for(size_t sh = 0; sh < vShape.size(); ++sh)
        {
          phiSideValue[i] += phi[sh] * vShape[sh];
        }
      }
      for (size_t ip=0;ip<nip;ip++){
        const typename DimCRFVGeometry<dim>::SCVF& scvf = geo.scvf(ip);
        number phiFrom = phiSideValue[scvf.from()];
        number phiTo = phiSideValue[scvf.to()];
        //                        UG_LOG("ip=" << ip << " inside=" << vValueInside[ip] << " outside=" << vValueOutside[ip] << "\n");
        if (phiFrom==0){
          if (phiTo<=0){
            vValue[ip] = vValueInside[ip];
          } else {
            vValue[ip] = vValueOutside[ip];
          }
        } else {
          number theta;
          if (phiFrom<0){
            if (phiTo<=0){
              vValue[ip] = vValueInside[ip];
            } else {
              theta = phiFrom/(phiFrom-phiTo);
              for (int d=0;d<dim;d++)
                vValue[ip][d] = vValueInside[ip][d]*vValueOutside[ip][d]/(theta*vValueInside[ip][d]+(1-theta)*vValueOutside[ip][d]);
              //                                    UG_LOG("theta=" << theta << " vValue=" << vValue[ip] << "\n");
            }
          } else {
            if (phiTo>=0){
              vValue[ip] = vValueOutside[ip];
            } else {
              theta = phiFrom/(phiFrom-phiTo);
              for (int d=0;d<dim;d++)
                vValue[ip][d] = vValueOutside[ip][d]*vValueInside[ip][d]/(theta*vValueOutside[ip][d]+(1-theta)*vValueInside[ip][d]);
              //                                    UG_LOG("theta=" << theta << " vValue=" << vValue[ip] << "\n");
            }
          }
        }
      }
    }; // ip-cr-average
  }; // evaluate
 
  void update(){}
 
    private:
  static const size_t max_number_of_ips = 20;
  };
 
 
template<typename TElem,typename TGrid>
void collect_finest_level_children(std::vector<TElem*>& childElemVector,TElem* elem,TGrid* grid){
  size_t numChildren = grid->template num_children<TElem>(elem);
  if (numChildren==0){
    childElemVector.push_back(elem);
  } else {
    for (size_t i=0;i<numChildren;i++){
      TElem* childElem = grid->template get_child<TElem>(elem,i);
      size_t nc=grid->template num_children<TElem>(childElem);
      if (nc==0){
        childElemVector.push_back(childElem);
      } else {
        collect_finest_level_children<TElem,TGrid>(childElemVector,childElem,grid);
      };
    };
  }
};
 
template<typename TElem,typename TPos,int dim>
void computeElemBarycenter(MathVector<dim>& bary,TElem* elem,TPos posA){
  bary = 0;
  size_t noc=elem->num_vertices();
  for (size_t i=0;i<noc;i++){
    bary+=posA[elem->vertex(i)];
  }
  bary/=noc;
}
 
template <typename TGridFunction>
class CRTwoPhaseSource
: public LevelSetUserDataBase<MathVector<TGridFunction::dim>, TGridFunction::dim,
  CRTwoPhaseSource<TGridFunction> >, virtual public INewtonUpdate
  {
  typedef typename TGridFunction::domain_type domain_type;
 
  typedef typename TGridFunction::algebra_type algebra_type;
 
  typedef typename domain_type::position_accessor_type position_accessor_type;
 
  static const int dim = domain_type::dim;
 
  typedef typename domain_type::grid_type grid_type;
 
  typedef typename TGridFunction::template dim_traits<dim>::grid_base_object elem_type;
 
    private:
  // level set grid function
  SmartPtr<TGridFunction> m_phi;
 
  //  approximation space for level and surface grid
  SmartPtr<ApproximationSpace<domain_type> > m_spApproxSpace;
 
  //  grid
  grid_type* m_grid;
 
    private:
 
  // gravitational constant
  number m_gravitational_constant;
 
  // surface tension factor
  number m_sigma;
 
  SmartPtr<CplUserData<MathVector<dim>,dim> > m_imSource;
 
    public:
  void set_gravitation(number gravityconst){
    m_gravitational_constant = gravityconst;
  }
 
  void set_sigma(number sigma){
    m_sigma = sigma;
  }
 
 
  void set_source(SmartPtr<CplUserData<MathVector<dim>, dim> > data)
  {
    m_imSource = data;
  }
 
  void set_source(number f_x)
  {
    SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
    for (int i=0;i<dim;i++){
      f->set_entry(i, f_x);
    }
    set_source(f);
  }
 
  void set_source(number f_x, number f_y)
  {
    if (dim!=2){
      UG_THROW("NavierStokes: Setting source vector of dimension 2"
          " to a Discretization for world dim " << dim);
    } else {
      SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
      f->set_entry(0, f_x);
      f->set_entry(1, f_y);
      set_source(f);
    }
  }
 
  void set_source(number f_x, number f_y, number f_z)
  {
    if (dim<3){
      UG_THROW("NavierStokes: Setting source vector of dimension 3"
          " to a Discretization for world dim " << dim);
    }
    else
    {
      SmartPtr<ConstUserVector<dim> > f(new ConstUserVector<dim>());
      f->set_entry(0, f_x);
      f->set_entry(1, f_y);
      f->set_entry(2, f_z);
      set_source(f);
    }
  }
 
#ifdef UG_FOR_LUA
  void set_source(const char* fctName)
  {
    set_source(LuaUserDataFactory<MathVector<dim>, dim>::create(fctName));
  }
#endif
 
    public:
  void set_density(SmartPtr<CplUserData<number, dim> > user){
    m_imDensity = user;
  }
  void set_density(number val){
    set_density(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_density(const char* fctName){
    set_density(LuaUserDataFactory<number, dim>::create(fctName));
  }
#endif
 
    private:
  SmartPtr<CplUserData<number,dim> > m_imDensity;
 
    public:
  CRTwoPhaseSource(SmartPtr<ApproximationSpace<domain_type> > approxSpace,SmartPtr<TGridFunction> spGridFct){
    m_phi = spGridFct;
    for (int si=0;si<m_phi->num_subsets();++si){
      if (m_phi->num_fct(si)<2){
        UG_THROW("No curvature component in approximation space.");
      }
      if (m_phi->local_finite_element_id(0) != LFEID(LFEID::LAGRANGE, dim,1)){
        UG_THROW("First component in approximation space must be of Lagrange 1 type.");
      }
      if (m_phi->local_finite_element_id(1) != LFEID(LFEID::PIECEWISE_CONSTANT,dim,0)){
        UG_THROW("Second component in approximation space must be of piecewise constant type.");
      }
    };
    domain_type& domain = *m_phi->domain().get();
    grid_type& grid = *domain.grid();
    m_grid = &grid;
    m_spApproxSpace = approxSpace;
    set_gravitation(980);
    set_source(0.0);
    set_density(1);
  }
 
  virtual ~CRTwoPhaseSource(){};
 
  template <int refDim>
  inline void evaluate(MathVector<dim> vValue[],
                       const MathVector<dim> vGlobIP[],
                       number time, int si,
                       GridObject* elem,
                       const MathVector<dim> vCornerCoords[],
                       const MathVector<refDim> vLocIP[],
                       const size_t nip,
                       LocalVector* u,
                       const MathMatrix<refDim, dim>* vJT = NULL) const
  {
    // evaluate source data
    //           for (size_t i=0;i<nip;i++) vValue[i]=0;
    (*m_imSource)(vValue,
        vGlobIP,
        time, si,
        elem,
        vCornerCoords,
        vLocIP,
        nip,
        u,
        vJT);
    // find corner values of level set function on finest level
    UG_ASSERT(dynamic_cast<elem_type*>(elem) != NULL, "Unsupported element type");
    elem_type* element = static_cast<elem_type*>(elem);
 
    const size_t numVertices = element->num_vertices();
    std::vector<Vertex*> childVertex(numVertices);
    std::vector<number> phi(numVertices);
 
    // find child vertices by injection
    for(size_t i = 0; i < numVertices; ++i){
      childVertex[i] = element->vertex(i);
    };
 
    bool onfinestlevel;
    // find out level, then get child vertices on finest level
    size_t numChildren = m_grid->template num_children<elem_type>(element);
    if (numChildren!=0){
      size_t lowerLevel=0;
      elem_type* childElem = m_grid->template get_child<elem_type>(element,0);
      numChildren = m_grid->template num_children<elem_type>(childElem);
      lowerLevel++;
      while (numChildren>0){
        childElem = m_grid->template get_child<elem_type>(childElem,0);
        numChildren = m_grid->template num_children<elem_type>(childElem);
        lowerLevel++;
      };
      for (size_t i=0;i<nip;i++){
        for (size_t j=0;j<lowerLevel;j++){
          childVertex[i]=m_grid->template get_child<Vertex>(childVertex[i],0);
        }
      }
    } else onfinestlevel=true;
 
    //    create Multiindex
    std::vector<DoFIndex> ind;
    for (size_t i=0;i<numVertices;i++){
      m_phi->dof_indices(childVertex[i], 0, ind);
      phi[i]=DoFRef(*m_phi, ind[0]);
    };
    bool onls=false;
    //bool inside=false; // inside seems to be unused (a.vogel)
    for (size_t i=0;i<numVertices;i++){
      //debug UG_LOG("phi(" << i+1 << ")=" << phi[i] << "\n");
      if (phi[i]==0){
        continue;
      };
      //if (phi[i]<0) inside=true;
      for (size_t j=i+1;j<numVertices;j++){
        //debug UG_LOG(" phi(" << j+1 << ")=" << phi[j] << "\n");
        if (phi[i]*phi[j]<0){
          onls = true;
          break;
        };
      }
      if (onls==true) break;
    }
    if (m_gravitational_constant!=0){
      number densityValue[max_number_of_ips];
      (*m_imDensity)(densityValue,
          vGlobIP,
          time, si,
          elem,
          vCornerCoords,
          vLocIP,
          nip,
          u,
          vJT);
      for (int d=0;d<dim;d++){
        for (size_t i=0;i<nip;i++){
          vValue[i][d] += m_gravitational_constant * densityValue[i];
        }
      }
    }
    if (onls==true){
      //  reference object id
      ReferenceObjectID roid = elem->reference_object_id();
 
      //  memory for shapes
      std::vector<number> vShape;
 
      //    get domain of grid function
      const domain_type& domain = *m_phi->domain().get();
 
      //    get position accessor
      typedef typename domain_type::position_accessor_type position_accessor_type;
      const position_accessor_type& posAcc = domain.position_accessor();
 
      position_accessor_type aaPos = m_phi->domain()->position_accessor();
 
      //    coord and vertex array
      MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
      Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
      elem_type* evaluationElem;
      MathVector<dim> surftensSource[max_number_of_ips];
      for (size_t i=0;i<nip;i++){
        surftensSource[i]=0;
      }
      if (onfinestlevel==true){
        evaluationElem = element;
      } else {
        // compute barycenter
        MathVector<dim> coarseElementBary;
        computeElemBarycenter<elem_type,position_accessor_type,dim>(coarseElementBary,element,aaPos);
        // collect child elements
        std::vector<elem_type*> childElems;
        collect_finest_level_children<elem_type,grid_type>(childElems,element,m_grid);
        evaluationElem = childElems[0];
        number dist=1e+9;
        for (size_t ce=0;ce<childElems.size();ce++){
          MathVector<dim> bary;
          elem_type* currentElem = childElems[ce];
          size_t elnoc = currentElem->num_vertices();
          bool currentElemOnLS=false;
          std::vector<number> currentPhi(elnoc);
          for (size_t j=0;j<elnoc;j++){
            m_phi->dof_indices(currentElem->vertex(j),0,ind);
            currentPhi[j]=DoFRef(*m_phi,ind[0]);
          }
          for (size_t i=0;i<numVertices;i++){
            if (currentPhi[i]==0){
              continue;
            };
            for (size_t j=i+1;j<numVertices;j++){
              if (currentPhi[i]*currentPhi[j]<0){
                currentElemOnLS = true;
                break;
              }
            }
            if (currentElemOnLS==true) break;
          }
          if (currentElemOnLS==false) continue;
          computeElemBarycenter<elem_type,position_accessor_type,dim>(bary,currentElem,aaPos);
          number localdist  = VecDistance(bary,coarseElementBary);
          if (localdist<dist){
            dist = localdist;
            evaluationElem = childElems[ce];
          }
        }
      }
      m_phi->dof_indices(evaluationElem,1, ind);
      // evaluate curvature
      number element_curvature=DoFRef(*m_phi, ind[0]);
      // UG_LOG(element_curvature << "\n");
      // scale with surface tension factor
      element_curvature*=m_sigma;
      DimCRFVGeometry<dim> geo;
      const size_t eeNumVertices = evaluationElem->num_vertices();
      for(size_t i = 0; i < eeNumVertices; ++i){
        vVrt[i] = evaluationElem->vertex(i);
        coCoord[i] = posAcc[vVrt[i]];
        // UG_LOG("co(" << i+1 << ",:)=[" << coCoord[i][0] << "," << coCoord[i][1] << "];\n");
      };
      // evaluate finite volume geometry
      geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
      std::vector<number> phiSideValue(geo.num_sh());
      // compute interpolated lsf values in cr dofs
      for (size_t i=0;i<geo.num_scv();i++){
        // compute lsf value in ip (from Lagrange-1 shape function) and set value according to position (inside or outside)
        const LocalShapeFunctionSet<dim>& rTrialSpace =
            LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim,1));
 
        //  evaluate shapes at ip
        rTrialSpace.shapes(vShape, geo.scv(i).local_ip());
 
        //    get multiindices of element
        std::vector<DoFIndex> ind;
        m_phi->dof_indices(elem, 0, ind);
 
        //     compute lsf at integration point
        phiSideValue[i] = 0.0;
        for(size_t sh = 0; sh < vShape.size(); ++sh)
        {
          phiSideValue[i] += phi[sh] * vShape[sh];
        }
      }
      //for (size_t i=0;i<geo.num_scv();i++){//debug
      //debug  UG_LOG("phi(" << i << ")=" << phiSideValue[i] << "\n");//debug
      //}//debug
      size_t nscvf = geo.num_scvf();
      for (size_t ip=0;ip<nscvf;ip++){
        const typename DimCRFVGeometry<dim>::SCVF& scvf = geo.scvf(ip);
        number phiFrom = phiSideValue[scvf.from()];
        number phiTo = phiSideValue[scvf.to()];
        if (phiFrom<0)   for (int d=0;d<dim;d++) surftensSource[scvf.from()][d] += scvf.normal()[d] * element_curvature;
        if (phiTo<0) for (int d=0;d<dim;d++) surftensSource[scvf.to()][d] -= scvf.normal()[d] * element_curvature;
      }
      for (size_t i=0;i<geo.num_scv();i++){
        surftensSource[i]/=geo.scv(i).volume();
        vValue[i] += surftensSource[i];
      }
      //for (size_t i=0;i<nip;i++){
        //  UG_LOG("rhs(" << i << ")=" << vValue[i] << "\n");//debug
      //};
      //UG_LOG("##############\n");//debug
    }
  }; // evaluate
 
  void update(){}
 
    private:
  static const size_t max_number_of_ips = 20;
  };
 
 
}; // end namespace levelset
}; // end namespace ug
 
#endif /* __LEVEL_SET__LEVEL_SET_USER_DATA__ */