dedicated header for solver-related algorithms

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);
       }
     }
   }

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;

Solver/solverAlgorithms.h

+//
+// 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

Solver/terms.h

@@ -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_