diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index 8a83e0154087c925c05372182eeef07af3218c14..4046dadca29841b98aeba1f7872470dcac9ab0ea 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -61,6 +61,7 @@ typedef unsigned long intptr_t;
#include "HighOrder.h"
#include "OS.h"
#include "onelabUtils.h"
+#include "gmshCrossFields.h"
#if defined(HAVE_3M)
#include "3M.h"
#endif
@@ -2117,6 +2118,12 @@ static void mesh_optimize_cb(Fl_Widget *w, void *data)
drawContext::global()->draw();
}
+static void mesh_cross_compute_cb(Fl_Widget *w, void *data)
+{
+ computeCrossField (GModel::current());
+ drawContext::global()->draw();
+}
+
static void mesh_refine_cb(Fl_Widget *w, void *data)
{
RefineMesh(GModel::current(), CTX::instance()->mesh.secondOrderLinear);
@@ -4298,6 +4305,8 @@ static menuItem static_modules[] = {
(Fl_Callback *)mesh_inspect_cb} ,
{"0Modules/Mesh/Refine by splitting",
(Fl_Callback *)mesh_refine_cb} ,
+ {"0Modules/Mesh/Compute Cross Field",
+ (Fl_Callback *)mesh_cross_compute_cb} ,
#if defined(HAVE_METIS)
{"0Modules/Mesh/Partition",
(Fl_Callback *)mesh_partition_cb} ,
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index eba8bb620e78bc59488d80d83be99b0bc8fa6aad..e1ceb7f1b460c4300a13f4b6f67709421f0fa9e8 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2951,7 +2951,6 @@ static void recurClassifyEdges(MTri3 *t, std::map<MTriangle *, GFace *> &reverse
while (!_stack.empty()){
t = _stack.top();
_stack.pop();
- // printf("%d\n",rec);
if(!t->isDeleted()) {
trisTouched.erase(t);
t->setDeleted(true);
@@ -3123,8 +3122,8 @@ void GModel::classifyFaces()
for(std::map<std::pair<int, int>, GEdge *>::iterator ite = newEdges.begin();
ite != newEdges.end(); ++ite) {
- // printf("%d %d -->
- // %d\n",ite->first.first,ite->first.second,ite->second->tag());
+
+ printf("%d %d --> %d\n",ite->first.first,ite->first.second,ite->second->tag());
std::list<MLine *> allSegments;
for(std::size_t i = 0; i < ite->second->lines.size(); i++)
allSegments.push_back(ite->second->lines[i]);
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 5b24751f858738fde1212d24ebc78e484f125d78..e35ef4c5995df9766d6a86cf4a626b04fd92c3e2 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -15,7 +15,7 @@ set(SRC
meshGRegionBoundaryRecovery.cpp
meshGRegionHxt.cpp meshGRegionNetgen.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
- meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
+ meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
meshGRegionLocalMeshMod.cpp meshGRegionMMG3D.cpp
meshGRegionBoundaryLayer.cpp
meshRelocateVertex.cpp
@@ -31,6 +31,7 @@ set(SRC
DivideAndConquer.cpp
Field.cpp
filterElements.cpp
+ gmshCrossFields.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Mesh/gmshCrossFields.cpp b/Mesh/gmshCrossFields.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d260ec3606f17dd720b02f6d7c5454d1024aecd0
--- /dev/null
+++ b/Mesh/gmshCrossFields.cpp
@@ -0,0 +1,248 @@
+#include <vector>
+#include "gmshCrossFields.h"
+#include "GModel.h"
+#include "GFace.h"
+#include "MEdge.h"
+#include "MTriangle.h"
+#include "GmshMessage.h"
+#include "Context.h"
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
+class cross2d {
+public:
+ MEdge _e;
+ std::vector<MEdge> _neighbors;
+ std::vector<cross2d*> _cneighbors;
+ // euler angles
+ double _a,_b,_c;
+ double _atemp,_btemp,_ctemp;
+ std::vector<MTriangle *> _t;
+ cross2d (MEdge &e, MTriangle *r, MEdge &e1, MEdge &e2) : _e(e) , _a (0), _b(0), _c(0)
+ {
+ _t.push_back(r);
+ _neighbors.push_back(e1);
+ _neighbors.push_back(e2);
+ }
+ void normalize(double &a){
+ double D = M_PI*.5;
+ if (a < 0 ) while (a < 0) a += D;
+ if (a >= D) while (a >= D) a -= D;
+ }
+ void finish ( std::map<MEdge,cross2d,Less_Edge > &C){
+ for (size_t i=0;i<_neighbors.size();i++){
+ std::map<MEdge,cross2d,Less_Edge >::iterator it = C.find(_neighbors[i]);
+ if (it == C.end())Msg::Error("impossible situation");
+ else _cneighbors.push_back(&(it->second));
+ }
+ if (_cneighbors.size() != 4){
+ _a = atan2 (_e.getVertex(1)->y()-_e.getVertex(0)->y(),_e.getVertex(1)->x()-_e.getVertex(0)->x());
+ normalize(_a);
+ _atemp = _a;
+ }
+ _neighbors.clear();
+ _b = _btemp = sin (4*_a);
+ _c = _ctemp = cos (4*_a);
+ }
+ double average_init (){
+ if (_cneighbors.size() == 4){
+
+ _btemp = 0.25*(_cneighbors[0]->_b+_cneighbors[1]->_b+_cneighbors[2]->_b+_cneighbors[3]->_b);
+ _ctemp = 0.25*(_cneighbors[0]->_c+_cneighbors[1]->_c+_cneighbors[2]->_c+_cneighbors[3]->_c);
+ _atemp = 0.25*atan2(_b,_c);
+ normalize(_atemp);
+ return 1;
+ }
+ return 1;
+ }
+
+ void update () {
+ _a = _atemp;
+ _b = _btemp;
+ _c = _ctemp;
+ }
+
+ double grad () {
+ if (_cneighbors.size() == 4){
+ double D = M_PI*.5;
+ double a[4] = {_cneighbors[0]->_a,_cneighbors[1]->_a,
+ _cneighbors[2]->_a,_cneighbors[3]->_a};
+ double b[4];
+ for (int i=0;i<4;i++){
+ if (fabs(_a - a[i]) < fabs(_a - (a[i]+D)) &&
+ fabs(_a - a[i]) < fabs(_a - (a[i]-D))) {b[i] = a[i];}
+ else if (fabs(_a - (a[i]+D)) < fabs(_a - (a[i])) &&
+ fabs(_a - (a[i]+D)) < fabs(_a - (a[i]-D))) {b[i] = a[i]+D;}
+ else {b[i] = a[i]-D;}
+ }
+ return fabs(_a-b[0]) + fabs(_a-b[1]) + fabs(_a-b[2]) + fabs(_a-b[3]) ;
+ }
+ return 0;
+ }
+
+ double average (){
+ if (_cneighbors.size() == 4){
+ double D = M_PI*.5;
+ double a[4] = {_cneighbors[0]->_a,_cneighbors[1]->_a,
+ _cneighbors[2]->_a,_cneighbors[3]->_a};
+ double b[4];
+ double avg = 0.0;
+ for (int i=0;i<4;i++){
+ if (fabs(_a - a[i]) < fabs(_a - (a[i]+D)) &&
+ fabs(_a - a[i]) < fabs(_a - (a[i]-D))) {b[i] = a[i];}
+ else if (fabs(_a - (a[i]+D)) < fabs(_a - (a[i])) &&
+ fabs(_a - (a[i]+D)) < fabs(_a - (a[i]-D))) {b[i] = a[i]+D;}
+ else {b[i] = a[i]-D;}
+ }
+ avg = 0.25*(b[0]+b[1]+b[2]+b[3]);
+
+ normalize(avg);
+
+ double d = fabs(_a-avg);
+ _atemp = avg;
+ return d;
+ }
+ return 0;
+ }
+};
+// ---------------------------------------------
+// TODO : MAKE IT PARALLEL AND SUPERFAST
+// DO IT ON SURFACES
+// ---------------------------------------------
+
+int computeCrossField2dTheta (GModel *gm,std::vector<GFace *> &f, const char * outputName){
+ Msg::SetNumThreads(Msg::GetMaxThreads());
+
+ std::map<MEdge,cross2d,Less_Edge > C;
+ for (size_t i=0;i<f.size();i++){
+ for (size_t j=0;j<f[i]->triangles.size();j++){
+ MTriangle *t = f[i]->triangles[j];
+ for (size_t k=0;k<3;k++){
+ MEdge e = t->getEdge(k);
+ MEdge e1 = t->getEdge((k+1)%3);
+ MEdge e2 = t->getEdge((k+2)%3);
+ cross2d c (e,t,e1,e2);
+ std::map<MEdge,cross2d,Less_Edge >::iterator it = C.find(e);
+ if (it == C.end())C.insert(std::make_pair(e,c));
+ else {
+ it->second._t.push_back(t);
+ it->second._neighbors.push_back(e1);
+ it->second._neighbors.push_back(e2);
+ }
+ }
+ }
+ }
+
+ std::map<MEdge,cross2d,Less_Edge >::iterator it = C.begin();
+ for (; it !=C.end();++it)it->second.finish(C);
+
+ std::vector<cross2d*> pc;
+ for (it = C.begin(); it !=C.end();++it)pc.push_back(&(it->second));
+
+ const int MAXITER = 20000;
+ int ITER = 0;
+ while (ITER++ < MAXITER){
+ double DELTA = 0;
+#if defined(_OPENMP)
+#pragma omp parallel for schedule(dynamic)
+#endif
+ for (size_t i=0;i<pc.size(); i++){
+ if (ITER<2000)DELTA += pc[i]->average_init();
+ else DELTA += pc[i]->average();
+ }
+
+#if defined(_OPENMP)
+#pragma omp parallel for schedule(dynamic)
+#endif
+ for (size_t i=0;i<pc.size(); i++)pc[i]->update();
+
+ if (ITER%100 == 0)printf("DELTA = %12.5E\n",DELTA);
+ if (DELTA < 1.e-12)break;
+ // getchar();
+ }
+
+ FILE *of = fopen ("cross.pos","w");
+ fprintf(of,"View \"\"{\n");
+ for (it = C.begin(); it !=C.end();++it){
+ double a = it->second.grad();
+ MEdge e0 = it->second._e;
+ fprintf(of,"SP(%g,%g,%g){%g};",
+ 0.5*(e0.getVertex(0)->x()+e0.getVertex(1)->x()),
+ 0.5*(e0.getVertex(0)->y()+e0.getVertex(1)->y()),
+ 0.5*(e0.getVertex(0)->z()+e0.getVertex(1)->z()), a);
+ }
+
+ for (size_t i=0;i<f.size();i++){
+ for (size_t j=0;j<f[i]->triangles.size();j++){
+ MTriangle *t = f[i]->triangles[j];
+ MEdge e0 = t->getEdge(0);
+ MEdge e1 = t->getEdge(1);
+ MEdge e2 = t->getEdge(2);
+ std::map<MEdge,cross2d,Less_Edge >::iterator it0 = C.find(e0);
+ std::map<MEdge,cross2d,Less_Edge >::iterator it1 = C.find(e1);
+ std::map<MEdge,cross2d,Less_Edge >::iterator it2 = C.find(e2);
+ double a0 = it0->second._a;
+ double a1 = it1->second._a;
+ double a2 = it2->second._a;
+ // double avg = atan2 (sin(a1)+sin(a2)+sin(a0),
+ // cos(a1)+cos(a2)+cos(a0));
+ fprintf(of,"VP(%g,%g,%g){%g,%g,%g};",
+ 0.5*(e0.getVertex(0)->x()+e0.getVertex(1)->x()),
+ 0.5*(e0.getVertex(0)->y()+e0.getVertex(1)->y()),
+ 0.5*(e0.getVertex(0)->z()+e0.getVertex(1)->z()),
+ cos (a0), sin(a0),0.0);
+ fprintf(of,"VP(%g,%g,%g){%g,%g,%g};",
+ 0.5*(e1.getVertex(0)->x()+e1.getVertex(1)->x()),
+ 0.5*(e1.getVertex(0)->y()+e1.getVertex(1)->y()),
+ 0.5*(e1.getVertex(0)->z()+e1.getVertex(1)->z()),
+ cos (a1), sin(a1),0.0);
+ fprintf(of,"VP(%g,%g,%g){%g,%g,%g};",
+ 0.5*(e2.getVertex(0)->x()+e2.getVertex(1)->x()),
+ 0.5*(e2.getVertex(0)->y()+e2.getVertex(1)->y()),
+ 0.5*(e2.getVertex(0)->z()+e2.getVertex(1)->z()),
+ cos (a2), sin(a2),0.0);
+ fprintf(of,"VP(%g,%g,%g){%g,%g,%g};",
+ 0.5*(e0.getVertex(0)->x()+e0.getVertex(1)->x()),
+ 0.5*(e0.getVertex(0)->y()+e0.getVertex(1)->y()),
+ 0.5*(e0.getVertex(0)->z()+e0.getVertex(1)->z()),
+ -sin (a0), cos(a0),0.0);
+ fprintf(of,"VP(%g,%g,%g){%g,%g,%g};",
+ 0.5*(e1.getVertex(0)->x()+e1.getVertex(1)->x()),
+ 0.5*(e1.getVertex(0)->y()+e1.getVertex(1)->y()),
+ 0.5*(e1.getVertex(0)->z()+e1.getVertex(1)->z()),
+ -sin (a1), cos(a1),0.0);
+ fprintf(of,"VP(%g,%g,%g){%g,%g,%g};",
+ 0.5*(e2.getVertex(0)->x()+e2.getVertex(1)->x()),
+ 0.5*(e2.getVertex(0)->y()+e2.getVertex(1)->y()),
+ 0.5*(e2.getVertex(0)->z()+e2.getVertex(1)->z()),
+ -sin (a2), cos(a2),0.0);
+
+ /*
+ fprintf(of,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};",
+ t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z(),
+ t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z(),
+ t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z(),
+ cos(avg),sin(avg),0.,cos(avg),sin(avg),0.,cos(avg),sin(avg),0.);
+ */
+
+ }
+ }
+ fprintf(of,"};\n");
+ fclose(of);
+ return 0;
+
+}
+
+int computeCrossField (GModel *gm)
+{
+ std::vector<GFace *> f ;
+ for (GModel::fiter it = gm->firstFace() ; it != gm->lastFace();++it){
+ GFace *gf = *it;
+ f.push_back(gf);
+ }
+ // if (!strcmp(method,"theta"))
+ return computeCrossField2dTheta (gm,f,"toto");
+ // else
+ // Msg::Error("Unknown Cross Field Computation Method %s",method);
+}
diff --git a/Mesh/gmshCrossFields.h b/Mesh/gmshCrossFields.h
new file mode 100644
index 0000000000000000000000000000000000000000..d1d43379169a97886334f087cb93526aa8949e60
--- /dev/null
+++ b/Mesh/gmshCrossFields.h
@@ -0,0 +1,6 @@
+#ifndef _CROSS_FIELDS_H_
+#define _CROSS_FIELDS_H_
+class GModel;
+int computeCrossField (GModel*);
+//int computeCrossField2d (const char * method, const char * outputName, std::vector<int> &faces);
+#endif