dune-pdelab  2.0.0
linearelasticity.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; c-basic-offset: 2; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LINEARELASTICITY_HH
3 #define DUNE_PDELAB_LINEARELASTICITY_HH
4 
5 #include <vector>
6 
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/static_assert.hh>
10 
11 #include <dune/geometry/type.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules.hh>
14 
15 #include "defaultimp.hh"
16 #include "pattern.hh"
17 #include "flags.hh"
18 #include "idefault.hh"
19 
21 
22 namespace Dune {
23  namespace PDELab {
27 
37  template<typename T>
40  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::DomainType>
41  {
42  public:
43 
44  typedef T ParameterType;
45 
46  // pattern assembly flags
47  enum { doPatternVolume = true };
48 
49  // residual assembly flags
50  enum { doAlphaVolume = true };
51  enum { doLambdaVolume = true };
52  enum { doLambdaBoundary = true };
53 
54  LinearElasticity (const ParameterType & p, int intorder=4)
55  : intorder_(intorder), param_(p)
56  {}
57 
58  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
59  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & mat) const
60  {
61  // extract local function spaces
62  typedef typename LFSU::template Child<0>::Type LFSU_SUB;
63 
64  // domain and range field type
65  typedef typename LFSU_SUB::Traits::FiniteElementType::
66  Traits::LocalBasisType::Traits::DomainFieldType DF;
67  typedef typename LFSU_SUB::Traits::FiniteElementType::
68  Traits::LocalBasisType::Traits::RangeFieldType RF;
69  typedef typename LFSU_SUB::Traits::FiniteElementType::
70  Traits::LocalBasisType::Traits::JacobianType JacobianType;
71 
72  typedef typename LFSU_SUB::Traits::SizeType size_type;
73 
74  // dimensions
75  const int dim = EG::Geometry::dimension;
76  const int dimw = EG::Geometry::dimensionworld;
77  dune_static_assert(dim == dimw, "doesn't work on manifolds");
78 
79  // select quadrature rule
80  GeometryType gt = eg.geometry().type();
81  const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,intorder_);
82 
83  // loop over quadrature points
84  for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
85  {
86  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
87  std::vector<JacobianType> js(lfsu.child(0).size());
88  lfsu.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
89 
90  // transform gradient to real element
91  const typename EG::Geometry::JacobianInverseTransposed jac
92  = eg.geometry().jacobianInverseTransposed(it->position());
93  std::vector<FieldVector<RF,dim> > gradphi(lfsu.child(0).size());
94  for (size_type i=0; i<lfsu.child(0).size(); i++)
95  {
96  gradphi[i] = 0.0;
97  jac.umv(js[i][0],gradphi[i]);
98  }
99 
100  // material parameters
101  RF mu = param_.mu(eg.entity(),it->position());
102  RF lambda = param_.lambda(eg.entity(),it->position());
103 
104  // geometric weight
105  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
106 
107  for(int d=0; d<dim; ++d)
108  {
109  for (size_type i=0; i<lfsu.child(0).size(); i++)
110  {
111  for (int k=0; k<dim; k++)
112  {
113  for (size_type j=0; j<lfsv.child(k).size(); j++)
114  {
115  // integrate \mu (grad u + (grad u)^T) * (grad phi_i + (grad phi_i)^T)
116  mat.accumulate(lfsv.child(k),j,lfsu.child(k),i,
117  // mu (d u_k / d x_d) (d v_k / d x_d)
118  mu * gradphi[i][d] * gradphi[j][d]
119  *factor);
120  mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
121  // mu (d u_d / d x_k) (d v_k / d x_d)
122  mu * gradphi[i][k] * gradphi[j][d]
123  *factor);
124  // integrate \lambda sum_(k=0..dim) (d u_d / d x_d) * (d v_k / d x_k)
125  mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
126  lambda * gradphi[i][d]*gradphi[j][k]
127  *factor);
128  }
129  }
130  }
131  }
132  }
133  }
134 
135  // volume integral depending on test and ansatz functions
136  template<typename EG, typename LFSU_HAT, typename X, typename LFSV, typename R>
137  void alpha_volume (const EG& eg, const LFSU_HAT& lfsu_hat, const X& x, const LFSV& lfsv, R& r) const
138  {
139  // extract local function spaces
140  typedef typename LFSU_HAT::template Child<0>::Type LFSU;
141 
142  // domain and range field type
143  typedef typename LFSU::Traits::FiniteElementType::
144  Traits::LocalBasisType::Traits::DomainFieldType DF;
145  typedef typename LFSU::Traits::FiniteElementType::
146  Traits::LocalBasisType::Traits::RangeFieldType RF;
147  typedef typename LFSU::Traits::FiniteElementType::
148  Traits::LocalBasisType::Traits::JacobianType JacobianType;
149 
150  typedef typename LFSU::Traits::SizeType size_type;
151 
152  // dimensions
153  const int dim = EG::Geometry::dimension;
154  const int dimw = EG::Geometry::dimensionworld;
155  dune_static_assert(dim == dimw, "doesn't work on manifolds");
156 
157  // select quadrature rule
158  GeometryType gt = eg.geometry().type();
159  const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,intorder_);
160 
161  // loop over quadrature points
162  for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
163  {
164  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
165  std::vector<JacobianType> js(lfsu_hat.child(0).size());
166  lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
167 
168  // transform gradient to real element
169  const typename EG::Geometry::JacobianInverseTransposed jac
170  = eg.geometry().jacobianInverseTransposed(it->position());
171  std::vector<FieldVector<RF,dim> > gradphi(lfsu_hat.child(0).size());
172  for (size_type i=0; i<lfsu_hat.child(0).size(); i++)
173  {
174  gradphi[i] = 0.0;
175  jac.umv(js[i][0],gradphi[i]);
176  }
177 
178  // material parameters
179  RF mu = param_.mu(eg.entity(),it->position());
180  RF lambda = param_.lambda(eg.entity(),it->position());
181 
182  // geometric weight
183  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
184 
185  for(int d=0; d<dim; ++d)
186  {
187  const LFSU & lfsu = lfsu_hat.child(d);
188 
189  // compute gradient of u
190  Dune::FieldVector<RF,dim> gradu(0.0);
191  for (size_t i=0; i<lfsu.size(); i++)
192  {
193  gradu.axpy(x(lfsu,i),gradphi[i]);
194  }
195 
196  for (size_type i=0; i<lfsv.child(d).size(); i++)
197  {
198  for (int k=0; k<dim; k++)
199  {
200  // integrate \mu (grad u + (grad u)^T) * (grad phi_i + (grad phi_i)^T)
201  r.accumulate(lfsv.child(d),i,
202  // mu (d u_d / d x_k) (d phi_i_d / d x_k)
203  mu * gradu[k] * gradphi[i][k]
204  *factor);
205  r.accumulate(lfsv.child(k),i,
206  // mu (d u_d / d x_k) (d phi_i_k / d x_d)
207  mu * gradu[k] * gradphi[i][d]
208  *factor);
209  // integrate \lambda sum_(k=0..dim) (d u / d x_d) * (d phi_i / d x_k)
210  r.accumulate(lfsv.child(k),i,
211  lambda * gradu[d]*gradphi[i][k]
212  *factor);
213  }
214  }
215  }
216  }
217  }
218 
219  // volume integral depending only on test functions
220  template<typename EG, typename LFSV_HAT, typename R>
221  void lambda_volume (const EG& eg, const LFSV_HAT& lfsv_hat, R& r) const
222  {
223  // extract local function spaces
224  typedef typename LFSV_HAT::template Child<0>::Type LFSV;
225 
226  // domain and range field type
227  typedef typename LFSV::Traits::FiniteElementType::
228  Traits::LocalBasisType::Traits::DomainFieldType DF;
229  typedef typename LFSV::Traits::FiniteElementType::
230  Traits::LocalBasisType::Traits::RangeFieldType RF;
231  typedef typename LFSV::Traits::FiniteElementType::
232  Traits::LocalBasisType::Traits::RangeType RangeType;
233 
234  typedef typename LFSV::Traits::SizeType size_type;
235 
236  // dimensions
237  const int dim = EG::Geometry::dimension;
238 
239  // select quadrature rule
240  GeometryType gt = eg.geometry().type();
241  const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,intorder_);
242 
243  // loop over quadrature points
244  for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
245  {
246  // evaluate shape functions
247  std::vector<RangeType> phi(lfsv_hat.child(0).size());
248  lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
249 
250  // evaluate right hand side parameter function
251  FieldVector<RF,dim> y(0.0);
252  param_.f(eg.entity(),it->position(),y);
253 
254  // weight
255  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
256 
257  for(int d=0; d<dim; ++d)
258  {
259  const LFSV & lfsv = lfsv_hat.child(d);
260 
261  // integrate f
262  for (size_type i=0; i<lfsv.size(); i++)
263  r.accumulate(lfsv,i, -y[d]*phi[i] * factor);
264  }
265  }
266  }
267 
268  // jacobian of boundary term
269  template<typename IG, typename LFSV_HAT, typename R>
270  void lambda_boundary (const IG& ig, const LFSV_HAT& lfsv_hat, R& r) const
271  {
272  // extract local function spaces
273  typedef typename LFSV_HAT::template Child<0>::Type LFSV;
274 
275  // domain and range field type
276  typedef typename LFSV::Traits::FiniteElementType::
277  Traits::LocalBasisType::Traits::DomainFieldType DF;
278  typedef typename LFSV::Traits::FiniteElementType::
279  Traits::LocalBasisType::Traits::RangeFieldType RF;
280  typedef typename LFSV::Traits::FiniteElementType::
281  Traits::LocalBasisType::Traits::RangeType RangeType;
282 
283  typedef typename LFSV::Traits::SizeType size_type;
284 
285  // dimensions
286  const int dim = IG::Geometry::dimension;
287 
288  // select quadrature rule
289  GeometryType gt = ig.geometry().type();
290  const QuadratureRule<DF,dim-1>& rule = QuadratureRules<DF,dim-1>::rule(gt,intorder_);
291 
292  // loop over quadrature points
293  for (typename QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
294  {
295  // position of quadrature point in local coordinates of element
296  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
297 
298  // evaluate boundary condition type
299  // skip rest if we are on Dirichlet boundary
300  if( param_.isDirichlet( ig.intersection(), it->position() ) )
301  continue;
302 
303  // evaluate shape functions
304  std::vector<RangeType> phi(lfsv_hat.child(0).size());
305  lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
306 
307  // evaluate surface force
308  FieldVector<RF,dim> y(0.0);
309  // currently we only implement homogeneous Neumann (e.g. Stress) BC
310  // param_.g(eg.entity(),it->position(),y);
311 
312  // weight
313  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
314 
315  for(int d=0; d<dim; ++d)
316  {
317  const LFSV & lfsv = lfsv_hat.child(d);
318 
319  // integrate f
320  for (size_type i=0; i<lfsv.size(); i++)
321  r.accumulate(lfsv,i, y[d]*phi[i] * factor);
322  }
323  }
324  }
325 
326  protected:
329  };
330 
332  } // namespace PDELab
333 } // namespace Dune
334 
335 #endif
Definition: linearelasticity.hh:51
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearelasticity.hh:59
Definition: linearelasticity.hh:47
T ParameterType
Definition: linearelasticity.hh:44
Default flags for all local operators.
Definition: flags.hh:18
static const int dim
Definition: adaptivity.hh:82
Definition: linearelasticity.hh:50
void lambda_boundary(const IG &ig, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:270
const IG & ig
Definition: common/constraints.hh:146
int intorder_
Definition: linearelasticity.hh:327
LinearElasticity(const ParameterType &p, int intorder=4)
Definition: linearelasticity.hh:54
const ParameterType & param_
Definition: linearelasticity.hh:328
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: linearelasticity.hh:52
sparsity pattern generator
Definition: pattern.hh:14
const EG & eg
Definition: common/constraints.hh:277
void alpha_volume(const EG &eg, const LFSU_HAT &lfsu_hat, const X &x, const LFSV &lfsv, R &r) const
Definition: linearelasticity.hh:137
R p(int i, D x)
Definition: qkdg.hh:62
void lambda_volume(const EG &eg, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:221
Definition: linearelasticity.hh:38