Skip to content
Snippets Groups Projects
Commit 00d385e8 authored by Éric Béchet's avatar Éric Béchet
Browse files

elastitysolver working w/ new "functionspaces"

parent 4469c334
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/elasticitySolver.cpp
+ 9
− 24
View file @ 00d385e8
...@@ -14,9 +14,7 @@ ...@@ -14,9 +14,7 @@
#include "linearSystemGMM.h" #include "linearSystemGMM.h"
#include "Numeric.h" #include "Numeric.h"
#include "functionSpace.h" #include "functionSpace.h"
#ifdef THIS_FILE_IS_MISSING
#include "terms.h" #include "terms.h"
#endif
#if defined(HAVE_POST) #if defined(HAVE_POST)
#include "PView.h" #include "PView.h"
...@@ -324,6 +322,7 @@ void elasticitySolver::solve() ...@@ -324,6 +322,7 @@ void elasticitySolver::solve()
// solving // solving
lsys->systemSolve(); lsys->systemSolve();
printf("-- done solving!\n");
} }
...@@ -444,6 +443,7 @@ void MyelasticitySolver::solve() ...@@ -444,6 +443,7 @@ void MyelasticitySolver::solve()
} }
} }
/*
// assembling the stifness matrix // assembling the stifness matrix
for (unsigned int i = 0; i < elasticFields.size(); i++){ for (unsigned int i = 0; i < elasticFields.size(); i++){
SElement::setShapeEnrichement(elasticFields[i]._enrichment); SElement::setShapeEnrichement(elasticFields[i]._enrichment);
...@@ -455,22 +455,15 @@ void MyelasticitySolver::solve() ...@@ -455,22 +455,15 @@ void MyelasticitySolver::solve()
El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g); El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);
} }
} }
*/
/* VectorLagrangeFunctionSpace L1(Dof::createTypeWithTwoInts(0,_tag),VectorLagrangeFunctionSpace::VECTOR_X); VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
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);
for (unsigned int i = 0; i < elasticFields.size(); i++) 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); ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
fullMatrix<double> localMatrix(12,12); fullMatrix<double> localMatrix;
std::vector<Dof> R;R.reserve(100); std::vector<Dof> R;
for ( groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end() ; ++it) for ( groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end() ; ++it)
{ {
MElement *e = *it; MElement *e = *it;
...@@ -480,19 +473,11 @@ throw; ...@@ -480,19 +473,11 @@ throw;
IntPt *GP; IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP); e->getIntegrationPoints(integrationOrder, &npts, &GP);
Eterm.get(e,npts,GP,localMatrix); Eterm.get(e,npts,GP,localMatrix);
El.addToMatrix(*pAssembler,localMatrix,R); P123.getKeys(e,R);
pAssembler->assemble(R, localMatrix);
} }
#endif
} }
printf("-- done assembling!\n"); printf("-- done assembling!\n");
// for (int i=0;i<pAssembler->sizeOfR();i++){
// printf("%g ",lsys->getFromRightHandSide(i));
// }
// printf("\n");
// solving
lsys->systemSolve(); lsys->systemSolve();
printf("-- done solving!\n");
} }
Solver/elasticityTerm.cpp
+ 10
− 17
View file @ 00d385e8
...@@ -59,19 +59,20 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const ...@@ -59,19 +59,20 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
if (se->getShapeEnrichement() == se->getTestEnrichement()){ if (se->getShapeEnrichement() == se->getTestEnrichement()){
for (int j = 0; j < nbNodes; j++){ 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, 0) = B(0, j) = Grads[j][0];
BT(j, 3) = B(3, j) = Grads[j][1]; 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, 1) = B(1, j + nbNodes) = Grads[j][1];
BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0]; 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, 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, 4) = B(4, j + 2 * nbNodes) = Grads[j][1];
BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1]; BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][0];
} }
} }
else{ else{
...@@ -80,15 +81,15 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const ...@@ -80,15 +81,15 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
for (int j = 0; j < nbNodes; j++){ for (int j = 0; j < nbNodes; j++){
BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0]; 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, 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, 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, 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, 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, 4) = Grads[j][1]; B(4, j + 2 * nbNodes) = GradsT[j][1];
BT(j + 2 * nbNodes, 5) = Grads[j][1]; B(5, 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 ...@@ -97,14 +98,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
m.gemm(BTH, B, weight * detJ, 1.); m.gemm(BTH, B, weight * detJ, 1.);
} }
return; 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 void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
......
Solver/femTerm.h
+ 2
− 27
View file @ 00d385e8
...@@ -101,22 +101,12 @@ class femTerm { ...@@ -101,22 +101,12 @@ class femTerm {
for (int k = 0; k < nbC; k++) for (int k = 0; k < nbC; k++)
C.push_back(getLocalDofC(se, 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) if (!sym)
dm.assemble(R, C, localMatrix); dm.assemble(R, C, localMatrix);
else else
dm.assemble(R, localMatrix); dm.assemble(R, localMatrix);
} }
void dirichletNodalBC(int physical, int dim, int comp, int field, void dirichletNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &e, const simpleFunction<dataVec> &e,
dofManager<dataVec> &dm) dofManager<dataVec> &dm)
...@@ -127,6 +117,7 @@ class femTerm { ...@@ -127,6 +117,7 @@ class femTerm {
for (unsigned int i = 0; i < v.size(); i++) 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())); 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, void neumannNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &fct, const simpleFunction<dataVec> &fct,
dofManager<dataVec> &dm) dofManager<dataVec> &dm)
...@@ -195,22 +186,6 @@ class DummyfemTerm : public femTerm<double> ...@@ -195,22 +186,6 @@ class DummyfemTerm : public femTerm<double>
virtual Dof getLocalDofC(SElement *se, int iCol) const {return Dof(0,0);} 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 elementMatrix(SElement *se, fullMatrix<dataMat> &m) const {m.scale(0.);}
virtual void elementVector(SElement *se, fullVector<dataVec> &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);
}
}; };
......
Solver/functionSpace.h
+ 23
− 88
View file @ 00d385e8
...@@ -5,11 +5,11 @@ ...@@ -5,11 +5,11 @@
#include "STensor3.h" #include "STensor3.h"
#include <vector> #include <vector>
#include <iterator> #include <iterator>
#include <iostream>
#include "Numeric.h" #include "Numeric.h"
#include "MElement.h" #include "MElement.h"
#include "dofManager.h" #include "dofManager.h"
//class STensor3{};
class SVoid{}; class SVoid{};
class basisFunction{ class basisFunction{
...@@ -79,7 +79,6 @@ class FunctionSpace { ...@@ -79,7 +79,6 @@ class FunctionSpace {
virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs); 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 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 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; virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
}; };
...@@ -97,8 +96,6 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> ...@@ -97,8 +96,6 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
Dof getLocalDof(MElement *ele, int i) const Dof getLocalDof(MElement *ele, int i) const
{ {
// int iComp = i / ele->getNumVertices();
// int ithLocalVertex = i % ele->getNumVertices();
return Dof(ele->getVertex(i)->getNum(), _iField); return Dof(ele->getVertex(i)->getNum(), _iField);
} }
...@@ -109,18 +106,18 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> ...@@ -109,18 +106,18 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
{ {
int ndofs= ele->getNumVertices(); int ndofs= ele->getNumVertices();
int curpos=vals.size(); int curpos=vals.size();
vals.resize(curpos+ndofs); vals.reserve(curpos+ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos])); ele->getShapeFunctions(u, v, w, &(vals[curpos]));
}; };
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{ {
int ndofs= ele->getNumVertices(); int ndofs= ele->getNumVertices();
// grads.reserve(grads.size()+ndofs); grads.reserve(grads.size()+ndofs);
double gradsuvw[256][3]; double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw); ele->getGradShapeFunctions(u, v, w, gradsuvw);
double jac[3][3]; double jac[3][3];
double invjac[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); inv3x3(jac, invjac);
for (int i=0;i<ndofs;++i) for (int i=0;i<ndofs;++i)
grads.push_back(GradType( grads.push_back(GradType(
...@@ -130,61 +127,16 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> ...@@ -130,61 +127,16 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
)); ));
}; };
virtual int getNumKeys(MElement *ele) {return ele->getNumVertices();} 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 ... virtual int getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{ {
int ndofs= ele->getNumVertices(); int ndofs= ele->getNumVertices();
// keys.reserve(keys.size()+ndofs); keys.reserve(keys.size()+ndofs);
for (int i=0;i<ndofs;++i) for (int i=0;i<ndofs;++i)
keys.push_back(getLocalDof(ele,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> template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T>
{ {
...@@ -224,6 +176,8 @@ public : ...@@ -224,6 +176,8 @@ public :
ScalarFS->f(ele,u,v,w,valsd); ScalarFS->f(ele,u,v,w,valsd);
int nbdofs=valsd.size(); int nbdofs=valsd.size();
int nbcomp=comp.size(); int nbcomp=comp.size();
int curpos=vals.size();
vals.reserve(curpos+nbcomp*nbdofs);
for (int j=0;j<nbcomp;++j) for (int j=0;j<nbcomp;++j)
{ {
for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]); for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
...@@ -236,6 +190,8 @@ public : ...@@ -236,6 +190,8 @@ public :
ScalarFS->gradf(ele,u,v,w,gradsd); ScalarFS->gradf(ele,u,v,w,gradsd);
int nbdofs=gradsd.size(); int nbdofs=gradsd.size();
int nbcomp=comp.size(); int nbcomp=comp.size();
int curpos=grads.size();
grads.reserve(curpos+nbcomp*nbdofs);
GradType val; GradType val;
for (int j=0;j<nbcomp;++j) for (int j=0;j<nbcomp;++j)
{ {
...@@ -245,47 +201,28 @@ public : ...@@ -245,47 +201,28 @@ public :
grads.push_back(val); 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) virtual int getKeys(MElement *ele, std::vector<Dof> &keys)
{ {
int nk=ScalarFS->getNumKeys(ele); int nk=ScalarFS->getNumKeys(ele);
std::vector<Dof> bufk; std::vector<Dof> bufk;
bufk.reserve(nk); bufk.reserve(nk);
ScalarFS->getKeys(ele,bufk); ScalarFS->getKeys(ele,bufk);
int nbdofs=bufk.size();
int nbcomp=comp.size(); int nbcomp=comp.size();
int curpos=keys.size();
keys.reserve(curpos+nbcomp*nbdofs);
for (int j=0;j<nbcomp;++j) for (int j=0;j<nbcomp;++j)
{ {
for (int i=0;i<nk;++i) for (int i=0;i<nk;++i)
{ {
int i1,i2; int i1,i2;
Dof::getTwoIntsFromType(bufk[i].getType(), 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> ...@@ -316,7 +253,7 @@ class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
{} {}
VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2, Along comp3) : VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2, Along comp3) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id), 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> ...@@ -406,10 +343,9 @@ class xFemFunctionSpace : public FunctionSpace<T>
FunctionSpace<T>* _spacebase; FunctionSpace<T>* _spacebase;
// Function<double>* enrichment; // Function<double>* enrichment;
public: public:
ValType f(MElement *ele, double u, double v, double w); virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
GradType gradf(MElement *ele, double u, double v, double w); virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
int getNumKeys(MElement *ele); int getNumKeys(MElement *ele);
void getKeys(MElement *ele, Dof *keys);
void getKeys(MElement *ele, std::vector<Dof> &keys); void getKeys(MElement *ele, std::vector<Dof> &keys);
}; };
...@@ -427,10 +363,9 @@ class FilteredFunctionSpace : public FunctionSpace<T> ...@@ -427,10 +363,9 @@ class FilteredFunctionSpace : public FunctionSpace<T>
F &_filter; F &_filter;
public: public:
ValType f(MElement *ele, double u, double v, double w); virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
GradType gradf(MElement *ele, double u, double v, double w); virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
int getNumKeys(MElement *ele); int getNumKeys(MElement *ele);
void getKeys(MElement *ele, Dof *keys);
void getKeys(MElement *ele, std::vector<Dof> &keys); void getKeys(MElement *ele, std::vector<Dof> &keys);
}; };
......
Solver/terms.h 0 → 100644
+ 227
− 0
View file @ 00d385e8
//
// 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_
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment