From 609106dced9e83f99f315114d66ef63ee19d7cf3 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 20 Sep 2009 18:07:57 +0000
Subject: [PATCH] * ported all old solver code to new api * removed old code
---
CMakeLists.txt | 7 +-
Geo/GFaceCompound.cpp | 120 ++++++-----
Geo/GFaceCompound.h | 6 +-
Numeric/CMakeLists.txt | 5 -
Numeric/gmshAssembler.h | 194 ------------------
Numeric/gmshConvexCombination.cpp | 27 ---
Numeric/gmshConvexCombination.h | 36 ----
Numeric/gmshCrossConf.cpp | 58 ------
Numeric/gmshCrossConf.h | 56 ------
Numeric/gmshDistance.cpp | 85 --------
Numeric/gmshDistance.h | 51 -----
Numeric/gmshLaplace.cpp | 59 ------
Numeric/gmshLaplace.h | 49 -----
Numeric/gmshLinearSystem.h | 31 ---
Numeric/gmshLinearSystemCSR.cpp | 324 ------------------------------
Numeric/gmshLinearSystemCSR.h | 151 --------------
Numeric/gmshLinearSystemFull.h | 121 -----------
Numeric/gmshLinearSystemGmm.h | 123 ------------
Numeric/gmshTermOfFormulation.h | 188 -----------------
Solver/convexCombinationTerm.h | 2 +-
Solver/crossConfTerm.h | 6 +-
Solver/helmholtzTerm.h | 14 +-
Solver/laplaceTerm.h | 19 +-
23 files changed, 79 insertions(+), 1653 deletions(-)
delete mode 100644 Numeric/gmshAssembler.h
delete mode 100644 Numeric/gmshConvexCombination.cpp
delete mode 100644 Numeric/gmshConvexCombination.h
delete mode 100644 Numeric/gmshCrossConf.cpp
delete mode 100644 Numeric/gmshCrossConf.h
delete mode 100644 Numeric/gmshDistance.cpp
delete mode 100644 Numeric/gmshDistance.h
delete mode 100644 Numeric/gmshLaplace.cpp
delete mode 100644 Numeric/gmshLaplace.h
delete mode 100644 Numeric/gmshLinearSystem.h
delete mode 100644 Numeric/gmshLinearSystemCSR.cpp
delete mode 100644 Numeric/gmshLinearSystemCSR.h
delete mode 100644 Numeric/gmshLinearSystemFull.h
delete mode 100644 Numeric/gmshLinearSystemGmm.h
delete mode 100644 Numeric/gmshTermOfFormulation.h
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f445c50b5d..25842832fd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -68,13 +68,10 @@ set(GMSH_API
Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h
Mesh/meshGFaceDelaunayInsertion.h
+ Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
+ Solver/linearSystem.h Solver/linearSystemGmm.h Solver/linearSystemFull.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
Graphics/drawContext.h
- Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
- Numeric/gmshAssembler.h Numeric/gmshTermOfFormulation.h
- Numeric/gmshLaplace.h Numeric/gmshLinearSystem.h
- Numeric/gmshLinearSystemGmm.h Numeric/gmshLinearSystemFull.h
- Numeric/gmshLinearSystemCSR.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
contrib/kbipack/gmp_blas.h
contrib/DiscreteIntegration/DILevelset.h)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 414ffb0fea..f1444dc57d 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -10,23 +10,23 @@
#include "MTriangle.h"
#include "Numeric.h"
#include "Octree.h"
-#include "gmshAssembler.h"
-#include "gmshLaplace.h"
-#include "gmshDistance.h"
#include "SBoundingBox3d.h"
#include "SPoint3.h"
-#include "gmshCrossConf.h"
-#include "gmshConvexCombination.h"
-#include "gmshLinearSystemGmm.h"
-#include "gmshLinearSystemCSR.h"
-#include "gmshLinearSystemFull.h"
#include "functionSpace.h"
#include "robustPredicates.h"
#include "meshGFaceOptimize.h"
#include "MElementCut.h"
#include "GEntity.h"
-
-static void fixEdgeToValue(GEdge *ed, double value, gmshAssembler<double> &myAssembler)
+#include "dofManager.h"
+#include "laplaceTerm.h"
+#include "distanceTerm.h"
+#include "crossConfTerm.h"
+#include "convexCombinationTerm.h"
+#include "linearSystemGmm.h"
+#include "linearSystemCSR.h"
+#include "linearSystemFull.h"
+
+static void fixEdgeToValue(GEdge *ed, double value, dofManager<double, double> &myAssembler)
{
myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], 0, 1, value);
@@ -64,7 +64,7 @@ static void computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
double x1, x2, y1, y2;
x1 = u(i, 0); y1 = u(i, 1);
x2 = u(next, 0); y2 = u(next, 1);
- for( int j = i + 2; j < nbPts; j++){
+ for(int j = i + 2; j < nbPts; j++){
int jnext = j + 1;
if(j == nbPts - 1) jnext = 0;
double x3, x4, y3, y4, x,y;
@@ -161,11 +161,11 @@ static void myPolygon(std::vector<MElement*> &vTri, std::vector<MVertex*> &vPoly
while(ite != ePoly.end()){
MVertex *vB = ite->getVertex(0);
MVertex *vE = ite->getVertex(1);
- if( vB == vPoly.back()){
+ if(vB == vPoly.back()){
if(vE != vINIT) vPoly.push_back(vE);
ePoly.erase(ite);
}
- else if( vE == vPoly.back()){
+ else if(vE == vPoly.back()){
if(vB != vINIT) vPoly.push_back(vB);
ePoly.erase(ite);
}
@@ -356,7 +356,7 @@ void GFaceCompound::getBoundingEdges()
//in case the bounding edges are explicitely given
if(_U0.size()){
- std::list<GEdge*> :: const_iterator it = _U0.begin();
+ std::list<GEdge*>::const_iterator it = _U0.begin();
for( ; it != _U0.end() ; ++it){
l_edges.push_back(*it);
(*it)->addFace(this);
@@ -377,10 +377,10 @@ void GFaceCompound::getBoundingEdges()
(*itf)->addFace(this);
}
- computeALoop(_unique,_U0);
+ computeALoop(_unique, _U0);
while(!_unique.empty()){
- computeALoop(_unique,_U1);
- printf("%d in unique\n",_unique.size());
+ computeALoop(_unique, _U1);
+ printf("%d in unique\n", (int)_unique.size());
_interior_loops.push_back(_U1);
}
}
@@ -393,7 +393,7 @@ void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
std::list<GFace*>::iterator it = _compound.begin();
for( ; it != _compound.end(); ++it){
std::list<GEdge*> ed = (*it)->edges();
- std::list<GEdge*> :: iterator ite = ed.begin();
+ std::list<GEdge*>::iterator ite = ed.begin();
for( ; ite != ed.end(); ++ite){
_touched.insert(*ite);
}
@@ -401,7 +401,7 @@ void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
it = _compound.begin();
for( ; it != _compound.end(); ++it){
std::list<GEdge*> ed = (*it)->edges();
- std::list<GEdge*> :: iterator ite = ed.begin();
+ std::list<GEdge*>::iterator ite = ed.begin();
for( ; ite != ed.end() ; ++ite){
if(!(*ite)->degenerate(0) && _touched.count(*ite) == 1) {
_unique.insert(*ite); }
@@ -474,21 +474,21 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
std::list<GEdge*> &U1, std::list<GEdge*> &V1,
- gmshLinearSystem<double> *lsys)
+ linearSystem<double> *lsys)
: GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0),
_lsys(lsys)
{
if (!_lsys) {
#if defined(HAVE_GMM)
- gmshLinearSystemGmm<double> *_lsysb = new gmshLinearSystemGmm<double>;
- //gmshLinearSystemCSRGmm<double> lsys;
+ linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
+ //linearSystemCSRGmm<double> lsys;
_lsysb->setPrec(1.e-15);
_lsysb->setGmres(1);
if(Msg::GetVerbosity() == 99) _lsysb->setNoisy(2);
_lsys = _lsysb;
#else
- _lsys = new gmshLinearSystemFull<double>;
+ _lsys = new linearSystemFull<double>;
#endif
}
@@ -699,7 +699,7 @@ void GFaceCompound::parametrize(iterationStep step) const
{
Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
- gmshAssembler<double> myAssembler(_lsys);
+ dofManager<double, double> myAssembler(_lsys);
simpleFunction<double> ONE(1.0);
if(_type == UNITCIRCLE){
@@ -721,7 +721,7 @@ void GFaceCompound::parametrize(iterationStep step) const
}
else if(_type == SQUARE){
if(step == ITERU){
- std::list<GEdge*> :: const_iterator it = _U0.begin();
+ std::list<GEdge*>::const_iterator it = _U0.begin();
for( ; it != _U0.end(); ++it){
fixEdgeToValue(*it, 0.0, myAssembler);
}
@@ -731,7 +731,7 @@ void GFaceCompound::parametrize(iterationStep step) const
}
}
else if(step == ITERV){
- std::list<GEdge*> :: const_iterator it = _V0.begin();
+ std::list<GEdge*>::const_iterator it = _V0.begin();
for( ; it != _V0.end(); ++it){
fixEdgeToValue(*it, 0.0, myAssembler);
}
@@ -756,13 +756,13 @@ void GFaceCompound::parametrize(iterationStep step) const
Msg::Debug("Creating term %d dofs numbered %d fixed",
myAssembler.sizeOfR(), myAssembler.sizeOfF());
- //gmshConvexCombinationTerm laplace(model(), &ONE, 1);
- gmshLaplaceTerm laplace(model(), &ONE, 1);
+ //convexCombinationTerm laplace(model(), 1, &ONE);
+ laplaceTerm<double> laplace(model(), 1, &ONE);
it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- laplace.addToMatrix(myAssembler, t);
+ SElement se((*it)->triangles[i]);
+ laplace.addToMatrix(myAssembler, &se);
}
}
@@ -772,7 +772,7 @@ void GFaceCompound::parametrize(iterationStep step) const
it = _compound.begin();
- std::set<MVertex *>allNodes;
+ std::set<MVertex *> allNodes;
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
@@ -784,10 +784,10 @@ void GFaceCompound::parametrize(iterationStep step) const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value = myAssembler.getDofValue(v,0,1);
+ double value = myAssembler.getDofValue(v, 0, 1);
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if(itf == coordinates.end()){
- SPoint3 p(0,0,0);
+ SPoint3 p(0, 0, 0);
p[step] = value;
coordinates[v] = p;
}
@@ -804,7 +804,7 @@ void GFaceCompound::parametrize_conformal() const
Msg::Info("Parametrizing Surface %d", tag());
- gmshAssembler<double> myAssembler(_lsys);
+ dofManager<double, double> myAssembler(_lsys);
std::vector<MVertex*> ordered;
std::vector<double> coords;
@@ -851,18 +851,18 @@ void GFaceCompound::parametrize_conformal() const
simpleFunction<double> ONE(1.0);
simpleFunction<double> MONE(-1.0 );
- gmshLaplaceTerm laplace1(model(), &ONE, 1);
- gmshLaplaceTerm laplace2(model(), &ONE, 2);
- gmshCrossConfTerm cross12(model(), &ONE, 1 , 2);
- gmshCrossConfTerm cross21(model(), &MONE, 2 , 1);
+ laplaceTerm<double> laplace1(model(), 1, &ONE);
+ laplaceTerm<double> laplace2(model(), 2, &ONE);
+ crossConfTerm cross12(model(), 1, 2, &ONE);
+ crossConfTerm cross21(model(), 2, 1, &MONE);
it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- laplace1.addToMatrix(myAssembler, t);
- laplace2.addToMatrix(myAssembler, t);
- cross12.addToMatrix(myAssembler, t);
- cross21.addToMatrix(myAssembler, t);
+ SElement se((*it)->triangles[i]);
+ laplace1.addToMatrix(myAssembler, &se);
+ laplace2.addToMatrix(myAssembler, &se);
+ cross12.addToMatrix(myAssembler, &se);
+ cross21.addToMatrix(myAssembler, &se);
}
}
@@ -897,14 +897,13 @@ void GFaceCompound::parametrize_conformal() const
void GFaceCompound::compute_distance() const
{
SBoundingBox3d bbox = GModel::current()->bounds();
- double L = norm(SVector3(bbox.max(), bbox.min()));
- //printf("L=%g \n", L);
- double mu = L/28;
- simpleFunction<double> DIFF(mu*mu);
- gmshAssembler<double> myAssembler(_lsys);
- gmshDistanceTerm distance(model(), &DIFF, 1);
-
- std::vector<MVertex*> ordered;
+ double L = norm(SVector3(bbox.max(), bbox.min()));
+ double mu = L / 28;
+ simpleFunction<double> DIFF(mu * mu), ONE(1.0);
+ dofManager<double, double> myAssembler(_lsys);
+ distanceTerm distance(model(), 1, 1, &DIFF, &ONE);
+
+ std::vector<MVertex*> ordered;
boundVertices(_U0, ordered);
for(unsigned int i = 0; i < ordered.size(); i++)
myAssembler.fixVertex(ordered[i], 0, 1, 0.0);
@@ -922,10 +921,10 @@ void GFaceCompound::compute_distance() const
it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- distance.addToMatrix(myAssembler, t);
+ SElement se((*it)->triangles[i]);
+ distance.addToMatrix(myAssembler, &se);
}
- distance.addToRightHandSide(myAssembler,*it);
+ distance.addToRightHandSide(myAssembler, *it);
}
Msg::Info("Distance Computation: Assembly done");
@@ -945,9 +944,9 @@ void GFaceCompound::compute_distance() const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value = std::min(0.9999, myAssembler.getDofValue(v,0,1));
- double dist = -mu*log(1.- value);
- coordinates[v] = SPoint3(dist,0.0,0.0);
+ double value = std::min(0.9999, myAssembler.getDofValue(v, 0, 1));
+ double dist = -mu * log(1. - value);
+ coordinates[v] = SPoint3(dist, 0.0, 0.0);
}
_lsys->clear();
@@ -957,7 +956,6 @@ void GFaceCompound::compute_distance() const
printf("--> Write distance function in file: XYZU-*.pos\n");
printf("--> Exit\n");
exit(1);
-
}
void GFaceCompound::computeNormals(std::map<MVertex*,SVector3> &normals) const
@@ -1220,7 +1218,7 @@ void GFaceCompound::buildOct() const
SBoundingBox3d bb;
int count = 0;
- std::list<GFace*> :: const_iterator it = _compound.begin();
+ std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -1280,7 +1278,7 @@ void GFaceCompound::buildOct() const
int GFaceCompound::genusGeom()
{
- std::list<GFace*> :: const_iterator it = _compound.begin();
+ std::list<GFace*>::const_iterator it = _compound.begin();
std::set<MEdge, Less_Edge> es;
std::set<MVertex*> vs;
int N = 0;
@@ -1299,7 +1297,7 @@ int GFaceCompound::genusGeom()
void GFaceCompound::printStuff() const
{
- std::list<GFace*> :: const_iterator it = _compound.begin();
+ std::list<GFace*>::const_iterator it = _compound.begin();
char name1[256], name2[256], name3[256];
char name4[256], name5[256], name6[256];
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index d22969ffe5..0a95efa9e3 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -12,7 +12,7 @@
#include "GEdge.h"
#include "GEdgeCompound.h"
#include "meshGFaceOptimize.h"
-#include "gmshLinearSystem.h"
+#include "linearSystem.h"
/*
A GFaceCompound is a model face that is the compound of model faces.
@@ -75,13 +75,13 @@ class GFaceCompound : public GFace {
virtual double curvature(MTriangle *t, double u, double v) const;
void printStuff() const;
bool trivial() const ;
- gmshLinearSystem <double> *_lsys;
+ linearSystem <double> *_lsys;
public:
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &U1,
std::list<GEdge*> &V0, std::list<GEdge*> &V1,
- gmshLinearSystem<double>* lsys =0);
+ linearSystem<double>* lsys =0);
virtual ~GFaceCompound();
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 85ac70f7dd..5ac9a6be5e 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -17,11 +17,6 @@ set(SRC
robustPredicates.cpp
mathEvaluator.cpp
Iso.cpp
- gmshLaplace.cpp
- gmshCrossConf.cpp
- gmshDistance.cpp
- gmshConvexCombination.cpp
- gmshLinearSystemCSR.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Numeric *.h)
diff --git a/Numeric/gmshAssembler.h b/Numeric/gmshAssembler.h
deleted file mode 100644
index 5bd43b57d8..0000000000
--- a/Numeric/gmshAssembler.h
+++ /dev/null
@@ -1,194 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_ASSEMBLER_H_
-#define _GMSH_ASSEMBLER_H_
-
-#include <map>
-#include <vector>
-#include "gmshLinearSystem.h"
-#include <stdlib.h>
-
-class MVertex;
-class MElement;
-
-struct gmshDofKey{
- MVertex *v;
- int comp;
- int field;
- gmshDofKey (MVertex *V, int iComp, int iField)
- : v(V), comp(iComp), field(iField) {}
- bool operator < (const gmshDofKey& d) const
- {
- if (v < d.v) return true;
- if (v > d.v) return false;
- if (comp < d.comp) return true;
- if (comp > d.comp) return false;
- if (field < d.field) return true;
- return false;
- }
-};
-
-template<class scalar>
-class gmshAssembler {
- private:
- std::map<gmshDofKey, int> numbering;
- std::map<gmshDofKey, scalar> fixed;
- std::map<gmshDofKey, scalar> fixedFULL;
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > > constraints;
- gmshLinearSystem<scalar> *lsys;
- public:
- gmshAssembler(gmshLinearSystem<scalar> *l) : lsys(l) {}
- inline void constraintVertex(MVertex*v, int iComp, int iField,
- std::vector<MVertex*> &verts,
- std::vector<scalar> &coeffs)
- {
- std::vector<std::pair<gmshDofKey, scalar> > constraint;
- gmshDofKey key(v, iComp, iField);
- for (unsigned int i = 0; i < verts.size(); i++){
- gmshDofKey key2(verts[i], iComp, iField);
- constraint.push_back(std::make_pair(key2, coeffs[i]));
- }
- constraints[key] = constraint;
- }
- inline void numberVertex(MVertex*v, int iComp, int iField)
- {
- gmshDofKey key(v, iComp, iField);
- if (fixed.find(key) != fixed.end()) return;
- if (constraints.find(key) != constraints.end()) return;
- std::map<gmshDofKey, int>::iterator it = numbering.find(key);
- if (it == numbering.end()){
- int size = numbering.size();
- numbering[key] = size;
- }
- }
- inline void numberVertex(MVertex*v, int iComp, int iField, int iField2)
- {
- gmshDofKey key(v, iComp, iField);
- gmshDofKey key2(v, iComp, iField2);
- if (fixed.find(key) != fixed.end()) return;
- if (fixed.find(key2) != fixed.end()) return;
- if (constraints.find(key) != constraints.end()) return;
- std::map<gmshDofKey, int> :: iterator it = numbering.find(key);
- if (it == numbering.end()){
- int size = numbering.size();
- numbering[key] = size;
- }
- }
- inline void fixVertex(MVertex*v, int iComp, int iField, scalar val)
- {
- fixed[gmshDofKey(v, iComp, iField)] = val;
- }
- inline void fixVertexFULL(MVertex*v, int iComp, int iField, scalar val)
- {
- fixedFULL[gmshDofKey(v, iComp, iField)] = val;
- }
- inline scalar getDofValue(MVertex *v, int iComp, int iField) const
- {
- gmshDofKey key(v, iComp, iField);
- {
- typename std::map<gmshDofKey, scalar>::const_iterator it = fixed.find(key);
- if (it != fixed.end()) return it->second;
- }
- {
- std::map<gmshDofKey, int>::const_iterator it = numbering.find(key);
- if (it != numbering.end())
- return lsys->getFromSolution(it->second);
- }
- {
- typename std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > >::
- const_iterator itConstr = constraints.find(key);
- if (itConstr != constraints.end()){
- scalar val = 0.;
- for (unsigned int i = 0; i < itConstr->second.size(); i++){
- const gmshDofKey &dofKeyConstr = itConstr->second[i].first;
- scalar valConstr = itConstr->second[i].second;
- val += getDofValue(dofKeyConstr.v, dofKeyConstr.comp, dofKeyConstr.field)
- * valConstr;
- }
- return val;
- }
- }
- return 0.0;
- }
- void assemble(MVertex *vR, int iCompR, int iFieldR,
- MVertex *vC, int iCompC, int iFieldC, scalar val)
- {
- if (!lsys->isAllocated()) lsys->allocate(numbering.size());
-
- std::map<gmshDofKey, int>::iterator
- itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itR != numbering.end()){
- std::map<gmshDofKey, int>::iterator
- itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
- typename std::map<gmshDofKey, scalar>::iterator
- itFF = fixedFULL.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itC != numbering.end() && itFF == fixedFULL.end()){
- lsys->addToMatrix(itR->second, itC->second, val);
- }
- else if (itFF != fixedFULL.end()){
- //printf("RHS = %g, ligne=%d \n", itFF->second, itR->second);
- lsys->addToMatrix(itR->second,itR->second, 1. );
- lsys->addToRightHandSide(itR->second, itFF->second);
- }
- else {
- typename std::map<gmshDofKey, scalar>::iterator
- itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
- if (itF != fixed.end()){
- //printf("RHS = val%g itF=%g \n", val, itF->second);
- lsys->addToRightHandSide(itR->second, -val*itF->second);
- }
- else{
- typename std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > >::
- iterator itConstrC = constraints.find(gmshDofKey(vC, iCompC, iFieldC));
- if (itConstrC != constraints.end()){
- for (unsigned int i = 0; i < itConstrC->second.size(); i++){
- gmshDofKey &dofKeyConstrC = itConstrC->second[i].first;
- scalar valConstrC = itConstrC->second[i].second;
- assemble(vR, iCompR, iFieldR,
- dofKeyConstrC.v, dofKeyConstrC.comp, dofKeyConstrC.field,
- val * valConstrC);
- }
- }
- }
- }
- }
- else{
- typename std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > >::
- iterator itConstrR = constraints.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itConstrR != constraints.end()){
- for (unsigned int i = 0; i < itConstrR->second.size(); i++){
- gmshDofKey &dofKeyConstrR = itConstrR->second[i].first;
- scalar valConstrR = itConstrR->second[i].second;
- assemble(dofKeyConstrR.v, dofKeyConstrR.comp, dofKeyConstrR.field,
- vC, iCompC, iFieldC,
- val * valConstrR);
- }
- }
- }
- }
- void assemble(MVertex *vR, int iCompR, int iFieldR, scalar val)
- {
- if (!lsys->isAllocated())lsys->allocate(numbering.size());
- std::map<gmshDofKey, int>::iterator
- itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itR != numbering.end()){
- lsys->addToRightHandSide(itR->second, val);
- }
- }
- int sizeOfR() const { return numbering.size(); }
- int sizeOfF() const { return fixed.size(); }
- void systemSolve(){
- lsys->systemSolve();
- }
- void systemClear(){
- lsys->zeroMatrix();
- lsys->zeroRightHandSide();
- }
-};
-
-#endif
diff --git a/Numeric/gmshConvexCombination.cpp b/Numeric/gmshConvexCombination.cpp
deleted file mode 100644
index 8f0f2d9315..0000000000
--- a/Numeric/gmshConvexCombination.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "gmshConvexCombination.h"
-#include "Numeric.h"
-
-void gmshConvexCombinationTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
-{
- m.set_all(0.);
- int nbNodes = e->getNumVertices();
- //double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
- //double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
- //double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
- const double _diff = 1.0;
-
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k < nbNodes; k++){
- m(j,k) = -1.*_diff;
- }
- m(j,j) = (nbNodes-1)*_diff;
- }
-
-}
diff --git a/Numeric/gmshConvexCombination.h b/Numeric/gmshConvexCombination.h
deleted file mode 100644
index d3fc41a948..0000000000
--- a/Numeric/gmshConvexCombination.h
+++ /dev/null
@@ -1,36 +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 _GMSH_CONVEX_H_
-#define _GMSH_CONVEX_H_
-
-#include "gmshTermOfFormulation.h"
-#include "simpleFunction.h"
-#include "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "fullMatrix.h"
-
-class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
- protected:
- const simpleFunction<double> *_diffusivity;
- const int _iField ;
- public:
- virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
- void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iField;
- }
- public:
- gmshConvexCombinationTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0) :
- gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField){}
- virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
-};
-
-
-#endif
diff --git a/Numeric/gmshCrossConf.cpp b/Numeric/gmshCrossConf.cpp
deleted file mode 100644
index f81d5d3562..0000000000
--- a/Numeric/gmshCrossConf.cpp
+++ /dev/null
@@ -1,58 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "gmshCrossConf.h"
-#include "Numeric.h"
-
-void gmshCrossConfTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
-{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
- int npts;
- IntPt *GP;
- double jac[3][3];
- double invjac[3][3];
- SVector3 Grads [256];
- double grads[256][3];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.set_all(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);
- SPoint3 p; e->pnt(u, v, w, p);
- const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
- inv3x3(jac, invjac) ;
- e->getGradShapeFunctions(u, v, w, grads);
- for (int j = 0; j < nbNodes; j++){
- Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
- invjac[0][2] * grads[j][2],
- invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
- invjac[1][2] * grads[j][2],
- invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
- invjac[2][2] * grads[j][2]);
- }
-
- SVector3 N (jac[2][0],jac[2][1],jac[2][2]);
- //printf("nx=%g, ny=%g, nz=%g norm =%g\n", N[0], N[1], N[2], sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2] ));
-
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k <= j; k++){
- m(j, k) += dot(crossprod(Grads[j],Grads[k]),N) * weight * detJ * _diff;
- }
- }
- }
-
- for (int j = 0; j < nbNodes; j++)
- for (int k = 0; k < j; k++)
- m(k, j) = -1.*m(j, k);
-}
-
diff --git a/Numeric/gmshCrossConf.h b/Numeric/gmshCrossConf.h
deleted file mode 100644
index dfdbe4d47c..0000000000
--- a/Numeric/gmshCrossConf.h
+++ /dev/null
@@ -1,56 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_CROSSCONF_H_
-#define _GMSH_CROSSCONF_H_
-
-#include "gmshTermOfFormulation.h"
-#include "simpleFunction.h"
-#include "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "fullMatrix.h"
-
-class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
- protected:
- const simpleFunction<double> *_diffusivity;
- const int _iField;
- int _iFieldB ;
- public:
- virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
- void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iField;
- }
- void getLocalDofC(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iFieldB;
- }
- public:
- gmshCrossConfTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) :
- gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
- {
- _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
- }
- virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
-};
-
-class gmshCrossConfTerm2DParametric : gmshCrossConfTerm {
- GFace *_gf;
- public:
- gmshCrossConfTerm2DParametric(GFace *gf, simpleFunction<double> *diffusivity, int iField = 0) :
- gmshCrossConfTerm(gf->model(), diffusivity, iField), _gf(gf) {}
- virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
-};
-
-
-#endif
diff --git a/Numeric/gmshDistance.cpp b/Numeric/gmshDistance.cpp
deleted file mode 100644
index 506d777210..0000000000
--- a/Numeric/gmshDistance.cpp
+++ /dev/null
@@ -1,85 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "gmshDistance.h"
-#include "Numeric.h"
-
-void gmshDistanceTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
-{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
- int npts;
- IntPt *GP;
- double jac[3][3];
- double invjac[3][3];
- double Grads[256][3], grads[256][3];
- double sf[256];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.set_all(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);
- SPoint3 p; e->pnt(u, v, w, p);
- const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
- inv3x3(jac, invjac) ;
- e->getShapeFunctions(u, v, w, sf);
- e->getGradShapeFunctions(u, v, w, grads);
- 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];
- }
-
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k <= j; k++){
- m(j, k) += _diff*(Grads[j][0] * Grads[k][0] +
- Grads[j][1] * Grads[k][1] +
- Grads[j][2] * Grads[k][2] ) * weight * detJ
- + (sf[j]*sf[k]) * weight * detJ ;
- }
- }
- }
-
- for (int j = 0; j < nbNodes; j++)
- for (int k = 0; k < j; k++)
- m(k, j) = m(j, k);
-}
-
-
-void gmshDistanceTerm::elementVector(MElement *e, fullVector<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) += ff[j] * weight * detJ;
- }
- }
-
-}
diff --git a/Numeric/gmshDistance.h b/Numeric/gmshDistance.h
deleted file mode 100644
index af55a6fc67..0000000000
--- a/Numeric/gmshDistance.h
+++ /dev/null
@@ -1,51 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_DISTANCE_H_
-#define _GMSH_DISTANCE_H_
-
-#include "gmshTermOfFormulation.h"
-#include "simpleFunction.h"
-#include "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "fullMatrix.h"
-
-class gmshDistanceTerm : public gmshNodalFemTerm<double> {
- protected:
- const simpleFunction<double> *_diffusivity;
- const int _iField;
- int _iFieldB ;
- public:
- virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
- void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iField;
- }
- void getLocalDofC(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iFieldB;
- }
- public:
- gmshDistanceTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) :
- gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
- {
- _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
- }
- void elementMatrix(MElement *e, fullMatrix<double> &m) const;
- void elementVector(MElement *e, fullVector<double> &m) const;
-
-};
-
-
-
-#endif
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
deleted file mode 100644
index 122b301211..0000000000
--- a/Numeric/gmshLaplace.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "gmshLaplace.h"
-#include "Numeric.h"
-
-void gmshLaplaceTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
-{
- int nbNodes = e->getNumVertices();
- int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
- int npts;
- IntPt *GP;
- double jac[3][3];
- double invjac[3][3];
- double Grads[256][3], grads[256][3];
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
- m.set_all(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);
- SPoint3 p; e->pnt(u, v, w, p);
- const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
- inv3x3(jac, invjac) ;
- e->getGradShapeFunctions(u, v, w, grads);
- 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];
- }
- double K_x=1.0;
- double K_y=1.0;
- double K_z=1.0;
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k <= j; k++){
- m(j, k) += (K_x*Grads[j][0] * Grads[k][0] +
- K_y*Grads[j][1] * Grads[k][1] +
- K_z*Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
- }
- }
- }
-
- for (int j = 0; j < nbNodes; j++)
- for (int k = 0; k < j; k++)
- m(k, j) = m(j, k);
-}
-
-
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
deleted file mode 100644
index cee00d054d..0000000000
--- a/Numeric/gmshLaplace.h
+++ /dev/null
@@ -1,49 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_LAPLACE_H_
-#define _GMSH_LAPLACE_H_
-
-#include "gmshTermOfFormulation.h"
-#include "simpleFunction.h"
-#include "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "fullMatrix.h"
-
-class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
- protected:
- const simpleFunction<double> *_diffusivity;
- const int _iField;
- int _iFieldB ;
- public:
- virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
- void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iField;
- }
- void getLocalDofC(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
- {
- *vR = e->getVertex(iRow);
- *iCompR = 0;
- *iFieldR = _iFieldB;
- }
- public:
- gmshLaplaceTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) :
- gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
- {
- _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
- }
- virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
-};
-
-
-
-#endif
diff --git a/Numeric/gmshLinearSystem.h b/Numeric/gmshLinearSystem.h
deleted file mode 100644
index b4be6df56b..0000000000
--- a/Numeric/gmshLinearSystem.h
+++ /dev/null
@@ -1,31 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_LINEAR_SYSTEM_H_
-#define _GMSH_LINEAR_SYSTEM_H_
-
-// A class that encapsulates a linear system solver interface :
-// building a sparse matrix, solving a linear system
-
-template <class scalar>
-class gmshLinearSystem {
- public :
- gmshLinearSystem (){}
- virtual ~gmshLinearSystem (){}
- virtual bool isAllocated() const = 0;
- virtual void allocate(int nbRows) = 0;
- virtual void clear() = 0;
- virtual void addToMatrix(int _row, int _col, scalar val) = 0;
- virtual scalar getFromMatrix(int _row, int _col) const = 0;
- virtual void addToRightHandSide(int _row, scalar val) = 0;
- 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/Numeric/gmshLinearSystemCSR.cpp b/Numeric/gmshLinearSystemCSR.cpp
deleted file mode 100644
index c2ddd50e6c..0000000000
--- a/Numeric/gmshLinearSystemCSR.cpp
+++ /dev/null
@@ -1,324 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include "GmshConfig.h"
-#include "GmshMessage.h"
-#include "gmshLinearSystemCSR.h"
-
-#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
-#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
-
-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);
- }
- 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;
- }
- }
-}
-
-
-static 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);
-}
-
-static void CSRList_Delete(CSRList_T *liste)
-{
- if (liste != 0) {
- free(liste->array);
- free(liste);
- }
-}
-
-template<>
-void gmshLinearSystemCSR<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);
-}
-
-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;
-
- 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 {
- 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];
-
- 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++){
- for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
- ptr[j] = j+1;
- }
- ptr[jptr[i+1]] = 0;
- }
-
- delete[] count;
-}
-
-#if defined(HAVE_GMM)
-#include "gmm.h"
-template<>
-int gmshLinearSystemCSRGmm<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
-
-#if defined(HAVE_GMM)
-#include "gmm.h"
-template<>
-int gmshLinearSystemCSRGmm<double>::checkSystem()
-{
- 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;
-
- printf("Coucou check system \n");
- for (unsigned int i = 0; i < _b->size(); i++)
- printf("%d ",((INDEX_TYPE *) jptr_->array)[i]);
- printf("\n");
-
- return 1;
-}
-#endif
-
diff --git a/Numeric/gmshLinearSystemCSR.h b/Numeric/gmshLinearSystemCSR.h
deleted file mode 100644
index 3e6b433584..0000000000
--- a/Numeric/gmshLinearSystemCSR.h
+++ /dev/null
@@ -1,151 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_LINEAR_SYSTEM_CSR_H_
-#define _GMSH_LINEAR_SYSTEM_CSR_H_
-
-#include <vector>
-#include "GmshConfig.h"
-#include "GmshMessage.h"
-#include "gmshLinearSystem.h"
-#include "linearSystemCSR.h"
-
-template <class scalar>
-class gmshLinearSystemCSR : public gmshLinearSystem<scalar> {
- protected:
- bool sorted;
- char *something;
- CSRList_T *a_,*ai_,*ptr_,*jptr_;
- std::vector<scalar> *_b, *_x;
- public:
- gmshLinearSystemCSR()
- : sorted(false), a_(0) {}
- virtual bool isAllocated() const { return a_ != 0; }
- virtual void allocate(int) ;
- virtual void clear()
- {
- allocate(0);
- }
- virtual ~gmshLinearSystemCSR()
- {
- 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 (il == 0) printf("NOT FOUND %d %d %d\n",il,ic,n);
-
- if(!something[il]) {
- jptr[il] = n;
- something[il] = 1;
- }
- else ptr[position_] = n;
-
-/* for (int i=0;i<_b->size();i++)printf("%d ",something[i]?jptr[i]:-1); */
-/* printf("\n"); */
-/* for (int i=0;i<CSRList_Nbr(ai_);i++)printf("(%d %d %g)",ai[i],ptr[i],a[i]); */
-/* printf("\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 gmshLinearSystemCSRGmm : public gmshLinearSystemCSR<scalar> {
- private:
- double _prec;
- int _noisy, _gmres;
- public:
- gmshLinearSystemCSRGmm()
- : _prec(1.e-8), _noisy(0), _gmres(0) {}
- virtual ~gmshLinearSystemCSRGmm()
- {}
- 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;
- }
-#endif
- virtual int checkSystem()
-#if defined(HAVE_GMM)
- ;
-#else
- {
- Msg::Error("Gmm++ is not available in this version of Gmsh");
- return 0;
- }
-#endif
-
-};
-
-#endif
diff --git a/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
deleted file mode 100644
index 55059b7a2e..0000000000
--- a/Numeric/gmshLinearSystemFull.h
+++ /dev/null
@@ -1,121 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_LINEAR_SYSTEM_FULL_H_
-#define _GMSH_LINEAR_SYSTEM_FULL_H_
-
-// Interface to a linear system with a FULL matrix
-
-#include "GmshMessage.h"
-#include "gmshLinearSystem.h"
-#include "fullMatrix.h"
-#include <stdlib.h>
-#include <set>
-
-template <class scalar>
-class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
- private:
- fullMatrix<scalar> *_a;
- fullVector<scalar> *_b, *_x;
- public :
- gmshLinearSystemFull() : _a(0), _b(0), _x(0){}
- virtual bool isAllocated() const { return _a != 0; }
- virtual void allocate(int _nbRows)
- {
- clear();
- _a = new fullMatrix<scalar>(_nbRows, _nbRows);
- _b = new fullVector<scalar>(_nbRows);
- _x = new fullVector<scalar>(_nbRows);
- }
-
- virtual ~gmshLinearSystemFull()
- {
- 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;
- }
- virtual int checkSystem()
- {
- //a->print("myMatrix");
- //_b->print("myVector");
-
- int nbNonConvex=0;
- for(int i = 0; i < _b->size(); i++){
-
- double diag = (*_a)(i, i);
- double offDiag = 0.0;
- bool convex = true;
- std::set<int> Ni;
- for(int j = 0; j < _b->size(); j++){
- if ( (j !=i) && (*_a)(i, j) > 0.0 ) convex = false;
- if ( j !=i && (*_a)(i, j) != 0.0){
- offDiag += (*_a)(i, j);
- Ni.insert(j);
- }
- }
- if (fabs(offDiag+diag) > 1.e-10 && fabs(offDiag) > 1.e-10 ) convex= false;
-
- if (convex == false){
- nbNonConvex+=1;
- printf("*** WARNING NON CONVEX LINE !!!!!\n");
- printf("*** i=%d : diag=%g offDiag=%g diff=%g convex=%s size Ni=%d\n", i , diag, offDiag, fabs(offDiag+diag), (convex)?"true":"false", Ni.size());
- for (std::set<int>::iterator it = Ni.begin(); it != Ni.end(); ++it) (*_a)(i, *it) = -1.0;
- (*_a)(i, i) = Ni.size();
- }
- }
-
- printf("nonConvex=%d total=%d ratio=%g\n ",nbNonConvex, _b->size(), nbNonConvex/(double)_b->size());
- //_a->print("myMatrix AFTER !!! ");
-
- return 1;
- }
-};
-
-#endif
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
deleted file mode 100644
index 372e1d605b..0000000000
--- a/Numeric/gmshLinearSystemGmm.h
+++ /dev/null
@@ -1,123 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_LINEAR_SYSTEM_GMM_H_
-#define _GMSH_LINEAR_SYSTEM_GMM_H_
-
-// Interface to GMM++
-
-#include "GmshConfig.h"
-#include "GmshMessage.h"
-#include "gmshLinearSystem.h"
-
-#if defined(HAVE_GMM)
-
-#include <gmm.h>
-
-template <class scalar>
-class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
- private:
- gmm::row_matrix<gmm::wsvector<scalar> > *_a;
- std::vector<scalar> *_b, *_x;
- double _prec;
- int _noisy, _gmres;
- public:
- gmshLinearSystemGmm()
- : _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 ~gmshLinearSystemGmm()
- {
- 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
-
-template <class scalar>
-class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
-public :
- gmshLinearSystemGmm()
- {
- 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/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
deleted file mode 100644
index 721413c5e3..0000000000
--- a/Numeric/gmshTermOfFormulation.h
+++ /dev/null
@@ -1,188 +0,0 @@
-// 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
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _GMSH_TERM_OF_FORMULATION_H_
-#define _GMSH_TERM_OF_FORMULATION_H_
-
-#include <math.h>
-#include <map>
-#include <vector>
-#include "fullMatrix.h"
-#include "simpleFunction.h"
-#include "gmshAssembler.h"
-#include "GModel.h"
-#include "MElement.h"
-
-template<class scalar>
-class gmshTermOfFormulation {
- protected:
- GModel *_gm;
- public:
- gmshTermOfFormulation(GModel *gm) : _gm(gm) {}
- virtual ~gmshTermOfFormulation(){}
- virtual void addToMatrix(gmshAssembler<scalar> &) const = 0;
-};
-
-// a nodal finite element term : variables are always defined at nodes
-// of the mesh
-template<class scalar>
-class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
- 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 key associated to a given row
- // of the local element matrix
- virtual void getLocalDofR(MElement *e, int iRow, MVertex **vR,
- int *iCompR, int *iFieldR) const = 0;
- virtual void getLocalDofC(MElement *e, int iCol, MVertex **vC,
- int *iCompC, int *iFieldC) const
- {
- getLocalDofR(e, iCol, vC, iCompC, iFieldC);
- }
- public:
- gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation<scalar>(gm) {}
- virtual ~gmshNodalFemTerm (){}
- virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const = 0;
- virtual void elementVector(MElement *e, fullVector<scalar> &m) const {
- m.scale(0.0);
- }
- void addToMatrix(gmshAssembler<scalar> &lsys) const
- {
- GModel *m = gmshTermOfFormulation<scalar>::_gm;
- if (m->getNumRegions()){
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
- addToMatrix(lsys, *it);
- }
- }
- else if(m->getNumFaces()){
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- addToMatrix(lsys, *it);
- }
- }
- }
- void addToMatrix(gmshAssembler<scalar> &lsys, GEntity *ge) const
- {
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- addToMatrix(lsys, e);
- }
- }
- void addToMatrix(gmshAssembler<scalar> &lsys, MElement *e) const
- {
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- fullMatrix<scalar> localMatrix(nbR, nbC);
- elementMatrix(e, localMatrix);
- addToMatrix(lsys, localMatrix, e);
- }
- void addToMatrix(gmshAssembler<scalar> &lsys, const std::vector<MElement*> &v) const
- {
- for (unsigned int i = 0; i < v.size(); i++)
- addToMatrix(lsys, v[i]);
- }
- void addToMatrix(gmshAssembler<scalar> &lsys, fullMatrix<scalar> &localMatrix,
- MElement *e) const
- {
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- for (int j = 0; j < nbR; j++){
- MVertex *vR;
- int iCompR, iFieldR;
- getLocalDofR(e, j, &vR, &iCompR, &iFieldR);
- for (int k = 0; k < nbC; k++){
- MVertex *vC;
- int iCompC, iFieldC;
- getLocalDofC(e, k, &vC, &iCompC, &iFieldC);
- lsys.assemble(vR, iCompR, iFieldR, vC, iCompC, iFieldC, localMatrix(j, k));
- }
- }
- }
- void addDirichlet(int physical,
- int dim,
- int comp,
- int field,
- const simpleFunction<scalar> &e,
- gmshAssembler<scalar> &lsys)
- {
- std::vector<MVertex *> v;
- GModel *m = gmshTermOfFormulation<scalar>::_gm;
- m->getMeshVerticesForPhysicalGroup(dim, physical, v);
- for (unsigned int i = 0; i < v.size(); i++)
- lsys.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
- }
- void addNeumann(int physical,
- int dim,
- int comp,
- int field,
- const simpleFunction<scalar> &fct,
- gmshAssembler<scalar> &lsys)
- {
- std::map<int, std::vector<GEntity*> > groups[4];
- GModel *m = gmshTermOfFormulation<scalar>::_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 scalar FCT = fct(p.x(), p.y(), p.z());
- for (int k = 0; k < nbNodes; k++){
- lsys.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
- }
- }
- }
- }
- }
- void addToRightHandSide (gmshAssembler<scalar> &lsys) const {
- GModel *m = gmshTermOfFormulation<scalar>::_gm;
- if (m->getNumRegions()){
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
- addToRightHandSide(lsys,*it);
- }
- }
- else if(m->getNumFaces()){
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- addToRightHandSide(lsys,*it);
- }
- }
- }
- void addToRightHandSide (gmshAssembler<scalar> &lsys, GEntity *ge) const {
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement (i);
- int nbR = sizeOfR(e);
- fullVector<scalar> V (nbR);
- elementVector (e, V);
- // assembly
- for (int j=0;j<nbR;j++){
- MVertex *vR;int iCompR,iFieldR;
- getLocalDofR (e,j,&vR,&iCompR,&iFieldR);
- lsys.assemble(vR,iCompR,iFieldR,V(j));
- }
- }
- }
-};
-
-#endif
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index 57ef4af802..1f77c8137a 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -28,7 +28,7 @@ class convexCombinationTerm : public femTerm<double, double> {
{
return se->getMeshElement()->getNumVertices();
}
- void getLocalDofR(SElement *se, int iRow) const
+ Dof getLocalDofR(SElement *se, int iRow) const
{
return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iField);
}
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index 1565addd4e..c8d313d7cf 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -36,11 +36,11 @@ class crossConfTerm : public femTerm<double, double> {
{
return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldR);
}
- void getLocalDofC(SElement *se, int iRow) const
+ Dof getLocalDofC(SElement *se, int iRow) const
{
return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldC);
}
- virtual void elementMatrix(SElement *e, fullMatrix<double> &m) const
+ virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
{
MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
@@ -82,7 +82,7 @@ class crossConfTerm : public femTerm<double, double> {
}
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k < j; k++)
- m(k, j) = -1.*m(j, k);
+ m(k, j) = -1. * m(j, k);
}
};
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index fa709a8a4b..2942f9901a 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -17,14 +17,14 @@
// \nabla \cdot k \nabla U + a U
template<class scalar>
-class helmoltzTerm : public femTerm<scalar, scalar> {
+class helmholtzTerm : public femTerm<scalar, scalar> {
protected:
const simpleFunction<scalar> *_k, *_a;
const int _iFieldR;
int _iFieldC ;
public:
- helmoltzTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k,
- simpleFunction<scalar> *a)
+ helmholtzTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k,
+ simpleFunction<scalar> *a)
: femTerm<scalar, scalar>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
_k(k), _a(a) {}
// one dof per vertex (nodal fem)
@@ -38,11 +38,13 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldR);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ Dof::createTypeWithTwoInts(0, _iFieldR));
}
Dof getLocalDofC(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldC);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ Dof::createTypeWithTwoInts(0, _iFieldC));
}
virtual void elementMatrix(SElement *se, fullMatrix<scalar> &m) const
{
@@ -50,7 +52,7 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
// compute integration rule
const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
2 * (e->getPolynomialOrder() - 1);
- int npts;IntPt *GP;
+ int npts; IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP);
// get the number of nodes
const int nbNodes = e->getNumVertices();
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index f712884401..f97c8d0bfc 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -9,24 +9,11 @@
#include "helmholtzTerm.h"
// \nabla \cdot k \nabla U
-template<class scalar>
+template<class scalar>
class laplaceTerm : public helmholtzTerm<scalar> {
- simpleFunction<scalar> *ONE;
public:
- laplaceTerm(GModel *gm, int iFieldR) :
- helmholtzTerm<scalar>(gm, iFieldR, iFieldR, 0, 0)
- {
- ONE = new simpleFunction<scalar>(1.0);
- _k = ONE;
- }
- laplaceTerm(GModel *gm, int iFieldR, simpleFunction<scalar> *k) :
- helmholtzTerm<scalar>(gm, iFieldR, iFieldR, k, 0), ONE(0) {}
- laplaceTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k) :
- helmholtzTerm<scalar>(gm, iFieldR, iFieldC, k, 0), ONE(0) {}
- virtual ~laplaceTerm()
- {
- if(ONE) delete ONE;
- }
+ laplaceTerm(GModel *gm, int iFieldR, simpleFunction<scalar> *k)
+ : helmholtzTerm<scalar>(gm, iFieldR, iFieldR, k, 0) {}
};
#endif
--
GitLab