lpcvt

CMakeLists.txt

@@ -43,7 +43,7 @@ option(ENABLE_OCC "Enable Open CASCADE geometrical models" ON)
 option(ENABLE_ACIS "Enable ACIS geometrical models" ON)
 option(ENABLE_OSMESA "Use OSMesa for offscreen rendering" OFF)
 option(ENABLE_PARSER "Build the GEO file parser" ON)
-option(ENABLE_PETSC "Enable PETSc linear algebra solvers" ON)
+option(ENABLE_PETSC "Enable PETSc linear algebra solvers" OFF)
 option(ENABLE_POST "Build the post-processing module" ON)
 option(ENABLE_PLUGINS "Build the post-processing plugins" ON)
 option(ENABLE_QT "Build QT GUI" OFF)
@@ -439,6 +439,9 @@ if(ENABLE_CHACO)
   set_config_option(HAVE_CHACO "Chaco")
 endif(ENABLE_CHACO)
 
+add_subdirectory(contrib/lbfgs)
+include_directories(contrib/lbfgs)
+
 if(ENABLE_DINTEGRATION)
   add_subdirectory(contrib/DiscreteIntegration)
   include_directories(contrib/DiscreteIntegration)

Common/LuaBindings.cpp

@@ -33,6 +33,8 @@
 #include "polynomialBasis.h"
 #include "Gauss.h"
 #include "meshPartitionOptions.h"
+#include "meshGFaceOptimize.h"
+#include "meshGFaceLloyd.h"
 
 #if defined(HAVE_OPENGL)
 #include "drawContext.h"
@@ -418,6 +420,7 @@ binding::binding()
   polynomialBasis::registerBindings(this);
   gaussIntegration::registerBindings(this);
   meshPartitionOptions::registerBindings(this);
+  Temporary::registerBindings(this);
 #if defined(HAVE_SOLVER)
   function::registerBindings(this);
   linearSystem<double>::registerBindings(this);

Common/gmshpy.i

@@ -31,6 +31,7 @@
   #include "meshPartitionOptions.h"
   #include "linearSystemCSR.h"
   #include "elasticitySolver.h"
+  #include "meshGFaceLloyd.h"
   #include "PView.h"
   #include "PViewData.h"
   #include "PViewFactory.h"
