fract_dimadapt_impl.h

Go to the documentation of this file.

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Sabine Stichel
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
 
/*
 * Implementation of the dimension-adaptive fracture representation
 */
 
// ug4 headers
#include "lib_disc/domain_traits.h"
#include"lib_grid/grid/grid.h"
#include "lib_grid/lg_base.h"
 
//  Only required for debug saves
#include "lib_grid/file_io/file_io.h"
 
namespace ug {
namespace d3f {
 
template <typename TDomain,typename TAlgebra>
FractDimadapt<TDomain,TAlgebra>::FractDimadapt
(
  SmartPtr<domain_type> domain, 
  const char* ss_name, 
  const std::vector<number> normal, 
  number aperture, 
  bool shift  
)
{
// Check the domain:
  if (domain.invalid ())
    UG_THROW ("FractDimadapt: Invalid domain specified!");
 
  m_dom = domain;
 
  m_shift = shift;
 
//  Initialize the data:
  m_spSH = domain->subset_handler ();
  m_fractSsGrp.set_subset_handler (m_spSH);
  
//initialize
  m_FirstFullLevel = -1;
  m_FractureIsFullDimensional = -1;
  m_density_type = "NONE";
  m_viscosity_type = "NONE";
  m_permeability_f = 0.0;
  m_permeability_m = 0.0;
//set default gravity
  m_Gravity[0] = 0.0;
  m_Gravity[1] = 0.0;
  m_Gravity[dim-1] = -9.81;
 
//  Initialize the attachment:
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
  pMG->template attach_to_dv<Vertex> (m_aVertexAuxMarks, F_NONE);
  m_aaVertAuxMarks.access (*pMG, m_aVertexAuxMarks);
  
  pMG->template attach_to_dv<Edge> (m_aEdgeAuxMarks, F_NONE);
  m_aaEdgeAuxMarks.access (*pMG, m_aEdgeAuxMarks);
 
  MathVector<dim> initial;
  for(int i=0;i<dim;i++)
    initial[i]=0.0;
  pMG->template attach_to_dv<Vertex> (m_aVertexShift, initial);
  m_aaVertexShift.access (*pMG, m_aVertexShift);
 
//  Register the callbacks in the message hub:
  m_spGridAdaptionCallbackID =
    pMG->message_hub()->register_class_callback
      (this, & FractDimadapt<TDomain,TAlgebra>::grid_adaption_callback);
  m_spGridDistributionCallbackID =
    pMG->message_hub()->register_class_callback
      (this, & FractDimadapt<TDomain,TAlgebra>::grid_distribution_callback);
 
  //  Add the subsets to the group:
  m_fractSsGrp.add (ss_name);
  for (int d=0;d < dim;d++)
    m_fractNormal[d] = normal[d];
  m_aperture = aperture;
 
  //  Mark the vertices:
  mark_vertices ();
}
 
template <typename TDomain,typename TAlgebra>
FractDimadapt<TDomain,TAlgebra>::~FractDimadapt ()
{
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
  m_aaVertAuxMarks.invalidate ();
  pMG->template detach_from<Vertex> (m_aVertexAuxMarks);
  m_aaEdgeAuxMarks.invalidate ();
  pMG->template detach_from<side_type> (m_aEdgeAuxMarks);
}
 
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::grid_adaption_callback
(
  const GridMessage_Adaption & msg 
)
{
  if (msg.adaption_ends ())
    mark_vertices ();
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::grid_distribution_callback
(
  const GridMessage_Distribution & msg 
)
{
  if (msg.msg () == GMDT_DISTRIBUTION_STOPS)
    mark_vertices ();
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::mark_vertices ()
{
  typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
 
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
  m_NumberLevels = pMG->num_levels ();
  if(m_NumberLevels==(int)max_number_levels)
    UG_THROW ("FractDimadapt: Not enough Levels allocated!");
  
  typename Grid::traits<side_type>::secure_container side_list;
  typename Grid::traits<Edge>::secure_container edge_list;
 
  const position_accessor_type& posAcc = m_dom->position_accessor();
 
  number eps = 1E-10;
 
  UG_LOG("mark vertices\n");
 
 
  for (int lev = 0; lev < m_NumberLevels; lev++)
  {
    //First find nodes on fracture mean plane
    //To this end. Mark Fracture Vertices with F_INNER
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
    int si = m_fractSsGrp [i];
 
    t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
    for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != list_end; ++iElem)
    {
      element_type * elem = *iElem;
      for (size_t j = 0; j < elem->num_vertices (); j++)
      {
        Vertex * vert = elem->vertex(j);
        m_aaVertAuxMarks [vert] = (! vert->is_constrained ())? F_INNER : F_OUTER;
          // REMARK: Constrained vertices cannot be inner for fractures!
          // We cannot catch this situation somewhere else, so we do it here.
      }
    }
    }
 
    //  Unmark all vertices in the bulk medium elements:
    SubsetGroup bulkSsGrp (m_spSH);
    for (int si = 0; si < m_spSH->num_subsets (); si++)
      if (DimensionOfSubset (*m_spSH, si) == dim && ! m_fractSsGrp.contains (si))
      {
        t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
        for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
            iElem != list_end; ++iElem)
        {
          element_type * elem = *iElem;
          for (size_t i = 0; i < elem->num_vertices (); i++)
            m_aaVertAuxMarks [elem->vertex(i)] = F_OUTER;
        }
      }
  }
 
  //Then add marked nodes to list
  //First find inner nodes on level 0 (on lev 0 all inner are mid)
  NumberFractNodes[0] = 0.0;
  for (size_t i = 0; i < m_fractSsGrp.size (); i++)
  {
    int si = m_fractSsGrp [i];
    t_elem_iter list_end = m_spSH->template end<element_type> (si, 0);
    for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, 0);
          iElem != list_end; ++iElem)
    {
      element_type * elem = *iElem;
      for (size_t j = 0; j < elem->num_vertices (); j++)
      {
        Vertex * vert = elem->vertex(j);
        int mark = m_aaVertAuxMarks [vert];
        if(mark == F_INNER)
        {
          FractNode[0][NumberFractNodes[0]][0] = vert;
          NumberFractNodes[0]++;
          if(NumberFractNodes[0]==(int)max_number_fracture_midnodes)
            UG_THROW ("FractDimadapt: Not enough Fracture nodes allocated!");
          m_aaVertAuxMarks [vert] = F_MID; //to avoid double counting
        }
      }
    }
  }
  //then on other levels find inner nodes that lie on fracture mean plane
  //given by line through mid with normal
  number c = VecDot(posAcc[FractNode[0][0][0]],m_fractNormal);
  for (int lev = 1; lev < m_NumberLevels; lev++)
  {
    NumberFractNodes[lev] = 0.0;
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
      int si = m_fractSsGrp [i];
      t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
            iElem != list_end; ++iElem)
      {
        element_type * elem = *iElem;
        for (size_t j = 0; j < elem->num_vertices (); j++)
        {
          Vertex * vert = elem->vertex(j);
          int mark = m_aaVertAuxMarks [vert];
          if(mark == F_INNER)
          {
            number c_a = VecDot(posAcc[vert],m_fractNormal);
            if(fabs(c-c_a)<eps)
            {
              FractNode[lev][NumberFractNodes[lev]][0] = vert;
              NumberFractVertNodes[lev][NumberFractNodes[lev]] = 0.0;
              NumberFractNodes[lev]++;
              if(NumberFractNodes[lev]==(int)max_number_fracture_midnodes)
                UG_THROW ("FractDimadapt: Not enough Fracture nodes allocated!");
              m_aaVertAuxMarks [vert] = F_MID; //to avoid double counting
            }
          }
        }
      }
    }
  }
 
  //now find corresponding nodes along the line given by normal
  //(normal to the fracture is parallel to this line)
  MathVector<dim> tangential;
  tangential[0]= m_fractNormal[1];
  tangential[1]= -m_fractNormal[0];
  //todo: 3d tangential berechnen!!
  for (int lev = 0; lev < m_NumberLevels; lev++)
  {
    int TriangleCount = 0;
    std::vector<int> bnd_count(NumberFractNodes[lev]);
    for (int k = 0; k < NumberFractNodes[lev]; k++)
      bnd_count[k] = 0;
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
      int si = m_fractSsGrp [i];
      t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != list_end; ++iElem)
      {
        element_type * elem = *iElem;
        for (size_t j = 0; j < elem->num_vertices (); j++)
        {
          Vertex * vert = elem->vertex(j);
          int mark = m_aaVertAuxMarks [vert];
          if(mark == F_MID || mark == F_BND || mark == F_VERT)
            continue;
          if(mark == F_NONE)
            UG_THROW ("FractDimadapt:Fracture marks not set correctly!");
          //UG_LOG("\t teste " <<  posAcc[vert]<<" "<<mark<<"\n");
          number c_a = VecDot(posAcc[vert],tangential);
          number c_b;
          int corner = 1;
          for (int k = 0; k < NumberFractNodes[lev]; k++)
          {
            Vertex * midvert = FractNode[lev][k][0];
            //UG_LOG("\t mid " <<  posAcc[midvert]<<"\n");
            c_b = VecDot(posAcc[midvert],tangential);
            if(fabs(c_b-c_a)<eps)
            {
              if(mark == F_OUTER)
              {
                //UG_LOG("\t OUTER "<<"\n");
                FractNode[lev][k][bnd_count[k]+1] = vert;
                bnd_count[k]++;
                m_aaVertAuxMarks [vert] = F_BND; //to avoid double counting
                k=NumberFractNodes[lev];
              }
              else if(mark == F_INNER)
              {
                //UG_LOG("\t INNER "<<"\n");
                FractNode[lev][k][3+NumberFractVertNodes[lev][k]] = vert;
                NumberFractVertNodes[lev][k]++;
                if(NumberFractVertNodes[lev][k]== (int)max_number_fracture_layers)
                  UG_THROW ("FractDimadapt: Not enough Fracture vert nodes allocated!");
                m_FractureIsFullDimensional = 1;
                m_aaVertAuxMarks [vert] = F_VERT; //to avoid double counting
                k=NumberFractNodes[lev];
                if(m_FirstFullLevel < 0)
                {
                  m_FirstFullLevel = lev;
                }
              }
              else
                UG_THROW ("FractDimadapt:Fracture marks maybe not set correctly!");
              corner = 0;
            }
          }
          if(corner == 1)
          {
            m_aaVertAuxMarks [vert] = F_CORNER;
            CornerTriangleList[lev][TriangleCount] = elem;
            TriangleCount++;
            //UG_LOG("\t CORNER "<<" "<<posAcc[vert]<<"\n");
            if(TriangleCount > 4)
              UG_THROW ("FractDimadapt: Something wrong with corner triangle identification"<<TriangleCount);
 
          }
        }
      }
