From c091d76a9615259e015be8ce6852d97333ec1f95 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 16 Feb 2012 21:09:50 +0000
Subject: [PATCH] added MEANPLANE for harmonic map, works well :-)
---
Geo/GFace.h | 8 +--
Geo/GFaceCompound.cpp | 35 ++++++++++-
Geo/GFaceCompound.h | 5 +-
Geo/SVector3.h | 58 ++++++++++++++++++
Mesh/BackgroundMesh.cpp | 6 +-
Mesh/CenterlineField.cpp | 54 +++-------------
Numeric/Numeric.cpp | 129 +++++++++++++++++++++++++++++++++++++++
Numeric/Numeric.h | 14 +++++
8 files changed, 245 insertions(+), 64 deletions(-)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index fabe4085a1..0d430ff287 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -16,6 +16,7 @@
#include "SPoint2.h"
#include "SVector3.h"
#include "Pair.h"
+#include "Numeric.h"
class MElement;
class MTriangle;
@@ -24,13 +25,6 @@ class MPolygon;
class ExtrudeParams;
class GFaceCompound;
-struct mean_plane
-{
- double plan[3][3];
- double a, b, c, d;
- double x, y, z;
-};
-
struct surface_params
{
double radius, radius2, height, cx, cy, cz;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e88a83ff0c..a0806068da 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -41,6 +41,7 @@
#include "GRbf.h"
#include "Curvature.h"
#include "MPoint.h"
+#include "Numeric.h"
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
@@ -210,7 +211,7 @@ static void computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
vcg/=nbFinal;
}
else{
- Msg::Debug("----> No Kernel for polygon: place point at CG polygon");
+ //Msg::Debug("----> No Kernel for polygon: place point at CG polygon");
//place at CG polygon
for(std::vector<MVertex*>::iterator it = cavV.begin() ; it != cavV.end() ; ++it){
SPoint3 vsp = coordinates[*it];
@@ -644,10 +645,13 @@ bool GFaceCompound::parametrize() const
// Laplace parametrization
if (_mapping == HARMONIC){
- Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag());
+ Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
+ printStuff(111);
+ checkOrientation(0, true);
+ printStuff(222);
}
// Multiscale Laplace parametrization
else if (_mapping == MULTISCALE){
@@ -660,7 +664,7 @@ bool GFaceCompound::parametrize() const
}
// Conformal map parametrization
else if (_mapping == CONFORMAL){
- Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
+ Msg::Info("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
std::vector<MVertex *> vert;
bool oriented, hasOverlap;
@@ -971,6 +975,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
getBoundingEdges();
_type = UNITCIRCLE;
+ //_type = MEANPLANE;
nbSplit = 0;
fillTris.clear();
@@ -1096,6 +1101,30 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, sin(theta));
}
}
+ else if (_type == MEANPLANE){
+
+ std::vector<SPoint3> points, pointsProj, pointsUV;
+ SPoint3 ptCG(0.0, 0.0, 0.0);
+ for(unsigned int i = 0; i < _ordered.size(); i++){
+ MVertex *v = _ordered[i];
+ points.push_back(SPoint3(v->x(), v->y(), v->z()));
+ ptCG[0] += v->x();
+ ptCG[1] += v->y();
+ ptCG[2] += v->z();
+ }
+ ptCG /= points.size();
+ mean_plane meanPlane;
+ computeMeanPlaneSimple(points, meanPlane);
+ projectPointsToPlane(points, pointsProj, meanPlane);
+ transformPointsIntoOrthoBasis(pointsProj, pointsUV, ptCG, meanPlane);
+
+ for(unsigned int i = 0; i < pointsUV.size(); i++){
+ MVertex *v = _ordered[i];
+ if(step == ITERU) myAssembler.fixVertex(v, 0, 1, pointsUV[i][0]);
+ else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, pointsUV[i][1]);
+ }
+
+ }
else{
Msg::Error("Unknown type of parametrization");
return;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 41a239b8f3..6f92c9691d 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -57,7 +57,7 @@ class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4, RBF=5} typeOfMapping;
- typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
+ typedef enum {UNITCIRCLE, MEANPLANE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
protected:
mutable std::set<MVertex *> ov;
@@ -134,6 +134,7 @@ class GFaceCompound : public GFace {
mutable int nbSplit;
int getNbSplit() const { return nbSplit; }
int allowPartition() const{ return _allowPartition; }
+ void setType(typeOfIsomorphism type){ _type=type;}
private:
typeOfIsomorphism _type;
mutable typeOfMapping _mapping;
@@ -147,7 +148,7 @@ class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4, RBF=5} typeOfMapping;
- typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
+ typedef enum {UNITCIRCLE, MEANPLANE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC, int allowPartition=1)
: GFace(m, tag)
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 54ca826d80..db3c6e86d2 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -9,6 +9,7 @@
#include "SPoint3.h"
#include <string>
#include <stdio.h>
+#include "GmshMessage.h"
// concrete class for vector of size 3
class SVector3 {
@@ -71,6 +72,7 @@ class SVector3 {
// Needed to allow the initialization of a SPoint3 from a SPoint3, a distance and a direction
SPoint3 point() const{return P;}
+
};
inline double dot(const SVector3 &a, const SVector3 &b)
@@ -102,4 +104,60 @@ inline SVector3 operator+(const SVector3 &a,const SVector3 &b)
inline SVector3 operator-(const SVector3 &a,const SVector3 &b)
{ return SVector3(a[0] - b[0], a[1] - b[1], a[2] - b[2]); }
+inline void buildOrthoBasis_naive(SVector3 &dir, SVector3 &dir1, SVector3 &dir2)
+{
+ dir.normalize();
+ if (dir[1]!=0.0 && dir[2]!=0.0){
+ dir1 = SVector3(1.0, 0.0, -dir[0]/dir[2]);
+ dir2 = SVector3 (dir[0]/dir[2], -(dir[0]*dir[0]+dir[2]*dir[2])/(dir[1]*dir[2]), 1.0);
+ }
+ else if (dir[0]!=0.0 && dir[2]!=0.0){
+ dir1 = SVector3(-dir[1]/dir[0], 1.0, 0.0);
+ dir2 = SVector3(1.0, dir[1]/dir[0], -(dir[1]*dir[1]+dir[0]*dir[0])/(dir[0]*dir[2]));
+ }
+ else if (dir[0]!=0.0 && dir[1]!=0.0){
+ dir1 = SVector3(0.0, -dir[2]/dir[1], 1.0);
+ dir2 = SVector3(-(dir[1]*dir[1]+dir[2]*dir[2])/(dir[0]*dir[1]), 1.0, dir[2]/dir[1]);
+ }
+ else if (dir[0]==0.0 && dir[1]==0.0){
+ dir1 = SVector3(0.0, 1.0, 0.0);
+ dir2 = SVector3(1.0, 0.0, 0.0);
+ }
+ else if (dir[1]==0.0 && dir[2]==0.0){
+ dir1 = SVector3(0.0, 1.0, 0.0);
+ dir2 = SVector3(0.0, 0.0, 1.0);
+ }
+ else if (dir[0]==0.0 && dir[2]==0.0){
+ dir1 = SVector3(1.0, 0.0, 0.0);
+ dir2 = SVector3(0.0, 0.0, 1.0);
+ }
+ else{
+ Msg::Error("Problem with computing orthoBasis");
+ Msg::Exit(1);
+ }
+ // printf("XYZ =%g %g %g r=%g dir0 = %g %g %g dir1 = %g %g %g dir2 =%g %g %g\n",
+ // x,y,z,d, dir[0], dir[1], dir[2], dir1[0], dir1[1], dir1[2], dir2[0], dir2[1], dir2[2] );
+ // printf("0x1 =%g 1x2=%g 2x1=%g \n", dot(dir, dir1), dot(dir1,dir2), dot(dir2,dir));
+ dir1.normalize();
+ dir2.normalize();
+}
+
+inline void buildOrthoBasis(SVector3 &normal, SVector3 &tangent, SVector3 &binormal)
+{
+ //pick any Unit-Vector that's not parallel to normal:
+ normal.normalize();
+ if (normal[0] > normal[1] )
+ tangent = SVector3(0.0, 1.0, 0.0);
+ else
+ tangent = SVector3(1.0, 0.0, 0.0);
+
+ //build a binormal from tangent and normal:
+ binormal = crossprod(tangent, normal);
+ binormal.normalize();
+
+ //and correct the tangent from the binormal and the normal.
+ tangent = crossprod(normal, binormal);
+ tangent.normalize();
+}
+
#endif
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index fc8fffdf66..ce0cd09f6f 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -246,17 +246,13 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
}
break;
}
-
double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : MAX_LC;
return lc;
}
static SMetric3 LC_MVertex_CURV_ANISO(GEntity *ge, double U, double V)
{
- //std::cout << "I'm in LC_MVertex_CURV_ANISO" << std::endl;
- //std::cout << "The dimension of the entity is: "<< ge->dim() << std::endl;
switch(ge->dim()){
- //std::cout << "The dimension of the entity is: "<< ge->dim() << std::endl;
case 0: return metric_based_on_surface_curvature((const GVertex *)ge);
case 1: return metric_based_on_surface_curvature((const GEdge *)ge, U);
case 2: return metric_based_on_surface_curvature((const GFace *)ge, U, V);
@@ -397,7 +393,7 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
SMetric3 LC(1./(lc*lc));
SMetric3 m = intersection(intersection (l4, LC),l3);
- // printf("%g %g %g %g %g %g\n",m(0,0),m(1,1),m(2,2),m(0,1),m(0,2),m(1,2));
+ //printf("%g %g %g %g %g %g\n",m(0,0),m(1,1),m(2,2),m(0,1),m(0,2),m(1,2));
return m;
// return lc * CTX::instance()->mesh.lcFactor;
}
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index edc5a94871..070c25392c 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -698,20 +698,14 @@ void Centerline::createSplitCompounds(){
int num_gfc = NF+i+1;
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
+ //GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
+ GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
gfc->meshAttributes.recombine = recombine;
gfc->addPhysicalEntity(i+1);
current->add(gfc);
- //gfc->parametrize();
- // for(unsigned int j = 0; j < pf->triangles.size(); ++j){
- // MTriangle *t = pf->triangles[j];
- // for(int k = 0; k < 3; k++){
- // MVertex *v = t->getVertex(k);
- // v->setEntity(pf);
- // }
- //}
+
}
}
@@ -867,7 +861,7 @@ void Centerline::cutMesh(){
for(unsigned int i = 0; i < edges.size(); i++){
std::vector<MLine*> lines = edges[i].lines;
double L = edges[i].length;
- double R = edges[i].minRad;
+ double R = 0.5*(edges[i].minRad+edges[i].maxRad);
double AR = L/R;
printf("*** branch =%d \n", i, AR);
if( AR > 4.5){
@@ -1044,48 +1038,14 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
double lc = 2*M_PI*d/nbPoints;
SVector3 p0(nodes[index2[0]][0], nodes[index2[0]][1], nodes[index2[0]][2]);
SVector3 p1(nodes[index2[1]][0], nodes[index2[1]][1], nodes[index2[1]][2]);
- SVector3 dir = p1-p0;
- dir.normalize();
+ SVector3 dir0 = p1-p0;
SVector3 dir1, dir2;
- if (dir[1]!=0.0 && dir[2]!=0.0){
- dir1 = SVector3(1.0, 0.0, -dir[0]/dir[2]);
- dir2 = SVector3 (dir[0]/dir[2], -(dir[0]*dir[0]+dir[2]*dir[2])/(dir[1]*dir[2]), 1.0);
- }
- else if (dir[0]!=0.0 && dir[2]!=0.0){
- dir1 = SVector3(-dir[1]/dir[0], 1.0, 0.0);
- dir2 = SVector3(1.0, dir[1]/dir[0], -(dir[1]*dir[1]+dir[0]*dir[0])/(dir[0]*dir[2]));
- }
- else if (dir[0]!=0.0 && dir[1]!=0.0){
- dir1 = SVector3(0.0, -dir[2]/dir[1], 1.0);
- dir2 = SVector3(-(dir[1]*dir[1]+dir[2]*dir[2])/(dir[0]*dir[1]), 1.0, dir[2]/dir[1]);
- }
- else if (dir[0]==0.0 && dir[1]==0.0){
- dir1 = SVector3(0.0, 1.0, 0.0);
- dir2 = SVector3(1.0, 0.0, 0.0);
- }
- else if (dir[1]==0.0 && dir[2]==0.0){
- dir1 = SVector3(0.0, 1.0, 0.0);
- dir2 = SVector3(0.0, 0.0, 1.0);
- }
- else if (dir[0]==0.0 && dir[2]==0.0){
- dir1 = SVector3(1.0, 0.0, 0.0);
- dir2 = SVector3(0.0, 0.0, 1.0);
- }
- else {printf("ARGHH EMI DO STHG \n"); exit(1);}
- // printf("XYZ =%g %g %g r=%g dir0 = %g %g %g dir1 = %g %g %g dir2 =%g %g %g\n",
- // x,y,z,d, dir[0], dir[1], dir[2], dir1[0], dir1[1], dir1[2], dir2[0], dir2[1], dir2[2] );
- // printf("0x1 =%g 1x2=%g 2x1=%g \n", dot(dir, dir1), dot(dir1,dir2), dot(dir2,dir));
- dir1.normalize();
- dir2.normalize();
-
- //dir = SVector3(1.0, 0.0, 0.0);
- //dir1 = SVector3(0.0, 1.0, 0.0);
- //dir2 = SVector3(0.0, 0.0, 1.0);
+ buildOrthoBasis(dir0,dir1,dir2);
double lcA = 4.*lc;
double lam1 = 1./(lcA*lcA);
double lam2 = 1./(lc*lc);
- metr = SMetric3(lam1,lam2,lam2, dir, dir1, dir2);
+ metr = SMetric3(lam1,lam2,lam2, dir0, dir1, dir2);
delete[]index2;
delete[]dist2;
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 18192ca5a5..f1eabcea63 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1226,3 +1226,132 @@ int intersection_segments(SPoint3 &p1, SPoint3 &p2,
x[1] >= 0.0 && x[1] <= 1.);
}
}
+
+void computeMeanPlaneSimple(const std::vector<SPoint3> &points, mean_plane &meanPlane)
+{
+
+ double xm = 0., ym = 0., zm = 0.;
+ int ndata = points.size();
+ int na = 3;
+ for(int i = 0; i < ndata; i++) {
+ xm += points[i].x();
+ ym += points[i].y();
+ zm += points[i].z();
+ }
+ xm /= (double)ndata;
+ ym /= (double)ndata;
+ zm /= (double)ndata;
+
+ fullMatrix<double> U(ndata, na), V(na, na);
+ fullVector<double> sigma(na);
+ for(int i = 0; i < ndata; i++) {
+ U(i, 0) = points[i].x() - xm;
+ U(i, 1) = points[i].y() - ym;
+ U(i, 2) = points[i].z() - zm;
+ }
+ U.svd(V, sigma);
+ double res[4], svd[3];
+ svd[0] = sigma(0);
+ svd[1] = sigma(1);
+ svd[2] = sigma(2);
+ int min;
+ if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
+ min = 0;
+ else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
+ min = 1;
+ else
+ min = 2;
+ res[0] = V(0, min);
+ res[1] = V(1, min);
+ res[2] = V(2, min);
+ norme(res);
+
+ double ex[3], t1[3], t2[3];
+
+ ex[0] = ex[1] = ex[2] = 0.0;
+ if(res[0] == 0.)
+ ex[0] = 1.0;
+ else if(res[1] == 0.)
+ ex[1] = 1.0;
+ else
+ ex[2] = 1.0;
+
+ prodve(res, ex, t1);
+ norme(t1);
+ prodve(t1, res, t2);
+ norme(t2);
+
+ res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
+
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[0][i] = t1[i];
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[1][i] = t2[i];
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[2][i] = res[i];
+
+ meanPlane.a = res[0];
+ meanPlane.b = res[1];
+ meanPlane.c = res[2];
+ meanPlane.d = -res[3];//BUG HERE
+
+ meanPlane.x = meanPlane.y = meanPlane.z = 0.;
+ if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
+ fabs(meanPlane.a) >= fabs(meanPlane.c) ){
+ meanPlane.x = meanPlane.d / meanPlane.a;
+ }
+ else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
+ fabs(meanPlane.b) >= fabs(meanPlane.c)){
+ meanPlane.y = meanPlane.d / meanPlane.b;
+ }
+ else{
+ meanPlane.z = meanPlane.d / meanPlane.c;
+ }
+
+ }
+
+void projectPointToPlane(const SPoint3 &pt, SPoint3 &ptProj, const mean_plane &meanPlane)
+{
+
+ double u = pt.x();
+ double v = pt.y();
+ double w = pt.z();
+ double a = meanPlane.a;
+ double b = meanPlane.b;
+ double c = meanPlane.c;
+ double d = meanPlane.d;
+ double t0 = -(a*u+b*v+c*w+d)/(a*a+b*b+c*c);
+
+ ptProj[0] = u + a*t0;
+ ptProj[1] = v + b*t0;
+ ptProj[2] = w + c*t0;
+
+}
+
+void projectPointsToPlane(const std::vector<SPoint3> &pts, std::vector<SPoint3> &ptsProj, const mean_plane &meanPlane)
+{
+
+ ptsProj.resize(pts.size());
+ for (int i= 0; i< pts.size(); i++){
+ projectPointToPlane(pts[i],ptsProj[i], meanPlane);
+ }
+
+}
+
+void transformPointsIntoOrthoBasis(const std::vector<SPoint3> &ptsProj, std::vector<SPoint3> &pointsUV,
+ const SPoint3 &ptCG, const mean_plane &meanPlane)
+{
+
+ pointsUV.resize(ptsProj.size());
+ SVector3 normal(meanPlane.a, meanPlane.b, meanPlane.c);
+ SVector3 tangent, binormal;
+ buildOrthoBasis(normal, tangent, binormal);
+
+ for (int i= 0; i< ptsProj.size(); i++){
+ SVector3 pp(ptsProj[i][0]-ptCG[0],ptsProj[i][1]-ptCG[1],ptsProj[i][2]-ptCG[2]) ;
+ pointsUV[i][0] = dot(pp, tangent);
+ pointsUV[i][1] = dot(pp, binormal);
+ pointsUV[i][2] = dot(pp, normal);
+ }
+
+}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index b066b30a05..10ea14891b 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -15,6 +15,13 @@
#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
#define sign(x) (((x)>=0)?1:-1)
+struct mean_plane
+{
+ double plan[3][3];
+ double a, b, c, d;
+ double x, y, z;
+};
+
double myatan2(double a, double b);
double myasin(double a);
double myacos(double a);
@@ -128,4 +135,11 @@ int intersection_segments (SPoint3 &p1, SPoint3 &p2,
SPoint3 &q1, SPoint3 &q2,
double x[2]);
+//tools for projection onto plane
+void computeMeanPlaneSimple(const std::vector<SPoint3> &points, mean_plane &meanPlane);
+void projectPointToPlane(const SPoint3 &pt, SPoint3 &ptProj, const mean_plane &meanPlane);
+void projectPointsToPlane(const std::vector<SPoint3> &pts, std::vector<SPoint3> &ptsProj, const mean_plane &meanPlane);
+void transformPointsIntoOrthoBasis(const std::vector<SPoint3> &ptsProj, std::vector<SPoint3> &pointsUV,
+ const SPoint3 &ptCG, const mean_plane &meanPlane);
+
#endif
--
GitLab