diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 117005fdcfd439abaeaac523e5de690773dd17a9..1e25c25650253d07815f1de461393434ee8c5a64 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -15,6 +15,7 @@
#include "Numeric.h"
#include "functionSpace.h"
#include "terms.h"
+#include "solverAlgorithms.h"
#if defined(HAVE_POST)
#include "PView.h"
@@ -417,10 +418,6 @@ void MyelasticitySolver::solve()
LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
int dim=1;
- std::map<int, std::vector<GEntity*> > groups[4];
- pModel->getPhysicalGroups(groups);
- fullVector<double> localVector;
- std::vector<Dof> R;
std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iEdge);
if (itg == groups[dim].end()) {printf(" Nonexistent edge\n");break;}
for (unsigned int i = 0; i < itg->second.size(); ++i)
@@ -429,14 +426,7 @@ void MyelasticitySolver::solve()
for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
{
MElement *e = ge->getMeshElement(j);
- int integrationOrder = 2 * e->getPolynomialOrder();
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- R.clear();
- P123.getKeys(e,R);
- Lterm.get(e,npts,GP,localVector);
- pAssembler->assemble(R, localVector);
+ Assemble(Lterm,P123,e,*pAssembler);
}
}
}
@@ -451,10 +441,6 @@ void MyelasticitySolver::solve()
LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
int dim=2;
- std::map<int, std::vector<GEntity*> > groups[4];
- pModel->getPhysicalGroups(groups);
- fullVector<double> localVector;
- std::vector<Dof> R;
std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iFace);
if (itg == groups[dim].end()) {printf(" Nonexistent face\n");break;}
for (unsigned int i = 0; i < itg->second.size(); ++i)
@@ -463,14 +449,7 @@ void MyelasticitySolver::solve()
for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
{
MElement *e = ge->getMeshElement(j);
- int integrationOrder = 2 * e->getPolynomialOrder();
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- R.clear();
- P123.getKeys(e,R);
- Lterm.get(e,npts,GP,localVector);
- pAssembler->assemble(R, localVector);
+ Assemble(Lterm,P123,e,*pAssembler);
}
}
}
@@ -485,10 +464,6 @@ void MyelasticitySolver::solve()
LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
int dim=3;
- std::map<int, std::vector<GEntity*> > groups[4];
- pModel->getPhysicalGroups(groups);
- fullVector<double> localVector;
- std::vector<Dof> R;
std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iVolume);
if (itg == groups[dim].end()) {printf(" Nonexistent volume\n");break;}
for (unsigned int i = 0; i < itg->second.size(); ++i)
@@ -497,14 +472,7 @@ void MyelasticitySolver::solve()
for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
{
MElement *e = ge->getMeshElement(j);
- int integrationOrder = 2 * e->getPolynomialOrder();
- int npts;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- R.clear();
- P123.getKeys(e,R);
- Lterm.get(e,npts,GP,localVector);
- pAssembler->assemble(R, localVector);
+ Assemble(Lterm,P123,e,*pAssembler);
}
}
}
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index ca26866c02aa2ebf00d6078df84be8fab720b12c..105fea2127d7129472567354f7cadcad2e1e8457 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -61,8 +61,18 @@ template<> struct TensorialTraits<SVector3>
typedef SVector3 CurlType;
};
+
+
+class FunctionSpaceBase
+{
+ public:
+ virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
+ virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+};
+
template<class T>
-class FunctionSpace {
+class FunctionSpace : public FunctionSpaceBase
+{
protected:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
new file mode 100644
index 0000000000000000000000000000000000000000..2337d0afcb62de4402e816591bb6fdd4de371b8b
--- /dev/null
+++ b/Solver/solverAlgorithms.h
@@ -0,0 +1,85 @@
+//
+// C++ Interface: solverAlgorithms
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+#ifndef _SOLVERALGORITHMS_H_
+#define _SOLVERALGORITHMS_H_
+
+
+
+#include "terms.h"
+
+
+
+template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space1,Iterator itbegin,Iterator itend,Assembler &assembler) // symmetric
+{
+ fullMatrix<double> localMatrix;
+ std::vector<Dof> R;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MElement *e = *it;
+ R.clear();
+ int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ; // todo: to be given from outside
+ int npts=0;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ term.get(e,npts,GP,localMatrix);
+ space1.getKeys(e,R);
+ assembler.assemble(R, localMatrix);
+ }
+}
+
+template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space1,MElement *e,Assembler &assembler) // symmetric
+{
+ fullMatrix<double> localMatrix;
+ std::vector<Dof> R;
+ int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ; // todo: to be given from outside
+ int npts=0;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ term.get(e,npts,GP,localMatrix);
+ space1.getKeys(e,R);
+ assembler.assemble(R, localMatrix);
+}
+
+template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space1,Iterator itbegin,Iterator itend,Assembler &assembler)
+{
+ fullVector<double> localVector;
+ std::vector<Dof> R;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MElement *e = *it;
+ R.clear();
+ int integrationOrder = 2* (e->getPolynomialOrder()) ; // todo: to be given from outside
+ int npts=0;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ term.get(e,npts,GP,localVector);
+ space1.getKeys(e,R);
+ assembler.assemble(R, localVector);
+ }
+}
+
+template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space1,MElement *e,Assembler &assembler)
+{
+ fullVector<double> localVector;
+ std::vector<Dof> R;
+ int integrationOrder = 2* (e->getPolynomialOrder()) ; // todo: to be given from outside
+ int npts=0;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ term.get(e,npts,GP,localVector);
+ space1.getKeys(e,R);
+ assembler.assemble(R, localVector);
+}
+
+#endif// _SOLVERALGORITHMS_H_
\ No newline at end of file
diff --git a/Solver/terms.h b/Solver/terms.h
index 62161ba0612ad4f7972143b4411d9640396268b8..ab5ab4652fc9deb571a96948263fdebbee1cdaa0 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -61,7 +61,7 @@ class Formulation
class BilinearTermBase
{
- public :
+ public :
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
};
@@ -77,6 +77,7 @@ template<class S1,class S2> class BilinearTerm : public BilinearTermBase
class LinearTermBase
{
+ public:
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
};
@@ -259,22 +260,4 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
};
-template<class T,class S1, class I,class A> void Assemble(T& term,S1& space1,I itbegin,I itend,A& assembler) // symmetric
-{
- fullMatrix<double> localMatrix;
- std::vector<Dof> R;
- for (I it = itbegin;it!=itend; ++it)
- {
- MElement *e = *it;
- R.clear();
- int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
- int npts=0;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- term.get(e,npts,GP,localMatrix);
- space1.getKeys(e,R);
- assembler.assemble(R, localMatrix);
- }
-}
-
#endif// _TERMS_H_