@@ -125,5 +126,6 @@ namespace std {
 %include "SPoint2.h"
 %include "GPoint.h"  
 %include "functionPython.h"
+%include "meshGFaceLloyd.h"
 %include "DefaultOptions.h"
 

Mesh/meshGFaceLloyd.h

@@ -7,9 +7,17 @@
 #define _MESH_GFACE_LLOYD_H_
 #include "fullMatrix.h"
 #include "DivideAndConquer.h"
-#include "Bindings.h"
+#include "ap.h"
+#include "alglibinternal.h"
+#include "alglibmisc.h"
+#include "linalg.h"
+#include "optimization.h"
 
 class GFace;
+class boundary_edge;
+
+void function1_grad(const alglib::real_1d_array&,double&,alglib::real_1d_array&,void*);
+void topology(const alglib::real_1d_array &,double,void*);
 
 class lloydAlgorithm {
   int ITER_MAX;
@@ -20,19 +28,77 @@ class lloydAlgorithm {
   : ITER_MAX(itermax), infiniteNorm(infnorm) {}
   lloydAlgorithm() {}
   void operator () (GFace *);
-  double get_angle(double,double);
-  void get_rotation(double,double,double&,double&,double&,double&);
-  void get_unrotation(double,double,double&,double&,double&,double&);
-  void rotate(SPoint2,SPoint2&);
-  void unrotate(SPoint2,SPoint2&);
-  SPoint2 fixed_point(const std::vector<SPoint2>&);
-  double range(const std::vector<SPoint2>&,double,double,bool);
-  double gradient(const std::vector<SPoint2>&,double,bool);
-  bool point_side(SPoint2,double,bool);
-  double polygon_area(const std::vector<SPoint2>&);
-  SPoint2 line_intersection(SPoint2,SPoint2,double,bool,bool&);
   double optimize(int,int);
-  static void registerBindings(binding *b);
+  
+  void eval(DocRecord&,GFace*,std::vector<SVector3>&,std::vector<double>&,std::vector<SVector3>&,std::vector<double>&,int);
+  double total_energy(std::vector<double>);	
+  SVector3 numerical(DocRecord&,GFace*,int,int);
+  void write(DocRecord&,GFace*,int);
+  void write2(DocRecord&,GFace*);
+  SVector3 analytical(SPoint2,SVector3,double);
+  SVector3 area_centroid(SPoint2,SPoint2,SPoint2);
+	
+  double F(SPoint2,SPoint2,SPoint2,int);
+  SVector3 simple(SPoint2,SPoint2,SPoint2,int);
+  SVector3 dF_dC1(SPoint2,SPoint2,SPoint2,int);
+  SVector3 dF_dC2(SPoint2,SPoint2,SPoint2,int);
+  double f(SPoint2,SPoint2,int);
+  double df_dx(SPoint2,SPoint2,int);
+  double df_dy(SPoint2,SPoint2,int);
+  double L1(double,double);
+  double dL1_du();
+  double dL1_dv();
+  double L2(double,double);
+  double dL2_du();
+  double dL2_dv();
+  double L3(double,double);
+  double dL3_du();
+  double dL3_dv();
+  double Tx(double,double,SPoint2,SPoint2,SPoint2);
+  double dTx_dp1x(double,double);
+  double dTx_dp2x(double,double);
+  double dTx_dp3x(double,double);
+  double dTx_du(SPoint2,SPoint2,SPoint2);
+  double dTx_dv(SPoint2,SPoint2,SPoint2);
+  double Ty(double,double,SPoint2,SPoint2,SPoint2);
+  double dTy_dp1y(double,double);
+  double dTy_dp2y(double,double);
+  double dTy_dp3y(double,double);
+  double dTy_du(SPoint2,SPoint2,SPoint2);
+  double dTy_dv(SPoint2,SPoint2,SPoint2);
+  double J(SPoint2,SPoint2,SPoint2);
+  double dJ_dp1x(SPoint2,SPoint2,SPoint2);
+  double dJ_dp2x(SPoint2,SPoint2,SPoint2);
+  double dJ_dp3x(SPoint2,SPoint2,SPoint2);
+  double dJ_dp1y(SPoint2,SPoint2,SPoint2);
+  double dJ_dp2y(SPoint2,SPoint2,SPoint2);
+  double dJ_dp3y(SPoint2,SPoint2,SPoint2);
+  SVector3 inner_dFdx0(SVector3,SPoint2,SPoint2,SPoint2,SPoint2);
+  SVector3 boundary_dFdx0(SVector3,SPoint2,SPoint2,SPoint2,SVector3);
+  double exp(double,int);
+  
+  SVector3 normal(SPoint2,SPoint2);
+  SPoint2 mid(SPoint2,SPoint2);
+  void print_segment(SPoint2,SPoint2,std::ofstream&);
+  bool adjacent(const DocRecord&,int,int);
+  double triangle_area(SPoint2,SPoint2,SPoint2);
+  SPoint2 boundary_intersection(SPoint2,SPoint2,bool&,SVector3&);
+  bool inside_domain(double,double);
+  int get_number_boundary_edges();
+  boundary_edge get_boundary_edge(int);
+  SPoint2 intersection(SPoint2,SPoint2,SPoint2,SPoint2,bool&);
+};
+
+class boundary_edge{
+ private :
+  SPoint2 p1;
+  SPoint2 p2;
+ public :
+  boundary_edge(SPoint2,SPoint2);
+  boundary_edge();
+  ~boundary_edge();
+  SPoint2 get_p1();
+  SPoint2 get_p2();
 };
 
 #endif

Mesh/meshGFaceOptimize.cpp

@@ -19,6 +19,10 @@
 #include "Generator.h"
 #include "Context.h"
 #include "OS.h"
+#include "PView.h"
+#include "PViewData.h"
+#include "SVector3.h"
+#include "SPoint3.h"
 
 #ifdef HAVE_MATCH
 extern "C" int FAILED_NODE;
@@ -1512,6 +1516,9 @@ struct RecombineTriangle
 {
   MElement *t1, *t2;
   double angle;
+  double cost_quality; //addition for class Temporary
+  double cost_alignment; //addition for class Temporary
+  double total_cost; //addition for class Temporary
   MVertex *n1, *n2, *n3, *n4;
   RecombineTriangle(const MEdge &me, MElement *_t1, MElement *_t2)
     : t1(_t1), t2(_t2)
@@ -1534,10 +1541,15 @@ struct RecombineTriangle
     angle = std::max(fabs(90. - a2),angle);
     angle = std::max(fabs(90. - a3),angle);
     angle = std::max(fabs(90. - a4),angle);
+	cost_quality = 1.0 - std::max(1.0 - angle/90.0,0.0); //addition for class Temporary
+	cost_alignment = Temporary::compute_alignment(me,_t1,_t2); //addition for class Temporary
+	total_cost = Temporary::compute_total_cost(cost_quality,cost_alignment); //addition for class Temporary
+	total_cost = 100.0*cost_quality; //addition for class Temporary
   }
   bool operator < (const RecombineTriangle &other) const
   {
-    return angle < other.angle;
+    //return angle < other.angle;
+	return total_cost < other.total_cost; //addition for class Temporary
   }
 };
 
@@ -1640,7 +1652,8 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
       for (int i=0;i<pairs.size();++i){
 	elist[2*i]   = t2n[pairs[i].t1];
 	elist[2*i+1] = t2n[pairs[i].t2];
-	elen [i] =  (int) pairs[i].angle;
+	//elen [i] =  (int) pairs[i].angle;
+	elen [i] = (int) pairs[i].total_cost; //addition for class Temporary
 	double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
 			     pairs[i].n1->x()-pairs[i].n2->x());
 
@@ -2055,3 +2068,166 @@ void recombineIntoQuadsIterative(GFace *gf)
     if (COUNT == 5)break;
   }
 }
