diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 63a1aad70838ba979207647186fd4f72cec026dd..d57a348d5c7fec6513b1a35552f84dce2b7cba0d 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -41,7 +41,7 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound,
GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
: GEdge(m, tag, 0 , 0), _compound(compound)
{
- orderEdges ();
+ orderEdges();
int N = _compound.size();
v0 = _orientation[0] ? _compound[0]->getBeginVertex() : _compound[0]->getEndVertex();
v1 = _orientation[N-1] ? _compound[N-1]->getEndVertex() : _compound[N-1]->getBeginVertex();
@@ -56,6 +56,7 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
void GEdgeCompound::orderEdges()
{
+
std::vector<GEdge*> _c ;
std::list<GEdge*> edges ;
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 2ef74c6f4d62325b4f83290793f81f9f81dd5e57..31f5d3950c669ec39886ba8ec93d867e9adf24a7 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -199,6 +199,24 @@ class STensor3 {
static int _index[3][3] = {{0,1,2},{3,4,5},{6,7,8}};
return _index[i][j];
}
+ inline void set_m11(double x){ _val[0] = x; }
+ inline void set_m21(double x){ _val[1] = x; }
+ inline void set_m31(double x){ _val[2] = x; }
+ inline void set_m12(double x){ _val[3] = x; }
+ inline void set_m22(double x){ _val[4] = x; }
+ inline void set_m32(double x){ _val[5] = x; }
+ inline void set_m13(double x){ _val[6] = x; }
+ inline void set_m23(double x){ _val[7] = x; }
+ inline void set_m33(double x){ _val[8] = x; }
+ inline double get_m11(){ return _val[0]; }
+ inline double get_m21(){ return _val[1]; }
+ inline double get_m31(){ return _val[2]; }
+ inline double get_m12(){ return _val[3]; }
+ inline double get_m22(){ return _val[4]; }
+ inline double get_m32(){ return _val[5]; }
+ inline double get_m13(){ return _val[6]; }
+ inline double get_m23(){ return _val[7]; }
+ inline double get_m33(){ return _val[8]; }
void getMat(fullMatrix<double> &mat) const
{
for (int i = 0; i < 3; i++){
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 849b98494a510a57f2fbcfa0ec819cbc2d2d8354..c0c3cb5c2c2f095ff4c1b687bea48221659f8d1c 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -39,6 +39,8 @@
#include "Curvature.h"
#include "MElement.h"
#include "Context.h"
+#include "directions3D.h"
+#include "meshGRegion.h"
#if defined(HAVE_ANN)
#include <ANN/ANN.h>
@@ -769,7 +771,7 @@ void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundE
SVector3 pc(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
SVector3 nc = ps-pc;
if (dot(ne,nc) < 0) dir = 1;
- if (dir ==1 && hLayer > 0 ) hLayer *= -1.0;
+ if (dir == 1 && hLayer > 0 ) hLayer *= -1.0;
//int shift = 0;
//if(is_cut) shift = NE;
@@ -791,7 +793,7 @@ void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundE
//if double extruded layer
if (nbElemSecondLayer > 0){
std::vector<GEntity*> extrudedESec = current->extrudeBoundaryLayer
- (gfc, nbElemSecondLayer, hSecondLayer, dir ? 0 : 1, -5);
+ (eFace, nbElemSecondLayer, hSecondLayer, dir, -5);
GFace *eFaceSec = (GFace*) extrudedESec[0];
eFaceSec->addPhysicalEntity(9); //tag 9
current->setPhysicalName("outerSecondWall", 2, 9);//dim 2 tag 9
@@ -799,6 +801,7 @@ void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundE
eRegionSec->addPhysicalEntity(10); //tag 10
current->setPhysicalName("wallVolume", 3, 10);//dim 3 tag 10
}
+ //end double extrusion
for (unsigned int j = 2; j < extrudedE.size(); j++){
GFace *elFace = (GFace*) extrudedE[j];
@@ -1147,6 +1150,12 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
SVector3 dir_a = p1-p0; dir_a.normalize();
SVector3 dir_n(xyz[0]-p0.x(), xyz[1]-p0.y(), xyz[2]-p0.z()); dir_n.normalize();
SVector3 dir_cross = crossprod(dir_a,dir_n); dir_cross.normalize();
+ // if (ge->dim()==3){
+ // printf("coucou dim ==3 !!!!!!!!!!!!!!! \n");
+ // SVector3 d1,d2,d3;
+ // computeCrossField(x,y,z,d1,d2,d3);
+ // exit(1);
+ // }
//find discrete curvature at closest vertex
Curvature& curvature = Curvature::getInstance();
@@ -1239,6 +1248,20 @@ SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dir
return curvMetric;
}
+void Centerline::computeCrossField(double x,double y,double z,
+ SVector3 &d1, SVector3 &d2, SVector3 &d3)
+{
+
+ //Le code suivant permet d'obtenir les trois vecteurs orthogonaux en un point.
+ // int NumSmooth = 10;//CTX::instance()->mesh.smoothCrossField
+ // Matrix m2;
+ // if(NumSmooth)
+ // m2 = Frame_field::findCross(x,y,z);
+ // else
+ // m2 = Frame_field::search(x,y,z);
+
+}
+
void Centerline::printSplit() const
{
FILE * f = fopen("mySPLIT.pos","w");
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index 7c964863f57cba122a5278fce859a747eac5b443..abf5d57f9038ff58c3e2a61797a7612d76a6cbc7 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -121,6 +121,9 @@ class Centerline : public Field{
//anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
+ void computeCrossField(double x,double y,double z,
+ SVector3 &d1, SVector3 &d2, SVector3 &d3);
+
//import the 1D mesh of the centerlines (in vtk format)
//and fill the vector of lines
void importFile(std::string fileName);
diff --git a/Mesh/Levy3D.cpp b/Mesh/Levy3D.cpp
index 143f713ee209d8bfa843af4eedb9f0f9ce27345b..9159e439718c8b308160e8049da4e4a734f1406d 100644
--- a/Mesh/Levy3D.cpp
+++ b/Mesh/Levy3D.cpp
@@ -874,9 +874,7 @@ double LpCVT::get_size(double x,double y,double z){
Tensor LpCVT::get_tensor(double x,double y,double z){
Tensor t;
- Matrix m;
-
- m = Frame_field::search(x,y,z);
+ STensor3 m = Frame_field::search(x,y,z);
t.set_t11(m.get_m11());
t.set_t21(m.get_m12());
diff --git a/Mesh/Matrix.h b/Mesh/Matrix.h
index 81ef1227aac8e05f7aa3137e56f9ce36e8b7a3a6..83de7598a469e14855cb9af528511cb4b028f31a 100644
--- a/Mesh/Matrix.h
+++ b/Mesh/Matrix.h
@@ -11,75 +11,41 @@
class Matrix{
private:
- double m11,m21,m31,m12,m22,m32,m13,m23,m33;
- public:
+ double _val[9];
+ public:
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;
+ /* 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; */
+ _val[0] = _val[4] = _val[8] = 1.0;
+ _val[1] = _val[2] = _val[3] = 0.0;
+ _val[5] = _val[6] = _val[7] = 0.0;
}
~Matrix(){}
- void set_m11(double x){ m11 = x; }
- void set_m21(double x){ m21 = x; }
- void set_m31(double x){ m31 = x; }
- void set_m12(double x){ m12 = x; }
- void set_m22(double x){ m22 = x; }
- void set_m32(double x){ m32 = x; }
- void set_m13(double x){ m13 = x; }
- void set_m23(double x){ m23 = x; }
- void set_m33(double x){ m33 = x; }
- double get_m11(){ return m11; }
- double get_m21(){ return m21; }
- double get_m31(){ return m31; }
- double get_m12(){ return m12; }
- double get_m22(){ return m22; }
- double get_m32(){ return m32; }
- double get_m13(){ return m13; }
- double get_m23(){ return m23; }
- double get_m33(){ return m33; }
+ void set_m11(double x){ _val[0] = x; }
+ void set_m21(double x){ _val[1] = x; }
+ void set_m31(double x){ _val[2] = x; }
+ void set_m12(double x){ _val[3] = x; }
+ void set_m22(double x){ _val[4] = x; }
+ void set_m32(double x){ _val[5] = x; }
+ void set_m13(double x){ _val[6] = x; }
+ void set_m23(double x){ _val[7] = x; }
+ void set_m33(double x){ _val[8] = x; }
+ double get_m11(){ return _val[0]; }
+ double get_m21(){ return _val[1]; }
+ double get_m31(){ return _val[2]; }
+ double get_m12(){ return _val[3]; }
+ double get_m22(){ return _val[4]; }
+ double get_m32(){ return _val[5]; }
+ double get_m13(){ return _val[6]; }
+ double get_m23(){ return _val[7]; }
+ double get_m33(){ return _val[8]; }
};
-/*
-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; }
-
-
-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; }
-*/
-
-
#endif
diff --git a/Mesh/cross3D.h b/Mesh/cross3D.h
index ce0c88e424049181b574357feaf2cd0753b58adb..c3289e88a1ad68e12b811001a589a6f95f36afd4 100644
--- a/Mesh/cross3D.h
+++ b/Mesh/cross3D.h
@@ -12,8 +12,6 @@
#include <ostream>
using std::ostream;
-#include "Matrix.h"
-
/* Defined in SVector3.h */
/* double angle(const SVector3 v, const SVector3 w) { */
@@ -112,8 +110,8 @@ public:
b = (crossprod(a,Ex).norm() < 1e-2) ? Ey : Ex;
scnd = crossprod(crossprod(frst, b), frst).unit();
}
- cross3D(const Matrix &x){
- Matrix m = x;
+ cross3D(const STensor3 &x){
+ STensor3 m = x;
SVector3 a(m.get_m11(), m.get_m21(), m.get_m31());
SVector3 b(m.get_m12(), m.get_m22(), m.get_m32());
frst = a.unit();
@@ -318,8 +316,8 @@ Qtn cross3D::rotationToOnSurf(const cross3D &y) const{
return R2;
}
-Matrix convert(const cross3D x){
- Matrix m;
+STensor3 convert(const cross3D x){
+ STensor3 m;
SVector3 v;
v = x.getFrst() ; m.set_m11(v[0]); m.set_m21(v[1]); m.set_m31(v[2]);
v = x.getScnd() ; m.set_m12(v[0]); m.set_m22(v[1]); m.set_m32(v[2]);
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index b453e9a40d59b61211e32431a790aaa9a4705011..17f951430c1de952f682ced70bd2ad6c4b72724e 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -20,6 +20,8 @@
#include "linearSystemFull.h"
#endif
+
+
/****************class Frame_field****************/
Frame_field::Frame_field(){}
@@ -31,7 +33,7 @@ void Frame_field::init_region(GRegion* gr){
MVertex* vertex;
GFace* gf;
std::list<GFace*> faces;
- Matrix m;
+ STensor3 m(1.0);
Nearest_point::init_region(gr);
@@ -48,13 +50,13 @@ void Frame_field::init_region(GRegion* gr){
duplicate = annAllocPts(temp.size(),3);
index = 0;
- for(std::map<MVertex*,Matrix>::iterator it=temp.begin(); it != temp.end(); it++){
+ for(std::map<MVertex*,STensor3>::iterator it=temp.begin(); it != temp.end(); it++){
vertex = it->first;
m = it->second;
duplicate[index][0] = vertex->x();
duplicate[index][1] = vertex->y();
duplicate[index][2] = vertex->z();
- field.push_back(std::pair<MVertex*,Matrix>(vertex,m));
+ field.push_back(std::pair<MVertex*,STensor3>(vertex,m));
index++;
}
@@ -63,7 +65,7 @@ void Frame_field::init_region(GRegion* gr){
}
void Frame_field::init_face(GFace* gf){
- // Fills the auxiliary std::map "temp" with a pair <vertex, Matrix>
+ // Fills the auxiliary std::map "temp" with a pair <vertex, STensor3>
// for each vertex of the face gf.
unsigned int i;
int j;
@@ -74,7 +76,7 @@ void Frame_field::init_face(GFace* gf){
MVertex* vertex;
MElement* element;
MElementOctree* octree;
- Matrix m;
+ STensor3 m(1.0);
backgroundMesh::set(gf);
octree = backgroundMesh::current()->get_octree();
@@ -102,7 +104,7 @@ void Frame_field::init_face(GFace* gf){
ok = translate(gf,octree,vertex,SPoint2(average_x,average_y),v1,v2);
}
else{
- ok = improved_translate(gf,vertex,v1,v2);
+ ok = improved_translate(gf,vertex,v1,v2);
}
if(ok){
@@ -119,7 +121,7 @@ void Frame_field::init_face(GFace* gf){
m.set_m13(v3.x());
m.set_m23(v3.y());
m.set_m33(v3.z());
- temp.insert(std::pair<MVertex*,Matrix>(vertex,m));
+ temp.insert(std::pair<MVertex*,STensor3>(vertex,m));
}
}
}
@@ -212,7 +214,7 @@ bool Frame_field::improved_translate(GFace* gf,MVertex* vertex,SVector3& v1,SVec
return 1;
}
-Matrix Frame_field::search(double x,double y,double z){
+STensor3 Frame_field::search(double x,double y,double z){
// Determines the frame/cross at location (x,y,z)
int index=0;
#if defined(HAVE_ANN)
@@ -239,7 +241,7 @@ Matrix Frame_field::search(double x,double y,double z){
return field[index].second;
}
-Matrix Frame_field::combine(double x,double y,double z){
+STensor3 Frame_field::combine(double x,double y,double z){
// Determines the frame/cross at location (x,y,z)
// Alternative to Frame_field::search
bool ok;
@@ -247,7 +249,7 @@ Matrix Frame_field::combine(double x,double y,double z){
SVector3 vec,other;
SVector3 vec1,vec2,vec3;
SVector3 final1,final2;
- Matrix m,m2;
+ STensor3 m(1.0),m2(1.0);
m = search(x,y,z);
m2 = m;
@@ -323,8 +325,8 @@ void Frame_field::print_field1(){
SPoint3 p;
SPoint3 p1,p2,p3,p4,p5,p6;
MVertex* vertex;
- std::map<MVertex*,Matrix>::iterator it;
- Matrix m;
+ std::map<MVertex*,STensor3>::iterator it;
+ STensor3 m(1.0);
k = 0.05;
std::ofstream file("frame1.pos");
@@ -373,7 +375,7 @@ void Frame_field::print_field2(GRegion* gr){
SPoint3 p1,p2,p3,p4,p5,p6;
MVertex* vertex;
MElement* element;
- Matrix m;
+ STensor3 m(1.0);
k = 0.05;
std::ofstream file("frame2.pos");
@@ -456,7 +458,7 @@ int Frame_field::build_vertex_to_vertices(GEntity* gr, int onWhat, bool initiali
unsigned int n = pElem->getNumVertices();
for(unsigned int j=0; j<n; j++){
MVertex* pVertex = pElem->getVertex(j);
- if(pVertex->onWhat()->dim() != onWhat) continue;
+ if(onWhat >0 && pVertex->onWhat()->dim() != onWhat) continue;
std::map<MVertex*,std::set<MVertex*> >::iterator it = vertex_to_vertices.find(pVertex);
if(it != vertex_to_vertices.end()){
@@ -581,10 +583,10 @@ void Frame_field::initFace(GFace* gf){
Area += crossprod(V2-V0,V1-V0);
}
Area.normalize(); // average normal over neighbour face elements
- Matrix m = convert(cross3D(Area));
- std::map<MVertex*, Matrix>::iterator iter = crossField.find(it->first);
+ STensor3 m = convert(cross3D(Area));
+ std::map<MVertex*, STensor3>::iterator iter = crossField.find(it->first);
if(iter == crossField.end())
- crossField.insert(std::pair<MVertex*, Matrix>(it->first,m));
+ crossField.insert(std::pair<MVertex*, STensor3>(it->first,m));
else
crossField[it->first] = m;
}
@@ -594,8 +596,24 @@ void Frame_field::initFace(GFace* gf){
// compute cumulative cross-data vertices x elements for the whole contour of gf
std::list<GEdge*> edges = gf->edges();
vertex_to_elements.clear();
- for( std::list<GEdge*>::const_iterator it=edges.begin(); it!=edges.end(); it++)
+ //Replace edges by their compounds
+ std::set<GEdge*> mySet;
+ std::list<GEdge*>::iterator it = edges.begin();
+ while(it != edges.end()){
+ if((*it)->getCompound()){
+ GEdge *gec = (GEdge*)(*it)->getCompound();
+ mySet.insert(gec);
+ }
+ else{
+ mySet.insert(*it);
+ }
+ ++it;
+ }
+ edges.clear();
+ edges.insert(edges.begin(), mySet.begin(), mySet.end());
+ for( std::list<GEdge*>::const_iterator it=edges.begin(); it!=edges.end(); it++){
build_vertex_to_elements(*it,false);
+ }
// align crosses of the contour of "gf" with the tangent to the contour
for(std::map<MVertex*, std::set<MElement*> >::const_iterator it
@@ -611,13 +629,13 @@ void Frame_field::initFace(GFace* gf){
SVector3 tangent = edg1.scaledTangent() + edg2.scaledTangent();
tangent.normalize();
- std::map<MVertex*, Matrix>::iterator iter = crossField.find(pVertex);
+ std::map<MVertex*, STensor3>::iterator iter = crossField.find(pVertex);
if(iter == crossField.end()) {
std::cout << "This should not happen: cross not found 1" << std::endl;
exit(1);
}
cross3D x(cross3D((*iter).second));
- Matrix m = convert(cross3D(x.getFrst(), tangent));
+ STensor3 m = convert(cross3D(x.getFrst(), tangent));
crossField[pVertex] = m;
}
@@ -642,7 +660,7 @@ void Frame_field::initFace(GFace* gf){
MVertex* pVertex0 = it->first;
if(pVertex0->onWhat()->dim() != 2) continue;
- std::map<MVertex*, Matrix>::iterator iter = crossField.find(pVertex0);
+ std::map<MVertex*, STensor3>::iterator iter = crossField.find(pVertex0);
cross3D y;
if(iter == crossField.end()){
std::cout << "This should not happen: cross not found 2" << std::endl;
@@ -663,9 +681,9 @@ void Frame_field::initFace(GFace* gf){
}
cross3D x = cross3D((*iter).second); //nearest cross on contour, x.getFrst() is the normal
- //Matrix m = convert(cross3D(x.getFrst(),y.getScnd()));
+ //STensor3 m = convert(cross3D(x.getFrst(),y.getScnd()));
SVector3 v1 = y.getFrst();
- Matrix m = convert( y.rotate(conj(x.rotationToOnSurf(y)))); //rotation around the normal
+ STensor3 m = convert( y.rotate(conj(x.rotationToOnSurf(y)))); //rotation around the normal
SVector3 v2 = y.getFrst();
if(fabs(angle(v1,v2)) > 1e-8){
std::cout << "This should not happen: rotation affected the normal" << std::endl;
@@ -694,22 +712,22 @@ void Frame_field::initRegion(GRegion* gr, int n){
// Find the index of the nearest cross on the contour, using annTree
int index = findAnnIndex(pVertex0->point());
MVertex* pVertex = listVertices[index];
- Matrix m = crossField[pVertex];
- crossField.insert(std::pair<MVertex*, Matrix>(pVertex0,m));
+ STensor3 m = crossField[pVertex];
+ crossField.insert(std::pair<MVertex*, STensor3>(pVertex0,m));
}
deleteAnnData();
buildAnnData(gr,3);
}
-Matrix Frame_field::findCross(double x, double y, double z){
+STensor3 Frame_field::findCross(double x, double y, double z){
int index = Frame_field::findAnnIndex(SPoint3(x,y,z));
MVertex* pVertex = Frame_field::listVertices[index];
return crossField[pVertex];
}
-double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter, Matrix &m0){
+double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter, STensor3 &m0){
double theta, gradient, energy;
- Matrix m;
+ STensor3 m;
SVector3 axis;
bool debug = false;
@@ -726,7 +744,7 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
MVertex *pVertex = *it;
if(pVertex->getNum() == pVertex0->getNum())
std::cout << "This should not happen!" << std::endl;
- std::map<MVertex*, Matrix>::const_iterator itB = crossField.find(pVertex);
+ std::map<MVertex*, STensor3>::const_iterator itB = crossField.find(pVertex);
if(itB == crossField.end()){ // not found
std::cout << "This should not happen: cross not found 4" << std::endl; exit(1);
}
@@ -785,7 +803,7 @@ double Frame_field::smoothRegion(GRegion *gr, int n){
}
double Frame_field::smooth(){
- Matrix m, m0;
+ STensor3 m(1.0), m0(1.0);
double enew, eold;
double energy = 0;
@@ -794,7 +812,7 @@ double Frame_field::smooth(){
//MVertex* pVertex0 = iter->first;
SVector3 T(0, 0, 0);
- std::map<MVertex*, Matrix>::iterator itA = crossField.find(iter->first);
+ std::map<MVertex*, STensor3>::iterator itA = crossField.find(iter->first);
if(itA == crossField.end()){
std::cout << "This should not happen" << std::endl; exit(1);
}
@@ -814,7 +832,7 @@ double Frame_field::smooth(){
return energy;
}
-void Frame_field::save(const std::vector<std::pair<SPoint3, Matrix> > data,
+void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
const std::string& filename){
const cross3D origin(SVector3(1,0,0), SVector3(0,1,0));
SPoint3 p1;
@@ -823,7 +841,7 @@ void Frame_field::save(const std::vector<std::pair<SPoint3, Matrix> > data,
file << "View \"cross field\" {\n";
for(unsigned int i=0; i<data.size(); i++){
SPoint3 p = data[i].first;
- Matrix m = data[i].second;
+ STensor3 m = data[i].second;
double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin)), val2=val1;
p1 = SPoint3(p.x() + k*m.get_m11(),
p.y() + k*m.get_m21(),
@@ -854,25 +872,133 @@ void Frame_field::save(const std::vector<std::pair<SPoint3, Matrix> > data,
file.close();
}
+void Frame_field::recur_connect_vert(MVertex *v,
+ STensor3 &cross,
+ std::multimap<MVertex*,MVertex*> &v2v,
+ std::set<MVertex*> &touched){
+
+ if (touched.find(v) != touched.end()) return;
+ touched.insert(v);
+ //std::map<MVertex*, STensor3>::iterator iterv = crossField.find(v);
+ //STensor3 cross = iterv->second;
+
+ for (std::multimap <MVertex*,MVertex*>::iterator it = v2v.lower_bound(v);
+ it != v2v.upper_bound(v) ; ++it){
+ MVertex *nextV = it->second;
+
+ //change orientation
+ ///------------------
+ printf("*********************** Changing orientation \n");
+
+ //compute dot product (N0,R0,A0)^T dot (Ni,Ri,Ai)
+ //where N,R,A are the normal, radial and axial direction s
+ std::map<MVertex*, STensor3>::iterator iter = crossField.find(nextV);
+ STensor3 nextCross = iter->second;
+ STensor3 crossT = cross.transpose();
+ STensor3 prod = crossT.operator*=(nextCross);
+ cross.print("cross ");
+ nextCross.print("nextCross ");
+ prod.print("product");
+ fullMatrix<double> mat(3,3); prod.getMat(mat);
+
+ //find biggest dot product
+ fullVector<int> Id(3);
+ for (int i = 0; i < 3; i++){
+ double maxVal = 0.0;
+ for (int j = 0; j < 3; j++){
+ double val = fabs(mat(i,j));
+ if( val > maxVal ){
+ maxVal = val;
+ Id(i) = j;
+ }
+ }
+ }
+
+ if (Id(0) +Id(1)+ Id(2) != 3 || (Id(0) == 1 && Id(1)==1 && Id(2)==1)) {
+ std::cout << "This should not happen: sum should be 0+1+2" << std::endl;
+ printf("Id =%d %d %d \n", Id(0), Id(1), Id(2));
+ //exit(1);
+ return;
+ }
+
+ //create new cross
+ fullMatrix<double> newmat(3,3);
+ for (int i = 0; i < 3; i++){
+ for (int j = 0; j < 3; j++){
+ newmat(i,j) = nextCross(Id(j),i);
+ }
+ }
+
+ STensor3 newcross;
+ newcross.setMat(newmat);
+ crossField[iter->first] = newcross;
+ newcross.print("newcross");
+
+ //continue recursion
+ recur_connect_vert (nextV, newcross, v2v,touched);
+ }
+
+}
+void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
+
+ printf("continuous cross field \n");
+
+ //start from a vertex of a face
+ std::list<GFace*> faces = gr->faces();
+ bool foundFace = false;
+ for(std::list<GFace*>::const_iterator iter=faces.begin(); iter!=faces.end(); iter++){
+ if (*iter == gf){
+ foundFace = true;
+ break;
+ }
+ }
+ if (!foundFace){
+ std::cout << "This should not happen: face does not belong to region" << std::endl;
+ exit(1);
+ }
+
+ //build connectivity
+ build_vertex_to_vertices(gr, -1, true);
+ std::multimap<MVertex*,MVertex*> v2v;
+ for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
+ = vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
+ MVertex *v = iter->first;
+ std::set<MVertex*> mySet = iter->second;
+ for (std::set<MVertex*>::iterator it = mySet.begin(); it!=mySet.end(); ++it){
+ v2v.insert(std::make_pair(v,*it));
+ }
+ }
+
+ //recursive loop
+ MVertex *beginV = gf->getMeshVertex(0);
+ std::set<MVertex*> touched;
+ std::map<MVertex*, STensor3>::iterator iter = crossField.find(beginV);
+ STensor3 bCross = iter->second;
+ recur_connect_vert (beginV,bCross,v2v,touched);
+
+}
+
void Frame_field::saveCrossField(const std::string& filename, double scale, bool full){
const cross3D origin(SVector3(1,0,0), SVector3(0,1,0));
SPoint3 p1;
double k = scale;
std::ofstream file(filename.c_str());
file << "View \"cross field\" {\n";
- for(std::map<MVertex*, Matrix>::const_iterator it = crossField.begin();
+ for(std::map<MVertex*, STensor3>::const_iterator it = crossField.begin();
it != crossField.end(); it++){
SPoint3 p = it->first->point();
- Matrix m = it->second;
+ STensor3 m = it->second;
double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin))*180./M_PI, val2=val1;
+ //double val1=1.0, val2=1.0;
p1 = SPoint3(p.x() + k*m.get_m11(),
p.y() + k*m.get_m21(),
- p.z() + k*m.get_m31());
+ p.z() + k*m.get_m31());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() - k*m.get_m11(),
p.y() - k*m.get_m21(),
p.z() - k*m.get_m31());
if(full) print_segment(p,p1,val1,val2,file);
+ //val1=2.0; val2=2.0;
p1 = SPoint3(p.x() + k*m.get_m12(),
p.y() + k*m.get_m22(),
p.z() + k*m.get_m32());
@@ -881,6 +1007,7 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
p.y() - k*m.get_m22(),
p.z() - k*m.get_m32());
if(full) print_segment(p,p1,val1,val2,file);
+ //val1=3.0; val2=3.0;
p1 = SPoint3(p.x() + k*m.get_m13(),
p.y() + k*m.get_m23(),
p.z() + k*m.get_m33());
@@ -1563,9 +1690,9 @@ void Nearest_point::clear(){
/****************static declarations****************/
-std::map<MVertex*,Matrix> Frame_field::temp;
-std::vector<std::pair<MVertex*,Matrix> > Frame_field::field;
-std::map<MVertex*, Matrix> Frame_field::crossField;
+std::map<MVertex*,STensor3> Frame_field::temp;
+std::vector<std::pair<MVertex*,STensor3> > Frame_field::field;
+std::map<MVertex*, STensor3> Frame_field::crossField;
std::map<MEdge, double, Less_Edge> Frame_field::crossDist;
std::map<MVertex*,std::set<MVertex*> > Frame_field::vertex_to_vertices;
std::map<MVertex*,std::set<MElement*> > Frame_field::vertex_to_elements;
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index 6316bf4fdb256b354f322379bac9a5a9f2e71810..ce418ef83456a78b7ee3090dbb09199ad51d4394 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -16,18 +16,18 @@
#include <ANN/ANN.h>
#endif
#include "yamakawa.h"
-#include "Matrix.h"
+#include "STensor3.h"
struct lowerThan {
- bool operator() (const std::pair<int, Matrix>& lhs, const std::pair<int, Matrix>& rhs) const
+ bool operator() (const std::pair<int, STensor3>& lhs, const std::pair<int, STensor3>& rhs) const
{return lhs.first < rhs.first;}
};
class Frame_field{
private:
- static std::map<MVertex*, Matrix> temp;
- static std::vector<std::pair<MVertex*, Matrix> > field;
- static std::map<MVertex*, Matrix> crossField;
+ static std::map<MVertex*, STensor3> temp;
+ static std::vector<std::pair<MVertex*, STensor3> > field;
+ static std::map<MVertex*, STensor3> crossField;
static std::map<MEdge, double, Less_Edge> crossDist;
static std::vector<MVertex*> listVertices;
#if defined(HAVE_ANN)
@@ -42,8 +42,8 @@ class Frame_field{
static void init_face(GFace*);
static bool translate(GFace*,MElementOctree*,MVertex*,SPoint2,SVector3&,SVector3&);
static bool improved_translate(GFace*,MVertex*,SVector3&,SVector3&);
- static Matrix search(double,double,double);
- static Matrix combine(double,double,double);
+ static STensor3 search(double,double,double);
+ static STensor3 combine(double,double,double);
static bool inside_domain(MElementOctree*,double,double);
static double get_ratio(GFace*,SPoint2);
static void print_field1();
@@ -57,16 +57,18 @@ class Frame_field{
static int buildAnnData(GEntity* ge, int dim);
static void deleteAnnData();
static int findAnnIndex(SPoint3 p);
- static Matrix findCross(double x, double y, double z);
+ static STensor3 findCross(double x, double y, double z);
static void initFace(GFace* gf);
static void initRegion(GRegion* gr, int n);
static double smoothFace(GFace *gf, int n);
static double smoothRegion(GRegion *gr, int n);
static double smooth();
- static double findBarycenter(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter, Matrix &m0);
- static void save(const std::vector<std::pair<SPoint3, Matrix> > data,
+ static double findBarycenter(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter, STensor3 &m0);
+ static void save(const std::vector<std::pair<SPoint3, STensor3> > data,
const std::string& filename);
static void saveCrossField(const std::string& filename, double scale, bool full=true);
+ static void continuousCrossField(GRegion *gr, GFace *gf);
+ static void recur_connect_vert(MVertex *v,STensor3 &cross, std::multimap<MVertex*,MVertex*> &v2v, std::set<MVertex*> &touched);
static void save_energy(GRegion* gr, const std::string& filename);
static void save_dist(const std::string& filename);
static void checkAnnData(GEntity* ge, const std::string& filename);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 3c9fc1b91e9630a152c9060b346f421d8ed1277c..6e254cf5f4ef7732835888aa1a021e0e20b677e0 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -27,6 +27,7 @@
#include "meshGRegionMMG3D.h"
#include "simple3D.h"
#include "Levy3D.h"
+#include "directions3D.h"
#if defined(HAVE_ANN)
#include "ANN/ANN.h"
@@ -205,125 +206,6 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
}
-void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles)
-{
- std::vector<MTetrahedron*> elements = gr->tetrahedra;
- std::list<GFace*> allFaces = gr->faces();
-
- std::set<SPoint3> candidates;
- std::set<SPoint3> seeds;
-
- printf("computing box and barycenters\n");
- double Dmax = -1.0;
- std::list<GFace*>::iterator itf = allFaces.begin();
- for(; itf != allFaces.end(); itf ++){
- SBoundingBox3d bbox = (*itf)->bounds();
- SVector3 dd(bbox.max(), bbox.min());
- //if( (*itf)->tag()==5904 || (*itf)->tag()==5902 || (*itf)->tag()==5906) {
- //if( (*itf)->tag()==6 || (*itf)->tag()==28){
- if( (*itf)->tag() !=1 ) {
- Dmax = std::max(Dmax,norm(dd));
- seeds.insert(bbox.center());
- }
- }
- printf("Dmax =%g \n", Dmax);
-
- printf("printing skeleton nodes \n");
- FILE *outfile;
- outfile = fopen("skeletonNodes.pos", "w");
- fprintf(outfile, "View \"skeleton nodes\" {\n");
- for(unsigned int i = 0; i < elements.size(); i++){
- MTetrahedron *ele = elements[i];
- SPoint3 pc = ele->circumcenter();
- double x[3] = {pc.x(), pc.y(), pc.z()};
- double uvw[3];
- ele->xyz2uvw(x, uvw);
-
- std::set<SPoint3>::const_iterator it2 = voronoiPoles.find(pc);
-
- if(ele->isInside(uvw[0], uvw[1], uvw[2]) &&
- it2 != voronoiPoles.end()){
- double radius = ele->getCircumRadius();
- if(radius > Dmax/10.) {
- candidates.insert(pc);
- fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
- pc.x(), pc.y(), pc.z(), radius);
- }
- }
- }
- fprintf(outfile,"};\n");
- fclose(outfile);
-
- printf("Ann computation of neighbours and writing edges\n");
-
-#if defined(HAVE_ANN)
- FILE *outfile2;
- outfile2 = fopen("skeletonEdges.pos", "w");
- fprintf(outfile2, "View \"skeleton edges\" {\n");
-
- ANNkd_tree *_kdtree;
- ANNpointArray _zeronodes;
- ANNidxArray _index;
- ANNdistArray _dist;
-
- std::set<SPoint3>::iterator itseed = seeds.begin();
- SPoint3 beginPt=*itseed;
- seeds.erase(itseed);
- itseed = seeds.begin();
- for(; itseed != seeds.end(); itseed++){
- printf("seed =%g %g %g \n", (*itseed).x(), (*itseed).y(), (*itseed).z());
- candidates.insert(*itseed);
- }
- printf("begin seed =%g %g %g \n", beginPt.x(), beginPt.y(), beginPt.z());
-
- double color = 1.;
- while(candidates.size()>0){
-
- _zeronodes = annAllocPts(candidates.size(), 3);
- std::set<SPoint3>::iterator itset = candidates.begin();
- int i=0;
- for(; itset != candidates.end(); itset++){
- _zeronodes[i][0] = (*itset).x();
- _zeronodes[i][1] = (*itset).y();
- _zeronodes[i][2] = (*itset).z();
- i++;
- }
- _kdtree = new ANNkd_tree(_zeronodes, candidates.size(), 3);
- _index = new ANNidx[1];
- _dist = new ANNdist[1];
-
- double xyz[3] = {beginPt.x(), beginPt.y(), beginPt.z()};
- _kdtree->annkSearch(xyz, 1, _index, _dist);
- SPoint3 endPt( _zeronodes[_index[0]][0], _zeronodes[_index[0]][1], _zeronodes[_index[0]][2]);
- fprintf(outfile2,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
- beginPt.x(), beginPt.y(), beginPt.z(),
- endPt.x(), endPt.y(), endPt.z(),
- color, color);
-
- std::set<SPoint3>::iterator itse=seeds.find(endPt);
- std::set<SPoint3>::iterator its=candidates.find(endPt);
- if(itse != seeds.end()){
- printf("found seed =%g %g %g \n", endPt.x(), endPt.y(), endPt.z());
- seeds.erase(itse);
- beginPt = *(seeds.begin());
- std::set<SPoint3>::iterator itsee=candidates.find(beginPt);
- if (itsee != candidates.end()) candidates.erase(itsee);
- color=color*2.;
- }
- else beginPt=endPt;
-
- if(its != candidates.end()) candidates.erase(its);
-
- delete _kdtree;
- annDeallocPts(_zeronodes);
- delete [] _index;
- delete [] _dist;
- }
-
- fprintf(outfile2,"};\n");
- fclose(outfile2);
-#endif
-}
void getBoundingInfoAndSplitQuads(GRegion *gr,
std::map<MFace,GEntity*,Less_Face> &allBoundingFaces,
@@ -734,15 +616,9 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
- // EMI Voronoi for centerlines
- if (Msg::GetVerbosity() == 20) {
- std::set<SPoint3> candidates;
- printVoronoi(gr, candidates);
- skeletonFromVoronoi(gr, candidates);
- }
-
modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr);
+
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
bowyerWatsonFrontalLayers(gr, false);
@@ -756,6 +632,24 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
insertVerticesInRegion(gr);
}
+ // //emi test frame field
+ // int NumSmooth = 10;//CTX::instance()->mesh.smoothCrossField
+ // std::cout << "NumSmooth = " << NumSmooth << std::endl;
+ // if(NumSmooth && (gr->dim() == 3)){
+ // double scale = gr->bounds().diag()*1e-2;
+ // Frame_field::initRegion(gr,NumSmooth);
+ // Frame_field::saveCrossField("cross0.pos",scale);
+ // Frame_field::smoothRegion(gr,NumSmooth);
+ // Frame_field::saveCrossField("cross1.pos",scale);
+ // GFace *gf = GModel::current()->getFaceByTag(2);
+ // Frame_field::continuousCrossField(gr,gf);
+ // Frame_field::saveCrossField("cross2.pos",scale);
+ // }
+ // Frame_field::init_region(gr);
+ // Frame_field::clear();
+ // exit(1);
+ // //fin test emi
+
if (sqr.buildPyramids (gr->model())){
// relocate vertices if pyramids
for(unsigned int k = 0; k < regions.size(); k++){
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index 612bc149f87603496d2c5d32906cab51ca157149..e729a105101d391bbb018ce213b65689118769d4 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -500,7 +500,7 @@ void Filler::treat_region(GRegion* gr){
Metric Filler::get_metric(double x,double y,double z){
Metric m;
- Matrix m2;
+ STensor3 m2;
if(CTX::instance()->mesh.smoothCrossField){
m2 = Frame_field::findCross(x,y,z);
}
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index dafab665c4533c3d17554c149ec7df4e358488c5..7b84f7f498e61d78818163961dd3800ec3b0c710 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -1,3 +1,4 @@
+Mesh.LcIntegrationPrecision = 1.e-2;
lc = .033 ;
lc2 = 2.2 ;
lc3 = .033 ;
@@ -223,21 +224,19 @@ Line(8) = {1003,1000};
Line Loop(9) = {6,7,8,5};
Line Loop(10) = {2,3,4,1};
Plane Surface(11) = {9,10};
-
-//Physical Surface(11)={11};
-//Point(9999) = {0.6,0,0,1};
+//Recombine Surface {11};
Field[2] = BoundaryLayer;
//Field[2].NodesList = {1};
//Field[2].EdgesList = {1,2,3,4};
Field[2].EdgesList = {1,2,3,4};
Field[2].hfar = 1.5;
-Field[2].hwall_n = 0.001;
+Field[2].hwall_n = 0.0008;
Field[2].hwall_t = 0.01;
-Field[2].ratio = 1.3;
+Field[2].ratio = 1.2;
Field[2].thickness = .03;
Field[2].IntersectMetrics = 1;
-//Background Field = 2;
+Background Field = 2;
Field[1] = Box;
Field[1].VIn = 0.01;
@@ -251,6 +250,6 @@ Field[1].ZMin = -1;
Field[3] = MinAniso;
Field[3].FieldsList = {1, 2};
-Background Field = 2;
-//BoundaryLayer Field = 2;
-//Recombine Surface {11};
+//Background Field = 2;
+BoundaryLayer Field = 2;
+
diff --git a/benchmarks/centerlines/cyl_centerlines.geo b/benchmarks/centerlines/cyl_centerlines.geo
index 09309bff8c20d3dd7e6d8c70984327040b9aeda8..b7f307b0634cd3781a787d6f6395ceed39d7f991 100644
--- a/benchmarks/centerlines/cyl_centerlines.geo
+++ b/benchmarks/centerlines/cyl_centerlines.geo
@@ -1,4 +1,4 @@
-Mesh.Algorithm = 7; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 1; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D
Mesh.LcIntegrationPrecision = 1.e-3;