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