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 
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  using ParameterType = T;
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  // Define types
62  using LFSU_SUB = TypeTree::Child<LFSU,0>;
63  using RF = typename M::value_type;
64  using JacobianType = typename LFSU_SUB::Traits::FiniteElementType::
65  Traits::LocalBasisType::Traits::JacobianType;
66  using size_type = typename LFSU_SUB::Traits::SizeType;
67 
68  // dimensions
69  const int dim = EG::Entity::dimension;
70  const int dimw = EG::Geometry::coorddimension;
71  static_assert(dim == dimw, "doesn't work on manifolds");
72 
73  // Reference to cell
74  const auto& cell = eg.entity();
75 
76  // get geometry
77  auto geo = eg.geometry();
78 
79  // Transformation
80  typename EG::Geometry::JacobianInverseTransposed jac;
81 
82  // Initialize vectors outside for loop
83  std::vector<JacobianType> js(lfsu.child(0).size());
84  std::vector<FieldVector<RF,dim> > gradphi(lfsu.child(0).size());
85 
86  // loop over quadrature points
87  for (const auto& qp : quadratureRule(geo,intorder_))
88  {
89  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
90  lfsu.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
91 
92  // transform gradient to real element
93  jac = geo.jacobianInverseTransposed(qp.position());
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  auto mu = param_.mu(cell,qp.position());
102  auto lambda = param_.lambda(cell,qp.position());
103 
104  // geometric weight
105  auto factor = qp.weight() * geo.integrationElement(qp.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  // Define types
140  using LFSU = TypeTree::Child<LFSU_HAT,0>;
141  using RF = typename R::value_type;
142  using JacobianType = typename LFSU::Traits::FiniteElementType::
143  Traits::LocalBasisType::Traits::JacobianType;
144  using size_type = typename LFSU::Traits::SizeType;
145 
146  // dimensions
147  const int dim = EG::Entity::dimension;
148  const int dimw = EG::Geometry::coorddimension;
149  static_assert(dim == dimw, "doesn't work on manifolds");
150 
151  // Reference to cell
152  const auto& cell = eg.entity();
153 
154  // Get geometry
155  auto geo = eg.geometry();
156 
157  // Transformation
158  typename EG::Geometry::JacobianInverseTransposed jac;
159 
160  // Initialize vectors outside for loop
161  std::vector<JacobianType> js(lfsu_hat.child(0).size());
162  std::vector<FieldVector<RF,dim> > gradphi(lfsu_hat.child(0).size());
163  Dune::FieldVector<RF,dim> gradu(0.0);
164 
165  // loop over quadrature points
166  for (const auto& qp : quadratureRule(geo,intorder_))
167  {
168  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
169  lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
170 
171  // transform gradient to real element
172  jac = geo.jacobianInverseTransposed(qp.position());
173  for (size_type i=0; i<lfsu_hat.child(0).size(); i++)
174  {
175  gradphi[i] = 0.0;
176  jac.umv(js[i][0],gradphi[i]);
177  }
178 
179  // material parameters
180  auto mu = param_.mu(cell,qp.position());
181  auto lambda = param_.lambda(cell,qp.position());
182 
183  // geometric weight
184  auto factor = qp.weight() * geo.integrationElement(qp.position());
185 
186  for(int d=0; d<dim; ++d)
187  {
188  const LFSU & lfsu = lfsu_hat.child(d);
189 
190  // compute gradient of u
191  gradu = 0.0;
192  for (size_t i=0; i<lfsu.size(); i++)
193  {
194  gradu.axpy(x(lfsu,i),gradphi[i]);
195  }
196 
197  for (size_type i=0; i<lfsv.child(d).size(); i++)
198  {
199  for (int k=0; k<dim; k++)
200  {
201  // integrate \mu (grad u + (grad u)^T) * (grad phi_i + (grad phi_i)^T)
202  r.accumulate(lfsv.child(d),i,
203  // mu (d u_d / d x_k) (d phi_i_d / d x_k)
204  mu * gradu[k] * gradphi[i][k]
205  *factor);
206  r.accumulate(lfsv.child(k),i,
207  // mu (d u_d / d x_k) (d phi_i_k / d x_d)
208  mu * gradu[k] * gradphi[i][d]
209  *factor);
210  // integrate \lambda sum_(k=0..dim) (d u / d x_d) * (d phi_i / d x_k)
211  r.accumulate(lfsv.child(k),i,
212  lambda * gradu[d]*gradphi[i][k]
213  *factor);
214  }
215  }
216  }
217  }
218  }
219 
220  // volume integral depending only on test functions
221  template<typename EG, typename LFSV_HAT, typename R>
222  void lambda_volume (const EG& eg, const LFSV_HAT& lfsv_hat, R& r) const
223  {
224  // Define types
225  using LFSV = TypeTree::Child<LFSV_HAT,0>;
226  using RF = typename R::value_type;
227  using RangeType = typename LFSV::Traits::FiniteElementType::
228  Traits::LocalBasisType::Traits::RangeType;
229  using size_type = typename LFSV::Traits::SizeType;
230 
231  // dimensions
232  const int dim = EG::Entity::dimension;
233 
234  // Reference to cell
235  const auto& cell = eg.entity();
236 
237  // Get geometry
238  auto geo = eg.geometry();
239 
240  // Initialize vectors outside for loop
241  std::vector<RangeType> phi(lfsv_hat.child(0).size());
242  FieldVector<RF,dim> y(0.0);
243 
244  // loop over quadrature points
245  for (const auto& qp : quadratureRule(geo,intorder_))
246  {
247  // evaluate shape functions
248  lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(qp.position(),phi);
249 
250  // evaluate right hand side parameter function
251  y = 0.0;
252  param_.f(cell,qp.position(),y);
253 
254  // weight
255  auto factor = qp.weight() * geo.integrationElement(qp.position());
256 
257  for(int d=0; d<dim; ++d)
258  {
259  const auto& 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  // Define types
273  using namespace Indices;
274  using LFSV = TypeTree::Child<LFSV_HAT,0>;
275  using RF = typename R::value_type;
276  using RangeType = typename LFSV::Traits::FiniteElementType::
277  Traits::LocalBasisType::Traits::RangeType;
278  using size_type = typename LFSV::Traits::SizeType;
279 
280  // dimensions
281  const int dim = IG::Entity::dimension;
282 
283  // get geometry
284  auto geo = ig.geometry();
285 
286  // Get geometry of intersection in local coordinates of inside cell
287  auto geo_in_inside = ig.geometryInInside();
288 
289  // Initialize vectors outside for loop
290  std::vector<RangeType> phi(lfsv_hat.child(0).size());
291  FieldVector<RF,dim> y(0.0);
292 
293  // loop over quadrature points
294  for (const auto& qp : quadratureRule(geo,intorder_))
295  {
296  // position of quadrature point in local coordinates of element
297  auto local = geo_in_inside.global(qp.position());
298 
299  // evaluate boundary condition type
300  // skip rest if we are on Dirichlet boundary
301  if( param_.isDirichlet( ig.intersection(), qp.position() ) )
302  continue;
303 
304  // evaluate shape functions
305  lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
306 
307  // evaluate surface force
308  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  auto factor = qp.weight() * geo.integrationElement(qp.position());
314 
315  for(int d=0; d<dim; ++d)
316  {
317  const auto& 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 // DUNE_PDELAB_LOCALOPERATOR_LINEARELASTICITY_HH
Dune::PDELab::LinearElasticity::intorder_
int intorder_
Definition: linearelasticity.hh:327
flags.hh
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
dim
static const int dim
Definition: adaptivity.hh:84
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:18
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
defaultimp.hh
Dune::PDELab::LinearElasticity::param_
const ParameterType & param_
Definition: linearelasticity.hh:328
Dune::PDELab::LinearElasticity::doPatternVolume
@ doPatternVolume
Definition: linearelasticity.hh:47
idefault.hh
Dune::PDELab::LinearElasticity::doAlphaVolume
@ doAlphaVolume
Definition: linearelasticity.hh:50
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:13
Dune::PDELab::LinearElasticity::alpha_volume
void alpha_volume(const EG &eg, const LFSU_HAT &lfsu_hat, const X &x, const LFSV &lfsv, R &r) const
Definition: linearelasticity.hh:137
Dune::PDELab::LinearElasticity::jacobian_volume
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearelasticity.hh:59
ig
const IG & ig
Definition: constraints.hh:148
Dune::PDELab::LinearElasticity
Definition: linearelasticity.hh:38
Dune::PDELab::quadratureRule
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
pattern.hh
Dune::PDELab::LinearElasticity::lambda_boundary
void lambda_boundary(const IG &ig, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:270
linearelasticityparameter.hh
quadraturerules.hh
Dune::PDELab::LinearElasticity::lambda_volume
void lambda_volume(const EG &eg, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:222
Dune::PDELab::LinearElasticity::ParameterType
T ParameterType
Definition: linearelasticity.hh:44
Dune::PDELab::LinearElasticity::doLambdaBoundary
@ doLambdaBoundary
Definition: linearelasticity.hh:52
Dune::PDELab::LinearElasticity::LinearElasticity
LinearElasticity(const ParameterType &p, int intorder=4)
Definition: linearelasticity.hh:54
Dune::PDELab::LinearElasticity::doLambdaVolume
@ doLambdaVolume
Definition: linearelasticity.hh:51
p
const P & p
Definition: constraints.hh:147