//      if(TriangleCount != 4) //only for full dim
//        UG_THROW ("FractDimadapt: Something wrong with corner triangle identification"<<TriangleCount);
    }
 
    int deg_count = 0; //for degenerated quadrilaterals that are triangles in full dim
    for (int k = 0; k < NumberFractNodes[lev]; k++)
    {
      if(bnd_count[k]!=2)
      {
        if(m_shift == true)
        {
          if(bnd_count[k]==1) //degenerated end quadri sets corner node as bnd node to triquadnode
            {
              deg_count++;
              FractNode[lev][k][1] = NULL;
              FractNode[lev][k][2] = NULL;
            }
          //else //in debug case possible
          //  UG_THROW ("FractDimadapt: Wrong number fract boundary nodes!"<<bnd_count[k]<<" "<<posAcc[FractNode[lev][k][0]]);
        }
        else
        {
          if(bnd_count[k]==0) //degenerated end quadri without shift: quadtrinode has no bnd nodes
            deg_count++;
          else
            UG_THROW ("FractDimadapt: Wrong number fract boundary nodes!"<<bnd_count[k]);
        }
      }
    }
    if(deg_count!=0 && deg_count!=2)
      UG_THROW ("FractDimadapt: Wrong number fract boundary nodes (degcount)!");
  }
 
  //set FractureIsFullDimensional also for cases without F_VERT nodes in fracture
  //i.e. for full-dim fractures with only 2 layers
  //check only if FractureIsFullDimensional is not set before
  if(m_FractureIsFullDimensional < 0)
  {
    //check shift on toplevel
    int lev = m_NumberLevels-1;
    //check only for one point
    int k = 1;
    MathVector<dim> coords_1,coords_2,shift;
    coords_1 = posAcc[FractNode[lev][k][0]];
    coords_2 = posAcc[FractNode[lev][k][1]];
    VecSubtract(shift,coords_1,coords_2);
    if(fabs(VecLength(shift))>1e-6)
    {
      m_FractureIsFullDimensional = 1;
      m_FirstFullLevel = m_NumberLevels;
    }
  }
 
  //set shift of boundary nodes (can only be set for full dim fracture)
  if(m_FractureIsFullDimensional == 1)
  {
    for (int lev = 0; lev < m_NumberLevels; lev++)
    {
      for (int k = 0; k < NumberFractNodes[lev]; k++)
      {
        MathVector<dim> coords_1,coords_2,shift;
        coords_1 = posAcc[FractNode[lev][k][0]];
        coords_2 = posAcc[FractNode[lev][k][1]];
        VecSubtract(shift,coords_1,coords_2);
        m_aaVertexShift [FractNode[lev][k][1]] = shift;
        coords_2 = posAcc[FractNode[lev][k][2]];
        VecSubtract(shift,coords_1,coords_2);
        m_aaVertexShift [FractNode[lev][k][2]] = shift;
      }
    }
  }
 
 
  //print list
