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