From c4d326325697b2fd9a4269cd120f4014a8765e62 Mon Sep 17 00:00:00 2001
From: Francois Henrotte <francois.henrotte@ulg.ac.be>
Date: Mon, 18 Feb 2013 14:14:20 +0000
Subject: [PATCH] cross fields
---
Mesh/cross3D.h | 11 +++
Mesh/directions3D.cpp | 215 ++++++++++++++++--------------------------
Mesh/directions3D.h | 43 +++------
Mesh/simple3D.cpp | 1 +
Mesh/yamakawa.h | 6 ++
5 files changed, 112 insertions(+), 164 deletions(-)
diff --git a/Mesh/cross3D.h b/Mesh/cross3D.h
index d1ebcc59af..8c4a3553ae 100644
--- a/Mesh/cross3D.h
+++ b/Mesh/cross3D.h
@@ -1,3 +1,13 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributor(s):
+// François Henrotte
+
+#ifndef _CROSS_3D_H_
+#define _CROSS_3D_H_
#include "stat_cwmtx.h"
#include "stat_quatern.h"
@@ -196,3 +206,4 @@ Matrix convert(const cross3D x){
+#endif
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index d15e6f2b86..78c7cdd7cf 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -6,12 +6,12 @@
// Contributor(s):
// Tristan Carrier
-#include "directions3D.h"
+#include <fstream>
#include "GModel.h"
#include "BackgroundMesh.h"
#include "meshGFaceDelaunayInsertion.h"
-#include <fstream>
#include "MTetrahedron.h"
+#include "directions3D.h"
#if defined(HAVE_PETSC)
#include "dofManager.h"
@@ -20,94 +20,6 @@
#include "linearSystemFull.h"
#endif
-/****************class Matrix****************/
-
-Matrix::Matrix(){
- m11 = 1.0;
- m21 = 0.0;
- m31 = 0.0;
- m12 = 0.0;
- m22 = 1.0;
- m32 = 0.0;
- m13 = 0.0;
- m23 = 0.0;
- m33 = 1.0;
-}
-
-Matrix::~Matrix(){}
-
-void Matrix::set_m11(double new_m11){
- m11 = new_m11;
-}
-
-void Matrix::set_m21(double new_m21){
- m21 = new_m21;
-}
-
-void Matrix::set_m31(double new_m31){
- m31 = new_m31;
-}
-
-void Matrix::set_m12(double new_m12){
- m12 = new_m12;
-}
-
-void Matrix::set_m22(double new_m22){
- m22 = new_m22;
-}
-
-void Matrix::set_m32(double new_m32){
- m32 = new_m32;
-}
-
-void Matrix::set_m13(double new_m13){
- m13 = new_m13;
-}
-
-void Matrix::set_m23(double new_m23){
- m23 = new_m23;
-}
-
-void Matrix::set_m33(double new_m33){
- m33 = new_m33;
-}
-
-double Matrix::get_m11(){
- return m11;
-}
-
-double Matrix::get_m21(){
- return m21;
-}
-
-double Matrix::get_m31(){
- return m31;
-}
-
-double Matrix::get_m12(){
- return m12;
-}
-
-double Matrix::get_m22(){
- return m22;
-}
-
-double Matrix::get_m32(){
- return m32;
-}
-
-double Matrix::get_m13(){
- return m13;
-}
-
-double Matrix::get_m23(){
- return m23;
-}
-
-double Matrix::get_m33(){
- return m33;
-}
-
/****************class Frame_field****************/
Frame_field::Frame_field(){}
@@ -285,7 +197,6 @@ Matrix Frame_field::search(double x,double y,double z){
delete[] indices;
delete[] distances;
#endif
-
return field[index].second;
}
@@ -359,11 +270,11 @@ double Frame_field::get_ratio(GFace* gf,SPoint2 point){
return val;
}
-void Frame_field::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
+void Frame_field::print_segment(SPoint3 p1,SPoint3 p2,double val1, double val2, std::ofstream& file){
file << "SL ("
<< p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
<< p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
- << "{" << p1.z() << "," << p1.z() << "};\n";
+ << "{" << val1 << "," << val2 << "};\n";
}
void Frame_field::print_field1(){
@@ -403,12 +314,13 @@ void Frame_field::print_field1(){
p6 = SPoint3(vertex->x() - k*m.get_m13(),
vertex->y() - k*m.get_m23(),
vertex->z() - k*m.get_m33());
- print_segment(p,p1,file);
- print_segment(p,p2,file);
- print_segment(p,p3,file);
- print_segment(p,p4,file);
- print_segment(p,p5,file);
- print_segment(p,p6,file);
+ double val1=10, val2=20;
+ print_segment(p,p1,val1,val2,file);
+ print_segment(p,p2,val1,val2,file);
+ print_segment(p,p3,val1,val2,file);
+ print_segment(p,p4,val1,val2,file);
+ print_segment(p,p5,val1,val2,file);
+ print_segment(p,p6,val1,val2,file);
}
file << "};\n";
}
@@ -434,8 +346,6 @@ void Frame_field::print_field2(GRegion* gr){
vertex = element->getVertex(j);
if(vertex->onWhat()->dim()>2){
m = search(vertex->x(),vertex->y(),vertex->z());
- //m = combine(vertex->x(),vertex->y(),vertex->z());
-
p = SPoint3(vertex->x(),vertex->y(),vertex->z());
p1 = SPoint3(vertex->x() + k*m.get_m11(),
vertex->y() + k*m.get_m21(),
@@ -455,12 +365,13 @@ void Frame_field::print_field2(GRegion* gr){
p6 = SPoint3(vertex->x() - k*m.get_m13(),
vertex->y() - k*m.get_m23(),
vertex->z() - k*m.get_m33());
- print_segment(p,p1,file);
- print_segment(p,p2,file);
- print_segment(p,p3,file);
- print_segment(p,p4,file);
- print_segment(p,p5,file);
- print_segment(p,p6,file);
+ double val1=10, val2=20;
+ print_segment(p,p1,val1,val2,file);
+ print_segment(p,p2,val1,val2,file);
+ print_segment(p,p3,val1,val2,file);
+ print_segment(p,p4,val1,val2,file);
+ print_segment(p,p5,val1,val2,file);
+ print_segment(p,p6,val1,val2,file);
}
}
}
@@ -470,9 +381,7 @@ void Frame_field::print_field2(GRegion* gr){
GRegion* Frame_field::test(){
GRegion* gr;
GModel* model = GModel::current();
-
gr = *(model->firstRegion());
-
return gr;
}
@@ -480,15 +389,13 @@ void Frame_field::clear(){
Nearest_point::clear();
temp.clear();
field.clear();
- #if defined(HAVE_ANN)
+#if defined(HAVE_ANN)
delete duplicate;
delete kd_tree;
annClose();
- #endif
+#endif
}
-
-#include "yamakawa.h"
#include "cross3D.h"
CWTSSpaceVector<> convert(const SVector3 v){
@@ -504,27 +411,21 @@ ostream& operator<< (ostream &os, SVector3 &v) {
return os;
}
-void Frame_field::init(GRegion* gr){
- //
- MVertex* pVertex;
- Matrix m;
+void Frame_field::init(GRegion* gr){
CWTSSpaceVector<> Ex(1,0,0), Ey(0,1,0), Ez(0,0,1);
- cross3D x(Ez, Ex);
CWTSSpaceVector<> v(1,2,0), w(0,1,1.3);
+ cross3D x(Ez, Ex);
cross3D y = cross3D(v, w);
- Recombinator crossData;
+ //Recombinator crossData;
crossData.build_vertex_to_vertices(gr);
for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
= crossData.vertex_to_vertices.begin();
iter != crossData.vertex_to_vertices.end(); ++iter){
- pVertex = iter->first;
- m = search(pVertex->x(), pVertex->y(), pVertex->z());
+ MVertex* pVertex = iter->first;
+ Matrix m = search(pVertex->x(), pVertex->y(), pVertex->z());
// if(pVertex->onWhat()->dim() <= 2)
- // m = search(pVertex->x(), pVertex->y(), pVertex->z());
- // else
- // m = convert(y);
crossField[pVertex->getNum()] = m;
}
}
@@ -593,8 +494,8 @@ double Frame_field::smooth(GRegion* gr){
Matrix m, m0;
double enew, eold;
- Recombinator crossData;
- crossData.build_vertex_to_vertices(gr);
+ // Recombinator crossData;
+ // crossData.build_vertex_to_vertices(gr);
//std::cout << "NbVert= " << crossData.vertex_to_vertices.size() << std::endl ;
double energy = 0;
@@ -632,8 +533,8 @@ void Frame_field::save(GRegion* gr, const std::string& filename){
MVertex* pVertex;
std::map<MVertex*,Matrix>::iterator it;
Matrix m;
- Recombinator crossData;
- crossData.build_vertex_to_vertices(gr);
+ // Recombinator crossData;
+ // crossData.build_vertex_to_vertices(gr);
double k = 0.04;
std::ofstream file(filename.c_str());
@@ -652,32 +553,71 @@ void Frame_field::save(GRegion* gr, const std::string& filename){
p1 = SPoint3(pVertex->x() + k*m.get_m11(),
pVertex->y() + k*m.get_m21(),
pVertex->z() + k*m.get_m31());
- print_segment(p,p1,file);
+ print_segment(p,p1,pVertex->z(),pVertex->z(),file);
p1 = SPoint3(pVertex->x() - k*m.get_m11(),
pVertex->y() - k*m.get_m21(),
pVertex->z() - k*m.get_m31());
- print_segment(p,p1,file);
+ print_segment(p,p1,pVertex->z(),pVertex->z(),file);
p1 = SPoint3(pVertex->x() + k*m.get_m12(),
pVertex->y() + k*m.get_m22(),
pVertex->z() + k*m.get_m32());
- print_segment(p,p1,file);
+ print_segment(p,p1,pVertex->z(),pVertex->z(),file);
p1 = SPoint3(pVertex->x() - k*m.get_m12(),
pVertex->y() - k*m.get_m22(),
pVertex->z() - k*m.get_m32());
- print_segment(p,p1,file);
+ print_segment(p,p1,pVertex->z(),pVertex->z(),file);
p1 = SPoint3(pVertex->x() + k*m.get_m13(),
pVertex->y() + k*m.get_m23(),
pVertex->z() + k*m.get_m33());
- print_segment(p,p1,file);
+ print_segment(p,p1,pVertex->z(),pVertex->z(),file);
p1 = SPoint3(pVertex->x() - k*m.get_m13(),
pVertex->y() - k*m.get_m23(),
pVertex->z() - k*m.get_m33());
- print_segment(p,p1,file);
+ print_segment(p,p1,pVertex->z(),pVertex->z(),file);
}
file << "};\n";
file.close();
}
+void Frame_field::fillTreeVolume(GRegion* gr){
+#if defined(HAVE_ANN)
+ int n = crossData.vertex_to_vertices.size();
+ std::cout << "Fillin ANN tree with " << n << " vertices" << std::endl;
+ annTreeData = annAllocPts(n,3);
+ int index=0;
+ for(std::map<MVertex*, std::set<MVertex*> >::iterator iter = crossData.vertex_to_vertices.begin();
+ iter != crossData.vertex_to_vertices.end(); ++iter){
+ MVertex* pVertex = iter->first;
+ annTreeData[index][0] = pVertex->x();
+ annTreeData[index][1] = pVertex->y();
+ annTreeData[index][2] = pVertex->z();
+ vertIndices[index++] = pVertex->getNum();
+ }
+ annTree = new ANNkd_tree(annTreeData,n,3);
+#endif
+}
+
+Matrix Frame_field::findNearestCross(double x,double y,double z){
+#if defined(HAVE_ANN)
+ ANNpoint query = annAllocPt(3);
+ ANNidxArray indices = new ANNidx[1];
+ ANNdistArray distances = new ANNdist[1];
+ query[0] = x;
+ query[1] = y;
+ query[2] = z;
+ double e = 0.0;
+ annTree->annkSearch(query,1,indices,distances,e);
+ annDeallocPt(query);
+ delete[] indices;
+ delete[] distances;
+ std::map<int, Matrix>::const_iterator it = crossField.find(vertIndices[indices[0]]);
+ //annTreeData[indices[0]] gives the coordinates
+ return it->second;
+#else
+ return Matrix();
+#endif
+}
+
void Frame_field::save_dist(const std::string& filename){
std::ofstream file(filename.c_str());
file << "View \"Distance\" {\n";
@@ -1291,7 +1231,8 @@ void Nearest_point::print_field(GRegion* gr){
z = vertex->z();
val = search(x,y,z,vec);
if(val){
- print_segment(SPoint3(x+k*vec.x(),y+k*vec.y(),z+k*vec.z()),SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
+ print_segment(SPoint3(x+k*vec.x(),y+k*vec.y(),z+k*vec.z()),
+ SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
}
}
}
@@ -1331,10 +1272,14 @@ std::map<MVertex*,Matrix> Frame_field::temp;
std::vector<std::pair<MVertex*,Matrix> > Frame_field::field;
std::map<int, Matrix> Frame_field::crossField;
std::map<MEdge, double, Less_Edge> Frame_field::crossDist;
+Recombinator Frame_field::crossData;
#if defined(HAVE_ANN)
ANNpointArray Frame_field::duplicate;
ANNkd_tree* Frame_field::kd_tree;
+ANNpointArray Frame_field::annTreeData;
+ANNkd_tree* Frame_field::annTree;
+std::vector<int> Frame_field::vertIndices;
#endif
std::map<MVertex*,double> Size_field::boundary;
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index 87af72dce4..0f7a39ceca 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -6,39 +6,17 @@
// Contributor(s):
// Tristan Carrier
+#ifndef _DIRECTION_3D_H_
+#define _DIRECTION_3D_H_
+
#include "GFace.h"
#include "MEdge.h"
#include "MElementOctree.h"
#if defined(HAVE_ANN)
#include <ANN/ANN.h>
#endif
-
-class Matrix{
- private:
- double m11,m21,m31,m12,m22,m32,m13,m23,m33;
- public:
- Matrix();
- ~Matrix();
- void set_m11(double);
- void set_m21(double);
- void set_m31(double);
- void set_m12(double);
- void set_m22(double);
- void set_m32(double);
- void set_m13(double);
- void set_m23(double);
- void set_m33(double);
- double get_m11();
- double get_m21();
- double get_m31();
- double get_m12();
- double get_m22();
- double get_m32();
- double get_m13();
- double get_m23();
- double get_m33();
-};
-
+#include "yamakawa.h"
+#include "Matrix.h"
struct lowerThan {
bool operator() (const std::pair<int, Matrix>& lhs, const std::pair<int, Matrix>& rhs) const
@@ -49,12 +27,15 @@ class Frame_field{
private:
static std::map<MVertex*, Matrix> temp;
static std::vector<std::pair<MVertex*, Matrix> > field;
- //static std::set<std::pair<int, Matrix>, lowerThan> crossField;
static std::map<int, Matrix> crossField;
+ //static std::map<MVertex*, Matrix, MVertexLessThanNum> crossField;
static std::map<MEdge, double, Less_Edge> crossDist;
#if defined(HAVE_ANN)
static ANNpointArray duplicate;
static ANNkd_tree* kd_tree;
+ static ANNpointArray annTreeData;
+ static ANNkd_tree* annTree;
+ static std::vector<int> vertIndices;
#endif
Frame_field();
public:
@@ -67,10 +48,13 @@ class Frame_field{
static double get_ratio(GFace*,SPoint2);
static void print_field1();
static void print_field2(GRegion*);
- static void print_segment(SPoint3,SPoint3,std::ofstream&);
+ static void print_segment(SPoint3,SPoint3,double,double,std::ofstream&);
+ static Recombinator crossData;
static void init(GRegion* gr);
static double smooth(GRegion* gr);
static double findBarycenter(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter, Matrix &m0);
+ static void fillTreeVolume(GRegion* gr);
+ static Matrix findNearestCross(double x,double y,double z);
static void save(GRegion* gr, const std::string &filename);
static void save_energy(GRegion* gr, const std::string& filename);
static void save_dist(const std::string& filename);
@@ -114,3 +98,4 @@ class Nearest_point{
static void clear();
};
+#endif
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index 485c2c3ea6..25390cba6d 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -347,6 +347,7 @@ void Filler::treat_region(GRegion* gr){
RTree<Node*,double,3,double> rtree;
Frame_field::init_region(gr);
+ std::cout << "Smoothing: energy = " << Frame_field::smooth(gr) << std::endl;
Size_field::init_region(gr);
Size_field::solve(gr);
octree = new MElementOctree(gr->model());
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index 5d609616b6..541fa26529 100755
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -6,6 +6,10 @@
// Contributor(s):
// Tristan Carrier
+#ifndef _YAMAKAWA_H_
+#define _YAMAKAWA_H_
+
+
#include "GRegion.h"
#include <set>
@@ -316,3 +320,5 @@ class PostOp{
void build_vertex_to_pyramids(MElement*);
void erase_vertex_to_pyramids(MElement*);
};
+
+#endif
--
GitLab