From 514669beaf7438b352702e5d19e1381b5eb7e7ff Mon Sep 17 00:00:00 2001
From: Tristan Carrier Baudouin <tristan.carrier@uclouvain.be>
Date: Mon, 8 Oct 2012 15:17:30 +0000
Subject: [PATCH] 3D fields
---
Mesh/directions3D.cpp | 124 +++++++++++++++++++++++++++++++++++++++---
Mesh/directions3D.h | 3 +
2 files changed, 120 insertions(+), 7 deletions(-)
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index ca7bcab8f1..c5d155d59f 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -451,9 +451,7 @@ GRegion* Frame_field::test(){
}
void Frame_field::clear(){
- #if defined(HAVE_ANN)
Nearest_point::clear();
- #endif
field.clear();
random.clear();
#if defined(HAVE_ANN)
@@ -670,12 +668,13 @@ void Nearest_point::init_region(GRegion* gr){
fullMatrix<double> gauss_points;
fullVector<double> gauss_weights;
- gaussIntegration::getTriangle(12,gauss_points,gauss_weights);
+ gaussIntegration::getTriangle(8,gauss_points,gauss_weights);
gauss_num = gauss_points.size1();
faces = gr->faces();
random.clear();
+ vicinity.clear();
vertices.clear();
for(it=faces.begin();it!=faces.end();it++){
@@ -704,6 +703,7 @@ void Nearest_point::init_region(GRegion* gr){
z = T(u,v,z1,z2,z3);
random.push_back(SPoint3(x,y,z));
+ vicinity.push_back(element);
}
vertices.insert(element->getVertex(0));
@@ -716,7 +716,8 @@ void Nearest_point::init_region(GRegion* gr){
x = (*it2)->x();
y = (*it2)->y();
z = (*it2)->z();
- random.push_back(SPoint3(x,y,z));
+ //random.push_back(SPoint3(x,y,z));
+ //vicinity.push_back(NULL);
}
duplicate = annAllocPts(random.size(),3);
@@ -757,9 +758,15 @@ bool Nearest_point::search(double x,double y,double z,SVector3& vec){
annDeallocPt(query);
delete[] indices;
delete[] distances;
-
- e2 = 0.000001;
- point = random[index];
+
+ if(vicinity[index]!=NULL){
+ point = closest(vicinity[index],SPoint3(x,y,z));
+ }
+ else{
+ point = random[index];
+ }
+
+ e2 = 0.000001;
if(fabs(point.x()-x)>e2 || fabs(point.y()-y)>e2 || fabs(point.z()-z)>e2){
vec = SVector3(point.x()-x,point.y()-y,point.z()-z);
@@ -779,6 +786,107 @@ double Nearest_point::T(double u,double v,double val1,double val2,double val3){
return (1.0-u-v)*val1 + u*val2 + v*val3;
}
+//The following method comes from this page : gamedev.net/topic/552906-closest-point-on-triangle
+//It can also be found on this page : geometrictools.com/LibMathematics/Distance/Distance.html
+SPoint3 Nearest_point::closest(MElement* element,SPoint3 point){
+ SVector3 edge0 = SVector3(element->getVertex(1)->x()-element->getVertex(0)->x(),element->getVertex(1)->y()-element->getVertex(0)->y(),
+ element->getVertex(1)->z()-element->getVertex(0)->z());
+ SVector3 edge1 = SVector3(element->getVertex(2)->x()-element->getVertex(0)->x(),element->getVertex(2)->y()-element->getVertex(0)->y(),
+ element->getVertex(2)->z()-element->getVertex(0)->z());
+ SVector3 v0 = SVector3(element->getVertex(0)->x()-point.x(),element->getVertex(0)->y()-point.y(),
+ element->getVertex(0)->z()-point.z());
+
+ double a = dot(edge0,edge0);
+ double b = dot(edge0,edge1);
+ double c = dot(edge1,edge1);
+ double d = dot(edge0,v0);
+ double e = dot(edge1,v0);
+
+ double det = a*c-b*b;
+ double s = b*e-c*d;
+ double t = b*d-a*e;
+
+ if(s+t<det){
+ if(s<0.0){
+ if(t<0.0){
+ if(d<0.0){
+ s = clamp(-d/a,0.0,1.0);
+ t = 0.0;
+ }
+ else{
+ s = 0.0;
+ t = clamp(-e/c,0.0,1.0);
+ }
+ }
+ else{
+ s = 0.0;
+ t = clamp(-e/c,0.0,1.0);
+ }
+ }
+ else if(t<0.0){
+ s = clamp(-d/a,0.0,1.0);
+ t = 0.0;
+ }
+ else{
+ double invDet = 1.0/det;
+ s *= invDet;
+ t *= invDet;
+ }
+ }
+ else{
+ if(s<0.0){
+ double tmp0 = b+d;
+ double tmp1 = c+e;
+ if(tmp1>tmp0){
+ double numer = tmp1-tmp0;
+ double denom = a-2.0*b+c;
+ s = clamp(numer/denom,0.0,1.0);
+ t = 1.0-s;
+ }
+ else{
+ t = clamp(-e/c,0.0,1.0);
+ s = 0.0;
+ }
+ }
+ else if(t<0.0){
+ if(a+d>b+e){
+ double numer = c+e-b-d;
+ double denom = a-2.0*b+c;
+ s = clamp(numer/denom,0.0,1.0);
+ t = 1.0-s;
+ }
+ else{
+ s = clamp(-e/c,0.0,1.0);
+ t = 0.0;
+ }
+ }
+ else{
+ double numer = c+e-b-d;
+ double denom = a-2.0*b+c;
+ s = clamp(numer/denom,0.0,1.0);
+ t = 1.0-s;
+ }
+ }
+
+ return SPoint3(element->getVertex(0)->x()+s*edge0.x()+t*edge1.x(),element->getVertex(0)->y()+s*edge0.y()+t*edge1.y(),
+ element->getVertex(0)->z()+s*edge0.z()+t*edge1.z());
+}
+
+double Nearest_point::clamp(double x,double min,double max){
+ double val;
+
+ val = x;
+
+ if(val<min){
+ val = min;
+ }
+ else if(val>max){
+ val = max;
+ }
+
+ return val;
+}
+
void Nearest_point::print_field(GRegion* gr){
unsigned int i;
int j;
@@ -828,6 +936,7 @@ GRegion* Nearest_point::test(){
void Nearest_point::clear(){
random.clear();
+ vicinity.clear();
#if defined(HAVE_ANN)
delete duplicate;
delete kd_tree;
@@ -848,6 +957,7 @@ std::map<MVertex*,double> Size_field::boundary;
MElementOctree* Size_field::octree;
std::vector<SPoint3> Nearest_point::random;
+std::vector<MElement*> Nearest_point::vicinity;
#if defined(HAVE_ANN)
ANNpointArray Nearest_point::duplicate;
ANNkd_tree* Nearest_point::kd_tree;
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index 2b583e166b..95600ab751 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -80,6 +80,7 @@ class Size_field{
class Nearest_point{
private:
static std::vector<SPoint3> random;
+ static std::vector<MElement*> vicinity;
#if defined(HAVE_ANN)
static ANNpointArray duplicate;
static ANNkd_tree* kd_tree;
@@ -89,6 +90,8 @@ class Nearest_point{
static void init_region(GRegion*);
static bool search(double,double,double,SVector3&);
static double T(double,double,double,double,double);
+ static SPoint3 closest(MElement*,SPoint3);
+ static double clamp(double,double,double);
static void print_field(GRegion*);
static void print_segment(SPoint3,SPoint3,std::ofstream&);
static GRegion* test();
--
GitLab