+
+
+
+#include "Bindings.h"
+
+double Temporary::w1,Temporary::w2,Temporary::w3;
+std::vector<SVector3> Temporary::gradients;
+
+Temporary::Temporary(){}
+
+Temporary::~Temporary(){}
+
+SVector3 Temporary::compute_gradient(MElement*element){
+  double x1,y1,z1;
+  double x2,y2,z2;
+  double x3,y3,z3;
+  double x,y,z;
+  MVertex*vertex1 = element->getVertex(0);
+  MVertex*vertex2 = element->getVertex(1);
+  MVertex*vertex3 = element->getVertex(2);
+  x1 = vertex1->x();
+  y1 = vertex1->y();
+  z1 = vertex1->z();
+  x2 = vertex2->x();
+  y2 = vertex2->y();
+  z2 = vertex2->z();
+  x3 = vertex3->x();
+  y3 = vertex3->y();
+  z3 = vertex3->z();
+  x = (x1+x2+x3)/3.0;
+  y = (y1+y2+y3)/3.0;
+  z = (z1+z2+z3)/3.0;
+  return SVector3(0.0,1.0,0.0);
+}
+
+void Temporary::select_weights(double new_w1,double new_w2,double new_w3){
+  w1 = new_w1;
+  w2 = new_w2;
+  w3 = new_w3;
+}
+
+double Temporary::compute_total_cost(double f1,double f2){
+  double cost;
+  cost = w1*f1 + w2*f2 + w3*f1*f2;
+  return cost;
+}
+
+SVector3 Temporary::compute_normal(MElement*element){
+  double x1,y1,z1;
+  double x2,y2,z2;
+  double x3,y3,z3;
+  double length;
+  SVector3 vectorA;
+  SVector3 vectorB;
+  SVector3 normal;
+  MVertex*vertex1 = element->getVertex(0);
+  MVertex*vertex2 = element->getVertex(1);
+  MVertex*vertex3 = element->getVertex(2);
+  x1 = vertex1->x();
+  y1 = vertex1->y();
+  z1 = vertex1->z();
+  x2 = vertex2->x();
+  y2 = vertex2->y();
+  z2 = vertex2->z();
+  x3 = vertex3->x();
+  y3 = vertex3->y();
+  z3 = vertex3->z();
+  vectorA = SVector3(x2-x1,y2-y1,z2-z1);
+  vectorB = SVector3(x3-x1,y3-y1,z3-z1);
+  normal = crossprod(vectorA,vectorB);
+  length = norm(normal);
+  return SVector3(normal.x()/length,normal.y()/length,normal.z()/length);
+}
+
+SVector3 Temporary::compute_other_vector(MElement*element){
+  int number;
+  double length;
+  SVector3 normal;
+  SVector3 gradient;
+  SVector3 other_vector;
+  number = element->getNum();
+  normal = Temporary::compute_normal(element);
+  gradient = Temporary::compute_gradient(element);//gradients[number];
+  other_vector = crossprod(gradient,normal);
+  length = norm(other_vector);
+  return SVector3(other_vector.x()/length,other_vector.y()/length,other_vector.z()/length);
+}
+
+double Temporary::compute_alignment(const MEdge&_edge, MElement*element1, MElement*element2){
+  int number;
+  double scalar_productA,scalar_productB;
+  double alignment;
+  SVector3 gradient;
+  SVector3 other_vector;
+  SVector3 edge;
+  MVertex*vertexA;
+  MVertex*vertexB;
+  number = element1->getNum();
+  gradient = Temporary::compute_gradient(element1);//gradients[number];
+  other_vector = Temporary::compute_other_vector(element1);
+  vertexA = _edge.getVertex(0);
+  vertexB = _edge.getVertex(1);
+  edge = SVector3(vertexB->x()-vertexA->x(),vertexB->y()-vertexA->y(),vertexB->z()-vertexA->z());
+  scalar_productA = fabs(dot(gradient,edge));
+  scalar_productB = fabs(dot(other_vector,edge));
+  alignment = std::max(scalar_productA,scalar_productB) - sqrt(2.0)/2.0;
+  alignment = alignment/(1.0-sqrt(2.0)/2.0);	
+  return alignment;
+}
+
+void Temporary::read_data(std::string file_name){
+  int i,j,number;
+  double x,y,z;
+  MElement*element;
+  PView*view;
+  PViewData*data;	
+  PView::readMSH(file_name,-1);
+  std::vector<PView*> list = PView::list;
+  data = list[0]->getData();
+  for(i = 0;i < data->getNumEntities(0);i++)
+  {
+    if(data->skipEntity(0,i)) continue;
+    for(j = 0;j < data->getNumElements(0,i);j++)
+	{
+	  if(data->skipElement(0,i,j)) continue;
+	  element = data->getElement(0,i,j);
+	  number = element->getNum();
+	  data->getValue(0,i,j,0,x);
+	  data->getValue(0,i,j,1,y);
+	  data->getValue(0,i,j,2,z);
+	  gradients[number] = SVector3(x,y,z);
+	}
+  }
+}
+
+void Temporary::quadrilaterize(std::string file_name,double _w1,double _w2,double _w3){
+  GFace*face;
+  GModel*model = GModel::current();
+  GModel::fiter iterator;
+  gradients.resize(model->getNumMeshElements() + 1);
+  w1 = _w1;
+  w2 = _w2;
+  w3 = _w3;
+  Temporary::read_data(file_name);
+  for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
+  {
+    face = *iterator;
+	_recombineIntoQuads(face,1,1);
+  }
+}
+
+void Temporary::registerBindings(binding *b){
+  classBinding*cb = b->addClass<Temporary>("Temporary");
+  cb->setDescription("This class is used to create quad meshes from a script.");
+  methodBinding*cm;
+  cm = cb->setConstructor<Temporary>();
+  cm->setDescription("This is the constructor.");
+  cm = cb->addMethod("quadrilaterize",&Temporary::quadrilaterize);
+  cm->setDescription("This function creates the quad mesh.");
+  cm->setArgNames("file_name","w1","w2","w3",NULL);
+}
+	
+

