ls_volume_impl.h

Go to the documentation of this file.

/*
 * Copyright (c) 2015:  G-CSC, Goethe University Frankfurt
 * Authors: Christian Wehner, Dmitry Logashenko
 * 
 * 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.
 */
 
#include "lib_disc/reference_element/reference_element_traits.h"
#ifdef UG_PARALLEL
#include "pcl/pcl_process_communicator.h"
#endif
 
namespace ug{
namespace LevelSet{
 
/*---- Class 'LSVolume': ----*/
 
template <typename TGridFunc>
void LSVolume<TGridFunc>::on_subsets
(
  const char * ss_names 
)
{
  m_ssGrp.set_subset_handler (m_spLSF->domain()->subset_handler ());
  m_ssGrp.add (TokenizeString (ss_names));
}
 
template <typename TGridFunc>
void LSVolume<TGridFunc>::compute ()
{
//  The full-dim. grid element types for this dimension:
  typedef typename domain_traits<dim>::DimElemList ElemList;
  
  const int n_ss = m_spLSF->num_subsets ();
 
//  sum up all the volumes
  m_volume_plus = m_volume_minus = 0;
  m_ss_vol_plus.resize (n_ss); m_ss_vol_minus.resize (n_ss);
  for (int si = 0; si < n_ss; si++) m_ss_vol_plus[si] = m_ss_vol_minus[si] = 0;
  boost::mpl::for_each<ElemList> (AddVolumes (this));
  
#ifdef UG_PARALLEL
  allreduce_volumes ();
#endif
}
 
template <typename TGridFunc>
void LSVolume<TGridFunc>::compute
(
  SmartPtr<UserData<number, dim> > spDensity, 
  number time 
)
{
//  The full-dim. grid element types for this dimension:
  typedef typename domain_traits<dim>::DimElemList ElemList;
  
  const int n_ss = m_spLSF->num_subsets ();
 
//  sum up all the volumes
  m_volume_plus = m_volume_minus = 0;
  m_ss_vol_plus.resize (n_ss); m_ss_vol_minus.resize (n_ss);
  for (int si = 0; si < n_ss; si++) m_ss_vol_plus[si] = m_ss_vol_minus[si] = 0;
  boost::mpl::for_each<ElemList> (AddWeightedVolumes (this, spDensity.get (), time));
  
#ifdef UG_PARALLEL
  allreduce_volumes ();
#endif
}
 
template <typename TGridFunc>
template <typename TElem>
void LSVolume<TGridFunc>::add_volumes_of_all ()
{
  typedef typename grid_func_type::template traits<TElem>::const_iterator ElemIter;
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_t;
 
  static const size_t num_corners = ref_elem_t::numCorners;
  
  const grid_func_type & lsf = * m_spLSF;
  const position_accessor_type & aaPos = lsf.domain()->position_accessor ();
  std::vector<DoFIndex> ind (1);
  MathVector<dim> corners [num_corners];
  number lsf_values [num_corners];
  
  for (int si = 0; si < lsf.num_subsets (); si++)
  {
    if (m_ssGrp.subset_handler().valid () && ! m_ssGrp.contains (si))
      continue; // skip this subset: it is not mentioned in the specified list
    
    number ss_vol_plus = 0, ss_vol_minus = 0;
    ElemIter iterEnd = lsf.template end<TElem> (si);
    for (ElemIter iter = lsf.template begin<TElem> (si); iter != iterEnd; ++iter)
    {
      TElem * elem = *iter;
      
    //  get the corner coordinates ans the values of the LSF
      for (size_t i = 0; i < num_corners; i++)
      {
        Vertex * vrt = elem->vertex (i);
        corners [i] = aaPos [vrt];
        if (lsf.inner_dof_indices (vrt, 0, ind) != 1)
          UG_THROW ("LSVolume: Not a scalar grid function for the LSF!");
        lsf_values [i] = DoFRef (lsf, ind [0]);
      }
      
    //  compute the volumes
      number vol_plus, vol_minus;
      LSElementSize<ref_elem_t, dim>::compute (corners, lsf_values, vol_plus, vol_minus);
      
    //  add them to the total volumes
      m_volume_plus += vol_plus; m_volume_minus += vol_minus;
      
    //  add them to the subset
      ss_vol_plus += vol_plus; ss_vol_minus += vol_minus;
      
    //  check the positivity (if requested)
      if (m_check_positivity && (vol_plus < 0 || vol_minus < 0))
      {
        UG_LOG ("Warning in LSVolume: V+ = " << vol_plus << " < 0 or V- = "
            << vol_minus << " < 0\n");
        for (size_t i = 0; i < num_corners; i++)
        {
          UG_LOG (" - corner[" << i << "] = " << corners[i]
            << ", grid_data_index = " << elem->vertex(i)->grid_data_index ()
            << ", lsf = " << lsf_values[i] << '\n');
        }
      }
      
      /*-- For debugging only: --*
      number test_vol_plus, test_vol_minus;
      for (size_t i = 0; i < num_corners; i++)
        lsf_values [i] = - lsf_values [i];
      LSElementSize<ref_elem_t, dim>::compute (corners, lsf_values, test_vol_minus, test_vol_plus);
      if (test_vol_minus < 0 || test_vol_plus < 0 || vol_minus < 0 || vol_plus < 0
        || std::abs (test_vol_minus - vol_minus) / (vol_plus + vol_minus) >= 1e-12
          || std::abs (test_vol_plus - vol_plus) / (vol_plus + vol_minus) >= 1e-12)
      {
        UG_LOG ("---- LSVolume: Inconsistent values of the volumes for elem id " << elem->grid_data_index () << ":\n");
        UG_LOG ("-- V_(-) = " << vol_minus << " or " << test_vol_minus << ", diff = " << vol_minus - test_vol_minus << "\n");
        UG_LOG ("-- V_(+) = " << vol_plus << " or " << test_vol_plus << ", diff = " << vol_plus - test_vol_plus << "\n");
        for (size_t i = 0; i < num_corners; i++)
        {
          UG_LOG ("-- co " << i << ": " << corners[i] << ", lsf = " << - lsf_values [i] << "\n");
        }
        UG_THROW ("LSVolume: Inaccurate computation of the volumes.");
      }
       *--*/
    }
    m_ss_vol_plus[si] += ss_vol_plus; m_ss_vol_minus[si] += ss_vol_minus;
  }
}
 
template <typename TGridFunc>
template <typename TElem>
void LSVolume<TGridFunc>::add_volumes_of_all
(
  UserData<number, dim> * pDensity, 
  number time 
)
{
  typedef typename grid_func_type::template traits<TElem>::const_iterator ElemIter;
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_t;
 
  static const size_t num_corners = ref_elem_t::numCorners;
  
  const grid_func_type & lsf = * m_spLSF;
  const position_accessor_type & aaPos = lsf.domain()->position_accessor ();
  std::vector<DoFIndex> ind (1);
  MathVector<dim> corners [num_corners], midpoint;
  number lsf_values [num_corners];
  
  for (int si = 0; si < lsf.num_subsets (); si++)
  {
    if (m_ssGrp.subset_handler().valid () && ! m_ssGrp.contains (si))
      continue; // skip this subset: it is not mentioned in the specified list
    
    number ss_vol_plus = 0, ss_vol_minus = 0;
    ElemIter iterEnd = lsf.template end<TElem> (si);
    for (ElemIter iter = lsf.template begin<TElem> (si); iter != iterEnd; ++iter)
    {
      TElem * elem = *iter;
      
    //  get the corner coordinates ans the values of the LSF
      midpoint = 0;
      for (size_t i = 0; i < num_corners; i++)
      {
        Vertex * vrt = elem->vertex (i);
        midpoint += (corners [i] = aaPos [vrt]);
        if (lsf.inner_dof_indices (vrt, 0, ind) != 1)
          UG_THROW ("LSVolume: Not a scalar grid function for the LSF!");
        lsf_values [i] = DoFRef (lsf, ind [0]);
      }
      midpoint /= num_corners;
      
    //  compute the volumes
      number vol_plus, vol_minus;
      LSElementSize<ref_elem_t, dim>::compute (corners, lsf_values, vol_plus, vol_minus);
    
    //  compute the weight
      number density;
      (* pDensity) (density, midpoint, time, si);
      vol_plus *= density; vol_minus *= density;
      
    //  add them to the total volumes
      m_volume_plus += vol_plus; m_volume_minus += vol_minus;
      
    //  add them to the subset
      ss_vol_plus += vol_plus; ss_vol_minus += vol_minus;
      
    //  check the positivity (if requested)
      if (m_check_positivity && (vol_plus < 0 || vol_minus < 0))
      {
        UG_LOG ("Warning in LSVolume: V+ = " << vol_plus << " < 0 or V- = "
            << vol_minus << " < 0\n");
        for (size_t i = 0; i < num_corners; i++)
        {
          UG_LOG (" - corner[" << i << "] = " << corners[i]
            << ", grid_data_index = " << elem->vertex(i)->grid_data_index ()
            << ", lsf = " << lsf_values[i] << '\n');
        }
      }
    }
    m_ss_vol_plus[si] += ss_vol_plus; m_ss_vol_minus[si] += ss_vol_minus;
  }
}
 
#ifdef UG_PARALLEL
template <typename TGridFunc>
void LSVolume<TGridFunc>::allreduce_volumes ()
{
  pcl::ProcessCommunicator procComm;
  m_volume_minus = procComm.allreduce (m_volume_minus, PCL_RO_SUM);
  m_volume_plus = procComm.allreduce (m_volume_plus, PCL_RO_SUM);
  
  const int n_ss = m_spLSF->num_subsets ();
  for (int si = 0; si < n_ss; si++)
  {
    m_ss_vol_minus[si] = procComm.allreduce (m_ss_vol_minus[si], PCL_RO_SUM);
    m_ss_vol_plus[si] = procComm.allreduce (m_ss_vol_plus[si], PCL_RO_SUM);
  }
}
#endif
 
template <typename TGridFunc>
void LSVolume<TGridFunc>::volume_in_subsets
(
  const char * ss_names, 
  number & vol_plus, 
  number & vol_minus 
) const
{
  SubsetGroup ss_grp (m_spLSF->domain()->subset_handler ());
  ss_grp.add (TokenizeString (ss_names));
  
  vol_plus = vol_minus = 0;
  for (size_t i = 0; i < ss_grp.size (); i++)
  {
    int si = ss_grp [i];
    vol_plus += m_ss_vol_plus[si];
    vol_minus += m_ss_vol_minus[si];
  }
}
 
template <typename TGridFunc>
void LSVolume<TGridFunc>::print_details () const
{
  ConstSmartPtr<domain_type> spDom = m_spLSF->domain ();
  
  UG_LOG ("Volumes of the subset:\n" << std::fixed);
  number vol_plus = 0, vol_minus = 0;
  for (int si = 0; si < m_spLSF->num_subsets (); si++)
  {
    vol_plus += m_ss_vol_plus[si]; vol_minus += m_ss_vol_minus[si];
    UG_LOG (spDom->subset_handler()->get_subset_name (si) << " (" << si
          << "): V(+) = " << m_ss_vol_plus[si] << ", V(-) = " << m_ss_vol_minus[si] << "\n");
  }
  UG_LOG ("Totally: V(+) = " << vol_plus << ", V(-) = " << vol_minus << "\n");
}
 
/*---- Class 'LSElementSize': ----*/
 
template <int WDim>
int LSElementSize<ReferenceEdge, WDim>::compute
(
  const MathVector<WDim> * corner, 
  const number * lsf, 
  number & vol_plus, 
  number & vol_minus 
)
{
  number vol_0, vol_1;
  int r;
  
  if (lsf[0] * lsf[1] > 0)
  {
    vol_0 = ElementSize<ref_element_type, dim> (corner);
    vol_1 = 0;
    r = 1;
  }
  else
  {
    number t = lsf[1] / (lsf[1] - lsf[0]);
    MathVector<dim> pnt[3];
    pnt[0] = corner[0];
    VecScaleAdd (pnt[1], t, corner[0], 1 - t, corner[1]);
    pnt[2] = corner[1];
    vol_0 = ElementSize<ref_element_type, dim> (pnt);
    vol_1 = ElementSize<ref_element_type, dim> (pnt + 1);
    r = 0;
  }
  
  if (lsf[0] >= 0)
  {
    vol_plus = vol_0; vol_minus = vol_1;
  }
  else
  {
    vol_minus = vol_0; vol_plus = vol_1; r = - r;
  }
  return r;
}
 
template <int WDim>
int LSElementSize<ReferenceTriangle, WDim>::compute
(
  const MathVector<WDim> * corner, 
  const number * lsf, 
  number & vol_plus, 
  number & vol_minus 
)
{
  number vol = ElementSize<ref_element_type, dim> (corner);
  
//  look for 'positive' and 'negative' corners
  int i_pos = -1, i_neg = -1;
  size_t num_neg = 0;
  for (int i = 0; i < 3; i++)
    if (lsf[i] < 0)
    {
      i_neg = i;
      num_neg++;
    }
    else
      i_pos = i;
    
//  check if all are positive or all are negative
  if (num_neg == 0)
  {
    vol_plus = vol; vol_minus = 0;
    return 1;
  }
  if (num_neg == 3)
  {
    vol_plus = 0; vol_minus = vol;
    return -1;
  }
  
//  look for the cut corner
  int i_0 = (num_neg == 1)? i_neg : i_pos;
  int i_1 = (i_0 + 1) % 3, i_2 = (i_0 + 2) % 3;
 
//  compose a new triangle
  number t;
  MathVector<dim> cut_tri [3];
  cut_tri [0] = corner [i_0];
  t = lsf[i_1] / (lsf[i_1] - lsf[i_0]);
  VecScaleAdd (cut_tri [1], t, corner[i_0], 1 - t, corner[i_1]);
  t = lsf[i_2] / (lsf[i_2] - lsf[i_0]);
  VecScaleAdd (cut_tri [2], t, corner[i_0], 1 - t, corner[i_2]);
  
//  get the cut volume
  number cut_vol = ElementSize<ReferenceTriangle, dim> (cut_tri);
  
//  get the volumes
  if (i_0 == i_neg) // i_0 is either i_pos or i_neg
  {
    vol_minus = cut_vol; vol_plus = vol - cut_vol;
  }
  else
  {
    vol_plus = cut_vol; vol_minus = vol - cut_vol;
  }
  return 0;
}
 
template <int WDim>
int LSElementSize<ReferenceTetrahedron, WDim>::compute
(
  const MathVector<WDim> * corner, 
  const number * lsf, 
  number & vol_plus, 
  number & vol_minus 
)
{
  number vol = ElementSize<ref_element_type, dim> (corner);
  
//  look for 'positive' and 'negative' corners
  int i_pos[2], i_neg[2];
  i_pos[0] = -1; i_pos[1] = -1;
  i_neg[0] = -1; i_neg[1] = -1;
  size_t num_neg = 0;
  for (int i = 0; i < 4; i++)
    if (lsf[i] < 0)
    {
      if (i_neg[0] < 0) i_neg[0] = i; else i_neg[1] = i;
      num_neg++;
    }
    else
    {
      if (i_pos[0] < 0) i_pos[0] = i; else i_pos[1] = i;
    }
  
//  check if all are positive or all are negative
  if (num_neg == 0)
  {
    vol_plus = vol; vol_minus = 0;
    return 1;
  }
  if (num_neg == 4)
  {
    vol_plus = 0; vol_minus = vol;
    return -1;
  }
  
//  if only one corner is cut out
  if (num_neg == 1 || num_neg == 3)
  {
    int i_0 = (num_neg == 1)? i_neg[0] : i_pos[0];
    MathVector<dim> cut_tet [4];
    for (int i = 0; i < 4; i++)
      if (i == i_0)
        cut_tet [i] = corner [i];
      else
      {
        number t = lsf[i] / (lsf[i] - lsf[i_0]);
        VecScaleAdd (cut_tet [i], t, corner[i_0], 1 - t, corner[i]);
      }
    // Note: The orientation of cut_tet should be the same as for the original shape, i.e. correct.
    number cut_vol = ElementSize<ReferenceTetrahedron, dim> (cut_tet);
    if (i_0 == i_neg[0]) // i_0 is either i_pos[0] or i_neg[0]
    {
      vol_minus = cut_vol; vol_plus = vol - cut_vol;
    }
    else
    {
      vol_plus = cut_vol; vol_minus = vol - cut_vol;
    }
    
    return 0;
  }
  
//  if two corners at every side: two prisms
  // We make the triangle containing i_neg[0] to be the botton.
  // To this end, we order i_pos[0] and i_pos[1] in such a way that
  // the new prism is properly oriented.
  int side_idx = tet_rules::FACE_FROM_VRTS [i_neg[0]] [i_pos[0]] [i_pos[1]];
  const int * side_co = tet_rules::FACE_VRT_INDS [side_idx];
  int k;
  for (k = 0; k < 3; k++)
    if (side_co [k] == i_neg[0]) break;
  UG_ASSERT (k < 3, "LSElementSize<ReferenceTetrahedron, WDim>::compute: internal error.");
  i_pos[0] = side_co [(k + 1) % 3];
  i_pos[1] = side_co [(k + 2) % 3];
  
  // Now compute the coordinates of the corners of the prism
  MathVector<dim> neg_prism [6];
  for (int k = 0; k < 2; k++)
  {
    int i_0 = i_neg [k];
    neg_prism [3 * k] = corner [i_0];
    for (int l = 0; l < 2; l++)
    {
      int i = i_pos [l];
      number t = lsf[i] / (lsf[i] - lsf[i_0]);
      VecScaleAdd (neg_prism [3 * k + l + 1], t, corner[i_0], 1 - t, corner[i]);
    }
  }
  vol_minus = ElementSize<ReferencePrism, dim> (neg_prism);
  vol_plus = vol - vol_minus;
  return 0;
}
 
template <int WDim>
int LSElementSize<ReferencePrism, WDim>::compute
(
  const MathVector<WDim> * corner, 
  const number * lsf, 
  number & vol_plus, 
  number & vol_minus 
)
{
  number vol = ElementSize<ref_element_type, dim> (corner);
  
//  look for 'positive' and 'negative' corners
  int i_pos[3], i_neg[3];
  i_pos[0] = -1; i_pos[1] = -1; i_pos[2] = -1;
  i_neg[0] = -1; i_neg[1] = -1; i_neg[2] = -1;
  size_t num_neg = 0;
  for (int i = 0; i < 6; i++)
    if (lsf[i] < 0) // note that we consider 0 as a positive value! This is used below!
    {
      if (i_neg[0] < 0) i_neg[0] = i;
      else if (i_neg[1] < 0) i_neg[1] = i;
      else i_neg[2] = i;
      num_neg++;
    }
    else
    {
      if (i_pos[0] < 0) i_pos[0] = i;
      else if (i_pos[1] < 0) i_pos[1] = i;
      else i_pos[2] = i;
    }
  /*
   * Note that i_neg[] and i_pos[] are now ordered in ascending order.
   */
  
//  check if all are positive or all are negative
  if (num_neg == 0)
  {
    vol_plus = vol; vol_minus = 0;
    return 1;
  }
  if (num_neg == 6)
  {
    vol_plus = 0; vol_minus = vol;
    return -1;
  }
  
//  if only one corner is cut out, i.e. a tetrahedron is cut out
  if (num_neg == 1 || num_neg == 5)
  {
    int i_0 = (num_neg == 1)? i_neg[0] : i_pos[0];
    int shift = (i_0 < 3)? 0 : 3;
    MathVector<dim> cut_tet [4];
    for (int k = 0; k < 3; k++) // loop the corners of the base where i_0 lies
    {
      int i = k + shift;
      if (i == i_0)
        cut_tet [k] = corner [i];
      else
      {
        number t = lsf[i] / (lsf[i] - lsf[i_0]);
        VecScaleAdd (cut_tet [k], t, corner[i_0], 1 - t, corner[i]);
      }
    }
    // For the base of the tetrahedron, we used the same orientation
    // as for the base of the prism. This is OK if that was the bottom
    // base of the prism. For the top base, we should invert the orientation:
    if (shift != 0)
    {
      cut_tet[3] = cut_tet[2]; // cut_tet[3] is here a temporary variable
      cut_tet[2] = cut_tet[1];
      cut_tet[1] = cut_tet[3];
    }
    int i_1 = (i_0 + 3) % 6; // connected corner of the opposite base
    number t_1 = lsf[i_1] / (lsf[i_1] - lsf[i_0]);
    VecScaleAdd (cut_tet [3], t_1, corner[i_0], 1 - t_1, corner[i_1]);
    number cut_vol = ElementSize<ReferenceTetrahedron, dim> (cut_tet);
    if (i_0 == i_neg[0]) // i_0 is either i_pos[0] or i_neg[0]
    {
      vol_minus = cut_vol; vol_plus = vol - cut_vol;
    }
    else
    {
      vol_plus = cut_vol; vol_minus = vol - cut_vol;
    }
    
    return 0;
  }
  
//  if two corners at every side: 2 tetrahedra or a prism cut out
  if (num_neg == 2 || num_neg == 4)
  {
    int i_0, i_1;
    if (num_neg == 2)
    {
      i_0 = i_neg[0]; i_1 = i_neg[1];
    }
    else
    {
      i_0 = i_pos[0]; i_1 = i_pos[1];
    }
    /*  Here, i_0 < i_1 (s. the remark above). */
    
  // Case 1: 2 corners of one base are cut out; this means, a prism is cut out
    if (i_1 < 3 || i_0 >= 3)
    {
      int shift = (i_0 < 3)? 0 : 3;
      int i_2 = -1; // the 3rd corner of the base
      for (int k = 0; k < 3; k++)
      {
        i_2 = k + shift;
        if (i_2 != i_0 && i_2 != i_1) break;
      }
      int i_0a = (i_0 + 3) % 6; // corner connected to i_0 on the opposite base
      int i_1a = (i_1 + 3) % 6; // corner connected to i_1 on the opposite base
      
      MathVector<dim> cut_prism [6];
      number t;
      
      // Remark: We construct the prism in such a way that the 0th
      // side (with corners 0, 1, 4, 3) is planar (to be consistent
      // with the computation of the volume in geometry_util.h).
      // This can be achieved by taking this side from the triangular
      // base of the original prism.
      // Note however that we should provide the correct orientation
      // of the prism. Thus we should provide (i_0 + 1) % 3 == i_1
      // if the bottom base of the original prism is cut, and
      // i_0a == (i_1a + 1) % 3 for the top base. To this end, we swap
      // i_0 with i_1 (and i_0a with i_1a respectively) if necessary.
      if ((i_1 < 3 && (i_0 + 1) % 3 != i_1) || (i_0 >= 3 && i_0a != (i_1a + 1) % 3))
      {
        int i;
        i = i_0; i_0 = i_1; i_1 = i;
        i = i_0a; i_0a = i_1a; i_1a = i;
      }
      
      cut_prism [0] = corner[i_0];
      t = lsf[i_2] / (lsf[i_2] - lsf[i_0]);
      VecScaleAdd (cut_prism [1], t, corner[i_0], 1 - t, corner[i_2]);
      t = lsf[i_0a] / (lsf[i_0a] - lsf[i_0]);
      VecScaleAdd (cut_prism [2], t, corner[i_0], 1 - t, corner[i_0a]);
      
      cut_prism [3] = corner[i_1];
      t = lsf[i_2] / (lsf[i_2] - lsf[i_1]);
      VecScaleAdd (cut_prism [4], t, corner[i_1], 1 - t, corner[i_2]);
      t = lsf[i_1a] / (lsf[i_1a] - lsf[i_1]);
      VecScaleAdd (cut_prism [5], t, corner[i_1], 1 - t, corner[i_1a]);
      
      number cut_vol = ElementSize<ReferencePrism, dim> (cut_prism);
      if (num_neg == 2)
      {
        vol_minus = cut_vol;
        vol_plus = vol - cut_vol;
      }
      else
      {
        vol_plus = cut_vol;
        vol_minus = vol - cut_vol;
      }
      
      return 0;
    }
    
  // Case 2: one corner of each base are cut out, and they are connected; a prism is cut out
    if (i_1 - i_0 == 3)
    {
      // We construct a new (cut out) prism and try to keep the 0th
      // side as planar as possible: see the previous remark. For this,
      // we take a part of a side of the original prism as a side of
      // the new one.
      MathVector<dim> cut_prism [6];
      cut_prism [0] = corner [i_0];
      cut_prism [3] = corner [i_1];
      {
        int i = (i_0 + 1) % 3; int j = i + 3;
        number t;
        t = lsf[i] / (lsf[i] - lsf[i_0]);
        VecScaleAdd (cut_prism [1], t, corner[i_0], 1 - t, corner[i]);
        t = lsf[j] / (lsf[j] - lsf[i_1]);
        VecScaleAdd (cut_prism [4], t, corner[i_1], 1 - t, corner[j]);
      }
      {
        int i = (i_0 + 2) % 3; int j = i + 3;
        number t;
        t = lsf[i] / (lsf[i] - lsf[i_0]);
        VecScaleAdd (cut_prism [2], t, corner[i_0], 1 - t, corner[i]);
        t = lsf[j] / (lsf[j] - lsf[i_1]);
        VecScaleAdd (cut_prism [5], t, corner[i_1], 1 - t, corner[j]);
      }
      number cut_vol = ElementSize<ReferencePrism, dim> (cut_prism);
      if (num_neg == 2)
      {
        vol_minus = cut_vol;
        vol_plus = vol - cut_vol;
      }
      else
      {
        vol_plus = cut_vol;
        vol_minus = vol - cut_vol;
      }
      
      return 0;
    }
    
  //  Case 3: The corners belong to different bases but are not connected.
  //  We consider two separated tetrahedra. They do not hinder each other.
  //  We call the function recursively.
    number tmp_lsf [6], tmp_vol;
    memcpy (tmp_lsf, lsf, 6 * sizeof (number));
    if (num_neg == 2)
    {
    //  2 negative corners
      number neg_vol;
      
      tmp_lsf [i_neg[0]] = 1; // correct the first negative corner
      LSElementSize<ReferencePrism, dim>::compute (corner, tmp_lsf, tmp_vol, neg_vol);
      vol_minus = neg_vol;
      
      tmp_lsf [i_neg[0]] = lsf [i_neg[0]];
      tmp_lsf [i_neg[1]] = 1; // correct the second negative corner
      LSElementSize<ReferencePrism, dim>::compute (corner, tmp_lsf, tmp_vol, neg_vol);
      vol_minus += neg_vol;
      
      vol_plus = vol - vol_minus;
      return 0;
    }
    if (num_neg == 4)
    {
    //  2 positive corners
      number pos_vol;
      
      tmp_lsf [i_pos[0]] = -1; // correct the first positive corner
      LSElementSize<ReferencePrism, dim>::compute (corner, tmp_lsf, pos_vol, tmp_vol);
      vol_plus = pos_vol;
      
      tmp_lsf [i_pos[0]] = lsf [i_pos[0]];
      tmp_lsf [i_pos[1]] = -1; // correct the second positive corner
      LSElementSize<ReferencePrism, dim>::compute (corner, tmp_lsf, pos_vol, tmp_vol);
      vol_plus += pos_vol;
      
      vol_minus = vol - vol_plus;
      return 0;
    }
  }
  
//  if three corners at every side: two prisms or "a very complicated case"
  if ((i_neg[0] == 0 && i_neg[2] == 2) || (i_neg[0] == 3 && i_neg[2] == 5))
  {
  //  Case 1: 2 prisms: every base has its own sign
    number cut_vol;
    MathVector<dim> cut_prism [6];
    for (int i = 0; i < 3; i++) cut_prism [i] = corner [i];
    for (int i = 0; i < 3; i++)
    {
      int j = i + 3;
      number t = lsf[j] / (lsf[j] - lsf[i]);
      VecScaleAdd (cut_prism [j], t, corner [i], 1 - t, corner [j]);
    }
    cut_vol = ElementSize<ReferencePrism, dim> (cut_prism);
    
    if (i_neg[0] == 0)
    {
      vol_minus = cut_vol;
      vol_plus = vol - cut_vol;
    }
    else
    {
      vol_plus = cut_vol;
      vol_minus = vol - cut_vol;
    }
    
    return 0;
  }
  
  // Case 2: The "very complicated case": Two corners of one base and one
  // of the other one are cut.
  // Only one of the three sides connecting the bases is shared by
  // a negative and a positive corner. Denote this side by [A, B], and its
  // point, where the LSF is 0, by S. (Without loss of generality, we
  // assume that LSF[A] < 0.) We divide the whole prism into two parts:
  // a smaller prism and a tetrahedron. The smaller prism is obtained from
  // the original one by shifting the corner A into S. The tetrahedron is
  // the difference of these prisms. After that, we call the function
  // recursively to compute the volumes for the smaller prism. It has
  // now one negative corner less than the original one (becase 0 is
  // considered as a positive number here), so that this case is implemented
  // above.
  MathVector<dim> part_prism [6];
  number part_lsf [6];
  int base_i = -1;
  
  for (int i = 0; i < 3; i++)
  {
    int j = i + 3;
    part_prism [i] = corner [i]; part_lsf [i] = lsf [i];
    part_prism [j] = corner [j]; part_lsf [j] = lsf [j];
    if ((lsf[i] >= 0 && lsf[j] < 0)
        || (lsf[i] < 0 && lsf[j] >= 0)) // note: 0 is considered as a positive value
    {
      base_i = (lsf[i] < 0)? i : j; // the corner to shift
      part_lsf[base_i] = 0; // we make it to be a positive corner!
      number t = lsf[j] / (lsf[j] - lsf[i]);
      VecScaleAdd (part_prism [base_i], t, corner [i], 1 - t, corner [j]);
    }
  }
  UG_ASSERT (base_i >= 0, "LSElementSize<ReferencePrism, dim>::compute: internal error");
  number part_prism_vol_plus, part_prism_vol_minus;
  LSElementSize<ReferencePrism, dim>::compute (part_prism, part_lsf,
    part_prism_vol_plus, part_prism_vol_minus);
  
  MathVector<dim> part_tet [4];
  int shift = (base_i < 3)? 0 : 3;
  for (int k = 0; k < 3; k++)
  {
    part_tet [k] = corner [k + shift];
    part_lsf [k] = lsf [k + shift];
  }
  // As before, we must provide the correct orientation of the tetrahedron.
  // To this end, we change the current orientation, if base_i is on the
  // top base of the prism. (Cf. the case of one cut corner.)
  if (shift != 0)
  {
    part_tet[3] = part_tet[2]; part_lsf[3] = part_lsf[2];
    part_tet[2] = part_tet[1]; part_lsf[2] = part_lsf[1];
    part_tet[1] = part_tet[3]; part_lsf[1] = part_lsf[3];
  }
  part_tet [3] = part_prism [base_i];
  part_lsf [3] = 0;
  number part_tet_vol_plus, part_tet_vol_minus;
  LSElementSize<ReferenceTetrahedron, dim>::compute (part_tet, part_lsf,
    part_tet_vol_plus, part_tet_vol_minus);
  
  vol_plus = part_prism_vol_plus + part_tet_vol_plus;
  vol_minus = part_prism_vol_minus + part_tet_vol_minus;
  
  return 0;
}
 
} // end namespace LevelSet
} // end namespace ug
 
/* End of File */