-added assembly for linear terms in dofmanager

-done BC with linear terms in elasiticiy

Solver/dofManager.h

@@ -193,25 +193,49 @@ class dofManager{
       }
     }
   }
-  inline void assemble(std::vector<Dof> &R, fullMatrix<dataMat> &m)
+
+  inline void assemble(std::vector<Dof> &R, fullVector<dataMat> &m) // for linear forms
   {
     if (!_current->isAllocated()) _current->allocate(unknown.size());
-
     std::vector<int> NR(R.size());
-
-    for (unsigned int i = 0; i < R.size(); i++){
+    for (unsigned int i = 0; i < R.size(); i++)
+    {
       std::map<Dof, int>::iterator itR = unknown.find(R[i]);
       if (itR != unknown.end()) NR[i] = itR->second;
       else NR[i] = -1;
     }
+    for (unsigned int i = 0; i < R.size(); i++)
+    {
+      if (NR[i] != -1)
+      {
+        _current->addToRightHandSide(NR[i], m(i));
+      }
+    }
+  }
 
-    for (unsigned int i = 0; i < R.size(); i++){
-      if (NR[i] != -1){
-	for (unsigned int j = 0; j < R.size(); j++){
-	  if (NR[j] != -1){
+
+  inline void assemble(std::vector<Dof> &R, fullMatrix<dataMat> &m)
+  {
+    if (!_current->isAllocated()) _current->allocate(unknown.size());
+    std::vector<int> NR(R.size());
+    for (unsigned int i = 0; i < R.size(); i++)
+    {
+      std::map<Dof, int>::iterator itR = unknown.find(R[i]);
+      if (itR != unknown.end()) NR[i] = itR->second;
+      else NR[i] = -1;
+    }
+    for (unsigned int i = 0; i < R.size(); i++)
+    {
+      if (NR[i] != -1)
+      {
+        for (unsigned int j = 0; j < R.size(); j++)
+        {
+          if (NR[j] != -1)
+          {
             _current->addToMatrix(NR[i], NR[j], m(i, j));
           }
-	  else{
+          else
+          {
             typename std::map<Dof,  dataVec>::iterator itFixed = fixed.find(R[j]);
             if (itFixed != fixed.end()){
               _current->addToRightHandSide(NR[i], -m(i, j) * itFixed->second);
@@ -221,6 +245,7 @@ class dofManager{
       }
     }
   }
+
   inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value)
   {
     assemble(Dof(entR, typeR), Dof(entC, typeC), value);

Solver/elasticitySolver.cpp

@@ -313,13 +313,15 @@ void elasticitySolver::solve()
       El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);      
     }
   }
-
+/*
   printf("-- done assembling!\n");
-  //  for (int i=0;i<pAssembler->sizeOfR();i++){
-    //    printf("%g ",lsys->getFromRightHandSide(i));
-  //  }
-  //  printf("\n");
-
+    for (int i=0;i<pAssembler->sizeOfR();i++)
+    {
+        if (fabs(lsys->getFromRightHandSide(i))>0.000001)
+        printf("%g ",lsys->getFromRightHandSide(i));
+    }
+    printf("\n");
+*/
   // solving
   lsys->systemSolve();
   printf("-- done solving!\n");
@@ -404,60 +406,109 @@ void MyelasticitySolver::solve()
     }
   }
 
+  VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
   //line forces
   for (std::map<int, SVector3 >::iterator it =lineForces.begin();
-       it != lineForces.end(); ++it){
-    DummyfemTerm El(pModel);
+       it != lineForces.end(); ++it)
+  {
     int iEdge = it->first;
     SVector3 f = it->second;
-    El.neumannNodalBC(iEdge, 1, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
-    El.neumannNodalBC(iEdge, 1, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
-    El.neumannNodalBC(iEdge, 1, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+    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)
+    {
+      GEntity *ge = itg->second[i];
+      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);
+      }
     }
+  }
+
 
 //face forces
   for (std::map<int, SVector3 >::iterator it = faceForces.begin();
-       it != faceForces.end(); ++it){
-    DummyfemTerm El(pModel);
+       it != faceForces.end(); ++it)
+  {
     int iFace = it->first;
     SVector3 f = it->second;
-    El.neumannNodalBC(iFace, 2, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
-    El.neumannNodalBC(iFace, 2, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
-    El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+    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)
+    {
+      GEntity *ge = itg->second[i];
+      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);
+      }
     }
+  }
+
 
 //volume forces
   for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
-       it != volumeForces.end(); ++it){
-    DummyfemTerm El(pModel);
+       it != volumeForces.end(); ++it)
+  {
     int iVolume = it->first;
     SVector3 f = it->second;
-//    El.setVector(f);
+    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
     printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
-    std::vector<GEntity*> ent = groups[_dim][iVolume];
-    for (unsigned int i = 0; i < ent.size(); i++){
-      //   to do
-      //      El.addToRightHandSide(*pAssembler, ent[i]);
-    }
+    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)
+    {
+      GEntity *ge = itg->second[i];
+      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);
       }
-
-/*
-  // assembling the stifness matrix
-  for (unsigned int i = 0; i < elasticFields.size(); i++){
-    SElement::setShapeEnrichement(elasticFields[i]._enrichment);
-    for (unsigned int j = 0; j < elasticFields.size(); j++){
-      elasticityTerm El(0, elasticFields[i]._E, elasticFields[i]._nu,
-                        elasticFields[i]._tag, elasticFields[j]._tag);
-      SElement::setTestEnrichement(elasticFields[j]._enrichment);
-      //      printf("%d %d\n",elasticFields[i]._tag,elasticFields[j]._tag);
-      El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);      
     }
   }
-*/
 
-  VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
 
   for (unsigned int i = 0; i < elasticFields.size(); i++)
   {
@@ -477,6 +528,31 @@ void MyelasticitySolver::solve()
       pAssembler->assemble(R, localMatrix);
     }
   }
+
+
+/*
+
+  for (std::map<int, SVector3 >::iterator it = volumeForces.begin();it != volumeForces.end(); ++it)
+  {
+    int iVolume = it->first;
+    SVector3 f = it->second;
+    printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
+    std::vector<GEntity*> ent = groups[_dim][iVolume];
+    for (unsigned int i = 0; i < ent.size(); i++)
+    {
+      //   to do
+      //      El.addToRightHandSide(*pAssembler, ent[i]);
+    }
+  }
+*/
+/*
+    for (int i=0;i<pAssembler->sizeOfR();i++)
+    {
+        if (fabs(lsys->getFromRightHandSide(i))>0.000001)
+        printf("%g ",lsys->getFromRightHandSide(i));
+    }
+    printf("\n");
+*/
   printf("-- done assembling!\n");
   lsys->systemSolve();
   printf("-- done solving!\n");

Solver/functionSpace.h

@@ -106,7 +106,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
   {
     int ndofs= ele->getNumVertices();
     int curpos=vals.size();
-    vals.reserve(curpos+ndofs);
+    vals.resize(curpos+ndofs);
     ele->getShapeFunctions(u, v, w, &(vals[curpos]));
   }; 
   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)

Solver/terms.h

@@ -82,8 +82,11 @@ class  LinearTermBase
 
 template<class S1> class LinearTerm : public LinearTermBase
 {
+ protected :
+  S1& space1;
  public : 
-  virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+  LinearTerm(S1& space1_) : space1(space1_) {}
+  virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) =0;
 };
 
 class  ScalarTermBase
@@ -168,6 +171,9 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no
 
 
 
+
+
+
 template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
 {
  protected : 
@@ -224,4 +230,34 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
 }; // class
 
 
+inline double dot(const double &a, const double &b)
+{ return a*b; }
+
+template<class S1> class LoadTerm : public LinearTerm<S1>
+{
+  typename S1::ValType Load;
+ public : 
+  LoadTerm(S1& space1_,typename S1::ValType Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+  virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
+  {
+    double nbFF=LinearTerm<S1>::space1.getNumKeys(ele);
+    double jac[3][3];
+    m.resize(nbFF);
+    m.scale(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::ValType> Vals;
+      LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
+      for (int j = 0; j < nbFF ; ++j)
+      {
+        m(j)+=dot(Vals[j],Load)*weight*detJ;
+      }
+    }
+  }
+};
+
+
+
 #endif// _TERMS_H_