From 8a29a1a70c730b2dcb138847c323ab7d6d43eca7 Mon Sep 17 00:00:00 2001
From: Eric Bechet <eric.bechet@ulg.ac.be>
Date: Mon, 4 Jan 2010 20:22:58 +0000
Subject: [PATCH] less template, more inheritance
---
Solver/elasticitySolver.cpp | 75 +++++++----
Solver/elasticitySolver.h | 5 +-
Solver/materialLaw.h | 25 ++++
Solver/solverAlgorithms.h | 24 ++++
Solver/terms.h | 254 +++++++++++++++++++-----------------
5 files changed, 237 insertions(+), 146 deletions(-)
create mode 100644 Solver/materialLaw.h
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index b69c5ef747..35051215ff 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -188,7 +188,7 @@ void elasticitySolver::solve()
for (unsigned int i = 0; i < allNeumann.size(); i++)
{
- LoadTerm<FunctionSpace<SVector3> > Lterm(*LagSpace,allNeumann[i]._f);
+ LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
if (allNeumann[i].onWhat==BoundaryCondition::ON_VERTEX)
Assemble(Lterm,*LagSpace,allNeumann[i].g->vbegin(),allNeumann[i].g->vend(),*pAssembler);
else
@@ -200,13 +200,23 @@ void elasticitySolver::solve()
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
- IsotropicElasticTerm<FunctionSpace<SVector3> ,FunctionSpace<SVector3> > Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
+ IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
+// LaplaceTerm<SVector3,SVector3> Eterm(*LagSpace);
Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
}
printf("-- done assembling!\n");
lsys->systemSolve();
printf("-- done solving!\n");
+ double energ=0;
+ for (unsigned int i = 0; i < elasticFields.size(); i++)
+ {
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+ }
+ printf("elastic energy=%f\n",energ);
}
@@ -264,21 +274,12 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
}
}
- std::map<int, std::vector<double> > data;
+ std::map<int, std::vector<double> > data;
+ SolverField<SVector3> Field(pAssembler, LagSpace);
for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
{
- SVector3 val(0,0,0);
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
- std::vector<Dof> D;
- std::vector<SVector3> SFVals;
- std::vector<double> Vals;
- LagSpace->getKeys(*it,D);
- pAssembler->getDofValue(D,Vals);
- LagSpace->f(*it,SFVals);
- for (int i=0;i<D.size();++i)
- val+=SFVals[i]*Vals[i];
- }
+ SVector3 val;
+ Field.f(*it,val);
std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
data[(*it)->getNum()]=vec;
}
@@ -290,30 +291,54 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
{
std::map<int, std::vector<double> > data;
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- double energ=0;
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
SolverField<SVector3> Field(pAssembler, LagSpace);
- IsotropicElasticTerm<SolverField<SVector3> ,SolverField<SVector3> > Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- ScalarTermOne One;
+ IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ ScalarTermConstant One(1.0);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
{
MElement *e=*it;
- fullMatrix<double> localMatrix;
+ double energ;
+ double vol;
IntPt *GP;
int npts=Integ_Bulk.getIntPoints(e,&GP);
- Eterm.get(e,npts,GP,localMatrix);
- double vol=0;
+ Elastic_Energy_Term.get(e,npts,GP,energ);
One.get(e,npts,GP,vol);
std::vector<double> vec;
- vec.push_back(localMatrix(0,0)/vol);
- energ+=localMatrix(0,0);
+ vec.push_back(energ/vol);
data[(*it)->getNum()]=vec;
}
}
- std::cout<< "elastic energy=" << energ << std::endl;
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
- return pv;
+ return pv;
+}
+
+PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+ std::map<int, std::vector<double> > data;
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+ {
+ MElement *e=*it;
+ double energ;
+ double vol;
+ IntPt *GP;
+ int npts=Integ_Bulk.getIntPoints(e,&GP);
+ Elastic_Energy_Term.get(e,npts,GP,energ);
+ std::vector<double> vec;
+ vec.push_back(energ);
+ data[(*it)->getNum()]=vec;
+ }
+ }
+ PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
+ return pv;
}
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index ef5109181f..68ada04d27 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -20,9 +20,8 @@ class groupOfElements;
struct elasticField {
int _tag; // tag for the dofManager
groupOfElements *g; // support for this field
- simpleFunction<double> *_enrichment; // XFEM enrichment
double _E, _nu; // specific elastic datas (should be somewhere else)
- elasticField () : g(0), _enrichment(0),_tag(0){}
+ elasticField () : g(0),_tag(0){}
};
struct BoundaryCondition
@@ -75,7 +74,7 @@ class elasticitySolver
virtual void solve();
virtual PView *buildDisplacementView(const std::string &postFileName);
virtual PView *buildElasticEnergyView(const std::string &postFileName);
- // virtual PView *buildVonMisesView(const std::string &postFileName);
+ virtual PView *buildVonMisesView(const std::string &postFileName);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
// std::pair<PView *, PView*> buildErrorEstimateView
diff --git a/Solver/materialLaw.h b/Solver/materialLaw.h
new file mode 100644
index 0000000000..05863c58b1
--- /dev/null
+++ b/Solver/materialLaw.h
@@ -0,0 +1,25 @@
+//
+// C++ Interface: materialLaw
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _MATERIALLAW_H_
+#define _MATERIALLAW_H_
+
+class Material
+{
+ public:
+ virtual ~Material() {}
+
+};
+
+
+
+#endif //_MATERIALLAW_H_
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 0afc37ef51..1e0cdc5590 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -90,6 +90,30 @@ template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &
assembler.assemble(R, localVector);
}
+template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,Iterator itbegin,Iterator itend,QuadratureBase &integrator,dataMat & val)
+{
+ dataMat localval;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MElement *e = *it;
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localval);
+ val+=localval;
+ }
+}
+
+template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,MElement *e,QuadratureBase &integrator,dataMat & val)
+{
+ dataMat localval;
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localval);
+ val+=localval;
+}
+
+
+
template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
{
int nbff=dofs.size();
diff --git a/Solver/terms.h b/Solver/terms.h
index dff4c69fd0..de33311651 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -20,51 +20,77 @@
#include "Numeric.h"
#include "functionSpace.h"
#include "groupOfElements.h"
+#include "materialLaw.h"
class BilinearTermBase
{
public :
+ virtual ~BilinearTermBase() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
};
-template<class S1,class S2> class BilinearTerm : public BilinearTermBase
+template<class T1,class T2> class BilinearTerm : public BilinearTermBase
{
protected :
- S1& space1;
- S2& space2;
+ FunctionSpace<T1>& space1;
+ FunctionSpace<T2>& space2;
public :
- BilinearTerm(S1& space1_,S2& space2_) : space1(space1_),space2(space2_) {}
+ BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T1>& space2_) : space1(space1_),space2(space2_) {}
+ virtual ~BilinearTerm() {}
};
class LinearTermBase
{
public:
+ virtual ~LinearTermBase() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
virtual void get(MVertex *ver,fullVector<double> &m) =0;
};
-template<class S1> class LinearTerm : public LinearTermBase
+template<class T1> class LinearTerm : public LinearTermBase
{
protected :
- S1& space1;
- public :
- LinearTerm(S1& space1_) : space1(space1_) {}
+ FunctionSpace<T1>& space1;
+ public :
+ LinearTerm(FunctionSpace<T1>& space1_) : space1(space1_) {}
+ virtual ~LinearTerm() {}
};
class ScalarTermBase
{
public :
+ virtual ~ScalarTermBase() {}
virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
};
class ScalarTerm : public ScalarTermBase
{
- public :
+ public :
+ virtual ~ScalarTerm() {}
};
-class ScalarTermOne : public ScalarTerm
+
+template<class T1,class T2> class BilinearTermToScalarTerm : public ScalarTerm
{
- public :
+ BilinearTerm<T1,T2> &bilterm;
+ public :
+ BilinearTermToScalarTerm(BilinearTerm<T1,T2> &bilterm_): bilterm(bilterm_){}
+ virtual ~BilinearTermToScalarTerm() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+ {
+ fullMatrix<double> localMatrix;
+ bilterm.get(ele,npts,GP,localMatrix);
+ val=localMatrix(0,0);
+ }
+};
+
+
+class ScalarTermConstant : public ScalarTerm
+{
+ double val;
+ public :
+ ScalarTermConstant(double val_=1.0): val(val_) {}
+ virtual ~ScalarTermConstant() {}
virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
{
double jac[3][3];
@@ -76,16 +102,20 @@ class ScalarTermOne : public ScalarTerm
val+=weight*detJ;
}
}
+ virtual void get(MVertex *ver,double &val)
+ {
+ val=1;
+ }
};
-template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2>
+template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
{
public :
- LaplaceTerm(S1& space1_,S2& space2_) : BilinearTerm<S1,S2>(space1_,space2_)
+ LaplaceTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : BilinearTerm<T1,T2>(space1_,space2_)
{}
-
+ virtual ~LaplaceTerm() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
{
Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
@@ -98,15 +128,15 @@ template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2>
-template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
+template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
{
public :
- LaplaceTerm(S1& space1_) : BilinearTerm<S1,S1>(space1_,space1_)
+ LaplaceTerm(FunctionSpace<T1>& space1_) : BilinearTerm<T1,T1>(space1_,space1_)
{}
-
+ virtual ~LaplaceTerm() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
{
- int nbFF = BilinearTerm<S1,S1>::space1.getNumKeys(ele);
+ int nbFF = BilinearTerm<T1,T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF, nbFF);
m.setAll(0.);
@@ -114,8 +144,8 @@ template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
{
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 = ele->getJacobian(u, v, w, jac);
- std::vector<typename S1::GradType> Grads;
- BilinearTerm<S1,S1>::space1.gradf(ele,u, v, w, Grads);
+ std::vector<typename TensorialTraits<T1>::GradType> Grads;
+ BilinearTerm<T1,T1>::space1.gradf(ele,u, v, w, Grads);
for (int j = 0; j < nbFF; j++)
{
for (int k = j; k < nbFF; k++)
@@ -134,10 +164,11 @@ template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
-template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,S2> // non symmetric
+class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{
protected :
double E,nu;
+ bool sym;
fullMatrix<double> H;/* =
{ {C11, C12, C12, 0, 0, 0},
{C12, C11, C12, 0, 0, 0},
@@ -147,7 +178,7 @@ template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,
{ 0, 0, 0, 0, 0, C44} };*/
public :
- IsotropicElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6)
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),H(6,6)
{
double FACT = E / (1 + nu);
double C11 = FACT * (1 - nu) / (1 - 2 * nu);
@@ -156,69 +187,9 @@ template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,
H.scale(0.);
for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
+ sym=(&space1_==&space2_);
}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
- {
- int nbFF1 = BilinearTerm<S1,S2>::space1.getNumKeys(ele);
- int nbFF2 = BilinearTerm<S1,S2>::space2.getNumKeys(ele);
- double jac[3][3];
- fullMatrix<double> B(6, nbFF2);
- fullMatrix<double> BTH(nbFF2, 6);
- fullMatrix<double> BT(nbFF1, 6);
- m.resize(nbFF1, nbFF2);
- m.setAll(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 = ele->getJacobian(u, v, w, jac);
- std::vector<typename S1::GradType> Grads;// tableau de matrices...
- std::vector<typename S2::GradType> GradsT;// tableau de matrices...
- BilinearTerm<S1,S2>::space1.gradf(ele,u, v, w, Grads);
- BilinearTerm<S1,S2>::space2.gradf(ele,u, v, w, GradsT);
- for (int j = 0; j < nbFF1; j++)
- {
- BT(j, 0) = Grads[j](0,0);
- BT(j, 1) = Grads[j](1,1);
- BT(j, 2) = Grads[j](2,2);
- BT(j, 3) = Grads[j](0,1)+Grads[j](1,0);
- BT(j, 4) = Grads[j](1,2)+Grads[j](2,1);
- BT(j, 5) = Grads[j](0,2)+Grads[j](2,0);
- }
- for (int j = 0; j < nbFF2; j++)
- {
- B(0, j) = GradsT[j](0,0);
- B(1, j) = GradsT[j](1,1);
- B(2, j) = GradsT[j](2,2);
- B(3, j) = GradsT[j](0,1)+GradsT[j](1,0);
- B(4, j) = GradsT[j](1,2)+GradsT[j](2,1);
- B(5, j) = GradsT[j](0,2)+GradsT[j](2,0);
- }
- BTH.setAll(0.);
- BTH.gemm(BT, H);
- m.gemm(BTH, B, weight * detJ, 1.);
- }
- }
-}; // class
-
-
-
-
-
-
-template<class S1> class IsotropicElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
-{
- protected :
- double E,nu;
- fullMatrix<double> H;/* =
- { {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} };*/
-
- public :
- IsotropicElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6)
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space1_),E(E_),nu(nu_),H(6,6)
{
double FACT = E / (1 + nu);
double C11 = FACT * (1 - nu) / (1 - 2 * nu);
@@ -227,35 +198,80 @@ template<class S1> class IsotropicElasticTerm<S1,S1> : public BilinearTerm<S1,S1
H.scale(0.);
for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
+ sym=true;
}
-
+ virtual ~IsotropicElasticTerm() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
{
- int nbFF = BilinearTerm<S1,S1>::space1.getNumKeys(ele);
- double jac[3][3];
- fullMatrix<double> B(6, nbFF);
- fullMatrix<double> BTH(nbFF, 6);
- fullMatrix<double> BT(nbFF, 6);
- m.resize(nbFF, nbFF);
- m.setAll(0.);
- for (int i = 0; i < npts; i++)
+ if (sym)
{
- 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 = ele->getJacobian(u, v, w, jac);
- std::vector<typename S1::GradType> Grads;
- BilinearTerm<S1,S1>::space1.gradf(ele,u, v, w, Grads); // a optimiser ??
- for (int j = 0; j < nbFF; j++)
+ int nbFF = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+ double jac[3][3];
+ fullMatrix<double> B(6, nbFF);
+ fullMatrix<double> BTH(nbFF, 6);
+ fullMatrix<double> BT(nbFF, 6);
+ m.resize(nbFF, nbFF);
+ m.setAll(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 = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::GradType> Grads;
+ BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads); // a optimiser ??
+ for (int j = 0; j < nbFF; j++)
+ {
+ BT(j, 0) = B(0, j) = Grads[j](0,0);
+ BT(j, 1) = B(1, j) = Grads[j](1,1);
+ BT(j, 2) = B(2, j) = Grads[j](2,2);
+ BT(j, 3) = B(3, j) = Grads[j](0,1)+Grads[j](1,0);
+ BT(j, 4) = B(4, j) = Grads[j](1,2)+Grads[j](2,1);
+ BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
+ }
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+ m.gemm(BTH, B, weight * detJ, 1.);
+ }
+ }
+ else
+ {
+ int nbFF1 = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+ int nbFF2 = BilinearTerm<SVector3,SVector3>::space2.getNumKeys(ele);
+ double jac[3][3];
+ fullMatrix<double> B(6, nbFF2);
+ fullMatrix<double> BTH(nbFF2, 6);
+ fullMatrix<double> BT(nbFF1, 6);
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ for (int i = 0; i < npts; i++)
{
- BT(j, 0) = B(0, j) = Grads[j](0,0);
- BT(j, 1) = B(1, j) = Grads[j](1,1);
- BT(j, 2) = B(2, j) = Grads[j](2,2);
- BT(j, 3) = B(3, j) = Grads[j](0,1)+Grads[j](1,0);
- BT(j, 4) = B(4, j) = Grads[j](1,2)+Grads[j](2,1);
- BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
+ 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 = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::GradType> Grads;// tableau de matrices...
+ std::vector<TensorialTraits<SVector3>::GradType> GradsT;// tableau de matrices...
+ BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads);
+ BilinearTerm<SVector3,SVector3>::space2.gradf(ele,u, v, w, GradsT);
+ for (int j = 0; j < nbFF1; j++)
+ {
+ BT(j, 0) = Grads[j](0,0);
+ BT(j, 1) = Grads[j](1,1);
+ BT(j, 2) = Grads[j](2,2);
+ BT(j, 3) = Grads[j](0,1)+Grads[j](1,0);
+ BT(j, 4) = Grads[j](1,2)+Grads[j](2,1);
+ BT(j, 5) = Grads[j](0,2)+Grads[j](2,0);
+ }
+ for (int j = 0; j < nbFF2; j++)
+ {
+ B(0, j) = GradsT[j](0,0);
+ B(1, j) = GradsT[j](1,1);
+ B(2, j) = GradsT[j](2,2);
+ B(3, j) = GradsT[j](0,1)+GradsT[j](1,0);
+ B(4, j) = GradsT[j](1,2)+GradsT[j](2,1);
+ B(5, j) = GradsT[j](0,2)+GradsT[j](2,0);
+ }
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+ m.gemm(BTH, B, weight * detJ, 1.);
}
- BTH.setAll(0.);
- BTH.gemm(BT, H);
- m.gemm(BTH, B, weight * detJ, 1.);
}
}
}; // class
@@ -265,14 +281,16 @@ inline double dot(const double &a, const double &b)
{ return a*b; }
-template<class S1> class LoadTerm : public LinearTerm<S1>
+template<class T1> class LoadTerm : public LinearTerm<T1>
{
- simpleFunction<typename S1::ValType> &Load;
+ simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
public :
- LoadTerm(S1& space1_,simpleFunction<typename S1::ValType> &Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+ LoadTerm(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :LinearTerm<T1>(space1_),Load(Load_) {}
+ virtual ~LoadTerm() {}
+
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
{
- double nbFF=LinearTerm<S1>::space1.getNumKeys(ele);
+ double nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF);
m.scale(0.);
@@ -280,11 +298,11 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
{
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 = ele->getJacobian(u, v, w, jac);
- std::vector<typename S1::ValType> Vals;
- LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
+ std::vector<typename TensorialTraits<T1>::ValType> Vals;
+ LinearTerm<T1>::space1.f(ele,u, v, w, Vals);
SPoint3 p;
ele->pnt(u, v, w, p);
- typename S1::ValType load=Load(p.x(),p.y(),p.z());
+ typename TensorialTraits<T1>::ValType load=Load(p.x(),p.y(),p.z());
for (int j = 0; j < nbFF ; ++j)
{
m(j)+=dot(Vals[j],load)*weight*detJ;
@@ -294,12 +312,12 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
virtual void get(MVertex *ver,fullVector<double> &m)
{
- double nbFF=LinearTerm<S1>::space1.getNumKeys(ver);
+ double nbFF=LinearTerm<T1>::space1.getNumKeys(ver);
double jac[3][3];
m.resize(nbFF);
- std::vector<typename S1::ValType> Vals;
- LinearTerm<S1>::space1.f(ver, Vals);
- typename S1::ValType load=Load(ver->x(),ver->y(),ver->z());
+ std::vector<typename TensorialTraits<T1>::ValType> Vals;
+ LinearTerm<T1>::space1.f(ver, Vals);
+ typename TensorialTraits<T1>::ValType load=Load(ver->x(),ver->y(),ver->z());
for (int j = 0; j < nbFF ; ++j)
{
m(j)=dot(Vals[j],load);
--
GitLab