You have to provide material parameters to instiantiate the class now

[[Imported from SVN: r4182]]

disc/elasticity/linearelasticityassembler.hh

@@ -24,7 +24,7 @@
 
 /**
  * @file
- * @brief  compute local stiffness matrix for conforming finite elements for linear elasticity equation
+ * @brief Compute local stiffness matrix for conforming finite elements for linear elasticity equation
  * @author Oliver Sander
  */
 
@@ -43,11 +43,11 @@ namespace Dune
 
   //! A class for computing local stiffness matrices
   /*! A class for computing local stiffness matrix for the
-        diffusion equation
+        linear elasticity equation
 
-            div j = q; j = -K grad u; in Omega
+        - div \sigma = f in Omega
 
-                u = g on Gamma1; j*n = J on Gamma2.
+        u = g on Gamma1; j*n = J on Gamma2.
 
         Uses conforming finite elements with the Lagrange shape functions.
         It should work for all dimensions and element types.
@@ -70,7 +70,6 @@ namespace Dune
 
     // some other sizes
     enum {dim=GridType::dimension};
-    //enum {SIZE=Dune::LagrangeShapeFunctionSetContainer<DT,RT,dim>::maxsize};
 
     //! The engineers' way of writing a symmetric second-order tensor
     typedef FieldVector<double, (dim+1)*dim/2> SymmTensor;
@@ -95,13 +94,9 @@ namespace Dune
     double nu_;
 
     //! Constructor
-    LinearElasticityLocalStiffness (bool procBoundaryAsDirichlet_=true)
-      : E_(2.5e5)
+    LinearElasticityLocalStiffness (double E, double nu, bool procBoundaryAsDirichlet_=true)
+      : E_(E), nu_(nu)
     {
-
-      //E_ = 2.5e5;
-      nu_ = 0.3;
-
       this->currentsize_ = 0;
       procBoundaryAsDirichlet = procBoundaryAsDirichlet_;
 
@@ -111,6 +106,13 @@ namespace Dune
         this->bctype[i] = BoundaryConditions::neumann;
     }
 
+    /** \brief Set material parameters */
+    void setEandNu(double E, double nu)
+    {
+      E_  = E;
+      nu_ = nu;
+    }
+
 
     //! assemble local stiffness matrix for given element and order
     /*! On exit the following things have been done:
@@ -131,10 +133,6 @@ namespace Dune
 
       this->currentsize_ = sfs.size();
 
-      //          std::cout << "Entity:" << std::endl;
-      //          for (int i=0; i<e.template count<n>(); i++)
-      //                std::cout << "corner i=" << i << " pos=" << e.geometry()[i] << std::endl;
-
       // clear assemble data
       for (int i=0; i<sfs.size(); i++)
       {
@@ -144,8 +142,6 @@ namespace Dune
           this->A[i][j] = 0;
       }
 
-      // Loop over all quadrature points and assemble matrix and right hand side
-
       // Compute suitable quadrature order
       int p=2;
       if (gt.isSimplex())
@@ -165,22 +161,12 @@ namespace Dune
         // eval jacobian inverse
         const Dune::FieldMatrix<DT,dim,dim> jac = e.geometry().jacobianInverseTransposed(local);
 
-#if 0
-        // eval diffusion tensor
-        const Dune::FieldMatrix<DT,dim,dim> K = problem.K(global,e,local);
-#endif
-
         // weight of quadrature point
         double weight = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].weight();
 
         // determinant of jacobian
         DT detjac = e.geometry().integrationElement(local);
 
-#if 0
-        // Source term;
-        RT q = problem.q(global,e,local);
-#endif
-
         RT factor = weight*detjac;
 
         // evaluate gradients at Gauss points
@@ -206,7 +192,6 @@ namespace Dune
             FieldMatrix<double, dim, dim> Grad(0);
             Grad[k] = grad[i];
 
-            //std::cout << Grad << std::endl;
             /* Computes the linear strain tensor from the deformation gradient*/
             computeStrain(Grad,strain[i*dim + k]);
 
@@ -225,28 +210,10 @@ namespace Dune
               for (int ccomp=0; ccomp<dim; ccomp++) {
 
                 this->A[row][col][rcomp][ccomp] += stress*strain[col*dim + ccomp] * weight * detjac;
-                //                               printf("adding %g to %d %d %d %d\n", stress*strain[col*dim + ccomp] * weight * detjac,
-                //                                      row, col, rcomp, ccomp);
 
               }
             }
           }
-#if 0
-          // rhs
-          b[i] += q*sfs[i].evaluateFunction(0,local)*factor;
-
-          // matrix
-          gv = 0;
-          K.umv(grad[i],gv);         // multiply with diffusion tensor
-          A[i][i] += (grad[i]*gv)*factor;
-          for (int j=0; j<i; j++)
-          {
-            RT t = (grad[j]*gv)*factor;
-            A[i][j] += t;
-            A[j][i] += t;
-          }
-
-#endif
 
         }