diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 5bf38803ce60b4caa38fa9f21b3914090d0a2ce2..9748f815cb31785141932c54d9a5e5bb5401b288 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -108,7 +108,7 @@ struct p2data{
gmshHighOrderSmoother *_s)
: gf(_gf), t1(_t1), t2(_t2), n12(_n12), s(_s)
{
- gsolver::elasticityTerm el(0, 1.e3, .3333,1);
+ elasticityTerm el(0, 1.e3, .3333,1);
s->moveToStraightSidedLocation(t1);
s->moveToStraightSidedLocation(t2);
m1 = new gmshMatrix<double>(3 * t1->getNumVertices(),
@@ -136,7 +136,7 @@ struct pNdata{
gmshHighOrderSmoother *_s)
: gf(_gf), t1(_t1), t2(_t2), n(_n), s(_s)
{
- gsolver::elasticityTerm el(0, 1.e3, .3333,1);
+ elasticityTerm el(0, 1.e3, .3333,1);
s->moveToStraightSidedLocation(t1);
s->moveToStraightSidedLocation(t2);
m1 = new gmshMatrix<double>(3 * t1->getNumVertices(),
@@ -429,15 +429,15 @@ void gmshHighOrderSmoother::computeMetricVector(GFace *gf,
void gmshHighOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
{
#ifdef HAVE_TAUCS__
- gsolver::linearSystemCSRTaucs<double> *lsys = new gsolver::linearSystemCSRTaucs<double>;
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
- gsolver::linearSystemCSRGmm<double> *lsys = new gsolver::linearSystemCSRGmm<double>;
+ linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
- gsolver::dofManager<double,double> myAssembler(lsys);
- gsolver::elasticityTerm El(0, 1.0, .333, getTag());
+ dofManager<double,double> myAssembler(lsys);
+ elasticityTerm El(0, 1.0, .333, getTag());
std::vector<MElement*> layer, v;
@@ -542,9 +542,9 @@ void gmshHighOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *
double gmshHighOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
GFace *gf,
- gsolver::dofManager<double,double> &myAssembler,
+ dofManager<double,double> &myAssembler,
std::set<MVertex*> &verticesToMove,
- gsolver::elasticityTerm &El)
+ elasticityTerm &El)
{
std::set<MVertex*>::iterator it;
@@ -572,12 +572,12 @@ double gmshHighOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
for (int j = 0; j < n2; j++){
MVertex *vR;
int iCompR, iFieldR;
- gsolver::Dof RDOF = El.getLocalDofR(e, j);
+ Dof RDOF = El.getLocalDofR(e, j);
myAssembler.assemble(RDOF,-R2(j));
for (int k = 0; k < n2; k++){
MVertex *vC;
int iCompC, iFieldC;
- gsolver::Dof CDOF = El.getLocalDofC(e, k);
+ Dof CDOF = El.getLocalDofC(e, k);
myAssembler.assemble(RDOF,CDOF,K22(j, k));
}
}
@@ -607,15 +607,15 @@ double gmshHighOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
void gmshHighOrderSmoother::smooth(std::vector<MElement*> &all)
{
#ifdef HAVE_TAUCS
- gsolver::linearSystemCSRTaucs<double> *lsys = new gsolver::linearSystemCSRTaucs<double>;
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
- gsolver::linearSystemCSRGmm<double> *lsys = new gsolver::linearSystemCSRGmm<double>;
+ linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
- gsolver::dofManager<double,double> myAssembler(lsys);
- gsolver::elasticityTerm El(0, 1.0, .333, getTag());
+ dofManager<double,double> myAssembler(lsys);
+ elasticityTerm El(0, 1.0, .333, getTag());
std::vector<MElement*> layer, v;
double minD;
diff --git a/Mesh/gmshSmoothHighOrder.h b/Mesh/gmshSmoothHighOrder.h
index 00d549f03cf157f9de0559da5ca97118601cd473..184bd82d0939fe5dab22f455b92724ddff08d6a2 100644
--- a/Mesh/gmshSmoothHighOrder.h
+++ b/Mesh/gmshSmoothHighOrder.h
@@ -24,51 +24,54 @@ class gmshHighOrderSmoother
std::map<MVertex*,SVector3> _straightSidedLocation;
std::map<MVertex*,SVector3> _targetLocation;
int _dim;
- void _clean () {
+ void _clean()
+ {
_straightSidedLocation.clear();
_targetLocation.clear();
}
double _MIDDLE;
- void moveTo (MVertex *v, const std::map<MVertex*,SVector3> &) const;
+ void moveTo(MVertex *v, const std::map<MVertex*,SVector3> &) const;
public:
- gmshHighOrderSmoother (int dim) : _tag(111), _dim(dim) {_clean();}
- void add ( MVertex * v, const SVector3 &d ) {
+ gmshHighOrderSmoother(int dim) : _tag(111), _dim(dim) {_clean();}
+ void add(MVertex * v, const SVector3 &d ) {
_straightSidedLocation[v] = d;
_targetLocation[v] = SPoint3(v->x(),v->y(),v->z());
}
- void smooth ( std::vector<MElement*> & );
- double smooth_metric_ ( std::vector<MElement*> &, GFace *gf,
- gsolver::dofManager<double,double> &myAssembler,
- std::set<MVertex*> &verticesToMove,
- gsolver::elasticityTerm &El);
- void smooth_metric ( std::vector<MElement*> &, GFace *gf );
- void smooth ( GFace* , bool metric = false);
- void smooth_p2point ( GFace* );
- void smooth_pNpoint ( GFace* );
- void smooth ( GRegion* );
- int getTag() const {return _tag;}
+ void smooth(std::vector<MElement*> & );
+ double smooth_metric_(std::vector<MElement*> &, GFace *gf,
+ dofManager<double,double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ elasticityTerm &El);
+ void smooth_metric(std::vector<MElement*> &, GFace *gf );
+ void smooth(GFace *, bool metric = false);
+ void smooth_p2point(GFace *);
+ void smooth_pNpoint(GFace *);
+ void smooth(GRegion*);
+ int getTag() const { return _tag; }
void swap(GFace *,
edgeContainer &edgeVertices,
faceContainer &faceVertices);
void optimize(GFace *,
edgeContainer &edgeVertices,
faceContainer &faceVertices);
- void computeMetricVector (GFace *gf,
- MElement *e,
- gmshMatrix<double> &J,
- gmshMatrix<double> &JT,
- gmshVector<double> &D);
- void moveToStraightSidedLocation (MVertex *) const;
- void moveToTargetLocation (MVertex *) const;
- void moveToStraightSidedLocation (MElement *) const;
- void moveToTargetLocation (MElement *) const;
- void updateTargetLocation (MVertex*, const SPoint3 &, const SPoint2 &) ;
- inline SVector3 getSSL(MVertex *v) const{
+ void computeMetricVector(GFace *gf,
+ MElement *e,
+ gmshMatrix<double> &J,
+ gmshMatrix<double> &JT,
+ gmshVector<double> &D);
+ void moveToStraightSidedLocation(MVertex *) const;
+ void moveToTargetLocation(MVertex *) const;
+ void moveToStraightSidedLocation(MElement *) const;
+ void moveToTargetLocation(MElement *) const;
+ void updateTargetLocation(MVertex*, const SPoint3 &, const SPoint2 &) ;
+ inline SVector3 getSSL(MVertex *v) const
+ {
std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
if (it != _straightSidedLocation.end())return it->second;
else return SVector3(v->x(),v->y(),v->z());
}
- inline SVector3 getDisplacement(MVertex *v) const{
+ inline SVector3 getDisplacement(MVertex *v) const
+ {
std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
std::map<MVertex*,SVector3>::const_iterator itt = _targetLocation.find(v);
if (it == _straightSidedLocation.end())
diff --git a/Numeric/gmshAssembler.h b/Numeric/gmshAssembler.h
index 3a6710ba3bd3384c496351f8639b520a36867766..5bd43b57d8ad0a85ef2123d2213cce98b739a4d5 100644
--- a/Numeric/gmshAssembler.h
+++ b/Numeric/gmshAssembler.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshConvexCombination.cpp b/Numeric/gmshConvexCombination.cpp
index f960e4bed014dc6c36e587bb1749833fb423fc38..814bb2a3189ddcf5bcd1178c2d7495e8b9158150 100644
--- a/Numeric/gmshConvexCombination.cpp
+++ b/Numeric/gmshConvexCombination.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshCrossConf.cpp b/Numeric/gmshCrossConf.cpp
index 33d26de5b078f53b7aa931d675f5fc3adfffb404..938bafbf0c7410df7b7e20033850e52b0a4b06ec 100644
--- a/Numeric/gmshCrossConf.cpp
+++ b/Numeric/gmshCrossConf.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshCrossConf.h b/Numeric/gmshCrossConf.h
index 28ba7cc173daaded4e8941e99b9f4f9d386a4479..801b52ba084697e6c4ba19e9899eb6ae0cc22db0 100644
--- a/Numeric/gmshCrossConf.h
+++ b/Numeric/gmshCrossConf.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshDistance.cpp b/Numeric/gmshDistance.cpp
index 6ce1e563707392800be349455eb14f47827f547c..ec11b25abc954623a8c050f48628ab6db6e8decc 100644
--- a/Numeric/gmshDistance.cpp
+++ b/Numeric/gmshDistance.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshDistance.h b/Numeric/gmshDistance.h
index 18b7720d6136e1076d4d62f000e9bd3024212dcc..588cdaec207a048144ad3ac738ede4b1631f6517 100644
--- a/Numeric/gmshDistance.h
+++ b/Numeric/gmshDistance.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index 05d4b83a01b1228c2ee33cd74b051fbc17b89fb2..f9e1990283d925d2ae93548f69d4218b321deb6a 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index ac4c81268c2c7d0fbda8cc00cd2306d24b0529ab..7c291ac53fa1d5ec0528d27c5422954aff778c20 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshFunction.h b/Numeric/gmshFunction.h
index 04478989ed3b6b410187f0aebf3b8f7c9f8fdfbd..f8d4fa5797621e0a8e5ff6152a8425329508c32f 100644
--- a/Numeric/gmshFunction.h
+++ b/Numeric/gmshFunction.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshHelmholtz.cpp b/Numeric/gmshHelmholtz.cpp
index d78bb53b6ab0fa7552365e841c7d7d2105b3530e..661f40ac7c1d3810c22f7d4f726e56ccde01ff02 100644
--- a/Numeric/gmshHelmholtz.cpp
+++ b/Numeric/gmshHelmholtz.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshHelmholtz.h b/Numeric/gmshHelmholtz.h
index f767f24a31a09ebc3338e8fe320c771a7e9eb5e0..5e08e804d8ac5b173081e63814df32b339c1dc84 100644
--- a/Numeric/gmshHelmholtz.h
+++ b/Numeric/gmshHelmholtz.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 9495e77183ce8cac03e8a7ed8c7717f6f8fc36b7..83fc8daa913be8b31ee7a695c1429bea79b8ffee 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index fb88b0b61af902f75bf2ad243195251b0db9784e..ae63420ddf7cf82c5c17fbc197f0799b4a0a9026 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshLinearSystem.h b/Numeric/gmshLinearSystem.h
index 6e544ac8c399a579640addd542bc9207e361a82a..b4be6df56bf902ce2d049f7d9f7d25ae85deeaa9 100644
--- a/Numeric/gmshLinearSystem.h
+++ b/Numeric/gmshLinearSystem.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshLinearSystemCSR.cpp b/Numeric/gmshLinearSystemCSR.cpp
index a435c03ebc0bd64d20752b4b8e37bda648b39f6e..c2ddd50e6c3a5d91edddab865e611dab3c06d2ad 100644
--- a/Numeric/gmshLinearSystemCSR.cpp
+++ b/Numeric/gmshLinearSystemCSR.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
@@ -13,7 +15,7 @@
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
-void *CSRMalloc(size_t size)
+static void *CSRMalloc(size_t size)
{
void *ptr;
if (!size) return(NULL);
@@ -21,14 +23,14 @@ void *CSRMalloc(size_t size)
return(ptr);
}
-void *CSRRealloc(void *ptr, size_t size)
+static void *CSRRealloc(void *ptr, size_t size)
{
if (!size) return(NULL);
ptr = realloc(ptr,size);
return(ptr);
}
-void CSRList_Realloc(CSRList_T *liste,int n)
+static void CSRList_Realloc(CSRList_T *liste,int n)
{
char* temp;
if (n <= 0) return;
@@ -46,7 +48,7 @@ void CSRList_Realloc(CSRList_T *liste,int n)
}
-CSRList_T *CSRList_Create(int n, int incr, int size)
+static CSRList_T *CSRList_Create(int n, int incr, int size)
{
CSRList_T *liste;
@@ -66,7 +68,7 @@ CSRList_T *CSRList_Create(int n, int incr, int size)
return(liste);
}
-void CSRList_Delete(CSRList_T *liste)
+static void CSRList_Delete(CSRList_T *liste)
{
if (liste != 0) {
free(liste->array);
@@ -74,20 +76,6 @@ void CSRList_Delete(CSRList_T *liste)
}
}
-void CSRList_Add(CSRList_T *liste, void *data)
-{
- liste->n++;
-
- CSRList_Realloc(liste,liste->n);
- liste->isorder = 0;
- memcpy(&liste->array[(liste->n - 1) * liste->size],data,liste->size);
-}
-
-int CSRList_Nbr(CSRList_T *liste)
-{
- return(liste->n);
-}
-
template<>
void gmshLinearSystemCSR<double>::allocate(int _nbRows)
{
@@ -153,7 +141,7 @@ static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj)
}
template <class scalar>
-void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
+static void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
{
unsigned long i,ir=n,j,k,l=1;
int *istack,jstack=0;
@@ -245,12 +233,12 @@ void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[]
}
template <class scalar>
-void sortColumns(int NbLines,
- int nnz,
- INDEX_TYPE *ptr,
- INDEX_TYPE *jptr,
- INDEX_TYPE *ai,
- scalar *a)
+static void sortColumns(int NbLines,
+ int nnz,
+ INDEX_TYPE *ptr,
+ INDEX_TYPE *jptr,
+ INDEX_TYPE *ai,
+ scalar *a)
{
// replace pointers by lines
int *count = new int [NbLines];
diff --git a/Numeric/gmshLinearSystemCSR.h b/Numeric/gmshLinearSystemCSR.h
index 0256182a76c14962af4474584aeb0da9359a045f..3e6b43358437a7830e49e7ee3a92a31635184e76 100644
--- a/Numeric/gmshLinearSystemCSR.h
+++ b/Numeric/gmshLinearSystemCSR.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
@@ -10,20 +12,7 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "gmshLinearSystem.h"
-
-typedef int INDEX_TYPE ;
-
-typedef struct {
- int nmax;
- int size;
- int incr;
- int n;
- int isorder;
- char *array;
-} CSRList_T;
-
-void CSRList_Add(CSRList_T *liste, void *data);
-int CSRList_Nbr(CSRList_T *liste);
+#include "linearSystemCSR.h"
template <class scalar>
class gmshLinearSystemCSR : public gmshLinearSystem<scalar> {
diff --git a/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
index ffc3292a75cc1442e6f6fe0a1982dd6122318ec9..a3c8b129efa72663e8badf93621d99cc2ead45ae 100644
--- a/Numeric/gmshLinearSystemFull.h
+++ b/Numeric/gmshLinearSystemFull.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
index afffe4ce9c85f1ebaa231b9e46d6442ca63fdc6d..372e1d605bb5e2bfa8b21fa467f1e866323d0be1 100644
--- a/Numeric/gmshLinearSystemGmm.h
+++ b/Numeric/gmshLinearSystemGmm.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshProjection.cpp b/Numeric/gmshProjection.cpp
index 231e50b15b3379572b7061eb7ff5919d718b6948..5875c8e668df017e1c02be08e9f1868d0685cb2c 100644
--- a/Numeric/gmshProjection.cpp
+++ b/Numeric/gmshProjection.cpp
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshProjection.h b/Numeric/gmshProjection.h
index 9df8cbb408e3054c3ef36d23f2a00a7e327463be..d47435ae912a0e8eb7491ee4b1e007cfb00239a6 100644
--- a/Numeric/gmshProjection.h
+++ b/Numeric/gmshProjection.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index b2f5d4170c0b8e3b37c76b54a9c0c56d54b298da..5728de1b190602669a89219899e6524541fd0a3a 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -1,3 +1,5 @@
+// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index be55ad186c0ddc6c70ab2e963858008e8846c18a..c85d0f53b9234b0f48a609c8427920ca9d383437 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -12,167 +12,162 @@
#include "MVertex.h"
#include "linearSystem.h"
-namespace gsolver {
-
- class Dof{
- private:
- // v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
- long int _entity; // "i": node, edge, group, etc.
- int _type; // "f": basis function type index, etc.
-
- public:
+class Dof{
+ private:
+ // v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
+ long int _entity; // "i": node, edge, group, etc.
+ int _type; // "f": basis function type index, etc.
+ public:
Dof(long int entity, int type) : _entity(entity), _type(type) {}
- inline long int getEntity() const { return _entity; }
- inline int getType() const { return _type; }
- static int createTypeWithTwoInts( int i1, int i2 ){
- return i1 + 10000 * i2;
- }
- static void getTwoIntsFromType( int t, int &i1, int &i2){
- i1 = t % 10000;
- i2 = t / 10000;
- }
- bool operator < ( const Dof & other) const{
- if (_entity < other._entity)return true;
- if (_entity > other._entity)return false;
- if (_type < other._type)return true;
- return false;
- }
- };
-
+ inline long int getEntity() const { return _entity; }
+ inline int getType() const { return _type; }
+ static int createTypeWithTwoInts(int i1, int i2)
+ {
+ return i1 + 10000 * i2;
+ }
+ static void getTwoIntsFromType(int t, int &i1, int &i2)
+ {
+ i1 = t % 10000;
+ i2 = t / 10000;
+ }
+ bool operator < ( const Dof & other) const
+ {
+ if(_entity < other._entity) return true;
+ if(_entity > other._entity) return false;
+ if(_type < other._type) return true;
+ return false;
+ }
+};
- template<class dataVec, class dataMat>
- class DofAffineConstraint{
- public:
- std::vector<std::pair<Dof, dataMat> > linear;
- dataVec shift;
- };
+template<class dataVec, class dataMat>
+class DofAffineConstraint{
+ public:
+ std::vector<std::pair<Dof, dataMat> > linear;
+ dataVec shift;
+};
- // data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
- template <class dataVec, class dataMat>
- class dofManager{
- private:
- // general affine constraint on sub-blocks, treated by adding
- // equations:
- // dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
- //std::map<std::pair<Dof, dataMat>, DofAffineConstraint> constraints;
-
- // fixations on full blocks, treated by eliminating equations:
- // DofVec = dataVec
- std::map<Dof, dataVec> fixed, initial;
+// data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
+template <class dataVec, class dataMat>
+class dofManager{
+ private:
+ // general affine constraint on sub-blocks, treated by adding
+ // equations:
+ // dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
+ std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec, dataMat> > constraints;
- // numbering of unknown dof blocks
- std::map<Dof, int> unknown;
-
- // associatations
- std::map<Dof, Dof> associatedWith;
+ // fixations on full blocks, treated by eliminating equations:
+ // DofVec = dataVec
+ std::map<Dof, dataVec> fixed;
+
+ // initial conditions
+ std::map<Dof, std::vector<dataVec> > initial;
- // linearSystems
- std::map<const std::string, linearSystem<dataMat>* > _linearSystems;
- linearSystem<dataMat>* _current;
+ // numbering of unknown dof blocks
+ std::map<Dof, int> unknown;
+
+ // associatations
+ std::map<Dof, Dof> associatedWith;
- private:
+ // linearSystems
+ std::map<const std::string, linearSystem<dataMat>* > _linearSystems;
+ linearSystem<dataMat>* _current;
+
public:
- dofManager(linearSystem<dataMat> *l) : _current(l){
- _linearSystems["A"]= l;
- }
- inline void fixDof(long int ent, int type, const dataVec & value)
- {
- fixed [Dof(ent, type)] = value;
- }
- inline void fixVertex(MVertex*v, int iComp, int iField, const dataVec & value)
- {
- fixDof (v->getNum(), Dof::createTypeWithTwoInts(iComp, iField ),value);
- }
- inline void numberDof(long int ent, int type){
- Dof key (ent, type);
- if (fixed.find(key) != fixed.end()) return;
- // if (constraints.find(key) != constraints.end()) return;
+ dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
+ inline void fixDof(long int ent, int type, const dataVec & value)
+ {
+ fixed[Dof(ent, type)] = value;
+ }
+ inline void fixVertex(MVertex*v, int iComp, int iField, const dataVec & value)
+ {
+ fixDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField ), value);
+ }
+ inline void numberDof(long int ent, int type)
+ {
+ Dof key(ent, type);
+ if(fixed.find(key) != fixed.end()) return;
+ // if (constraints.find(key) != constraints.end()) return;
std::map<Dof, int> :: iterator it = unknown.find(key);
if (it == unknown.end()){
unsigned int size = unknown.size();
unknown[key] = size;
}
- }
- inline void numberVertex(MVertex*v, int iComp, int iField){
- numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp,iField ));
- }
- inline dataVec getDofValue(int ent, int type) const
- {
- Dof key (ent, type);
- {
- typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
- if (it != fixed.end()) return it->second;
- }
- {
- std::map<Dof, int>::const_iterator it = unknown.find(key);
- if (it != unknown.end())
- return _current->getFromSolution(it->second);
- }
- return dataVec(0.0);
- }
- inline dataVec getDofValue(MVertex *v, int iComp, int iField) const
+ }
+ inline void numberVertex(MVertex*v, int iComp, int iField)
+ {
+ numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
+ }
+ inline dataVec getDofValue(int ent, int type) const
+ {
+ Dof key(ent, type);
{
- return getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp,iField ));
+ typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
+ if (it != fixed.end()) return it->second;
}
-
- inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
{
- if (!_current->isAllocated()) _current->allocate(unknown.size());
-
- std::map<Dof, int>::iterator itR = unknown.find(R);
- if (itR != unknown.end()){
- std::map<Dof, int>::iterator itC = unknown.find(C);
- if (itC != unknown.end()){
- _current->addToMatrix(itR->second, itC->second, value);
- }
- else{
- typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
- if (itFixed != fixed.end()){
- _current->addToRightHandSide(itR->second, -value*itFixed->second);
- }
- }
- }
+ std::map<Dof, int>::const_iterator it = unknown.find(key);
+ if (it != unknown.end())
+ return _current->getFromSolution(it->second);
}
+ return dataVec(0.0);
+ }
+ inline dataVec getDofValue(MVertex *v, int iComp, int iField) const
+ {
+ return getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
+ }
+ inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
+ {
+ if (!_current->isAllocated()) _current->allocate(unknown.size());
- inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value)
- {
- assemble (Dof(entR,typeR),Dof(entC,typeC), value);
- }
-
- inline void assemble(MVertex *vR, int iCompR, int iFieldR,
- MVertex *vC, int iCompC, int iFieldC,
- const dataMat &value){
- assemble (vR->getNum(), Dof::createTypeWithTwoInts(iCompR,iFieldR),
- vC->getNum(), Dof::createTypeWithTwoInts(iCompC,iFieldC),
- value);
- }
-
- inline void assemble(const Dof &R, const dataMat &value)
- {
- if (!_current->isAllocated())_current->allocate(unknown.size());
- std::map<Dof, int>::iterator itR = unknown.find(R);
- if (itR != unknown.end()){
- _current->addToRightHandSide(itR->second, value);
+ std::map<Dof, int>::iterator itR = unknown.find(R);
+ if (itR != unknown.end()){
+ std::map<Dof, int>::iterator itC = unknown.find(C);
+ if (itC != unknown.end()){
+ _current->addToMatrix(itR->second, itC->second, value);
+ }
+ else{
+ typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
+ if (itFixed != fixed.end()){
+ _current->addToRightHandSide(itR->second, -value * itFixed->second);
+ }
}
}
- inline void assemble(int entR, int typeR, const dataMat &value)
- {
- assemble(Dof(entR,typeR),value);
+ }
+ inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value)
+ {
+ assemble(Dof(entR, typeR), Dof(entC, typeC), value);
+ }
+ inline void assemble(MVertex *vR, int iCompR, int iFieldR,
+ MVertex *vC, int iCompC, int iFieldC,
+ const dataMat &value){
+ assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR),
+ vC->getNum(), Dof::createTypeWithTwoInts(iCompC, iFieldC),
+ value);
+ }
+ inline void assemble(const Dof &R, const dataMat &value)
+ {
+ if(!_current->isAllocated()) _current->allocate(unknown.size());
+ std::map<Dof, int>::iterator itR = unknown.find(R);
+ if(itR != unknown.end()){
+ _current->addToRightHandSide(itR->second, value);
}
- inline void assemble(MVertex *vR, int iCompR, int iFieldR,
- const dataMat &value){
- assemble (vR->getNum(), Dof::createTypeWithTwoInts(iCompR,iFieldR),
- value);
- }
- int sizeOfR() const { return unknown.size(); }
- int sizeOfF() const { return fixed.size(); }
- void systemSolve(){
- _current->systemSolve();
- }
- void systemClear(){
- _current->zeroMatrix();
- _current->zeroRightHandSide();
- }
- };
-}
+ }
+ inline void assemble(int entR, int typeR, const dataMat &value)
+ {
+ assemble(Dof(entR, typeR), value);
+ }
+ inline void assemble(MVertex *vR, int iCompR, int iFieldR,
+ const dataMat &value){
+ assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR), value);
+ }
+ int sizeOfR() const { return unknown.size(); }
+ int sizeOfF() const { return fixed.size(); }
+ void systemSolve(){ _current->systemSolve(); }
+ void systemClear()
+ {
+ _current->zeroMatrix();
+ _current->zeroRightHandSide();
+ }
+};
+
#endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index c81b232895fbb69d271f822a20d38831e442f7f6..e47cef39aaa389ffd4330d61e5d57bb7513151b1 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -1,244 +1,249 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include <string.h>
#include "elasticityTerm.h"
#include "MTriangle.h"
#include "MTetrahedron.h"
#include "elasticitySolver.h"
#include "linearSystemTAUCS.h"
+#include "Numeric.h"
-extern double ComputeVonMises(double *V);
-namespace gsolver {
- double vonMises (dofManager<double,double> *a ,
- MElement *e,
- double u, double v, double w,
- double E, double nu, int _tag){
- double valx[256];
- double valy[256];
- double valz[256];
- for (int k=0;k<e->getNumVertices();k++){
- valx[k]=a->getDofValue(e->getVertex(k),0,_tag);
- valy[k]=a->getDofValue(e->getVertex(k),1,_tag);
- valz[k]=a->getDofValue(e->getVertex(k),2,_tag);
- }
- double gradux[3];
- double graduy[3];
- double graduz[3];
- e->interpolateGrad(valx, u, v, w, gradux);
- e->interpolateGrad(valy, u, v, w, graduy);
- e->interpolateGrad(valz, u, v, w, graduz);
-
- double eps[6] = {gradux[0],graduy[1],graduz[2],
- 0.5*(gradux[1]+graduy[0]),
- 0.5*(gradux[2]+graduz[0]),
- 0.5*(graduy[2]+graduz[1])};
- double A = E/(1.+nu);
- double B = A * (nu/(1.-2*nu));
- double trace = eps[0] + eps[1] + eps[2] ;
- double sxx = A * eps[0] + B * trace;
- double syy = A * eps[1] + B * trace;
- double szz = A * eps[2] + B * trace;
- double sxy = A * eps[3];
- double sxz = A * eps[4];
- double syz = A * eps[5];
-
- // return fabs(sxx);
-
- double s[9] = {sxx,sxy,sxz,
- sxy,syy,syz,
- sxz,syz,szz};
-
- return ComputeVonMises(s);
+static double vonMises(dofManager<double,double> *a, MElement *e,
+ double u, double v, double w,
+ double E, double nu, int _tag)
+{
+ double valx[256];
+ double valy[256];
+ double valz[256];
+ for (int k = 0; k < e->getNumVertices(); k++){
+ valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
+ valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
+ valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
}
+ double gradux[3];
+ double graduy[3];
+ double graduz[3];
+ e->interpolateGrad(valx, u, v, w, gradux);
+ e->interpolateGrad(valy, u, v, w, graduy);
+ e->interpolateGrad(valz, u, v, w, graduz);
- void elasticitySolver::setMesh (const std::string &meshFileName){
- pModel = new GModel();
- pModel->readMSH(meshFileName.c_str());
- _dim = pModel->getNumRegions() ? 3 : 2;
- }
+ double eps[6] = {gradux[0], graduy[1], graduz[2],
+ 0.5 * (gradux[1] + graduy[0]),
+ 0.5 * (gradux[2] + graduz[0]),
+ 0.5 * (graduy[2] + graduz[1])};
+ double A = E / (1. + nu);
+ double B = A * (nu / (1. - 2 * nu));
+ double trace = eps[0] + eps[1] + eps[2] ;
+ double sxx = A * eps[0] + B * trace;
+ double syy = A * eps[1] + B * trace;
+ double szz = A * eps[2] + B * trace;
+ double sxy = A * eps[3];
+ double sxz = A * eps[4];
+ double syz = A * eps[5];
+
+ double s[9] = {sxx, sxy, sxz,
+ sxy, syy, syz,
+ sxz, syz, szz};
- void elasticitySolver::readInputFile (const std::string &fn){
- FILE *f = fopen (fn.c_str(),"r");
- char what[256];
- while(!feof(f)){
- fscanf(f,"%s",what);
- // printf("%s\n",what);
- if (!strcmp(what,"ElasticMaterial")){
- double E,nu;
- int physical;
- fscanf(f,"%d %lf %lf",&physical,&E,&nu);
- elasticConstants[physical] = std::make_pair(E,nu);
- }
- else if (!strcmp(what,"NodalDisplacement")){
- double val;
- int node, comp;
- fscanf(f,"%d %d %lf",&node,&comp,&val);
- nodalDisplacements[ std::make_pair(node,comp) ] = val;
- }
- else if (!strcmp(what,"EdgeDisplacement")){
- double val;
- int edge, comp;
- fscanf(f,"%d %d %lf",&edge,&comp,&val);
- edgeDisplacements[ std::make_pair(edge,comp) ] = val;
- }
- else if (!strcmp(what,"FaceDisplacement")){
- double val;
- int face, comp;
- fscanf(f,"%d %d %lf",&face,&comp,&val);
- faceDisplacements[ std::make_pair(face,comp) ] = val;
- }
- else if (!strcmp(what,"NodalForce")){
- double val1,val2,val3;
- int node;
- fscanf(f,"%d %lf %lf %lf",&node,&val1,&val2,&val3);
- nodalForces[node] = SVector3(val1,val2,val3);
- }
- else if (!strcmp(what,"LineForce")){
- double val1,val2,val3;
- int node;
- fscanf(f,"%d %lf %lf %lf",&node,&val1,&val2,&val3);
- //printf("%d %lf %lf %lf\n",node,val1,val2,val3);
- lineForces[node] = SVector3(val1,val2,val3);
- }
- else if (!strcmp(what,"FaceForce")){
- double val1,val2,val3;
- int face;
- fscanf(f,"%d %lf %lf %lf",&face,&val1,&val2,&val3);
- faceForces[face] = SVector3(val1,val2,val3);
- }
- else if (!strcmp(what,"VolumeForce")){
- double val1,val2,val3;
- int volume;
- fscanf(f,"%d %lf %lf %lf",&volume,&val1,&val2,&val3);
- volumeForces[volume] = SVector3(val1,val2,val3);
- }
- else if (!strcmp(what,"MeshFile")){
- char name[245];
- fscanf(f,"%s",name);
- setMesh(name);
- }
- }
- fclose(f);
- }
+ return ComputeVonMises(s);
+}
- void elasticitySolver::solve (){
-#ifdef HAVE_TAUCS
- linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
-#else
- linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
-#endif
- pAssembler = new dofManager<double,double>(lsys);
-
- std::map<int, std::vector<GEntity*> > groups[4];
- pModel->getPhysicalGroups(groups);
-
-
- // we first do all fixations. the behavior of the dofManager is to
- // give priority to fixations : when a dof is fixed, it cannot be
- // numbered afterwards
-
- for (std::map<std::pair<int,int>, double>::iterator it = nodalDisplacements.begin();
- it != nodalDisplacements.end(); ++it){
- MVertex *v = pModel->getMeshVertexByTag(it->first.first);
- pAssembler-> fixVertex(v , it->first.second, _tag, it->second);
- printf("-- Fixing node %3d comp %1d to %8.5f\n",it->first.first,it->first.second,it->second);
+void elasticitySolver::setMesh(const std::string &meshFileName)
+{
+ pModel = new GModel();
+ pModel->readMSH(meshFileName.c_str());
+ _dim = pModel->getNumRegions() ? 3 : 2;
+}
+
+void elasticitySolver::readInputFile(const std::string &fn)
+{
+ FILE *f = fopen (fn.c_str(),"r");
+ char what[256];
+ while(!feof(f)){
+ fscanf(f, "%s", what);
+ // printf("%s\n", what);
+ if (!strcmp(what,"ElasticMaterial")){
+ double E, nu;
+ int physical;
+ fscanf(f, "%d %lf %lf", &physical, &E, &nu);
+ elasticConstants[physical] = std::make_pair(E, nu);
}
-
- for (std::map<std::pair<int,int>, double>::iterator it = edgeDisplacements.begin();
- it != edgeDisplacements.end(); ++it){
- elasticityTerm El(pModel,1,1,_tag);
- El.dirichletNodalBC(it->first.first,1,it->first.second,_tag,
- gmshFunction<double>(it->second),*pAssembler);
- printf("-- Fixing edge %3d comp %1d to %8.5f\n",it->first.first,it->first.second,it->second);
+ else if (!strcmp(what, "NodalDisplacement")){
+ double val;
+ int node, comp;
+ fscanf(f, "%d %d %lf", &node, &comp, &val);
+ nodalDisplacements[ std::make_pair(node, comp) ] = val;
}
-
- for (std::map<std::pair<int,int>, double>::iterator it = faceDisplacements.begin();
- it != faceDisplacements.end(); ++it){
- elasticityTerm El(pModel,1,1,_tag);
- El.dirichletNodalBC(it->first.first,2,it->first.second,_tag,
- gmshFunction<double>(it->second),*pAssembler);
- printf("-- Fixing face %3d comp %1d to %8.5f\n",it->first.first,it->first.second,it->second);
+ else if (!strcmp(what, "EdgeDisplacement")){
+ double val;
+ int edge, comp;
+ fscanf(f, "%d %d %lf", &edge, &comp, &val);
+ edgeDisplacements[ std::make_pair(edge, comp) ] = val;
}
-
- // we number the nodes : when a node is numbered, it cannot be numbered
- // again with another number. so we loop over elements
- for (std::map<int, std::pair<double,double> >::iterator it = elasticConstants.begin();
- it != elasticConstants.end() ; it++){
- int physical = it->first;
- std::vector<MVertex *> v;
- pModel->getMeshVerticesForPhysicalGroup(_dim,physical,v);
- printf("Physical %d, dim: %d, nb vert: %d\n",physical,_dim,v.size());
- for (int i=0;i<v.size();i++){
- pAssembler->numberVertex (v[i], 0, _tag);
- pAssembler->numberVertex (v[i], 1, _tag);
- pAssembler->numberVertex (v[i], 2, _tag);
- }
+ else if (!strcmp(what, "FaceDisplacement")){
+ double val;
+ int face, comp;
+ fscanf(f, "%d %d %lf", &face, &comp, &val);
+ faceDisplacements[ std::make_pair(face, comp) ] = val;
}
-
- // Now we start the assembly process
- // First build the force vector
- // Nodes at geometrical vertices
- for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
- it != nodalForces.end(); ++it){
- int iVertex = it->first;
- SVector3 f = it->second;
- std::vector<GEntity*> ent = groups[0][iVertex];
- for (int i=0;i<ent.size();i++){
- pAssembler-> assemble(ent[i]->mesh_vertices[0] , 0, _tag, f.x());
- pAssembler-> assemble(ent[i]->mesh_vertices[0] , 1, _tag, f.y());
- pAssembler-> assemble(ent[i]->mesh_vertices[0] , 2, _tag, f.z());
- printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",ent[i]->mesh_vertices[0]->getNum(),iVertex,f.x(),f.y(),f.z());
- }
+ else if (!strcmp(what, "NodalForce")){
+ double val1, val2, val3;
+ int node;
+ fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3);
+ nodalForces[node] = SVector3(val1, val2, val3);
}
- /// line forces
- for (std::map<int, SVector3 >::iterator it = lineForces.begin();
- it != lineForces.end(); ++it){
- elasticityTerm El(pModel,1,1,_tag);
- int iEdge = it->first;
- SVector3 f = it->second;
- El.neumannNodalBC(iEdge,1,0,_tag,gmshFunction<double>(f.x()),*pAssembler);
- El.neumannNodalBC(iEdge,1,1,_tag,gmshFunction<double>(f.y()),*pAssembler);
- El.neumannNodalBC(iEdge,1,2,_tag,gmshFunction<double>(f.z()),*pAssembler);
- printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n",iEdge,f.x(),f.y(),f.z());
+ else if (!strcmp(what, "LineForce")){
+ double val1, val2, val3;
+ int node;
+ fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3);
+ //printf("%d %lf %lf %lf\n", node, val1, val2, val3);
+ lineForces[node] = SVector3(val1, val2, val3);
}
- /// face forces
- for (std::map<int, SVector3 >::iterator it = faceForces.begin();
- it != faceForces.end(); ++it){
- elasticityTerm El(pModel,1,1,_tag);
- int iFace = it->first;
- SVector3 f = it->second;
- El.neumannNodalBC(iFace,2,0,_tag,gmshFunction<double>(f.x()),*pAssembler);
- El.neumannNodalBC(iFace,2,1,_tag,gmshFunction<double>(f.y()),*pAssembler);
- El.neumannNodalBC(iFace,2,2,_tag,gmshFunction<double>(f.z()),*pAssembler);
- printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n",iFace,f.x(),f.y(),f.z());
+ else if (!strcmp(what, "FaceForce")){
+ double val1, val2, val3;
+ int face;
+ fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3);
+ faceForces[face] = SVector3(val1, val2, val3);
}
-
- // volume forces
- for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
- it != volumeForces.end(); ++it){
- elasticityTerm El(pModel,1,1,_tag);
- int iVolume = it->first;
- SVector3 f = it->second;
- El.setVector(f);
- printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n",iVolume,f.x(),f.y(),f.z());
- std::vector<GEntity*> ent = groups[_dim][iVolume];
- for (int i=0;i<ent.size();i++){
- El.addToRightHandSide(*pAssembler,ent[i]);
- }
+ else if (!strcmp(what, "VolumeForce")){
+ double val1, val2, val3;
+ int volume;
+ fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3);
+ volumeForces[volume] = SVector3(val1, val2, val3);
}
-
- // assembling the stifness matrix
- for (std::map<int, std::pair<double,double> > :: iterator it = elasticConstants.begin();
- it != elasticConstants.end() ; it++){
- int physical = it->first;
- double E = it->second.first;
- double nu = it->second.second;
- elasticityTerm El(0,E,nu,_tag);
- std::vector<GEntity*> ent = groups[_dim][physical];
- for (int i=0;i<ent.size();i++){
- El.addToMatrix(*pAssembler,ent[i]);
- }
+ else if (!strcmp(what, "MeshFile")){
+ char name[245];
+ fscanf(f, "%s", name);
+ setMesh(name);
}
- // solving
- lsys->systemSolve();
- }// end of solve
-}// end of namespace
+ }
+ fclose(f);
+}
+
+void elasticitySolver::solve()
+{
+#ifdef HAVE_TAUCS
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#else
+ linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
+#endif
+ pAssembler = new dofManager<double, double>(lsys);
+ std::map<int, std::vector<GEntity*> > groups[4];
+ pModel->getPhysicalGroups(groups);
+
+ // we first do all fixations. the behavior of the dofManager is to
+ // give priority to fixations : when a dof is fixed, it cannot be
+ // numbered afterwards
+
+ for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
+ it != nodalDisplacements.end(); ++it){
+ MVertex *v = pModel->getMeshVertexByTag(it->first.first);
+ pAssembler-> fixVertex(v , it->first.second, _tag, it->second);
+ printf("-- Fixing node %3d comp %1d to %8.5f\n", it->first.first, it->first.second, it->second);
+ }
+
+ for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
+ it != edgeDisplacements.end(); ++it){
+ elasticityTerm El(pModel, 1, 1, _tag);
+ El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
+ gmshFunction<double>(it->second), *pAssembler);
+ printf("-- Fixing edge %3d comp %1d to %8.5f\n", it->first.first, it->first.second, it->second);
+ }
+
+ for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
+ it != faceDisplacements.end(); ++it){
+ elasticityTerm El(pModel, 1, 1, _tag);
+ El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
+ gmshFunction<double>(it->second), *pAssembler);
+ printf("-- Fixing face %3d comp %1d to %8.5f\n", it->first.first, it->first.second, it->second);
+ }
+
+ // we number the nodes : when a node is numbered, it cannot be numbered
+ // again with another number. so we loop over elements
+ for (std::map<int, std::pair<double, double> >::iterator it = elasticConstants.begin();
+ it != elasticConstants.end() ; it++){
+ int physical = it->first;
+ std::vector<MVertex *> v;
+ pModel->getMeshVerticesForPhysicalGroup(_dim, physical, v);
+ printf("Physical %d, dim: %d, nb vert: %d\n", physical, _dim, v.size());
+ for (unsigned int i = 0; i < v.size(); i++){
+ pAssembler->numberVertex (v[i], 0, _tag);
+ pAssembler->numberVertex (v[i], 1, _tag);
+ pAssembler->numberVertex (v[i], 2, _tag);
+ }
+ }
+
+ // Now we start the assembly process
+ // First build the force vector
+ // Nodes at geometrical vertices
+ for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
+ it != nodalForces.end(); ++it){
+ int iVertex = it->first;
+ SVector3 f = it->second;
+ std::vector<GEntity*> ent = groups[0][iVertex];
+ for (unsigned int i = 0; i < ent.size(); i++){
+ pAssembler-> assemble(ent[i]->mesh_vertices[0] , 0, _tag, f.x());
+ pAssembler-> assemble(ent[i]->mesh_vertices[0] , 1, _tag, f.y());
+ pAssembler-> assemble(ent[i]->mesh_vertices[0] , 2, _tag, f.z());
+ printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n", ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
+ }
+ }
+
+ // line forces
+ for (std::map<int, SVector3 >::iterator it = lineForces.begin();
+ it != lineForces.end(); ++it){
+ elasticityTerm El(pModel, 1, 1, _tag);
+ int iEdge = it->first;
+ SVector3 f = it->second;
+ El.neumannNodalBC(iEdge, 1, 0, _tag, gmshFunction<double>(f.x()), *pAssembler);
+ El.neumannNodalBC(iEdge, 1, 1, _tag, gmshFunction<double>(f.y()), *pAssembler);
+ El.neumannNodalBC(iEdge, 1, 2, _tag, gmshFunction<double>(f.z()), *pAssembler);
+ printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
+ }
+
+ // face forces
+ for (std::map<int, SVector3 >::iterator it = faceForces.begin();
+ it != faceForces.end(); ++it){
+ elasticityTerm El(pModel, 1, 1, _tag);
+ int iFace = it->first;
+ SVector3 f = it->second;
+ El.neumannNodalBC(iFace, 2, 0, _tag, gmshFunction<double>(f.x()), *pAssembler);
+ El.neumannNodalBC(iFace, 2, 1, _tag, gmshFunction<double>(f.y()), *pAssembler);
+ El.neumannNodalBC(iFace, 2, 2, _tag, gmshFunction<double>(f.z()), *pAssembler);
+ printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
+ }
+
+ // volume forces
+ for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
+ it != volumeForces.end(); ++it){
+ elasticityTerm El(pModel, 1, 1, _tag);
+ int iVolume = it->first;
+ SVector3 f = it->second;
+ El.setVector(f);
+ printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
+ std::vector<GEntity*> ent = groups[_dim][iVolume];
+ for (unsigned int i = 0; i < ent.size(); i++){
+ El.addToRightHandSide(*pAssembler, ent[i]);
+ }
+ }
+
+ // assembling the stifness matrix
+ for (std::map<int, std::pair<double, double> > :: iterator it = elasticConstants.begin();
+ it != elasticConstants.end() ; it++){
+ int physical = it->first;
+ double E = it->second.first;
+ double nu = it->second.second;
+ elasticityTerm El(0, E, nu, _tag);
+ std::vector<GEntity*> ent = groups[_dim][physical];
+ for (unsigned int i = 0; i < ent.size(); i++){
+ El.addToMatrix(*pAssembler, ent[i]);
+ }
+ }
+
+ // solving
+ lsys->systemSolve();
+}
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 276d3db572a8725df055f7f81de5764c2dc19f35..29a8726d220650e2d426c788c8fc35972aa7086d 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _ELASTICITY_SOLVER_H_
#define _ELASTICITY_SOLVER_H_
@@ -9,46 +14,50 @@
class GModel;
// an elastic solver ...
-namespace gsolver {
- class elasticitySolver
+class elasticitySolver{
+ protected:
+ GModel *pModel;
+ int _dim, _tag;
+ dofManager<double, double> *pAssembler;
+ // young modulus and poisson coefficient per physical
+ std::map<int, std::pair<double, double> > elasticConstants;
+ // imposed nodal forces
+ std::map<int, SVector3> nodalForces;
+ // imposed line forces
+ std::map<int, SVector3> lineForces;
+ // imposed face forces
+ std::map<int, SVector3> faceForces;
+ // imposed volume forces
+ std::map<int, SVector3> volumeForces;
+ // imposed nodal displacements
+ std::map<std::pair<int,int>, double> nodalDisplacements;
+ // imposed edge displacements
+ std::map<std::pair<int,int>, double> edgeDisplacements;
+ // imposed face displacements
+ std::map<std::pair<int,int>, double> faceDisplacements;
+ public:
+ elasticitySolver(int tag) : _tag(tag) {}
+ void addNodalForces (int iNode, const SVector3 &f)
+ {
+ nodalForces[iNode] = f;
+ }
+ void addNodalDisplacement(int iNode, int dir, double val)
+ {
+ nodalDisplacements[std::make_pair(iNode, dir)] = val;
+ }
+ void addElasticConstants(double e, double nu, int physical)
{
- protected:
- GModel *pModel;
- int _dim,_tag;
- dofManager<double,double> *pAssembler;
- // young modulus and poisson coefficient per physical
- std::map<int, std::pair<double,double> > elasticConstants;
- // imposed nodal forces
- std::map<int, SVector3 > nodalForces;
- // imposed line forces
- std::map<int, SVector3 > lineForces;
- // imposed face forces
- std::map<int, SVector3 > faceForces;
- // imposed volume forces
- std::map<int, SVector3 > volumeForces;
- // imposed nodal displacements
- std::map<std::pair<int,int>, double> nodalDisplacements;
- // imposed edge displacements
- std::map<std::pair<int,int>, double> edgeDisplacements;
- // imposed face displacements
- std::map<std::pair<int,int>, double> faceDisplacements;
- public:
- elasticitySolver (int tag) : _tag(tag){}
- void addNodalForces (int iNode, const SVector3 & f )
- {nodalForces[iNode] = f;}
- void addNodalDisplacement (int iNode, int dir, double val )
- {nodalDisplacements[std::make_pair(iNode,dir)] = val;}
- void addElasticConstants (double e, double nu, int physical)
- {elasticConstants[physical] = std::make_pair(e,nu);}
- void setMesh (const std::string &meshFileName);
- virtual void solve ();
- void readInputFile (const std::string &meshFileName);
- // PView * buildDisplacementView (const std::string &postFileName);
- // PView * buildVonMisesView(const std::string &postFileName);
- // std::pair<PView *, PView*> buildErrorEstimateView (const std::string &errorFileName, double, int);
- // std::pair<PView *, PView*> buildErrorEstimateView (const std::string &errorFileName, const gmshElasticityData &ref, double, int);
- }; // end of class definition
-}// end of namespace
+ elasticConstants[physical] = std::make_pair(e, nu);
+ }
+ void setMesh(const std::string &meshFileName);
+ virtual void solve();
+ void readInputFile(const std::string &meshFileName);
+ // PView * buildDisplacementView (const std::string &postFileName);
+ // PView * buildVonMisesView(const std::string &postFileName);
+ // std::pair<PView *, PView*> buildErrorEstimateView (const std::string &errorFileName, double, int);
+ // std::pair<PView *, PView*> buildErrorEstimateView (const std::string &errorFileName, const gmshElasticityData &ref, double, int);
+};
+
#endif
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 1545d4cf00d88348b97be30040de3157ce176d99..06eec04a06f4e916faa1d0ed2664e4aeede9e0dd 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -6,99 +6,98 @@
#include "elasticityTerm.h"
#include "Numeric.h"
-namespace gsolver {
- void elasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
- {
- int nbNodes = e->getNumVertices();
- int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
- int npts;
- IntPt *GP;
- double jac[3][3];
- double invjac[3][3];
- double Grads[256][3], grads[100][3];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.set_all(0.);
-
- double FACT = _E / (1 + _nu);
- double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
- double C12 = FACT * _nu / (1 - 2 * _nu);
- double C44 = (C11 - C12) / 2;
- const double C[6][6] =
- { {C11, C12, C12, 0, 0, 0},
- {C12, C11, C12, 0, 0, 0},
- {C12, C12, C11, 0, 0, 0},
- { 0, 0, 0, C44, 0, 0},
- { 0, 0, 0, 0, C44, 0},
- { 0, 0, 0, 0, 0, C44} };
-
- gmshMatrix<double> H(6, 6);
- gmshMatrix<double> B(6, 3 * nbNodes);
- gmshMatrix<double> BTH(3 * nbNodes, 6);
- gmshMatrix<double> BT(3 * nbNodes, 6);
- for (int i = 0; i < 6; i++)
- for (int j = 0; j < 6; j++)
- H(i, j) = C[i][j];
-
- for (int i = 0; i < npts; i++){
- const double u = GP[i].pt[0];
- const double v = GP[i].pt[1];
- const double w = GP[i].pt[2];
- const double weight = GP[i].weight;
- const double detJ = e->getJacobian(u, v, w, jac);
- e->getGradShapeFunctions(u, v, w, grads);
- inv3x3(jac, invjac);
- B.set_all(0.);
- BT.set_all(0.);
- for (int j = 0; j < nbNodes; j++){
- Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
- invjac[0][2] * grads[j][2];
- Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
- invjac[1][2] * grads[j][2];
- Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
- invjac[2][2] * grads[j][2];
-
- BT(j, 0) = B(0, j) = Grads[j][0];
- BT(j, 3) = B(3, j) = Grads[j][1];
- BT(j, 4) = B(4, j) = Grads[j][2];
-
- BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
- BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
- BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
-
- BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
- BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
- BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
- }
- BTH.set_all(0.);
- BTH.gemm(BT, H);
- m.gemm(BTH, B, weight * detJ, 1.);
- }
- }
+void elasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+{
+ int nbNodes = e->getNumVertices();
+ int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ double invjac[3][3];
+ double Grads[256][3], grads[100][3];
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
- void elasticityTerm::elementVector(MElement *e, gmshVector<double> &m) const{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2 * e->getPolynomialOrder();
- int npts;
- IntPt *GP;
- double jac[3][3];
- double ff[256];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.scale(0.);
-
- for (int i = 0; i < npts; i++){
- const double u = GP[i].pt[0];
- const double v = GP[i].pt[1];
- const double w = GP[i].pt[2];
- const double weight = GP[i].weight;
- const double detJ = e->getJacobian(u, v, w, jac);
- e->getShapeFunctions(u, v, w, ff);
- for (int j = 0; j < nbNodes; j++){
- m(j) += _volumeForce.x() *ff[j] * weight * detJ *.5;
- m(j+nbNodes) += _volumeForce.y() *ff[j] * weight * detJ *.5;
- m(j+2*nbNodes) += _volumeForce.z() *ff[j] * weight * detJ *.5;
- }
- }
- }
+ m.set_all(0.);
+
+ double FACT = _E / (1 + _nu);
+ double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
+ double C12 = FACT * _nu / (1 - 2 * _nu);
+ double C44 = (C11 - C12) / 2;
+ const double C[6][6] =
+ { {C11, C12, C12, 0, 0, 0},
+ {C12, C11, C12, 0, 0, 0},
+ {C12, C12, C11, 0, 0, 0},
+ { 0, 0, 0, C44, 0, 0},
+ { 0, 0, 0, 0, C44, 0},
+ { 0, 0, 0, 0, 0, C44} };
+
+ gmshMatrix<double> H(6, 6);
+ gmshMatrix<double> B(6, 3 * nbNodes);
+ gmshMatrix<double> BTH(3 * nbNodes, 6);
+ gmshMatrix<double> BT(3 * nbNodes, 6);
+ for (int i = 0; i < 6; i++)
+ for (int j = 0; j < 6; j++)
+ H(i, j) = C[i][j];
+
+ for (int i = 0; i < npts; i++){
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
+ const double weight = GP[i].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ e->getGradShapeFunctions(u, v, w, grads);
+ inv3x3(jac, invjac);
+ B.set_all(0.);
+ BT.set_all(0.);
+ for (int j = 0; j < nbNodes; j++){
+ Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2];
+ Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2];
+ Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2];
+
+ BT(j, 0) = B(0, j) = Grads[j][0];
+ BT(j, 3) = B(3, j) = Grads[j][1];
+ BT(j, 4) = B(4, j) = Grads[j][2];
+
+ BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
+ BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
+ BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
+
+ BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
+ BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
+ BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
+ }
+ BTH.set_all(0.);
+ BTH.gemm(BT, H);
+ m.gemm(BTH, B, weight * detJ, 1.);
+ }
+}
+
+void elasticityTerm::elementVector(MElement *e, gmshVector<double> &m) const
+{
+ int nbNodes = e->getNumVertices();
+ int integrationOrder = 2 * e->getPolynomialOrder();
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ double ff[256];
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+
+ m.scale(0.);
+
+ for (int i = 0; i < npts; i++){
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
+ const double weight = GP[i].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ e->getShapeFunctions(u, v, w, ff);
+ for (int j = 0; j < nbNodes; j++){
+ m(j) += _volumeForce.x() *ff[j] * weight * detJ *.5;
+ m(j+nbNodes) += _volumeForce.y() *ff[j] * weight * detJ *.5;
+ m(j+2*nbNodes) += _volumeForce.z() *ff[j] * weight * detJ *.5;
+ }
+ }
}
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index 015a1b4962186518ad82f50204b8cda82e53a4dc..db1c4574586dbe11d421dc7b6912778af36639c2 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -3,41 +3,39 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _ELASTICITY_H_
-#define _ELASTICITY_H_
+#ifndef _ELASTICITY_TERM_H_
+#define _ELASTICITY_TERM_H_
-#include "termOfFormulation.h"
+#include "femTerm.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "GmshMatrix.h"
-namespace gsolver {
- class elasticityTerm : public femTerm<double,double> {
- protected:
- double _E, _nu;
- int _iField;
- SVector3 _volumeForce;
+class elasticityTerm : public femTerm<double,double> {
+ protected:
+ double _E, _nu;
+ int _iField;
+ SVector3 _volumeForce;
public:
- // element matrix size : 3 dofs per vertex
- virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
- // order dofs in the local matrix :
- // dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn
- Dof getLocalDofR(MElement *e, int iRow) const
- {
- int iCompR = iRow / e->getNumVertices();
- int ithLocalVertex = iRow % e->getNumVertices();
- return Dof (e->getVertex(ithLocalVertex)->getNum(),
- Dof::createTypeWithTwoInts(iCompR, _iField));
- }
+ // element matrix size : 3 dofs per vertex
+ virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
+ virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
+ // order dofs in the local matrix :
+ // dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn
+ Dof getLocalDofR(MElement *e, int iRow) const
+ {
+ int iCompR = iRow / e->getNumVertices();
+ int ithLocalVertex = iRow % e->getNumVertices();
+ return Dof(e->getVertex(ithLocalVertex)->getNum(),
+ Dof::createTypeWithTwoInts(iCompR, _iField));
+ }
public:
- elasticityTerm(GModel *gm, double E, double nu, int iField) :
- femTerm<double,double>(gm), _E(E), _nu(nu), _iField(iField){}
- void setVector(const SVector3 &f) {_volumeForce = f;}
- void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
- void elementVector(MElement *e, gmshVector<double> &m) const;
- };
-}
+ elasticityTerm(GModel *gm, double E, double nu, int iField) :
+ femTerm<double,double>(gm), _E(E), _nu(nu), _iField(iField){}
+ void setVector(const SVector3 &f) {_volumeForce = f;}
+ void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+ void elementVector(MElement *e, gmshVector<double> &m) const;
+};
#endif
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..001109e203ef1fa96f1382ee4fd3a6bf303db7f7
--- /dev/null
+++ b/Solver/femTerm.h
@@ -0,0 +1,133 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _FEM_TERM_H_
+#define _FEM_TERM_H_
+
+#include <math.h>
+#include <map>
+#include <vector>
+#include "GmshMatrix.h"
+#include "gmshFunction.h"
+#include "dofManager.h"
+#include "GModel.h"
+#include "MElement.h"
+
+// a nodal finite element term : variables are always defined at nodes
+// of the mesh
+template<class dataVec, class dataMat>
+class femTerm {
+ protected:
+ GModel *_gm;
+ public:
+ // return the number of columns of the element matrix
+ virtual int sizeOfC(MElement*) const = 0;
+ // return the number of rows of the element matrix
+ virtual int sizeOfR(MElement*) const = 0;
+ // in a given element, return the dof associated to a given row (column)
+ // of the local element matrix
+ virtual Dof getLocalDofR(MElement *e, int iRow) const = 0;
+ // default behavior : symmetric
+ virtual Dof getLocalDofC(MElement *e, int iCol) const
+ {
+ return getLocalDofR(e, iCol);
+ }
+ public:
+ femTerm(GModel *gm) : _gm(gm) {}
+ virtual ~femTerm (){}
+ virtual void elementMatrix(MElement *e, gmshMatrix<dataMat> &m) const = 0;
+ virtual void elementVector(MElement *e, gmshVector<dataVec> &m) const
+ {
+ m.scale(0.0);
+ }
+ void addToMatrix(dofManager<dataVec,dataMat> &dm, GEntity *ge) const
+ {
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement(i);
+ addToMatrix(dm, e);
+ }
+ }
+ void addToMatrix(dofManager<dataVec,dataMat> &dm, MElement *e) const
+ {
+ const int nbR = sizeOfR(e);
+ const int nbC = sizeOfC(e);
+ gmshMatrix<dataMat> localMatrix(nbR, nbC);
+ elementMatrix(e, localMatrix);
+ addToMatrix(dm, localMatrix, e);
+ }
+ void addToMatrix(dofManager<dataVec,dataMat> &dm,
+ gmshMatrix<dataMat> &localMatrix,
+ MElement *e) const
+ {
+ const int nbR = localMatrix.size1();
+ const int nbC = localMatrix.size2();
+ for (int j = 0; j < nbR; j++){
+ Dof R = getLocalDofR(e, j);
+ for (int k = 0; k < nbC; k++){
+ Dof C = getLocalDofC(e, k);
+ dm.assemble(R,C, localMatrix(j, k));
+ }
+ }
+ }
+ void dirichletNodalBC(int physical, int dim, int comp, int field,
+ const gmshFunction<dataVec> &e,
+ dofManager<dataVec,dataMat> &dm)
+ {
+ std::vector<MVertex *> v;
+ GModel *m = _gm;
+ m->getMeshVerticesForPhysicalGroup(dim, physical, v);
+ for (unsigned int i = 0; i < v.size(); i++)
+ dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
+ }
+ void neumannNodalBC(int physical, int dim, int comp,int field,
+ const gmshFunction<dataVec> &fct,
+ dofManager<dataVec,dataMat> &dm)
+ {
+ std::map<int, std::vector<GEntity*> > groups[4];
+ GModel *m = _gm;
+ m->getPhysicalGroups(groups);
+ std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
+ if (it == groups[dim].end()) return;
+ double jac[3][3];
+ double sf[256];
+ for (unsigned int i = 0; i < it->second.size(); ++i){
+ GEntity *ge = it->second[i];
+ for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
+ MElement *e = ge->getMeshElement(j);
+ int integrationOrder = 2 * e->getPolynomialOrder();
+ int nbNodes = e->getNumVertices();
+ int npts;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ for (int ip = 0; ip < npts; ip++){
+ const double u = GP[ip].pt[0];
+ const double v = GP[ip].pt[1];
+ const double w = GP[ip].pt[2];
+ const double weight = GP[ip].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ SPoint3 p; e->pnt(u, v, w, p);
+ e->getShapeFunctions(u, v, w, sf);
+ const dataVec FCT = fct(p.x(), p.y(), p.z());
+ for (int k = 0; k < nbNodes; k++){
+ dm.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
+ }
+ }
+ }
+ }
+ }
+ void addToRightHandSide(dofManager<dataVec,dataMat> &dm, GEntity *ge) const
+ {
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement (i);
+ int nbR = sizeOfR(e);
+ gmshVector<dataVec> V (nbR);
+ elementVector (e, V);
+ // assembly
+ for (int j=0;j<nbR;j++)dm.assemble(getLocalDofR(e,j),V(j));
+ }
+ }
+};
+
+#endif
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index c2a6c3179e41ea78aac558f70d45ccb3be38346c..e88cfe3e8e8ad4f6dc2f7fc52747bedc23e255af 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -3,11 +3,11 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _HELMHOLTZ_H_
-#define _HELMHOLTZ_H_
+#ifndef _HELMHOLTZ_TERM_H_
+#define _HELMHOLTZ_TERM_H_
#include <assert.h>
-#include "termOfFormulation.h"
+#include "femTerm.h"
#include "gmshFunction.h"
#include "Gmsh.h"
#include "GModel.h"
@@ -15,80 +15,82 @@
#include "GmshMatrix.h"
#include "Numeric.h"
-namespace gsolver {
- // \nabla \cdot k \nabla U + a U
- template<class scalar>
- class helmoltzTerm : public femTerm<scalar,scalar> {
- protected:
- const gmshFunction<scalar> *_k, *_a;
- const int _iFieldR;
- int _iFieldC ;
+// \nabla \cdot k \nabla U + a U
+template<class scalar>
+class helmoltzTerm : public femTerm<scalar,scalar> {
+ protected:
+ const gmshFunction<scalar> *_k, *_a;
+ const int _iFieldR;
+ int _iFieldC ;
+ public:
+ // one dof per vertex (nodal fem)
+ virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
+ virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
+ Dof getLocalDofR(MElement *e, int iRow) const
+ {
+ return Dof(e->getVertex(iRow)->getNum(), _iFieldR);
+ }
+ Dof getLocalDofC(MElement *e, int iRow) const
+ {
+ return Dof(e->getVertex(iRow)->getNum(), _iFieldC);
+ }
public:
- // one dof per vertex (nodal fem)
- virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
- Dof getLocalDofR(MElement *e, int iRow) const
- {return Dof (e->getVertex(iRow)->getNum(),_iFieldR);}
- Dof getLocalDofC(MElement *e, int iRow) const
- {return Dof (e->getVertex(iRow)->getNum(),_iFieldC);}
- public:
- helmoltzTerm(GModel *gm, int iFieldR, int iFieldC,
- gmshFunction<scalar> *k,
- gmshFunction<scalar> *a) :
- femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR),_iFieldC(iFieldC)
- {}
- virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const{
- // compute integration rule
- const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() : 2 * (e->getPolynomialOrder() - 1);
- int npts;IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- // get the number of nodes
- const int nbNodes = e->getNumVertices();
- // assume a maximum of 100 nodes
- assert (nbNodes < 100);
- double jac[3][3];
- double invjac[3][3];
- double Grads[100][3], grads[100][3];
- double sf[100];
- // set the local matrix to 0
- m.set_all(0.);
- // loop over integration points
- for (int i = 0; i < npts; i++){
- // compute stuff at this point
- const double u = GP[i].pt[0];
- const double v = GP[i].pt[1];
- const double w = GP[i].pt[2];
- const double weightDetJ = GP[i].weight * e->getJacobian(u, v, w, jac);
- SPoint3 p; e->pnt(u, v, w, p);
- const scalar K = (*_k)(p.x(), p.y(), p.z());
- const scalar A = _a ? (*_a)(p.x(), p.y(), p.z()) : 0.0;
- inv3x3(jac, invjac) ;
- e->getGradShapeFunctions(u, v, w, grads);
- if (_a) e->getShapeFunctions(u, v, w, sf);
- for (int j = 0; j < nbNodes; j++){
- Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
- invjac[0][2] * grads[j][2];
- Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
- invjac[1][2] * grads[j][2];
- Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
- invjac[2][2] * grads[j][2];
- if (!_a) sf[j] = 0;
- }
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k <= j; k++){
- m(j, k) += (K* (Grads[j][0] * Grads[k][0] +
- Grads[j][1] * Grads[k][1] +
- Grads[j][2] * Grads[k][2]) + A * sf[j] * sf[k]) *
- weightDetJ;
- }
- }
+ helmoltzTerm(GModel *gm, int iFieldR, int iFieldC,
+ gmshFunction<scalar> *k,
+ gmshFunction<scalar> *a) :
+ femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR), _iFieldC(iFieldC){}
+ virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const
+ {
+ // compute integration rule
+ const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
+ 2 * (e->getPolynomialOrder() - 1);
+ int npts;IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ // get the number of nodes
+ const int nbNodes = e->getNumVertices();
+ // assume a maximum of 100 nodes
+ assert (nbNodes < 100);
+ double jac[3][3];
+ double invjac[3][3];
+ double Grads[100][3], grads[100][3];
+ double sf[100];
+ // set the local matrix to 0
+ m.set_all(0.);
+ // loop over integration points
+ for (int i = 0; i < npts; i++){
+ // compute stuff at this point
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
+ const double weightDetJ = GP[i].weight * e->getJacobian(u, v, w, jac);
+ SPoint3 p; e->pnt(u, v, w, p);
+ const scalar K = (*_k)(p.x(), p.y(), p.z());
+ const scalar A = _a ? (*_a)(p.x(), p.y(), p.z()) : 0.0;
+ inv3x3(jac, invjac) ;
+ e->getGradShapeFunctions(u, v, w, grads);
+ if (_a) e->getShapeFunctions(u, v, w, sf);
+ for (int j = 0; j < nbNodes; j++){
+ Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2];
+ Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2];
+ Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2];
+ if (!_a) sf[j] = 0;
+ }
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k <= j; k++){
+ m(j, k) += (K* (Grads[j][0] * Grads[k][0] +
+ Grads[j][1] * Grads[k][1] +
+ Grads[j][2] * Grads[k][2]) + A * sf[j] * sf[k]) *
+ weightDetJ;
+ }
}
- for (int j = 0; j < nbNodes; j++)
- for (int k = 0; k < j; k++)
- m(k, j) = m(j, k);
-
- }// end of element matrix
- };// end of the class definition
-}// end of namespace
+ }
+ for (int j = 0; j < nbNodes; j++)
+ for (int k = 0; k < j; k++)
+ m(k, j) = m(j, k);
+ }
+};
#endif
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index a1c0394133c62257366a21188f17c1bb891c6666..e402d28c3cbc1a26c69f48c3533c113e1ad65f85 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -3,33 +3,32 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _LAPLACE_H_
-#define _LAPLACE_H_
+#ifndef _LAPLACE_TERM_H_
+#define _LAPLACE_TERM_H_
#include "helmholtzTerm.h"
-namespace gsolver {
- // \nabla \cdot k \nabla U
- template<class scalar>
- class laplaceTerm : public helmholtzTerm<scalar> {
- gmshFunction<scalar> *ONE;
- public:
- laplaceTerm(GModel *gm, int iFieldR) :
- helmholtzTerm<scalar>(gm,iFieldR,iFieldR,0,0)
- {
- ONE = new gmshFunction<scalar>(1.0);
- _k = ONE;
- }
- laplaceTerm(GModel *gm, int iFieldR, gmshFunction<scalar> *k) :
- helmholtzTerm<scalar>(gm,iFieldR,iFieldR,k,0),ONE(0)
- {}
- laplaceTerm(GModel *gm, int iFieldR, int iFieldC, gmshFunction<scalar> *k) :
- helmholtzTerm<scalar>(gm,iFieldR,iFieldC,k,0),ONE(0)
- {}
- virtual ~laplaceTerm(){
- if(ONE) delete ONE;
- }
- };// end of the class definition
-}// end of namespace
+// \nabla \cdot k \nabla U
+template<class scalar>
+class laplaceTerm : public helmholtzTerm<scalar> {
+ gmshFunction<scalar> *ONE;
+ public:
+ laplaceTerm(GModel *gm, int iFieldR) :
+ helmholtzTerm<scalar>(gm, iFieldR, iFieldR, 0, 0)
+ {
+ ONE = new gmshFunction<scalar>(1.0);
+ _k = ONE;
+ }
+ laplaceTerm(GModel *gm, int iFieldR, gmshFunction<scalar> *k) :
+ helmholtzTerm<scalar>(gm, iFieldR, iFieldR, k, 0), ONE(0)
+ {}
+ laplaceTerm(GModel *gm, int iFieldR, int iFieldC, gmshFunction<scalar> *k) :
+ helmholtzTerm<scalar>(gm, iFieldR, iFieldC, k, 0), ONE(0)
+ {}
+ virtual ~laplaceTerm()
+ {
+ if(ONE) delete ONE;
+ }
+};
#endif
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index e274702104443a949fe896bcfffb958c62fa7d98..fe6925ac228585595c33468cfbf7e1a0a8f32bad 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -8,23 +8,22 @@
// A class that encapsulates a linear system solver interface :
// building a sparse matrix, solving a linear system
-namespace gsolver {
- template <class scalar>
- class linearSystem {
- public :
- linearSystem (){}
- virtual ~linearSystem (){}
- 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;
- virtual scalar getFromRightHandSide(int _row) const = 0;
- virtual scalar getFromSolution(int _row) const = 0;
- virtual void zeroMatrix() = 0;
- virtual void zeroRightHandSide() = 0;
- virtual int systemSolve() = 0;
- };
-}
+template <class scalar>
+class linearSystem {
+ public :
+ linearSystem (){}
+ virtual ~linearSystem (){}
+ 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;
+ virtual scalar getFromRightHandSide(int _row) const = 0;
+ virtual scalar getFromSolution(int _row) const = 0;
+ virtual void zeroMatrix() = 0;
+ virtual void zeroRightHandSide() = 0;
+ virtual int systemSolve() = 0;
+};
+
#endif
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 5405999d17bcaa3ae5bdaeb7000c5af9f2b184bf..52b9d8faf942ed26a937a617f389ac1551840ed1 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -12,306 +12,305 @@
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
-namespace gsolver {
- void *CSRMalloc(size_t size)
- {
- void *ptr;
- if (!size) return(NULL);
- ptr = malloc(size);
- return(ptr);
- }
-
- void *CSRRealloc(void *ptr, size_t size)
- {
- if (!size) return(NULL);
- ptr = realloc(ptr,size);
- return(ptr);
+static void *CSRMalloc(size_t size)
+{
+ void *ptr;
+ if (!size) return(NULL);
+ ptr = malloc(size);
+ return(ptr);
+}
+
+static void *CSRRealloc(void *ptr, size_t size)
+{
+ if (!size) return(NULL);
+ ptr = realloc(ptr,size);
+ return(ptr);
+}
+
+static void CSRList_Realloc(CSRList_T *liste,int n)
+{
+ char* temp;
+ if (n <= 0) return;
+ if (liste->array == NULL) {
+ liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+ liste->array = (char *)CSRMalloc(liste->nmax * liste->size);
}
-
- void CSRList_Realloc(CSRList_T *liste,int n)
- {
- char* temp;
- if (n <= 0) return;
- if (liste->array == NULL) {
+ else {
+ if (n > liste->nmax) {
liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
- liste->array = (char *)CSRMalloc(liste->nmax * liste->size);
- }
- else {
- if (n > liste->nmax) {
- liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
- temp = (char *)CSRRealloc(liste->array, liste->nmax * liste->size);
- liste->array = temp;
- }
+ temp = (char *)CSRRealloc(liste->array, liste->nmax * liste->size);
+ liste->array = temp;
}
}
+}
+
+
+static CSRList_T *CSRList_Create(int n, int incr, int size)
+{
+ CSRList_T *liste;
+ if (n <= 0) n = 1 ;
+ if (incr <= 0) incr = 1;
- CSRList_T *CSRList_Create(int n, int incr, int size)
- {
- CSRList_T *liste;
-
- if (n <= 0) n = 1 ;
- if (incr <= 0) incr = 1;
-
- liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
-
- liste->nmax = 0;
- liste->incr = incr;
- liste->size = size;
- liste->n = 0;
- liste->isorder = 0;
- liste->array = NULL;
-
- CSRList_Realloc(liste,n);
- return(liste);
- }
+ liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
- void CSRList_Delete(CSRList_T *liste)
- {
- if (liste != 0) {
- free(liste->array);
- free(liste);
- }
- }
+ liste->nmax = 0;
+ liste->incr = incr;
+ liste->size = size;
+ liste->n = 0;
+ liste->isorder = 0;
+ liste->array = NULL;
- void CSRList_Add(gsolver::CSRList_T *liste, void *data)
- {
- liste->n++;
-
- CSRList_Realloc(liste,liste->n);
- liste->isorder = 0;
- memcpy(&liste->array[(liste->n - 1) * liste->size],data,liste->size);
+ CSRList_Realloc(liste,n);
+ return(liste);
+}
+
+static void CSRList_Delete(CSRList_T *liste)
+{
+ if (liste != 0) {
+ free(liste->array);
+ free(liste);
}
+}
+
+void CSRList_Add(CSRList_T *liste, void *data)
+{
+ liste->n++;
- int CSRList_Nbr(CSRList_T *liste)
- {
- return(liste->n);
+ CSRList_Realloc(liste,liste->n);
+ liste->isorder = 0;
+ memcpy(&liste->array[(liste->n - 1) * liste->size],data,liste->size);
+}
+
+int CSRList_Nbr(CSRList_T *liste)
+{
+ return(liste->n);
+}
+
+template<>
+void linearSystemCSR<double>::allocate(int _nbRows)
+{
+ if(a_) {
+ CSRList_Delete(a_);
+ CSRList_Delete(ai_);
+ CSRList_Delete(ptr_);
+ CSRList_Delete(jptr_);
+ delete _x;
+ delete _b;
+ delete something;
}
- template<>
- void linearSystemCSR<double>::allocate(int _nbRows)
- {
- if(a_) {
- CSRList_Delete(a_);
- CSRList_Delete(ai_);
- CSRList_Delete(ptr_);
- CSRList_Delete(jptr_);
- delete _x;
- delete _b;
- delete something;
- }
-
- if(_nbRows == 0){
- a_ = 0;
- ai_ = 0;
- ptr_ = 0;
- jptr_ = 0;
- _b = 0;
- _x = 0;
- something = 0;
- return;
- }
-
- a_ = CSRList_Create (_nbRows, _nbRows, sizeof(double));
- ai_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
- ptr_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
- jptr_ = CSRList_Create (_nbRows+1, _nbRows, sizeof(INDEX_TYPE));
-
- something = new char[_nbRows];
-
- for (int i=0;i<_nbRows;i++)something[i] = 0;
-
- _b = new std::vector<double>(_nbRows);
- _x = new std::vector<double>(_nbRows);
+ if(_nbRows == 0){
+ a_ = 0;
+ ai_ = 0;
+ ptr_ = 0;
+ jptr_ = 0;
+ _b = 0;
+ _x = 0;
+ something = 0;
+ return;
}
- const int NSTACK = 50;
- const unsigned int M_sort2 = 7;
+ a_ = CSRList_Create (_nbRows, _nbRows, sizeof(double));
+ ai_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
+ ptr_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
+ jptr_ = CSRList_Create (_nbRows+1, _nbRows, sizeof(INDEX_TYPE));
- static void free_ivector(int *v, long nl, long nh)
- {
- // free an int vector allocated with ivector()
- free((char*)(v+nl-1));
- }
+ something = new char[_nbRows];
- static int *ivector(long nl, long nh)
- {
- // allocate an int vector with subscript range v[nl..nh]
- int *v;
- v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
- if (!v) fprintf(stderr, "allocation failure in ivector()\n");
- return v-nl+1;
- }
+ for (int i=0;i<_nbRows;i++)something[i] = 0;
- static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj)
- {
- if(ai<bi)return -1;
- if(ai>bi)return 1;
- if(aj<bj)return -1;
- if(aj>bj)return 1;
- return 0;
- }
+ _b = new std::vector<double>(_nbRows);
+ _x = new std::vector<double>(_nbRows);
+}
+
+const int NSTACK = 50;
+const unsigned int M_sort2 = 7;
+
+static void free_ivector(int *v, long nl, long nh)
+{
+ // free an int vector allocated with ivector()
+ free((char*)(v+nl-1));
+}
+
+static int *ivector(long nl, long nh)
+{
+ // allocate an int vector with subscript range v[nl..nh]
+ int *v;
+ v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
+ if (!v) fprintf(stderr, "allocation failure in ivector()\n");
+ return v-nl+1;
+}
+
+static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj)
+{
+ if(ai<bi)return -1;
+ if(ai>bi)return 1;
+ if(aj<bj)return -1;
+ if(aj>bj)return 1;
+ return 0;
+}
+
+template <class scalar>
+static void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
+{
+ unsigned long i,ir=n,j,k,l=1;
+ int *istack,jstack=0;
+ INDEX_TYPE tempi;
+ scalar a,temp;
+ int b,c;
- template <class scalar>
- void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
- {
- unsigned long i,ir=n,j,k,l=1;
- int *istack,jstack=0;
- INDEX_TYPE tempi;
- scalar a,temp;
- int b,c;
-
- istack=ivector(1,NSTACK);
- for (;;) {
- if (ir-l < M_sort2) {
- for (j=l+1;j<=ir;j++) {
- a=arr[j -1];
- b=ai[j -1];
- c=aj[j -1];
- for (i=j-1;i>=1;i--) {
- if (cmpij(ai[i -1],aj[i -1],b,c) <= 0) break;
- arr[i+1 -1]=arr[i -1];
- ai[i+1 -1]=ai[i -1];
- aj[i+1 -1]=aj[i -1];
- }
- arr[i+1 -1]=a;
- ai[i+1 -1]=b;
- aj[i+1 -1]=c;
- }
- if (!jstack) {
- free_ivector(istack,1,NSTACK);
- return;
- }
- ir=istack[jstack];
- l=istack[jstack-1];
- jstack -= 2;
+ istack=ivector(1,NSTACK);
+ for (;;) {
+ if (ir-l < M_sort2) {
+ for (j=l+1;j<=ir;j++) {
+ a=arr[j -1];
+ b=ai[j -1];
+ c=aj[j -1];
+ for (i=j-1;i>=1;i--) {
+ if (cmpij(ai[i -1],aj[i -1],b,c) <= 0) break;
+ arr[i+1 -1]=arr[i -1];
+ ai[i+1 -1]=ai[i -1];
+ aj[i+1 -1]=aj[i -1];
+ }
+ arr[i+1 -1]=a;
+ ai[i+1 -1]=b;
+ aj[i+1 -1]=c;
+ }
+ if (!jstack) {
+ free_ivector(istack,1,NSTACK);
+ return;
+ }
+ ir=istack[jstack];
+ l=istack[jstack-1];
+ jstack -= 2;
+ }
+ else {
+ k=(l+ir) >> 1;
+ SWAP(arr[k -1],arr[l+1 -1])
+ SWAPI(ai[k -1],ai[l+1 -1])
+ SWAPI(aj[k -1],aj[l+1 -1])
+ if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[ir -1],aj[ir -1])>0){
+ SWAP(arr[l+1 -1],arr[ir -1])
+ SWAPI(ai[l+1 -1],ai[ir -1])
+ SWAPI(aj[l+1 -1],aj[ir -1])
+ }
+ if (cmpij(ai[l -1],aj[l -1],ai[ir -1],aj[ir -1])>0){
+ SWAP(arr[l -1],arr[ir -1])
+ SWAPI(ai[l -1],ai[ir -1])
+ SWAPI(aj[l -1],aj[ir -1])
+ }
+ if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[l -1],aj[l -1])>0){
+ SWAP(arr[l+1 -1],arr[l -1])
+ SWAPI(ai[l+1 -1],ai[l -1])
+ SWAPI(aj[l+1 -1],aj[l -1])
+ }
+ i=l+1;
+ j=ir;
+ a=arr[l -1];
+ b=ai[l -1];
+ c=aj[l -1];
+ for (;;) {
+ do i++; while (cmpij(ai[i -1],aj[i -1],b,c) < 0);
+ do j--; while (cmpij(ai[j -1],aj[j -1],b,c) > 0);
+ if (j < i) break;
+ SWAP(arr[i -1],arr[j -1])
+ SWAPI(ai[i -1],ai[j -1])
+ SWAPI(aj[i -1],aj[j -1])
+ }
+ arr[l -1]=arr[j -1];
+ arr[j -1]=a;
+ ai[l -1]=ai[j -1];
+ ai[j -1]=b;
+ aj[l -1]=aj[j -1];
+ aj[j -1]=c;
+ jstack += 2;
+ if (jstack > NSTACK) {
+ Msg::Fatal("NSTACK too small while sorting the columns of the matrix");
+ throw;
+ }
+ if (ir-i+1 >= j-l) {
+ istack[jstack]=ir;
+ istack[jstack-1]=i;
+ ir=j-1;
}
else {
- k=(l+ir) >> 1;
- SWAP(arr[k -1],arr[l+1 -1])
- SWAPI(ai[k -1],ai[l+1 -1])
- SWAPI(aj[k -1],aj[l+1 -1])
- if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[ir -1],aj[ir -1])>0){
- SWAP(arr[l+1 -1],arr[ir -1])
- SWAPI(ai[l+1 -1],ai[ir -1])
- SWAPI(aj[l+1 -1],aj[ir -1])
- }
- if (cmpij(ai[l -1],aj[l -1],ai[ir -1],aj[ir -1])>0){
- SWAP(arr[l -1],arr[ir -1])
- SWAPI(ai[l -1],ai[ir -1])
- SWAPI(aj[l -1],aj[ir -1])
- }
- if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[l -1],aj[l -1])>0){
- SWAP(arr[l+1 -1],arr[l -1])
- SWAPI(ai[l+1 -1],ai[l -1])
- SWAPI(aj[l+1 -1],aj[l -1])
- }
- i=l+1;
- j=ir;
- a=arr[l -1];
- b=ai[l -1];
- c=aj[l -1];
- for (;;) {
- do i++; while (cmpij(ai[i -1],aj[i -1],b,c) < 0);
- do j--; while (cmpij(ai[j -1],aj[j -1],b,c) > 0);
- if (j < i) break;
- SWAP(arr[i -1],arr[j -1])
- SWAPI(ai[i -1],ai[j -1])
- SWAPI(aj[i -1],aj[j -1])
- }
- arr[l -1]=arr[j -1];
- arr[j -1]=a;
- ai[l -1]=ai[j -1];
- ai[j -1]=b;
- aj[l -1]=aj[j -1];
- aj[j -1]=c;
- jstack += 2;
- if (jstack > NSTACK) {
- Msg::Fatal("NSTACK too small while sorting the columns of the matrix");
- throw;
- }
- if (ir-i+1 >= j-l) {
- istack[jstack]=ir;
- istack[jstack-1]=i;
- ir=j-1;
- }
- else {
- istack[jstack]=j-1;
- istack[jstack-1]=l;
- l=i;
- }
+ istack[jstack]=j-1;
+ istack[jstack-1]=l;
+ l=i;
}
}
}
+}
+
+template <class scalar>
+static void sortColumns(int NbLines,
+ int nnz,
+ INDEX_TYPE *ptr,
+ INDEX_TYPE *jptr,
+ INDEX_TYPE *ai,
+ scalar *a)
+{
+ // replace pointers by lines
+ int *count = new int [NbLines];
- template <class scalar>
- void sortColumns(int NbLines,
- int nnz,
- INDEX_TYPE *ptr,
- INDEX_TYPE *jptr,
- INDEX_TYPE *ai,
- scalar *a)
- {
- // replace pointers by lines
- int *count = new int [NbLines];
-
- for(int i=0;i<NbLines;i++){
- count[i] = 0;
- INDEX_TYPE _position = jptr[i];
- while(1){
- count[i]++;
- INDEX_TYPE _position_temp = _position;
- _position = ptr[_position];
- ptr[_position_temp] = i;
- if (_position == 0) break;
- }
- }
- _sort2_xkws<double>(nnz,a,ptr,ai);
- jptr[0] = 0;
- for(int i=1;i<=NbLines;i++){
- jptr[i] = jptr[i-1] + count[i-1];
+ for(int i=0;i<NbLines;i++){
+ count[i] = 0;
+ INDEX_TYPE _position = jptr[i];
+ while(1){
+ count[i]++;
+ INDEX_TYPE _position_temp = _position;
+ _position = ptr[_position];
+ ptr[_position_temp] = i;
+ if (_position == 0) break;
}
-
- for(int i=0;i<NbLines;i++){
- for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
- ptr[j] = j+1;
- }
- ptr[jptr[i+1]] = 0;
+ }
+ _sort2_xkws<double>(nnz,a,ptr,ai);
+ jptr[0] = 0;
+ for(int i=1;i<=NbLines;i++){
+ jptr[i] = jptr[i-1] + count[i-1];
+ }
+
+ for(int i=0;i<NbLines;i++){
+ for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
+ ptr[j] = j+1;
}
-
- delete[] count;
+ ptr[jptr[i+1]] = 0;
}
-}
+ delete[] count;
+}
+
#if defined(HAVE_GMM)
+
#include "gmm.h"
-namespace gsolver {
- template<>
- int linearSystemCSRGmm<double>::systemSolve()
- {
- if (!sorted)
- sortColumns(_b->size(),
- CSRList_Nbr(a_),
- (INDEX_TYPE *) ptr_->array,
- (INDEX_TYPE *) jptr_->array,
- (INDEX_TYPE *) ai_->array,
- (double*) a_->array);
- sorted = true;
-
- gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>
- ref((double*)a_->array, (INDEX_TYPE *) ai_->array,
- (INDEX_TYPE *)jptr_->array, _b->size(), _b->size());
- gmm::csr_matrix<double,0> M;
- M.init_with(ref);
-
- gmm::ildltt_precond<gmm::csr_matrix<double, 0> > P(M, 10, 1.e-10);
- gmm::iteration iter(_prec);
- iter.set_noisy(_noisy);
- if(_gmres) gmm::gmres(M, *_x, *_b, P, 100, iter);
- else gmm::cg(M, *_x, *_b, P, iter);
- return 1;
- }
+
+template<>
+int linearSystemCSRGmm<double>::systemSolve()
+{
+ if (!sorted)
+ sortColumns(_b->size(),
+ CSRList_Nbr(a_),
+ (INDEX_TYPE *) ptr_->array,
+ (INDEX_TYPE *) jptr_->array,
+ (INDEX_TYPE *) ai_->array,
+ (double*) a_->array);
+ sorted = true;
+
+ gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>
+ ref((double*)a_->array, (INDEX_TYPE *) ai_->array,
+ (INDEX_TYPE *)jptr_->array, _b->size(), _b->size());
+ gmm::csr_matrix<double,0> M;
+ M.init_with(ref);
+
+ gmm::ildltt_precond<gmm::csr_matrix<double, 0> > P(M, 10, 1.e-10);
+ gmm::iteration iter(_prec);
+ iter.set_noisy(_noisy);
+ if(_gmres) gmm::gmres(M, *_x, *_b, P, 100, iter);
+ else gmm::cg(M, *_x, *_b, P, iter);
+ return 1;
}
+
#endif
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 1dcaae4e714292fe419aa299ef9477afedfe487a..822a4b38fef3192e3d03c4f2222b3762d7b9b86c 100644
--- a/Solver/linearSystemCSR.h
+++ b/Solver/linearSystemCSR.h
@@ -11,134 +11,132 @@
#include "GmshMessage.h"
#include "linearSystem.h"
-namespace gsolver {
- typedef int INDEX_TYPE ;
- typedef struct {
- int nmax;
- int size;
- int incr;
- int n;
- int isorder;
- char *array;
- } CSRList_T;
+typedef int INDEX_TYPE ;
+typedef struct {
+ int nmax;
+ int size;
+ int incr;
+ int n;
+ int isorder;
+ char *array;
+} CSRList_T;
- void CSRList_Add(CSRList_T *liste, void *data);
- int CSRList_Nbr(CSRList_T *liste);
-
- template <class scalar>
- class linearSystemCSR : public linearSystem<scalar> {
- protected:
- bool sorted;
- char *something;
- CSRList_T *a_,*ai_,*ptr_,*jptr_;
- std::vector<scalar> *_b, *_x;
- public:
- linearSystemCSR()
- : sorted(false), a_(0) {}
- virtual bool isAllocated() const { return a_ != 0; }
- virtual void allocate(int) ;
- virtual void clear()
- {
- allocate(0);
- }
- virtual ~linearSystemCSR()
- {
- allocate(0);
- }
- virtual void addToMatrix ( int il, int ic, double val)
- {
- // if (sorted)throw;
-
- INDEX_TYPE *jptr = (INDEX_TYPE*) jptr_->array;
- INDEX_TYPE *ptr = (INDEX_TYPE*) ptr_->array;
- INDEX_TYPE *ai = (INDEX_TYPE*) ai_->array;
- scalar *a = ( scalar * ) a_->array;
-
- INDEX_TYPE position_ = jptr[il];
-
- if(something[il]) {
- while(1){
- if(ai[position_] == ic){
- a[position_] += val;
- // if (il == 0) printf("FOUND %d %d %d\n",il,ic,position_);
- return;
- }
- if (ptr[position_] == 0)break;
- position_ = ptr[position_];
- }
- }
-
- INDEX_TYPE zero = 0;
- CSRList_Add (a_, &val);
- CSRList_Add (ai_, &ic);
- CSRList_Add (ptr_, &zero);
- // The pointers may have been modified
- // if there has been a reallocation in CSRList_Add
-
- ptr = (INDEX_TYPE*) ptr_->array;
- ai = (INDEX_TYPE*) ai_->array;
- a = (scalar*) a_->array;
-
- INDEX_TYPE n = CSRList_Nbr(a_) - 1;
-
- if(!something[il]) {
- jptr[il] = n;
- something[il] = 1;
+void CSRList_Add(CSRList_T *liste, void *data);
+int CSRList_Nbr(CSRList_T *liste);
+
+template <class scalar>
+class linearSystemCSR : public linearSystem<scalar> {
+ protected:
+ bool sorted;
+ char *something;
+ CSRList_T *a_,*ai_,*ptr_,*jptr_;
+ std::vector<scalar> *_b, *_x;
+ public:
+ linearSystemCSR()
+ : sorted(false), a_(0) {}
+ virtual bool isAllocated() const { return a_ != 0; }
+ virtual void allocate(int) ;
+ virtual void clear()
+ {
+ allocate(0);
+ }
+ virtual ~linearSystemCSR()
+ {
+ allocate(0);
+ }
+ virtual void addToMatrix ( int il, int ic, double val)
+ {
+ // if (sorted)throw;
+
+ INDEX_TYPE *jptr = (INDEX_TYPE*) jptr_->array;
+ INDEX_TYPE *ptr = (INDEX_TYPE*) ptr_->array;
+ INDEX_TYPE *ai = (INDEX_TYPE*) ai_->array;
+ scalar *a = ( scalar * ) a_->array;
+
+ INDEX_TYPE position_ = jptr[il];
+
+ if(something[il]) {
+ while(1){
+ if(ai[position_] == ic){
+ a[position_] += val;
+ // if (il == 0) printf("FOUND %d %d %d\n",il,ic,position_);
+ return;
+ }
+ if (ptr[position_] == 0)break;
+ position_ = ptr[position_];
}
- else ptr[position_] = n;
- }
+ }
- virtual scalar getFromMatrix (int _row, int _col) const
- {
- throw;
- }
- virtual void addToRightHandSide(int _row, scalar _val)
- {
- if(_val != 0.0) (*_b)[_row] += _val;
- }
- virtual scalar getFromRightHandSide(int _row) const
- {
- return (*_b)[_row];
- }
- virtual scalar getFromSolution(int _row) const
- {
- return (*_x)[_row];
- }
- virtual void zeroMatrix()
- {
- int N=CSRList_Nbr(a_);
- scalar *a = (scalar*) a_->array;
- for (int i=0;i<N;i++)a[i]=0;
- }
- virtual void zeroRightHandSide()
- {
- for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+ INDEX_TYPE zero = 0;
+ CSRList_Add (a_, &val);
+ CSRList_Add (ai_, &ic);
+ CSRList_Add (ptr_, &zero);
+ // The pointers may have been modified
+ // if there has been a reallocation in CSRList_Add
+
+ ptr = (INDEX_TYPE*) ptr_->array;
+ ai = (INDEX_TYPE*) ai_->array;
+ a = (scalar*) a_->array;
+
+ INDEX_TYPE n = CSRList_Nbr(a_) - 1;
+
+ if(!something[il]) {
+ jptr[il] = n;
+ something[il] = 1;
}
- };
+ else ptr[position_] = n;
+ }
+
+ virtual scalar getFromMatrix (int _row, int _col) const
+ {
+ throw;
+ }
+ virtual void addToRightHandSide(int _row, scalar _val)
+ {
+ if(_val != 0.0) (*_b)[_row] += _val;
+ }
+ virtual scalar getFromRightHandSide(int _row) const
+ {
+ return (*_b)[_row];
+ }
+ virtual scalar getFromSolution(int _row) const
+ {
+ return (*_x)[_row];
+ }
+ virtual void zeroMatrix()
+ {
+ int N=CSRList_Nbr(a_);
+ scalar *a = (scalar*) a_->array;
+ for (int i=0;i<N;i++)a[i]=0;
+ }
+ virtual void zeroRightHandSide()
+ {
+ for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+ }
+};
- template <class scalar>
- class linearSystemCSRGmm : public linearSystemCSR<scalar> {
- private:
- double _prec;
- int _noisy, _gmres;
- public:
- linearSystemCSRGmm()
- : _prec(1.e-8), _noisy(0), _gmres(0) {}
- virtual ~linearSystemCSRGmm()
- {}
- void setPrec(double p){ _prec = p; }
- void setNoisy(int n){ _noisy = n; }
- void setGmres(int n){ _gmres = n; }
- virtual int systemSolve()
+template <class scalar>
+class linearSystemCSRGmm : public linearSystemCSR<scalar> {
+ private:
+ double _prec;
+ int _noisy, _gmres;
+ public:
+ linearSystemCSRGmm()
+ : _prec(1.e-8), _noisy(0), _gmres(0) {}
+ virtual ~linearSystemCSRGmm()
+ {}
+ void setPrec(double p){ _prec = p; }
+ void setNoisy(int n){ _noisy = n; }
+ void setGmres(int n){ _gmres = n; }
+ virtual int systemSolve()
#if defined(HAVE_GMM)
- ;
+ ;
#else
- {
- Msg::Error("Gmm++ is not available in this version of Gmsh");
- return 0;
- }
+ {
+ Msg::Error("Gmm++ is not available in this version of Gmsh");
+ return 0;
+ }
#endif
- };
-}
+};
#endif
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index 85ab46429434241497842b5a14839346a4532114..831fc60110485d8ea3910c9f3ba0ebdc811594de 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -14,72 +14,69 @@
#include <stdlib.h>
#include <set>
-namespace gsolver {
- template <class scalar>
- class linearSystemFull : public linearSystem<scalar> {
- private:
- gmshMatrix<scalar> *_a;
- gmshVector<scalar> *_b, *_x;
- public :
- linearSystemFull() : _a(0), _b(0), _x(0){}
- virtual bool isAllocated() const { return _a != 0; }
- virtual void allocate(int _nbRows)
- {
- clear();
- _a = new gmshMatrix<scalar>(_nbRows, _nbRows);
- _b = new gmshVector<scalar>(_nbRows);
- _x = new gmshVector<scalar>(_nbRows);
+template <class scalar>
+class linearSystemFull : public linearSystem<scalar> {
+ private:
+ gmshMatrix<scalar> *_a;
+ gmshVector<scalar> *_b, *_x;
+ public :
+ linearSystemFull() : _a(0), _b(0), _x(0){}
+ virtual bool isAllocated() const { return _a != 0; }
+ virtual void allocate(int _nbRows)
+ {
+ clear();
+ _a = new gmshMatrix<scalar>(_nbRows, _nbRows);
+ _b = new gmshVector<scalar>(_nbRows);
+ _x = new gmshVector<scalar>(_nbRows);
+ }
+ virtual ~linearSystemFull()
+ {
+ clear();
+ }
+ virtual void clear()
+ {
+ if(_a){
+ delete _a;
+ delete _b;
+ delete _x;
}
-
- virtual ~linearSystemFull()
- {
- clear();
- }
-
- virtual void clear()
- {
- if(_a){
- delete _a;
- delete _b;
- delete _x;
- }
- _a = 0;
- }
- virtual void addToMatrix(int _row, int _col, scalar _val)
- {
- if(_val != 0.0) (*_a)(_row, _col) += _val;
- }
- virtual scalar getFromMatrix(int _row, int _col) const
- {
- return (*_a)(_row, _col);
- }
- virtual void addToRightHandSide(int _row, scalar _val)
- {
- if(_val != 0.0) (*_b)(_row) += _val;
- }
- virtual scalar getFromRightHandSide(int _row) const
- {
- return (*_b)(_row);
- }
- virtual scalar getFromSolution(int _row) const
- {
- return (*_x)(_row);
- }
- virtual void zeroMatrix()
- {
- _a->set_all(0.);
- }
- virtual void zeroRightHandSide()
- {
- for(int i = 0; i < _b->size(); i++) (*_b)(i) = 0.;
- }
- virtual int systemSolve()
- {
- if (_b->size())
- _a->lu_solve(*_b, *_x);
- // _x->print("********* mySol");
- return 1;
- }
- };
-}
+ _a = 0;
+ }
+ virtual void addToMatrix(int _row, int _col, scalar _val)
+ {
+ if(_val != 0.0) (*_a)(_row, _col) += _val;
+ }
+ virtual scalar getFromMatrix(int _row, int _col) const
+ {
+ return (*_a)(_row, _col);
+ }
+ virtual void addToRightHandSide(int _row, scalar _val)
+ {
+ if(_val != 0.0) (*_b)(_row) += _val;
+ }
+ virtual scalar getFromRightHandSide(int _row) const
+ {
+ return (*_b)(_row);
+ }
+ virtual scalar getFromSolution(int _row) const
+ {
+ return (*_x)(_row);
+ }
+ virtual void zeroMatrix()
+ {
+ _a->set_all(0.);
+ }
+ virtual void zeroRightHandSide()
+ {
+ for(int i = 0; i < _b->size(); i++) (*_b)(i) = 0.;
+ }
+ virtual int systemSolve()
+ {
+ if (_b->size())
+ _a->lu_solve(*_b, *_x);
+ // _x->print("********* mySol");
+ return 1;
+ }
+};
+
#endif
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
index 822d28fc67079d164089b8a8a9d71a465c818076..3f0c81bf1a51f81b32f193e64befc3e8696fbe04 100644
--- a/Solver/linearSystemGMM.h
+++ b/Solver/linearSystemGMM.h
@@ -14,109 +14,106 @@
#if defined(HAVE_GMM)
#include <gmm.h>
-namespace gsolver {
- template <class scalar>
- class linearSystemGmm : public linearSystem<scalar> {
- private:
- gmm::row_matrix<gmm::wsvector<scalar> > *_a;
- std::vector<scalar> *_b, *_x;
- double _prec;
- int _noisy, _gmres;
- public:
- linearSystemGmm()
- : _a(0), _b(0), _x(0), _prec(1.e-8), _noisy(0), _gmres(0) {}
- virtual bool isAllocated() const { return _a != 0; }
- virtual void allocate(int _nbRows)
- {
- clear();
- _a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
- _b = new std::vector<scalar>(_nbRows);
- _x = new std::vector<scalar>(_nbRows);
- }
- virtual ~linearSystemGmm()
- {
- clear();
- }
-
- virtual void clear()
- {
- if (_a){
- delete _a;
- delete _b;
- delete _x;
- }
- _a = 0;
- }
- virtual void addToMatrix(int _row, int _col, scalar _val)
- {
- if(_val != 0.0) (*_a)(_row, _col) += _val;
- }
- virtual scalar getFromMatrix (int _row, int _col) const
- {
- return (*_a)(_row, _col);
- }
- virtual void addToRightHandSide(int _row, scalar _val)
- {
- if(_val != 0.0) (*_b)[_row] += _val;
- }
- virtual scalar getFromRightHandSide(int _row) const
- {
- return (*_b)[_row];
- }
- virtual scalar getFromSolution(int _row) const
- {
- return (*_x)[_row];
- }
- virtual void zeroMatrix()
- {
- gmm::clear(*_a);
- }
- virtual void zeroRightHandSide()
- {
- for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
- }
- void setPrec(double p){ _prec = p; }
- void setNoisy(int n){ _noisy = n; }
- void setGmres(int n){ _gmres = n; }
- virtual int systemSolve()
- {
- //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25, 0.);
- gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30, 1.e-10);
- gmm::iteration iter(_prec);
- iter.set_noisy(_noisy);
- if(_gmres) gmm::gmres(*_a, *_x, *_b, P, 100, iter);
- else gmm::cg(*_a, *_x, *_b, P, iter);
- return 1;
+
+template <class scalar>
+class linearSystemGmm : public linearSystem<scalar> {
+ private:
+ gmm::row_matrix<gmm::wsvector<scalar> > *_a;
+ std::vector<scalar> *_b, *_x;
+ double _prec;
+ int _noisy, _gmres;
+ public:
+ linearSystemGmm()
+ : _a(0), _b(0), _x(0), _prec(1.e-8), _noisy(0), _gmres(0) {}
+ virtual bool isAllocated() const { return _a != 0; }
+ virtual void allocate(int _nbRows)
+ {
+ clear();
+ _a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
+ _b = new std::vector<scalar>(_nbRows);
+ _x = new std::vector<scalar>(_nbRows);
+ }
+ virtual ~linearSystemGmm()
+ {
+ clear();
+ }
+ virtual void clear()
+ {
+ if (_a){
+ delete _a;
+ delete _b;
+ delete _x;
}
- };
-}
+ _a = 0;
+ }
+ virtual void addToMatrix(int _row, int _col, scalar _val)
+ {
+ if(_val != 0.0) (*_a)(_row, _col) += _val;
+ }
+ virtual scalar getFromMatrix (int _row, int _col) const
+ {
+ return (*_a)(_row, _col);
+ }
+ virtual void addToRightHandSide(int _row, scalar _val)
+ {
+ if(_val != 0.0) (*_b)[_row] += _val;
+ }
+ virtual scalar getFromRightHandSide(int _row) const
+ {
+ return (*_b)[_row];
+ }
+ virtual scalar getFromSolution(int _row) const
+ {
+ return (*_x)[_row];
+ }
+ virtual void zeroMatrix()
+ {
+ gmm::clear(*_a);
+ }
+ virtual void zeroRightHandSide()
+ {
+ for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+ }
+ void setPrec(double p){ _prec = p; }
+ void setNoisy(int n){ _noisy = n; }
+ void setGmres(int n){ _gmres = n; }
+ virtual int systemSolve()
+ {
+ //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25, 0.);
+ gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30, 1.e-10);
+ gmm::iteration iter(_prec);
+ iter.set_noisy(_noisy);
+ if(_gmres) gmm::gmres(*_a, *_x, *_b, P, 100, iter);
+ else gmm::cg(*_a, *_x, *_b, P, iter);
+ return 1;
+ }
+};
#else
-namespace gsolver {
- template <class scalar>
- class linearSystemGmm : public linearSystem<scalar> {
- public :
- linearSystemGmm()
- {
- Msg::Error("Gmm++ is not available in this version of Gmsh");
- }
- virtual bool isAllocated() const { return false; }
- virtual void allocate(int nbRows) {}
- virtual void addToMatrix(int _row, int _col, scalar val) {}
- virtual scalar getFromMatrix(int _row, int _col) const { return 0.; }
- virtual void addToRightHandSide(int _row, scalar val) {}
- virtual scalar getFromRightHandSide(int _row) const { return 0.; }
- virtual scalar getFromSolution(int _row) const { return 0.; }
- virtual void zeroMatrix() {}
- virtual void zeroRightHandSide() {}
- virtual int systemSolve() { return 0; }
- void setPrec(double p){}
- virtual void clear(){}
- void setNoisy(int n){}
- void setGmres(int n){}
- };
-}
+template <class scalar>
+class linearSystemGmm : public linearSystem<scalar> {
+ public :
+ linearSystemGmm()
+ {
+ Msg::Error("Gmm++ is not available in this version of Gmsh");
+ }
+ virtual bool isAllocated() const { return false; }
+ virtual void allocate(int nbRows) {}
+ virtual void addToMatrix(int _row, int _col, scalar val) {}
+ virtual scalar getFromMatrix(int _row, int _col) const { return 0.; }
+ virtual void addToRightHandSide(int _row, scalar val) {}
+ virtual scalar getFromRightHandSide(int _row) const { return 0.; }
+ virtual scalar getFromSolution(int _row) const { return 0.; }
+ virtual void zeroMatrix() {}
+ virtual void zeroRightHandSide() {}
+ virtual int systemSolve() { return 0; }
+ void setPrec(double p){}
+ virtual void clear(){}
+ void setNoisy(int n){}
+ void setGmres(int n){}
+};
+
#endif
#endif
diff --git a/Solver/linearSystemTAUCS.cpp b/Solver/linearSystemTAUCS.cpp
index b767a37dca9b003af66b219b218dc18cb2fa2250..f2806251e2a9d7492a183e32aedb261c2b9eb1be 100644
--- a/Solver/linearSystemTAUCS.cpp
+++ b/Solver/linearSystemTAUCS.cpp
@@ -3,7 +3,6 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
@@ -14,50 +13,49 @@ extern "C" {
#include "taucs.h"
}
-namespace gsolver {
- template <class scalar>
- void sortColumns(int NbLines,
- int nnz,
- INDEX_TYPE *ptr,
- INDEX_TYPE *jptr,
- INDEX_TYPE *ai,
- scalar *a);
-
- template<>
- int linearSystemCSRTaucs<double>::systemSolve()
- {
- if(!sorted){
- sortColumns(_b->size(),
- CSRList_Nbr(a_),
- (INDEX_TYPE *) ptr_->array,
- (INDEX_TYPE *) jptr_->array,
- (INDEX_TYPE *) ai_->array,
- (double*) a_->array);
- }
- sorted = true;
-
- taucs_ccs_matrix myVeryCuteTaucsMatrix;
- myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
- //myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ptr_->array;
- //myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)ai_->array;
- myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ai_->array;
- myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)jptr_->array;
- myVeryCuteTaucsMatrix.values.d = (double*) a_->array;
- myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
-
- char* options[] = { "taucs.factor.LLT=true", NULL };
- Msg::Info("TAUCS solves %d unknowns", _b->size());
- int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
- NULL,
- 1,
- &(*_x)[0],
- &(*_b)[0],
- options,
- NULL);
- if (result != TAUCS_SUCCESS){
- Msg::Error("Taucs Was Not Successfull %d",result);
- }
- return 1;
+template <class scalar>
+void sortColumns(int NbLines,
+ int nnz,
+ INDEX_TYPE *ptr,
+ INDEX_TYPE *jptr,
+ INDEX_TYPE *ai,
+ scalar *a);
+
+template<>
+int linearSystemCSRTaucs<double>::systemSolve()
+{
+ if(!sorted){
+ sortColumns(_b->size(),
+ CSRList_Nbr(a_),
+ (INDEX_TYPE *) ptr_->array,
+ (INDEX_TYPE *) jptr_->array,
+ (INDEX_TYPE *) ai_->array,
+ (double*) a_->array);
}
+ sorted = true;
+
+ taucs_ccs_matrix myVeryCuteTaucsMatrix;
+ myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
+ //myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ptr_->array;
+ //myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)ai_->array;
+ myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ai_->array;
+ myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)jptr_->array;
+ myVeryCuteTaucsMatrix.values.d = (double*) a_->array;
+ myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
+
+ char* options[] = { "taucs.factor.LLT=true", NULL };
+ Msg::Info("TAUCS solves %d unknowns", _b->size());
+ int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
+ NULL,
+ 1,
+ &(*_x)[0],
+ &(*_b)[0],
+ options,
+ NULL);
+ if (result != TAUCS_SUCCESS){
+ Msg::Error("Taucs Was Not Successfull %d",result);
+ }
+ return 1;
}
+
#endif
diff --git a/Solver/linearSystemTAUCS.h b/Solver/linearSystemTAUCS.h
index 73a43686bd1e5e382f797a664ce63fced449cce8..fd7eef60e1ae00f769b30b5f3fbeee908404cc06 100644
--- a/Solver/linearSystemTAUCS.h
+++ b/Solver/linearSystemTAUCS.h
@@ -8,27 +8,22 @@
#include "linearSystemCSR.h"
-namespace gsolver {
- template <class scalar>
- class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
- private:
- public:
- linearSystemCSRTaucs()
- {}
- virtual ~linearSystemCSRTaucs()
- {}
- virtual void addToMatrix ( int il, int ic, double val) {
- if (il <= ic)
- linearSystemCSR<scalar>::addToMatrix(il,ic,val);
- }
-
- virtual int systemSolve()
+template <class scalar>
+class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
+ public:
+ linearSystemCSRTaucs(){}
+ virtual ~linearSystemCSRTaucs(){}
+ virtual void addToMatrix(int il, int ic, double val)
+ {
+ if (il <= ic)
+ linearSystemCSR<scalar>::addToMatrix(il,ic,val);
+ }
+ virtual int systemSolve()
#if defined(HAVE_TAUCS)
-;
+ ;
#else
- {throw;}
+ { throw; }
#endif
- };
-}
+};
#endif
diff --git a/Solver/termOfFormulation.h b/Solver/termOfFormulation.h
deleted file mode 100644
index c1398e7675a402824e2e01d200604aa99d6ed50e..0000000000000000000000000000000000000000
--- a/Solver/termOfFormulation.h
+++ /dev/null
@@ -1,137 +0,0 @@
-// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _TERM_OF_FORMULATION_H_
-#define _TERM_OF_FORMULATION_H_
-
-#include <math.h>
-#include <map>
-#include <vector>
-#include "GmshMatrix.h"
-#include "gmshFunction.h"
-#include "dofManager.h"
-#include "GModel.h"
-#include "MElement.h"
-
-namespace gsolver {
- // a nodal finite element term : variables are always defined at nodes
- // of the mesh
- template<class dataVec, class dataMat>
- class femTerm {
- protected:
- GModel *_gm;
- public:
- // return the number of columns of the element matrix
- virtual int sizeOfC(MElement*) const = 0;
- // return the number of rows of the element matrix
- virtual int sizeOfR(MElement*) const = 0;
- // in a given element, return the dof associated to a given row (column)
- // of the local element matrix
- virtual Dof getLocalDofR(MElement *e, int iRow) const = 0;
- // default behavior : symmetric
- virtual Dof getLocalDofC(MElement *e, int iCol) const
- { return getLocalDofR(e, iCol); }
- public:
- femTerm(GModel *gm) : _gm(gm) {}
- virtual ~femTerm (){}
- virtual void elementMatrix(MElement *e, gmshMatrix<dataMat> &m) const = 0;
- virtual void elementVector(MElement *e, gmshVector<dataVec> &m) const {
- m.scale(0.0);
- }
- void addToMatrix(dofManager<dataVec,dataMat> &dm, GEntity *ge) const
- {
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- addToMatrix(dm, e);
- }
- }
- void addToMatrix(dofManager<dataVec,dataMat> &dm, MElement *e) const
- {
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- gmshMatrix<dataMat> localMatrix(nbR, nbC);
- elementMatrix(e, localMatrix);
- addToMatrix(dm, localMatrix, e);
- }
- void addToMatrix(dofManager<dataVec,dataMat> &dm,
- gmshMatrix<dataMat> &localMatrix,
- MElement *e) const
- {
- const int nbR = localMatrix.size1();
- const int nbC = localMatrix.size2();
- for (int j = 0; j < nbR; j++){
- Dof R = getLocalDofR(e, j);
- for (int k = 0; k < nbC; k++){
- Dof C = getLocalDofC(e, k);
- dm.assemble(R,C, localMatrix(j, k));
- }
- }
- }
- void dirichletNodalBC(int physical,
- int dim,
- int comp,
- int field,
- const gmshFunction<dataVec> &e,
- dofManager<dataVec,dataMat> &dm)
- {
- std::vector<MVertex *> v;
- GModel *m = _gm;
- m->getMeshVerticesForPhysicalGroup(dim, physical, v);
- for (unsigned int i = 0; i < v.size(); i++)
- dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
- }
- void neumannNodalBC(int physical,
- int dim,
- int comp,
- int field,
- const gmshFunction<dataVec> &fct,
- dofManager<dataVec,dataMat> &dm)
- {
- std::map<int, std::vector<GEntity*> > groups[4];
- GModel *m = _gm;
- m->getPhysicalGroups(groups);
- std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
- if (it == groups[dim].end()) return;
- double jac[3][3];
- double sf[256];
- for (unsigned int i = 0; i < it->second.size(); ++i){
- GEntity *ge = it->second[i];
- for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
- MElement *e = ge->getMeshElement(j);
- int integrationOrder = 2 * e->getPolynomialOrder();
- int nbNodes = e->getNumVertices();
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- for (int ip = 0; ip < npts; ip++){
- const double u = GP[ip].pt[0];
- const double v = GP[ip].pt[1];
- const double w = GP[ip].pt[2];
- const double weight = GP[ip].weight;
- const double detJ = e->getJacobian(u, v, w, jac);
- SPoint3 p; e->pnt(u, v, w, p);
- e->getShapeFunctions(u, v, w, sf);
- const dataVec FCT = fct(p.x(), p.y(), p.z());
- for (int k = 0; k < nbNodes; k++){
- dm.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
- }
- }
- }
- }
- }
- void addToRightHandSide (dofManager<dataVec,dataMat> &dm, GEntity *ge) const {
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement (i);
- int nbR = sizeOfR(e);
- gmshVector<dataVec> V (nbR);
- elementVector (e, V);
- // assembly
- for (int j=0;j<nbR;j++)dm.assemble(getLocalDofR(e,j),V(j));
- }
- }
- }; // end of class definition
-}// end of namespace
-
-#endif