//  for (int lev = 0; lev < m_NumberLevels; lev++)
//  {
//    UG_LOG(" level "<<lev<<" "<<NumberFractNodes[lev] <<" \n");
//    for (int k = 0; k < NumberFractNodes[lev]; k++)
//    {
//      UG_LOG("vert nodes "<<NumberFractVertNodes[lev][k]<<"\n ");
//      UG_LOG("\t vert " <<  posAcc[FractNode[lev][k][0]] <<" "<< posAcc[FractNode[lev][k][1]] <<" "<<posAcc[FractNode[lev][k][2]] );
//      UG_LOG("\t shift " << vert_shift(FractNode[lev][k][1]) <<" "<<vert_shift(FractNode[lev][k][2]) );
//      for (int i = 3; i < NumberFractVertNodes[lev][k]+3; i++)
//      {
//        UG_LOG(" vert " <<  posAcc[FractNode[lev][k][i]] );
//      }
//      UG_LOG(" \n ");
//    }
//  }
 
 
  //Mark Fracture Sides
  for (int lev = 0; lev < m_NumberLevels; lev++)
  {
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
      int si = m_fractSsGrp [i];
      t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != list_end; ++iElem)
      {
        element_type * elem = *iElem;
 
        //todo 3d faces
 
        //  Loop over the edges: We look for sides with only inner or only outer corners
        pMG->associated_elements_sorted (edge_list, elem);
        for (size_t j = 0; j < edge_list.size (); j++)
        {
          Edge * edge = edge_list [j];
          bool has_inner = false, has_outer = false;
          // bool has_corner = false; // currently unused, silencing warning (avogel)
          bool has_vert = false;
          //todo: check corner case in 3d
          size_t n_co = edge->num_vertices ();
          for (size_t co = 0; co < n_co; co++)
          {
            int mark = vert_mark (edge->vertex (co));
            if (mark == F_BND)
              has_outer = true;
            else if (mark == F_MID)
              has_inner = true;
            else if (mark == F_CORNER)
              /* has_corner = true */ ; // currently unused, silencing warning (avogel)
            else if (mark == F_VERT)
              has_vert = true;
            else
              UG_THROW ("FractDimadapt:Fracture marks maybe not set correctly!");
          }
 
          if (has_inner && has_outer)
            m_aaEdgeAuxMarks [edge] = F_VERT;
          else if (has_vert)
            m_aaEdgeAuxMarks [edge] = F_VERT;
          else if (has_inner)
          { // this is the inner side
            m_aaEdgeAuxMarks [edge] = F_MID;
          }
          else
          { // this is the outer side
            m_aaEdgeAuxMarks [edge] = F_BND;
          }
        }
      }
    }
  }
 
 
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::density_func (
  number c, 
  number &rho, 
  number &d_rho 
  )
{
  UG_COND_THROW(strcmp(m_density_type, "NONE") !=0, "FractDimadapt: density not set: Use 'add_density'!");
 
  if(strcmp(m_density_type,"linear")==0)
  {
    rho = m_rho_pw + ( m_rho_pb - m_rho_pw ) * c;
    d_rho = m_rho_pb - m_rho_pw;
  }
  else
    UG_THROW ("FractDimadapt: density type "<<m_density_type << "not yet considered: use 'linear'!");
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::viscosity_func (
  number c, 
  number &mu, 
  number &d_mu 
  )
{
  UG_COND_THROW(strcmp(m_viscosity_type, "NONE") !=0, "FractDimadapt: viscosity not set: Use 'add_viscosity'!");
 
  if(strcmp(m_viscosity_type,"const")==0)
  {
    mu = m_zero_viscosity;
    d_mu = 0.0;
  }
  else
    UG_THROW ("FractDimadapt: density type "<<m_viscosity_type << "not yet considered: use 'const'!");
}
 
 
template <typename TDomain,typename TAlgebra>
template <typename TGridFunction>
void FractDimadapt<TDomain,TAlgebra>::get_ave_val (
    const std::vector<number> pnt, 
    TGridFunction& u, 
    number time 
    )
{
  UG_LOG("get_ave_val\n");
 
  const position_accessor_type& posAcc = m_dom->position_accessor();
 
  //  create Multiindex
  std::vector<DoFIndex> multInd;
 
  number eps = 1E-10;
 
  number c_1,c_2,c_ave=0,c_jump=0;
  number p_1,p_2,p_ave=0,p_jump=0;
 
  FILE *file;
 
  MathVector<dim> eval_point;
  for (int d=0;d < dim;d++)
    eval_point[d] = pnt[d];
  int toplevel = m_NumberLevels-1;
  MathVector<dim> coords_1,coords_2,diff;
 
  if(m_FirstFullLevel < 0)
    UG_THROW ("FractDimadapt:get_ave_val: FirstFullLevel not set");
 
 
 
  //Find Fracture mid node corresponding to eval point
  int check = 0;
  for (int k = 0; k < NumberFractNodes[toplevel]; k++)
  {
    coords_1 = posAcc[FractNode[toplevel][k][0]];
    number dist = VecDistance(coords_1,eval_point);
    if(dist < eps)
    {
      //if fracture is full-dimensional
      //then integrate over fracture width.
      if(m_FractureIsFullDimensional == 1)
      {
        //the nodes over which is integrated are stored in FractNode
        Vertex * vert = FractNode[toplevel][k][0];
        u.dof_indices(vert, 0, multInd);
        c_ave = 2*DoFRef(u,multInd[0]);
        u.dof_indices(vert, 1, multInd);
        p_ave = 2*DoFRef(u,multInd[0]);
        for (int i = 3; i < NumberFractVertNodes[toplevel][k]+3; i++)
        {
          u.dof_indices(FractNode[toplevel][k][i], 0, multInd);
          c_ave += 2*DoFRef(u,multInd[0]);
          u.dof_indices(FractNode[toplevel][k][i], 1, multInd);
          p_ave += 2*DoFRef(u,multInd[0]);
        }
 
        coords_1 = posAcc[FractNode[toplevel][k][1]];
        u.dof_indices(FractNode[toplevel][k][1], 0, multInd);
        c_1=DoFRef(u,multInd[0]);
        u.dof_indices(FractNode[toplevel][k][1], 1, multInd);
        p_1=DoFRef(u,multInd[0]);
        coords_2 = posAcc[FractNode[toplevel][k][2]];
        u.dof_indices(FractNode[toplevel][k][2], 0, multInd);
        c_2=DoFRef(u,multInd[0]);
        u.dof_indices(FractNode[toplevel][k][2], 1, multInd);
        p_2=DoFRef(u,multInd[0]);
 
        // distinguish upper and lower boundary
        //jump: upper-lower
        VecSubtract(diff,coords_1,coords_2);
        if(diff[dim-1]>1e-10) //coords_1 is upper
        {
          c_jump = c_1 - c_2;
          p_jump = p_1 - p_2;
        }
        else if(diff[dim-1]<-1e-10) //coords_2 is upper
        {
          c_jump = c_2 - c_1;
          p_jump = p_2 - p_1;
        }
        else if(fabs(diff[dim-1])<1e-10)
        {
          if(fabs(diff[0])<1e-10 && fabs(diff[dim-2])<1e-10)
          {
            UG_THROW ("FractDimadapt:get_ave_val: full-dim fracture with Fracture width 0!");
          }
          UG_THROW ("FractDimadapt:get_ave_val: Horizontal fracture not yet considered for jump");
        }
 
        c_ave += c_1 + c_2;
        p_ave += p_1 + p_2;
        c_ave /= (NumberFractVertNodes[toplevel][k]+1)*2+2;
        p_ave /= (NumberFractVertNodes[toplevel][k]+1)*2+2;
      }
      else //fracture is low-dimensional
      {
        //average is stored in mid node
        Vertex * vert = FractNode[toplevel][k][0];
        u.dof_indices(vert, 0, multInd);
        c_ave = DoFRef(u,multInd[0]);
        u.dof_indices(vert, 1, multInd);
        p_ave = DoFRef(u,multInd[0]);
 
        u.dof_indices(FractNode[toplevel][k][1], 0, multInd);
        c_1=DoFRef(u,multInd[0]);
        u.dof_indices(FractNode[toplevel][k][1], 1, multInd);
        p_1=DoFRef(u,multInd[0]);
 
        u.dof_indices(FractNode[toplevel][k][2], 0, multInd);
        c_2=DoFRef(u,multInd[0]);
        u.dof_indices(FractNode[toplevel][k][2], 1, multInd);
        p_2=DoFRef(u,multInd[0]);
 
        //jump: upper-lower
        c_jump = c_1 - c_2;
        p_jump = p_1 - p_2;
        MathVector<dim> shift = vert_shift(FractNode[toplevel][k][1]);
        if(shift[dim-1]>1e-10) //c_1 is lower
        {
          c_jump = c_2 - c_1;
          p_jump = p_2 - p_1;
        }
        else if(shift[dim-1]<1e-10) //c_1 is upper
        {
          c_jump = c_1 - c_2;
          p_jump = p_1 - p_2;
        }
        else if(fabs(shift[dim-1])<1e-10)
        {
          if(fabs(shift[0])<1e-10 && fabs(shift[1])<1e-10)
          {
            UG_THROW ("FractDimadapt:get_ave_val: shift not set!");
          }
          UG_THROW ("FractDimadapt:get_ave_val: Horizontal fracture not yet considered for jump");
        }
        //add neglected pressure drop over fracture
        //cf. e.g. [Reiter, S. et al. Models and simulations of variable-density flow
        //in fractured porous media. IJSCE Vol. 9, Nos. 5/6, 2014]
        if (m_shift == true)
        {
          number buffer = 0.0;
          number density = 0.0;
          density_func(c_ave,density,buffer);
          number p_hydro = m_aperture*density*m_Gravity[dim-1];
          p_jump += p_hydro;
        }
      }
 
      check = 1;
      k=NumberFractNodes[toplevel];
    }
  }
  if(check==0)
    UG_THROW ("FractDimadapt:get_ave_val: No inner fracture node at this position found!");
 
  UG_LOG("\t t " << time<<" c_ave "<<c_ave<<" c_1 "<<c_1<<" c_2 "<<c_2<<" c_jump "<<c_jump <<"\n");
  UG_LOG("\t t " << time<<" p_ave "<<p_ave<<" p_1 "<<p_1<<" p_2 "<<p_2<<" p_jump "<<p_jump <<"\n");
 
  file = fopen("ave_c.dat","a");
  fprintf(file, "%g %g %g\n",time,c_ave,c_jump);
  fclose(file);
  file = fopen("ave_p.dat","a");
  fprintf(file, "%g %g %g\n",time,p_ave,p_jump);
  fclose(file);
}
 
template <typename TDomain,typename TAlgebra>
template <typename TGridFunction>
number FractDimadapt<TDomain,TAlgebra>::evaluate_criterion(
    TGridFunction& u,
    DegeneratedLayerManager<TDomain::dim>* fractManager,
    number omega_theta,
    number time
    )
{
  //DegeneratedLayerManager<TDomain::dim>* fractManager = NULL;
  UG_LOG("evaluate_criterion\n");
 
  if(m_permeability_f==0.0)
    UG_THROW ("FractDimadapt: Permeability not set. Use 'add_permeability'");
 
  if (! fractManager->is_closed ())
    UG_THROW ("FractDimadapt:evaluate_crit: Fracture manager not closed");
 
  typedef DegeneratedLayerManager<TDomain::dim> fract_manager_type;
  typedef typename fract_manager_type::side_type DFside_type;
  static const size_t maxLayerSideCorners = fract_manager_type::maxLayerSideCorners;
 
  typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
 
  const position_accessor_type& posAcc = m_dom->position_accessor();
 
  //finite volume geometry
  DimFV1Geometry<dim> fullgeo; //for full-dim
  TFractFVGeom *lowgeo;
  lowgeo = &GeomProvider<TFractFVGeom>::get (LFEID (LFEID::LAGRANGE, low_dim, 1), 1); //for low-dim
 
 
  MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
 
  //  create Multiindex
  std::vector<DoFIndex> multInd;
 
  number ViscosityIP;
  number DensityIP;
  number d_Density = 0.0;
  number d_Visc = 0.0;
 
  number theta_max = 0;
  number theta_inv_max = 0;
  number theta_min = 0;
  number q_max = 0;
  number q_rot_max = 0;
  number test_max = 0;
 
  number scale = m_permeability_f;
 
  FILE *file;
 
  //compute elementwise
  for (size_t j = 0; j < m_fractSsGrp.size (); j++)
  {
    int si = m_fractSsGrp [j];
    bool isFracture = fractManager->contains (si);
 
    for (t_elem_iter iElem = m_spSH->template begin<element_type> (si,m_NumberLevels-1);
        iElem != m_spSH->template end<element_type> (si,m_NumberLevels-1); ++iElem)
    {
      element_type * elem = *iElem;
 
      //  get vertices
      const size_t numVertices = elem->num_vertices();
 
      //extract corner coordinates
      for(size_t i = 0; i < numVertices; ++i)
        coCoord[i] = posAcc[elem->vertex(i)];
 
      if(isFracture) //low-dim
      {
        typedef StdLinConsistentGravity<low_dim> TConsGravity;
        TConsGravity ConsGravityMethod;
          MathVector<low_dim> vConsGravity [maxLayerSideCorners];
 
          size_t n_co, inner_side_idx, outer_side_idx;
        size_t inner_side_corners [maxLayerSideCorners];
        size_t outer_side_corners [maxLayerSideCorners];
        DFside_type * inner_side, * outer_side;
        size_t ass_co [2*maxLayerSideCorners];
 
        number DensityCo[maxLayerSideCorners];
        number ViscosityCo[maxLayerSideCorners];
        number cValue[maxLayerSideCorners];
        number pValue[maxLayerSideCorners];
 
        //  get the non-degenerated sides of the fracture element
        fractManager->get_layer_sides (elem, n_co,
          inner_side, inner_side_idx, inner_side_corners,
          outer_side, outer_side_idx, outer_side_corners, ass_co);
 
        //  compute the FV geometry of the inner side
        MathVector<dim> vSideCornerCoords [maxLayerSideCorners];
        try
        {
          for (size_t co = 0; co < n_co; co++)
            vSideCornerCoords [co] = posAcc [elem->vertex(inner_side_corners[co])];
          lowgeo->update (inner_side, vSideCornerCoords, m_spSH.get());
        }
        UG_CATCH_THROW("EvaluateCriterion: Cannot update the Finite Volume Geometry for a fracture element.");
          const size_t numSh = lowgeo->num_sh();
        const size_t numScvf = lowgeo->num_scvf();
 
        //get c,rho, mu and p in corners
        for (size_t sh=0; sh < numSh; sh++)
        {
          size_t co = inner_side_corners[sh];
          Vertex *vert = elem->vertex(co);
          u.dof_indices(vert, 0, multInd);
          cValue[sh] = DoFRef(u,multInd[0]);
          density_func(cValue[sh],DensityCo[sh],d_Density);
          viscosity_func(cValue[sh],ViscosityCo[sh],d_Visc);
          u.dof_indices(vert, 1, multInd);
          pValue[sh] = DoFRef(u,multInd[0]);
        }
 
        //  Prepare the consistent gravity
        try
        {
          ConsGravityMethod.template prepare<dim>
            (vConsGravity, numSh, lowgeo->corners(), DensityCo, m_Gravity);
        }
        UG_CATCH_THROW ("EvaluateCriterion: Cannot prepare Consistent Gravity.");
 
        for (size_t ip=0; ip < numScvf; ip++)
        {
          number absDarcy = 0;
          number vorticity, q_rot;
          number theta = 0;
          number theta_inv = 0;
          MathVector<dim> DarcyVel,Vel;
          MathVector<dim> grad_c_ip;
          MathVector<dim> grad_p_ip;
          number c_ip;
 
          const typename TFractFVGeom::SCVF& scvf = lowgeo->scvf(ip);
 
          VecSet (grad_c_ip, 0.0);
          VecSet (grad_p_ip, 0.0);
          c_ip = 0;
          for(size_t co = 0 ; co < scvf.num_sh(); ++co){
            VecScaleAppend(grad_c_ip,cValue[co],scvf.global_grad(co));
            VecScaleAppend(grad_p_ip,pValue[co],scvf.global_grad(co));
            c_ip += cValue[co] * scvf.shape(co);
          };
          density_func(c_ip,DensityIP,d_Density);
          viscosity_func(c_ip,ViscosityIP,d_Visc);
 
          // DARCY VEL
          //  Compute rho * g (as the consistent gravity force)
          ConsGravityMethod.template compute<dim>
            (Vel, scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(), vConsGravity);
 
          // The pressure-gradient part:
          VecScaleAppend (Vel, -1, grad_p_ip);
          // The viscosity and permeability factor:
          VecScale (DarcyVel, Vel, m_permeability_f / ViscosityIP);
          absDarcy = VecLength(DarcyVel);
          absDarcy = fabs(absDarcy);
 
          // Q_ROT
          //compute vorticity
          vorticity = m_permeability_f / ViscosityIP * d_Density * m_Gravity[1] * grad_c_ip[0];
 
          //compute q_rot
          q_rot = m_aperture/2 * m_permeability_f/m_permeability_m * vorticity
              * c_ip / omega_theta;
 
          //compute theta
          if(q_rot!=0)
            theta = absDarcy / fabs(q_rot);
          if(absDarcy != 0)
            theta_inv = fabs(q_rot) / absDarcy;
 
          //max
          if(theta > theta_max)
            theta_max = theta;
          if(theta_inv > theta_inv_max)
            theta_inv_max = theta_inv;
          if(absDarcy > q_max)
            q_max = absDarcy;
          if(q_rot > q_rot_max)
            q_rot_max = q_rot;
 
          //min
          if(theta < theta_min || theta_min==0)
            theta_min = theta;
        }
      }
      else //full-dim
      {
        std::vector<number> DensityCo(numVertices);
        std::vector<number> ViscosityCo(numVertices);
        std::vector<number> cValue(numVertices);
        std::vector<number> pValue(numVertices);
 
        //  Consistent gravity and its derivative at corners
        typedef StdLinConsistentGravityX<element_type::dim> TConsGravity;
        TConsGravity ConsGravityMethod;
        MathVector<dim> vConsGravity [domain_traits<dim>::MaxNumVerticesOfElem];
 
        //  evaluate finite volume geometry
        fullgeo.update(elem, &(coCoord[0]), u.domain()->subset_handler().get());
        const size_t numSh = fullgeo.num_sh();
        const size_t numScvf = fullgeo.num_scvf();
 
        for (size_t co=0; co < numVertices; co++)
        {
          Vertex *vert = elem->vertex(co);
          u.dof_indices(vert, 0, multInd);
          cValue[co] = DoFRef(u,multInd[0]);
          density_func(cValue[co],DensityCo[co],d_Density);
          viscosity_func(cValue[co],ViscosityCo[co],d_Visc);
          u.dof_indices(vert, 1, multInd);
          pValue[co] = DoFRef(u,multInd[0]);
        }
 
        //  Prepare the consistent gravity
        try
        {
          ConsGravityMethod.template prepare<dim>
          (vConsGravity, numSh, coCoord, &DensityCo[0], m_Gravity);
        }
        UG_CATCH_THROW ("EvaluateCriterion: Cannot prepare Consistent Gravity.");
 
        for (size_t ip=0;ip < numScvf;ip++)
        {
          number absDarcy = 0;
          number vorticity, q_rot;
          number theta = 0;
          number theta_inv = 0;
          MathVector<dim> DarcyVel,Vel;
          MathVector<dim> grad_c_ip;
          MathVector<dim> grad_p_ip;
          number c_ip;
          number vort_test,test;
 
          const typename DimFV1Geometry<dim>::SCVF& scvf = fullgeo.scvf(ip);
 
          VecSet (grad_c_ip, 0.0);
          VecSet (grad_p_ip, 0.0);
          c_ip = 0;
          for(size_t co = 0 ; co < scvf.num_sh(); ++co){
            VecScaleAppend(grad_c_ip,cValue[co],scvf.global_grad(co));
            VecScaleAppend(grad_p_ip,pValue[co],scvf.global_grad(co));
            c_ip += cValue[co] * scvf.shape(co);
          };
          density_func(c_ip,DensityIP,d_Density);
          viscosity_func(c_ip,ViscosityIP,d_Visc);
 
          // DARCY VEL
          //  Compute rho * g (as the consistent gravity force)
          ConsGravityMethod.template compute<dim>
            (Vel, scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(), vConsGravity);
 
          // The pressure-gradient part:
          VecScaleAppend (Vel, -1, grad_p_ip);
 
          // The viscosity and permeability factor:
          VecScale (DarcyVel, Vel, m_permeability_f / ViscosityIP); //MatVecMult if perm tensor
          absDarcy = VecLength(DarcyVel);
          absDarcy = fabs(absDarcy);
 
          // Q_ROT
          //compute vorticity
          vorticity = m_permeability_f / ViscosityIP * d_Density * m_Gravity[1] * grad_c_ip[0];
 
          vort_test = m_permeability_f / ViscosityIP * 25.0/1025.0 * 1000 / (1 - 25.0/1025.0 *c_ip) / (1 - 25.0/1025.0 *c_ip)* m_Gravity[1] * grad_c_ip[0];
 
          //compute q_rot
          q_rot = m_aperture/2 * m_permeability_f/m_permeability_m * vorticity * c_ip / omega_theta;
 
          test = m_aperture * m_permeability_f/m_permeability_m * vort_test * c_ip / 50.0 /10.0 /0.25;
 
          //compute theta
          if(q_rot!=0)
            theta = absDarcy / fabs(q_rot);
          if(absDarcy != 0)
            theta_inv = fabs(q_rot) / absDarcy;
 
          test = fabs(test);
 
          //max
          if(theta > theta_max)
            theta_max = theta;
          if(theta_inv > theta_inv_max)
            theta_inv_max = theta_inv;
          if(absDarcy > q_max)
            q_max = absDarcy;
          if(q_rot > q_rot_max)
            q_rot_max = q_rot;
          if(test > test_max)
            test_max = test;
 
          //min
          if(theta < theta_min || theta_min==0)
            theta_min = theta;
        }
      }
    }
  }
 
  //number theta =  q_max / q_rot_max;
  number test = q_max / test_max;
  //UG_LOG("\n\t Max: q_max: "<<  q_max<<" q_rot_max: " << q_rot_max<<" q_rot_UG3_max: " << test_max<<"\n");
  //UG_LOG("\n\t Max: theta_max: "<< theta_max <<" q_max / q_rot_max: " << theta <<" q_max / q_rot_UG3_max: " << test<<"\n");
 
  UG_LOG("\t Min: theta_min: "<<  theta_min<<" theta_inv_max: " << theta_inv_max<<"\n");
 
  file = fopen("crit.dat","a");
  fprintf(file, "%g: theta_min: %g theta_inv_max: %g theta_UG3 %g\n",time,theta_min/scale,theta_inv_max*scale,test);
  fclose(file);
 
  return theta_inv_max*scale;
}
 
 
template <typename TDomain,typename TAlgebra>
template <typename TGridFunction>
void FractDimadapt<TDomain,TAlgebra>::check_transfer(
    TGridFunction& u,
    DegeneratedLayerManager<TDomain::dim>* fractManager,
    number omega_theta,
    number delta,
    number time
    )
{
  UG_LOG("check_transfer\n");
 
  //  do we have fractures?
  if (fractManager == NULL) return;
 
  for (size_t i = 0; i < m_fractSsGrp.size (); i++)
  {
    int si = m_fractSsGrp [i];
 
    bool isDegenerated = fractManager->contains (si);
 
    UG_LOG("\n\t "<< isDegenerated << " "<< m_FractureIsFullDimensional<<"\n");
 
    if(m_FractureIsFullDimensional == 1 && isDegenerated == true)
      fractManager->remove(m_spSH->get_subset_name(si));
 
    isDegenerated = fractManager->contains (si);
    UG_LOG("\n\t "<< isDegenerated << " "<< m_FractureIsFullDimensional<<"\n");
 
    number theta;
    theta = evaluate_criterion(u,fractManager,omega_theta,time);
 
    if(theta < (1 - delta) )
    {
      UG_LOG("\n\t Should be low-dimensional\n");
      if(m_FractureIsFullDimensional == 1)
      {
        UG_LOG("\n\t switch\n");
        //first calculate average and store it in mid node
        transfer_full_to_low(u);
        //then transfer the grid (other inner fracture nodes are deleted)
        grid_transfer_full_to_low ();
        convert_end_triangles_to_quadris(u,fractManager);
        fractManager->add(m_spSH->get_subset_name(si));
        fractManager->close();
 
        isDegenerated = fractManager->contains (si);
        UG_LOG("\n\t "<< isDegenerated << " "<< m_FractureIsFullDimensional<<"\n");
      }
    }
    if(theta > (1 + delta) )
    {
      UG_LOG("\n\t Should be full-dimensional\n");
      if(m_FractureIsFullDimensional == 0)
      {
        UG_LOG("\n\t switch\n");
        //first transfer grid
        convert_end_quadris_to_triangles(fractManager);
        grid_transfer_low_to_full ();
        //then calculate value of unknown for the new inner fracture nodes
        transfer_low_to_full(u);
        fractManager->remove(m_spSH->get_subset_name(si));
        fractManager->close();
 
        isDegenerated = fractManager->contains (si);
        UG_LOG("\n\t "<< isDegenerated << " "<< m_FractureIsFullDimensional<<"\n");
      }
    }
  }
}
 
 
template <typename TDomain,typename TAlgebra>
template <typename TGridFunction>
void FractDimadapt<TDomain,TAlgebra>::transfer_full_to_low (
    TGridFunction& u
    )
{
  //  create Multiindex
  std::vector<DoFIndex> multInd;
 
  UG_LOG("transfer full to low \n");
 
  int lev = m_NumberLevels-1; //transfer is only needed on toplevel
  number c_ave,c_1,c_2;
  number p_ave,p_1,p_2;
  number p_1_korr,p_2_korr;
  number density = 0.0;
 
  for (int k = 0; k < NumberFractNodes[lev]; k++)
  {
    //integrate over fracture width.
    //the nodes over which is integrated are stored in FractNode
    Vertex * vert = FractNode[lev][k][0];
    u.dof_indices(vert, 0, multInd);
    c_ave = 2*DoFRef(u,multInd[0]);
    u.dof_indices(vert, 1, multInd);
    p_ave = 2*DoFRef(u,multInd[0]);
 
    for (int i = 3; i < NumberFractVertNodes[lev][k]+3; i++)
    {
      u.dof_indices(FractNode[lev][k][i], 0, multInd);
      c_ave += 2*DoFRef(u,multInd[0]);
      u.dof_indices(FractNode[lev][k][i], 1, multInd);
      p_ave += 2*DoFRef(u,multInd[0]);
    }
 
    u.dof_indices(FractNode[lev][k][1], 0, multInd);
    c_1=DoFRef(u,multInd[0]);
    u.dof_indices(FractNode[lev][k][1], 1, multInd);
    p_1=DoFRef(u,multInd[0]);
    u.dof_indices(FractNode[lev][k][2], 0, multInd);
    c_2=DoFRef(u,multInd[0]);
    u.dof_indices(FractNode[lev][k][2], 1, multInd);
    p_2=DoFRef(u,multInd[0]);
 
    c_ave += c_1 + c_2;
    p_ave += p_1 + p_2;
    c_ave /= (NumberFractVertNodes[lev][k]+1)*2+2;
    p_ave /= (NumberFractVertNodes[lev][k]+1)*2+2;
 
    //correct pressure (vertices are shifted, therefore hydrostatic pressure needs to be adapted)
    if(m_shift == true)
    {
      number buffer = 0.0;
      density_func(c_1,density,buffer);
      MathVector<dim> shift = vert_shift(FractNode[lev][k][1]);
      p_1_korr = p_1 + shift[dim-1]*density*m_Gravity[dim-1];
      density_func(c_2,density,buffer);
      shift = vert_shift(FractNode[lev][k][2]);
      p_2_korr = p_2 + shift[dim-1]*density*m_Gravity[dim-1];
 
      //write corrected pressure in bnd_verts
      u.dof_indices(FractNode[lev][k][1], 1, multInd);
      DoFRef(u,multInd[0]) = p_1_korr;
      u.dof_indices(FractNode[lev][k][2], 1, multInd);
      DoFRef(u,multInd[0]) = p_2_korr;
    }
 
    //write average in mid_vert
    u.dof_indices(vert, 0, multInd);
    DoFRef(u,multInd[0]) = c_ave;
    u.dof_indices(vert, 1, multInd);
    DoFRef(u,multInd[0]) = p_ave;
  }
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::calculate_interpolation (
  number bnd[2][2],
  number mid[2],
  MathVector<dim> bndpoint[2],
  number param[2][3]  //return value
  )
{
   number a[2][2],b[2][2],tem,temp,temp1,temp2,temp4,temp5;
   int i,j;
   int n=2,m=2;
   int p,q;
   number length;
   MathVector<dim> vec,x;
 
   //do not multiply mistakes!!
   for(i=0;i<2;i++)
   {
      if(fabs(bnd[i][0])<1E-10 && bnd[i][0]<0)
        bnd[i][0]=0;
   }
   if(fabs(mid[0])<1E-10 && mid[0]<0)
        mid[0]=0;
 
   //|x_2 - x_1|
   VecSubtract(x,bndpoint[1],bndpoint[0]);
   length = VecLength(x);
 
    a[0][0]=length*length;
    a[0][1]=length;
    a[1][0]=length*length*length/3;
    a[1][1]=length*length/2;
 
   //Calculate Inverse of A=a_ij
   for(i=0;i<n;i++)
   {
      for(j=0;j<m;j++)
    {
      if(i==j)
        b[i][j]=1;
      else
        b[i][j]=0;
    }
   }
   for(i=0;i<n;i++)
   {
      temp=a[i][i];
      if(temp<0)
        temp=temp*(-1);
      p=i;
      for(j=i+1;j<n;j++)
    {
        if(a[j][i]<0)
          tem=a[j][i]*(-1);
        else
          tem=a[j][i];
        if(temp<0)
          temp=temp*(-1);
        if(tem>temp)
          {
            p=j;
            temp=a[j][i];
          }
    }
      //row exchange in both the matrix
      for(j=0;j<n;j++)
      {
        temp1=a[i][j];
      a[i][j]=a[p][j];
      a[p][j]=temp1;
      temp2=b[i][j];
      b[i][j]=b[p][j];
      b[p][j]=temp2;
      }
      //dividing the row by a[i][i]
      temp4=a[i][i];
      for(j=0;j<n;j++)
    {
      a[i][j]=(number)a[i][j]/temp4;
      b[i][j]=(number)b[i][j]/temp4;
    }
        //making other elements 0 in order to make the matrix a[][] an indentity matrix and obtaining a inverse b[][] matrix
      for(q=0;q<n;q++)
    {
      if(q==i)
        continue;
      temp5=a[q][i];
      for(j=0;j<n;j++)
        {
          a[q][j]=a[q][j]-(temp5*a[i][j]);
          b[q][j]=b[q][j]-(temp5*b[i][j]);
        }
    }
    //b is inverse
  }
 
    for(i=0;i<2;i++)
  {
      vec[0] = (bnd[1][i]-bnd[0][i]); //u_2 - u_1
      vec[1] = (mid[i]-bnd[0][i])*length; // |x_2-x_1|*(u_mid - u_1)
      param[i][0] = b[0][0]*vec[0] + b[0][1]*vec[1];
      param[i][1] = b[1][0]*vec[0] + b[1][1]*vec[1];
      param[i][2] = bnd[0][i];  // u_1
  }
}
 
template <typename TDomain,typename TAlgebra>
template <typename TGridFunction>
void FractDimadapt<TDomain,TAlgebra>::transfer_low_to_full (
    TGridFunction& u 
    )
{
  //  create Multiindex
  std::vector<DoFIndex> multInd;
 
  const position_accessor_type& posAcc = m_dom->position_accessor();
 
  UG_LOG("transfer low to full \n");
 
  int lev = m_NumberLevels-1; //transfer only needed on toplevel
  for (int k = 0; k < NumberFractNodes[lev]; k++)
  {
    number mid[2],bnd[2][2];
    number inner[2];
    MathVector<dim> bndpoint[2];
    MathVector<dim> calcpoint;
    number param[2][3];
    number x;
    number density = 0.0;
    number buffer = 0.0; //d_dens not needed here
 
    //get values in mid node (=average)
    Vertex * vert = FractNode[lev][k][0];
    u.dof_indices(vert, 0, multInd);
    mid[0] = DoFRef(u,multInd[0]);
    u.dof_indices(vert, 1, multInd);
    mid[1] = DoFRef(u,multInd[0]);
 
    //get values in bnd nodes
    // and correct pressure (vertices are shifted, therefore hydrostatic pressure needs to be adapted)
    // bnd 1
    vert = FractNode[lev][k][1];
    u.dof_indices(vert, 0, multInd);
    bnd[0][0]=DoFRef(u,multInd[0]);
    u.dof_indices(vert, 1, multInd);
    bnd[0][1]=DoFRef(u,multInd[0]);
    bndpoint[0]=posAcc[vert];
    //correct p
    if(m_shift == true)
    {
      density_func(bnd[0][0],density,buffer);
      MathVector<dim> shift = vert_shift(vert);
      bnd[0][1] = bnd[0][1] - shift[dim-1]*density*m_Gravity[dim-1];
      DoFRef(u,multInd[0]) = bnd[0][1];
    }
 
    // bnd 2
    vert = FractNode[lev][k][2];
    u.dof_indices(vert, 0, multInd);
    bnd[1][0]=DoFRef(u,multInd[0]);
    u.dof_indices(vert, 1, multInd);
    bnd[1][1]=DoFRef(u,multInd[0]);
    bndpoint[1]=posAcc[vert];
    //correct p
    if(m_shift == true)
    {
      density_func(bnd[1][0],density,buffer);
      MathVector<dim> shift = vert_shift(vert);
      bnd[1][1] = bnd[1][1] - shift[dim-1]*density*m_Gravity[dim-1];
      DoFRef(u,multInd[0]) = bnd[1][1];
    }
 
 
    //assume quadratic behaviour of unknowns in fracture
    //calculate parameters a,b,c
    calculate_interpolation(bnd,mid,bndpoint,param);
    size_t use_linear = 0;
    for (int i = 0; i < NumberFractVertNodes[lev][k]+3; i++)
    {
      //boundary values are not changed
      // though, considering them here, should not change anything.
      if(i==1 || i==2)
        continue;
      vert = FractNode[lev][k][i];
      calcpoint = posAcc[vert];
      //|x - x_1|
      VecSubtract(calcpoint,bndpoint[0],calcpoint);
      x = VecLength(calcpoint);
      //calculate and set value
      for(size_t j=0; j<2; j++)
      {
        inner[j] = param[j][0]*x*x + param[j][1]*x + param[j][2]; //f(x)= ax^2+bx+c
//          UG_LOG("\t inner " <<  posAcc[vert] <<" "<< bnd[0][j] <<" "<<bnd[1][j] <<" "
//              <<mid[j]<<" "<<inner[j]<< " \n" );
        u.dof_indices(vert, j, multInd);
        DoFRef(u,multInd[0]) = inner[j];
      }
 
      //if computed unknown are out of bounds (e.g. negative concentration)
      //use a linear interpolation instead
      if(inner[0]<-1e-10 || inner[0]>1 || inner[1]<-1e-10)
      {
        UG_LOG ("FractDimadapt:transfer_low_to_full: c or p neg. Use linear fct instead");
        use_linear = 1;
        i=NumberFractVertNodes[lev][k]+3;
      }
    }
    if(use_linear == 1)
    {
      for (int i = 0; i < NumberFractVertNodes[lev][k]+3; i++)
      {
        //boundary values are not changed
        if(i==1 || i==2)
          continue;
        vert = FractNode[lev][k][i];
        //|x - x_1|
        VecSubtract(calcpoint,bndpoint[0],posAcc[vert]);
        x = VecLength(calcpoint);
        //|x2 - x_1|
        VecSubtract(calcpoint,bndpoint[0],bndpoint[1]);
        x = x / VecLength(calcpoint);
        //calculate and set value
        for(size_t j=0; j<2; j++)
        {
          inner[j] = (1-x) * bnd[0][j] + x * bnd[1][j];
//            UG_LOG("\t inner " <<  posAcc[vert] <<" "<< bnd[0][j] <<" "<<bnd[1][j] <<" "
//                <<mid[j]<<" "<<inner[j]<< " \n" );
          u.dof_indices(vert, j, multInd);
          DoFRef(u,multInd[0]) = inner[j];
        }
      }
    }
 
  }
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::grid_transfer_full_to_low ()
{
  typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
 
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
 
  typename Grid::traits<side_type>::secure_container side_list;
  typename Grid::traits<Edge>::secure_container edge_list;
 
  if(m_FirstFullLevel < 0)
    UG_THROW ("FractDimadapt:grid_transfer_full_to_low: FirstFullLevel not set");
 
  typedef typename domain_type::position_attachment_type pos_att_type;
  pos_att_type aPos = m_dom->position_attachment();
  Grid::VertexAttachmentAccessor<pos_att_type> aaPos;
  aaPos.access (*pMG, aPos);
 
  UG_LOG("grid transfer full to low \n");
 
  //begin of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_BEGINS));
 
  //delete triangles (because they have no vert fract node)
  //todo: check for 3d. edges/sides??
  for (int lev = m_FirstFullLevel; lev < m_NumberLevels; lev++)
  {
    for(int i=0; i<4; i++)
    {
      element_type * elem = CornerTriangleList[lev][i];
 
      pMG->associated_elements_sorted (edge_list, elem);
 
      std::vector<Edge*> edge_to_erase(edge_list.size ());
      int counter = 0;
      for (size_t j = 0; j < edge_list.size (); j++)
      {
        Edge * edge = edge_list [j];
        if(edge_mark(edge)==F_MID || edge_mark(edge)==F_BND)
          continue;
        edge_to_erase[counter]= edge;
        counter++;
      }
      for(int j=0; j<counter; j++)
      {
        pMG->erase(edge_to_erase[j]);
      }
    }
  }
 
  //delete vert fract nodes (then all edges and faces depending on them are deleted too)
  for (int lev = 0; lev < m_NumberLevels; lev++)
  {
    for (int k = 0; k < NumberFractNodes[lev]; k++)
    {
      for (int i = 3; i < NumberFractVertNodes[lev][k]+3; i++)
      {
        pMG->erase(FractNode[lev][k][i]);
      }
      NumberFractVertNodes[lev][k] = 0;
    }
  }
 
  //insert new elements
  int baselev = m_FirstFullLevel-1;
  int toplev = m_NumberLevels-1;
  for (int lev = baselev; lev < toplev; lev++)
  {
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
      int si = m_fractSsGrp [i];
      t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != list_end; ++iElem)
      {
        element_type * elem = *iElem;
 
        //Find the corner vertices for the new elements
        int numvert = elem->num_vertices ();
        Vertex *v[4];//MaxNumCornerElements
        Vertex *v2[4];
        pMG->associated_elements_sorted (edge_list, elem);
 
        //for each element two new elems need to be created
        int counter = 0;
        int counter2 = 0;
        int done = 0;
        int triangle_done = 0;
        for (size_t j = 0; j < edge_list.size (); j++)
        {
          Edge * edge = edge_list [j];
          //edge has no vertex children: both endpoints are corners of new elem
          if(pMG->num_child_vertices(edge)==0)
          {
            //Find edge co numbering in correspondence to elem co numbering
            int num0,num1;
            for(int k=0; k<numvert; k++)
            {
              if(elem->vertex(k)==edge->vertex(0))
              {
                if(elem->vertex((k+1)%numvert)==edge->vertex(1))//std
                {
                  num0 = 0;
                  num1 = 1;
                }
                else if(elem->vertex((numvert+k-1)%numvert)==edge->vertex(1))
                {
                  num0 = 1;
                  num1 = 0;
                }
                else
                  UG_THROW ("FractDimadapt:grid_transfer_full_to_low: Elem Edge numbering\n");
              }
            }
            if(done==0) //corners for the first elem
            {
              v[counter] = pMG->get_child_vertex(edge->vertex(num0));
              counter++;
              v[counter] = pMG->get_child_vertex(edge->vertex(num1));
              counter++;
              done = 1;
            }
            else //corners for the second elem
            {
              v2[counter2] = pMG->get_child_vertex(edge->vertex(num0));
              counter2++;
              v2[counter2] = pMG->get_child_vertex(edge->vertex(num1));
              counter2++;
              done = 1;
            }
          }
          else //edge has vertex child: this is corner of new elem
          {
            v[counter] = pMG->get_child_vertex(edge);
            counter++;
            if(numvert==3 && triangle_done ==0) //for triangle
            {
              for(int k=0; k<numvert; k++)
              {
                Vertex *vert =elem->vertex(k);
                if(vert_mark(vert)==F_CORNER)
                {
                  if(elem->vertex(k)==edge->vertex(0))
                  {
                    if(elem->vertex((k+1)%numvert)==edge->vertex(1))//std
                    {
                      v2[counter2] = pMG->get_child_vertex(vert);
                      counter2++;
                      v2[counter2] = pMG->get_child_vertex(edge);
                      counter2++;
                    }
                    else if(elem->vertex((numvert+k-1)%numvert)==edge->vertex(1))
                    {
                      v2[counter2] = pMG->get_child_vertex(edge);
                      counter2++;
                      v2[counter2] = pMG->get_child_vertex(vert);
                      counter2++;
                    }
                    else
                      UG_THROW ("FractDimadapt:grid_transfer_full_to_low: Triangle\n");
                  }
                  else if(elem->vertex(k)==edge->vertex(1))
                  {
                    if(elem->vertex((k+1)%numvert)==edge->vertex(0))
                    {
                      v2[counter2] = pMG->get_child_vertex(edge);
                      counter2++;
                      v2[counter2] = pMG->get_child_vertex(vert);
                      counter2++;
                    }
                    else if(elem->vertex((numvert+k-1)%numvert)==edge->vertex(0))//std
                    {
                      v2[counter2] = pMG->get_child_vertex(vert);
                      counter2++;
                      v2[counter2] = pMG->get_child_vertex(edge);
                      counter2++;
                    }
                    else
                      UG_THROW ("FractDimadapt:grid_transfer_full_to_low: Triangle\n");
                  }
                  else
                    UG_THROW ("FractDimadapt:grid_transfer_full_to_low:Triangle\n");
                }
              }
              triangle_done = 1;
            }
            else //quadrilateral and second triangle side
            {
              v2[counter2] = pMG->get_child_vertex(edge);
              counter2++;
            }
          }
        }
 
        //for quadirlaterals two new quadrilaterals are created
        if(numvert==4)
        {
          if(counter!=4)
            UG_THROW ("FractDimadapt:grid_transfer_full_to_low: wrong number corners"<<counter<<"\n");
          pMG->create<Quadrilateral>(QuadrilateralDescriptor(v[0],v[1],v[2],v[3]),elem);
          if(counter2!=4)
            UG_THROW ("FractDimadapt:grid_transfer_full_to_low: wrong number corners"<<counter2<<"\n");
          pMG->create<Quadrilateral>(QuadrilateralDescriptor(v2[0],v2[1],v2[2],v2[3]),elem);
        }
        else if(numvert==3)//for triangles one quadrilateral and one triangle are created
        {
          if(counter!=4)
            UG_THROW ("FractDimadapt:grid_transfer_full_to_low: wrong number corners"<<counter<<"\n");
          pMG->create<Quadrilateral>(QuadrilateralDescriptor(v[0],v[1],v[2],v[3]),elem);
          if(counter2!=3)
            UG_THROW ("FractDimadapt:grid_transfer_full_to_low: wrong number corners"<<counter2<<"\n");
          pMG->create<Triangle>(TriangleDescriptor(v2[0],v2[1],v2[2]),elem);
        }
      }
    }
 
  }
 
  //shift fracture boundary vertices so that fracture has width 0
  if(m_shift == true)
  {
    for (int lev = 0; lev < m_NumberLevels; lev++)
    {
      for (int k = 0; k < NumberFractNodes[lev]; k++)
      {
        MathVector<dim> newcoord;
        if(FractNode[lev][k][1]!=NULL && FractNode[lev][k][2]!=NULL)
        {
          MathVector<dim> shift = vert_shift(FractNode[lev][k][1]);
 
          //Todo SHIFT Debug
          number debug_shift = 1e-10;
          if(shift[1]>0)
            shift[1] += debug_shift;
          else if (shift[1]<0)
            shift[1] -= debug_shift;
 
          VecAdd(newcoord,aaPos[FractNode[lev][k][1]],shift);
          aaPos[FractNode[lev][k][1]] = newcoord;
 
 
          shift = vert_shift(FractNode[lev][k][2]);
 
          //Todo SHIFT Debug
          if(shift[1]>0)
            shift[1] += debug_shift;
          else if (shift[1]<0)
            shift[1] -= debug_shift;
 
          VecAdd(newcoord,aaPos[FractNode[lev][k][2]],shift);
          aaPos[FractNode[lev][k][2]] = newcoord;
        }
      }
    }
  }
 
 
  m_FractureIsFullDimensional = 0;
 
  //end of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_ENDS));
 
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::grid_transfer_low_to_full ()
{
  typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
 
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
 
  typename Grid::traits<side_type>::secure_container side_list;
  typename Grid::traits<Edge>::secure_container edge_list;
 
  typedef typename domain_type::position_attachment_type pos_att_type;
  pos_att_type aPos = m_dom->position_attachment();
  Grid::VertexAttachmentAccessor<pos_att_type> aaPos;
  aaPos.access (*pMG, aPos);
 
  if(m_FirstFullLevel < 0)
    UG_THROW ("FractDimadapt:grid_transfer_low_to_full: FirstFullLevel not set");
 
  UG_LOG("grid transfer low to full \n");
 
  //begin of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_BEGINS));
 
  //shift fracture boundary vertices so that fracture has no longer width 0
  if(m_shift == true)
  {
    for (int lev = 0; lev < m_NumberLevels; lev++)
    {
      for (int k = 0; k < NumberFractNodes[lev]; k++)
      {
        if(FractNode[lev][k][1]!=NULL && FractNode[lev][k][2]!=NULL)
        {
          MathVector<dim> newcoord;
          MathVector<dim> shift = vert_shift(FractNode[lev][k][1]);
 
          //Todo SHIFT Debug
          number debug_shift = 1e-10;
          if(shift[1]>0)
            shift[1] -= debug_shift;
          else if (shift[1]<0)
            shift[1] += debug_shift;
 
          VecSubtract(newcoord,aaPos[FractNode[lev][k][1]],shift);
          aaPos[FractNode[lev][k][1]] = newcoord;
          shift = vert_shift(FractNode[lev][k][2]);
 
          //Todo SHIFT Debug
          if(shift[1]>0)
            shift[1] -= debug_shift;
          else if (shift[1]<0)
            shift[1] += debug_shift;
 
          VecSubtract(newcoord,aaPos[FractNode[lev][k][2]],shift);
          aaPos[FractNode[lev][k][2]] = newcoord;
        }
      }
    }
  }
 
  //todo: check for 3d. edges/sides??
 
  //delete vert edges. Then all fracture elems are deleted as well
  for (int lev = m_FirstFullLevel; lev < m_NumberLevels; lev++)
  {
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
      int si = m_fractSsGrp [i];
      Edge * edge_to_erase[200];
      int counter = 0;
      t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != list_end; ++iElem)
      {
        element_type * elem = *iElem;
        pMG->associated_elements_sorted (edge_list, elem);
        for (size_t j = 0; j < edge_list.size (); j++)
        {
          Edge * edge = edge_list [j];
          if(edge_mark(edge)==F_MID || edge_mark(edge)==F_BND || edge_mark(edge)==F_NONE)
            continue;
          edge_to_erase[counter]= edge;
          m_aaEdgeAuxMarks [edge] = F_NONE; // mark as none to avoid double counting
          counter++;
          if(counter==200)
            UG_THROW ("FractDimadapt:grid_transfer_low_to_full: number of edges out of bounds");
        }
      }
      for(int j=0; j<counter; j++)
      {
        pMG->erase(edge_to_erase[j]);
      }
    }
  }
 
  //now insert new vertices
  int baselev = m_FirstFullLevel-1;
  int toplev = m_NumberLevels-1;
  for (int lev = baselev; lev < toplev; lev++)
  {
    for (size_t i = 0; i < m_fractSsGrp.size (); i++)
    {
      int si = m_fractSsGrp [i];
      t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != list_end; ++iElem)
      {
        element_type * elem = *iElem;
        Vertex* v[4][4];
        if(elem->num_vertices ()==4)
        {
          //insert vrt and calculate new position
          Vertex *newvrt = *pMG->create<RegularVertex>(elem);
          aaPos[newvrt] = CalculateCenter(elem, aaPos);
          //insert vertex recursively on finer grids.
          Vertex *newvrt_old = newvrt;
          for(int recursive_lev = lev; recursive_lev < toplev-1; recursive_lev++)
          {
            newvrt = *pMG->create<RegularVertex>(newvrt_old);
            aaPos[newvrt] = aaPos[newvrt_old];
            newvrt_old = newvrt;
          }
        }
 
        pMG->associated_elements_sorted (edge_list, elem);
        for (size_t j = 0; j < edge_list.size (); j++)
        {
          Edge * edge = edge_list [j];
          if(pMG->num_child_vertices(edge)==0)
          {
            //insert vrt and calculate new position
            Vertex *newvrt = *pMG->create<RegularVertex>(edge);
            aaPos[newvrt] = CalculateCenter(edge, aaPos);
            //insert corresponding edges
            pMG->create<RegularEdge>(EdgeDescriptor(pMG->get_child_vertex(edge->vertex(0)),newvrt),edge);
            pMG->create<RegularEdge>(EdgeDescriptor(newvrt,pMG->get_child_vertex(edge->vertex(1))),edge);
            //insert vertex recursively on finer grids.
            Vertex *newvrt_old = newvrt;
            for(int recursive_lev = lev; recursive_lev < toplev-1; recursive_lev++)
            {
              newvrt = *pMG->create<RegularVertex>(newvrt_old);
              aaPos[newvrt] = aaPos[newvrt_old];
              newvrt_old = newvrt;
            }
          }
        }
 
        //list corners of new elems
        if(elem->num_vertices ()==4)
        {
          for(size_t j = 0; j < 4; ++j)
          {
            v[j][0] = pMG->get_child_vertex(elem->vertex(j));
          //  v[j][1] = pMG->get_child_vertex(pMG->get_edge(elem, j));
            v[j][1] = pMG->get_child_vertex(edge_list [j]);
            v[j][2] = pMG->get_child_vertex(elem);
            v[j][3] = pMG->get_child_vertex(edge_list [(j+3)%edge_list.size ()]);
            //UG_LOG("vrtlist quad " <<posAcc[v[j][0]]<<" "<<posAcc[v[j][1]]<<" "<<posAcc[v[j][2]]<<" "<<posAcc[v[j][3]]<<"\n ");
            pMG->create<Quadrilateral>(QuadrilateralDescriptor(v[j][0],v[j][1],v[j][2],v[j][3]),elem);
          }
        }
        else if(elem->num_vertices ()==3)
        {
          for (size_t j = 0; j < edge_list.size (); j++)
          {
            Edge * edge = edge_list [j];
            v[j][0] = pMG->get_child_vertex(elem->vertex(j));
            v[j][1] = pMG->get_child_vertex(edge);
            v[j][2] = pMG->get_child_vertex(edge_list [(j+2)%edge_list.size ()]);
            //UG_LOG("vrtlist tri 1 " <<posAcc[v[j][0]]<<" "<<posAcc[v[j][1]]<<" "<<posAcc[v[j][2]]);
            pMG->create<Triangle>(TriangleDescriptor(v[j][0],v[j][1],v[j][2]),elem);
          }
          v[3][0] = pMG->get_child_vertex(edge_list [0]);
          v[3][1] = pMG->get_child_vertex(edge_list [1]);
          v[3][2] = pMG->get_child_vertex(edge_list [2]);
          //UG_LOG("vrtlist tri 2 " <<posAcc[v[3][0]]<<" "<<posAcc[v[3][1]]<<" "<<posAcc[v[3][2]]);
          pMG->create<Triangle>(TriangleDescriptor(v[3][0],v[3][1],v[3][2]),elem);
        }
        else
          UG_THROW ("FractDimadapt:grid_transfer_low_to_full: only triangles and quadrilaterals so far");
      }
    }
  }
 
  m_FractureIsFullDimensional = 1;
 
  //end of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_ENDS));
}
 
