diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index a68607dff5535aa5434db04d9c149ec2f102e7ff..1996015a457a3aaf1989cf8b65da1e97145c40d8 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -6,6 +6,8 @@
#include "Context.h"
#include "MVertex.h"
#include "MElement.h"
+#include "GFace.h"
+#include "DivideAndConquer.h"
#include "Bindings.h"
#include "dgSystemOfEquations.h"
#include "luaFunction.h"
@@ -234,6 +236,7 @@ binding::binding(){
}
_instance=this;
L = lua_open();
+ /*
luaopen_base(L);
luaopen_table(L);
luaopen_os(L);
@@ -241,6 +244,8 @@ binding::binding(){
luaopen_string(L);
luaopen_math(L);
luaopen_debug(L);
+ */
+ luaL_openlibs(L);
lua_register(L,"help",luaHelp);
#ifdef HAVE_READLINE
@@ -248,6 +253,11 @@ binding::binding(){
lua_register(L,"clearHistory",luaClear);
#endif
+ // lua_pushcfunction(L, luaopen_io);
+ // lua_call(L, 0, 0);
+
+
+
// Register Lua bindings
GModel::registerBindings(this);
fullMatrix<double>::registerBindings(this);
@@ -263,6 +273,9 @@ binding::binding(){
function::registerDefaultFunctions();
MVertex::registerBindings(this);
MElement::registerBindings(this);
+ DocRecord::registerBindings(this);
+ GEntity::registerBindings(this);
+ GFace::registerBindings(this);
}
binding *binding::_instance=NULL;
#endif
diff --git a/Geo/GEntity.cpp b/Geo/GEntity.cpp
index 39572e24ce6031557155ba13c27715d510f87292..17247cb3a65e769b737b0f384b870bdf717474df 100644
--- a/Geo/GEntity.cpp
+++ b/Geo/GEntity.cpp
@@ -64,3 +64,25 @@ std::string GEntity::getInfoString()
return sstream.str();
}
+
+#include "Bindings.h"
+
+void GEntity::registerBindings(binding *b)
+{
+ classBinding *cb = b->addClass<GEntity>("GEntity");
+ cb->setDescription("A GEntity is a geometrical entity of the model.");
+
+ methodBinding *cm;
+ cm = cb->addMethod("model", &GEntity::model);
+ cm->setDescription("returns the geometric model the entity belongs to.");
+
+ /*
+ cm = cb->addMethod("getNumMeshElements", (unsigned int (GEntity::*)() ) &GEntity::getNumMeshElements);
+ cm->setDescription("Return the number of elements of the mesh of the entity.");
+ cm = cb->addMethod("getMeshElement", &GEntity::getMeshElement);
+ cm->setDescription("returns the ith MElement.");
+ cm->setArgNames("i",NULL);
+ */
+
+}
+
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 90bc88fefec67dc03d4e676474f9270f1bdbc022..b24b6db0d7b9c8877ddebc087c7c75d2779c7a69 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -22,6 +22,7 @@ class GRegion;
class MVertex;
class MElement;
class VertexArray;
+class binding;
// A geometric model entity.
class GEntity {
@@ -266,6 +267,10 @@ class GEntity {
// get the mesh vertex at the given index
MVertex *getMeshVertex(unsigned int index) { return mesh_vertices[index]; }
+
+ // bindings
+ static void registerBindings(binding *b);
+
};
class GEntityLessThan {
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 803c7ef6d9682f69088593ab494c5c1e9e9331e6..7032c53b1a9990b75a0337af854d3eb8b7da3716 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -17,6 +17,7 @@
#include "Numeric.h"
#include "GaussLegendre1D.h"
#include "Context.h"
+#include "meshGFaceLloyd.h"
#define SQU(a) ((a)*(a))
@@ -1044,3 +1045,22 @@ bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
}
return true;
}
+
+void GFace::lloyd(int nbiter, int infn){
+ lloydAlgorithm algo(nbiter, infn);
+ algo(this);
+}
+
+#include "Bindings.h"
+
+void GFace::registerBindings(binding *b)
+{
+ classBinding *cb = b->addClass<GFace>("GFace");
+ cb->setParentClass<GEntity>();
+ cb->setDescription("A Geometrical Face.");
+ methodBinding *cm;
+ cm = cb->addMethod("lloyd", &GFace::lloyd);
+ cm->setDescription("do N iteration of Lloyd's algorithm using or not the infinite norm");
+ cm->setArgNames("N","infiniteNorm",NULL);
+}
+
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 768897e4a12545b4983d511084f764a16cae9be9..09a2f988e254857d5da33a5ea38c0d2a176e3b49 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -238,6 +238,9 @@ class GFace : public GEntity
bool fillPointCloud(double maxDist, std::vector<SPoint3> *points,
std::vector<SVector3> *normals=0);
+ // apply Lloyd's algorithm to the mesh
+ void lloyd (int nIter, int infNorm = 0);
+
struct {
// do we recombine the triangles of the mesh?
int recombine;
@@ -285,6 +288,16 @@ class GFace : public GEntity
std::vector<MTriangle*> triangles;
std::vector<MQuadrangle*> quadrangles;
std::vector<MPolygon*> polygons;
+
+ // an array with additional vertices that are supposed to exist
+ // in the final mesh of the model face. This can be used for
+ // boundary layer meshes or when using Lloyd-like smoothing algorithms
+ // those vertices are classifed on this GFace, their type is MFaceVertex.
+ // After mesh generation, those are moved to the mesh_vertices array
+ std::vector<MVertex*> _additional_vertices;
+
+
+ static void registerBindings(binding *b);
};
#endif
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 81c9fde60c4e1401664b3259eac110ef1951cf0d..82707ba69848a06e570ebcf336d507c267aea321 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1372,6 +1372,9 @@ void GModel::registerBindings(binding *b)
cm = cb->addMethod("getMeshVertexByTag",&GModel::getMeshVertexByTag);
cm->setDescription("access a mesh vertex by tag, using the vertex cache");
cm->setArgNames("tag",NULL);
+ cm = cb->addMethod("getFaceByTag",&GModel::getFaceByTag);
+ cm->setDescription("access a model face by its tag");
+ cm->setArgNames("tag",NULL);
cm = cb->setConstructor<GModel>();
cm->setDescription("Create an empty GModel.");
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 02c352e1e1d2adf2f0fd940b2b5a07c758d98ead..39dc4f2feddefe0b39f9f016bc534d0aff42b25b 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -818,6 +818,8 @@ void MElement::registerBindings(binding *b)
cm->setArgNames("i",NULL);
cm = cb->addMethod("getType", &MElement::getType);
cm->setDescription("get the type of the element");
+ cm = cb->addMethod("getTypeForMSH", &MElement::getTypeForMSH);
+ cm->setDescription("get the gmsh type of the element");
cm = cb->addMethod("getPartition", &MElement::getPartition);
cm->setDescription("get the partition to which the element belongs");
cm = cb->addMethod("getPolynomialOrder", &MElement::getPolynomialOrder);
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index b49f985f36202be68c20d83d5d65be364b834bcc..53395398219266715b500f2f2295a9e7b9dc008f 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -13,6 +13,7 @@ set(SRC
meshGFaceExtruded.cpp
meshGFaceBDS.cpp
meshGFaceDelaunayInsertion.cpp
+ meshGFaceLloyd.cpp
meshGFaceOptimize.cpp
meshGFaceQuadrilateralize.cpp
meshGRegion.cpp
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 818fb7eed069e6279a80b1cffce89b54c0bb52be..822e029a63c5c0d0acc42f0ae898f2835609fe1b 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -301,6 +301,11 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
++itvx;
}
+ // add _additional_vertices
+ all_vertices.insert(gf->_additional_vertices.begin(),
+ gf->_additional_vertices.end());
+
+
if(all_vertices.size() < 3){
Msg::Warning("Mesh Generation of Model Face %d Skipped: "
"Only %d Mesh Vertices on The Contours",
@@ -478,7 +483,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
if(RECUR_ITER > 0)
Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations",
- RECUR_ITER);
+ RECUR_ITER);
Msg::Debug("Boundary Edges recovered for surface %d", gf->tag());
@@ -1234,7 +1239,7 @@ void deMeshGFace::operator() (GFace *gf)
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
}
-int debugSurface = -1;
+int debugSurface = -100;
void meshGFace::operator() (GFace *gf)
{
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 3f79f9ef02dba2ffaf6c5eb771bf56995c184ad5..bd4ce25304e9f7c0a8f489716f6e482e6546e66f 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -19,7 +19,6 @@
#include "DivideAndConquer.h"
#include "Numeric.h"
#include "robustPredicates.h"
-#include "MallocUtils.h"
#define Pred(x) ((x)->prev)
#define Succ(x) ((x)->next)
@@ -334,7 +333,7 @@ int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPo
double alpha1, alpha, beta, xx, yy;
int first;
- newp = (DListRecord *) Malloc(sizeof(DListRecord));
+ newp = new DListRecord;
newp->point_num = newPoint;
if(*dlist == NULL) {
@@ -406,7 +405,7 @@ int DocRecord::DListDelete(DListPeek *dlist, PointNumero oldPoint)
return 0;
if(Succ(*dlist) == *dlist) {
if((*dlist)->point_num == oldPoint) {
- Free(*dlist);
+ delete *dlist;
*dlist = NULL;
return 1;
}
@@ -421,7 +420,7 @@ int DocRecord::DListDelete(DListPeek *dlist, PointNumero oldPoint)
if(p == *dlist) {
*dlist = Succ(p);
}
- Free(p);
+ delete p;
return 1;
}
p = Succ(p);
@@ -440,10 +439,10 @@ int DocRecord::Delete(PointNumero a, PointNumero b)
}
// compte les points sur le polygone convexe
-int DocRecord::CountPointsOnHull(int n)
+int DocRecord::CountPointsOnHull()
{
PointNumero p, p2, temp;
- int i;
+ int i, n=numPoints;
if(points[0].adjacent == NULL)
return 0;
@@ -459,6 +458,26 @@ int DocRecord::CountPointsOnHull(int n)
return (i <= n) ? i : -1;
}
+// compte les points sur le polygone convexe
+void DocRecord::ConvexHull()
+{
+ PointNumero p, p2, temp;
+
+ if(points[0].adjacent == NULL)
+ return;
+ int count = 0;
+ p = 0;
+ _hull[count++] = p;
+ p2 = First(0);
+ while((p2 != 0) && (count < numPoints)) {
+ temp = p2;
+ p2 = Successor(p2, p);
+ p = temp;
+ _hull[count++] = p;
+ }
+}
+
+
PointNumero *DocRecord::ConvertDlistToArray(DListPeek *dlist, int *n)
{
DListPeek p, temp;
@@ -470,7 +489,7 @@ PointNumero *DocRecord::ConvertDlistToArray(DListPeek *dlist, int *n)
max++;
p = Pred(p);
} while(p != *dlist);
- ptr = (PointNumero *) Malloc((max + 1) * sizeof(PointNumero));
+ ptr = new PointNumero[max + 1];
if(ptr == NULL)
return NULL;
p = *dlist;
@@ -478,7 +497,7 @@ PointNumero *DocRecord::ConvertDlistToArray(DListPeek *dlist, int *n)
ptr[i] = p->point_num;
temp = p;
p = Pred(p);
- Free(temp);
+ delete temp;
}
ptr[max] = ptr[0];
*dlist = NULL;
@@ -486,6 +505,185 @@ PointNumero *DocRecord::ConvertDlistToArray(DListPeek *dlist, int *n)
return ptr;
}
+// build ready to use voronoi datas
+void DocRecord::ConvertDListToVoronoiData()
+{
+ if (_adjacencies){
+ for(int i = 0; i < numPoints; i++)
+ if (_adjacencies[i].t)
+ delete [] _adjacencies[i].t;
+ delete [] _adjacencies;
+ }
+ if (_hull)delete [] _hull;
+ _hullSize = CountPointsOnHull ();
+ _hull = new PointNumero[_hullSize];
+ ConvexHull();
+ std::sort(_hull, _hull+_hullSize);
+ _adjacencies = new STriangle[numPoints];
+
+ for(PointNumero i = 0; i < numPoints; i++)
+ _adjacencies[i].t = ConvertDlistToArray(&points[i].adjacent,
+ &_adjacencies[i].t_length);
+
+}
+
+void DocRecord::voronoiCell (PointNumero pt, std::vector<SPoint2> &pts) const {
+
+ if (!_adjacencies){
+ printf("no adjacencies were created\n");
+ throw;
+ }
+ const int n = _adjacencies[pt].t_length;
+ for(int j = 0; j < n; j++) {
+ PointNumero a = _adjacencies[pt].t[j];
+ PointNumero b = _adjacencies[pt].t[(j+1) %n];
+ double pa[2] = {(double)points[a].where.h, (double)points[a].where.v};
+ double pb[2] = {(double)points[b].where.h, (double)points[b].where.v};
+ double pc[2] = {(double)points[pt].where.h, (double)points[pt].where.v};
+ double center[2];
+ circumCenterXY (pa,pb,pc,center);
+ pts.push_back(SPoint2(center[0],center[1]));
+ }
+}
+
+/*
+ Consider N points {X_1, \dots, X_N}
+ We want to find the point X_P that verifies
+
+ min max | X_i - X_P | , i=1,\dots,N
+
+ L2 norm
+
+ min \int_V || X - X_P||^2 dv
+
+ => 2 \int_V || X - X_P|| dv = 0 => X_P is the centroid
+
+ min \int_V || X - X_P||^{2m} dv
+
+ => 2 \int_V || X - X_P||^{2m-1} dv = 0 => X_P is somewhere ...
+
+ now in infinite norm, how to find X_P ?
+
+*/
+
+void DocRecord::makePosView(std::string fileName) {
+ FILE *f = fopen(fileName.c_str(),"w");
+ if (_adjacencies){
+ fprintf(f,"View \"voronoi\" {\n");
+ for(PointNumero i = 0; i < numPoints; i++) {
+ std::vector<SPoint2> pts;
+ double pc[2] = {(double)points[i].where.h, (double)points[i].where.v};
+ if (!onHull(i)){
+ fprintf(f,"SP(%g,%g,%g) {%g};\n",pc[0],pc[1],0.0,(double)i);
+ voronoiCell (i,pts);
+ for (int j=0;j<pts.size();j++){
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
+ pts[j].x(),pts[j].y(),0.0,
+ pts[(j+1)%pts.size()].x(),pts[(j+1)%pts.size()].y(),0.0,
+ (double)i,(double)i);
+ }
+ }
+ else {
+ fprintf(f,"SP(%g,%g,%g) {%g};\n",pc[0],pc[1],0.0,-(double)i);
+ }
+ }
+ fprintf(f,"};\n");
+ }
+ fclose(f);
+}
+
+void centroidOfOrientedBox(std::vector<SPoint2> &pts,
+ const double &angle,
+ double &xc,
+ double &yc,
+ double &inertia) {
+ const int N = pts.size();
+
+ double sina = sin(angle);
+ double cosa = cos(angle);
+
+ double xmin = cosa* pts[0].x()+ sina*pts[0].y();
+ double xmax = cosa* pts[0].x()+ sina*pts[0].y();
+ double ymin = -sina* pts[0].x()+ cosa*pts[0].y();
+ double ymax = -sina* pts[0].x()+ cosa*pts[0].y();
+ for (int j=1;j<N;j++){
+ xmin = std::min(cosa* pts[j].x()+ sina*pts[j].y(),xmin);
+ ymin = std::min(-sina* pts[j].x()+ cosa*pts[j].y(),ymin);
+ xmax = std::max(cosa* pts[j].x()+ sina*pts[j].y(),xmax);
+ ymax = std::max(-sina* pts[j].x()+ cosa*pts[j].y(),ymax);
+ }
+ double XC = 0.5*(xmax+xmin);
+ double YC = 0.5*(ymax+ymin);
+ xc = XC*cosa - YC *sina;
+ yc = XC*sina + YC *cosa;
+ inertia = std::max(xmax-xmin,ymax-ymin);
+}
+
+void centroidOfPolygon(SPoint2 &pc,
+ std::vector<SPoint2> &pts,
+ double &xc,
+ double &yc,
+ double &inertia) {
+ double area_tot = 0;
+ SPoint2 center(0,0);
+ for (int j=0;j<pts.size();j++){
+ SPoint2 &pa =pts[j];
+ SPoint2 &pb =pts[(j+1)%pts.size()];
+ const double area = triangle_area2d(pa,pb,pc);
+ area_tot += area;
+ center += ((pa+pb+pc) * (area/3.0));
+ }
+ SPoint2 x = center * (1.0/area_tot);
+ inertia = 0;
+ for (int j=0;j<pts.size();j++){
+ SPoint2 &pa =pts[j];
+ SPoint2 &pb =pts[(j+1)%pts.size()];
+ const double area = triangle_area2d(pa,pb,pc);
+
+ const double b = sqrt ( (pa.x()-pb.x())*(pa.x()-pb.x()) +
+ (pa.y()-pb.y())*(pa.y()-pb.y()) );
+ const double h = 2.0 * area / b;
+ const double a = fabs((pb.x()-pa.x())*(pc.x()-pa.x())*
+ (pb.y()-pa.y())*(pc.y()-pa.y()))/b;
+
+ const double j2 = (h*b*b*b + h*a*b*b + h*a*a*b + b*h*h*h) / 12.0;
+
+ center = (pa+pb+pc) * (1.0/3.0);
+ const SPoint2 dx = x - center;
+ inertia += j2 + area*area*(dx.x()+dx.x()+dx.y()*dx.y());
+ }
+ xc = x.x();
+ yc = x.y();
+}
+
+double DocRecord::Lloyd(int type) {
+ fullMatrix<double> cgs(numPoints,2);
+ double inertia_tot;
+ for(PointNumero i = 0; i < numPoints; i++) {
+ PointRecord &pt = points[i];
+ std::vector<SPoint2> pts;
+ voronoiCell (i,pts);
+ double E;
+
+ if (!points[i].data){
+ SPoint2 p (pt.where.h,pt.where.v);
+ if (type == 0)
+ centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E);
+ else
+ centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E);
+ }
+ inertia_tot += E;
+ }
+
+ for(PointNumero i = 0; i < numPoints; i++) {
+ if (!points[i].data){
+ points[i].where.h = cgs(i,0);
+ points[i].where.v = cgs(i,1);
+ }
+ }
+ return inertia_tot;
+}
+
// Convertir les listes d'adjacence en triangles
int DocRecord::ConvertDListToTriangles()
{
@@ -495,14 +693,15 @@ int DocRecord::ConvertDListToTriangles()
int count = 0, count2 = 0;
PointNumero aa, bb, cc;
- STriangle *striangle = (STriangle *) Malloc(n * sizeof(STriangle));
+ STriangle *striangle = new STriangle[n];
- count2 = CountPointsOnHull(n);
+ count2 = CountPointsOnHull();
// nombre de triangles que l'on doit obtenir
count2 = 2 * (n - 1) - count2;
- triangles = (Triangle *) Malloc(2 * count2 * sizeof(Triangle));
+ // FIXME ::: WHY 2 * ???????????????????
+ triangles = new Triangle[2 * count2];
for(i = 0; i < n; i++) {
// on cree une liste de points connectes au point i (t) + nombre
@@ -530,8 +729,8 @@ int DocRecord::ConvertDListToTriangles()
numTriangles = count2;
for(int i = 0; i < n; i++)
- Free(striangle[i].t);
- Free(striangle);
+ delete [] striangle[i].t;
+ delete striangle;
return 1;
}
@@ -548,23 +747,29 @@ void DocRecord::RemoveAllDList()
do {
temp = p;
p = Pred(p);
- Free(temp);
+ delete temp;
} while(p != points[i].adjacent);
points[i].adjacent = NULL;
}
}
DocRecord::DocRecord(int n)
- : numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
+ : numPoints(n), points(NULL), numTriangles(0), triangles(NULL), _hullSize(0), _hull(NULL),_adjacencies(NULL)
{
if(numPoints)
- points = (PointRecord*)Malloc(numPoints * sizeof(PointRecord));
+ points = new PointRecord[numPoints];
}
DocRecord::~DocRecord()
{
- if(points) Free(points);
- if(triangles) Free(triangles);
+ if(points) delete [] points;
+ if(triangles) delete []triangles;
+ if(_hull) delete [] _hull;
+ if (_adjacencies){
+ for(int i = 0; i < numPoints; i++)
+ delete [] _adjacencies[i].t;
+ delete _adjacencies;
+ }
}
void DocRecord::MakeMeshWithPoints()
@@ -575,3 +780,49 @@ void DocRecord::MakeMeshWithPoints()
RemoveAllDList();
}
+void DocRecord::Voronoi()
+{
+ if(numPoints < 3) return;
+ BuildDelaunay();
+ ConvertDListToVoronoiData();
+}
+
+void DocRecord::setPoints(fullMatrix<double> *p){
+ if (numPoints != p->size1())throw;
+ for (int i=0;i<p->size1();i++){
+ x(i) = (*p)(i,0);
+ y(i) = (*p)(i,1);
+ data(i) = (*p)(i,2) < 0 ? (void *) 1 : NULL;
+ }
+}
+
+
+#include "Bindings.h"
+
+void DocRecord::registerBindings(binding *b)
+{
+ classBinding *cb = b->addClass<DocRecord>("Triangulator");
+ cb->setDescription("A class that does 2D delaunay triangulation (JF's SANDBOX for the moment)");
+ methodBinding *cm;
+
+ cm = cb->addMethod("setPoints", &DocRecord::setPoints);
+ cm->setDescription("Set the NumPoints points of the triangulation (x,y,fixed)");
+ cm->setArgNames("points",NULL);
+ cm = cb->addMethod("Triangulate", &DocRecord::MakeMeshWithPoints);
+ cm->setDescription("Compute the Delaunay triangulation");
+ cm = cb->addMethod("Voronoi", &DocRecord::Voronoi);
+ cm->setDescription("Compute the Voronoi cells");
+ cm = cb->addMethod("hullSize", &DocRecord::hullSize);
+ cm->setDescription("returns the size of the hull");
+ cm = cb->addMethod("makePosView", &DocRecord::makePosView);
+ cm->setDescription("save a .pos file with the voronoi");
+ cm->setArgNames("FileName",NULL);
+ cm = cb->addMethod("Lloyd", &DocRecord::Lloyd);
+ cm->setDescription("do one iteration of Lloyd's algorithm");
+ cm->setArgNames("Type",NULL);
+ cm = cb->setConstructor<DocRecord,int>();
+ cm->setDescription ("A Triangulator");
+ cm->setArgNames("NumPoints",NULL);
+}
+
+
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index 5f23ab104773420f66f624fabac294b9fcb935c9..90f04276d431fc5b063a826b44f7a4a14a70067c 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -5,7 +5,13 @@
#ifndef _DIVIDE_AND_CONQUER_H_
#define _DIVIDE_AND_CONQUER_H_
+#include <vector>
+#include <algorithm>
+#include "fullMatrix.h"
+#include "SPoint2.h"
+class binding;
+class GFace;
typedef struct _CDLIST DListRecord, *DListPeek;
typedef int PointNumero;
@@ -14,11 +20,12 @@ typedef struct{
double h;
}DPoint;
-typedef struct{
+struct PointRecord{
DPoint where;
DListPeek adjacent;
void *data;
-}PointRecord;
+ PointRecord () : data (0) , adjacent(0) {}
+};
struct _CDLIST{
PointNumero point_num;
@@ -46,6 +53,9 @@ typedef struct{
class DocRecord{
private:
+ int _hullSize;
+ PointNumero *_hull;
+ STriangle * _adjacencies;
PointNumero Predecessor(PointNumero a, PointNumero b);
PointNumero Successor(PointNumero a, PointNumero b);
int FixFirst(PointNumero x, PointNumero f);
@@ -57,15 +67,17 @@ class DocRecord{
int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k);
int Merge(DT vl, DT vr);
DT RecurTrig(PointNumero left, PointNumero right);
- int BuildDelaunay();
int DListInsert(DListRecord **dlist, DPoint center, PointNumero newPoint);
int Insert(PointNumero a, PointNumero b);
int DListDelete(DListPeek *dlist, PointNumero oldPoint);
int Delete(PointNumero a, PointNumero b);
- int CountPointsOnHull(int n);
PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n);
int ConvertDListToTriangles();
+ void ConvertDListToVoronoiData();
void RemoveAllDList();
+ int BuildDelaunay();
+ int CountPointsOnHull();
+ void ConvexHull();
public:
int numPoints; // number of points
PointRecord *points; // points to triangulate
@@ -74,8 +86,29 @@ class DocRecord{
DocRecord(int n);
double &x(int i){ return points[i].where.v; }
double &y(int i){ return points[i].where.h; }
+ void* &data(int i){ return points[i].data; }
+ void setPoints(fullMatrix<double> *p);
~DocRecord();
void MakeMeshWithPoints();
+ void Voronoi ();
+ int hullSize() {return _hullSize;}
+ int onHull(PointNumero i) {return std::binary_search(_hull, _hull+_hullSize, i);}
+ void makePosView(std::string);
+ double Lloyd (int);
+ void voronoiCell (PointNumero pt, std::vector<SPoint2> &pts) const;
+ static void registerBindings(binding *b);
};
+void centroidOfOrientedBox(std::vector<SPoint2> &pts,
+ const double &angle,
+ double &xc,
+ double &yc,
+ double &inertia);
+void centroidOfPolygon(SPoint2 &pc,
+ std::vector<SPoint2> &pts,
+ double &xc,
+ double &yc,
+ double &inertia);
+
+
#endif
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 70ec4ba66228dcd37f63dc10a9dd11b102cca838..6166c77bce209cf36c97a3307dc35429780b4e4a 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -281,6 +281,21 @@ double triangle_area(double p0[3], double p1[3], double p2[3])
return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
}
+double triangle_area2d(double p0[2], double p1[2], double p2[2])
+{
+ const double c =
+ (p2[0] - p1[0])*(p0[1] - p1[1]) -
+ (p2[1] - p1[1])*(p0[0] - p1[0]);
+
+ return 0.5 * sqrt(c*c);
+}
+
+double triangle_polar_inertia(double p0[2], double p1[2], double p2[2])
+{
+ throw;
+}
+
+
void circumCenterXY(double *p1, double *p2, double *p3, double *res)
{
double d, a1, a2, a3;
@@ -294,7 +309,7 @@ void circumCenterXY(double *p1, double *p2, double *p3, double *res)
d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
if(d == 0.0) {
- Msg::Warning("Colinear points in circum circle computation");
+ // Msg::Warning("Colinear points in circum circle computation");
res[0] = res[1] = -99999.;
return ;
}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 631683a9c04b2e4d7653981aedd76c84923fd604..605b6688273a19cad90bdcf0ba84b812f79307f0 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -62,6 +62,8 @@ double inv2x2(double mat[2][2], double inv[2][2]);
double angle_02pi(double A3);
double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
double triangle_area(double p0[3], double p1[3], double p2[3]);
+double triangle_area2d(double p0[2], double p1[2], double p2[2]);
+double triangle_polar_inertia(double p0[2], double p1[2], double p2[2]);
void circumCenterXY(double *p1, double *p2, double *p3, double *res);
void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv=0);
char float2char(float f);
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 97f7a42375105eca659dc3d571d0fc0537317a35..f7d154e46f0af611d77f6785a0374be3ccceb46e 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -308,6 +308,9 @@ void fullMatrix<double>::registerBindings(binding *b)
cm = cb->addMethod("gemm", &fullMatrix<double>::gemm);
cm->setArgNames("A","B","alpha","beta",NULL);
cm->setDescription("this = beta*this + alpha * (A.B)");
+ cm = cb->addMethod("gemm_naive", &fullMatrix<double>::gemm_naive);
+ cm->setArgNames("A","B","alpha","beta",NULL);
+ cm->setDescription("this = beta*this + alpha * (A.B)");
cm = cb->setConstructor<fullMatrix<double>,int,int>();
cm->setDescription ("A new matrix of size 'nRows' x 'nColumns'");
cm->setArgNames("nRows","nColumns",NULL);
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 49e91e46b430fca1c5968db6770fa1ada50bfe24..dd6018c3209e9275bef5bd8045ac2e770c9a355e 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -197,8 +197,7 @@ class fullMatrix
for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
(*this)(desti, destj) = a(i, j);
}
- void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
-#if !defined(HAVE_BLAS)
+ void mult_naive(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
{
c.scale(0.);
for(int i = 0; i < _r; i++)
@@ -206,18 +205,30 @@ class fullMatrix
for(int k = 0; k < _c; k++)
c._data[i + _r * j] += (*this)(i, k) * b(k, j);
}
-#endif
;
- void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
- scalar alpha=1., scalar beta=1.)
+ void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
#if !defined(HAVE_BLAS)
+ {
+ mult_naive(b,c);
+ }
+#endif
+ ;
+ void gemm_naive(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
+ scalar alpha=1., scalar beta=1.)
{
fullMatrix<scalar> temp(a.size1(), b.size2());
- a.mult(b, temp);
+ a.mult_naive(b, temp);
temp.scale(alpha);
scale(beta);
add(temp);
}
+ ;
+ void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
+ scalar alpha=1., scalar beta=1.)
+#if !defined(HAVE_BLAS)
+ {
+ gemm_naive(anb,alpha,beta);
+ }
#endif
;
inline void setAll(const scalar &m)
diff --git a/Solver/TESTCASES/CylinderEddies.lua b/Solver/TESTCASES/CylinderEddies.lua
index 92987c61d022ac23b26d132629cfa003ab93186a..b35dc895fef726a621989fe53b4785792f5eded8 100644
--- a/Solver/TESTCASES/CylinderEddies.lua
+++ b/Solver/TESTCASES/CylinderEddies.lua
@@ -71,7 +71,7 @@ for i=1,2 do
end
if (i % 100 == 0) then
DG:exportSolution(string.format("output/cyl-%06d", i))
- DG:saveSolution(string.format("output/cyl-%06d.bin", i))
+-- DG:saveSolution(string.format("output/cyl-%06d.bin", i))
end
end
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index cec41f9134b95ab020b52d427c130682e5a30c6d..8d59bd00fab50acf1160bee46f4cb0d8f58364dc 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -14,7 +14,7 @@ r = 1.0 * unit ;
ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
ssin = ( 1.0 - ccos*ccos )^0.5 ;
-Lc1 = 0.01 ;
+Lc1 = 0.003 ;
Lc2 = 0.003 ;
Point(1) = { -e1-e2, 0.0 , 0.0 , Lc1};