mirror of https://github.com/CGAL/cgal
format
This commit is contained in:
Abdelkrim Mebarki
parent
01c80c0e4c
commit
5871846965
|
|
@ -27,9 +27,10 @@ CGAL_BEGIN_NAMESPACE
|
|||
|
||||
// The class Euler_integrator_2 is a model of the concept Integrator
|
||||
template <class VectorField_2>
|
||||
class Euler_integrator_2{
|
||||
class Euler_integrator_2
|
||||
{
|
||||
public:
|
||||
typedef Euler_integrator_2<VectorField_2> self;
|
||||
typedef Euler_integrator_2<VectorField_2> self;
|
||||
typedef typename VectorField_2::FT FT;
|
||||
typedef typename VectorField_2::Point_2 Point_2;
|
||||
typedef typename VectorField_2::Vector_2 Vector_2;
|
||||
|
|
@ -37,17 +38,32 @@ public:
|
|||
protected:
|
||||
FT default_integration_step;
|
||||
public:
|
||||
Euler_integrator_2();
|
||||
Euler_integrator_2();
|
||||
|
||||
Euler_integrator_2(const FT & integration_step);
|
||||
inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const bool & index) const;
|
||||
inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, const bool & index) const;
|
||||
inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, Vector_2 v, const bool & index) const;
|
||||
inline FT get_default_integration_step(){return default_integration_step;}
|
||||
// just for debugging
|
||||
inline FT distance(const Point_2 & p, const Point_2 & q)
|
||||
{
|
||||
return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y())));
|
||||
}
|
||||
|
||||
inline Point_2 operator()(const Point_2 & p, const Vector_field_2 &
|
||||
vector_field_2, const bool & index) const;
|
||||
|
||||
inline Point_2 operator()(const Point_2 & p, const Vector_field_2 &
|
||||
vector_field_2, const FT &
|
||||
integration_step, const bool & index)
|
||||
const;
|
||||
|
||||
inline Point_2 operator()(const Point_2 & p, const Vector_field_2 &
|
||||
vector_field_2, const FT &
|
||||
integration_step, Vector_2 v, const bool &
|
||||
index) const;
|
||||
|
||||
inline FT get_default_integration_step()
|
||||
{
|
||||
return default_integration_step;
|
||||
}
|
||||
// just for debugging
|
||||
inline FT distance(const Point_2 & p, const Point_2 & q)
|
||||
{
|
||||
return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y())));
|
||||
}
|
||||
};
|
||||
|
||||
template <class Vector_field>
|
||||
|
|
|
|||
|
|
@ -31,9 +31,10 @@ CGAL_BEGIN_NAMESPACE
|
|||
// bilinear interpolation to extract the vector field values
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
class Regular_grid_2{
|
||||
class Regular_grid_2
|
||||
{
|
||||
public:
|
||||
typedef Regular_grid_2<StreamLinesTraits_2> Vector_field_2;
|
||||
typedef Regular_grid_2<StreamLinesTraits_2> Vector_field_2;
|
||||
typedef StreamLinesTraits_2 Geom_traits;
|
||||
typedef typename StreamLinesTraits_2::FT FT;
|
||||
typedef typename StreamLinesTraits_2::Point_2 Point_2;
|
||||
|
|
@ -50,68 +51,107 @@ protected:
|
|||
bool is_in_samples(const int & i,const int & j) const;
|
||||
public:
|
||||
Regular_grid_2(const int & m, const int & n,const FT & x, const FT & y);
|
||||
// Regular_grid_2();
|
||||
~Regular_grid_2(){delete [] vector_field;}
|
||||
void fill(std::ifstream & f, const FT & x, const FT & y);
|
||||
// Regular_grid_2();
|
||||
~Regular_grid_2()
|
||||
{
|
||||
delete [] vector_field;
|
||||
}
|
||||
void fill(std::ifstream & f, const FT & x, const FT & y);
|
||||
inline typename Geom_traits::Iso_rectangle_2 bbox() const;
|
||||
|
||||
std::pair<Vector_2,FT>
|
||||
get_field(const Point_2 & p) const
|
||||
{
|
||||
CGAL_streamlines_precondition(is_in_domain(p));
|
||||
Vector_2 v = get_vector_field(p);
|
||||
FT density = get_density_field(p);
|
||||
return std::pair<Vector_2, FT>(v,density);
|
||||
}
|
||||
{
|
||||
CGAL_streamlines_precondition(is_in_domain(p));
|
||||
Vector_2 v = get_vector_field(p);
|
||||
FT density = get_density_field(p);
|
||||
return std::pair<Vector_2, FT>(v,density);
|
||||
}
|
||||
|
||||
inline bool is_in_domain(const Point_2 & p) const;
|
||||
inline FT get_integration_step(const Point_2 &) const;
|
||||
inline FT get_integration_step() const;
|
||||
inline FT get_integration_step() const;
|
||||
inline void set_field(const int & i, const int & j, const Vector_2 & v);
|
||||
inline std::pair<int, int> get_dimension(){return std::pair<int, int>(number_of_samples_x, number_of_samples_y);}
|
||||
inline std::pair<FT, FT> get_size(){return std::pair<FT, FT>(domain_size_x, domain_size_y);}};
|
||||
inline std::pair<int, int> get_dimension()
|
||||
{
|
||||
return std::pair<int, int>(number_of_samples_x, number_of_samples_y);
|
||||
}
|
||||
inline std::pair<FT, FT> get_size()
|
||||
{
|
||||
return std::pair<FT,
|
||||
FT>(domain_size_x,
|
||||
domain_size_y);
|
||||
}
|
||||
};
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
inline
|
||||
typename Regular_grid_2<StreamLinesTraits_2>::Geom_traits::Iso_rectangle_2
|
||||
Regular_grid_2<StreamLinesTraits_2>::bbox() const{
|
||||
return typename Geom_traits::Iso_rectangle_2(0.0, 0.0, domain_size_x, domain_size_y);}
|
||||
Regular_grid_2<StreamLinesTraits_2>::bbox() const
|
||||
{
|
||||
return typename Geom_traits::Iso_rectangle_2(0.0, 0.0,
|
||||
domain_size_x,
|
||||
domain_size_y);
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
inline int
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_index(const int & i, const int & j) const{
|
||||
return 2*(number_of_samples_x*j + i);}
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_index(const int & i, const
|
||||
int & j) const
|
||||
{
|
||||
return 2*(number_of_samples_x*j + i);
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
Regular_grid_2<StreamLinesTraits_2>::Regular_grid_2(const int & m, const int & n,const FT & x, const FT & y){
|
||||
Regular_grid_2<StreamLinesTraits_2>::Regular_grid_2(const int & m,
|
||||
const int &
|
||||
n,const FT & x,
|
||||
const FT & y)
|
||||
{
|
||||
number_of_samples_x = m;
|
||||
number_of_samples_y = n;
|
||||
domain_size_x = x;
|
||||
domain_size_y = y;
|
||||
vector_field = new FT[number_of_samples_x*number_of_samples_y* 2];}
|
||||
vector_field = new FT[number_of_samples_x*number_of_samples_y* 2];
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
inline void
|
||||
Regular_grid_2<StreamLinesTraits_2>::set_field(const int & i, const int & j, const Vector_2 & v){
|
||||
Regular_grid_2<StreamLinesTraits_2>::set_field(const int & i, const
|
||||
int & j, const Vector_2
|
||||
& v)
|
||||
{
|
||||
CGAL_streamlines_precondition(is_in_samples(i,j));
|
||||
int index = get_index(i,j);
|
||||
vector_field[index++] = v.x();
|
||||
vector_field[index] = v.y();}
|
||||
vector_field[index] = v.y();
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
inline bool
|
||||
Regular_grid_2<StreamLinesTraits_2>::is_in_domain(const Point_2 & p) const{
|
||||
return ((p.x()>=0.0) && (p.x()<=domain_size_x) && (p.y()>=0.0) && (p.y()<=domain_size_y));}
|
||||
Regular_grid_2<StreamLinesTraits_2>::is_in_domain(const Point_2 & p)
|
||||
const
|
||||
{
|
||||
return ((p.x()>=0.0) && (p.x()<=domain_size_x) && (p.y()>=0.0) &&
|
||||
(p.y()<=domain_size_y));
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
bool
|
||||
Regular_grid_2<StreamLinesTraits_2>::is_in_samples(const int & i, const int & j) const{
|
||||
return ((i>=0) && (i<=number_of_samples_x-1) && (j>=0) && (j<=number_of_samples_y-1));}
|
||||
Regular_grid_2<StreamLinesTraits_2>::is_in_samples(const int & i,
|
||||
const int & j)
|
||||
const
|
||||
{
|
||||
return ((i>=0) && (i<=number_of_samples_x-1) && (j>=0) &&
|
||||
(j<=number_of_samples_y-1));
|
||||
}
|
||||
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
typename Regular_grid_2<StreamLinesTraits_2>::Vector_2
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_vector_field(const Point_2 & p) const{
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_vector_field(const Point_2 &
|
||||
p) const
|
||||
{
|
||||
FT fXv,fYv;
|
||||
FT x = (p.x() / domain_size_x) * (number_of_samples_x-1);
|
||||
FT y = (p.y() / domain_size_y) * (number_of_samples_y-1);
|
||||
|
|
@ -146,24 +186,34 @@ Regular_grid_2<StreamLinesTraits_2>::get_vector_field(const Point_2 & p) const{
|
|||
fXv = fXv / normal;
|
||||
fYv = fYv / normal;
|
||||
Vector_2 v = Vector_2(fXv, fYv);
|
||||
return v;}
|
||||
return v;
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
typename Regular_grid_2<StreamLinesTraits_2>::FT
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_density_field(const Point_2 & p) const{
|
||||
return 1.0;}
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_density_field(const Point_2 &
|
||||
p) const
|
||||
{
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
template<class StreamLinesTraits_2>
|
||||
inline
|
||||
typename Regular_grid_2<StreamLinesTraits_2>::FT
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_integration_step(const Point_2 &) const{
|
||||
return 1.0;}
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_integration_step(const
|
||||
Point_2 &)
|
||||
const
|
||||
{
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
template<class StreamLinesTraits_2>
|
||||
inline
|
||||
typename Regular_grid_2<StreamLinesTraits_2>::FT
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_integration_step() const{
|
||||
return 1.0;}
|
||||
Regular_grid_2<StreamLinesTraits_2>::get_integration_step() const
|
||||
{
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
|
|
|
|||
|
|
@ -28,36 +28,55 @@ CGAL_BEGIN_NAMESPACE
|
|||
|
||||
// The class Runge_kutta_integrator_2 is a model of the concept Integrator
|
||||
template <class VectorField_2>
|
||||
class Runge_kutta_integrator_2{
|
||||
class Runge_kutta_integrator_2
|
||||
{
|
||||
public:
|
||||
typedef Runge_kutta_integrator_2<VectorField_2> self;
|
||||
typedef Runge_kutta_integrator_2<VectorField_2> self;
|
||||
typedef typename VectorField_2::FT FT;
|
||||
typedef typename VectorField_2::Point_2 Point_2;
|
||||
typedef typename VectorField_2::Vector_2 Vector_2;
|
||||
typedef typename VectorField_2::Vector_field_2 Vector_field_2;
|
||||
protected:
|
||||
typedef CGAL::Euler_integrator_2<VectorField_2> Euler_integrator_2;
|
||||
typedef CGAL::Euler_integrator_2<VectorField_2> Euler_integrator_2;
|
||||
|
||||
Euler_integrator_2 * euler_integrator_2;
|
||||
|
||||
FT default_integration_step;
|
||||
|
||||
public:
|
||||
Runge_kutta_integrator_2();
|
||||
|
||||
Runge_kutta_integrator_2(FT integration_step);
|
||||
|
||||
~Runge_kutta_integrator_2();
|
||||
|
||||
Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const bool & index) const;
|
||||
|
||||
Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, const bool & index) const;
|
||||
|
||||
Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, Vector_2 v, const bool & index) const;
|
||||
|
||||
inline FT get_default_integration_step(){return default_integration_step;}
|
||||
Euler_integrator_2 * get_euler_integrator_2(){return euler_integrator_2;}
|
||||
inline FT distance(const Point_2 & p, const Point_2 & q){
|
||||
return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y())));}
|
||||
inline FT get_default_integration_step()
|
||||
{
|
||||
return default_integration_step;
|
||||
}
|
||||
Euler_integrator_2 * get_euler_integrator_2()
|
||||
{
|
||||
return euler_integrator_2;
|
||||
}
|
||||
inline FT distance(const Point_2 & p, const Point_2 & q)
|
||||
{
|
||||
return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y()
|
||||
-
|
||||
q.y())*(p.y() - q.y())));
|
||||
}
|
||||
};
|
||||
|
||||
template <class VectorField_2>
|
||||
Runge_kutta_integrator_2<VectorField_2>::
|
||||
Runge_kutta_integrator_2()
|
||||
{
|
||||
euler_integrator_2 = new Euler_integrator_2();
|
||||
euler_integrator_2 = new Euler_integrator_2();
|
||||
default_integration_step =
|
||||
euler_integrator_2->get_default_integration_step();
|
||||
}
|
||||
|
|
@ -66,7 +85,7 @@ template <class VectorField_2>
|
|||
Runge_kutta_integrator_2<VectorField_2>::
|
||||
Runge_kutta_integrator_2(FT integration_step)
|
||||
{
|
||||
euler_integrator_2 = new Euler_integrator_2(integration_step);
|
||||
euler_integrator_2 = new Euler_integrator_2(integration_step);
|
||||
default_integration_step =
|
||||
euler_integrator_2->get_default_integration_step();
|
||||
}
|
||||
|
|
@ -75,7 +94,7 @@ template <class VectorField_2>
|
|||
Runge_kutta_integrator_2<VectorField_2>::
|
||||
~Runge_kutta_integrator_2()
|
||||
{
|
||||
delete euler_integrator_2;
|
||||
delete euler_integrator_2;
|
||||
}
|
||||
|
||||
template <class VectorField_2>
|
||||
|
|
@ -95,7 +114,9 @@ Point_2 Runge_kutta_integrator_2<VectorField_2>::operator()
|
|||
template <class VectorField_2>
|
||||
inline
|
||||
typename Runge_kutta_integrator_2<VectorField_2>::Point_2 Runge_kutta_integrator_2<VectorField_2>::operator()
|
||||
(const Point_2 & p, const Vector_field_2& vector_field_2, const FT & integration_step, const bool & index) const{
|
||||
(const Point_2 & p, const Vector_field_2& vector_field_2, const FT &
|
||||
integration_step, const bool & index) const
|
||||
{
|
||||
Vector_2 v;
|
||||
v = vector_field_2.get_field(p).first;
|
||||
Runge_kutta_integrator_2<VectorField_2>
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -20,19 +20,15 @@
|
|||
#ifndef CGAL_TRIANGULAR_FIELD_2_H_
|
||||
#define CGAL_TRIANGULAR_FIELD_2_H_
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/Cartesian.h>
|
||||
#include <CGAL/Delaunay_triangulation_2.h>
|
||||
#include <CGAL/Polygon_2_algorithms.h>
|
||||
|
||||
#include <float.h>
|
||||
// #include <strings.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <list>
|
||||
#include <queue>
|
||||
#include <cmath>
|
||||
#include <string>
|
||||
|
||||
|
|
@ -41,143 +37,190 @@
|
|||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
class Triangular_field_2{
|
||||
class Triangular_field_2
|
||||
{
|
||||
public:
|
||||
typedef Triangular_field_2<StreamLinesTraits_2> Vector_field_2;
|
||||
typedef StreamLinesTraits_2 Geom_traits;
|
||||
typedef typename StreamLinesTraits_2::FT FT;
|
||||
typedef typename StreamLinesTraits_2::Point_2 Point_2;
|
||||
typedef typename StreamLinesTraits_2::Vector_2 Vector_2;
|
||||
typedef Triangular_field_2<StreamLinesTraits_2> Vector_field_2;
|
||||
typedef StreamLinesTraits_2 Geom_traits;
|
||||
typedef typename StreamLinesTraits_2::FT FT;
|
||||
typedef typename StreamLinesTraits_2::Point_2 Point_2;
|
||||
typedef typename StreamLinesTraits_2::Vector_2 Vector_2;
|
||||
protected:
|
||||
typedef CGAL::Triangulation_vertex_base_2<StreamLinesTraits_2> Vb;
|
||||
typedef CGAL::Triangulation_face_base_with_info_2<Vector_2, StreamLinesTraits_2> Fb;
|
||||
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
|
||||
typedef CGAL::Delaunay_triangulation_2<StreamLinesTraits_2,TDS> D_Ttr;
|
||||
typedef typename D_Ttr::Vertex_handle Vertex_handle;
|
||||
typedef typename TDS::Face_handle Face_handle;
|
||||
typedef CGAL::Triangulation_vertex_base_2<StreamLinesTraits_2> Vb;
|
||||
typedef CGAL::Triangulation_face_base_with_info_2<Vector_2, StreamLinesTraits_2> Fb;
|
||||
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
|
||||
typedef CGAL::Delaunay_triangulation_2<StreamLinesTraits_2,TDS> D_Ttr;
|
||||
typedef typename D_Ttr::Vertex_handle Vertex_handle;
|
||||
typedef typename TDS::Face_handle Face_handle;
|
||||
public:
|
||||
typedef typename D_Ttr::Vertex_iterator Vertex_iterator;
|
||||
D_Ttr m_D_Ttr;
|
||||
typedef typename D_Ttr::Vertex_iterator Vertex_iterator;
|
||||
D_Ttr m_D_Ttr;
|
||||
protected:
|
||||
Vector_2 get_vector_field(const Point_2 & p) const;
|
||||
FT get_density_field(const Point_2 & p) const;
|
||||
template <class PointInputIterator, class VectorInputIterator>
|
||||
void fill(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin){
|
||||
std::cout << "reading file...\n";
|
||||
while(pBegin != pEnd){
|
||||
Point_2 p;
|
||||
Vector_2 v;
|
||||
p = (*pBegin);
|
||||
v = (*vBegin);
|
||||
Vertex_handle m_Vertex_handle = m_D_Ttr.insert(p);
|
||||
field_map[p] = v;
|
||||
if (m_D_Ttr.number_of_vertices() == 1){
|
||||
maxx = minx = p.x();
|
||||
maxy = miny = p.y();}
|
||||
if(p.x()<minx)
|
||||
minx = p.x();
|
||||
if(p.y()<miny)
|
||||
miny = p.y();
|
||||
if(p.x()>maxx)
|
||||
maxx = p.x();
|
||||
if(p.y()>maxy)
|
||||
maxy = p.y();
|
||||
pBegin++;
|
||||
vBegin++;}
|
||||
std::cout << "number of samples " << m_D_Ttr.number_of_vertices() << "\n";}
|
||||
public:
|
||||
template <class PointInputIterator, class VectorInputIterator>
|
||||
Triangular_field_2(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin){
|
||||
fill(pBegin, pEnd, vBegin);}
|
||||
inline typename Geom_traits::Iso_rectangle_2 bbox() const;
|
||||
Vector_2 get_vector_field(const Point_2 & p) const;
|
||||
|
||||
std::pair<Vector_2,FT> get_field(const Point_2 & p) const
|
||||
{
|
||||
assert(is_in_domain(p));
|
||||
Vector_2 v = get_vector_field(p);
|
||||
FT density = get_density_field(p);
|
||||
std::pair<Vector_2, FT> field_value(v,density);
|
||||
return field_value;
|
||||
FT get_density_field(const Point_2 & p) const;
|
||||
|
||||
template <class PointInputIterator, class VectorInputIterator>
|
||||
void fill(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin)
|
||||
{
|
||||
std::cout << "reading file...\n";
|
||||
while(pBegin != pEnd)
|
||||
{
|
||||
Point_2 p;
|
||||
Vector_2 v;
|
||||
p = (*pBegin);
|
||||
v = (*vBegin);
|
||||
Vertex_handle m_Vertex_handle = m_D_Ttr.insert(p);
|
||||
field_map[p] = v;
|
||||
if (m_D_Ttr.number_of_vertices() == 1)
|
||||
{
|
||||
maxx = minx = p.x();
|
||||
maxy = miny = p.y();
|
||||
}
|
||||
if(p.x()<minx)
|
||||
{
|
||||
minx = p.x();
|
||||
}
|
||||
if(p.y()<miny)
|
||||
{
|
||||
miny = p.y();
|
||||
}
|
||||
if(p.x()>maxx)
|
||||
{
|
||||
maxx = p.x();
|
||||
}
|
||||
if(p.y()>maxy)
|
||||
{
|
||||
maxy = p.y();
|
||||
}
|
||||
++pBegin;
|
||||
++vBegin;
|
||||
}
|
||||
std::cout << "number of samples " << m_D_Ttr.number_of_vertices() << "\n";
|
||||
}
|
||||
public:
|
||||
template <class PointInputIterator, class VectorInputIterator>
|
||||
Triangular_field_2(PointInputIterator pBegin, PointInputIterator
|
||||
pEnd, VectorInputIterator vBegin)
|
||||
{
|
||||
fill(pBegin, pEnd, vBegin);
|
||||
}
|
||||
|
||||
inline typename Geom_traits::Iso_rectangle_2 bbox() const;
|
||||
|
||||
std::pair<Vector_2,FT> get_field(const Point_2 & p) const
|
||||
{
|
||||
assert(is_in_domain(p));
|
||||
Vector_2 v = get_vector_field(p);
|
||||
FT density = get_density_field(p);
|
||||
return std::make_pair(v, density);
|
||||
}
|
||||
|
||||
bool is_in_domain(const Point_2 & p) const;
|
||||
FT get_integration_step(const Point_2 &) const;
|
||||
FT get_integration_step() const;
|
||||
bool is_in_domain(const Point_2 & p) const;
|
||||
|
||||
FT get_integration_step(const Point_2 &) const;
|
||||
|
||||
FT get_integration_step() const;
|
||||
protected:
|
||||
FT minx;
|
||||
FT miny;
|
||||
FT maxx;
|
||||
FT maxy;
|
||||
FT minx;
|
||||
FT miny;
|
||||
FT maxx;
|
||||
FT maxy;
|
||||
protected:
|
||||
mutable std::map<Point_2, Vector_2> field_map;
|
||||
FT distance(const Point_2 & p, const Point_2 & q){
|
||||
return sqrt(((p.x() - q.x()) * (p.x() - q.x())) + ((p.y() - q.y()) * (p.y() - q.y())));}
|
||||
mutable std::map<Point_2, Vector_2> field_map;
|
||||
|
||||
FT distance(const Point_2 & p, const Point_2 & q)
|
||||
{
|
||||
return sqrt( CGAL::squared_distance(p, q) );
|
||||
}
|
||||
};
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
inline
|
||||
typename Triangular_field_2<StreamLinesTraits_2>::Geom_traits::Iso_rectangle_2
|
||||
Triangular_field_2<StreamLinesTraits_2>::bbox() const{
|
||||
return typename Geom_traits::Iso_rectangle_2(minx, miny, maxx, maxy);}
|
||||
Triangular_field_2<StreamLinesTraits_2>::bbox() const
|
||||
{
|
||||
return typename Geom_traits::Iso_rectangle_2(minx, miny, maxx,
|
||||
maxy);
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
bool
|
||||
Triangular_field_2<StreamLinesTraits_2>::is_in_domain(const Point_2 & p) const{
|
||||
Face_handle f = m_D_Ttr.locate(p);
|
||||
return !m_D_Ttr.is_infinite(f);}
|
||||
Triangular_field_2<StreamLinesTraits_2>::is_in_domain(const Point_2 &
|
||||
p) const
|
||||
{
|
||||
Face_handle f = m_D_Ttr.locate(p);
|
||||
return !m_D_Ttr.is_infinite(f);
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
typename Triangular_field_2<StreamLinesTraits_2>::Vector_2
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_vector_field(const Point_2 & p) const{
|
||||
Face_handle m_Face_handle = m_D_Ttr.locate(p);
|
||||
assert(is_in_domain(p));
|
||||
Point_2 p0 = m_Face_handle->vertex(0)->point();
|
||||
Point_2 p1 = m_Face_handle->vertex(1)->point();
|
||||
Point_2 p2 = m_Face_handle->vertex(2)->point();
|
||||
Vertex_handle v0 = m_Face_handle->vertex(0);
|
||||
Vertex_handle v1 = m_Face_handle->vertex(1);
|
||||
Vertex_handle v2 = m_Face_handle->vertex(2);
|
||||
FT s0,s1,s2,s;
|
||||
std::vector<Point_2> vec;
|
||||
vec.push_back(p0); vec.push_back(p1); vec.push_back(p2);
|
||||
s = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
vec.push_back(p); vec.push_back(p1); vec.push_back(p2);
|
||||
s0 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
vec.push_back(p0); vec.push_back(p); vec.push_back(p2);
|
||||
s1 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
vec.push_back(p0); vec.push_back(p1); vec.push_back(p);
|
||||
s2 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
s0 = s0 / s; s1 = s1 / s; s2 = s2 / s;
|
||||
Vector_2 v_0 = field_map[p0];
|
||||
Vector_2 v_1 = field_map[p1];
|
||||
Vector_2 v_2 = field_map[p2];
|
||||
FT x = ((v_0.x()*s0)+(v_1.x()*s1)+(v_2.x()*s2));
|
||||
FT y = ((v_0.y()*s0)+(v_1.y()*s1)+(v_2.y()*s2));
|
||||
FT normal = sqrt((x)*(x) + (y)*(y));
|
||||
if (normal != 0){
|
||||
x = x / normal;
|
||||
y = y / normal;}
|
||||
Vector_2 v = Vector_2(x, y);
|
||||
return v;}
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_vector_field(const
|
||||
Point_2 & p)
|
||||
const
|
||||
{
|
||||
Face_handle m_Face_handle = m_D_Ttr.locate(p);
|
||||
assert(is_in_domain(p));
|
||||
Vertex_handle v0 = m_Face_handle->vertex(0);
|
||||
Vertex_handle v1 = m_Face_handle->vertex(1);
|
||||
Vertex_handle v2 = m_Face_handle->vertex(2);
|
||||
const Point_2 & p0 = v0->point();
|
||||
const Point_2 & p1 = v1->point();
|
||||
const Point_2 & p2 = v2->point();
|
||||
FT s0,s1,s2,s;
|
||||
std::vector<Point_2> vec;
|
||||
vec.push_back(p0); vec.push_back(p1); vec.push_back(p2);
|
||||
s = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
vec.push_back(p); vec.push_back(p1); vec.push_back(p2);
|
||||
s0 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
vec.push_back(p0); vec.push_back(p); vec.push_back(p2);
|
||||
s1 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
vec.push_back(p0); vec.push_back(p1); vec.push_back(p);
|
||||
s2 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits());
|
||||
vec.clear();
|
||||
s0 = s0 / s; s1 = s1 / s; s2 = s2 / s;
|
||||
Vector_2 v_0 = field_map[p0];
|
||||
Vector_2 v_1 = field_map[p1];
|
||||
Vector_2 v_2 = field_map[p2];
|
||||
FT x = ((v_0.x()*s0)+(v_1.x()*s1)+(v_2.x()*s2));
|
||||
FT y = ((v_0.y()*s0)+(v_1.y()*s1)+(v_2.y()*s2));
|
||||
FT normal = sqrt((x)*(x) + (y)*(y));
|
||||
if (normal != 0)
|
||||
{
|
||||
x = x / normal;
|
||||
y = y / normal;
|
||||
}
|
||||
Vector_2 v = Vector_2(x, y);
|
||||
return v;
|
||||
}
|
||||
|
||||
template <class StreamLinesTraits_2>
|
||||
typename Triangular_field_2<StreamLinesTraits_2>::FT
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_density_field(const Point_2 & p) const{
|
||||
return p.x();}
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_density_field(const
|
||||
Point_2 &
|
||||
p) const
|
||||
{
|
||||
return p.x();
|
||||
}
|
||||
|
||||
template<class StreamLinesTraits_2>
|
||||
typename Triangular_field_2<StreamLinesTraits_2>::FT
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_integration_step(const Point_2 &) const{
|
||||
return 1.0;}
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_integration_step(const
|
||||
Point_2
|
||||
&) const
|
||||
{
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
template<class StreamLinesTraits_2>
|
||||
typename Triangular_field_2<StreamLinesTraits_2>::FT
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_integration_step() const{
|
||||
return 1.0;}
|
||||
Triangular_field_2<StreamLinesTraits_2>::get_integration_step() const
|
||||
{
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue