diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 115bb1ffecdaab8ef2aa9bf1596a0ca772b33e88..d5c7870885bcdd5868bd6ee3c823e6e54de5ca84 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -2442,7 +2442,7 @@ void GUI::create_option_window()
mesh_choice[5]->align(FL_ALIGN_RIGHT);
mesh_choice[5]->callback(mesh_options_ok_cb);
// not reimplemented yet
- ((Fl_Menu_Item*)mesh_choice[5]->menu())[1].deactivate();
+ // ((Fl_Menu_Item*)mesh_choice[5]->menu())[1].deactivate();
mesh_value[0] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 4 * BH, IW, BH, "Smoothing steps");
mesh_value[0]->minimum(0);
diff --git a/Fltk/GUI_Classifier.cpp b/Fltk/GUI_Classifier.cpp
index f2779e955069e728492adb4547a3c61383a5abad..61c21b863ecafcafd69b835c85a715372413090e 100644
--- a/Fltk/GUI_Classifier.cpp
+++ b/Fltk/GUI_Classifier.cpp
@@ -288,7 +288,7 @@ void updateedges_cb(Fl_Widget* w, void* data)
Draw();
}
-edge_angle:: edge_angle ( MVertex *_v1, MVertex *_v2, MTriangle *t1, MTriangle *t2)
+edge_angle:: edge_angle ( MVertex *_v1, MVertex *_v2, MElement *t1, MElement *t2)
: v1(_v1), v2(_v2)
{
if (!t2) angle = 0;
diff --git a/Fltk/GUI_Classifier.h b/Fltk/GUI_Classifier.h
index d0084ed96ac38f4fdd4c0e6a3028a8420843f541..e6ac0694703bd851b778c80adf82b1b0913b6821 100644
--- a/Fltk/GUI_Classifier.h
+++ b/Fltk/GUI_Classifier.h
@@ -42,7 +42,7 @@ class edge_angle
public :
MVertex *v1, *v2;
double angle;
- edge_angle ( MVertex *_v1, MVertex *_v2, MTriangle *t1, MTriangle *t2);
+ edge_angle ( MVertex *_v1, MVertex *_v2, MElement *t1, MElement *t2);
bool operator < (const edge_angle & other) const
{
return other.angle < angle;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index c96dfd1c18cf0a010e7357a2ba5c3026d6be957f..0b1649a93e6c7c1bd19a77bd582a6ca71cc000f0 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -9,6 +9,7 @@
#include "Numeric.h"
#include "Octree.h"
#include "gmshLinearSystemGmm.h"
+#include "gmshLinearSystemFull.h"
class gmshGradientBasedDiffusivity : public gmshFunction
{
@@ -201,9 +202,13 @@ void GFaceCompound::parametrize (bool _isU, int ITER) const
{
Msg::Info("Parametrizing Surface %d coordinate %d (ITER %d)", tag(), _isU, ITER);
+#ifdef HAVE_GMM
gmshLinearSystemGmm lsys;
lsys.setPrec(1.e-10);
- lsys.setNoisy(2);
+ //lsys.setNoisy(2);
+#else
+ gmshLinearSystemFull lsys;
+#endif
gmshAssembler myAssembler(&lsys);
gmshGradientBasedDiffusivity diffusivity (coordinates);
if (ITER > 0) diffusivity.setComponent(_isU);
@@ -212,10 +217,10 @@ void GFaceCompound::parametrize (bool _isU, int ITER) const
// maps the boundary onto a circle
std::vector<MVertex*> ordered;
std::vector<double> coords;
- Msg::Info("%d edges on the contour", l_edges.size());
- for (std::list<GEdge*>::const_iterator it = l_edges.begin();
- it !=l_edges.end();++it)printf("%d ",(*it)->tag());
- printf("\n");
+ // Msg::Info("%d edges on the contour", l_edges.size());
+ // for (std::list<GEdge*>::const_iterator it = l_edges.begin();
+ // it !=l_edges.end();++it)printf("%d ",(*it)->tag());
+ // printf("\n");
bool success = orderVertices (l_edges, ordered, coords);
if (!success)throw;
for (int i=0; i<ordered.size();i++){
@@ -355,7 +360,8 @@ GPoint GFaceCompound::point(double par1, double par2) const
getTriangle (par1, par2, <, U,V);
SPoint3 p(0,0,0);
if (!lt){
- printf("coucou\n");
+ Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
+
return GPoint(p.x(),p.y(),p.z(),this);
}
p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 1c8507f66758185bd04054cf409357ccbd4f3c09..25a8cd271daaedf4af95168a2e6d0f542895c7e9 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1551,6 +1551,3 @@ void BDS_Mesh::recombineIntoQuads(const double angle_limit, GFace *gf)
}
}
-void FullQuadMesh()
-{
-}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 730f5f2b02d6169510fab47630b65905fabe3cdb..96c51568f57e4e11f0feacdb359d1f5f03f5cf23 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -17,6 +17,7 @@
#include "BackgroundMesh.h"
#include "BoundaryLayers.h"
#include "HighOrder.h"
+#include "Generator.h"
#if !defined(HAVE_NO_POST)
#include "PView.h"
@@ -406,6 +407,15 @@ static void Mesh2D(GModel *m)
if(nIter++ > 10) break;
}
}
+
+ // look if there a re model faces for which
+ // full quad algo is set ON
+ bool fullQuad = false;
+ for(GModel::fiter it = m->firstFace() ; it!=m->lastFace(); ++it)
+ if ( CTX.mesh.algo_recombine == 2 && (*it)->quadrangles.size())
+ fullQuad = true;
+ if (fullQuad)RefineMesh(m,false,true);
+
// gmshCollapseSmallEdges (*m);
diff --git a/Mesh/Generator.h b/Mesh/Generator.h
index 182b7032d88ea6d88ebffd5ea465ccd6d25d1060..b9eb74a9beb958c003d23772ea3233b71116114a 100644
--- a/Mesh/Generator.h
+++ b/Mesh/Generator.h
@@ -13,6 +13,6 @@ void AdaptMesh(GModel *m);
void GenerateMesh(GModel *m, int dimension);
void OptimizeMesh(GModel *m);
void OptimizeMeshNetgen(GModel *m);
-void RefineMesh(GModel *m, bool linear=true);
+void RefineMesh(GModel *m, bool linear=true, bool splitTrianglesIntoQuads = false);
#endif
diff --git a/Mesh/Refine.cpp b/Mesh/Refine.cpp
index db40254384690f533447abf3b55bfd24858f3049..2346b4258012d5115af9e2d49de68a82a78f5b47 100644
--- a/Mesh/Refine.cpp
+++ b/Mesh/Refine.cpp
@@ -44,25 +44,28 @@ static void Subdivide(GEdge *ge)
ge->deleteVertexArrays();
}
-static void Subdivide(GFace *gf)
+static void Subdivide(GFace *gf,bool splitTrianglesIntoQuads)
{
- std::vector<MTriangle*> triangles2;
- for(unsigned int i = 0; i < gf->triangles.size(); i++){
- MTriangle *t = gf->triangles[i];
- if(t->getNumVertices() == 6){
- triangles2.push_back
- (new MTriangle(t->getVertex(0), t->getVertex(3), t->getVertex(5)));
- triangles2.push_back
- (new MTriangle(t->getVertex(3), t->getVertex(4), t->getVertex(5)));
- triangles2.push_back
- (new MTriangle(t->getVertex(3), t->getVertex(1), t->getVertex(4)));
- triangles2.push_back
- (new MTriangle(t->getVertex(5), t->getVertex(4), t->getVertex(2)));
+
+ if (! splitTrianglesIntoQuads ){
+ std::vector<MTriangle*> triangles2;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ if(t->getNumVertices() == 6){
+ triangles2.push_back
+ (new MTriangle(t->getVertex(0), t->getVertex(3), t->getVertex(5)));
+ triangles2.push_back
+ (new MTriangle(t->getVertex(3), t->getVertex(4), t->getVertex(5)));
+ triangles2.push_back
+ (new MTriangle(t->getVertex(3), t->getVertex(1), t->getVertex(4)));
+ triangles2.push_back
+ (new MTriangle(t->getVertex(5), t->getVertex(4), t->getVertex(2)));
+ }
+ delete t;
}
- delete t;
+ gf->triangles = triangles2;
}
- gf->triangles = triangles2;
-
+
std::vector<MQuadrangle*> quadrangles2;
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
@@ -78,6 +81,31 @@ static void Subdivide(GFace *gf)
}
delete q;
}
+ if (splitTrianglesIntoQuads){
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ if(t->getNumVertices() == 6){
+ SPoint2 pt, temp;
+ bool reparamOK = true;
+ for(int k = 0; k<6; k++){
+ reparamMeshVertexOnFace(t->getVertex(k), gf, temp);
+ pt[0] += temp[0]/6.;
+ pt[1] += temp[1]/6.;
+ }
+ GPoint gp = gf->point(pt);
+ MFaceVertex *newv = new MFaceVertex (gp.x(),gp.y(),gp.z(),gf,pt[0],pt[1]);
+ gf->mesh_vertices.push_back(newv);
+ quadrangles2.push_back
+ (new MQuadrangle(t->getVertex(0), t->getVertex(3), newv, t->getVertex(5)));
+ quadrangles2.push_back
+ (new MQuadrangle(t->getVertex(3), t->getVertex(1), t->getVertex(4), newv));
+ quadrangles2.push_back
+ (new MQuadrangle(t->getVertex(5), newv,t->getVertex(4), t->getVertex(2)));
+ delete t;
+ }
+ }
+ gf->triangles.clear();
+ }
gf->quadrangles = quadrangles2;
for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
@@ -228,7 +256,7 @@ static void Subdivide(GRegion *gr)
gr->deleteVertexArrays();
}
-void RefineMesh(GModel *m, bool linear)
+void RefineMesh(GModel *m, bool linear, bool splitTrianglesIntoQuads)
{
Msg::StatusBar(1, true, "Refining mesh...");
double t1 = Cpu();
@@ -242,7 +270,7 @@ void RefineMesh(GModel *m, bool linear)
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
Subdivide(*it);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- Subdivide(*it);
+ Subdivide(*it,splitTrianglesIntoQuads);
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
Subdivide(*it);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index f6eb2c1a32b950550cfa56435e9ee157341d33a5..712fd40eaab0aed17552c98e8b0d8354fd7d8e50 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -752,6 +752,9 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
if(AlgoDelaunay2D(gf)){
gmshBowyerWatson(gf);
for (int i=0;i<CTX.mesh.nb_smoothing;i++)laplaceSmoothing(gf);
+ if (gf->meshAttributes.recombine){
+ gmshRecombineIntoQuads(gf);
+ }
}
else if (debug){
char name[256];
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 1287e2ced915e59556d0b7edce0a33f4567e9137..e8bb0bcbc831369ea0871132fa8a39647db0fc8c 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -12,6 +12,7 @@
#include "MElement.h"
#include "BackgroundMesh.h"
#include "GmshMessage.h"
+#include "Generator.h"
static void setLcsInit(MTriangle *t, std::map<MVertex*, double> &vSizes)
{
@@ -126,17 +127,17 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris)
AllTris.erase(AllTris.begin());
}
}
-
-void buildVertexToTriangle(std::vector<MTriangle*> &triangles, v2t_cont &adj)
+template <class T>
+void buildVertexToElement(std::vector<T*> &eles,
+ v2t_cont &adj)
{
- adj.clear();
- for (unsigned int i = 0; i < triangles.size(); i++){
- MTriangle *t = triangles[i];
- for (unsigned int j = 0; j < 3; j++){
+ for (unsigned int i = 0; i < eles.size(); i++){
+ T *t = eles[i];
+ for (unsigned int j = 0; j < t->getNumVertices(); j++){
MVertex *v = t->getVertex(j);
v2t_cont :: iterator it = adj.find(v);
if (it == adj.end()){
- std::vector<MTriangle*> one;
+ std::vector<MElement*> one;
one.push_back(t);
adj[v] = one;
}
@@ -147,23 +148,23 @@ void buildVertexToTriangle(std::vector<MTriangle*> &triangles, v2t_cont &adj)
}
}
-void buildEdgeToTriangle(GFace *gf, e2t_cont &adj)
-{
- buildEdgeToTriangle(gf->triangles, adj);
+
+void buildVertexToTriangle(std::vector<MTriangle*> &eles, v2t_cont &adj){
+ adj.clear();
+ buildVertexToElement(eles,adj);
}
-void buildEdgeToTriangle(std::vector<MTriangle*> &triangles, e2t_cont &adj)
+
+template <class T>
+void buildEdgeToElement(std::vector<T*> &triangles, e2t_cont &adj)
{
- adj.clear();
for (unsigned int i = 0; i < triangles.size(); i++){
- MTriangle *t = triangles[i];
- for (unsigned int j = 0; j < 3; j++){
- MVertex *v1 = t->getVertex(j);
- MVertex *v2 = t->getVertex((j + 1) % 3);
- MEdge e(v1, v2);
+ T *t = triangles[i];
+ for (unsigned int j = 0; j < t->getNumEdges(); j++){
+ MEdge e = t->getEdge(j);
e2t_cont::iterator it = adj.find(e);
if (it == adj.end()){
- std::pair<MTriangle*, MTriangle*> one = std::make_pair(t, (MTriangle*)0);
+ std::pair<MElement*, MElement*> one = std::make_pair(t, (MElement*)0);
adj[e] = one;
}
else
@@ -174,9 +175,23 @@ void buildEdgeToTriangle(std::vector<MTriangle*> &triangles, e2t_cont &adj)
}
}
-void parametricCoordinates(MTriangle *t, GFace *gf, double u[3], double v[3])
+void buildEdgeToElement(GFace *gf, e2t_cont &adj)
{
- for (unsigned int j = 0; j < 3; j++){
+ adj.clear();
+ buildEdgeToElement(gf->triangles, adj);
+ buildEdgeToElement(gf->quadrangles, adj);
+}
+
+void buildEdgeToTriangle(std::vector<MTriangle*> &tris, e2t_cont &adj){
+ adj.clear();
+ buildEdgeToElement(tris, adj);
+}
+
+
+
+void parametricCoordinates(MElement *t, GFace *gf, double u[4], double v[4])
+{
+ for (unsigned int j = 0; j < t->getNumVertices(); j++){
MVertex *ver = t->getVertex(j);
SPoint2 param;
reparamMeshVertexOnFace(ver, gf, param);
@@ -188,7 +203,8 @@ void parametricCoordinates(MTriangle *t, GFace *gf, double u[3], double v[3])
void laplaceSmoothing(GFace *gf)
{
v2t_cont adj;
- buildVertexToTriangle(gf->triangles, adj);
+ buildVertexToElement(gf->triangles, adj);
+ buildVertexToElement(gf->quadrangles, adj);
for (int i = 0; i < 5; i++){
v2t_cont :: iterator it = adj.begin();
@@ -200,22 +216,29 @@ void laplaceSmoothing(GFace *gf)
double initu,initv;
ver->getParameter(0, initu);
ver->getParameter(1, initv);
- const std::vector<MTriangle*> < = it->second;
- double fact = lt.size() ? 1. / (3. * lt.size()) : 0;
+ const std::vector<MElement*> < = it->second;
double cu = 0, cv = 0;
- double pu[3], pv[3];
+ double pu[4], pv[4];
+ double fact = 0.0;
for (unsigned int i = 0; i < lt.size(); i++){
parametricCoordinates(lt[i], gf, pu, pv);
- cu += fact * (pu[0] + pu[1] + pu[2]);
- cv += fact * (pv[0] + pv[1] + pv[2]);
+ cu += (pu[0] + pu[1] + pu[2]);
+ cv += (pv[0] + pv[1] + pv[2]);
+ if (lt[i]->getNumVertices() == 4){
+ cu += pu[3];
+ cv += pv[3];
+ }
+ fact += lt[i]->getNumVertices();
// have to test validity !
}
- ver->setParameter(0, cu);
- ver->setParameter(1, cv);
- GPoint pt = gf->point(SPoint2(cu, cv));
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
+ if (fact != 0.0){
+ ver->setParameter(0, cu/fact);
+ ver->setParameter(1, cv/fact);
+ GPoint pt = gf->point(SPoint2(cu/fact, cv/fact));
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
+ }
}
++it;
}
@@ -772,4 +795,102 @@ int edgeCollapsePass(double minLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &a
return nbCollapse;
}
+extern double angle3Points ( MVertex *p1, MVertex *p2, MVertex *p3 );
+
+struct recombine_triangle
+{
+ MElement *t1, *t2;
+ double angle;
+ MVertex *n1,*n2,*n3,*n4;
+ recombine_triangle(const MEdge &me, MElement *_t1, MElement *_t2)
+ : t1(_t1),t2(_t2)
+ {
+ n1 = me.getVertex(0);
+ n2 = me.getVertex(1);
+
+ if (t1->getVertex (0) != n1 && t1->getVertex (0) != n2)n3 = t1->getVertex(0);
+ else if (t1->getVertex (1) != n1 && t1->getVertex (1) != n2)n3 = t1->getVertex(1);
+ else if (t1->getVertex (2) != n1 && t1->getVertex (2) != n2)n3 = t1->getVertex(2);
+ if (t2->getVertex (0) != n1 && t2->getVertex (0) != n2)n4 = t2->getVertex(0);
+ else if (t2->getVertex (1) != n1 && t2->getVertex (1) != n2)n4 = t2->getVertex(1);
+ else if (t2->getVertex (2) != n1 && t2->getVertex (2) != n2)n4 = t2->getVertex(2);
+
+ double a1 = 180*angle3Points(n1,n4,n2)/M_PI;
+ double a2 = 180*angle3Points(n4,n2,n3)/M_PI;
+ double a3 = 180*angle3Points(n2,n3,n1)/M_PI;
+ double a4 = 180*angle3Points(n3,n1,n4)/M_PI;
+ // printf("%g %g %g %g\n",a1,a2,a3,a4);
+ angle = fabs(90. - a1);
+ angle = std::max(fabs(90. - a2),angle);
+ angle = std::max(fabs(90. - a3),angle);
+ angle = std::max(fabs(90. - a4),angle);
+ }
+ bool operator < (const recombine_triangle &other) const
+ {
+ return angle < other.angle;
+ }
+};
+
+
+void gmshQuadrilateralizeAlgoBrutal(GFace *gf){
+}
+
+void _gmshRecombineIntoQuads(GFace *gf)
+{
+ e2t_cont adj;
+ std::set<MElement*> _touched;
+ std::set<recombine_triangle> pairs;
+ buildEdgeToElement(gf->triangles, adj);
+ e2t_cont::iterator it = adj.begin();
+ for ( ; it!= adj.end(); ++it){
+ if (it->second.second && it->second.first->getNumVertices() == 3 &&
+ it->second.second->getNumVertices() == 3)
+ pairs.insert(recombine_triangle(it->first,
+ it->second.first,
+ it->second.second));
+ }
+ bool rec = false;
+ std::set<recombine_triangle>::iterator itp = pairs.begin();
+ while(itp != pairs.end()){
+ // recombine if difference between max quad angle and right
+ // angle is smaller than tol
+ if(itp->angle < gf->meshAttributes.recombineAngle){
+ MElement *t1 = itp->t1;
+ MElement *t2 = itp->t2;
+ if (_touched.find(t1) == _touched.end() &&
+ _touched.find(t2) == _touched.end()){
+ _touched.insert(t1);
+ _touched.insert(t2);
+ MQuadrangle *q = new MQuadrangle (itp->n1,itp->n3,itp->n2,itp->n4);
+ gf->quadrangles.push_back(q);
+ }
+ }
+ ++itp;
+ }
+
+ std::vector<MTriangle*> _newt;
+ for ( int i = 0 ; i<gf->triangles.size();i++){
+ if (_touched.find(gf->triangles[i]) == _touched.end()){
+ _newt.push_back(gf->triangles[i]);
+ }
+ else {
+ delete gf->triangles[i];
+ }
+ }
+ gf->triangles = _newt;
+
+}
+
+void gmshRecombineIntoQuads(GFace *gf){
+ _gmshRecombineIntoQuads (gf);
+ laplaceSmoothing(gf);
+ _gmshRecombineIntoQuads (gf);
+ laplaceSmoothing(gf);
+ _gmshRecombineIntoQuads (gf);
+ laplaceSmoothing(gf);
+}
+
+
+
+
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index dd25030ef73eafa0d65e62471c03fc88779ae316..24bbb5947d7040ac92058c7132cef21967214356 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -14,8 +14,8 @@
class GFace;
class MVertex;
-typedef std::map<MVertex*, std::vector<MTriangle*> > v2t_cont;
-typedef std::map<MEdge, std::pair<MTriangle*,MTriangle*>, Less_Edge> e2t_cont;
+typedef std::map<MVertex*, std::vector<MElement*> > v2t_cont;
+typedef std::map<MEdge, std::pair<MElement*,MElement*>, Less_Edge> e2t_cont;
void buildVertexToTriangle(std::vector<MTriangle*> &, v2t_cont &adj);
void buildEdgeToTriangle(std::vector<MTriangle*> &, e2t_cont &adj);
void laplaceSmoothing(GFace *gf);
@@ -47,6 +47,7 @@ void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr
std::vector<double> &Us,
std::vector<double> &Vs);
void transferDataStructure (GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris);
+void gmshRecombineIntoQuads(GFace *gf);
struct swapquad{
int v[4];
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index d228aac8272535604c21de33ea90b5428e4cfd36..8cedf26a0ebcff18c2fe2a35e8e8bd41ad6800b3 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -540,9 +540,12 @@ Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial
}
// check for independence
-
+
+
double det = Vandermonde.determinant();
+ // printf("coucou2 %g\n",det);
+
if (det == 0.0){
Msg::Fatal("Vandermonde matrix has zero determinant!?");
return Double_Matrix(1, 1);
@@ -557,6 +560,8 @@ Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial
coefficient(i, j) = f * cofactor.determinant() / det;
}
}
+ // printf("coucou3 %g\n",det);
+
Vandermonde.set_all(0.);
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index d8a67ed2c8e4be7499a8f4de9c7e41f9a991d1f2..f4e3b37eaf50f1bb518300c3b954f8b81a3901d1 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -35,6 +35,14 @@ class Gmsh_Vector
{
return data[i];
}
+ inline SCALAR operator [] (int i) const
+ {
+ return data[i];
+ }
+ inline SCALAR & operator [] (int i)
+ {
+ return data[i];
+ }
inline double norm()
{
double n = 0.;
@@ -52,7 +60,86 @@ class Gmsh_Matrix
{
private:
int r, c;
+
public:
+ void _back_substitution(int *indx, SCALAR *b)
+ {
+ int i,ii=-1,ip,j;
+ SCALAR sum;
+
+ for(i=0; i< c; i++){
+ ip=indx[i];
+ sum=b[ip];
+ b[ip]=b[i];
+ if(ii != -1)
+ for(j=ii; j <= i-1; j++) sum -= (*this)(i,j)*b[j];
+ else if (sum) ii=i;
+ b[i]=sum;
+ }
+ for(i=c-1; i>=0; i--){
+ sum=b[i];
+ for(j=i+1;j<c;j++) sum -= (*this)(i,j)*b[j];
+ b[i]=sum/(*this)(i,i);
+ }
+ }
+
+
+ bool _lu_decomposition(int *indx , SCALAR & determinant)
+ {
+ if (r != c)
+ Msg::Fatal("Gmsh_Matrix::_lu_decomposition : cannot lu decompose a non-square matrix");
+ int i, imax, j,k;
+ SCALAR big,dum,sum,temp;
+ SCALAR *vv = new SCALAR [c];
+
+ determinant=1.0;
+ for(i=0; i<c; i++){
+ big=0.0;
+ for(j=0;j<c; j++)
+ if((temp=fabs((*this)(i,j))) > big) big=temp;
+ if(big==0.0) {
+ return false;
+ big = 1.e-12;
+ }
+ vv[i]=1.0/big;
+ }
+ for(j=0; j<c;j++){
+ for(i=0;i<j;i++){
+ sum=(*this)(i,j);
+ for(k=0;k<i;k++) sum -= (*this)(i,k)*(*this)(k,j);
+ (*this)(i,j) = sum;
+ }
+ big=0.0;
+ for(i=j; i<c;i++){
+ sum=(*this)(i,j);
+ for(k=0;k<j;k++)
+ sum -= (*this)(i,k)*(*this)(k,j);
+ (*this)(i,j)=sum;
+ if((dum=vv[i]*fabs(sum)) >= big){
+ big=dum;
+ imax=i;
+ }
+ }
+ if(j != imax){
+ for(k=0; k <c; k++){
+ dum=(*this)(imax,k);
+ (*this)(imax,k)=(*this)(j,k);
+ (*this)(j,k) = dum;
+ }
+ determinant = -(determinant);
+ vv[imax]=vv[j];
+ }
+ indx[j]=imax;
+ if( (*this)(j,j) == 0.0) (*this)(j,j) = 1.e-20;
+ if(j !=c){
+ dum=1.0/((*this)(j,j));
+ for(i=j+1;i<c;i++) (*this)(i,j) *= dum;
+ }
+ }
+ delete [] vv;
+ return true;
+ }
+
inline int size1() const { return r; }
inline int size2() const { return c; }
SCALAR *data;
@@ -64,9 +151,18 @@ class Gmsh_Matrix
}
Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other) : r(other.r), c(other.c)
{
- data = new double[r * c];
+ data = new SCALAR[r * c];
memcpy(other);
}
+ Gmsh_Matrix & operator=(const Gmsh_Matrix<SCALAR> &other)
+ {
+ if (this != &other){
+ r = other.r; c=other.c;
+ data = new SCALAR[r * c];
+ memcpy(other);
+ }
+ return *this;
+ }
Gmsh_Matrix() : r(0), c(0), data(0) {}
void memcpy(const Gmsh_Matrix &other)
{
@@ -94,33 +190,70 @@ class Gmsh_Matrix
b.data[i] += (*this)(i, j) * x(j);
}
inline void blas_dgemm(const Gmsh_Matrix<SCALAR> & x, const Gmsh_Matrix<SCALAR> & b,
- const double c_a = 1.0, const double c_b = 1.0)
+ const SCALAR c_a = 1.0, const SCALAR c_b = 1.0)
{
- Msg::Fatal("Gmsh_Matrix::blas_dgemm is not implemented");
+ Gmsh_Matrix<SCALAR> temp (x.size1(),b.size2());
+ temp.mult(x,b);
+ scale(c_b);
+ temp.scale(c_a);
+ add (temp);
}
- inline void set_all(const double &m)
+ inline void set_all(const SCALAR &m)
{
for(int i = 0; i < r * c; i++) data[i] = m;
}
inline void lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
- Msg::Fatal("Gmsh_Matrix::lu_solve is not implemented");
- result.scale(0);
+ SCALAR d;
+ int i,*indx;
+ indx = new int [c];
+ if (!_lu_decomposition(indx,d))
+ Msg::Fatal("Singular matrix in Gmsh_Matrix::_lu_decomposition");
+ for(i=0; i < c; i++) result[i] = rhs[i];
+ _back_substitution(indx,result.data);
+ delete [] indx;
}
Gmsh_Matrix cofactor(int i, int j) const
{
- Msg::Fatal("Gmsh_Matrix::cofactor is not implemented");
- Gmsh_Matrix cof(size1() - 1, size2() - 1);
+ int ni = size1();
+ int nj = size2();
+ Gmsh_Matrix<SCALAR> cof(ni - 1, nj - 1);
+ for (int I=0;I<ni;I++){
+ for (int J=0;J<nj;J++){
+ if (J!=j && I!=i)
+ cof (I<i?I:I-1,J<j?J:J-1) = (*this)(I,J);
+ }
+ }
return cof;
}
- inline void invert()
- {
- Msg::Fatal("Gmsh_Matrix::invert is not implemented");
- }
- double determinant() const
- {
- Msg::Fatal("Gmsh_Matrix::determinant is not implemented");
- return 0.;
+ inline void invert(Gmsh_Matrix& y)
+ {
+ SCALAR d,*col;
+ int i,j,*indx;
+
+ col = new SCALAR [c];
+ indx = new int [c];
+ if (!_lu_decomposition(indx,d))
+ Msg::Fatal("Singular matrix in Gmsh_Matrix::_lu_decomposition");
+ for(j=0;j<c;j++){
+ for(i=0; i < c; i++) col[i] = 0.0;
+ col[j]=1.0;
+ _back_substitution(indx,col);
+ for(i=0;i<c;i++) y(i,j)=col[i];
+ }
+ delete [] col;
+ delete [] indx;
+ }
+ SCALAR determinant() const
+ {
+ Gmsh_Matrix copy = *this;
+ SCALAR factor = 1.0;
+ int *indx = new int [c];
+ if (!copy._lu_decomposition(indx, factor))return 0.0;
+ SCALAR det = factor;
+ for (int i=0;i<c;i++)det *= copy(i,i);
+ delete [] indx;
+ return det;
}
inline Gmsh_Matrix touchSubmatrix(int i0, int ni, int j0, int nj)
{
@@ -128,7 +261,7 @@ class Gmsh_Matrix
Gmsh_Matrix subm(ni, nj);
return subm;
}
- inline void scale(const SCALAR s)
+ inline void scale(const double s)
{
for (int i = 0; i < r * c; ++i) data[i] *= s;
}
diff --git a/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
new file mode 100644
index 0000000000000000000000000000000000000000..c20fb7348ce54bd42faadfdbd0018b7d991495d4
--- /dev/null
+++ b/Numeric/gmshLinearSystemFull.h
@@ -0,0 +1,72 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _GMSH_LINEAR_SYSTEM_FULL_H_
+#define _GMSH_LINEAR_SYSTEM_FULL_H_
+
+// Interface to a linear system with a FULL matrix
+
+#include "GmshMessage.h"
+#include "gmshLinearSystem.h"
+
+#include "GmshMatrix.h"
+
+class gmshLinearSystemFull : public gmshLinearSystem {
+ Double_Matrix *_a;
+ Double_Vector *_b, *_x;
+public :
+ gmshLinearSystemFull () : _a(0), _b(0), _x(0){}
+ virtual bool isAllocated () const {return _a != 0;}
+ virtual void allocate (int _nbRows)
+ {
+ if (_a) delete _a;
+ if (_x) delete _x;
+ if (_b) delete _b;
+ _a = new Double_Matrix(_nbRows,_nbRows);
+ _b = new Double_Vector(_nbRows);
+ _x = new Double_Vector(_nbRows);
+ }
+ virtual ~gmshLinearSystemFull ()
+ {
+ delete _a;
+ delete _b;
+ delete _x;
+ }
+ virtual void addToMatrix (int _row, int _col, double _val)
+ {
+ if (_val != 0.0) (*_a)(_row, _col) += _val;
+ }
+ virtual double getFromMatrix (int _row, int _col) const
+ {
+ return (*_a)(_row, _col);
+ }
+ virtual void addToRightHandSide (int _row, double _val)
+ {
+ if (_val != 0.0) (*_b)(_row)+=_val;
+ }
+ virtual double getFromRightHandSide (int _row) const
+ {
+ return (*_b)(_row);
+ }
+ virtual double getFromSolution (int _row) const
+ {
+ return (*_x)(_row);
+ }
+ virtual void zeroMatrix ()
+ {
+ _a->set_all(0.0);
+ }
+ virtual void zeroRightHandSide ()
+ {
+ for (unsigned int i = 0; i < _b->size(); i++) (*_b)(i) = 0;
+ }
+ virtual int systemSolve ()
+ {
+ _a->lu_solve(*_b,*_x);
+ return 1;
+ }
+};
+
+#endif