template <typename TDomain,typename TAlgebra>
template <typename TGridFunction>
void FractDimadapt<TDomain,TAlgebra>::convert_end_triangles_to_quadris (
    TGridFunction& u,
    DegeneratedLayerManager<TDomain::dim>* fractManager
)
{
  UG_LOG("convert_tri_to_quad \n");
 
  //  create Multiindex
  std::vector<DoFIndex> multInd;
 
  typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
 
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
 
  typename Grid::traits<Edge>::secure_container edge_list;
 
 
  typedef typename domain_type::position_attachment_type pos_att_type;
  pos_att_type aPos = m_dom->position_attachment();
  Grid::VertexAttachmentAccessor<pos_att_type> aaPos;
  aaPos.access (*pMG, aPos);
 
 
  //begin of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_BEGINS));
 
  //each new inserted vertex is saved to be used on the next level as parent
  const size_t max_number_deg_quadris = 2;
  Vertex * new_vrt_father[max_number_deg_quadris];
  Edge * new_degedge_father[max_number_deg_quadris];
  Edge * new_fulledge_father[max_number_deg_quadris];
  Face* new_face1_father[max_number_deg_quadris];
  size_t number_new_vrts = 0;
 
  size_t max_verts = max_number_deg_quadris * m_NumberLevels;
  std::vector<std::vector<Vertex*> > vert_pairs(max_verts);
  for(size_t i = 0; i < vert_pairs.size(); ++i) vert_pairs[i].resize(2);
//  Vertex * vert_pairs[max_verts][2];
  size_t vert_count = 0;
 
  //degenerated fractures in UG4 can only end with quadrilaterals
  //in the full-dim case the fractures end with triangles
  for (int lev = 0; lev < m_NumberLevels; lev++)
  {
    size_t triangle_counter = 0;
    Edge * corneredge[4];
    Face* cornerface[4];
    Vertex * cornervrt[4];
    Vertex * vrtlist[4][4];
    number sign[4];
    size_t partner[4];
    size_t edgenum[4]; //for ordering
 
 
 
    //find end elements of fracture (triangle + cornernode)
    for (size_t h = 0; h < m_fractSsGrp.size (); h++)
    {
      int si = m_fractSsGrp [h];
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != m_spSH->template end<element_type> (si, lev); ++iElem)
      {
        element_type * elem = *iElem;
 
        if(elem->num_vertices()==3)
        {
          pMG->associated_elements_sorted (edge_list, elem);
          for (size_t j = 0; j < edge_list.size (); j++)
          {
            Edge * edge = edge_list [j];
            if(edge_mark(edge)==F_MID)
            {
              if(vert_mark(edge->vertex(0))==F_CORNER)
              {
                cornervrt[triangle_counter] = edge->vertex(0);
                number s = aaPos[edge->vertex(1)][0] - aaPos[edge->vertex(0)][0];
                sign[triangle_counter] = s/(fabs(s));
              }
              else if(vert_mark(edge->vertex(1))==F_CORNER)
              {
                cornervrt[triangle_counter] = edge->vertex(1);
                number s = aaPos[edge->vertex(0)][0] - aaPos[edge->vertex(1)][0];
                sign[triangle_counter] = s/(fabs(s));
              }
              else
              {
                UG_LOG("triangle mid edge with no corner node \n ");
                continue; //oder Fehler bei deg case
              }
              corneredge[triangle_counter] = edge;
              edgenum[triangle_counter] = j;
 
              Face *nbfaces[2];
              if(GetAssociatedFaces(nbfaces,*pMG,edge,2)!=2)
                UG_THROW("FractDimadapt:tri_to_quad: edge has more than two faces\n");
 
              int equal[2];
              for(size_t j=0; j< 2; j++)
              {
                equal[j]=1;
                for(size_t k = 0; k < nbfaces[j]->num_vertices(); k++)
                  if(nbfaces[j]->vertex(k) != elem->vertex(k))
                  {
                    k = nbfaces[j]->num_vertices();
                    equal[j] = 0;
                  }
              }
              if(equal[0]==0 && equal[1]==0)
                UG_THROW("FractDimadapt:tri_to_quad: no elem found\n");
              for(size_t j=0; j< 2; j++)
                if(equal[j]==1)
                  cornerface[triangle_counter] = nbfaces[j];
 
              //list corners for new quadrilateral
              size_t count = 0;
              for(size_t k=0; k< elem->num_vertices(); k++)
              {
                vrtlist[triangle_counter][count+k] = elem->vertex(k);
                if(k==edgenum[triangle_counter])
                {
                  count=1;
                  vrtlist[triangle_counter][k+count] = NULL; //to be set later (newvrt)
                }
              }
 
              triangle_counter++;
              if(triangle_counter>(2*max_number_deg_quadris))
                UG_THROW ("FractDimadapt:tri_to_quad: too many corner triangles");
            }
          }
        }
      }
 
      //find corresponding pairs. Two triangles at each end.
      for (size_t i = 0; i < triangle_counter; i++)
      {
        for (size_t j = 0; j < triangle_counter; j++)
        {
          if(i==j)
            continue;
          if(cornervrt[i] == cornervrt[j])
            partner[i]=j;
        }
      }
 
      //insert new vertices, edges and elems.
      //correct child structure and erase old parts
      int new_count = 0;
      for (size_t i = 0; i < triangle_counter; i++)
      {
        if(partner[i]<i)
          continue;
 
        //special treatment for baselevel
        //as there are no parents here, subdomain info needs to be set
        if(lev == 0)
        {
          pMG->enable_hierarchical_insertion(false); //use only for base level! to set subset info for new vertices...
 
          //insert vrt
          Vertex *newvrt = *pMG->create<RegularVertex>(corneredge[i]);
          if(m_shift == true)
          {
            //Todo SHIFT Debug
            number debug_shift = 1e-10;
 
            aaPos[newvrt][0] = aaPos[cornervrt[i]][0] + sign[i]*debug_shift;
            aaPos[newvrt][1] = aaPos[cornervrt[i]][1];
          }
          else
          {
            MathVector<dim> tangential;
            tangential[0]= m_fractNormal[1];
            tangential[1]= -m_fractNormal[0];
            number step_x = fabs(aaPos[corneredge[i]->vertex(0)][0] - aaPos[corneredge[i]->vertex(1)][0]);
            number step_y= fabs(aaPos[corneredge[i]->vertex(0)][1] - aaPos[corneredge[i]->vertex(1)][1]);
            for(int l=0; l<m_NumberLevels; l++)
            {
              step_x /= 2;
              step_y /= 2;
            }
            step_x /= 5;
            step_y /= 5;
            aaPos[newvrt][0] = aaPos[cornervrt[i]][0] + fabs(tangential[0]) * sign[i] * step_x;
            aaPos[newvrt][1] = aaPos[cornervrt[i]][1] + fabs(tangential[1]) * sign[i] * step_y;
          }
 
          new_vrt_father[new_count] = newvrt;
          //UG_LOG("new vrt "<<aaPos[newvrt] <<" \n ");
 
          //save new vert and corner vert to assign c and p values later (after vector has new length)
          vert_pairs[vert_count][0] = newvrt;
          vert_pairs[vert_count][1] = cornervrt[i];
          vert_count++;
 
          //insert corresponding edges
          Edge *newedge1 = *pMG->create<RegularEdge>(EdgeDescriptor(corneredge[i]->vertex(0),newvrt),corneredge[i]);
          Edge *newedge2 = *pMG->create<RegularEdge>(EdgeDescriptor(newvrt,corneredge[i]->vertex(1)),corneredge[i]);
 
          pMG->enable_hierarchical_insertion(true);
 
          Edge *newparent;
          if(corneredge[i]->vertex(0)==cornervrt[i])
          {
            newparent = newedge2;
            new_fulledge_father[new_count] = newedge2;
            new_degedge_father[new_count] = newedge1;
          }
          else if(corneredge[i]->vertex(1)==cornervrt[i])
          {
            newparent = newedge1;
            new_fulledge_father[new_count] = newedge1;
            new_degedge_father[new_count] = newedge2;
          }
          else
            UG_THROW("FractDimadapt:tri_to_quad: some error in collected data\n");
 
          //set new parents of child_edges and child_vertices
          for(size_t j = 0; j < pMG->num_child_edges(corneredge[i]); j++)
            pMG->associate_parent(pMG->get_child_edge(corneredge[i],j),newparent);
          for(size_t j = 0; j < pMG->num_child_vertices(corneredge[i]); j++)
            pMG->associate_parent(pMG->get_child_vertex(corneredge[i]),newparent);
 
          //add newrt to list corners for new quadrilaterals.
          vrtlist[i][edgenum[i]+1] = newvrt;
          vrtlist[partner[i]][edgenum[partner[i]]+1] = newvrt;
 
//          UG_LOG("vrtlist " <<aaPos[vrtlist[i][0]]<<" "<<aaPos[vrtlist[i][1]]<<" "
//              <<aaPos[vrtlist[i][2]]<<" "<<aaPos[vrtlist[i][3]]<<"\n ");
//          UG_LOG("vrtlist " <<aaPos[vrtlist[partner[i]][0]]<<" "<<aaPos[vrtlist[partner[i]][1]]<<" "
//              <<aaPos[vrtlist[partner[i]][2]]<<" "<<aaPos[vrtlist[partner[i]][3]]<<"\n ");
 
          Face* face = NULL;
          //insert new elements
          //set subdoms
          // set parents of children of old elements
          Face* Fchildren_list[4]; //max_num_children macro rausfinden
          Edge* Echildren_list[2];
          size_t num_children = 0;
 
          face = *pMG->create<Quadrilateral>(QuadrilateralDescriptor(vrtlist[i][0],vrtlist[i][1],
            vrtlist[i][2],vrtlist[i][3]));
          m_spSH->assign_subset(face, si);
 
          num_children = pMG->num_child_faces(cornerface[i]);
          for(size_t k = 0; k < num_children; k++)
            Fchildren_list[k] = pMG->get_child_face(cornerface[i],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Fchildren_list[k],face);
 
          num_children = pMG->num_child_edges(cornerface[i]);
          for(size_t k = 0; k < num_children; k++)
            Echildren_list[k] = pMG->get_child_edge(cornerface[i],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Echildren_list[k],face);
 
          new_face1_father[new_count] = face;
 
          face = *pMG->create<Quadrilateral>(QuadrilateralDescriptor(vrtlist[partner[i]][0],vrtlist[partner[i]][1],
              vrtlist[partner[i]][2],vrtlist[partner[i]][3]));
          m_spSH->assign_subset(face, si);
 
          num_children = pMG->num_child_faces(cornerface[partner[i]]);
          for(size_t k = 0; k < num_children; k++)
            Fchildren_list[k] = pMG->get_child_face(cornerface[partner[i]],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Fchildren_list[k],face);
 
          num_children = pMG->num_child_edges(cornerface[partner[i]]);
          for(size_t k = 0; k < num_children; k++)
            Echildren_list[k] = pMG->get_child_edge(cornerface[partner[i]],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Echildren_list[k],face);
 
          new_count++;
          number_new_vrts++;
 
          //erase old edge and therefore old triangles
          pMG->erase(corneredge[i]);
        }
        else //lev >0
        {
          //UG_LOG("lev "<<lev<<"\n");
          //get parent elem
          Face* parent1 = dynamic_cast<Face*>(pMG->get_parent(cornerface[i]));
          Face* parent2 = dynamic_cast<Face*>(pMG->get_parent(cornerface[partner[i]]));
 
          //find out which saved new_vrt is the parent of the vrt to be inserted
          //by comparing the parent face with the saved face
          size_t this_id = max_number_deg_quadris+10;
          for(size_t k = 0; k < number_new_vrts; k++)
          {
            if(parent1 == new_face1_father[k] || parent2 == new_face1_father[k])
              this_id = k;
          }
          if(this_id == max_number_deg_quadris+10)
            UG_THROW("FractDimadapt:tri_to_quad: parent not found in list (lev "<<lev<<")\n");
 
//          //insert vrt as child of the one on the coarser lev
          Vertex *newvrt = *pMG->create<RegularVertex>(new_vrt_father[this_id]);
          aaPos[newvrt][0] = aaPos[new_vrt_father[this_id]][0];
          aaPos[newvrt][1] = aaPos[new_vrt_father[this_id]][1];
 
          //save new vert and corner vert to assign c and p values later (after vector has new length)
          vert_pairs[vert_count][0] = newvrt;
          vert_pairs[vert_count][1] = cornervrt[i];
          vert_count++;
 
          new_vrt_father[new_count] = newvrt;
          //UG_LOG("lev "<<lev<<" new vrt "<<aaPos[newvrt] <<" \n ");
 
          //insert corresponding edges
          Edge *newdegedge = NULL;
          Edge *newfulledge = NULL;
 
          if(corneredge[i]->vertex(0)==cornervrt[i])
          {
            newdegedge = *pMG->create<RegularEdge>(EdgeDescriptor(corneredge[i]->vertex(0),newvrt),new_degedge_father[this_id]);
            newfulledge = *pMG->create<RegularEdge>(EdgeDescriptor(newvrt,corneredge[i]->vertex(1)),new_fulledge_father[this_id]);
          }
          else if(corneredge[i]->vertex(1)==cornervrt[i])
          {
            newfulledge = *pMG->create<RegularEdge>(EdgeDescriptor(corneredge[i]->vertex(0),newvrt),new_fulledge_father[this_id]);
            newdegedge = *pMG->create<RegularEdge>(EdgeDescriptor(newvrt,corneredge[i]->vertex(1)),new_degedge_father[this_id]);
          }
          else
            UG_THROW("FractDimadapt:tri_to_quad: some error in collected data\n");
 
          new_fulledge_father[new_count] = newfulledge;
          new_degedge_father[new_count] = newdegedge;
 
          //set new parents of children
          for(size_t j = 0; j < pMG->num_child_edges(corneredge[i]); j++)
            pMG->associate_parent(pMG->get_child_edge(corneredge[i],j),newfulledge);
          for(size_t j = 0; j < pMG->num_child_vertices(corneredge[i]); j++)
            pMG->associate_parent(pMG->get_child_vertex(corneredge[i]),newfulledge);
 
          //add newrt to list corners for new quadrilaterals.
          vrtlist[i][edgenum[i]+1] = newvrt;
          vrtlist[partner[i]][edgenum[partner[i]]+1] = newvrt;
 
//          UG_LOG("vrtlist " <<aaPos[vrtlist[i][0]]<<" "<<aaPos[vrtlist[i][1]]<<" "
//              <<aaPos[vrtlist[i][2]]<<" "<<aaPos[vrtlist[i][3]]<<"\n ");
//          UG_LOG("vrtlist " <<aaPos[vrtlist[partner[i]][0]]<<" "<<aaPos[vrtlist[partner[i]][1]]<<" "
//              <<aaPos[vrtlist[partner[i]][2]]<<" "<<aaPos[vrtlist[partner[i]][3]]<<"\n ");
 
          //insert new elements
          // and set parents of children of old elements
          Face* face = NULL;
          Face* Fchildren_list[4]; //max_num_children macro rausfinden
          Edge* Echildren_list[2];
          size_t num_children = 0;
 
          face = *pMG->create<Quadrilateral>(QuadrilateralDescriptor(vrtlist[i][0],vrtlist[i][1],
            vrtlist[i][2],vrtlist[i][3]),parent1);
 
          num_children = pMG->num_child_faces(cornerface[i]);
          for(size_t k = 0; k < num_children; k++)
            Fchildren_list[k] = pMG->get_child_face(cornerface[i],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Fchildren_list[k],face);
 
          num_children = pMG->num_child_edges(cornerface[i]);
          for(size_t k = 0; k < num_children; k++)
            Echildren_list[k] = pMG->get_child_edge(cornerface[i],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Echildren_list[k],face);
          new_face1_father[new_count] = face;
 
 
          face = *pMG->create<Quadrilateral>(QuadrilateralDescriptor(vrtlist[partner[i]][0],vrtlist[partner[i]][1],
              vrtlist[partner[i]][2],vrtlist[partner[i]][3]),parent2);
 
          num_children = pMG->num_child_faces(cornerface[partner[i]]);
          for(size_t k = 0; k < num_children; k++)
            Fchildren_list[k] = pMG->get_child_face(cornerface[partner[i]],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Fchildren_list[k],face);
 
          num_children = pMG->num_child_edges(cornerface[partner[i]]);
          for(size_t k = 0; k < num_children; k++)
            Echildren_list[k] = pMG->get_child_edge(cornerface[partner[i]],k);
          for(size_t k = 0; k < num_children; k++)
            pMG->associate_parent(Echildren_list[k],face);
 
          new_count++;
 
          //erase old edge and therefore old triangles
          pMG->erase(corneredge[i]);
        }
//        SubsetHandler &sh = *m_spSH.get();
//        SaveGridHierarchyTransformed(*pMG, m_spSH.get(), "hallo.ugx", 1);
      }
    }
  }
 
  //debug
