DUNE PDELab (git)

linearelasticity.hh
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
9
10#include <dune/geometry/type.hh>
11#include <dune/geometry/referenceelements.hh>
12
13#include<dune/pdelab/common/quadraturerules.hh>
14
15#include "defaultimp.hh"
16#include "pattern.hh"
17#include "flags.hh"
18#include "idefault.hh"
19
20#include "linearelasticityparameter.hh"
21
22namespace 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:
327 int intorder_;
328 const ParameterType & param_;
329 };
330
332 } // namespace PDELab
333} // namespace Dune
334
335#endif // DUNE_PDELAB_LOCALOPERATOR_LINEARELASTICITY_HH
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
vector space out of a tensor product of fields.
Definition: fvector.hh:91
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Definition: linearelasticity.hh:41
Default flags for all local operators.
Definition: flags.hh:19
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
typename impl::_Child< Node, indices... >::type Child
Template alias for the type of a child node given by a list of child indices.
Definition: childextraction.hh:225
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)