Skip to content
Snippets Groups Projects
Commit 2c463548 authored by Gauthier Becker's avatar Gauthier Becker
Browse files

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)
parent b6a93443
Branches
Tags
No related merge requests found
......@@ -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
......@@ -27,6 +27,7 @@ class linearSystem {
virtual void zeroRightHandSide() = 0;
virtual int systemSolve() = 0;
static void registerBindings (binding*);
virtual double normInfRightHandSide() const = 0;
};
#endif
......@@ -113,6 +113,17 @@ class linearSystemCSR : 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;
}
};
template <class scalar>
......
......@@ -70,6 +70,16 @@ class linearSystemFull : public linearSystem<scalar> {
{
for(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;
}
virtual int systemSolve()
{
if (_b->size())
......
......@@ -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; }
......
......@@ -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)
{
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment