4 #ifndef DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
11 #include <dune/common/fvector.hh>
12 #include <dune/common/static_assert.hh>
14 #include <dune/geometry/type.hh>
15 #include <dune/geometry/quadraturerules.hh>
30 namespace VectorWave {
46 template<
class Imp,
class GV,
class RF =
double,
class Time_ =
double>
57 typedef FieldVector<DomainField, dimension>
Domain;
62 typedef FieldVector<RangeField, dimension>
Range;
65 typedef typename GV::template Codim<0>::Entity
Element;
77 {
return asImp().epsilonGlobal(e.geometry().global(xl)); }
94 {
return asImp().muGlobal(e.geometry().global(xl)); }
110 const Imp &asImp()
const {
return *
static_cast<const Imp*
>(
this); }
123 template<
class GV,
class RF =
double,
class Time =
double>
125 public Parameters<ConstantParameters<GV, RF, Time>, GV, RF, Time>
133 template<
typename Domain>
136 template<
typename Domain>
166 template<
class Params>
192 R0(Params ¶ms_, std::size_t qorder_ = 0) :
193 params(params_), qorder(qorder_)
196 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
199 const LFSV& lfsv, M& mat)
const
202 typedef typename LFSU::Traits::FiniteElementType FEU;
203 typedef typename LFSV::Traits::FiniteElementType FEV;
205 const FEU &feU = lfsu.finiteElement();
206 const FEV &feV = lfsv.finiteElement();
208 typedef typename FEU::Traits::Basis BasisU;
209 typedef typename FEV::Traits::Basis BasisV;
211 const BasisU &basisU = feU.basis();
212 const BasisV &basisV = feV.basis();
214 typedef typename Params::RangeField RF;
216 typedef typename BasisU::Traits::DomainField DF;
217 static const std::size_t dimDLocal = BasisU::Traits::dimDomainLocal;
219 typedef typename BasisU::Traits::Jacobian JacobianU;
220 typedef typename BasisV::Traits::Jacobian JacobianV;
225 static const J2CU &j2CU = J2CU();
226 static const J2CV &j2CV = J2CV();
228 typedef typename J2CU::Curl CurlU;
229 typedef typename J2CV::Curl CurlV;
231 dune_static_assert(J2CU::dimCurl == J2CV::dimCurl,
"Curl dimension "
232 "of ansatz and test functions must match in "
236 typedef QuadratureRules<DF,dimDLocal> QRs;
237 typedef QuadratureRule<DF,dimDLocal> QR;
238 typedef typename QR::const_iterator QRIterator;
239 GeometryType gt = eg.geometry().type();
240 const QR& rule = QRs::rule(gt,qorder);
242 std::vector<JacobianU> jacobianU(lfsu.size());
243 std::vector<JacobianV> jacobianV(lfsv.size());
245 std::vector<CurlU> curlU(lfsu.size());
246 std::vector<CurlV> curlV(lfsv.size());
249 for(QRIterator it=rule.begin(); it != rule.end(); ++it) {
251 basisU.evaluateJacobian(it->position(), jacobianU);
252 basisV.evaluateJacobian(it->position(), jacobianV);
254 for(std::size_t j = 0; j < lfsu.size(); ++j)
255 j2CU(jacobianU[j], curlU[j]);
256 for(std::size_t i = 0; i < lfsv.size(); ++i)
257 j2CV(jacobianV[i], curlV[i]);
259 RF factor = it->weight()
260 * eg.geometry().integrationElement(it->position())
261 / params.mu(eg.entity(), it->position());
263 for(std::size_t j = 0; j < lfsu.size(); ++j)
264 for(std::size_t i = 0; i < lfsv.size(); ++i)
265 mat(i,j) += factor * (curlU[j] * curlV[i]);
271 params.setTime(time);
295 template<
class Params>
309 R1(Params ¶ms_, std::size_t qorder_ = 0) {}
335 template<
class Params>
361 R2(Params ¶ms_, std::size_t qorder_ = 2)
362 : params(params_), qorder(qorder_)
365 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
368 const LFSV& lfsv, M& mat)
const
371 typedef typename LFSU::Traits::FiniteElementType FEU;
372 typedef typename LFSV::Traits::FiniteElementType FEV;
374 const FEU &feU = lfsu.finiteElement();
375 const FEV &feV = lfsv.finiteElement();
377 typedef typename FEU::Traits::Basis BasisU;
378 typedef typename FEV::Traits::Basis BasisV;
380 const BasisU &basisU = feU.basis();
381 const BasisV &basisV = feV.basis();
383 typedef typename Params::RangeField RF;
385 typedef typename BasisU::Traits::DomainField DF;
386 static const std::size_t dimDLocal = BasisU::Traits::dimDomainLocal;
388 typedef typename BasisU::Traits::Range RangeU;
389 typedef typename BasisV::Traits::Range RangeV;
392 typedef QuadratureRules<DF,dimDLocal> QRs;
393 typedef QuadratureRule<DF,dimDLocal> QR;
394 typedef typename QR::const_iterator QRIterator;
395 GeometryType gt = eg.geometry().type();
396 const QR& rule = QRs::rule(gt,qorder);
398 std::vector<RangeU> phiU(lfsu.size());
399 std::vector<RangeV> phiV(lfsv.size());
402 for(QRIterator it=rule.begin(); it != rule.end(); ++it) {
404 basisU.evaluateFunction(it->position(), phiU);
405 basisV.evaluateFunction(it->position(), phiV);
407 RF factor = it->weight()
408 * eg.geometry().integrationElement(it->position())
409 * params.epsilon(eg.entity(), it->position());
411 for(std::size_t j = 0; j < lfsu.size(); ++j)
412 for(std::size_t i = 0; i < lfsv.size(); ++i)
413 mat(i,j) += factor * (phiU[j] * phiV[i]);
419 params.setTime(time);
425 template<
class Params,
class Step =
int,
class Time =
double>
430 const Params ¶ms;
442 virtual bool isAffine(std::size_t order, Step step)
const
451 Step requested, Step available)
const
455 if(t1 > t2) std::swap(t1, t2);
458 case 0:
return !params.muChanged(t1, t2);
460 case 2:
return !params.epsilonChanged(t1, t2);
463 DUNE_THROW(InvalidStateException,
464 "VectorWave::CachePolicy::canReuseJacobian(): Invalid "
465 "temporal derivative order=" << order);
468 Step requested, Step available)
const
471 Step available)
const
475 if(t1 > t2) std::swap(t1, t2);
477 return !params.muChanged(t1, t2) && !params.epsilonChanged(t1, t2);
487 #endif // DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
Time_ Time
export type of temporal values
Definition: vectorwave.hh:49
RangeField mu(const Element &e, const Domain &xl) const
evaluate magnetic permeability
Definition: vectorwave.hh:93
bool epsilonChanged(Time t1, Time t2) const
Definition: vectorwave.hh:135
Local operator for the vector wave problem, one-temporal-derivative part.
Definition: vectorwave.hh:296
R1(Params ¶ms_, std::size_t qorder_=0)
Construct a local operator object.
Definition: vectorwave.hh:309
ConstantParameters(RF epsilon, RF mu)
Definition: vectorwave.hh:131
Default flags for all local operators.
Definition: flags.hh:18
virtual bool canReuseJacobian(std::size_t order, Step requested, Step available) const
Whether a jacobian can be reused.
Definition: vectorwave.hh:450
CachePolicy(const Params ¶ms_)
construct a CachePolicy object
Definition: vectorwave.hh:439
void setTime(typename Params::Time time)
set time on the function object
Definition: vectorwave.hh:418
void setTime(typename Params::Time time)
set time on the parameter object
Definition: vectorwave.hh:270
Definition: vectorwave.hh:350
Time dt
Time step size.
Definition: cache.hh:75
Parameter class interface for the vector wave local operators.
Definition: vectorwave.hh:47
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
Definition: vectorwave.hh:182
RF RangeField
field type of range
Definition: vectorwave.hh:60
Step currentStep
The number of the step set via preStep().
Definition: cache.hh:61
virtual bool isComposedAffine(Step step) const
The composed operator is affine.
Definition: vectorwave.hh:445
bool muChanged(Time t1, Time t2) const
Definition: vectorwave.hh:138
MultiStepCachePolicy for VectorWave operators.
Definition: vectorwave.hh:426
Local operator for the vector wave problem, no-temporal-derivatives part.
Definition: vectorwave.hh:167
RF epsilonGlobal(const Domain &) const
Definition: vectorwave.hh:134
bool muChanged(Time t1, Time t2) const
Whether mu changes between the given time steps.
Definition: vectorwave.hh:101
FieldVector< RangeField, dimension > Range
vector type of range
Definition: vectorwave.hh:62
Policy class for the MultiStepCache.
Definition: cache.hh:46
virtual bool hasPureLinearAlpha(std::size_t order, Step step) const
All local operators have pure linear alpha_*()
Definition: vectorwave.hh:448
extract the curl of a function from the jacobian of that function
Definition: jacobiantocurl.hh:28
R0(Params ¶ms_, std::size_t qorder_=0)
Construct a local operator object.
Definition: vectorwave.hh:192
GV::ctype DomainField
field type of domain
Definition: vectorwave.hh:55
bool epsilonChanged(Time t1, Time t2) const
Whether epsilon changes between the given time steps.
Definition: vectorwave.hh:84
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Time endTime
The end-time of the step set via preStep().
Definition: cache.hh:68
GV::template Codim< 0 >::Entity Element
element type of GridView
Definition: vectorwave.hh:65
RangeField epsilon(const Element &e, const Domain &xl) const
evaluate dielectric permittivity
Definition: vectorwave.hh:76
Local operator for the vector wave problem, second-temporal-derivatives part.
Definition: vectorwave.hh:336
R2(Params ¶ms_, std::size_t qorder_=2)
Construct a local operator object.
Definition: vectorwave.hh:361
static const std::size_t dimension
export dimension (both domain and range)
Definition: vectorwave.hh:52
sparsity pattern generator
Definition: pattern.hh:14
virtual bool canReuseZeroResidual(std::size_t order, Step requested, Step available) const
Whether a zero-residual can be reused.
Definition: vectorwave.hh:467
void setTime(const Time &time)
set the time for subsequent evaluation
Definition: vectorwave.hh:107
const EG & eg
Definition: common/constraints.hh:277
Homogenous parameter class for the vector wave local operators.
Definition: vectorwave.hh:124
Definition: vectorwave.hh:351
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: vectorwave.hh:198
RF muGlobal(const Domain &) const
Definition: vectorwave.hh:137
const E & e
Definition: interpolate.hh:172
FieldVector< DomainField, dimension > Domain
vector type of domain
Definition: vectorwave.hh:57
virtual bool canReuseComposedJacobian(Step requested, Step available) const
Whether a composed jacobian can be reused.
Definition: vectorwave.hh:470
Definition: vectorwave.hh:181
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: vectorwave.hh:367
virtual bool isAffine(std::size_t order, Step step) const
All component operators are affine.
Definition: vectorwave.hh:442
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104