2 #ifndef DUNE_PDELAB_LINEARELASTICITY_HH
3 #define DUNE_PDELAB_LINEARELASTICITY_HH
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/static_assert.hh>
11 #include <dune/geometry/type.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules.hh>
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
62 typedef typename LFSU::template Child<0>::Type LFSU_SUB;
65 typedef typename LFSU_SUB::Traits::FiniteElementType::
66 Traits::LocalBasisType::Traits::DomainFieldType DF;
67 typedef typename LFSU_SUB::Traits::FiniteElementType::
68 Traits::LocalBasisType::Traits::RangeFieldType RF;
69 typedef typename LFSU_SUB::Traits::FiniteElementType::
70 Traits::LocalBasisType::Traits::JacobianType JacobianType;
72 typedef typename LFSU_SUB::Traits::SizeType size_type;
75 const int dim = EG::Geometry::dimension;
76 const int dimw = EG::Geometry::dimensionworld;
77 dune_static_assert(dim == dimw,
"doesn't work on manifolds");
80 GeometryType gt = eg.geometry().type();
81 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,
intorder_);
84 for (
typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
87 std::vector<JacobianType> js(lfsu.child(0).size());
88 lfsu.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
91 const typename EG::Geometry::JacobianInverseTransposed jac
92 = eg.geometry().jacobianInverseTransposed(it->position());
93 std::vector<FieldVector<RF,dim> > gradphi(lfsu.child(0).size());
94 for (size_type i=0; i<lfsu.child(0).size(); i++)
97 jac.umv(js[i][0],gradphi[i]);
101 RF mu =
param_.mu(eg.entity(),it->position());
102 RF lambda =
param_.lambda(eg.entity(),it->position());
105 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
107 for(
int d=0; d<
dim; ++d)
109 for (size_type i=0; i<lfsu.child(0).size(); i++)
111 for (
int k=0; k<
dim; k++)
113 for (size_type j=0; j<lfsv.child(k).size(); j++)
116 mat.accumulate(lfsv.child(k),j,lfsu.child(k),i,
118 mu * gradphi[i][d] * gradphi[j][d]
120 mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
122 mu * gradphi[i][k] * gradphi[j][d]
125 mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
126 lambda * gradphi[i][d]*gradphi[j][k]
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
140 typedef typename LFSU_HAT::template Child<0>::Type LFSU;
143 typedef typename LFSU::Traits::FiniteElementType::
144 Traits::LocalBasisType::Traits::DomainFieldType DF;
145 typedef typename LFSU::Traits::FiniteElementType::
146 Traits::LocalBasisType::Traits::RangeFieldType RF;
147 typedef typename LFSU::Traits::FiniteElementType::
148 Traits::LocalBasisType::Traits::JacobianType JacobianType;
150 typedef typename LFSU::Traits::SizeType size_type;
153 const int dim = EG::Geometry::dimension;
154 const int dimw = EG::Geometry::dimensionworld;
155 dune_static_assert(dim == dimw,
"doesn't work on manifolds");
158 GeometryType gt = eg.geometry().type();
159 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,
intorder_);
162 for (
typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
165 std::vector<JacobianType> js(lfsu_hat.child(0).size());
166 lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
169 const typename EG::Geometry::JacobianInverseTransposed jac
170 = eg.geometry().jacobianInverseTransposed(it->position());
171 std::vector<FieldVector<RF,dim> > gradphi(lfsu_hat.child(0).size());
172 for (size_type i=0; i<lfsu_hat.child(0).size(); i++)
175 jac.umv(js[i][0],gradphi[i]);
179 RF mu =
param_.mu(eg.entity(),it->position());
180 RF lambda =
param_.lambda(eg.entity(),it->position());
183 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
185 for(
int d=0; d<
dim; ++d)
187 const LFSU & lfsu = lfsu_hat.child(d);
190 Dune::FieldVector<RF,dim> gradu(0.0);
191 for (
size_t i=0; i<lfsu.size(); i++)
193 gradu.axpy(x(lfsu,i),gradphi[i]);
196 for (size_type i=0; i<lfsv.child(d).size(); i++)
198 for (
int k=0; k<
dim; k++)
201 r.accumulate(lfsv.child(d),i,
203 mu * gradu[k] * gradphi[i][k]
205 r.accumulate(lfsv.child(k),i,
207 mu * gradu[k] * gradphi[i][d]
210 r.accumulate(lfsv.child(k),i,
211 lambda * gradu[d]*gradphi[i][k]
220 template<
typename EG,
typename LFSV_HAT,
typename R>
224 typedef typename LFSV_HAT::template Child<0>::Type LFSV;
227 typedef typename LFSV::Traits::FiniteElementType::
228 Traits::LocalBasisType::Traits::DomainFieldType DF;
229 typedef typename LFSV::Traits::FiniteElementType::
230 Traits::LocalBasisType::Traits::RangeFieldType RF;
231 typedef typename LFSV::Traits::FiniteElementType::
232 Traits::LocalBasisType::Traits::RangeType RangeType;
234 typedef typename LFSV::Traits::SizeType size_type;
237 const int dim = EG::Geometry::dimension;
240 GeometryType gt = eg.geometry().type();
241 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,
intorder_);
244 for (
typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
247 std::vector<RangeType> phi(lfsv_hat.child(0).size());
248 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
251 FieldVector<RF,dim> y(0.0);
252 param_.f(eg.entity(),it->position(),y);
255 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
257 for(
int d=0; d<
dim; ++d)
259 const LFSV & lfsv = lfsv_hat.child(d);
262 for (size_type i=0; i<lfsv.size(); i++)
263 r.accumulate(lfsv,i, -y[d]*phi[i] * factor);
269 template<
typename IG,
typename LFSV_HAT,
typename R>
273 typedef typename LFSV_HAT::template Child<0>::Type LFSV;
276 typedef typename LFSV::Traits::FiniteElementType::
277 Traits::LocalBasisType::Traits::DomainFieldType DF;
278 typedef typename LFSV::Traits::FiniteElementType::
279 Traits::LocalBasisType::Traits::RangeFieldType RF;
280 typedef typename LFSV::Traits::FiniteElementType::
281 Traits::LocalBasisType::Traits::RangeType RangeType;
283 typedef typename LFSV::Traits::SizeType size_type;
286 const int dim = IG::Geometry::dimension;
289 GeometryType gt = ig.geometry().type();
290 const QuadratureRule<DF,dim-1>& rule = QuadratureRules<DF,dim-1>::rule(gt,
intorder_);
293 for (
typename QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
296 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
300 if(
param_.isDirichlet( ig.intersection(), it->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 = it->weight() * ig.geometry().integrationElement(it->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:51
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearelasticity.hh:59
Definition: linearelasticity.hh:47
T ParameterType
Definition: linearelasticity.hh:44
Default flags for all local operators.
Definition: flags.hh:18
static const int dim
Definition: adaptivity.hh:82
Definition: linearelasticity.hh:50
void lambda_boundary(const IG &ig, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:270
const IG & ig
Definition: common/constraints.hh:146
int intorder_
Definition: linearelasticity.hh:327
LinearElasticity(const ParameterType &p, int intorder=4)
Definition: linearelasticity.hh:54
const ParameterType & param_
Definition: linearelasticity.hh:328
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: linearelasticity.hh:52
sparsity pattern generator
Definition: pattern.hh:14
const EG & eg
Definition: common/constraints.hh:277
void alpha_volume(const EG &eg, const LFSU_HAT &lfsu_hat, const X &x, const LFSV &lfsv, R &r) const
Definition: linearelasticity.hh:137
R p(int i, D x)
Definition: qkdg.hh:62
void lambda_volume(const EG &eg, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:221
Definition: linearelasticity.hh:38