Commit c8b7f315 authored by Sauvage's avatar Sauvage

Normal cycles : calcul du tenseur déplacé dans le collecteur

parent 2aaa9d77
......@@ -27,6 +27,8 @@
#include "Geometry/basic.h"
#include "Utils/convertType.h"
#include "OpenNL/linear_solver.h"
#include "OpenNL/sparse_matrix.h"
#include "OpenNL/full_vector.h"
......
......@@ -310,43 +310,23 @@ void computeCurvatureVertex_NormalCycles(
VertexAttribute<typename PFP::VEC3>& Kmin,
VertexAttribute<typename PFP::VEC3>& Knormal, unsigned int thread)
{
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ;
typedef typename PFP::VEC3 VEC3 ;
typedef Geom::Matrix<3,3,REAL> MATRIX;
typedef Eigen::Matrix<REAL,3,1> E_VEC3;
typedef Eigen::Matrix<REAL,3,3> E_MATRIX;
// collect the normal cycle tensor
Algo::Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
neigh.collectAll(dart) ;
neigh.computeArea() ;
VEC3 center = position[dart] ;
typename PFP::MATRIX33 tensor(0) ;
// collect edges inside the neighborhood
const std::vector<Dart>& vd1 = neigh.getInsideEdges() ;
for (std::vector<Dart>::const_iterator it = vd1.begin(); it != vd1.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(map, *it, position) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
}
// collect edges crossing the neighborhood's border
const std::vector<Dart>& vd2 = neigh.getBorder() ;
for (std::vector<Dart>::const_iterator it = vd2.begin(); it != vd2.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(map, *it, position) ;
REAL alpha ;
Algo::Geometry::intersectionSphereEdge<PFP>(map, center, radius, *it, position, alpha) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) * alpha ;
}
tensor /= neigh.getArea() ;
MATRIX tensor(0) ;
neigh.computeNormalCyclesTensor(edgeangle,tensor);
// solve eigen problem
Eigen::Matrix3f m3 ;
m3 << tensor(0,0) , tensor(0,1) , tensor(0,2) , tensor(1,0) , tensor(1,1) , tensor(1,2) , tensor(2,0) , tensor(2,1) , tensor(2,2) ;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> solver (m3);
Eigen::Vector3f ev = solver.eigenvalues();
Eigen::Matrix3f evec = solver.eigenvectors();
Eigen::SelfAdjointEigenSolver<E_MATRIX> solver (Utils::convertRef<E_MATRIX>(tensor));
const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());
// sort eigen components : ev[s[0]] has minimal absolute value ; kmin = ev[s[1]] <= ev[s[2]] = kmax
int s[3] = {0, 1, 2} ;
......
......@@ -112,6 +112,13 @@ public:
template <typename PPFP>
friend std::ostream& operator<<(std::ostream &out, const Collector<PPFP>& c);
REAL computeArea () {
assert(!"Warning: Collector<PFP>::computeArea() should be overloaded in non-virtual derived classes");
return 0.0;
}
void computeNormalCyclesTensor (const EdgeAttribute<REAL>&, typename PFP::MATRIX33&) {assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes"); }
};
/*********************************************************
......@@ -149,26 +156,27 @@ template <typename PFP>
class Collector_WithinSphere : public Collector<PFP>
{
protected:
const VertexAttribute<typename PFP::VEC3>& position;
typename PFP::REAL radius;
typename PFP::REAL area;
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ;
const VertexAttribute<VEC3>& position;
REAL radius;
public:
Collector_WithinSphere(typename PFP::MAP& m, const VertexAttribute<typename PFP::VEC3>& p, typename PFP::REAL r = 0, unsigned int thread=0) :
Collector_WithinSphere(typename PFP::MAP& m, const VertexAttribute<VEC3>& p, REAL r = 0, unsigned int thread=0) :
Collector<PFP>(m, thread),
position(p),
radius(r),
area(0)
radius(r)
{}
inline void setRadius(typename PFP::REAL r) { radius = r; }
inline typename PFP::REAL getRadius() const { return radius; }
inline const VertexAttribute<typename PFP::VEC3>& getPosition() const { return position; }
inline void setRadius(REAL r) { radius = r; }
inline REAL getRadius() const { return radius; }
inline const VertexAttribute<VEC3>& getPosition() const { return position; }
void collectAll(Dart d);
void collectBorder(Dart d);
void computeArea();
inline typename PFP::REAL getArea() const { return area; }
REAL computeArea();
void computeNormalCyclesTensor (const EdgeAttribute<REAL>&, typename PFP::MATRIX33&);
};
/*********************************************************
......
......@@ -256,12 +256,14 @@ void Collector_WithinSphere<PFP>::collectBorder(Dart d)
this->insideVertices.clear();
}
template <typename PFP>
void Collector_WithinSphere<PFP>::computeArea()
typename PFP::REAL Collector_WithinSphere<PFP>::computeArea()
{
assert(this->isInsideCollected || !"computeArea: inside cells have not been collected.") ;
area = 0;
typename PFP::VEC3 centerPosition = this->position[this->centerDart];
VEC3 centerPosition = position[this->centerDart];
REAL area = 0;
for (std::vector<Dart>::const_iterator it = this->insideFaces.begin(); it != this->insideFaces.end(); ++it)
area += Algo::Geometry::triangleArea<PFP>(this->map, *it, this->position);
......@@ -285,8 +287,35 @@ void Collector_WithinSphere<PFP>::computeArea()
area += alpha * beta * Algo::Geometry::triangleArea<PFP>(this->map, *it, this->position);
}
}
return area;
}
template <typename PFP>
void Collector_WithinSphere<PFP>::computeNormalCyclesTensor (const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33& tensor){
assert(this->isInsideCollected || !"computeNormalCyclesTensor: inside cells have not been collected.") ;
VEC3 centerPosition = position[this->centerDart];
tensor.zero() ;
// collect edges inside the neighborhood
for (std::vector<Dart>::const_iterator it = this->insideEdges.begin(); it != this->insideEdges.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, position) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
}
// collect edges crossing the neighborhood's border
for (std::vector<Dart>::const_iterator it = this->border.begin(); it != this->border.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, position) ;
REAL alpha ;
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, radius, *it, position, alpha) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) * alpha ;
}
tensor /= computeArea() ;
}
/*********************************************************
* Collector Normal Angle (Vertices)
*********************************************************/
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment