geo_seq_check.cc

Go to the documentation of this file.

#include "rheolef/geo.h"
#include "rheolef/geo_domain.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/space_numbering.h"
 
#include "geo_header.h"
 
namespace rheolef {
 
// TODO: quadrangle may also be convex
// ----------------------------------------------------------------------------
// compute measure
// ----------------------------------------------------------------------------
template <class T>
static
T
meas_t (
  const point_basic<T>& a,
  const point_basic<T>& b,
  const point_basic<T>& c)
{
  return vect2d(b-a, c-a)/2.0;
}
template <class T>
static
T
meas_T (
  const point_basic<T>& a,
  const point_basic<T>& b,
  const point_basic<T>& c,
  const point_basic<T>& d)
{
  return mixt(b-a, c-a, d-a)/6;
}
template <class T>
static
void
meas_P (
  const point_basic<T>& a0,
  const point_basic<T>& a1,
  const point_basic<T>& a2,
  const point_basic<T>& a3,
  const point_basic<T>& a4,
  const point_basic<T>& a5,
  T& meas1,
  T& meas2,
  T& meas3)
{
  meas1 = meas_T (a0, a1, a2, a3);
  meas2 = meas_T (a1, a2, a3, a4);
  meas3 = meas_T (a2, a3, a4, a5);
}
template <class T>
static
bool
check_element (const geo_element& K, const disarray<point_basic<T>,sequential>& p, bool verbose = false)
{
  std::cerr << std::setprecision(std::numeric_limits<T>::digits10);
  bool orient = true;
  switch (K.variant()) {
    case reference_element::p:
      return true;
    case reference_element::e: {
      T meas = p[K[1]][0] - p[K[0]][0];
      orient = (meas > 0);
      if (verbose && !orient) {
        warning_macro("meas(K)="<<meas<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
      }
      break;
    }
    case reference_element::t: {
      T meas = meas_t (p[K[0]], p[K[1]], p[K[2]]);
      orient = (meas > 0);
      if (verbose && !orient) {
        warning_macro("meas(K)="<<meas<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
      }
      break;
    }
    case reference_element::q: {
      T meas1 = meas_t (p[K[0]], p[K[1]], p[K[2]]);  // t=012
      T meas2 = meas_t (p[K[0]], p[K[2]], p[K[3]]);  // t=023
      orient = (meas1 > 0 && meas2 > 0);
      if (verbose && !orient) {
        warning_macro("meas(K)="<<meas1<<"+"<<meas2<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
      }
      break;
    }
    case reference_element::T: {
      T meas = meas_T (p[K[0]], p[K[1]], p[K[2]], p[K[3]]);
      orient = (meas > 0);
      if (verbose && !orient) {
        warning_macro("meas(K)="<<meas<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
      }
      break;
    }
    case reference_element::P: {
      T meas1, meas2, meas3;
      meas_P (p[K[0]], p[K[1]], p[K[2]], p[K[3]], p[K[4]], p[K[5]], meas1, meas2, meas3);
      orient = (meas1 > 0 && meas2 > 0 && meas3 > 0);
      if (verbose && !orient) {
        warning_macro("meas(K)="<<meas1<<"+"<<meas2<<"+"<<meas3
                <<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
      }
      break;
    }
    case reference_element::H: {
      T meas1, meas2, meas3, meas4, meas5, meas6;
      meas_P (p[K[0]], p[K[1]], p[K[2]], p[K[4]], p[K[5]], p[K[6]], meas1, meas2, meas3);
      meas_P (p[K[0]], p[K[2]], p[K[3]], p[K[4]], p[K[6]], p[K[7]], meas4, meas5, meas6);
      orient = (meas1 > 0 && meas2 > 0 && meas3 > 0 && meas4 > 0 && meas5 > 0 && meas6 > 0);
      if (verbose && !orient) {
        warning_macro("meas(K)="<<meas1<<"+"<<meas2<<"+"<<meas3
                <<"+"<<meas4<<"+"<<meas5<<"+"<<meas6
                <<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
      }
      break;
    }
    default:
      error_macro ("unexpected element variant `"<<K.variant()<<"'");
  }
  bool is_convex = true; // TODO: q,P,H
  return (orient && is_convex);
}
// ----------------------------------------------------------------------------
// check element orientation
// ----------------------------------------------------------------------------
template <class T>
bool
geo_rep<T,sequential>::check(bool verbose) const
{
  if (base::_dimension != base::_gs._map_dimension) return true;
  bool status = true;
  for (const_iterator iter = begin(base::_gs._map_dimension), last = end(base::_gs._map_dimension); iter != last; iter++) {
    const geo_element& K = *iter;
    status = status && check_element (K, base::_node, verbose);
  }
  return status;
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
bool
geo_rep<T,distributed>::check(bool verbose) const
{
  // TODO
  return true;
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class geo_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class geo_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef