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