//  for (int lev = 0; lev < m_NumberLevels; lev++)
//  {
//    for (int k = 0; k < NumberFractNodes[lev]; k++)
//    {
//      MathVector<dim> newcoord;
//      MathVector<dim> shift = vert_shift(FractNode[lev][k][1]);
//      VecAdd(newcoord,aaPos[FractNode[lev][k][1]],shift);
//      aaPos[FractNode[lev][k][1]] = newcoord;
//      shift = vert_shift(FractNode[lev][k][2]);
//      VecAdd(newcoord,aaPos[FractNode[lev][k][2]],shift);
//      aaPos[FractNode[lev][k][2]] = newcoord;
//    }
//  }
 
 
 
  //update fract manager
  fractManager->close();
 
  //end of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_ENDS));
 
 
  //assign c and p values at new node (copy those of cornernode)
  for(size_t i=vert_count-2; i < vert_count; i++)
  {
    number c,p;
    u.dof_indices(vert_pairs[i][1], 0, multInd);
    c = DoFRef(u,multInd[0]);
    u.dof_indices(vert_pairs[i][0], 0, multInd);
    DoFRef(u,multInd[0]) = c;
 
    u.dof_indices(vert_pairs[i][1], 1, multInd);
    p = DoFRef(u,multInd[0]);
    u.dof_indices(vert_pairs[i][0], 1, multInd);
    DoFRef(u,multInd[0]) = p;
  }
 
}
 
