diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index eec404cde3ac0a82aebb9fbd2260110db5fef0e4..35a5c62715a1ebd8e9f0a317e93e4ceb70eda321 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -14,9 +14,7 @@
#include "linearSystemGMM.h"
#include "Numeric.h"
#include "functionSpace.h"
-#ifdef THIS_FILE_IS_MISSING
#include "terms.h"
-#endif
#if defined(HAVE_POST)
#include "PView.h"
@@ -324,6 +322,7 @@ void elasticitySolver::solve()
// solving
lsys->systemSolve();
+ printf("-- done solving!\n");
}
@@ -444,6 +443,7 @@ void MyelasticitySolver::solve()
}
}
+/*
// assembling the stifness matrix
for (unsigned int i = 0; i < elasticFields.size(); i++){
SElement::setShapeEnrichement(elasticFields[i]._enrichment);
@@ -455,22 +455,15 @@ void MyelasticitySolver::solve()
El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);
}
}
+*/
-/* VectorLagrangeFunctionSpace L1(Dof::createTypeWithTwoInts(0,_tag),VectorLagrangeFunctionSpace::VECTOR_X);
- VectorLagrangeFunctionSpace L2(Dof::createTypeWithTwoInts(1,_tag),VectorLagrangeFunctionSpace::VECTOR_Y);
- VectorLagrangeFunctionSpace L3(Dof::createTypeWithTwoInts(2,_tag),VectorLagrangeFunctionSpace::VECTOR_Z);
- CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType> P123(L1,L2,L3);*/
- VectorLagrangeFunctionSpace P123(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+ VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
- DummyfemTerm El(pModel);
-// ElasticTerm<CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType>,CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType> > Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
-throw;
-#ifdef DEBUG_ME_ELASTIC_TERM_IS_NOT_DEFINE_IT_IS_PROBABLY_IN_TERM_H_BUT_THIS_FILE_IS_NOT_IN_THE_SVN
ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
- fullMatrix<double> localMatrix(12,12);
- std::vector<Dof> R;R.reserve(100);
+ fullMatrix<double> localMatrix;
+ std::vector<Dof> R;
for ( groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end() ; ++it)
{
MElement *e = *it;
@@ -480,19 +473,11 @@ throw;
IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP);
Eterm.get(e,npts,GP,localMatrix);
- El.addToMatrix(*pAssembler,localMatrix,R);
+ P123.getKeys(e,R);
+ pAssembler->assemble(R, localMatrix);
}
-#endif
}
-
printf("-- done assembling!\n");
- // for (int i=0;i<pAssembler->sizeOfR();i++){
- // printf("%g ",lsys->getFromRightHandSide(i));
- // }
- // printf("\n");
-
- // solving
lsys->systemSolve();
+ printf("-- done solving!\n");
}
-
-
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 7fd8bbec30707dd574b744fcd6b3c47c64137f2a..e60656d05707d30dcebe3e5694263a36322e68ac 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -59,19 +59,20 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
if (se->getShapeEnrichement() == se->getTestEnrichement()){
for (int j = 0; j < nbNodes; j++){
- //printf(" GR(j) = %12.5E,%12.5E,%12.5E\n", Grads[j][0],Grads[j][1],Grads[j][2]);
+// printf(" GR(j) = %12.5E,%12.5E,%12.5E\n", Grads[j][0],Grads[j][1],Grads[j][2]);
BT(j, 0) = B(0, j) = Grads[j][0];
BT(j, 3) = B(3, j) = Grads[j][1];
- BT(j, 4) = B(4, j) = Grads[j][2];
+ BT(j, 5) = B(5, j) = Grads[j][2];
BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
- BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
+ BT(j + nbNodes, 4) = B(4, j + nbNodes) = Grads[j][2];
BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
- BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
- BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
+ BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][1];
+ BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][0];
+
}
}
else{
@@ -80,15 +81,15 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
for (int j = 0; j < nbNodes; j++){
BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1];
- BT(j, 4) = Grads[j][2]; B(4, j) = GradsT[j][2];
+ BT(j, 5) = Grads[j][2]; B(5, j) = GradsT[j][2];
BT(j + nbNodes, 1) = Grads[j][1]; B(1, j + nbNodes) = GradsT[j][1];
BT(j + nbNodes, 3) = Grads[j][0]; B(3, j + nbNodes) = GradsT[j][0];
- BT(j + nbNodes, 5) = Grads[j][2]; B(5, j + nbNodes) = GradsT[j][2];
+ BT(j + nbNodes, 4) = Grads[j][2]; B(4, j + nbNodes) = GradsT[j][2];
BT(j + 2 * nbNodes, 2) = Grads[j][2]; B(2, j + 2 * nbNodes) = GradsT[j][2];
- BT(j + 2 * nbNodes, 4) = Grads[j][0]; B(4, j + 2 * nbNodes) = GradsT[j][0];
- BT(j + 2 * nbNodes, 5) = Grads[j][1]; B(5, j + 2 * nbNodes) = GradsT[j][1];
+ BT(j + 2 * nbNodes, 4) = Grads[j][1]; B(4, j + 2 * nbNodes) = GradsT[j][1];
+ BT(j + 2 * nbNodes, 5) = Grads[j][0]; B(5, j + 2 * nbNodes) = GradsT[j][0];
}
}
@@ -97,14 +98,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
m.gemm(BTH, B, weight * detJ, 1.);
}
return;
- for (int i = 0; i < 3 * nbNodes; i++){
- for (int j = 0; j < 3 * nbNodes; j++){
- printf("%g ", m(i, j));
- }
- printf("\n");
- }
- printf("\n");
- printf("\n");
}
void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index ee759aa75ca90841282cabba8f773b8baff34700..5a16cac81c8059091f52f9eac1c9d860a14d610c 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -101,22 +101,12 @@ class femTerm {
for (int k = 0; k < nbC; k++)
C.push_back(getLocalDofC(se, k));
}
-/* for (int i = 0; i < nbC; i++) */
-/* for (int j = 0; j < nbR; j++) */
-/* dm.assemble(getLocalDofR(se, i), getLocalDofC(se, j), localMatrix(i,j)); */
-/* return; */
-
-
-// if (nbR == nbC){
-// for (int i=0;i<nbR;i++)
-// if (!(C[i] == R[i]))sym = false;
-// }
-// else sym = false;
if (!sym)
dm.assemble(R, C, localMatrix);
else
dm.assemble(R, localMatrix);
}
+
void dirichletNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &e,
dofManager<dataVec> &dm)
@@ -127,6 +117,7 @@ class femTerm {
for (unsigned int i = 0; i < v.size(); i++)
dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
}
+
void neumannNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &fct,
dofManager<dataVec> &dm)
@@ -195,22 +186,6 @@ class DummyfemTerm : public femTerm<double>
virtual Dof getLocalDofC(SElement *se, int iCol) const {return Dof(0,0);}
virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const {m.scale(0.);}
virtual void elementVector(SElement *se, fullVector<dataVec> &m) const {m.scale(0.);}
- public:
- void addToMatrix(dofManager<dataVec> &dm,
- fullMatrix<dataMat> &localMatrix,
- std::vector<Dof> R,std::vector<Dof> C) const
- {
- dm.assemble(R, C, localMatrix);
- }
-
-
- void addToMatrix(dofManager<dataVec> &dm,
- fullMatrix<dataMat> &localMatrix,
- std::vector<Dof> R) const // symmetric version.
- {
- dm.assemble(R, localMatrix);
- }
-
};
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 3a540f69ec38f24062ee27220896ca1aac28585b..0130c38804694d1f607680646ec7e4302d358556 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -5,11 +5,11 @@
#include "STensor3.h"
#include <vector>
#include <iterator>
+#include <iostream>
#include "Numeric.h"
#include "MElement.h"
#include "dofManager.h"
-//class STensor3{};
class SVoid{};
class basisFunction{
@@ -79,7 +79,6 @@ class FunctionSpace {
virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs);
virtual int curlf(MElement *ele, double u, double v, double w,std::vector<CurlType> &curls);*/
virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
- virtual int getKeys(MElement *ele, Dof *keys)=0; // may be faster once the number of dofs is known
virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
};
@@ -97,8 +96,6 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
Dof getLocalDof(MElement *ele, int i) const
{
-// int iComp = i / ele->getNumVertices();
-// int ithLocalVertex = i % ele->getNumVertices();
return Dof(ele->getVertex(i)->getNum(), _iField);
}
@@ -109,18 +106,18 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
{
int ndofs= ele->getNumVertices();
int curpos=vals.size();
- vals.resize(curpos+ndofs);
+ vals.reserve(curpos+ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
};
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
int ndofs= ele->getNumVertices();
-// grads.reserve(grads.size()+ndofs);
+ grads.reserve(grads.size()+ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
double jac[3][3];
double invjac[3][3];
- const double detJ = ele->getJacobian(u, v, w, jac); // a faire une fois pour toute ??
+ const double detJ = ele->getJacobian(u, v, w, jac); // redondant : on fait cet appel a l'exterieur
inv3x3(jac, invjac);
for (int i=0;i<ndofs;++i)
grads.push_back(GradType(
@@ -130,61 +127,16 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
));
};
virtual int getNumKeys(MElement *ele) {return ele->getNumVertices();}
- virtual int getKeys(MElement *ele, Dof *keys)
- {
- int ndofs= ele->getNumVertices();
- for (int i=0;i<ndofs;++i)
- keys[i]=getLocalDof(ele,i);
- }
+
virtual int getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{
int ndofs= ele->getNumVertices();
-// keys.reserve(keys.size()+ndofs);
+ keys.reserve(keys.size()+ndofs);
for (int i=0;i<ndofs;++i)
keys.push_back(getLocalDof(ele,i));
}
};
-/*
-template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T> // scalarFS * const vector (avec vecteur non const, peut etre utilise pour xfem directement
-{
-public :
- typedef typename TensorialTraits<T>::ValType ValType;
- typedef typename TensorialTraits<T>::GradType GradType;
- typedef typename TensorialTraits<T>::HessType HessType;
- typedef typename TensorialTraits<T>::DivType DivType;
- typedef typename TensorialTraits<T>::CurlType CurlType;
-protected :
- T multiplier;
- FunctionSpace<double> *ScalarFS;
-public :
- template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult): multiplier(mult),ScalarFS(new T2(SFS)) {}
- virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
- virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
- {
- std::vector<double> valsd;
- ScalarFS->f(ele,u,v,w,valsd);
- int nbdofs=valsd.size();
- for (int i=0;i<nbdofs;++i) vals.push_back(multiplier*valsd[i]);
- }
-
- virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
- {
- std::vector<SVector3> gradsd;
- ScalarFS->gradf(ele,u,v,w,gradsd);
- int nbdofs=gradsd.size();
- GradType val;
- for (int i=0;i<nbdofs;++i)
- {
- tensprod(multiplier,gradsd[i],val);
- grads.push_back(val);
- }
- };
- virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele);}
- virtual int getKeys(MElement *ele, Dof *keys){ ScalarFS->getKeys(ele,keys);}
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys) {ScalarFS->getKeys(ele,keys);}
-};
-*/
template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T>
{
@@ -224,6 +176,8 @@ public :
ScalarFS->f(ele,u,v,w,valsd);
int nbdofs=valsd.size();
int nbcomp=comp.size();
+ int curpos=vals.size();
+ vals.reserve(curpos+nbcomp*nbdofs);
for (int j=0;j<nbcomp;++j)
{
for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
@@ -236,6 +190,8 @@ public :
ScalarFS->gradf(ele,u,v,w,gradsd);
int nbdofs=gradsd.size();
int nbcomp=comp.size();
+ int curpos=grads.size();
+ grads.reserve(curpos+nbcomp*nbdofs);
GradType val;
for (int j=0;j<nbcomp;++j)
{
@@ -245,47 +201,28 @@ public :
grads.push_back(val);
}
}
- };
- virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele)*comp.size();}
- virtual int getKeys(MElement *ele, Dof *keys)
- {
- int nbcomp=comp.size();
- int nk=ScalarFS->getNumKeys(ele);
- Dof *kptr=keys;
- ScalarFS->getKeys(ele,kptr);
- kptr+=nk;
- for (int j=1;j<nbcomp;++j)
- {
- for (int i=0;i<nk;++i)
- kptr[i]=kptr[i-nk];
- kptr+=nk;
- }
- kptr=keys;
- for (int j=0;j<nbcomp;++j)
- {
- for (int i=0;i<nk;++i)
- {
- int i1,i2;
- Dof::getTwoIntsFromType(kptr[i].getType(), i1,i2);
- kptr[i]=Dof(kptr[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i2));
- }
- kptr+=nk;
- }
}
+
+ virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele)*comp.size();}
+
virtual int getKeys(MElement *ele, std::vector<Dof> &keys)
{
+
int nk=ScalarFS->getNumKeys(ele);
std::vector<Dof> bufk;
bufk.reserve(nk);
ScalarFS->getKeys(ele,bufk);
+ int nbdofs=bufk.size();
int nbcomp=comp.size();
+ int curpos=keys.size();
+ keys.reserve(curpos+nbcomp*nbdofs);
for (int j=0;j<nbcomp;++j)
{
for (int i=0;i<nk;++i)
{
int i1,i2;
Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
- keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i2)));
+ keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i1)));
}
}
}
@@ -316,7 +253,7 @@ class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
{}
VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2, Along comp3) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
- BasisVectors[comp1],comp3, BasisVectors[comp2],comp2, BasisVectors[comp3],comp3)
+ BasisVectors[comp1],comp1, BasisVectors[comp2],comp2, BasisVectors[comp3],comp3)
{}
};
@@ -406,10 +343,9 @@ class xFemFunctionSpace : public FunctionSpace<T>
FunctionSpace<T>* _spacebase;
// Function<double>* enrichment;
public:
- ValType f(MElement *ele, double u, double v, double w);
- GradType gradf(MElement *ele, double u, double v, double w);
+ virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
+ virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
int getNumKeys(MElement *ele);
- void getKeys(MElement *ele, Dof *keys);
void getKeys(MElement *ele, std::vector<Dof> &keys);
};
@@ -427,10 +363,9 @@ class FilteredFunctionSpace : public FunctionSpace<T>
F &_filter;
public:
- ValType f(MElement *ele, double u, double v, double w);
- GradType gradf(MElement *ele, double u, double v, double w);
+ virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
+ virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
int getNumKeys(MElement *ele);
- void getKeys(MElement *ele, Dof *keys);
void getKeys(MElement *ele, std::vector<Dof> &keys);
};
diff --git a/Solver/terms.h b/Solver/terms.h
new file mode 100644
index 0000000000000000000000000000000000000000..df76f2987d57c76aa361ba5060c9b01b056b02d4
--- /dev/null
+++ b/Solver/terms.h
@@ -0,0 +1,227 @@
+//
+// C++ Interface: terms
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+#ifndef _TERMS_H_
+#define _TERMS_H_
+
+#include "SVector3.h"
+#include <vector>
+#include <iterator>
+#include "Numeric.h"
+#include "functionSpace.h"
+
+
+// evaluation of a field ???
+template<class T>
+class Field {
+ protected:
+ typedef typename TensorialTraits<T>::ValType ValType;
+ typedef typename TensorialTraits<T>::GradType GradType;
+ dofManager<double> *dm; //
+/* typedef typename TensorialTraits<T>::HessType HessType;
+ typedef typename TensorialTraits<T>::DivType DivType;
+ typedef typename TensorialTraits<T>::CurlType CurlType;*/
+
+ public:
+ virtual int f(MElement *ele, double u, double v, double w, ValType &val)=0;
+ virtual int gradf(MElement *ele, double u, double v, double w,GradType &grad)=0;
+// virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads, STensor3 &invjac)=0;// on passe le jacobien que l'on veut ...
+/* virtual int hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hesss);
+ virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs);
+ virtual int curlf(MElement *ele, double u, double v, double w,std::vector<CurlType> &curls);*/
+ virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
+ virtual int getKeys(MElement *ele, Dof *keys)=0; // may be faster once the number of dofs is known
+ virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+};
+
+
+
+
+
+
+class Formulation
+{
+ std::vector<FunctionSpace<double>* > scalarfs;
+ std::vector<FunctionSpace<SVector3>* > vectorfs;
+ std::vector<groupOfElements* > groups;
+ std::vector<std::pair<MElement*,std::vector<groupOfElements&> > > links;
+ dofManager<double> *dm; //
+
+};
+
+class BilinearTermBase
+{
+ public :
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+};
+
+template<class S1,class S2> class BilinearTerm : public BilinearTermBase
+{
+ protected :
+ S1& space1;
+ S2& space2;
+ public :
+ BilinearTerm(S1& space1_,S2& space2_) : space1(space1_),space2(space2_) {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+};
+
+class LinearTermBase
+{
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
+};
+
+template<class S1> class LinearTerm : public LinearTermBase
+{
+ public :
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+};
+
+class ScalarTermBase
+{
+ public :
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &v) =0;
+};
+
+class ScalarTerm : public ScalarTermBase
+{
+ public :
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &v) =0;
+};
+
+
+
+template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // non 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 :
+ ElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6)
+ {
+ double FACT = E / (1 + nu);
+ double C11 = FACT * (1 - nu) / (1 - 2 * nu);
+ double C12 = FACT * nu / (1 - 2 * nu);
+ double C44 = (C11 - C12) / 2;
+ 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;
+ }
+ 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 ElasticTerm<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 :
+ ElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6)
+ {
+ double FACT = E / (1 + nu);
+ double C11 = FACT * (1 - nu) / (1 - 2 * nu);
+ double C12 = FACT * nu / (1 - 2 * nu);
+ double C44 = (C11 - C12) / 2;
+ 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;
+ }
+
+ 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++)
+ {
+ 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++)
+ {
+ 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.);
+ }
+ }
+}; // class
+
+
+#endif// _TERMS_H_