Mesh/meshGFaceOptimize.h

@@ -126,4 +126,22 @@ struct swapquad{
   }
 };
 
+class Temporary{
+  private :
+	static double w1,w2,w3;
+	static std::vector<SVector3> gradients;
+	void read_data(std::string);
+	static SVector3 compute_normal(MElement*);
+	static SVector3 compute_other_vector(MElement*);
+	static SVector3 compute_gradient(MElement*);
+  public :
+	Temporary();
+	~Temporary();
+	void quadrilaterize(std::string,double,double,double);
+	static double compute_total_cost(double,double);
+	static void select_weights(double,double,double);
+	static double compute_alignment(const MEdge&,MElement*,MElement*);
+	static void registerBindings(binding *b);
+};
+
 #endif

Numeric/DivideAndConquer.cpp

@@ -28,6 +28,30 @@
 #define Pred(x) ((x)->prev)
 #define Succ(x) ((x)->next)
 
+int DocRecord::get_p(){
+  return p;
+}
+
+void DocRecord::set_p(int new_p){
+  p = new_p;
+}
+
+int DocRecord::get_dimension(){
+  return dimension;
+}
+
+void DocRecord::set_dimension(int new_dimension){
+  dimension = new_dimension;
+}
+
+GFace* DocRecord::get_face(){
+  return gf;
+}
+
+void DocRecord::set_face(GFace* new_gf){
+  gf = new_gf;
+}
+
 PointNumero DocRecord::Predecessor(PointNumero a, PointNumero b)
 {
   DListPeek p = points[a].adjacent;
@@ -824,7 +848,7 @@ DocRecord::DocRecord(int n)
     numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
 {
   if(numPoints)
-    points = new PointRecord[numPoints+1000];
+    points = new PointRecord[numPoints+3000];
 }
 
 DocRecord::~DocRecord()
@@ -898,6 +922,7 @@ void DocRecord::remove_all(){
 	  points2[index].where.v = points[i].where.v;
 	  points2[index].data = points[i].data;
 	  points2[index].flag = points[i].flag;
+	  points2[index].identificator = points[i].identificator;
 	  index++;
 	}
   }
