From 0b4486fd49c382560ead8d639647cb3391acc10f Mon Sep 17 00:00:00 2001
From: Bruno Seny <bruno.seny@student.uclouvain.be>
Date: Mon, 20 Dec 2010 16:35:26 +0000
Subject: [PATCH] Several things for multi-constraint mesh partitioning + start
 circumscribed circle of MElement (need to be finished)

---
 Geo/MElement.h              |  5 +++++
 Geo/MQuadrangle.cpp         |  6 +++++-
 Geo/MQuadrangle.h           |  2 ++
 Geo/MTriangle.cpp           | 16 +++++++++++++++-
 Geo/MTriangle.h             |  1 +
 Mesh/meshPartition.cpp      | 15 +++++++++++++--
 Mesh/meshPartitionOptions.h |  2 ++
 7 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/Geo/MElement.h b/Geo/MElement.h
index 0dd0b769be..38ba9eb16c 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -192,6 +192,11 @@ class MElement
   // otherwise get the minimum radius of all the circles/spheres
   // tangent to the most boundaries of the element.
   virtual double getInnerRadius(){ return 0.; }
+  
+  // get the radius of the circumscribed circle/sphere if it exists,
+  // otherwise get the maximum radius of all the circles/spheres
+  // tangent to the most boundaries of the element.
+  virtual double getOuterRadius(){ return 0.; }
 
   // compute the barycenter
   virtual SPoint3 barycenter();
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index bc9aa6f41b..ff011fea99 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -267,7 +267,11 @@ double MQuadrangle::angleShapeMeasure()
   return 1.;
 #endif
 }
-
+double MQuadrangle::getOuterRadius()
+{
+  // TO DO!!!!!!!!!!!!! (BRUNO SENY)
+  return 1.0;
+}
 double MQuadrangle::getInnerRadius()
 {
   // get the coordinates (x, y, z) of the 4 points defining the Quad
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 412224a8ff..8de5e45fcb 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -39,6 +39,7 @@ class MQuadrangle : public MElement {
     v[2] = _v[2];
     v[3] = _v[3];
   }
+void projectInMeanPlane(double *xn, double *yn);
  public :
   MQuadrangle(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0)
     : MElement(num, part)
@@ -141,6 +142,7 @@ class MQuadrangle : public MElement {
   // planar, we compute the mean plane due to the least-square
   // criterion.
   virtual double getInnerRadius();
+  virtual double getOuterRadius();
  private:
   int edges_quad(const int edge, const int vert) const
   {
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 92d0aab1fa..21dc2f496b 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -44,7 +44,21 @@ double MTriangle::getInnerRadius()
     dist[i] = e.getVertex(0)->distance(e.getVertex(1));
     k += 0.5 * dist[i];
   }
-  return sqrt(k * (k - dist[0]) * (k - dist[1]) * (k - dist[2])) / k;
+  double area = sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]));
+  return area/k;
+}
+
+double MTriangle::getOuterRadius()
+{
+  // radius of circle circumscribing a triangle
+  double dist[3], k = 0.;
+  for (int i = 0; i < 3; i++){
+    MEdge e = getEdge(i);
+    dist[i] = e.getVertex(0)->distance(e.getVertex(1));
+    k += 0.5 * dist[i];
+  }
+  double area = sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]));
+  return dist[0]*dist[1]*dist[2]/(4*area);
 }
 
 double MTriangle::angleShapeMeasure()
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 87150f25bc..fdd769b25e 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -54,6 +54,7 @@ class MTriangle : public MElement {
   virtual double gammaShapeMeasure();
   virtual double distoShapeMeasure();
   virtual double getInnerRadius();
+  virtual double getOuterRadius();
   virtual double angleShapeMeasure();
   virtual int getNumVertices() const { return 3; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 7b1283f913..95ee43016a 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -59,7 +59,7 @@ extern "C" void METIS_mCPartGraphRecursive
 (int *, int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *,
  int *, idxtype *);
 extern "C" void METIS_mCPartGraphKway
-(int *, int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *,
+(int *, int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, float *, int *,
  int *, idxtype *);
 
 
@@ -437,6 +437,7 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
         int numflag = 0;
         // if metisOptions[0]=0 then default options
         int metisOptions[5];
+        float ubvec[options.ncon];
         int edgeCut;
         const int iSec = 0;
         switch(options.algorithm) {
@@ -507,11 +508,21 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
           metisOptions[2] = 1;
           metisOptions[3] = options.refine_algorithm;
           metisOptions[4] = 0;
+          printf("Tolerance for Constraints:[");
+          for(int u=0;u<options.ncon;u++){
+           ubvec[u]=1.0;
+           if(options.tolerance[u]%options.num_partitions>0){
+             ubvec[u] = (float) ceil((float)options.tolerance[u]/options.num_partitions)/((float)options.tolerance[u]/options.num_partitions);
+           }
+           printf(" %f", ubvec[u]);
+          }
+          printf("] \n");
+          graph.fillWithMultipleWeights(options.ncon,options.getWeightMap());
           if (options.num_partitions > 1) {
             METIS_mCPartGraphKway
               (&n,&options.ncon,&graph.xadj[graph.section[iSec]],
                &graph.adjncy[graph.section[iSec]], &graph.vwgts[graph.section[iSec]], NULL, &wgtflag, &numflag,
-               &options.num_partitions, metisOptions, &edgeCut,
+               &options.num_partitions,ubvec, metisOptions, &edgeCut,
                &graph.partition[graph.section[iSec]]);
           }
           break;
diff --git a/Mesh/meshPartitionOptions.h b/Mesh/meshPartitionOptions.h
index b90bdd9ecc..549cb0af45 100644
--- a/Mesh/meshPartitionOptions.h
+++ b/Mesh/meshPartitionOptions.h
@@ -96,6 +96,8 @@ class meshPartitionOptions
     std::vector<int> nodalWeights;
 
     std::map<int, std::vector<int> > weightMap;
+
+    std::vector<int> tolerance;
   
   
   public:
-- 
GitLab