template <typename TDomain,typename TAlgebra>
void FractDimadapt<TDomain,TAlgebra>::convert_end_quadris_to_triangles (
    DegeneratedLayerManager<TDomain::dim>* fractManager
)
{
  UG_LOG("convert_quad_to_tri \n");
 
  //  create Multiindex
  std::vector<DoFIndex> multInd;
 
  typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
 
  MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
 
  typename Grid::traits<Edge>::secure_container edge_list;
 
 
  typedef typename domain_type::position_attachment_type pos_att_type;
  pos_att_type aPos = m_dom->position_attachment();
  Grid::VertexAttachmentAccessor<pos_att_type> aaPos;
  aaPos.access (*pMG, aPos);
 
 
  //begin of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_BEGINS));
 
  size_t num_vert = 2 * m_NumberLevels;
  std::vector<Vertex*> vertex_to_erase(num_vert);
  size_t vert_count = 0;
 
  for (int lev = 0; lev < m_NumberLevels; lev++)
  {
    size_t triangle_counter = 0;
 
    Edge * corneredge[4];
    Edge * midedge[4];
    Face* cornerface[4];
 
    Vertex * cornervrt[4];
    Vertex * triquadvrt[4];
    Vertex * vrtlist[4][3];
    size_t partner[4];
    size_t edge_num[4]; //corneredge->vertex(edge_num) is cornervert
 
    for (size_t h = 0; h < m_fractSsGrp.size (); h++)
    {
      int si = m_fractSsGrp [h];
 
      //find the end quadrilaterals.
      //store all relevant information in arrays
      for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
          iElem != m_spSH->template end<element_type> (si, lev); ++iElem)
      {
        element_type * elem = *iElem;
        for (size_t j = 0; j < elem->num_vertices (); j++)
        {
          Vertex * vert = elem->vertex(j);
          //find the corner-quadri
          if(vert_mark(vert) == F_CORNER)
          {
            if(elem->num_vertices ()!=4)
              UG_THROW("FractDimadapt:quad_to_tri: corner elem is no quadrilateral\n");
 
            //save vertices and edges
            cornervrt[triangle_counter] = vert;
 
            pMG->associated_elements_sorted (edge_list, elem);
            for (size_t i = 0; i < edge_list.size (); i++)
            {
              Edge * edge = edge_list [i];
              if(vert_mark(edge->vertex(0)) == F_MID)
              {
                if(vert_mark(edge->vertex(1)) == F_MID)
                {
                  midedge[triangle_counter] = edge;
                }
                else if(vert_mark(edge->vertex(1)) == F_CORNER)
                {
                  corneredge[triangle_counter] = edge;
                  triquadvrt[triangle_counter] = edge->vertex(0);
                  edge_num[triangle_counter] = 1;
                }
              }
              else if(vert_mark(edge->vertex(0)) == F_CORNER)
              {
                if(vert_mark(edge->vertex(1)) == F_MID)
                {
                  corneredge[triangle_counter] = edge;
                  triquadvrt[triangle_counter] = edge->vertex(1);
                  edge_num[triangle_counter] = 0;
                }
              }
            }
 
            //list of vertices for new element
            size_t list_num = 0;
            for (size_t i = 0; i < elem->num_vertices (); i++)
            {
              Vertex * covert = elem->vertex(i);
              if(covert != triquadvrt[triangle_counter])
              {
                vrtlist[triangle_counter][list_num] = covert;
                list_num++;
              }
            }
 
            //save the face
            Face *nbfaces[2];
            if(GetAssociatedFaces(nbfaces,*pMG,corneredge[triangle_counter],2)!=2)
              UG_THROW("FractDimadapt:quad_to_tri: edge has more than two faces\n");
            int equal[2];
            for(size_t i=0; i< 2; i++)
            {
              equal[i]=1;
              for(size_t k = 0; k < nbfaces[i]->num_vertices(); k++)
                if(nbfaces[i]->vertex(k) != elem->vertex(k))
                {
                  k = nbfaces[i]->num_vertices();
                  equal[i] = 0;
                }
            }
            if(equal[0]==0 && equal[1]==0)
              UG_THROW("FractDimadapt:quad_to_tri: no elem found\n");
            for(size_t k=0; k < 2; k++)
              if(equal[k]==1)
                cornerface[triangle_counter] = nbfaces[k];
 
            triangle_counter++;
            j=elem->num_vertices ();
          }
        }
      }
 
      //find partner
      for (size_t i = 0; i < triangle_counter; i++)
      {
        for (size_t j = 0; j < triangle_counter; j++)
        {
          if(i==j)
            continue;
          if(cornervrt[i] == cornervrt[j])
            partner[i]=j;
        }
      }
 
      //create new elems
      //reset parent child connections
      for (size_t i = 0; i < triangle_counter; i++)
      {
        if(partner[i]<i)
          continue;
 
        //store vert in array (erased later after parent-child-connections are all clear)
        vertex_to_erase[vert_count] = triquadvrt[i];
        vert_count++;
 
        //insert new edges
        Edge *newedge = NULL;
        if(lev==0)
        {
          pMG->enable_hierarchical_insertion(false); //use only for base level! to set subset info for new vertices...
          if(edge_num[i]==0)
            newedge = *pMG->create<RegularEdge>(EdgeDescriptor(cornervrt[i],midedge[i]->vertex(1)),midedge[i]);
          else
            newedge = *pMG->create<RegularEdge>(EdgeDescriptor(midedge[i]->vertex(0),cornervrt[i]),midedge[i]);
          pMG->enable_hierarchical_insertion(true);
        }
        else
        {
          Edge * parentedge = dynamic_cast<Edge*>(pMG->get_parent(midedge[i]));
          if(edge_num[i]==0)
            newedge = *pMG->create<RegularEdge>(EdgeDescriptor(cornervrt[i],midedge[i]->vertex(1)),parentedge);
          else
            newedge = *pMG->create<RegularEdge>(EdgeDescriptor(midedge[i]->vertex(0),cornervrt[i]),parentedge);
        }
        //set new parents of child_edges and child_vertices
        for(size_t j = 0; j < pMG->num_child_edges(midedge[i]); j++)
          pMG->associate_parent(pMG->get_child_edge(midedge[i],j),newedge);
        for(size_t j = 0; j < pMG->num_child_vertices(midedge[i]); j++)
          pMG->associate_parent(pMG->get_child_vertex(midedge[i]),newedge);
 
//        UG_LOG("vrtlist " <<aaPos[vrtlist[i][0]]<<" "<<aaPos[vrtlist[i][1]]<<" "
//              <<aaPos[vrtlist[i][2]]<<"\n ");
//        UG_LOG("vrtlist " <<aaPos[vrtlist[partner[i]][0]]<<" "<<aaPos[vrtlist[partner[i]][1]]<<" "
//              <<aaPos[vrtlist[partner[i]][2]]<<"\n ");
 
        //insert new elements
        // and set parents of children of old elements
        Face* face = NULL;
        Face* Fchildren_list[4]; //max_num_children macro rausfinden
        Edge* Echildren_list[2];
        size_t num_children = 0;
 
        //first elem
        if(lev==0)
        {
          face = *pMG->create<Triangle>(TriangleDescriptor(vrtlist[i][0],vrtlist[i][1],
              vrtlist[i][2]));
          m_spSH->assign_subset(face, si);
        }
        else
        {
          Face * parentface = dynamic_cast<Face*>(pMG->get_parent(cornerface[i]));
          face = *pMG->create<Triangle>(TriangleDescriptor(vrtlist[i][0],vrtlist[i][1],
            vrtlist[i][2]),parentface);
        }
        num_children = pMG->num_child_faces(cornerface[i]);
        for(size_t k = 0; k < num_children; k++)
          Fchildren_list[k] = pMG->get_child_face(cornerface[i],k);
        for(size_t k = 0; k < num_children; k++)
          pMG->associate_parent(Fchildren_list[k],face);
 
        num_children = pMG->num_child_edges(cornerface[i]);
        for(size_t k = 0; k < num_children; k++)
          Echildren_list[k] = pMG->get_child_edge(cornerface[i],k);
        for(size_t k = 0; k < num_children; k++)
          pMG->associate_parent(Echildren_list[k],face);
 
        //second elem
        if(lev==0)
        {
          face = *pMG->create<Triangle>(TriangleDescriptor(vrtlist[partner[i]][0],vrtlist[partner[i]][1],
              vrtlist[partner[i]][2]));
          m_spSH->assign_subset(face, si);
        }
        else
        {
          Face * parentface = dynamic_cast<Face*>(pMG->get_parent(cornerface[partner[i]]));
          face = *pMG->create<Triangle>(TriangleDescriptor(vrtlist[partner[i]][0],vrtlist[partner[i]][1],
            vrtlist[partner[i]][2]),parentface);
        }
        num_children = pMG->num_child_faces(cornerface[partner[i]]);
        for(size_t k = 0; k < num_children; k++)
          Fchildren_list[k] = pMG->get_child_face(cornerface[partner[i]],k);
        for(size_t k = 0; k < num_children; k++)
          pMG->associate_parent(Fchildren_list[k],face);
 
        num_children = pMG->num_child_edges(cornerface[partner[i]]);
        for(size_t k = 0; k < num_children; k++)
          Echildren_list[k] = pMG->get_child_edge(cornerface[partner[i]],k);
        for(size_t k = 0; k < num_children; k++)
          pMG->associate_parent(Echildren_list[k],face);
      }
    }
  }
 
  //erase vrt and therefore corresponding edges/elems
  for (size_t i = 0; i < vert_count; i++)
  {
    pMG->erase(vertex_to_erase[i]);
  }
 
 
  //update fract manager
  fractManager->close();
 
  //end of grid adaption noticed by messenger
  pMG->message_hub()->post_message(GridMessage_Adaption(GMAT_HNODE_ADAPTION_ENDS));
}
 
 
 
} // namespace d3f
} // end namespace ug
 
/* End of File */