diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 13e5142eb461f86036048456083f1eb1a20498d1..9b5d159d24361cf98f8bf903ba2111050bb9116f 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -34,6 +34,7 @@
#undef HAVE_OCC
#undef HAVE_OCC_MESH_CONSTRAINTS
#undef HAVE_OSMESA
+#undef HAVE_TAUCS
#undef HAVE_TETGEN
#undef HAVE_TREE_BROWSER
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 29fb4b56cf243ba16c4b9727b2f60c6449e34d94..1082434cb1ba987c7e72f6f45932f8ab4cbadd6f 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -147,9 +147,8 @@ void GEdgeCompound::orderEdges()
int GEdgeCompound::minimumMeshSegments() const
{
int N = 0;
- for (unsigned int i = 0; i < _compound.size(); i++)
- N +=_compound[i]->minimumMeshSegments();
- return N;
+ for (int i=0;i<_compound.size();i++)N +=_compound[i]->minimumMeshSegments();
+ return 3;
}
int GEdgeCompound::minimumDrawSegments() const
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 101f7fb45a44890f897f5dd9157a30cc08003178..6575cbb00a20c7a0c1414e235b09ed25b4ff91d4 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -964,3 +964,32 @@ double GFace::length(const SPoint2 &pt1, const SPoint2 &pt2, int nbQuadPoints)
}
return L;
}
+
+int GFace::poincareMesh()
+{
+ std::set<MEdge, Less_Edge> es;
+ std::set<MVertex*> vs;
+ for (int i=0; i<getNumMeshElements() ; i++){
+ MElement *e = getMeshElement(i);
+ for (int j=0;j<e->getNumVertices();j++)vs.insert(e->getVertex(j));
+ for (int j=0;j<e->getNumEdges();j++)es.insert(e->getEdge(j));
+ }
+ return vs.size() - es.size() + getNumMeshElements();
+}
+
+int GFace::genusGeom()
+{
+ int G = 0;
+ int nSeams = 0;
+ std::set<GEdge*> single_seams;
+ for (std::list<GEdge*>::iterator it = l_edges.begin();
+ it!=l_edges.end();++it){
+ if ((*it)->isSeam(this)){
+ nSeams++;
+ std::set<GEdge*>::iterator it2 = single_seams.find(*it);
+ if (it2 != single_seams.end())single_seams.erase(it2);
+ else single_seams.insert(*it);
+ }
+ }
+ return nSeams - single_seams.size();
+}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 312536b8b74b0532c411566b1b0a14ccca7ebf73..7bfaeaec820653e15cebc02eb8cb72aab6621532 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -49,13 +49,8 @@ class GFace : public GEntity
mean_plane meanPlane;
std::list<GEdge *> embedded_edges;
std::list<GVertex *> embedded_vertices;
- // given a list of GEdges, the function builds a list of wires,
- // i.e. closed edge loops. the first wire is the one that is the
- // outer contour of the face.
- void resolveWires();
GFaceCompound *compound; // this model edge belongs to a compound
-
public: // this will become protected or private
std::list<GEdgeLoop> edgeLoops;
@@ -117,6 +112,19 @@ class GFace : public GEntity
// retrieve surface params
virtual surface_params getSurfaceParams() const;
+ // compute the genus G of the surface
+ // we have the poincare constant CHI = #V-#E+#F
+ // where #V #E and #F are the number of vertices, edges
+ // and faces of the surface
+ // Then, CHI = 2 G + 2 - B
+ // where B is the number of boundaries (edge loops) of the
+ // surface. This topological constant can be computed using both the
+ // geometry and the mesh. Both approaches should give the same result ;-)
+ // by default, genus is ZERO
+ int poincareMesh ();
+ int genusMesh () {return (poincareMesh() + edgeLoops.size() - 2)/2;}
+ virtual int genusGeom ();
+
// return the point on the face corresponding to the given parameter
virtual GPoint point(double par1, double par2) const = 0;
virtual GPoint point(const SPoint2 &pt) const { return point(pt.x(), pt.y()); }
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 1b9375be288ec841404251b409abd790c97f5b59..f1aaf1a622566a84cf9f079bd46e95084e55b7cb 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -13,6 +13,7 @@
#include "gmshLaplace.h"
#include "gmshConvexCombination.h"
#include "gmshLinearSystemGmm.h"
+#include "gmshLinearSystemCSR.h"
#include "gmshLinearSystemFull.h"
#include "FunctionSpace.h"
#include "GmshDefines.h"
@@ -86,12 +87,21 @@ static void fixEdgeToValue(GEdge *ed, double value, gmshAssembler<double> &myAss
}
}
+bool GFaceCompound::trivial() const {
+ if (_compound.size() == 1 &&
+ (*(_compound.begin()))->getNativeType()==GEntity::OpenCascadeModel &&
+ (*(_compound.begin()))->geomType() != GEntity::DiscreteSurface ){
+ return true;
+ }
+ return false;
+}
+
void GFaceCompound::parametrize() const
{
+ if (trivial())return;
if (!oct){
coordinates.clear();
- //parametrize(ITERD);
parametrize(ITERU);
parametrize(ITERV);
computeNormals();
@@ -101,7 +111,6 @@ void GFaceCompound::parametrize() const
void GFaceCompound::getBoundingEdges()
{
-
for (std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
(*it)->setCompound(this);
@@ -148,7 +157,6 @@ void GFaceCompound::getBoundingEdges()
std::set<GEdge*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end(); ++itf){
l_edges.push_back(*itf);
- //printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
(*itf)->addFace(this);
}
@@ -157,40 +165,35 @@ void GFaceCompound::getBoundingEdges()
std::set<GEdge*>::iterator it = _unique.begin();
GVertex *vB = (*it)->getBeginVertex();
GVertex *vE = (*it)->getEndVertex();
- //printf("boundary add edge=%d \n", (*it)->tag());
+ // printf("boundary add edge=%d (%d,%d)\n", (*it)->tag(),vB->tag(),vE->tag());
_loop.push_back(*it);
_unique.erase(it);
it++;
bool found = false;
- for (int i=0; i<2; i++) {
-
- for (std::set<GEdge*>::iterator it = _unique.begin() ; it != _unique.end(); ++it){
- GVertex *v1 = (*it)->getBeginVertex();
- GVertex *v2 = (*it)->getEndVertex();
+ for (int i=0; i<2; i++) {
+ for (std::set<GEdge*>::iterator itx = _unique.begin() ; itx != _unique.end(); ++itx){
+ GVertex *v1 = (*itx)->getBeginVertex();
+ GVertex *v2 = (*itx)->getEndVertex();
std::set<GEdge*>::iterator itp;
if ( v1 == vE ){
- //printf("boundary add edge=%d \n", (*it)->tag());
- _loop.push_back(*it);
- itp = it;
- it++;
+ _loop.push_back(*itx);
+ itp = itx;
+ itx++;
_unique.erase(itp);
vE = v2;
i = -1;
}
else if ( v2 == vE){
- //printf("boundary add edge=%d \n", (*it)->tag());
- _loop.push_back(*it);
- itp = it;
- it++;
+ _loop.push_back(*itx);
+ itp = itx;
+ itx++;
_unique.erase(itp);
vE = v1;
i=-1;
}
-
- if (it == _unique.end()) break;
-
+ if (itx == _unique.end()) break;
}
if (vB == vE) {
@@ -208,9 +211,7 @@ void GFaceCompound::getBoundingEdges()
if (found == true) break;
}
-
_U0 = _loop;
-
}
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
@@ -250,6 +251,9 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
{
l.clear();
coord.clear();
+
+ // printf("%d lines\n",e.size());
+
std::list<GEdge*>::const_iterator it = e.begin();
std::list<MLine*> temp;
double tot_length = 0;
@@ -273,7 +277,7 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
temp.erase(temp.begin());
while (temp.size()){
- bool found = 0;
+ bool found = false;
for (std::list<MLine*>::iterator itl = temp.begin(); itl != temp.end(); ++itl){
MLine *ll = *itl;
MVertex *v0 = ll->getVertex(0);
@@ -308,6 +312,11 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
{
+ if (trivial()){
+ SPoint2 param;
+ reparamMeshVertexOnFace(v, (*(_compound.begin())),param);
+ return param;
+ }
parametrize() ;
std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
@@ -403,7 +412,9 @@ void GFaceCompound::parametrize(iterationStep step) const
Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
#ifdef HAVE_GMM
- gmshLinearSystemGmm<double> lsys;
+
+ // gmshLinearSystemGmm<double> lsys;
+ gmshLinearSystemCSRGmm<double> lsys;
lsys.setPrec(1.e-15);
if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
#else
@@ -425,6 +436,9 @@ void GFaceCompound::parametrize(iterationStep step) const
Msg::Error("Could not order vertices on boundary");
return;
}
+
+ // printf("%d v on the boundary\n", ordered.size());
+
for (unsigned int i = 0; i < ordered.size(); i++){
MVertex *v = ordered[i];
const double theta = 2 * M_PI * coords[i];
@@ -487,7 +501,7 @@ void GFaceCompound::parametrize(iterationStep step) const
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
- diffusivity.setCurrent(t);
+ // diffusivity.setCurrent(t);
laplace.addToMatrix(myAssembler, t);
}
}
@@ -500,8 +514,6 @@ void GFaceCompound::parametrize(iterationStep step) const
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
- //diffusivity.setCompound((*it)->tag());
- //diffusivity.setCurrent(t);
laplace.addToMatrix(myAssembler, t);
}
}
@@ -575,12 +587,16 @@ void GFaceCompound::computeNormals () const
double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
{
+ if (trivial()){
+ return (*(_compound.begin()))->curvatureMax(param);
+ }
+
parametrize();
double U,V;
GFaceCompoundTriangle *lt;
getTriangle (param.x(),param.y(), <, U,V);
if (!lt){
- printf("oops\n");
+ // printf("oops\n");
return 0.0;
}
if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
@@ -592,6 +608,7 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
double GFaceCompound::curvature(MTriangle *t, double u, double v) const
{
+
SVector3 n1 = _normals[t->getVertex(0)];
SVector3 n2 = _normals[t->getVertex(1)];
SVector3 n3 = _normals[t->getVertex(2)];
@@ -604,6 +621,10 @@ double GFaceCompound::curvature(MTriangle *t, double u, double v) const
GPoint GFaceCompound::point(double par1, double par2) const
{
+ if (trivial()){
+ return (*(_compound.begin()))->point(par1,par2);
+ }
+
parametrize();
double U,V;
double par[2] = {par1,par2};
@@ -615,7 +636,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
gp.setNoSuccess();
return gp;
}
- if (0 && lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
+ if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->point(pv.x(),pv.y());
}
@@ -670,13 +691,15 @@ GPoint GFaceCompound::point(double par1, double par2) const
SPoint3 PP(point.x(), point.y(), point.z());
return GPoint(PP.x(),PP.y(),PP.z(),this,par);
-
}
-
}
Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
{
+ if (trivial()){
+ return (*(_compound.begin()))->firstDer(param);
+ }
+
parametrize();
double U,V,UDU,UDV,VDU,VDV;
GFaceCompoundTriangle *lt;
@@ -849,6 +872,26 @@ void GFaceCompound::buildOct() const
printStuff();
}
+int GFaceCompound::genusGeom ()
+{
+ std::list<GFace*> :: const_iterator it = _compound.begin();
+ std::set<MEdge, Less_Edge> es;
+ std::set<MVertex*> vs;
+ int N = 0;
+ for ( ; it != _compound.end() ; ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ N++;
+ MTriangle *e = (*it)->triangles[i];
+ for (int j=0;j<e->getNumVertices();j++)vs.insert(e->getVertex(j));
+ for (int j=0;j<e->getNumEdges();j++)es.insert(e->getEdge(j));
+ }
+ }
+ int poincare = vs.size() - es.size() + N;// = 2g + 2 - b
+
+ return (poincare - 2 + edgeLoops.size())/2;
+}
+
+
void GFaceCompound::printStuff() const
{
std::list<GFace*> :: const_iterator it = _compound.begin();
@@ -892,11 +935,6 @@ void GFaceCompound::printStuff() const
coordinates.find(t->getVertex(1));
std::map<MVertex*,SPoint3>::const_iterator it2 =
coordinates.find(t->getVertex(2));
- fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
- t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
- t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
- t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- it0->second.x(),it1->second.x(),it2->second.x());
// fprintf(xyzu,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
// t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
// t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
@@ -909,11 +947,11 @@ void GFaceCompound::printStuff() const
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
it0->second.y(),it1->second.y(),it2->second.y());
- /* fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- it0->second.x(),it1->second.x(),it2->second.x());*/
+ it0->second.x(),it1->second.x(),it2->second.x());
double K1 = curvature(t,it0->second.x(),it0->second.y());
double K2 = curvature(t,it1->second.x(),it1->second.y());
double K3 = curvature(t,it2->second.x(),it2->second.y());
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index f1a17fc1069ad4d806a16192e41f8e0f196be304..8a302980728585790979d3b5cacc9ea149baa29e 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -63,6 +63,7 @@ class GFaceCompound : public GFace {
double &_u, double &_v) const;
virtual double curvature(MTriangle *t, double u, double v) const;
void printStuff() const;
+ bool trivial() const ;
public:
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag,
@@ -72,19 +73,17 @@ public:
std::list<GEdge*> &V0,
std::list<GEdge*> &V1);
virtual ~GFaceCompound();
- Range<double> parBounds(int i) const { return Range<double>(0, 1); }
+ Range<double> parBounds(int i) const { return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
virtual GPoint point(double par1, double par2) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
-
-
virtual SPoint2 getCoordinates(MVertex *v) const;
-
virtual bool buildRepresentationCross(){ return false; }
virtual double curvatureMax(const SPoint2 ¶m) const;
+ virtual int genusGeom ();
private:
typeOfIsomorphism _type;
};
diff --git a/Geo/Makefile b/Geo/Makefile
index c0c0a6cdbd0ecb050014bc5dbcfce73fb6659740..f1f17316d3775ea90ddffd0b28049c18b1d7ed0c 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -140,6 +140,7 @@ GFaceCompound${OBJEXT}: GFaceCompound.cpp ../Common/GmshConfig.h GFaceCompound.h
../Numeric/gmshConvexCombination.h ../Numeric/gmshTermOfFormulation.h \
../Numeric/gmshFunction.h ../Numeric/GmshMatrix.h \
../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
+ ../Numeric/gmshLinearSystemCSR.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLinearSystemFull.h ../Numeric/gmshLinearSystem.h \
../Numeric/GmshMatrix.h
GRegionCompound${OBJEXT}: GRegionCompound.cpp ../Common/GmshConfig.h \
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 923e1e87334a50acf7d1e0b1ab4d1da21d3a0239..60cdc2c6b9a79fca2f884bf6d8ec2d7434e6ed6a 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -17,12 +17,23 @@ gmshFace::gmshFace(GModel *m, Surface *face)
{
resetMeshAttributes();
+ std::list<GEdge*> l_wire;
+ GVertex *first = 0;
for(int i = 0; i < List_Nbr(s->Generatrices); i++){
Curve *c;
List_Read(s->Generatrices, i, &c);
GEdge *e = m->getEdgeByTag(abs(c->Num));
-
if(e){
+ GVertex *start = (c->Num > 0) ? e->getBeginVertex() : e->getEndVertex();
+ GVertex *next = (c->Num > 0) ? e->getEndVertex() : e->getBeginVertex();
+ if ( ! first ) first = start;
+ l_wire.push_back(e);
+ if ( next == first ){
+ edgeLoops.push_back(GEdgeLoop(l_wire));
+ l_wire.clear();
+ first = 0;
+ }
+
l_edges.push_back(e);
e->addFace(this);
l_dirs.push_back((c->Num > 0) ? 1 : -1);
@@ -62,6 +73,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
}
}
isSphere = iSRuledSurfaceASphere(s, center, radius);
+
}
double gmshFace::getMetricEigenvalue(const SPoint2 &pt)
diff --git a/Numeric/Makefile b/Numeric/Makefile
index d79c3d228e968d7e8492254b01aac8471aa6a5b0..e3cf870075fa646fa4e5fb5d338361b43dae2627 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -7,7 +7,7 @@ include ../variables
LIB = ../lib/libGmshNumeric${LIBEXT}
-INC = ${DASH}I../Common ${DASH}I../Numeric ${DASH}I../Geo
+INC = ${DASH}I../Common ${DASH}I../Numeric ${DASH}I../Geo ${DASH}I../contrib/gmm
CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
@@ -18,6 +18,7 @@ SRC = Numeric.cpp\
gmshHelmholtz.cpp\
gmshElasticity.cpp\
gmshProjection.cpp\
+ gmshLinearSystemCSR.cpp\
EigSolve.cpp\
FunctionSpace.cpp\
GaussQuadratureLin.cpp\
@@ -149,6 +150,8 @@ gmshProjection${OBJEXT}: gmshProjection.cpp gmshProjection.h \
../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/MEdge.h \
../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
../Numeric/Gauss.h ../Common/Gmsh.h ../Common/GmshMessage.h Numeric.h
+gmshLinearSystemCSR${OBJEXT}: gmshLinearSystemCSR.cpp ../Common/GmshConfig.h \
+ ../Common/GmshMessage.h gmshLinearSystemCSR.h gmshLinearSystem.h
EigSolve${OBJEXT}: EigSolve.cpp
FunctionSpace${OBJEXT}: FunctionSpace.cpp FunctionSpace.h GmshMatrix.h \
../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/GmshDefines.h
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 5a70e6736a1f573150a2cb2159fd80dfecc40446..14004b233ed7b26bd76428626cc9b94501e9e967 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -549,7 +549,6 @@ bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
for (int iter = 0; iter < MAXIT; iter++){
func(x, f, data);
- printf("coucou av break !! \n");
bool isZero = false;
for (int k=0; k<N; k++) {
if (f(k) == 0. ) isZero = true;
@@ -557,7 +556,6 @@ bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
if (isZero == false) break;
}
if (isZero) break;
- printf("coucou ap break !! \n");
//printf("**** fffffff0=%g %g %g %g %g %g %g %g %g\n", f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8));
for (int j = 0; j < N; j++){
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index da5a790994e7623ecddad123617bfa65e0e0fa43..05d4b83a01b1228c2ee33cd74b051fbc17b89fb2 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -6,6 +6,12 @@
#include "gmshElasticity.h"
#include "Numeric.h"
+/*
+Non linear version :
+ -) \sigma_{}
+
+*/
+
void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
@@ -71,7 +77,6 @@ void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
}
BTH.set_all(0.);
BTH.gemm(BT, H);
- // printf("detJ = %22.15E\n",detJ);
m.gemm(BTH, B, weight * detJ, 1.);
}
}
diff --git a/configure b/configure
index 7effc75829b3983370496e4609689cfef400ff00..02e0f7f0e434f7f13939551a8fd6758aa1c6f573 100755
--- a/configure
+++ b/configure
@@ -1285,6 +1285,7 @@ Optional Features:
--enable-osmesa use OSMesa for offscreen rendering (default=no)
--enable-cgns enable CGNS output (default=no)
--enable-occ enable OpenCascade support (default=no)
+ --enable-taucs enable Taucs support (default=no)
--enable-hdf5 enable HDF5 support (default=no)
--enable-med enable MED support (default=yes)
--enable-fm enable support for FourierModel (default=yes)
@@ -1300,6 +1301,7 @@ Optional Features:
Optional Packages:
--with-PACKAGE[=ARG] use PACKAGE [ARG=yes]
--without-PACKAGE do not use PACKAGE (same as --with-PACKAGE=no)
+ --with-taucs-prefix=PFX prefix where TAUCS is installed
--with-fltk-prefix=PFX prefix where FLTK is installed
--with-jpeg-prefix=PFX prefix where the JPEG library and includes are
installed
@@ -1761,6 +1763,12 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu
+# Check whether --with-taucs-prefix was given.
+if test "${with_taucs_prefix+set}" = set; then
+ withval=$with_taucs_prefix; TAUCS_PREFIX=$withval
+fi
+
+
# Check whether --with-fltk-prefix was given.
if test "${with_fltk_prefix+set}" = set; then
withval=$with_fltk_prefix; FLTK_PREFIX=$withval
@@ -1945,6 +1953,11 @@ if test "${enable_occ+set}" = set; then
enableval=$enable_occ;
fi
+# Check whether --enable-taucs was given.
+if test "${enable_taucs+set}" = set; then
+ enableval=$enable_taucs;
+fi
+
# Check whether --enable-hdf5 was given.
if test "${enable_hdf5+set}" = set; then
enableval=$enable_hdf5;
@@ -2013,6 +2026,7 @@ if test "x$enable_minimal" = "xyes"; then
enable_cgns=no;
enable_hdf5=no;
enable_zlib=no;
+ enable_taucs=no;
if test "x$enable_post" != "xyes"; then
enable_post=no;
fi
@@ -5110,6 +5124,54 @@ _ACEOF
fi
fi
+if test "x${TAUCS_PREFIX}" != "x"; then
+ if test "x$enable_taucs" != "xno"; then
+ enable_taucs="yes"
+ fi
+fi
+if test "x$enable_taucs" = "xyes"; then
+ if test "x${TAUCS_PREFIX}" != "x"; then
+ LDFLAGS="-L${TAUCS_PREFIX}/lib ${LDFLAGS}"
+ fi
+ as_ac_File=`echo "ac_cv_file_${TAUCS_PREFIX}/src/taucs.h" | $as_tr_sh`
+{ echo "$as_me:$LINENO: checking for ${TAUCS_PREFIX}/src/taucs.h" >&5
+echo $ECHO_N "checking for ${TAUCS_PREFIX}/src/taucs.h... $ECHO_C" >&6; }
+if { as_var=$as_ac_File; eval "test \"\${$as_var+set}\" = set"; }; then
+ echo $ECHO_N "(cached) $ECHO_C" >&6
+else
+ test "$cross_compiling" = yes &&
+ { { echo "$as_me:$LINENO: error: cannot check for file existence when cross compiling" >&5
+echo "$as_me: error: cannot check for file existence when cross compiling" >&2;}
+ { (exit 1); exit 1; }; }
+if test -r "${TAUCS_PREFIX}/src/taucs.h"; then
+ eval "$as_ac_File=yes"
+else
+ eval "$as_ac_File=no"
+fi
+fi
+ac_res=`eval echo '${'$as_ac_File'}'`
+ { echo "$as_me:$LINENO: result: $ac_res" >&5
+echo "${ECHO_T}$ac_res" >&6; }
+if test `eval echo '${'$as_ac_File'}'` = yes; then
+ TAUCS="yes"
+fi
+
+ if test "x${TAUCS}" = "xyes"; then
+ TAUCS_LIBS="-ltaucs"
+ cat >>confdefs.h <<\_ACEOF
+#define HAVE_TAUCS 1
+_ACEOF
+
+ BO="${BO} Taucs"
+ if test "x${TAUCS_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} ${TAUCS_LIBS}"
+ else
+ GMSH_LIBS="${GMSH_LIBS} -L${TAUCS_PREFIX}/lib ${TAUCS_LIBS}"
+ FLAGS="${FLAGS} -I${TAUCS_PREFIX}/src -I${TAUCS_PREFIX}"
+ fi
+ fi
+fi
+
if test "x${OCC}" = "xyes"; then
if test "x${OCC_MESH_CONTRAINTS_PREFIX}" != "x"; then
as_ac_File=`echo "ac_cv_file_${OCC_MESH_CONTRAINTS_PREFIX}/MeshGmsh_Constrain.hxx" | $as_tr_sh`
diff --git a/configure.in b/configure.in
index fad1a4bd955c1c91cfbece378133142b5b27df90..c3187ba066ae8da9843227d20b9d100df6b1c490 100644
--- a/configure.in
+++ b/configure.in
@@ -9,6 +9,10 @@ dnl Check that this is the gmsh source tree
AC_INIT(Geo/GModel.h)
dnl Parse '--with' command-line options
+AC_ARG_WITH(taucs-prefix,
+ AC_HELP_STRING([--with-taucs-prefix=PFX],
+ [prefix where TAUCS is installed]),
+ [TAUCS_PREFIX=$withval])
AC_ARG_WITH(fltk-prefix,
AC_HELP_STRING([--with-fltk-prefix=PFX],
[prefix where FLTK is installed]),
@@ -127,6 +131,9 @@ AC_ARG_ENABLE(cgns,
AC_ARG_ENABLE(occ,
AC_HELP_STRING([--enable-occ],
[enable OpenCascade support (default=no)]))
+AC_ARG_ENABLE(taucs,
+ AC_HELP_STRING([--enable-taucs],
+ [enable Taucs support (default=no)]))
AC_ARG_ENABLE(hdf5,
AC_HELP_STRING([--enable-hdf5],
[enable HDF5 support (default=no)]))
@@ -171,6 +178,7 @@ if test "x$enable_minimal" = "xyes"; then
enable_cgns=no;
enable_hdf5=no;
enable_zlib=no;
+ enable_taucs=no;
if test "x$enable_post" != "xyes"; then
enable_post=no;
fi
@@ -589,6 +597,30 @@ if test "x$enable_occ" = "xyes"; then
fi
fi
+dnl Check for Taucs
+if test "x${TAUCS_PREFIX}" != "x"; then
+ if test "x$enable_taucs" != "xno"; then
+ enable_taucs="yes"
+ fi
+fi
+if test "x$enable_taucs" = "xyes"; then
+ if test "x${TAUCS_PREFIX}" != "x"; then
+ LDFLAGS="-L${TAUCS_PREFIX}/lib ${LDFLAGS}"
+ fi
+ AC_CHECK_FILE(${TAUCS_PREFIX}/src/taucs.h,TAUCS="yes")
+ if test "x${TAUCS}" = "xyes"; then
+ TAUCS_LIBS="-ltaucs"
+ AC_DEFINE(HAVE_TAUCS)
+ BO="${BO} Taucs"
+ if test "x${TAUCS_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} ${TAUCS_LIBS}"
+ else
+ GMSH_LIBS="${GMSH_LIBS} -L${TAUCS_PREFIX}/lib ${TAUCS_LIBS}"
+ FLAGS="${FLAGS} -I${TAUCS_PREFIX}/src -I${TAUCS_PREFIX}"
+ fi
+ fi
+fi
+
dnl Check for OpenCascade mesh constraints
if test "x${OCC}" = "xyes"; then
if test "x${OCC_MESH_CONTRAINTS_PREFIX}" != "x"; then