diff --git a/Mesh/cross3D.h b/Mesh/cross3D.h
index 26648e250c0c2b08369b8b8a77ee187119aed029..ce0c88e424049181b574357feaf2cd0753b58adb 100644
--- a/Mesh/cross3D.h
+++ b/Mesh/cross3D.h
@@ -14,9 +14,6 @@ using std::ostream;
#include "Matrix.h"
-//double max(const double a, const double b) { return (b>a)?b:a;}
-double min(const double a, const double b) { return (b<a)?b:a; }
-double squ(const double a) { return a*a; }
/* Defined in SVector3.h */
/* double angle(const SVector3 v, const SVector3 w) { */
@@ -35,7 +32,7 @@ public:
Qtn(){};
~Qtn(){};
Qtn(const SVector3 axis, const double theta = M_PI ){
- double temp = sin(0.5*theta);
+ double temp = sin(0.5 * theta);
v[0] = axis[0] * temp;
v[1] = axis[1] * temp;
v[2] = axis[2] * temp;
@@ -63,15 +60,18 @@ void Qtn::storeProduct(const Qtn &x, const Qtn &y) {
v[2] = a0*b1 - a1*b0 + a2*b3 + a3*b2;
v[3] =-a0*b0 - a1*b1 - a2*b2 + a3*b3;
}
-Qtn Qtn::operator *(const Qtn &x) const {
- Qtn y;
- y.storeProduct(*this, x);
- return y;
+Qtn Qtn::operator *(const Qtn &y) const {
+ Qtn x;
+ x.storeProduct(*this, y);
+ return x;
}
SVector3 eulerAxisFromQtn(const Qtn &x){
double temp = sqrt(1. - x[3]*x[3]);
- return SVector3(x[0]/temp, x[1]/temp, x[2]/temp);
+ if(temp< 1e-10)
+ return SVector3(1,0,0);
+ else
+ return SVector3(x[0]/temp, x[1]/temp, x[2]/temp);
}
double eulerAngleFromQtn(const Qtn &x){
@@ -97,11 +97,23 @@ class cross3D{
private:
SVector3 frst, scnd;
public:
+ cross3D() {
+ frst = SVector3(1,0,0);
+ scnd = SVector3(0,1,0);
+ }
cross3D(const SVector3 &a, const SVector3 &b){
frst = a.unit();
scnd = crossprod(crossprod(frst, b), frst).unit();
}
- cross3D(Matrix &m){
+ cross3D(const SVector3 &a) {
+ //if only a is given, b is an arbitrary vector not parallel to a
+ SVector3 b, Ex(1,0,0), Ey(0,1,0);
+ frst = a.unit();
+ b = (crossprod(a,Ex).norm() < 1e-2) ? Ey : Ex;
+ scnd = crossprod(crossprod(frst, b), frst).unit();
+ }
+ cross3D(const Matrix &x){
+ Matrix 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();
@@ -127,7 +139,8 @@ public:
scnd = im(R*scnd*conj(R));
return *this;
}
- Qtn rotationTo(const cross3D &x) const;
+ Qtn rotationTo(const cross3D &y) const;
+ Qtn rotationToOnSurf(const cross3D &y) const;
};
ostream& operator<< (ostream &os, const cross3D &x) {
@@ -135,7 +148,7 @@ ostream& operator<< (ostream &os, const cross3D &x) {
return os;
}
-cross3D cross3D::get(const int k) const{
+cross3D cross3D::get(int k) const{
SVector3 a, b;
switch(k){
case 0: a = getFrst() ; b = getScnd() ; break;
@@ -175,81 +188,136 @@ Qtn cross3D::rotationTo(const cross3D &y) const{
/* x.rotationTo(y) returns a quaternion R representing the rotation
such that y = R x, x = conj(R) y
eulerAngleFromQtn(R) = distance(x, y)
+ if onFace is true, only the rotation around y.frst (which is the normal) is returned
*/
- int ind1, ind2;
- double d, dmin, th1, th2;
- SVector3 n, w, axis;
+ double d, dmin, jmin, kmin, th1, th2;
+ SVector3 axis;
Qtn Rxy1, Rxy2;
- cross3D xx = get(0);
- dmin = angle(xx.frst, y.get(0).frst); ind1 = 0;
- for(unsigned int k=4 ; k<24; k=k+4){
- if((d = angle(xx.frst, y.get(k).frst)) < dmin) {
- ind1 = k;
- dmin = d;
- }
- }
-
- // y.get(ind1) is the attitude of y whose y.first minimizes the angle with x.first
+ cross3D xx = *this;
+ cross3D yy = y;
- axis = crossprod(xx.frst, y.get(ind1).frst);
- if(axis.norm() > 1e-8)
- axis.normalize();
- else {
- axis = SVector3(1,0,0);
- dmin = 0.;
+ //ind1 = 0; jmin=0; dmin = angle(xx.get(kmin).frst, vect[jmin]);
+ dmin = M_PI; jmin = kmin = 0;
+ for(int j=0 ; j<24; j=j+4){
+ for(int k=0 ; k<12; k=k+4){
+ if((d = angle(xx.get(j).frst, yy.get(k).frst)) < dmin) {
+ kmin = k;
+ jmin = j;
+ dmin = d;
+ }
+ }
}
+ xx = xx.get(jmin);
+ yy = yy.get(kmin);
+ // xx.frst and yy.frst now form the smallest angle among the 6^2 possible.
th1 = dmin;
if(th1 > 1.00001*acos(1./sqrt(3.))){
std::cout << "This should not happen: th1 = " << th1 << std::endl;
exit(1);
}
+ if(th1 > 1e-8)
+ axis = crossprod(xx.frst, yy.frst).unit();
+ else {
+ axis = SVector3(1,0,0);
+ th1 = 0.;
+ }
- Rxy1 = Qtn(axis, dmin);
+ Rxy1 = Qtn(axis, th1);
xx.rotate(Rxy1);
- dmin = angle(xx.scnd, y.get(ind1).scnd); ind2 = ind1;
- for(unsigned int k=1 ; k<4; k++){
- if((d = angle(xx.scnd, y.get(ind1 + k).scnd)) < dmin) {
- ind2 = ind1 + k;
+ dmin = M_PI;
+ for(int j=0 ; j<4; j++){
+ if((d = angle(xx.get(j).scnd, yy.scnd)) < dmin) {
+ jmin = j;
dmin = d;
}
}
-
- // y.get(ind2) is the attitude of y whose y.second minimizes the angle with xx.second
-
- axis = crossprod(xx.scnd, y.get(ind2).scnd);
- //axis = xx.first;
- if(axis.norm() > 1e-8)
- axis.normalize();
- else {
- axis = SVector3(1,0,0);
- dmin = 0.;
- }
+ xx = xx.get(jmin);
+ // xx.scnd and yy.scnd now form the smallest angle among the 4^2 possible.
th2 = dmin;
- if(th2 > M_PI/2.){
+ if(th2 > M_PI/4.){
std::cout << "This should not happen: th2 = " << th2 << std::endl;
exit(1);
}
- Rxy2 = Qtn(axis, dmin);
+ if(th2 > 1e-8)
+ axis = crossprod(xx.scnd, yy.scnd).unit();
+ else {
+ axis = SVector3(1,0,0);
+ th2 = 0.;
+ }
- xx.rotate(Rxy2);
+ Rxy2 = Qtn(axis, th2);
Qtn R = Rxy2*Rxy1;
+ // test
double theta = eulerAngleFromQtn(R);
- if(theta > 1.11){
+ if(theta > 1.07 /*0.988*/){ //
std::cout << "Ouch! th1 = " << th1 << " th2 = " << th2 << std::endl;
std::cout << "x = " << *this << std::endl;
std::cout << "y = " << y << std::endl;
- /* std::cout << "R = " << R << std::endl; */
+ std::cout << "R = " << R << std::endl;
std::cout << "u = " << eulerAngleFromQtn(R) << std::endl;
std::cout << "axis = " << eulerAxisFromQtn(R) << std::endl;
}
+
return R;
}
+Qtn cross3D::rotationToOnSurf(const cross3D &y) const{
+ /* this->frst and y.frst are assumed to be the normal to the face
+ R1 is the rotation such that R1(this->frst) = y.frst.
+ R2 is then the rotation such that R2 o R1(this->scnd) = y.scnd
+ */
+ double d, dmin, jmin, th1, th2;
+ SVector3 axis;
+
+ cross3D xx = *this;
+ cross3D yy = y;
+
+ th1 = angle(xx.frst, yy.frst);
+ if(th1 > 1e-8)
+ axis = crossprod(xx.frst, yy.frst).unit();
+ else {
+ axis = SVector3(1,0,0);
+ th1 = 0.;
+ }
+
+ Qtn R1 = Qtn(axis, th1);
+ xx.rotate(R1);
+ if(fabs(angle(xx.getFrst(), yy.getFrst())) > 1e-8){
+ std::cout << "This should not happen: not aligned "<< std::endl;
+ exit(1);
+ }
+
+ dmin = M_PI; jmin = 0;
+ for(int j=0 ; j<4; j++){
+ if((d = angle(xx.get(j).scnd, yy.scnd)) < dmin) {
+ jmin = j;
+ dmin = d;
+ }
+ }
+ xx = xx.get(jmin);
+ // xx.scnd and yy.scnd now form the smallest angle among the 4^2 possible.
+
+ th2 = dmin;
+ if(th2 > M_PI/4.){
+ std::cout << "This should not happen: th2 = " << th2 << std::endl;
+ exit(1);
+ }
+ if(th2 > 1e-8)
+ axis = crossprod(xx.scnd, yy.scnd).unit();
+ else {
+ axis = SVector3(1,0,0);
+ th2 = 0.;
+ }
+
+ Qtn R2 = Qtn(axis, th2);
+ return R2;
+}
+
Matrix convert(const cross3D x){
Matrix m;
SVector3 v;
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 15fe08aac4fb66c40a4465ead68ed3bdc31007a8..b453e9a40d59b61211e32431a790aaa9a4705011 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -25,7 +25,7 @@
Frame_field::Frame_field(){}
void Frame_field::init_region(GRegion* gr){
- // Fill in the ANN tree with the bondary cross field of region gr
+ // Fill in a ANN tree with the bondary cross field of region gr
#if defined(HAVE_ANN)
int index;
MVertex* vertex;
@@ -63,6 +63,8 @@ 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>
+ // for each vertex of the face gf.
unsigned int i;
int j;
bool ok;
@@ -104,20 +106,20 @@ void Frame_field::init_face(GFace* gf){
}
if(ok){
- v1.normalize();
- v2.normalize();
- v3 = crossprod(v1,v2);
- v3.normalize();
- m.set_m11(v1.x());
- m.set_m21(v1.y());
- m.set_m31(v1.z());
- m.set_m12(v2.x());
- m.set_m22(v2.y());
- m.set_m32(v2.z());
- m.set_m13(v3.x());
- m.set_m23(v3.y());
- m.set_m33(v3.z());
- temp.insert(std::pair<MVertex*,Matrix>(vertex,m));
+ v1.normalize();
+ v2.normalize();
+ v3 = crossprod(v1,v2);
+ v3.normalize();
+ m.set_m11(v1.x());
+ m.set_m21(v1.y());
+ m.set_m31(v1.z());
+ m.set_m12(v2.x());
+ m.set_m22(v2.y());
+ m.set_m32(v2.z());
+ m.set_m13(v3.x());
+ m.set_m23(v3.y());
+ m.set_m33(v3.z());
+ temp.insert(std::pair<MVertex*,Matrix>(vertex,m));
}
}
}
@@ -431,40 +433,278 @@ void Frame_field::clear(){
delete kd_tree;
annClose();
#endif
+#if defined(HAVE_ANN)
+ if(annTreeData) delete annTreeData;
+ if(annTree) delete annTree;
+#endif
}
+// BARYCENTRIC
+
#include "cross3D.h"
-// CWTSSpaceVector<> convert(const SVector3 v){
-// return CWTSSpaceVector<> (v.x(), v.y(), v.z());
-// }
+//double max(const double a, const double b) { return (b>a)?b:a;}
+double min(const double a, const double b) { return (b<a)?b:a; }
+double squ(const double a) { return a*a; }
+
+int Frame_field::build_vertex_to_vertices(GEntity* gr, int onWhat, bool initialize){
+ // build vertex to vertices data
+ std::set<MVertex*> vertices;
+ if(initialize) vertex_to_vertices.clear();
+ for(unsigned int i=0; i<gr->getNumMeshElements(); i++){
+ MElement* pElem = gr->getMeshElement(i);
+ unsigned int n = pElem->getNumVertices();
+ for(unsigned int j=0; j<n; j++){
+ MVertex* pVertex = pElem->getVertex(j);
+ if(pVertex->onWhat()->dim() != onWhat) continue;
+
+ std::map<MVertex*,std::set<MVertex*> >::iterator it = vertex_to_vertices.find(pVertex);
+ if(it != vertex_to_vertices.end()){
+ for(unsigned int k=1; k<n; k++)
+ it->second.insert(pElem->getVertex((j+k) % n));
+ }
+ else{
+ vertices.clear();
+ for(unsigned int k=1; k<n; k++)
+ vertices.insert(pElem->getVertex((j+k) % n));
+ vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(pVertex,vertices));
+ }
+
+ }
+ }
+ return vertex_to_vertices.size();
+}
-// SVector3 convert(const CWTSSpaceVector<> x){
-// return SVector3 (x[0], x[1], x[2]);
-// }
+int Frame_field::build_vertex_to_elements(GEntity* gr, bool initialize){
+ std::set<MElement*> elements;
+ if(initialize) vertex_to_elements.clear();
+ for(unsigned int i=0; i<gr->getNumMeshElements(); i++){
+ MElement* pElem = gr->getMeshElement(i);
+ unsigned int n = pElem->getNumVertices();
+ for(unsigned int j=0; j<n; j++){
+ MVertex* pVertex = pElem->getVertex(j);
+ std::map<MVertex*,std::set<MElement*> >::iterator it = vertex_to_elements.find(pVertex);
+ if(it != vertex_to_elements.end())
+ it->second.insert(pElem);
+ else{
+ elements.clear();
+ elements.insert(pElem);
+ vertex_to_elements.insert(std::pair<MVertex*,std::set<MElement*> >(pVertex,elements));
+ }
+ }
+ }
+ return vertex_to_elements.size();
+}
+
+void Frame_field::build_listVertices(GEntity* gr, int dim, bool initialize){
+ std::set<MVertex*> list;
+ for(unsigned int i=0; i<gr->getNumMeshElements(); i++){
+ MElement* pElem = gr->getMeshElement(i);
+ for(int j=0; j<pElem->getNumVertices(); j++){
+ MVertex * pVertex = pElem->getVertex(j);
+ if(pVertex->onWhat()->dim() == dim)
+ list.insert(pVertex);
+ }
+ }
+ if(initialize) listVertices.clear();
+ for(std::set<MVertex*>::const_iterator it=list.begin(); it!=list.end(); it++)
+ listVertices.push_back(*it);
+}
-// ostream& operator<< (ostream &os, SVector3 &v) {
-// os << "(" << v.x() << ", " << v.y() << ", " << v.z() << ")";
-// return os;
-// }
+int Frame_field::buildAnnData(GEntity* ge, int dim){
+ build_listVertices(ge,dim);
+ int n = listVertices.size();
+#if defined(HAVE_ANN)
+ annTreeData = annAllocPts(n,3);
+ for(int i=0; i<n; i++){
+ MVertex* pVertex = listVertices[i];
+ annTreeData[i][0] = pVertex->x();
+ annTreeData[i][1] = pVertex->y();
+ annTreeData[i][2] = pVertex->z();
+ }
+ annTree = new ANNkd_tree(annTreeData,n,3);
+#endif
+ std::cout << "ANN data for " << ge->tag() << "(" << dim << ") contains " << n << " vertices" << std::endl;
+ return n;
+}
-void Frame_field::init(GRegion* gr){
- // SVector3 Ex(1,0,0), Ey(0,1,0), Ez(0,0,1);
- // SVector3 v(1,2,0), w(0,1,1.3);
- // cross3D x(Ez, Ex);
- // cross3D y = cross3D(v, w);
+void Frame_field::deleteAnnData(){
+#if defined(HAVE_ANN)
+ if(annTreeData) delete annTreeData;
+ if(annTree) delete annTree;
+ annTreeData = NULL;
+ annTree = NULL;
+#endif
+}
- //Recombinator crossData;
- crossData.build_vertex_to_vertices(gr);
- std::cout << "Vertices in crossData = " << crossData.vertex_to_vertices.size() << std::endl;
- for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
- = crossData.vertex_to_vertices.begin();
- iter != crossData.vertex_to_vertices.end(); ++iter){
- MVertex* pVertex = iter->first;
- Matrix m = search(pVertex->x(), pVertex->y(), pVertex->z());
- // if(pVertex->onWhat()->dim() <= 2)
- crossField[pVertex->getNum()] = m;
+int Frame_field::findAnnIndex(SPoint3 p){
+ int index = 0;
+#if defined(HAVE_ANN)
+ ANNpoint query = annAllocPt(3);
+ ANNidxArray indices = new ANNidx[1];
+ ANNdistArray distances = new ANNdist[1];
+ query[0] = p.x();
+ query[1] = p.y();
+ query[2] = p.z();
+ double e = 0.;
+ annTree->annkSearch(query,1,indices,distances,e);
+ annDeallocPt(query);
+ index = indices[0];
+ delete[] indices;
+ delete[] distances;
+#endif
+ return index;
+}
+
+void Frame_field::initFace(GFace* gf){
+ // align crosses of gf with the normal (average on neighbour elements)
+ // at all vertices of the GFace gf
+ build_vertex_to_elements(gf);
+ for(std::map<MVertex*, std::set<MElement*> >::const_iterator it
+ = vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
+ std::set<MElement*> elements = it->second;
+ SVector3 Area = SVector3(0,0,0);
+ for(std::set<MElement*>::const_iterator iter=elements.begin();
+ iter != elements.end(); iter++){
+ MElement *pElem = *iter;
+ int n = pElem->getNumVertices();
+ int i;
+ for(i=0; i<n; i++){
+ if(pElem->getVertex(i) == it->first) break;
+ if(i == n-1) {
+ std::cout << "This should not happen" << std:: endl; exit(1);
+ }
+ }
+ SVector3 V0 = pElem->getVertex(i)->point();
+ SVector3 V1 = pElem->getVertex((i+n-1)%n)->point(); // previous vertex
+ SVector3 V2 = pElem->getVertex((i+1)%n)->point(); // next vertex
+ 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);
+ if(iter == crossField.end())
+ crossField.insert(std::pair<MVertex*, Matrix>(it->first,m));
+ else
+ crossField[it->first] = m;
+ }
+
+ std::cout << "Nodes in face = " << vertex_to_elements.size() << std::endl;
+
+ // 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++)
+ 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
+ = vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
+ MVertex* pVertex = it->first;
+ std::set<MElement*> elements = it->second;
+ if(elements.size() != 2){
+ std::cout << "The face has an open boundary" << std::endl;
+ exit(1);
+ }
+ MEdge edg1 = (*elements.begin())->getEdge(0);
+ MEdge edg2 = (*elements.rbegin())->getEdge(0);
+ SVector3 tangent = edg1.scaledTangent() + edg2.scaledTangent();
+ tangent.normalize();
+
+ std::map<MVertex*, Matrix>::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));
+ crossField[pVertex] = m;
+ }
+
+ // fills ANN data with the vertices of the contour (dim=1) of "gf"
+ buildAnnData(gf,1);
+
+ std::cout << "Nodes on contour = " << vertex_to_elements.size() << std::endl;
+ std::cout << "crossField = " << crossField.size() << std::endl;
+ std::cout << "region = " << gf->getNumMeshVertices() << std::endl;
+
+ // align crosses.scnd with the nearest cross of the contour
+ // fixme: we have to rebuild the vertex_to_elements list
+ // to have a working iterator oin the verticies of gf
+ // The natural iterator gf->getMeshVertex(i) seems not to work
+ // when the option PackingOfParallelogram is on.
+ build_vertex_to_elements(gf);
+ for(std::map<MVertex*, std::set<MElement*> >::const_iterator it
+ = vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
+ //for(unsigned int i=0; i<gf->getNumMeshVertices(); i++){
+ //MVertex* pVertex0 = gf->getMeshVertex(i);
+
+ MVertex* pVertex0 = it->first;
+ if(pVertex0->onWhat()->dim() != 2) continue;
+
+ std::map<MVertex*, Matrix>::iterator iter = crossField.find(pVertex0);
+ cross3D y;
+ if(iter == crossField.end()){
+ std::cout << "This should not happen: cross not found 2" << std::endl;
+ std::cout << pVertex0->x() << "," << pVertex0->y() << "," << pVertex0->z() << std::endl;
+ //exit(1);
+ y = cross3D();
+ }
+ else
+ y = cross3D((*iter).second); //local cross y.getFirst() is the normal
+
+ // Find the index of the nearest cross on the contour, using annTree
+ int index = findAnnIndex(pVertex0->point());
+ MVertex* pVertex = listVertices[index]; //nearest vertex on contour
+ iter = crossField.find(pVertex);
+ if(iter == crossField.end()){
+ std::cout << "This should not happen: cross not found 3" << std::endl;
+ exit(1);
+ }
+ cross3D x = cross3D((*iter).second); //nearest cross on contour, x.getFrst() is the normal
+
+ //Matrix m = convert(cross3D(x.getFrst(),y.getScnd()));
+ SVector3 v1 = y.getFrst();
+ Matrix 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;
+ exit(1);
+ }
+ crossField[pVertex0] = m;
}
+ checkAnnData(gf, "zzz.pos");
+ deleteAnnData();
+}
+
+void Frame_field::initRegion(GRegion* gr, int n){
+ std::list<GFace*> faces = gr->faces();
+ for(std::list<GFace*>::const_iterator iter=faces.begin(); iter!=faces.end(); iter++){
+ initFace(*iter);
+ // smoothing must be done immediately because crosses on the contour vertices
+ // are now initialized with the normal to THIS face.
+ smoothFace(*iter, n);
+ }
+
+ // Fills ANN data with the vertices of the surface (dim=2) of "gr"
+ buildAnnData(gr,2);
+ for(unsigned int i=0; i<gr->getNumMeshVertices(); i++){
+ MVertex* pVertex0 = gr->getMeshVertex(i);
+ if(pVertex0->onWhat()->dim() != 3) continue;
+ // 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));
+ }
+ deleteAnnData();
+ buildAnnData(gr,3);
+}
+
+Matrix 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){
@@ -486,9 +726,10 @@ 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<int, Matrix>::const_iterator itB = crossField.find(pVertex->getNum());
- if(itB == crossField.end()) // not found
- m = search(pVertex->x(), pVertex->y(), pVertex->z());
+ std::map<MVertex*, Matrix>::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);
+ }
else
m = itB->second;
cross3D y(m);
@@ -501,22 +742,15 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
crossDist[edge] = theta;
if (fabs(theta) > 1e-10) {
axis = eulerAxisFromQtn(Rxy); // undefined if theta==0
- // dT = convert(axis) * (theta / edge.length() / edge.length()) ;
dT = axis * (theta / edge.length() / edge.length()) ;
T += dT;
}
temp += 1. / edge.length() / edge.length();
- //std::cout << " " << pVertex->getNum() << " : " << theta << dT << std::endl;
}
- if(list.size()) T *= 1.0/temp; // average rotation vector
- double a = T.norm();
- SVector3 b = T; b.normalize();
- if(debug) std::cout << "energy= " << energy << " T= " << a << b << std::endl;
+ if(temp) T *= 1.0/temp; // average rotation vector
- //std::cout << "xold=> " << x << std::endl;
theta = T.norm();
if(fabs(theta) > 1e-10){
- //axis = convert(T);
axis = T;
if(debug) std::cout << "-> " << pVertex0->getNum() << " : " << T << std::endl;
Qtn R(axis.unit(),theta);
@@ -526,135 +760,138 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
return energy;
}
-double Frame_field::smooth(GRegion* gr){
- // loops on crossData.vertex_to_vertices
- //
+
+double Frame_field::smoothFace(GFace *gf, int n){
+ double energy=0;
+ build_vertex_to_vertices(gf, 2);
+ build_vertex_to_vertices(gf, 0, false);
+ for(int i=0; i<n; i++){
+ energy = smooth();
+ std::cout << "energy = " << energy << std::endl;
+ }
+ return energy;
+}
+
+double Frame_field::smoothRegion(GRegion *gr, int n){
+ double energy=0;
+ build_vertex_to_vertices(gr, 3);
+ //build_vertex_to_vertices(gr, 1, false);
+ //build_vertex_to_vertices(gr, 0, false);
+ for(int i=0; i<n; i++){
+ energy = smooth();
+ std::cout << "energy = " << energy << std::endl;
+ }
+ return energy;
+}
+
+double Frame_field::smooth(){
Matrix m, m0;
double enew, eold;
- // Recombinator crossData;
- // crossData.build_vertex_to_vertices(gr);
- //std::cout << "NbVert= " << crossData.vertex_to_vertices.size() << std::endl ;
-
double energy = 0;
for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
- = crossData.vertex_to_vertices.begin();
- iter != crossData.vertex_to_vertices.end(); ++iter){
-
- MVertex* pVertex0 = iter->first;
- if(pVertex0->onWhat()->dim()>2){
- SVector3 T(0, 0, 0);
- std::map<int, Matrix>::iterator itA = crossField.find(iter->first->getNum());
- if(itA == crossField.end())
- m0 = search(pVertex0->x(), pVertex0->y(), pVertex0->z());
- else
- m0 = itA->second;
-
- m = m0;
- unsigned int NbIter = 0;
- enew = findBarycenter(iter, m); // initial energy in cell
- do{
- eold = enew;
- crossField[itA->first] = m;
- enew = findBarycenter(iter, m);
- } while ((enew < eold) && (++NbIter < 10));
- energy += eold;
+ = vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
+
+ //MVertex* pVertex0 = iter->first;
+ SVector3 T(0, 0, 0);
+ std::map<MVertex*, Matrix>::iterator itA = crossField.find(iter->first);
+ if(itA == crossField.end()){
+ std::cout << "This should not happen" << std::endl; exit(1);
}
+ else
+ m0 = itA->second;
+
+ m = m0;
+ unsigned int NbIter = 0;
+ enew = findBarycenter(iter, m); // initial energy in cell
+ do{
+ eold = enew;
+ crossField[itA->first] = m;
+ enew = findBarycenter(iter, m);
+ } while ((enew < eold) && (++NbIter < 10));
+ energy += eold;
}
return energy;
}
-void Frame_field::save(GRegion* gr, const std::string& filename){
- SPoint3 p, p1;
- MVertex* pVertex;
- std::map<MVertex*,Matrix>::iterator it;
- Matrix m;
- // Recombinator crossData;
- // crossData.build_vertex_to_vertices(gr);
-
- double k = 0.04;
+void Frame_field::save(const std::vector<std::pair<SPoint3, Matrix> > data,
+ const std::string& filename){
+ const cross3D origin(SVector3(1,0,0), SVector3(0,1,0));
+ SPoint3 p1;
+ double k = 0.1;
std::ofstream file(filename.c_str());
file << "View \"cross field\" {\n";
-
- for(std::map<MVertex*, std::set<MVertex*> >::iterator iter = crossData.vertex_to_vertices.begin();
- iter != crossData.vertex_to_vertices.end(); ++iter){
- pVertex = iter->first;
- std::map<int, Matrix>::iterator itA = crossField.find(pVertex->getNum());
- if(itA == crossField.end())
- m = search(pVertex->x(), pVertex->y(), pVertex->z());
- else
- m = itA->second;
-
- p = SPoint3(pVertex->x(),pVertex->y(),pVertex->z());
- double val1=0, val2=0;
- p1 = SPoint3(pVertex->x() + k*m.get_m11(),
- pVertex->y() + k*m.get_m21(),
- pVertex->z() + k*m.get_m31());
+ for(unsigned int i=0; i<data.size(); i++){
+ SPoint3 p = data[i].first;
+ Matrix 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(),
+ p.z() + k*m.get_m31());
print_segment(p,p1,val1,val2,file);
- p1 = SPoint3(pVertex->x() - k*m.get_m11(),
- pVertex->y() - k*m.get_m21(),
- pVertex->z() - k*m.get_m31());
+ p1 = SPoint3(p.x() - k*m.get_m11(),
+ p.y() - k*m.get_m21(),
+ p.z() - k*m.get_m31());
print_segment(p,p1,val1,val2,file);
- p1 = SPoint3(pVertex->x() + k*m.get_m12(),
- pVertex->y() + k*m.get_m22(),
- pVertex->z() + k*m.get_m32());
+ p1 = SPoint3(p.x() + k*m.get_m12(),
+ p.y() + k*m.get_m22(),
+ p.z() + k*m.get_m32());
print_segment(p,p1,val1,val2,file);
- p1 = SPoint3(pVertex->x() - k*m.get_m12(),
- pVertex->y() - k*m.get_m22(),
- pVertex->z() - k*m.get_m32());
+ p1 = SPoint3(p.x() - k*m.get_m12(),
+ p.y() - k*m.get_m22(),
+ p.z() - k*m.get_m32());
print_segment(p,p1,val1,val2,file);
- p1 = SPoint3(pVertex->x() + k*m.get_m13(),
- pVertex->y() + k*m.get_m23(),
- pVertex->z() + k*m.get_m33());
+ p1 = SPoint3(p.x() + k*m.get_m13(),
+ p.y() + k*m.get_m23(),
+ p.z() + k*m.get_m33());
print_segment(p,p1,val1,val2,file);
- p1 = SPoint3(pVertex->x() - k*m.get_m13(),
- pVertex->y() - k*m.get_m23(),
- pVertex->z() - k*m.get_m33());
+ p1 = SPoint3(p.x() - k*m.get_m13(),
+ p.y() - k*m.get_m23(),
+ p.z() - k*m.get_m33());
print_segment(p,p1,val1,val2,file);
}
file << "};\n";
file.close();
}
-void Frame_field::fillTreeVolume(GRegion* gr){
-#if defined(HAVE_ANN)
- int n = crossData.vertex_to_vertices.size();
- std::cout << "Filling ANN tree with " << n << " vertices" << std::endl;
- annTreeData = annAllocPts(n,3);
- int index=0;
- vertIndices.clear();
- 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.push_back(pVertex->getNum());
- index++;
+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();
+ it != crossField.end(); it++){
+ SPoint3 p = it->first->point();
+ Matrix m = it->second;
+ double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin))*180./M_PI, val2=val1;
+ p1 = SPoint3(p.x() + k*m.get_m11(),
+ p.y() + k*m.get_m21(),
+ 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);
+ p1 = SPoint3(p.x() + k*m.get_m12(),
+ p.y() + k*m.get_m22(),
+ p.z() + k*m.get_m32());
+ print_segment(p,p1,val1,val2,file);
+ p1 = SPoint3(p.x() - k*m.get_m12(),
+ p.y() - k*m.get_m22(),
+ p.z() - k*m.get_m32());
+ if(full) print_segment(p,p1,val1,val2,file);
+ p1 = SPoint3(p.x() + k*m.get_m13(),
+ p.y() + k*m.get_m23(),
+ p.z() + k*m.get_m33());
+ if(full) print_segment(p,p1,val1,val2,file);
+ p1 = SPoint3(p.x() - k*m.get_m13(),
+ p.y() - k*m.get_m23(),
+ p.z() - k*m.get_m33());
+ if(full) print_segment(p,p1,val1,val2,file);
}
- 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);
- std::map<int, Matrix>::const_iterator it = crossField.find(vertIndices[indices[0]]);
- //annTreeData[indices[0]] gives the coordinates
- delete[] indices;
- delete[] distances;
- return it->second;
-#else
- return Matrix();
-#endif
+ file << "};\n";
+ file.close();
}
void Frame_field::save_dist(const std::string& filename){
@@ -665,7 +902,7 @@ void Frame_field::save_dist(const std::string& filename){
it != crossDist.end(); it++){
MVertex* pVerta = it->first.getVertex(0);
MVertex* pVertb = it->first.getVertex(1);
- double value = it->second;
+ double value = it->second*180./M_PI;
if(it->first.length())
value /= it->first.length();
file << "SL ("
@@ -677,6 +914,24 @@ void Frame_field::save_dist(const std::string& filename){
file.close();
}
+void Frame_field::checkAnnData(GEntity* ge, const std::string& filename){
+#if defined(HAVE_ANN)
+ std::ofstream file(filename.c_str());
+ file << "View \"ANN pairing\" {\n";
+ for(unsigned int i=0; i<ge->getNumMeshVertices(); i++){
+ MVertex* pVerta = ge->getMeshVertex(i);
+ MVertex* pVertb = listVertices[findAnnIndex(pVerta->point())];
+ double value = pVerta->distance(pVertb);
+ file << "SL ("
+ << pVerta->x() << ", " << pVerta->y() << ", " << pVerta->z() << ", "
+ << pVertb->x() << ", " << pVertb->y() << ", " << pVertb->z() << ")"
+ << "{" << value << "," << value << "};\n";
+ }
+ file << "};\n";
+ file.close();
+#endif
+}
+
#include "PView.h"
#include "PViewDataList.h"
void Frame_field::save_energy(GRegion* gr, const std::string& filename){
@@ -1310,16 +1565,16 @@ void Nearest_point::clear(){
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<MVertex*, Matrix> Frame_field::crossField;
std::map<MEdge, double, Less_Edge> Frame_field::crossDist;
-Recombinator Frame_field::crossData;
-
+std::map<MVertex*,std::set<MVertex*> > Frame_field::vertex_to_vertices;
+std::map<MVertex*,std::set<MElement*> > Frame_field::vertex_to_elements;
+std::vector<MVertex*> Frame_field::listVertices;
#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 888ceafcbb258e719a75c213f1464c4e73dcdb34..6316bf4fdb256b354f322379bac9a5a9f2e71810 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -27,15 +27,14 @@ class Frame_field{
private:
static std::map<MVertex*, Matrix> temp;
static std::vector<std::pair<MVertex*, Matrix> > field;
- static std::map<int, Matrix> crossField;
- //static std::map<MVertex*, Matrix, MVertexLessThanNum> crossField;
+ static std::map<MVertex*, Matrix> crossField;
static std::map<MEdge, double, Less_Edge> crossDist;
+ static std::vector<MVertex*> listVertices;
#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:
@@ -50,15 +49,27 @@ class Frame_field{
static void print_field1();
static void print_field2(GRegion*);
static void print_segment(SPoint3,SPoint3,double,double,std::ofstream&);
- static Recombinator crossData;
- static void init(GRegion* gr);
- static double smooth(GRegion* gr);
+ static std::map<MVertex*,std::set<MVertex*> > vertex_to_vertices;
+ static std::map<MVertex*,std::set<MElement*> > vertex_to_elements;
+ static int build_vertex_to_vertices(GEntity* gr, int onWhat, bool initialize=true);
+ static int build_vertex_to_elements(GEntity* gr, bool initialize=true);
+ static void build_listVertices(GEntity* gr, int dim, bool initialize=true);
+ 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 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 fillTreeVolume(GRegion* gr);
- static Matrix findNearestCross(double x,double y,double z);
- static void save(GRegion* gr, const std::string &filename);
+ static void save(const std::vector<std::pair<SPoint3, Matrix> > data,
+ const std::string& filename);
+ static void saveCrossField(const std::string& filename, double scale, bool full=true);
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);
static GRegion* test();
static void clear();
};
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index 87e81067654ce30a4c74aed3eb5bbfb12e118dce..838ba5fe278393b0e2a5ff0195b0c1c7b91bdb38 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -325,7 +325,19 @@ void Filler::treat_model(){
}
void Filler::treat_region(GRegion* gr){
- #if defined(HAVE_RTREE)
+
+ int NumSmooth = 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);
+ }
+
+#if defined(HAVE_RTREE)
unsigned int i;
int j;
int count;
@@ -347,23 +359,6 @@ void Filler::treat_region(GRegion* gr){
RTree<Node*,double,3,double> rtree;
Frame_field::init_region(gr);
-
- int NumSmooth = CTX::instance()->mesh.smoothCrossField;
- if(NumSmooth){
- Frame_field::init(gr);
- Frame_field::save(gr,"cross_init.pos");
- double enew = Frame_field::smooth(gr);
- int NumIter = 0;
- double eold = 0;
- do{
- std::cout << "Smoothing: energy(" << NumIter << ") = " << enew << std::endl;
- eold = enew;
- enew = Frame_field::smooth(gr);
- } while((eold > enew) && (NumIter++ < NumSmooth));
- Frame_field::save(gr,"cross_smooth.pos");
- Frame_field::fillTreeVolume(gr);
- }
-
Size_field::init_region(gr);
Size_field::solve(gr);
octree = new MElementOctree(gr->model());
@@ -374,16 +369,16 @@ void Filler::treat_region(GRegion* gr){
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- temp.insert(vertex);
- }
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ temp.insert(vertex);
+ }
}
for(it=temp.begin();it!=temp.end();it++){
- if((*it)->onWhat()->dim()<3){
- boundary_vertices.push_back(*it);
- }
+ if((*it)->onWhat()->dim()<3){
+ boundary_vertices.push_back(*it);
+ }
}
//std::ofstream file("nodes.pos");
//file << "View \"test\" {\n";
@@ -393,11 +388,11 @@ void Filler::treat_region(GRegion* gr){
y = boundary_vertices[i]->y();
z = boundary_vertices[i]->z();
- node = new Node(SPoint3(x,y,z));
- compute_parameters(node,gr);
- rtree.Insert(node->min,node->max,node);
- fifo.push(node);
- //print_node(node,file);
+ node = new Node(SPoint3(x,y,z));
+ compute_parameters(node,gr);
+ rtree.Insert(node->min,node->max,node);
+ fifo.push(node);
+ //print_node(node,file);
}
count = 1;
@@ -469,16 +464,17 @@ void Filler::treat_region(GRegion* gr){
delete octree;
Size_field::clear();
Frame_field::clear();
- #endif
+#endif
}
Metric Filler::get_metric(double x,double y,double z){
Metric m;
Matrix m2;
- if(!CTX::instance()->mesh.smoothCrossField)
- m2 = Frame_field::search(x,y,z);
+ if(CTX::instance()->mesh.smoothCrossField){
+ m2 = Frame_field::findCross(x,y,z);
+ }
else
- m2 = Frame_field::findNearestCross(x,y,z);
+ m2 = Frame_field::search(x,y,z);
m.set_m11(m2.get_m11());
m.set_m21(m2.get_m21());
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 239a4db3eba8a0d1be6178b2a706c65a7d336993..598fe227298e94bade2bba11bb9c67380a29ca7d 100755
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -1887,26 +1887,26 @@ void Recombinator::build_vertex_to_vertices(GRegion* gr){
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
+ for(j=0;j<element->getNumVertices();j++){
a = element->getVertex(j);
- b = element->getVertex((j+1)%4);
- c = element->getVertex((j+2)%4);
- d = element->getVertex((j+3)%4);
-
- it = vertex_to_vertices.find(a);
- if(it!=vertex_to_vertices.end()){
- it->second.insert(b);
- it->second.insert(c);
- it->second.insert(d);
- }
- else{
- bin.clear();
- bin.insert(b);
- bin.insert(c);
- bin.insert(d);
- vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(a,bin));
- }
- }
+ b = element->getVertex((j+1)%4);
+ c = element->getVertex((j+2)%4);
+ d = element->getVertex((j+3)%4);
+
+ it = vertex_to_vertices.find(a);
+ if(it!=vertex_to_vertices.end()){
+ it->second.insert(b);
+ it->second.insert(c);
+ it->second.insert(d);
+ }
+ else{
+ bin.clear();
+ bin.insert(b);
+ bin.insert(c);
+ bin.insert(d);
+ vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(a,bin));
+ }
+ }
}
}