@@ -906,6 +931,15 @@ void DocRecord::remove_all(){
   numPoints = numPoints2;
 }
 
+void DocRecord::add_point(double x,double y,GFace*face){
+  PointRecord point;
+  point.where.h = x;
+  point.where.v = y; 
+  point.data = new MVertex(x,y,0.0,(GEntity*)face,2);
+  points[numPoints] = point;
+  numPoints = numPoints+1;
+}
+
 #include "Bindings.h"
 
 void DocRecord::registerBindings(binding *b)

Numeric/DivideAndConquer.h

@@ -30,6 +30,7 @@ struct PointRecord {
   DListPeek adjacent;
   void *data;
   int flag; //0:to be kept, 1:to be removed
+  int identificator;
   PointRecord() : adjacent(0), data (0) {}
 };
 
@@ -59,6 +60,9 @@ typedef struct{
 
 class DocRecord{
  private:
+  int p;
+  int dimension;
+  GFace* gf;
   int _hullSize;
   PointNumero *_hull;
   PointNumero Predecessor(PointNumero a, PointNumero b);
@@ -76,15 +80,22 @@ class DocRecord{
   int Insert(PointNumero a, PointNumero b);
   int DListDelete(DListPeek *dlist, PointNumero oldPoint);
   int Delete(PointNumero a, PointNumero b);
-  int ConvertDListToTriangles();
   void ConvertDListToVoronoiData();
+  int ConvertDListToTriangles();
   void RemoveAllDList();
   int BuildDelaunay();
   int CountPointsOnHull();
   void ConvexHull();
  public:
+  int get_p();
+  void set_p(int);
+  int get_dimension();
+  void set_dimension(int);
+  GFace* get_face();
+  void set_face(GFace*);
   STriangle *_adjacencies;
   int numPoints;        // number of points
+  int size_points;
   PointRecord *points;  // points to triangulate
   int numTriangles;     // number of triangles
   Triangle *triangles;  // 2D results
@@ -105,6 +116,7 @@ class DocRecord{
   void initialize();
   bool remove_point(int);
   void remove_all();
+  void add_point(double,double,GFace*);
   PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n);
   static void registerBindings(binding *b);
 };