2 #ifndef DUNE_PDELAB_LINEARELASTICITY_HH
3 #define DUNE_PDELAB_LINEARELASTICITY_HH
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
10 #include <dune/geometry/type.hh>
11 #include <dune/geometry/referenceelements.hh>
12 #include <dune/geometry/quadraturerules.hh>
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
61 typedef typename LFSU::template Child<0>::Type LFSU_SUB;
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;
70 typedef typename LFSU_SUB::Traits::SizeType size_type;
73 const int dim = EG::Entity::dimension;
74 const int dimw = EG::Geometry::coorddimension;
75 static_assert(dim == dimw,
"doesn't work on manifolds");
78 const auto& geometry = eg.geometry();
79 GeometryType gt = geometry.type();
80 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,
intorder_);
83 for (
const auto& qp : rule)
86 std::vector<JacobianType> js(lfsu.child(0).size());
87 lfsu.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
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++)
96 jac.umv(js[i][0],gradphi[i]);
100 RF mu =
param_.mu(eg.entity(),qp.position());
101 RF lambda =
param_.lambda(eg.entity(),qp.position());
104 RF factor = qp.weight() * geometry.integrationElement(qp.position());
106 for(
int d=0; d<
dim; ++d)
108 for (size_type i=0; i<lfsu.child(0).size(); i++)
110 for (
int k=0; k<
dim; k++)
112 for (size_type j=0; j<lfsv.child(k).size(); j++)
115 mat.accumulate(lfsv.child(k),j,lfsu.child(k),i,
117 mu * gradphi[i][d] * gradphi[j][d]
119 mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
121 mu * gradphi[i][k] * gradphi[j][d]
124 mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
125 lambda * gradphi[i][d]*gradphi[j][k]
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
139 typedef typename LFSU_HAT::template Child<0>::Type LFSU;
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;
148 typedef typename LFSU::Traits::SizeType size_type;
151 const int dim = EG::Entity::dimension;
152 const int dimw = EG::Geometry::coorddimension;
153 static_assert(dim == dimw,
"doesn't work on manifolds");
156 const auto& geometry = eg.geometry();
157 GeometryType gt = geometry.type();
158 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,
intorder_);
161 for (
const auto& qp : rule)
164 std::vector<JacobianType> js(lfsu_hat.child(0).size());
165 lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
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++)
174 jac.umv(js[i][0],gradphi[i]);
178 RF mu =
param_.mu(eg.entity(),qp.position());
179 RF lambda =
param_.lambda(eg.entity(),qp.position());
182 RF factor = qp.weight() * geometry.integrationElement(qp.position());
184 for(
int d=0; d<
dim; ++d)
186 const LFSU & lfsu = lfsu_hat.child(d);
189 Dune::FieldVector<RF,dim> gradu(0.0);
190 for (
size_t i=0; i<lfsu.size(); i++)
192 gradu.axpy(x(lfsu,i),gradphi[i]);
195 for (size_type i=0; i<lfsv.child(d).size(); i++)
197 for (
int k=0; k<
dim; k++)
200 r.accumulate(lfsv.child(d),i,
202 mu * gradu[k] * gradphi[i][k]
204 r.accumulate(lfsv.child(k),i,
206 mu * gradu[k] * gradphi[i][d]
209 r.accumulate(lfsv.child(k),i,
210 lambda * gradu[d]*gradphi[i][k]
219 template<
typename EG,
typename LFSV_HAT,
typename R>
223 typedef typename LFSV_HAT::template Child<0>::Type LFSV;
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;
232 typedef typename LFSV::Traits::SizeType size_type;
235 const int dim = EG::Entity::dimension;
238 const auto& geometry = eg.geometry();
239 GeometryType gt = geometry.type();
240 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,
intorder_);
243 for (
const auto& qp : rule)
246 std::vector<RangeType> phi(lfsv_hat.child(0).size());
247 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(qp.position(),phi);
250 FieldVector<RF,dim> y(0.0);
251 param_.f(eg.entity(),qp.position(),y);
254 RF factor = qp.weight() * geometry.integrationElement(qp.position());
256 for(
int d=0; d<
dim; ++d)
258 const LFSV & lfsv = lfsv_hat.child(d);
261 for (size_type i=0; i<lfsv.size(); i++)
262 r.accumulate(lfsv,i, -y[d]*phi[i] * factor);
268 template<
typename IG,
typename LFSV_HAT,
typename R>
272 typedef typename LFSV_HAT::template Child<0>::Type LFSV;
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;
281 typedef typename LFSV::Traits::SizeType size_type;
284 const int dim = IG::Entity::dimension;
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_);
293 for (
const auto& qp : rule)
296 Dune::FieldVector<DF,dim> local = geometryInInside.global(qp.position());
300 if(
param_.isDirichlet( ig.intersection(), qp.position() ) )
304 std::vector<RangeType> phi(lfsv_hat.child(0).size());
305 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
308 FieldVector<RF,dim> y(0.0);
313 RF factor = qp.weight() * geometry.integrationElement(qp.position());
315 for(
int d=0; d<
dim; ++d)
317 const LFSV & lfsv = lfsv_hat.child(d);
320 for (size_type i=0; i<lfsv.size(); i++)
321 r.accumulate(lfsv,i, y[d]*phi[i] * factor);
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