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> 12 #include <dune/geometry/quadraturerules.hh> 60 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
61 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, M & mat)
const 64 using LFSU_SUB = TypeTree::Child<LFSU,0>;
65 using RF =
typename M::value_type;
66 using JacobianType =
typename LFSU_SUB::Traits::FiniteElementType::
67 Traits::LocalBasisType::Traits::JacobianType;
68 using size_type =
typename LFSU_SUB::Traits::SizeType;
71 const int dim = EG::Entity::dimension;
72 const int dimw = EG::Geometry::coorddimension;
73 static_assert(dim == dimw,
"doesn't work on manifolds");
76 const auto& cell = eg.entity();
79 auto geo = eg.geometry();
82 typename EG::Geometry::JacobianInverseTransposed jac;
85 std::vector<JacobianType> js(lfsu.child(0).size());
86 std::vector<FieldVector<RF,dim> > gradphi(lfsu.child(0).size());
92 lfsu.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
95 jac = geo.jacobianInverseTransposed(qp.position());
96 for (size_type i=0; i<lfsu.child(0).size(); i++)
99 jac.umv(js[i][0],gradphi[i]);
103 auto mu =
param_.mu(cell,qp.position());
104 auto lambda =
param_.lambda(cell,qp.position());
107 auto factor = qp.weight() * geo.integrationElement(qp.position());
109 for(
int d=0; d<
dim; ++d)
111 for (size_type i=0; i<lfsu.child(0).size(); i++)
113 for (
int k=0; k<
dim; k++)
115 for (size_type j=0; j<lfsv.child(k).size(); j++)
118 mat.accumulate(lfsv.child(k),j,lfsu.child(k),i,
120 mu * gradphi[i][d] * gradphi[j][d]
122 mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
124 mu * gradphi[i][k] * gradphi[j][d]
127 mat.accumulate(lfsv.child(k),j,lfsu.child(d),i,
128 lambda * gradphi[i][d]*gradphi[j][k]
138 template<
typename EG,
typename LFSU_HAT,
typename X,
typename LFSV,
typename R>
139 void alpha_volume (
const EG& eg,
const LFSU_HAT& lfsu_hat,
const X& x,
const LFSV& lfsv, R& r)
const 142 using LFSU = TypeTree::Child<LFSU_HAT,0>;
143 using RF =
typename R::value_type;
144 using JacobianType =
typename LFSU::Traits::FiniteElementType::
145 Traits::LocalBasisType::Traits::JacobianType;
146 using size_type =
typename LFSU::Traits::SizeType;
149 const int dim = EG::Entity::dimension;
150 const int dimw = EG::Geometry::coorddimension;
151 static_assert(dim == dimw,
"doesn't work on manifolds");
154 const auto& cell = eg.entity();
157 auto geo = eg.geometry();
160 typename EG::Geometry::JacobianInverseTransposed jac;
163 std::vector<JacobianType> js(lfsu_hat.child(0).size());
164 std::vector<FieldVector<RF,dim> > gradphi(lfsu_hat.child(0).size());
165 Dune::FieldVector<RF,dim> gradu(0.0);
171 lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(qp.position(),js);
174 jac = geo.jacobianInverseTransposed(qp.position());
175 for (size_type i=0; i<lfsu_hat.child(0).size(); i++)
178 jac.umv(js[i][0],gradphi[i]);
182 auto mu =
param_.mu(cell,qp.position());
183 auto lambda =
param_.lambda(cell,qp.position());
186 auto factor = qp.weight() * geo.integrationElement(qp.position());
188 for(
int d=0; d<
dim; ++d)
190 const LFSU & lfsu = lfsu_hat.child(d);
194 for (
size_t i=0; i<lfsu.size(); i++)
196 gradu.axpy(x(lfsu,i),gradphi[i]);
199 for (size_type i=0; i<lfsv.child(d).size(); i++)
201 for (
int k=0; k<
dim; k++)
204 r.accumulate(lfsv.child(d),i,
206 mu * gradu[k] * gradphi[i][k]
208 r.accumulate(lfsv.child(k),i,
210 mu * gradu[k] * gradphi[i][d]
213 r.accumulate(lfsv.child(k),i,
214 lambda * gradu[d]*gradphi[i][k]
223 template<
typename EG,
typename LFSV_HAT,
typename R>
227 using LFSV = TypeTree::Child<LFSV_HAT,0>;
228 using RF =
typename R::value_type;
229 using RangeType =
typename LFSV::Traits::FiniteElementType::
230 Traits::LocalBasisType::Traits::RangeType;
231 using size_type =
typename LFSV::Traits::SizeType;
234 const int dim = EG::Entity::dimension;
237 const auto& cell = eg.entity();
240 auto geo = eg.geometry();
243 std::vector<RangeType> phi(lfsv_hat.child(0).size());
244 FieldVector<RF,dim> y(0.0);
250 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(qp.position(),phi);
254 param_.f(cell,qp.position(),y);
257 auto factor = qp.weight() * geo.integrationElement(qp.position());
259 for(
int d=0; d<
dim; ++d)
261 const auto& lfsv = lfsv_hat.child(d);
264 for (size_type i=0; i<lfsv.size(); i++)
265 r.accumulate(lfsv,i, -y[d]*phi[i] * factor);
271 template<
typename IG,
typename LFSV_HAT,
typename R>
275 using namespace Indices;
276 using LFSV = TypeTree::Child<LFSV_HAT,0>;
277 using RF =
typename R::value_type;
278 using RangeType =
typename LFSV::Traits::FiniteElementType::
279 Traits::LocalBasisType::Traits::RangeType;
280 using size_type =
typename LFSV::Traits::SizeType;
283 const int dim = IG::Entity::dimension;
286 auto geo = ig.geometry();
289 auto geo_in_inside = ig.geometryInInside();
292 std::vector<RangeType> phi(lfsv_hat.child(0).size());
293 FieldVector<RF,dim> y(0.0);
299 auto local = geo_in_inside.global(qp.position());
303 if(
param_.isDirichlet( ig.intersection(), qp.position() ) )
307 lfsv_hat.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
315 auto factor = qp.weight() * geo.integrationElement(qp.position());
317 for(
int d=0; d<
dim; ++d)
319 const auto& lfsv = lfsv_hat.child(d);
322 for (size_type i=0; i<lfsv.size(); i++)
323 r.accumulate(lfsv,i, y[d]*phi[i] * factor);
337 #endif // DUNE_PDELAB_LOCALOPERATOR_LINEARELASTICITY_HH void alpha_volume(const EG &eg, const LFSU_HAT &lfsu_hat, const X &x, const LFSV &lfsv, R &r) const
Definition: linearelasticity.hh:139
const IG & ig
Definition: constraints.hh:148
Definition: linearelasticity.hh:40
static const int dim
Definition: adaptivity.hh:83
Definition: linearelasticity.hh:53
T ParameterType
Definition: linearelasticity.hh:46
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearelasticity.hh:61
void lambda_boundary(const IG &ig, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:272
int intorder_
Definition: linearelasticity.hh:329
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const P & p
Definition: constraints.hh:147
Definition: linearelasticity.hh:49
void lambda_volume(const EG &eg, const LFSV_HAT &lfsv_hat, R &r) const
Definition: linearelasticity.hh:224
Definition: linearelasticity.hh:54
Definition: linearelasticity.hh:52
sparsity pattern generator
Definition: pattern.hh:13
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:18
LinearElasticity(const ParameterType &p, int intorder=4)
Definition: linearelasticity.hh:56
const ParameterType & param_
Definition: linearelasticity.hh:330