diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 3ed8303ef171e68693929c573f1d2788fe776f80..2304d1f6a1024e905463052094688e57ebfd8775 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -457,10 +457,24 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
- std::list<GEdge*> &U1, std::list<GEdge*> &V1)
- : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0)
+ std::list<GEdge*> &U1, std::list<GEdge*> &V1,
+ gmshLinearSystem<double> * lsys)
+ : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0),
+ _lsys(lsys)
{
+#ifdef HAVE_GMM
+ if (!_lsys) {
+ gmshLinearSystemGmm<double> *_lsysb = new gmshLinearSystemGmm<double>;
+ //gmshLinearSystemCSRGmm<double> lsys;
+ _lsysb->setPrec(1.e-15);
+ if(Msg::GetVerbosity() == 99) _lsysb->setNoisy(2);
+ _lsys = _lsysb;
+ }
+#else
+ _lsys = new gmshLinearSystemFull<double>;
+#endif
+
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
if(!(*it)){
Msg::Error("Incorrect face in compound surface %d\n", tag);
@@ -651,15 +665,7 @@ void GFaceCompound::parametrize(iterationStep step) const
{
Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
-#ifdef HAVE_GMM
- gmshLinearSystemGmm<double> lsys;
- //gmshLinearSystemCSRGmm<double> lsys;
- lsys.setPrec(1.e-15);
- if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
-#else
- gmshLinearSystemFull<double> lsys;
-#endif
- gmshAssembler<double> myAssembler(&lsys);
+ gmshAssembler<double> myAssembler(_lsys);
gmshFunction<double> ONE(1.0);
@@ -728,7 +734,7 @@ void GFaceCompound::parametrize(iterationStep step) const
}
Msg::Debug("Assembly done");
- lsys.systemSolve();
+ _lsys->systemSolve();
Msg::Debug("System solved");
@@ -755,8 +761,8 @@ void GFaceCompound::parametrize(iterationStep step) const
else{
itf->second[step]= value;
}
-
}
+ _lsys->clear();
}
void GFaceCompound::parametrize_conformal() const
@@ -765,16 +771,8 @@ void GFaceCompound::parametrize_conformal() const
Msg::Info("Parametrizing Surface %d", tag());
-#ifdef HAVE_GMM
- gmshLinearSystemGmm<double> lsys;
- //gmshLinearSystemCSRGmm<double> lsys;
- lsys.setPrec(1.e-15);
- lsys.setGmres(1);
- if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
-#else
- gmshLinearSystemFull<double> lsys;
-#endif
- gmshAssembler<double> myAssembler(&lsys);
+ gmshAssembler<double> myAssembler(_lsys);
+
std::vector<MVertex*> ordered;
std::vector<double> coords;
bool success = orderVertices(_U0, ordered, coords);
@@ -863,7 +861,7 @@ void GFaceCompound::parametrize_conformal() const
}
Msg::Debug("Assembly done");
- lsys.systemSolve();
+ _lsys->systemSolve();
Msg::Debug("System solved");
it = _compound.begin();
@@ -884,10 +882,10 @@ void GFaceCompound::parametrize_conformal() const
coordinates[v] = SPoint3(value1,value2,0.0);
}
- // checkOrientation();
-
+ // checkOrientation();
computeNormals();
buildOct();
+ _lsys->clear();
}
void GFaceCompound::computeNormals(std::map<MVertex*,SVector3> &normals) const
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 376c9cabae1de7e098821a92ae2ccc1c5e6de807..69755e9c5cf9bed11a813de85c44766a7ae8eec0 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -12,6 +12,7 @@
#include "GEdge.h"
#include "GEdgeCompound.h"
#include "meshGFaceOptimize.h"
+#include "gmshLinearSystem.h"
/*
A GFaceCompound is a model face that is the compound of model faces.
@@ -73,11 +74,13 @@ class GFaceCompound : public GFace {
virtual double curvature(MTriangle *t, double u, double v) const;
void printStuff() const;
bool trivial() const ;
+ gmshLinearSystem <double> *_lsys;
public:
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &U1,
- std::list<GEdge*> &V0, std::list<GEdge*> &V1);
+ std::list<GEdge*> &V0, std::list<GEdge*> &V1,
+ gmshLinearSystem<double>* lsys =0);
virtual ~GFaceCompound();
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index d9522cb04d21fecb8ca16993b2b5c4838959d996..20675b626c5dc252a4a1d447f051349d85be1efd 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -188,7 +188,7 @@ static void basisFuns(double u, int i, int deg, float *U, double *N)
static Vertex InterpolateNurbs(Curve *Curve, double u, int derivee)
{
- static double Nb[1000];
+ double Nb[1000];
int span = findSpan(u, Curve->degre, List_Nbr(Curve->Control_Points), Curve->k);
basisFuns(u, span, Curve->degre, Curve->k, Nb);
Vertex p;
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 078c84216347f01fb03b28b9039a9fa7302f9eb8..a0703b31677f2501933347d99b4e23faaba04e1d 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -13,6 +13,10 @@
#include "meshPartition.h"
#include "meshPartitionObjects.h"
#include "meshPartitionOptions.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MElementCut.h"
+#include "partitionEdge.h"
//--Prototype for Chaco interface
@@ -643,6 +647,87 @@ struct MatchBoElemToGrVertex
}
};
+template <class ITERATOR>
+void fillit_ ( std::multimap<MEdge,MElement*,Less_Edge> &edgeToElement,
+ ITERATOR it_beg, ITERATOR it_end){
+ for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
+ MElement *el = *IT;
+ for(int j=0;j < el->getNumEdges();j++) {
+ MEdge e = el->getEdge(j);
+ edgeToElement.insert(std::make_pair(e,el));
+ }// for every edge of the elem
+ }// for every elem of the face
+}
+
+template <class ITERATOR>
+void fillit_ ( std::multimap<MVertex*,MElement*> &vertexToElement,
+ ITERATOR it_beg, ITERATOR it_end){
+ for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
+ MElement *el = *IT;
+ for(int j=0;j < el->getNumVertices();j++) {
+ MVertex* e = el->getVertex(j);
+ vertexToElement.insert(std::make_pair(e,el));
+ }// for every edge of the elem
+ }// for every elem of the face
+}
+
+
+void assignPartitionBoundary (GModel *model,
+ MEdge &me,
+ std::set<partitionEdge*, Less_partitionEdge> &pedges,
+ std::vector<MElement*> &v){
+
+ std::vector<int> v2;
+ v2.push_back(v[0]->getPartition());
+ for (int i=1;i<v.size();i++){
+ bool found = false;
+ for (int j=0;j<v2.size();j++){
+ if (v[i]->getPartition() == v2[j]){
+ found = true;
+ break;
+ }
+ }
+ if (!found)v2.push_back(v[i]->getPartition());
+ }
+ if (v2.size() < 2)return;
+
+ partitionEdge pe (model, 1, 0, 0, v2);
+ std::set<partitionEdge*, Less_partitionEdge>::iterator it = pedges.find(&pe);
+ partitionEdge *ppe;
+ if (it == pedges.end()){
+ ppe = new partitionEdge(model, -pedges.size()-1, 0, 0, v2);
+ pedges.insert (ppe);
+ model->add(ppe);
+ }
+ else ppe = *it;
+ ppe->lines.push_back(new MLine (me.getVertex(0),me.getVertex(1)));
+}
+
+int CreatePartitionBoundaries (GModel *model) {
+ std::multimap<MEdge,MElement*,Less_Edge> edgeToElement;
+ std::set<partitionEdge*, Less_partitionEdge> pedges;
+ for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
+ fillit_ ( edgeToElement,(*it)->triangles.begin(),(*it)->triangles.end());
+ fillit_ ( edgeToElement,(*it)->quadrangles.begin(),(*it)->quadrangles.end());
+ fillit_ ( edgeToElement,(*it)->polygons.begin(),(*it)->polygons.end());
+ }
+
+ {
+ std::multimap<MEdge,MElement*,Less_Edge>::iterator it = edgeToElement.begin();
+ Equal_Edge oper;
+ while ( it != edgeToElement.end() ){
+ MEdge e = it->first;
+ std::vector<MElement*> voe;
+ do {
+ voe.push_back(it->second);
+ ++it;
+ }while (! oper (e,it->first));
+ assignPartitionBoundary (model,e,pedges,voe);
+ }
+ }
+}
+
+
/*******************************************************************************
*
diff --git a/Mesh/meshPartition.h b/Mesh/meshPartition.h
index 2c81367f4673d7b13d63ef55be2b66cd00fbc734..b9d40cdea91736285cfc11fa5a7c12df7ebb70de 100644
--- a/Mesh/meshPartition.h
+++ b/Mesh/meshPartition.h
@@ -25,5 +25,6 @@ int MakeGraph(GModel *const model, Graph &graph,
BoElemGrVec *const boElemGrVec = 0);
int PartitionGraph(Graph &graph, meshPartitionOptions &options);
int PartitionMesh(GModel *const model, meshPartitionOptions &options);
+int CreatePartitionBoundaries (GModel *model);
#endif
diff --git a/Numeric/gmshLinearSystem.h b/Numeric/gmshLinearSystem.h
index caf623604c0bc868a6adca7c1f0fa68b35e2f997..6e544ac8c399a579640addd542bc9207e361a82a 100644
--- a/Numeric/gmshLinearSystem.h
+++ b/Numeric/gmshLinearSystem.h
@@ -16,6 +16,7 @@ class gmshLinearSystem {
virtual ~gmshLinearSystem (){}
virtual bool isAllocated() const = 0;
virtual void allocate(int nbRows) = 0;
+ virtual void clear() = 0;
virtual void addToMatrix(int _row, int _col, scalar val) = 0;
virtual scalar getFromMatrix(int _row, int _col) const = 0;
virtual void addToRightHandSide(int _row, scalar val) = 0;
diff --git a/Numeric/gmshLinearSystemCSR.h b/Numeric/gmshLinearSystemCSR.h
index 8bf38ed74d3acfb4fefaef236db7df75b747dece..63d536dd0e19e2a1500dedf18c7e47b191aef0ec 100644
--- a/Numeric/gmshLinearSystemCSR.h
+++ b/Numeric/gmshLinearSystemCSR.h
@@ -37,6 +37,10 @@ class gmshLinearSystemCSR : public gmshLinearSystem<scalar> {
: sorted(false), a_(0) {}
virtual bool isAllocated() const { return a_ != 0; }
virtual void allocate(int) ;
+ virtual void clear()
+ {
+ allocate(0);
+ }
virtual ~gmshLinearSystemCSR()
{
allocate(0);
diff --git a/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
index ef3fda9a838cf83ec5d952d77e9a5f9e625518ff..ffc3292a75cc1442e6f6fe0a1982dd6122318ec9 100644
--- a/Numeric/gmshLinearSystemFull.h
+++ b/Numeric/gmshLinearSystemFull.h
@@ -24,18 +24,25 @@ class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
{
- if(_a) delete _a;
- if(_x) delete _x;
- if(_b) delete _b;
+ clear();
_a = new gmshMatrix<scalar>(_nbRows, _nbRows);
_b = new gmshVector<scalar>(_nbRows);
_x = new gmshVector<scalar>(_nbRows);
}
+
virtual ~gmshLinearSystemFull()
{
- delete _a;
- delete _b;
- delete _x;
+ clear();
+ }
+
+ virtual void clear()
+ {
+ if(_a){
+ delete _a;
+ delete _b;
+ delete _x;
+ }
+ _a = 0;
}
virtual void addToMatrix(int _row, int _col, scalar _val)
{
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
index 1d58acfdd9270cb59dd4b43b15bd5e82f8242ca3..83ee875c871c4190eb51e8625d67aa35691c2735 100644
--- a/Numeric/gmshLinearSystemGmm.h
+++ b/Numeric/gmshLinearSystemGmm.h
@@ -29,18 +29,24 @@ class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
{
- if(_a) delete _a;
- if(_x) delete _x;
- if(_b) delete _b;
+ clear();
_a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
_b = new std::vector<scalar>(_nbRows);
_x = new std::vector<scalar>(_nbRows);
}
virtual ~gmshLinearSystemGmm()
{
- delete _a;
- delete _b;
- delete _x;
+ clear();
+ }
+
+ virtual void clear()
+ {
+ if (_a){
+ delete _a;
+ delete _b;
+ delete _x;
+ }
+ _a = 0;
}
virtual void addToMatrix(int _row, int _col, scalar _val)
{
diff --git a/benchmarks/3d/Torus.geo b/benchmarks/3d/Torus.geo
index dccc9448e10231cdb9bcecd7f790a1fa3b81cf33..46241c987f6be8a8baf3a357eb48f10ca0422984 100644
--- a/benchmarks/3d/Torus.geo
+++ b/benchmarks/3d/Torus.geo
@@ -13,4 +13,4 @@ Line Loop(5) = {4,1,2,3};
Plane Surface(6) = {5};
Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2};
-Recombine Surface {6, 27, 23, 15, 19, 28};
+//Recombine Surface {6, 27, 23, 15, 19, 28};
diff --git a/benchmarks/boolean/FILLET.geo b/benchmarks/boolean/FILLET.geo
index 1a43702bac7f18fadb6bab37a335aeacf09f0441..b66b90f4c754f45b278d543d4c21ad3dac447da8 100644
--- a/benchmarks/boolean/FILLET.geo
+++ b/benchmarks/boolean/FILLET.geo
@@ -3,7 +3,9 @@ H = 30;
Z = 10;
OCCShape("Box",{0,0,0,L,H,Z},"none");
R = 4;
+
OCCShape("Fillet",{1:12,R},"none");
+/*
OCCShape("Cone",{0*L/2,H/2,-Z,0,0,1,.3*R,2*R,3*Z},"Fuse");
OCCShape("Fillet",{1,R},"none");
OCCShape("Fillet",{14,R/8},"none");
@@ -11,6 +13,3 @@ OCCShape("Cone",{0.99*L/2,H/2,-Z,0,0,1,.3*R,2*R,3*Z},"Fuse");
OCCShape("Fillet",{1,R},"none");
OCCShape("Fillet",{83,R/8},"none");
OCCShape("Cone",{0.99*L,H/2,-Z,0,0,1,3*R,.2*R,3*Z},"Cut");
-
-// not ready yet for seams!
-//Compound Surface(100000000) = {31, 40, 43, 13};