2 #ifndef DUNE_PDELAB_LOCALOPERATOR_LINEARELASTICITY_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_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>
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 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;
69 const int dim = EG::Entity::dimension;
70 const int dimw = EG::Geometry::coorddimension;
71 static_assert(
dim == dimw,
"doesn't work on manifolds");
74 const auto& cell = eg.entity();
77 auto geo = eg.geometry();
80 typename EG::Geometry::JacobianInverseTransposed jac;
83 std::vector<JacobianType> js(lfsu.child(0).size());
84 std::vector<FieldVector<RF,dim> > gradphi(lfsu.child(0).size());
90 lfsu.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
93 jac = geo.jacobianInverseTransposed(qp.position());
94 for (size_type i=0; i<lfsu.child(0).size(); i++)
97 jac.umv(js[i][0],gradphi[i]);
101 auto mu =
param_.mu(cell,qp.position());
102 auto lambda =
param_.lambda(cell,qp.position());
105 auto factor = qp.weight() * geo.integrationElement(qp.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 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;
147 const int dim = EG::Entity::dimension;
148 const int dimw = EG::Geometry::coorddimension;
149 static_assert(
dim == dimw,
"doesn't work on manifolds");
152 const auto& cell = eg.entity();
155 auto geo = eg.geometry();
158 typename EG::Geometry::JacobianInverseTransposed jac;
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);
169 lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
172 jac = geo.jacobianInverseTransposed(qp.position());
173 for (size_type i=0; i<lfsu_hat.child(0).size(); i++)
176 jac.umv(js[i][0],gradphi[i]);
180 auto mu =
param_.mu(cell,qp.position());
181 auto lambda =
param_.lambda(cell,qp.position());
184 auto factor = qp.weight() * geo.integrationElement(qp.position());
186 for(
int d=0; d<
dim; ++d)
188 const LFSU & lfsu = lfsu_hat.child(d);
192 for (
size_t i=0; i<lfsu.size(); i++)
194 gradu.axpy(x(lfsu,i),gradphi[i]);
197 for (size_type i=0; i<lfsv.child(d).size(); i++)
199 for (
int k=0; k<
dim; k++)
202 r.accumulate(lfsv.child(d),i,
204 mu * gradu[k] * gradphi[i][k]
206 r.accumulate(lfsv.child(k),i,
208 mu * gradu[k] * gradphi[i][d]
211 r.accumulate(lfsv.child(k),i,
212 lambda * gradu[d]*gradphi[i][k]
221 template<
typename EG,
typename LFSV_HAT,
typename R>
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;
232 const int dim = EG::Entity::dimension;
235 const auto& cell = eg.entity();
238 auto geo = eg.geometry();
241 std::vector<RangeType> phi(lfsv_hat.child(0).size());
242 FieldVector<RF,dim> y(0.0);
248 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(qp.position(),phi);
252 param_.f(cell,qp.position(),y);
255 auto factor = qp.weight() * geo.integrationElement(qp.position());
257 for(
int d=0; d<
dim; ++d)
259 const auto& 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 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;
281 const int dim = IG::Entity::dimension;
284 auto geo =
ig.geometry();
287 auto geo_in_inside =
ig.geometryInInside();
290 std::vector<RangeType> phi(lfsv_hat.child(0).size());
291 FieldVector<RF,dim> y(0.0);
297 auto local = geo_in_inside.global(qp.position());
301 if(
param_.isDirichlet(
ig.intersection(), qp.position() ) )
305 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
313 auto factor = qp.weight() * geo.integrationElement(qp.position());
315 for(
int d=0; d<
dim; ++d)
317 const auto& 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);
static const int dim
Definition: adaptivity.hh:84
const P & p
Definition: constraints.hh:148
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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:19
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Definition: linearelasticity.hh:41
const ParameterType & param_
Definition: linearelasticity.hh:328
@ doLambdaBoundary
Definition: linearelasticity.hh:52
void lambda_boundary(const IG &ig, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:270
@ doLambdaVolume
Definition: linearelasticity.hh:51
LinearElasticity(const ParameterType &p, int intorder=4)
Definition: linearelasticity.hh:54
@ doPatternVolume
Definition: linearelasticity.hh:47
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearelasticity.hh:59
@ doAlphaVolume
Definition: linearelasticity.hh:50
void alpha_volume(const EG &eg, const LFSU_HAT &lfsu_hat, const X &x, const LFSV &lfsv, R &r) const
Definition: linearelasticity.hh:137
void lambda_volume(const EG &eg, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:222
T ParameterType
Definition: linearelasticity.hh:44
int intorder_
Definition: linearelasticity.hh:327
sparsity pattern generator
Definition: pattern.hh:14