From 2c463548fca87084024e9a4c797d81fd0e325369 Mon Sep 17 00:00:00 2001
From: Gauthier Becker <gauthierbecker@gmail.com>
Date: Fri, 11 Jun 2010 15:18:28 +0000
Subject: [PATCH] Adds in dgplate : 1) Non linear solver (Newton-Raphson) 2)
cohesive law (fracture in bending OK) 3) lua version of code 4) a test case
with 3rd order quadrilateral element. Some modifications are made in
linearSystem (add a function to compute norm of RightHandSide) and in
dofManager (add a function to have the fixed dof outside dofManager)
---
Geo/MElement.cpp | 8 +++----
Post/PViewDataGModel.cpp | 46 +++++++++++++++++++-------------------
Solver/dofManager.h | 15 +++++++++----
Solver/linearSystem.h | 1 +
Solver/linearSystemCSR.h | 45 +++++++++++++++++++++++--------------
Solver/linearSystemFull.h | 18 +++++++++++----
Solver/linearSystemGMM.h | 17 +++++++++++---
Solver/linearSystemPETSc.h | 10 +++++++--
8 files changed, 103 insertions(+), 57 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index b62648a362..f885578e7d 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -156,7 +156,7 @@ int MElement::getVolumeSign()
}
bool MElement::setVolumePositive()
-{
+{
int s = getVolumeSign();
if(s < 0) revert();
if(!s) return false;
@@ -438,7 +438,7 @@ void MElement::interpolateCurl(double val[], double u, double v, double w, doubl
f[2] = fy[0] - fx[1];
}
-double MElement::interpolateDiv(double val[], double u, double v, double w,
+double MElement::interpolateDiv(double val[], double u, double v, double w,
int stride, int order)
{
double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
@@ -477,7 +477,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
fprintf(fp, " %d %d %d 1 %d", 4 + par, abs(physical), elementary, _partition);
else{
int numGhosts = ghosts->size();
- fprintf(fp, " %d %d %d %d %d", 4 + numGhosts + par, abs(physical), elementary,
+ fprintf(fp, " %d %d %d %d %d", 4 + numGhosts + par, abs(physical), elementary,
1 + numGhosts, _partition);
for(unsigned int i = 0; i < ghosts->size(); i++)
fprintf(fp, " %d", -(*ghosts)[i]);
@@ -501,7 +501,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
// tags change from element to element (third-party codes can
// still write MSH file optimized for reading speed, by grouping
// elements with the same number of tags in blobs)
- int blob[60] = {type, 1, numTags, num ? num : _num, abs(physical), elementary,
+ int blob[60] = {type, 1, numTags, num ? num : _num, abs(physical), elementary,
1 + numGhosts, _partition};
if(ghosts)
for(int i = 0; i < numGhosts; i++) blob[8 + i] = -(*ghosts)[i];
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 94fae4dda8..8b14f91532 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -15,7 +15,7 @@
#include "Numeric.h"
#include "GmshMessage.h"
-PViewDataGModel::PViewDataGModel(DataType type)
+PViewDataGModel::PViewDataGModel(DataType type)
: PViewData(), _min(VAL_INF), _max(-VAL_INF), _type(type)
{
}
@@ -69,25 +69,25 @@ bool PViewDataGModel::finalize(bool computeMinMax)
for(int step = 0; step < getNumTimeSteps(); step++){
GModel *m = _steps[step]->getModel();
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
- if((*it)->lines.size())
+ if((*it)->lines.size())
_addInterpolationMatricesForElement((*it)->lines[0]);
}
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
- if((*it)->triangles.size())
+ if((*it)->triangles.size())
_addInterpolationMatricesForElement((*it)->triangles[0]);
- if((*it)->quadrangles.size())
+ if((*it)->quadrangles.size())
_addInterpolationMatricesForElement((*it)->quadrangles[0]);
if((*it)->polygons.size())
_addInterpolationMatricesForElement((*it)->polygons[0]);
}
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
- if((*it)->tetrahedra.size())
+ if((*it)->tetrahedra.size())
_addInterpolationMatricesForElement((*it)->tetrahedra[0]);
- if((*it)->hexahedra.size())
+ if((*it)->hexahedra.size())
_addInterpolationMatricesForElement((*it)->hexahedra[0]);
- if((*it)->prisms.size())
+ if((*it)->prisms.size())
_addInterpolationMatricesForElement((*it)->prisms[0]);
- if((*it)->pyramids.size())
+ if((*it)->pyramids.size())
_addInterpolationMatricesForElement((*it)->pyramids[0]);
if((*it)->polyhedra.size())
_addInterpolationMatricesForElement((*it)->polyhedra[0]);
@@ -107,7 +107,7 @@ void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
fs->coefficients, fs->monomials);
else
setInterpolationMatrices(type, fs->coefficients, fs->monomials);
- }
+ }
}
MElement *PViewDataGModel::_getElement(int step, int ent, int ele)
@@ -149,7 +149,7 @@ double PViewDataGModel::getMax(int step)
}
SBoundingBox3d PViewDataGModel::getBoundingBox(int step)
-{
+{
if(step < 0 || _steps.empty()){
SBoundingBox3d tmp;
for(unsigned int i = 0; i < _steps.size(); i++)
@@ -297,7 +297,7 @@ int PViewDataGModel::getNumElements(int step, int ent)
// to generalize
if(step < 0 && ent < 0) return _steps[0]->getModel()->getNumMeshElements();
if(step < 0) return _steps[0]->getEntity(ent)->getNumMeshElements();
- if(ent < 0) return _steps[step]->getModel()->getNumMeshElements();
+ if(ent < 0) return _steps[step]->getModel()->getNumMeshElements();
return _steps[step]->getEntity(ent)->getNumMeshElements();
}
@@ -320,11 +320,11 @@ int PViewDataGModel::getNumNodes(int step, int ent, int ele)
}
}
-int PViewDataGModel::getNode(int step, int ent, int ele, int nod,
+int PViewDataGModel::getNode(int step, int ent, int ele, int nod,
double &x, double &y, double &z)
{
MElement *e = _getElement(step, ent, ele);
- if(_type == GaussPointData){
+ if(_type == GaussPointData){
std::vector<double> &p(_steps[step]->getGaussPoints(e->getTypeForMSH()));
if(p[0] == 1.e22){
// hack: the points are the element vertices
@@ -354,7 +354,7 @@ int PViewDataGModel::getNode(int step, int ent, int ele, int nod,
}
}
-void PViewDataGModel::setNode(int step, int ent, int ele, int nod,
+void PViewDataGModel::setNode(int step, int ent, int ele, int nod,
double x, double y, double z)
{
MVertex *v = _getElement(step, ent, ele)->getVertex(nod);
@@ -409,7 +409,7 @@ void PViewDataGModel::getValue(int step, int ent, int ele, int nod, int comp, do
{
MElement *e = _getElement(step, ent, ele);
switch(_type){
- case NodeData:
+ case NodeData:
val = _steps[step]->getData(e->getVertex(nod)->getNum())[comp];
break;
case ElementNodeData:
@@ -417,7 +417,7 @@ void PViewDataGModel::getValue(int step, int ent, int ele, int nod, int comp, do
val = _steps[step]->getData(e->getNum())[_steps[step]->getNumComponents() * nod + comp];
break;
case ElementData:
- default:
+ default:
val = _steps[step]->getData(e->getNum())[comp];
break;
}
@@ -427,22 +427,22 @@ void PViewDataGModel::setValue(int step, int ent, int ele, int nod, int comp, do
{
MElement *e = _getElement(step, ent, ele);
switch(_type){
- case NodeData:
+ case NodeData:
_steps[step]->getData(e->getVertex(nod)->getNum())[comp] = val;
break;
case ElementNodeData:
- case GaussPointData:
+ case GaussPointData:
_steps[step]->getData(e->getNum())[_steps[step]->getNumComponents() * nod + comp] = val;
break;
- case ElementData:
- default:
+ case ElementData:
+ default:
_steps[step]->getData(e->getNum())[comp] = val;
break;
}
}
int PViewDataGModel::getNumEdges(int step, int ent, int ele)
-{
+{
return _getElement(step, ent, ele)->getNumEdges();
}
@@ -511,7 +511,7 @@ bool PViewDataGModel::combineTime(nameData &nd)
}
// copy interpolation matrices
- for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it =
+ for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it =
data[0]->_interpolation.begin(); it != data[0]->_interpolation.end(); it++)
if(_interpolation[it->first].empty())
for(unsigned int i = 0; i < it->second.size(); i++)
@@ -522,7 +522,7 @@ bool PViewDataGModel::combineTime(nameData &nd)
for(unsigned int j = 0; j < data[i]->_steps.size(); j++)
if(data[i]->hasTimeStep(j))
_steps.push_back(new stepData<double>(*data[i]->_steps[j]));
-
+
std::string tmp;
if(nd.name == "__all__")
tmp = "all";
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index e2a4f4886f..d4682c6915 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -82,7 +82,7 @@ class dofManager{
public:
typedef typename dofTraits<T>::VecType dataVec;
typedef typename dofTraits<T>::MatType dataMat;
- private:
+ protected:
// general affine constraint on sub-blocks, treated by adding
// equations:
@@ -301,7 +301,7 @@ class dofManager{
}
}
- inline void assemble(std::vector<Dof> &R, const fullVector<dataMat> &m) // for linear forms
+ virtual inline void assemble(std::vector<Dof> &R, const fullVector<dataMat> &m) // for linear forms
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
std::vector<int> NR(R.size());
@@ -335,7 +335,7 @@ class dofManager{
}
- inline void assemble(std::vector<Dof> &R, const fullMatrix<dataMat> &m)
+ virtual inline void assemble(std::vector<Dof> &R, const fullMatrix<dataMat> &m)
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
std::vector<int> NR(R.size());
@@ -487,6 +487,13 @@ class dofManager{
}
}
}
-
+ void getFixedDof(std::vector<Dof> &R){
+ R.reserve(fixed.size());
+ std::map<Dof, double>::iterator it; // template problem ?? //TODO
+ for(it=fixed.begin(); it!=fixed.end();++it){
+ R.push_back((*it).first);
+ }
+ }
};
+
#endif
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index 5ff9c10a14..bec6ac2bd4 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -27,6 +27,7 @@ class linearSystem {
virtual void zeroRightHandSide() = 0;
virtual int systemSolve() = 0;
static void registerBindings (binding*);
+ virtual double normInfRightHandSide() const = 0;
};
#endif
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 5e9f325a35..a83561e326 100644
--- a/Solver/linearSystemCSR.h
+++ b/Solver/linearSystemCSR.h
@@ -13,7 +13,7 @@
class binding;
-typedef int INDEX_TYPE ;
+typedef int INDEX_TYPE ;
typedef struct {
int nmax;
int size;
@@ -22,7 +22,7 @@ typedef struct {
int isorder;
char *array;
} CSRList_T;
-
+
void CSRList_Add(CSRList_T *liste, const void *data);
int CSRList_Nbr(CSRList_T *liste);
@@ -31,7 +31,7 @@ class linearSystemCSR : public linearSystem<scalar> {
protected:
bool sorted;
char *something;
- CSRList_T *_a, *_ai, *_ptr, *_jptr;
+ CSRList_T *_a, *_ai, *_ptr, *_jptr;
std::vector<scalar> *_b, *_x;
public:
linearSystemCSR()
@@ -52,9 +52,9 @@ class linearSystemCSR : public linearSystem<scalar> {
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){
@@ -64,24 +64,24 @@ class linearSystemCSR : public linearSystem<scalar> {
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
-
+ // 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;
+ something[il] = 1;
}
else ptr[position] = n;
}
@@ -91,11 +91,11 @@ class linearSystemCSR : public linearSystem<scalar> {
{
Msg::Error("getFromMatrix not implemented for CSR");
}
- virtual void addToRightHandSide(int row, const scalar &val)
+ virtual void addToRightHandSide(int row, const scalar &val)
{
if(val != 0.0) (*_b)[row] += val;
}
- virtual void getFromRightHandSide(int row, scalar &val) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
val = (*_b)[row];
}
@@ -109,10 +109,21 @@ class linearSystemCSR : public linearSystem<scalar> {
scalar *a = (scalar*) _a->array;
for (int i = 0; i < N; i++) a[i] = 0;
}
- virtual void zeroRightHandSide()
+ virtual void zeroRightHandSide()
{
for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
}
+ virtual double normInfRightHandSide() const{
+ double nor = 0.;
+ double temp;
+ for(int i=0;i<_b->size();i++){
+ temp = (*_b)[i];
+ if(temp<0) temp = -temp;
+ if(nor<temp) nor=temp;
+
+ }
+ return nor;
+ }
};
template <class scalar>
@@ -126,7 +137,7 @@ class linearSystemCSRGmm : public linearSystemCSR<scalar> {
void setPrec(double p){ _prec = p; }
void setNoisy(int n){ _noisy = n; }
void setGmres(int n){ _gmres = n; }
- virtual int systemSolve()
+ virtual int systemSolve()
#if !defined(HAVE_GMM)
{
Msg::Error("Gmm++ is not available in this version of Gmsh");
@@ -147,7 +158,7 @@ class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
if (il <= ic)
linearSystemCSR<scalar>::addToMatrix(il, ic, val);
}
- virtual int systemSolve()
+ virtual int systemSolve()
#if !defined(HAVE_TAUCS)
{
Msg::Error("TAUCS is not available in this version of Gmsh");
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index 1c3cf865e4..eb8c17839b 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -54,15 +54,15 @@ class linearSystemFull : public linearSystem<scalar> {
{
if(val != 0.0) (*_b)(row) += val;
}
- virtual void getFromRightHandSide(int row, scalar &val) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
val = (*_b)(row);
}
- virtual void getFromSolution(int row, scalar &val) const
+ virtual void getFromSolution(int row, scalar &val) const
{
val = (*_x)(row);
}
- virtual void zeroMatrix()
+ virtual void zeroMatrix()
{
_a->setAll(0.);
}
@@ -70,7 +70,17 @@ class linearSystemFull : public linearSystem<scalar> {
{
for(int i = 0; i < _b->size(); i++) (*_b)(i) = 0.;
}
- virtual int systemSolve()
+ virtual double normInfRightHandSide() const{
+ double nor = 0.;
+ double temp;
+ for(int i=0;i<_b->size();i++){
+ temp = (*_b)(i);
+ if(temp<0) temp = -temp;
+ if(nor<temp) nor=temp;
+ }
+ return nor;
+ }
+ virtual int systemSolve()
{
if (_b->size())
_a->luSolve(*_b, *_x);
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
index 6198561a6a..8d8fabf5e6 100644
--- a/Solver/linearSystemGMM.h
+++ b/Solver/linearSystemGMM.h
@@ -47,7 +47,7 @@ class linearSystemGmm : public linearSystem<scalar> {
_a = 0;
}
- virtual void addToMatrix(int row, int col, const scalar &val)
+ virtual void addToMatrix(int row, int col, const scalar &val)
{
if(val != 0.0) (*_a)(row, col) += val;
}
@@ -56,12 +56,12 @@ class linearSystemGmm : public linearSystem<scalar> {
val = (*_a)(row, col);
}
- virtual void addToRightHandSide(int row, const scalar &val)
+ virtual void addToRightHandSide(int row, const scalar &val)
{
if(val != 0.0) (*_b)[row] += val;
}
- virtual void getFromRightHandSide(int row, scalar &val) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
val = (*_b)[row];
}
@@ -77,6 +77,17 @@ class linearSystemGmm : public linearSystem<scalar> {
{
for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
}
+ virtual double normInfRightHandSide() const {
+ double nor = 0.;
+ double temp;
+ for(int i=0;i<_b->size();i++){
+ temp = (*_b)[i];
+ if(temp<0) temp = -temp;
+ if(nor<temp) nor=temp;
+ }
+ return nor;
+ }
+
void setPrec(double p){ _prec = p; }
void setNoisy(int n){ _noisy = n; }
void setGmres(int n){ _gmres = n; }
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 4097a529bd..32cc3b0d1c 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -56,7 +56,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
virtual void allocate(int nbRows)
{
clear();
- _try(MatCreate(PETSC_COMM_WORLD, &_a));
+ _try(MatCreate(PETSC_COMM_WORLD, &_a));
_try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows, nbRows));
// override the default options with the ones from the option
// database (if any)
@@ -105,6 +105,12 @@ class linearSystemPETSc : public linearSystem<scalar> {
#else
val = s;
#endif
+ }
+ virtual double normInfRightHandSide() const {
+ PetscScalar nor;
+ _try(VecNorm(_b,NORM_INFINITY,&nor));
+ return nor;
+
}
virtual void addToMatrix(int row, int col, const scalar &val)
{
@@ -185,7 +191,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
virtual void zeroMatrix() {}
virtual void zeroRightHandSide() {}
virtual int systemSolve() { return 0; }
-};
+};
#endif
--
GitLab