From eb165f8cf8539903f46e44a1a9094d6d665367c6 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@dune-project.org>
Date: Thu, 23 Feb 2006 11:14:34 +0000
Subject: [PATCH] fixes and cleanup -- there'll be more
[[Imported from SVN: r4165]]
---
disc/elasticity/linearelasticityassembler.hh | 114 +++++--------------
1 file changed, 26 insertions(+), 88 deletions(-)
diff --git a/disc/elasticity/linearelasticityassembler.hh b/disc/elasticity/linearelasticityassembler.hh
index e17df3ef7..c420b16b5 100644
--- a/disc/elasticity/linearelasticityassembler.hh
+++ b/disc/elasticity/linearelasticityassembler.hh
@@ -17,14 +17,10 @@
#include "grid/common/referenceelements.hh"
#include "istl/operators.hh"
#include "istl/bvector.hh"
-#include "istl/bcrsmatrix.hh"
+#include <dune/disc/operators/localstiffness.hh>
#include "disc/shapefunctions/lagrangeshapefunctions.hh"
-#include "disc/operators/p1operator.hh"
#include "disc/operators/boundaryconditions.hh"
-//#include"disc/functions/p0function.hh"
-#include "disc/functions/p1function.hh"
-//#include"groundwater.hh"
/**
* @file
@@ -40,7 +36,7 @@ namespace Dune
* @{
*/
/**
- * @brief compute local stiffness matrix for conforming finite elements for diffusion equation
+ * @brief compute local stiffness matrix for the linear elasticity operator
*
*/
@@ -60,11 +56,11 @@ namespace Dune
Template parameters are:
- - Grid a DUNE grid type
+ - GridType a DUNE grid type
- RT type used for return values
*/
template<class GridType, class RT>
- class LinearElasticityLocalStiffness
+ class LinearElasticityLocalStiffness : public LocalStiffness<GridType,RT,GridType::dimension>
{
// grid types
typedef typename GridType::ctype DT;
@@ -74,7 +70,7 @@ namespace Dune
// some other sizes
enum {dim=GridType::dimension};
- enum {SIZE=Dune::LagrangeShapeFunctionSetContainer<DT,RT,dim>::maxsize};
+ //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;
@@ -99,21 +95,20 @@ namespace Dune
double nu_;
//! Constructor
- LinearElasticityLocalStiffness ( /*const GroundwaterEquationParameters<G,RT>& params, */
- bool procBoundaryAsDirichlet_=true)
- //: problem(params)
+ LinearElasticityLocalStiffness (bool procBoundaryAsDirichlet_=true)
+ : E_(2.5e5)
{
- E_ = 2.5e5;
+ //E_ = 2.5e5;
nu_ = 0.3;
- currentsize = 0;
+ this->currentsize_ = 0;
procBoundaryAsDirichlet = procBoundaryAsDirichlet_;
// For the time being: all boundary conditions are homogeneous Neumann
// This means no boundary condition handling is done at all
- for (int i=0; i<SIZE; i++)
- bctype[i] = BoundaryConditions::neumann;
+ for (int i=0; i<LocalStiffness<GridType,RT,m>::SIZE; i++)
+ this->bctype[i] = BoundaryConditions::neumann;
}
@@ -134,7 +129,7 @@ namespace Dune
const typename LagrangeShapeFunctionSetContainer<DT,RT,dim>::value_type& sfs
= LagrangeShapeFunctions<DT,RT,dim>::general(gt,k);
- currentsize = sfs.size();
+ this->currentsize_ = sfs.size();
// std::cout << "Entity:" << std::endl;
// for (int i=0; i<e.template count<n>(); i++)
@@ -143,10 +138,10 @@ namespace Dune
// clear assemble data
for (int i=0; i<sfs.size(); i++)
{
- b[i] = 0;
- bctype[i][0] = BoundaryConditions::neumann;
+ this->b[i] = 0;
+ this->bctype[i][0] = BoundaryConditions::neumann;
for (int j=0; j<sfs.size(); j++)
- A[i][j] = 0;
+ this->A[i][j] = 0;
}
// Loop over all quadrature points and assemble matrix and right hand side
@@ -188,7 +183,7 @@ namespace Dune
RT factor = weight*detjac;
// evaluate gradients at Gauss points
- Dune::FieldVector<DT,dim> grad[SIZE], temp, gv;
+ Dune::FieldVector<DT,dim> grad[LocalStiffness<GridType,RT,m>::SIZE], temp, gv;
for (int i=0; i<sfs.size(); i++) {
for (int l=0; l<dim; l++)
@@ -228,7 +223,7 @@ namespace Dune
for (int ccomp=0; ccomp<dim; ccomp++) {
- A[row][col][rcomp][ccomp] += stress*strain[col*dim + ccomp] * weight * detjac;
+ this->A[row][col][rcomp][ccomp] += stress*strain[col*dim + ccomp] * weight * detjac;
}
}
@@ -254,6 +249,7 @@ namespace Dune
}
+#if 0
for (int row=0; row<sfs.size(); row++)
for (int col=0; col<=row; col++) {
@@ -262,40 +258,12 @@ namespace Dune
if (row!=col) {
for (int rcomp=0; rcomp<dim; rcomp++)
for (int ccomp=0; ccomp<dim; ccomp++)
- A[col][row][ccomp][rcomp] = A[row][col][rcomp][ccomp];
- }
- }
-
-#if 0
- static int eleme = 0;
- printf("********* Element %d **********\n", eleme++);
- for (int row=0; row<sfs.size(); row++) {
-
- for (int rcomp=0; rcomp<dim; rcomp++) {
-
- for (int col=0; col<sfs.size(); col++) {
-
- for (int ccomp=0; ccomp<dim; ccomp++) {
-
- std::cout << A[row][col][rcomp][ccomp] << " ";
-
- }
-
- std::cout << " ";
-
+ this->A[col][row][ccomp][rcomp] = this->A[row][col][rcomp][ccomp];
}
- std::cout << std::endl;
-
-
-
}
-
- std::cout << std::endl;
-
- }
- //exit(0);
#endif
+ std::cout << this->A[0][0] << std::endl;
#if 0
// evaluate boundary conditions via intersection iterator
IntersectionIterator endit = e.iend();
@@ -484,41 +452,7 @@ namespace Dune
}
- // print contents of local stiffness matrix
- void print (std::ostream& s, int width, int precision)
- {
- // set the output format
- s.setf(std::ios_base::scientific, std::ios_base::floatfield);
- int oldprec = s.precision();
- s.precision(precision);
-
- for (int i=0; i<currentsize; i++)
- {
- s << "FEM"; // start a new row
- s << " "; // space in front of each entry
- s.width(4); // set width for counter
- s << i; // number of first entry in a line
- for (int j=0; j<currentsize; j++)
- {
- s << " "; // space in front of each entry
- s.width(width); // set width for each entry anew
- s << A[i][j]; // yeah, the number !
- }
- s << " "; // space in front of each entry
- s.width(width); // set width for each entry anew
- s << b[i];
- s << " "; // space in front of each entry
- s.width(width); // set width for each entry anew
- s << bctype[i][0];
- s << std::endl; // start a new line
- }
-
-
- // reset the output format
- s.precision(oldprec);
- s.setf(std::ios_base::fixed, std::ios_base::floatfield);
- }
-
+#if 0
//! access local stiffness matrix
/*! Access elements of the local stiffness matrix. Elements are
undefined without prior call to the assemble method.
@@ -545,17 +479,21 @@ namespace Dune
{
return bctype[i];
}
+#endif
+
private:
+
// parameters given in constructor
- //const GroundwaterEquationParameters<G,RT>& problem;
bool procBoundaryAsDirichlet;
+#if 0
// assembled data
int currentsize;
MBlockType A[SIZE][SIZE];
VBlockType b[SIZE];
BCBlockType bctype[SIZE];
+#endif
};
/** @} */
--
GitLab