dune-fem 2.8.0
Loading...
Searching...
No Matches
dgelliptic.hh
Go to the documentation of this file.
1/**************************************************************************
2
3 The dune-fem module is a module of DUNE (see www.dune-project.org).
4 It is based on the dune-grid interface library
5 extending the grid interface by a number of discretization algorithms
6 for solving non-linear systems of partial differential equations.
7
8 Copyright (C) 2003 - 2015 Robert Kloefkorn
9 Copyright (C) 2003 - 2010 Mario Ohlberger
10 Copyright (C) 2004 - 2015 Andreas Dedner
11 Copyright (C) 2005 Adrian Burri
12 Copyright (C) 2005 - 2015 Mirko Kraenkel
13 Copyright (C) 2006 - 2015 Christoph Gersbacher
14 Copyright (C) 2006 - 2015 Martin Nolte
15 Copyright (C) 2011 - 2015 Tobias Malkmus
16 Copyright (C) 2012 - 2015 Stefan Girke
17 Copyright (C) 2013 - 2015 Claus-Justus Heine
18 Copyright (C) 2013 - 2014 Janick Gerstenberger
19 Copyright (C) 2013 Sven Kaulman
20 Copyright (C) 2013 Tom Ranner
21 Copyright (C) 2015 Marco Agnese
22 Copyright (C) 2015 Martin Alkaemper
23
24
25 The dune-fem module is free software; you can redistribute it and/or
26 modify it under the terms of the GNU General Public License as
27 published by the Free Software Foundation; either version 2 of
28 the License, or (at your option) any later version.
29
30 The dune-fem module is distributed in the hope that it will be useful,
31 but WITHOUT ANY WARRANTY; without even the implied warranty of
32 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 GNU General Public License for more details.
34
35 You should have received a copy of the GNU General Public License along
36 with this program; if not, write to the Free Software Foundation, Inc.,
37 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
38
39**************************************************************************/
40#ifndef DGELLIPTIC_HH
41#define DGELLIPTIC_HH
42
43#include <dune/common/fmatrix.hh>
44
50
52
53// EllipticOperator
54// ----------------
55
56template <class DFSpace>
58{
59 DefaultPenalty(const DFSpace &space, double penalty)
60 : space_(space)
61 , penalty_(penalty)
62 {}
63 template <class Intersection>
64 double operator()(const Intersection &intersection,
65 double intersectionArea, double area, double nbArea) const
66 {
67 const double hInv = intersectionArea / std::min( area, nbArea );
68 return penalty_ * hInv;
69 }
70 const double &factor() const { return penalty_; }
71 private:
72 const DFSpace &space_;
73 double penalty_;
74};
75
76template< class DiscreteFunction, class Model,
79: public virtual Dune::Fem::Operator< DiscreteFunction >
80{
81 typedef DiscreteFunction DiscreteFunctionType;
82 typedef Model ModelType;
83
84 typedef DiscreteFunction DomainDiscreteFunctionType;
85 typedef DiscreteFunction RangeDiscreteFunctionType;
86protected:
87 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
88 typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
89 typedef typename LocalFunctionType::RangeType RangeType;
90 typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
91
92 typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
93 typedef typename IteratorType::Entity EntityType;
94 typedef typename EntityType::Geometry GeometryType;
95
96 typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
97
98 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
99 typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
100 typedef typename IntersectionIteratorType::Intersection IntersectionType;
101 typedef typename IntersectionType::Geometry IntersectionGeometryType;
102
105
106 static const int dimDomain = LocalFunctionType::dimDomain;
107 static const int dimRange = LocalFunctionType::dimRange;
108
109public:
112 const ModelType &model,
114 : model_( model ),
115 penalty_( space, parameter.getValue< double >( "penalty", 40 ) )
116 {
117 // std::cout << "dg operator with penalty:" << penalty_.factor() << std::endl;
118 }
119
121 const DiscreteFunctionSpaceType &rSpace,
122 const ModelType &model,
124 : model_( model ),
125 penalty_( dSpace, parameter.getValue< double >( "penalty", 40 ) )
126 {
127 // std::cout << "dg operator with penalty:" << penalty_.factor() << std::endl;
128 }
129
132 { apply(u,w); }
133 template <class GF>
134 void apply( const GF &u, RangeDiscreteFunctionType &w ) const;
135
136protected:
137 const ModelType &model () const { return model_; }
138 Penalty penalty() const { return penalty_; }
139
140private:
141 const ModelType &model_;
142 Penalty penalty_;
143};
144
145// DifferentiableDGEllipticOperator
146// ------------------------------
147
148template< class JacobianOperator, class Model,
151: public DGEllipticOperator< typename JacobianOperator::DomainFunctionType, Model, Penalty >,
152 public Dune::Fem::DifferentiableOperator< JacobianOperator >
153{
155
156 typedef JacobianOperator JacobianOperatorType;
157
162
163protected:
164 typedef typename DomainDiscreteFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
165 typedef typename DomainDiscreteFunctionType::LocalFunctionType DomainLocalFunctionType;
166 typedef typename DomainLocalFunctionType::RangeType DomainRangeType;
167 typedef typename DomainLocalFunctionType::JacobianRangeType DomainJacobianRangeType;
168 typedef typename RangeDiscreteFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
169 typedef typename RangeDiscreteFunctionType::LocalFunctionType RangeLocalFunctionType;
170 typedef typename RangeLocalFunctionType::RangeType RangeRangeType;
171 typedef typename RangeLocalFunctionType::JacobianRangeType RangeJacobianRangeType;
172
173 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
174 typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
175 typedef typename LocalFunctionType::RangeType RangeType;
176 typedef typename LocalFunctionType::RangeFieldType RangeFieldType;
177 typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
178
179 typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
180 typedef typename IteratorType::Entity EntityType;
181 typedef typename EntityType::Geometry GeometryType;
182
183 typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
184
185 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
186 typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
187 typedef typename IntersectionIteratorType::Intersection IntersectionType;
188 typedef typename IntersectionType::Geometry IntersectionGeometryType;
189
192
193 static const int dimDomain = LocalFunctionType::dimDomain;
194 static const int dimRange = LocalFunctionType::dimRange;
195
196public:
199 const ModelType &model,
201 : BaseType( space, model, parameter )
202 {}
204 const DiscreteFunctionSpaceType &rSpace,
205 const ModelType &model,
207 : BaseType( dSpace, rSpace, model, parameter )
208 {}
209
212 { apply(u,jOp); }
213 template <class GridFunctionType>
214 void apply ( const GridFunctionType &u, JacobianOperatorType &jOp ) const;
216
217protected:
220};
221
222// Implementation of DGEllipticOperator
223// ----------------------------------
224
225template< class RangeDiscreteFunction, class Model, class Penalty >
226template<class GF>
228 ::apply( const GF &u, RangeDiscreteFunctionType &w ) const
229{
230 // clear destination
231 w.clear();
232
233 // get discrete function space
234 const DiscreteFunctionSpaceType &dfSpace = w.space();
235
236 // iterate over grid
237 const IteratorType end = dfSpace.end();
238 for( IteratorType it = dfSpace.begin(); it != end; ++it )
239 {
240 // get entity (here element)
241 const EntityType &entity = *it;
242
243 bool needsCalculation = model().init( entity );
244 if (! needsCalculation )
245 continue;
246
247 // get elements geometry
248 const GeometryType &geometry = entity.geometry();
249
250 // get local representation of the discrete functions
251 const auto uLocal = u.localFunction( entity );
252 LocalFunctionType wLocal = w.localFunction( entity );
253
254 // obtain quadrature order
255 const int quadOrder = uLocal.order() + wLocal.order();
256
257 { // element integral
258 QuadratureType quadrature( entity, quadOrder );
259 const size_t numQuadraturePoints = quadrature.nop();
260 for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
261 {
262 const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
263 const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
264
265 RangeType vu;
266 uLocal.evaluate( quadrature[ pt ], vu );
268 uLocal.jacobian( quadrature[ pt ], du );
269
270 // compute mass contribution (studying linear case so linearizing around zero)
271 RangeType avu( 0 );
272 model().source( quadrature[ pt ], vu, du, avu );
273 avu *= weight;
274 // add to local functional wLocal.axpy( quadrature[ pt ], avu );
275
276 JacobianRangeType adu( 0 );
277 // apply diffusive flux
278 model().flux( quadrature[ pt ], vu, du, adu );
279 adu *= weight;
280
281 // add to local function
282 wLocal.axpy( quadrature[ pt ], avu, adu );
283 }
284 }
285 if ( ! dfSpace.continuous() )
286 {
287 const double area = entity.geometry().volume();
288 const IntersectionIteratorType iitend = dfSpace.gridPart().iend( entity );
289 for( IntersectionIteratorType iit = dfSpace.gridPart().ibegin( entity ); iit != iitend; ++iit ) // looping over intersections
290 {
292 const IntersectionType &intersection = *iit;
293 if ( intersection.neighbor() )
294 {
295 // const EntityType outside = Dune::Fem::make_entity( intersection.outside() );
296 const EntityType outside = intersection.outside() ;
297
298 typedef typename IntersectionType::Geometry IntersectionGeometryType;
299 const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
300
301 // compute penalty factor
302 const double intersectionArea = intersectionGeometry.volume();
303 const double beta = penalty()(intersection,intersectionArea,area,outside.geometry().volume());
304
305 auto uOutLocal = u.localFunction( outside ); // local u on outisde element
306
307 FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
308 FaceQuadratureType quadOutside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::OUTSIDE );
309 const size_t numQuadraturePoints = quadInside.nop();
311
312 for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
313 {
315 // get coordinate of quadrature point on the reference element of the intersection
316 const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
317 DomainType normal = intersection.integrationOuterNormal( x );
318 double faceVol = normal.two_norm();
319 normal /= faceVol; // make it into a unit normal
320 const double weight = quadInside.weight( pt ) * faceVol;
321
322 RangeType vuIn,vuOut;
323 JacobianRangeType duIn, duOut;
324 uLocal.evaluate( quadInside[ pt ], vuIn );
325 uLocal.jacobian( quadInside[ pt ], duIn );
326 uOutLocal.evaluate( quadOutside[ pt ], vuOut );
327 uOutLocal.jacobian( quadOutside[ pt ], duOut );
328 RangeType jump = vuIn - vuOut;
329 // compute -0.5 * [u] x normal
330 JacobianRangeType dvalue;
331 for (int r=0;r<dimRange;++r)
332 for (int d=0;d<dimDomain;++d)
333 dvalue[r][d] = -0.5 * normal[d] * jump[r];
334 JacobianRangeType aduIn,aduOut;
335 model().init( outside );
336 model().flux( quadOutside[ pt ], vuOut, duOut, aduOut );
337 auto pFactor = model().penalty( quadOutside[ pt ], vuOut );
338 model().init( entity );
339 model().flux( quadInside[ pt ], vuIn, duIn, aduIn );
340 JacobianRangeType affine;
341 model().flux( quadInside[ pt ], RangeType(0), JacobianRangeType(0), affine);
342 pFactor += model().penalty( quadInside[ pt ], vuIn );
344
346 RangeType value;
347 JacobianRangeType advalue;
348 // penalty term : beta [u] [phi] = beta (u+ - u-)(phi+ - phi-)=beta (u+ - u-)phi+
349 value = jump;
350 for (unsigned int r=0;r<dimRange;++r)
351 value[r] *= beta * pFactor[r]/2.;
352 // {A grad u}.[phi] = {A grad u}.phi+ n_+ = 0.5*(grad u+ + grad u-).n_+ phi+
353 aduIn += aduOut;
354 aduIn *= -0.5;
355 aduIn.umv(normal,value);
356 // [ u ] * { {A} grad phi_en } = -normal(u+ - u-) * 0.5 {A} grad phi_en
357 // we actually need G(u)tau with G(x,u) = d/sigma_j D_i(x,u,sigma)
358 // - we might need to assume D(x,u,sigma) = G(x,u)sigma + affine(x)
359 model().flux( quadInside[ pt ], vuIn, dvalue, advalue );
360 // advalue -= affine;
361
362 value *= weight;
363 advalue *= weight;
364 wLocal.axpy( quadInside[ pt ], value, advalue );
366 }
367 }
368 else if( intersection.boundary() )
369 {
370 std::array<int,dimRange> components;
371 components.fill(0);
372 model().isDirichletIntersection( intersection, components);
373
374 typedef typename IntersectionType::Geometry IntersectionGeometryType;
375 const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
376
377 // compute penalty factor
378 const double intersectionArea = intersectionGeometry.volume();
379 const double beta = penalty()(intersection,intersectionArea,area,area);
380
381 FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
382 const size_t numQuadraturePoints = quadInside.nop();
383
384 for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
385 {
386 const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
387 const DomainType normal = intersection.integrationOuterNormal( x );
388 const double weight = quadInside.weight( pt );
389
390 RangeType bndValue;
391 model().dirichlet(1, quadInside[pt], bndValue);
392
393 RangeType value;
394 JacobianRangeType dvalue,advalue;
395
396 RangeType vuIn,jump;
397 JacobianRangeType duIn, aduIn;
398 uLocal.evaluate( quadInside[ pt ], vuIn );
399 uLocal.jacobian( quadInside[ pt ], duIn );
400 jump = vuIn;
401 jump -= bndValue;
402
403 // penalty term : beta [u] [phi] = beta (u+ - u-)(phi+ - phi-)=beta (u+ - u-)phi+
404 auto pFactor = model().penalty( quadInside[ pt ], vuIn );
405 value = jump;
406 for (unsigned int r=0;r<dimRange;++r)
407 value[r] *= beta * pFactor[r] * intersectionGeometry.integrationElement( x );
408
409 // [ u ] * { grad phi_en } = -normal(u+ - u-) * 0.5 grad phi_en
410 // here we need a diadic product of u x n
411 for (int r=0;r<dimRange;++r)
412 for (int d=0;d<dimDomain;++d)
413 dvalue[r][d] = -0.5 * normal[d] * jump[r];
414
415 model().flux( quadInside[ pt ], jump, dvalue, advalue );
416
417 // consistency term
418 // {A grad u}.[phi] = {A grad u}.phi+ n_+ = 0.5*(grad u+ + grad u-).n_+ phi+
419 model().flux( quadInside[ pt ], vuIn, duIn, aduIn );
420 aduIn.mmv(normal,value);
421
422 for (int r=0;r<dimRange;++r)
423 if (!components[r]) // do not use dirichlet constraints here
424 {
425 value[r] = 0;
426 advalue[r] = 0;
427 }
428
429 value *= weight;
430 advalue *= weight;
431 wLocal.axpy( quadInside[ pt ], value, advalue );
432 }
433 }
434 }
435 }
436
437 }
438
439 // communicate data (in parallel runs)
440 w.communicate();
441}
442
443// Implementation of DifferentiableDGEllipticOperator
444// ------------------------------------------------
445template< class JacobianOperator, class Model, class Penalty >
446template<class GF>
448 ::apply ( const GF &u, JacobianOperator &jOp ) const
449{
450 // std::cout << "starting assembly\n";
451 // Dune::Timer timer;
452 typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
453 typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;
454
455 const DomainDiscreteFunctionSpaceType &domainSpace = jOp.domainSpace();
456 const RangeDiscreteFunctionSpaceType &rangeSpace = jOp.rangeSpace();
457
459 jOp.reserve(stencil);
460 jOp.clear();
461
462 const GridPartType& gridPart = rangeSpace.gridPart();
463
464 const unsigned int numDofs = rangeSpace.blockMapper().maxNumDofs() *
465 DiscreteFunctionSpaceType :: localBlockSize ;
466
467 std::vector< RangeType > phi( numDofs );
468 std::vector< JacobianRangeType > dphi( numDofs );
469
470 std::vector< RangeType > phiNb( numDofs );
471 std::vector< JacobianRangeType > dphiNb( numDofs );
472
473 // std::cout << " in assembly: start element loop size=" << rangeSpace.gridPart().grid().size(0) << " time= " << timer.elapsed() << std::endl;;
474 const IteratorType end = rangeSpace.end();
475 for( IteratorType it = rangeSpace.begin(); it != end; ++it )
476 {
477 const EntityType &entity = *it;
478
479 bool needsCalculation = model().init( entity );
480 if (! needsCalculation )
481 continue;
482
483 const GeometryType geometry = entity.geometry();
484
485 const auto uLocal = u.localFunction( entity );
486
487 LocalMatrixType jLocal = jOp.localMatrix( entity, entity );
488
489 const BasisFunctionSetType &baseSet = jLocal.domainBasisFunctionSet();
490 const unsigned int numBaseFunctions = baseSet.size();
491
492 QuadratureType quadrature( entity, 2*rangeSpace.order() );
493 const size_t numQuadraturePoints = quadrature.nop();
494 for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
495 {
496 const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
497 const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
498
499 // evaluate all basis functions at given quadrature point
500 baseSet.evaluateAll( quadrature[ pt ], phi );
501
502 // evaluate jacobians of all basis functions at given quadrature point
503 baseSet.jacobianAll( quadrature[ pt ], dphi );
504
505 // get value for linearization
506 RangeType u0;
507 JacobianRangeType jacU0;
508 uLocal.evaluate( quadrature[ pt ], u0 );
509 uLocal.jacobian( quadrature[ pt ], jacU0 );
510
511 RangeType aphi( 0 );
512 JacobianRangeType adphi( 0 );
513 for( unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
514 {
515 // if mass terms or right hand side is present
516 model().linSource( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], aphi );
517
518 // if gradient term is present
519 model().linFlux( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], adphi );
520
521 // get column object and call axpy method
522 jLocal.column( localCol ).axpy( phi, dphi, aphi, adphi, weight );
523 }
524 }
525 if ( rangeSpace.continuous() )
526 continue;
527
528 double area = geometry.volume();
529 const IntersectionIteratorType endiit = gridPart.iend( entity );
530 for ( IntersectionIteratorType iit = gridPart.ibegin( entity );
531 iit != endiit ; ++ iit )
532 {
533 const IntersectionType& intersection = *iit ;
534
535 if( intersection.neighbor() )
536 {
537 const EntityType neighbor = Dune::Fem::make_entity( intersection.outside() );
538
539 typedef typename IntersectionType::Geometry IntersectionGeometryType;
540 const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
541
543 // get local matrix for face entries
544 LocalMatrixType localOpNb = jOp.localMatrix( neighbor, entity );
545 // get neighbor's base function set
546 const BasisFunctionSetType &baseSetNb = localOpNb.domainBasisFunctionSet();
548
549 // compute penalty factor
550 const double intersectionArea = intersectionGeometry.volume();
551 const double beta = penalty()(intersection,intersectionArea,area,neighbor.geometry().volume());
552
553 // here we assume that the intersection is conforming
554 FaceQuadratureType faceQuadInside(gridPart, intersection, 2*rangeSpace.order() + 1,
555 FaceQuadratureType::INSIDE);
556 FaceQuadratureType faceQuadOutside(gridPart, intersection, 2*rangeSpace.order() + 1,
557 FaceQuadratureType::OUTSIDE);
558
559 const size_t numFaceQuadPoints = faceQuadInside.nop();
560 for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
561 {
562 const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
563 DomainType normal = intersection.integrationOuterNormal( x );
564 double faceVol = normal.two_norm();
565 normal /= faceVol; // make it into a unit normal
566
567 const double quadWeight = faceQuadInside.weight( pt );
568 const double weight = quadWeight * faceVol;
569
571 RangeType u0En;
572 JacobianRangeType u0EnJac;
573 uLocal.evaluate( faceQuadInside[ pt ], u0En );
574 uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );
575
577 // evaluate basis function of face inside E^- (entity)
579
580 // evaluate all basis functions for quadrature point pt
581 baseSet.evaluateAll( faceQuadInside[ pt ], phi );
582
583 // evaluate the jacobians of all basis functions
584 baseSet.jacobianAll( faceQuadInside[ pt ], dphi );
585
587 // evaluate basis function of face inside E^+ (neighbor)
589
590 // evaluate all basis functions for quadrature point pt on neighbor
591 baseSetNb.evaluateAll( faceQuadOutside[ pt ], phiNb );
592
593 // evaluate the jacobians of all basis functions on neighbor
594 baseSetNb.jacobianAll( faceQuadOutside[ pt ], dphiNb );
595
596 model().init( entity );
597 for( unsigned int i = 0; i < numBaseFunctions; ++i )
598 {
599 JacobianRangeType adphiEn = dphi[ i ];
600 model().linFlux( u0En, u0EnJac, faceQuadInside[ pt ], phi[i], adphiEn, dphi[ i ] );
601 }
602 auto pFactor = model().penalty( faceQuadInside[ pt ], u0En );
603 model().init( neighbor );
604 for( unsigned int i = 0; i < numBaseFunctions; ++i )
605 {
606 JacobianRangeType adphiNb = dphiNb[ i ];
607 model().linFlux( u0En, u0EnJac, faceQuadOutside[ pt ], phiNb[i], adphiNb, dphiNb[ i ] );
608 }
609 pFactor += model().penalty( faceQuadOutside[ pt ], u0En );
610 model().init( entity );
612
614 for( unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
615 {
616 RangeType valueEn(0), valueNb(0);
617 JacobianRangeType dvalueEn(0), dvalueNb(0);
618
619 // -{ A grad u } * [ phi_en ]
620 dphi[localCol].usmv( -0.5, normal, valueEn );
621
622 // -{ A grad u } * [ phi_en ]
623 dphiNb[localCol].usmv( -0.5, normal, valueNb );
624
625 // [ u ] * [ phi_en ] = u^- * phi_en^-
626 for (unsigned int r=0;r<dimRange;++r)
627 {
628 valueEn[r] += beta*pFactor[r]/2.*phi[localCol][r];
629 valueNb[r] -= beta*pFactor[r]/2.*phiNb[localCol][r];
630 }
631 // here we need a diadic product of u x n
632 for ( int r=0; r< dimRange; ++r )
633 for ( int d=0; d< dimDomain; ++d )
634 {
635 // [ u ] * { grad phi_en }
636 dvalueEn[r][d] = - 0.5 * normal[d] * phi[localCol][r];
637
638 // [ u ] * { grad phi_en }
639 dvalueNb[r][d] = 0.5 * normal[d] * phiNb[localCol][r];
640 }
641
642 jLocal.column( localCol ).axpy( phi, dphi, valueEn, dvalueEn, weight );
643 localOpNb.column( localCol ).axpy( phi, dphi, valueNb, dvalueNb, weight );
644 }
646 }
647 }
648 else if( intersection.boundary() )
649 {
650 std::array<int,dimRange> components;
651 components.fill(0);
652 model().isDirichletIntersection( intersection, components);
653
654 typedef typename IntersectionType::Geometry IntersectionGeometryType;
655 const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
656
657 // compute penalty factor
658 const double intersectionArea = intersectionGeometry.volume();
659 const double beta = penalty()(intersection,intersectionArea,area,area);
660
661 // here we assume that the intersection is conforming
662 FaceQuadratureType faceQuadInside(gridPart, intersection, 2*rangeSpace.order() + 1,
663 FaceQuadratureType::INSIDE);
664
665 const size_t numFaceQuadPoints = faceQuadInside.nop();
666 for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
667 {
668 const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
669 DomainType normal = intersection.integrationOuterNormal( x );
670 double faceVol = normal.two_norm();
671 normal /= faceVol; // make it into a unit normal
672
673 const double quadWeight = faceQuadInside.weight( pt );
674 const double weight = quadWeight * faceVol;
675
676 RangeType u0En;
677 JacobianRangeType u0EnJac;
678 uLocal.evaluate( faceQuadInside[ pt ], u0En );
679 uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );
680
682 // evaluate basis function of face inside E^- (entity)
684
685 // evaluate all basis functions for quadrature point pt
686 baseSet.evaluateAll( faceQuadInside[ pt ], phi );
687
688 // evaluate the jacobians of all basis functions
689 baseSet.jacobianAll( faceQuadInside[ pt ], dphi );
690
691 for( unsigned int i = 0; i < numBaseFunctions; ++i )
692 {
693 JacobianRangeType adphiEn = dphi[ i ];
694 model().linFlux( u0En, u0EnJac, faceQuadInside[ pt ], phi[i], adphiEn, dphi[ i ] );
695 }
696
697 auto pFactor = model().penalty( faceQuadInside[ pt ], u0En );
698
699 for( unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
700 {
701 RangeType valueEn(0);
702 JacobianRangeType dvalueEn(0);
703
704 // -{ A grad u } * [ phi_en ]
705 dphi[localCol].usmv( -1.0, normal, valueEn );
706
707 // [ u ] * [ phi_en ] = u^- * phi_en^-
708 for (unsigned int r=0;r<dimRange;++r)
709 valueEn[r] += beta*pFactor[r]*phi[ localCol ][r];
710
711 // here we need a diadic product of u x n
712 for ( int r=0; r< dimRange; ++r )
713 for ( int d=0; d< dimDomain; ++d )
714 {
715 // [ u ] * { grad phi_en }
716 dvalueEn[r][d] = - 0.5 * normal[d] * phi[localCol][r];
717 }
718
719 for (int r=0;r<dimRange;++r)
720 if (!components[r]) // do not use dirichlet constraints here
721 {
722 valueEn[r] = 0;
723 dvalueEn[r] = 0;
724 }
725
726 jLocal.column( localCol ).axpy( phi, dphi, valueEn, dvalueEn, weight );
727 }
728 }
729 }
730 }
731 } // end grid traversal
732 jOp.flushAssembly();
733 // std::cout << " in assembly: final " << timer.elapsed() << std::endl;;
734}
735
736#endif // ELLIPTIC_HH
double min(const Dune::Fem::Double &v, const double p)
Definition: double.hh:953
Dune::Entity< codim, dim, Grid, Implementation > make_entity(Dune::Entity< codim, dim, Grid, Implementation > entity)
Definition: compatibility.hh:21
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract differentiable operator
Definition: differentiableoperator.hh:29
abstract operator
Definition: operator.hh:34
Stencil contaning the entries (en,en) and (en,nb) for all entities en in the space and neighbors nb o...
Definition: stencil.hh:385
quadrature on the codim-0 reference element
Definition: elementquadrature.hh:58
Definition: dgelliptic.hh:58
const double & factor() const
Definition: dgelliptic.hh:70
DefaultPenalty(const DFSpace &space, double penalty)
Definition: dgelliptic.hh:59
double operator()(const Intersection &intersection, double intersectionArea, double area, double nbArea) const
Definition: dgelliptic.hh:64
Definition: dgelliptic.hh:80
DiscreteFunctionSpaceType::IteratorType IteratorType
Definition: dgelliptic.hh:92
DGEllipticOperator(const DiscreteFunctionSpaceType &space, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
constructor
Definition: dgelliptic.hh:111
DiscreteFunctionSpaceType::DomainType DomainType
Definition: dgelliptic.hh:96
GridPartType::IntersectionIteratorType IntersectionIteratorType
Definition: dgelliptic.hh:99
DGEllipticOperator(const DiscreteFunctionSpaceType &dSpace, const DiscreteFunctionSpaceType &rSpace, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: dgelliptic.hh:120
static const int dimDomain
Definition: dgelliptic.hh:106
IteratorType::Entity EntityType
Definition: dgelliptic.hh:93
EntityType::Geometry GeometryType
Definition: dgelliptic.hh:94
DiscreteFunctionSpaceType::GridPartType GridPartType
Definition: dgelliptic.hh:98
void apply(const GF &u, RangeDiscreteFunctionType &w) const
Definition: dgelliptic.hh:228
Penalty penalty() const
Definition: dgelliptic.hh:138
const ModelType & model() const
Definition: dgelliptic.hh:137
IntersectionType::Geometry IntersectionGeometryType
Definition: dgelliptic.hh:101
IntersectionIteratorType::Intersection IntersectionType
Definition: dgelliptic.hh:100
LocalFunctionType::JacobianRangeType JacobianRangeType
Definition: dgelliptic.hh:90
virtual void operator()(const DomainDiscreteFunctionType &u, RangeDiscreteFunctionType &w) const
application operator
Definition: dgelliptic.hh:131
static const int dimRange
Definition: dgelliptic.hh:107
DiscreteFunction DomainDiscreteFunctionType
Definition: dgelliptic.hh:84
Dune::Fem::ElementQuadrature< GridPartType, 1 > FaceQuadratureType
Definition: dgelliptic.hh:103
LocalFunctionType::RangeType RangeType
Definition: dgelliptic.hh:89
Dune::Fem::CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: dgelliptic.hh:104
DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: dgelliptic.hh:87
DiscreteFunction RangeDiscreteFunctionType
Definition: dgelliptic.hh:85
DiscreteFunctionType::LocalFunctionType LocalFunctionType
Definition: dgelliptic.hh:88
DiscreteFunction DiscreteFunctionType
Definition: dgelliptic.hh:81
Model ModelType
Definition: dgelliptic.hh:82
Definition: dgelliptic.hh:153
RangeLocalFunctionType::JacobianRangeType RangeJacobianRangeType
Definition: dgelliptic.hh:171
static const int dimRange
Definition: dgelliptic.hh:194
DiscreteFunctionSpaceType::IteratorType IteratorType
Definition: dgelliptic.hh:179
static const int dimDomain
Definition: dgelliptic.hh:193
DiscreteFunctionSpaceType::DomainType DomainType
Definition: dgelliptic.hh:183
DomainLocalFunctionType::JacobianRangeType DomainJacobianRangeType
Definition: dgelliptic.hh:167
BaseType::ModelType ModelType
Definition: dgelliptic.hh:159
DiscreteFunctionType::LocalFunctionType LocalFunctionType
Definition: dgelliptic.hh:174
DifferentiableDGEllipticOperator(const DiscreteFunctionSpaceType &dSpace, const DiscreteFunctionSpaceType &rSpace, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: dgelliptic.hh:203
DomainDiscreteFunctionType::LocalFunctionType DomainLocalFunctionType
Definition: dgelliptic.hh:165
IntersectionType::Geometry IntersectionGeometryType
Definition: dgelliptic.hh:188
Dune::Fem::ElementQuadrature< GridPartType, 1 > FaceQuadratureType
Definition: dgelliptic.hh:190
void apply(const GridFunctionType &u, JacobianOperatorType &jOp) const
IteratorType::Entity EntityType
Definition: dgelliptic.hh:180
DiscreteFunctionSpaceType::GridPartType GridPartType
Definition: dgelliptic.hh:185
DomainDiscreteFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: dgelliptic.hh:164
LocalFunctionType::JacobianRangeType JacobianRangeType
Definition: dgelliptic.hh:177
DGEllipticOperator< typename JacobianOperator::DomainFunctionType, Model, Penalty > BaseType
Definition: dgelliptic.hh:154
LocalFunctionType::RangeType RangeType
Definition: dgelliptic.hh:175
const ModelType & model() const
Definition: dgelliptic.hh:137
DifferentiableDGEllipticOperator(const DiscreteFunctionSpaceType &space, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
constructor
Definition: dgelliptic.hh:198
RangeDiscreteFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: dgelliptic.hh:168
GridPartType::IntersectionIteratorType IntersectionIteratorType
Definition: dgelliptic.hh:186
DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: dgelliptic.hh:173
RangeDiscreteFunctionType::LocalFunctionType RangeLocalFunctionType
Definition: dgelliptic.hh:169
BaseType::DiscreteFunctionType DiscreteFunctionType
Definition: dgelliptic.hh:158
BaseType::DomainDiscreteFunctionType DomainDiscreteFunctionType
Definition: dgelliptic.hh:160
Dune::Fem::CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: dgelliptic.hh:191
RangeLocalFunctionType::RangeType RangeRangeType
Definition: dgelliptic.hh:170
void jacobian(const DomainDiscreteFunctionType &u, JacobianOperatorType &jOp) const
method to setup the jacobian of the operator for storage in a matrix
Definition: dgelliptic.hh:211
IntersectionIteratorType::Intersection IntersectionType
Definition: dgelliptic.hh:187
JacobianOperator JacobianOperatorType
Definition: dgelliptic.hh:156
BaseType::RangeDiscreteFunctionType RangeDiscreteFunctionType
Definition: dgelliptic.hh:161
LocalFunctionType::RangeFieldType RangeFieldType
Definition: dgelliptic.hh:176
DomainLocalFunctionType::RangeType DomainRangeType
Definition: dgelliptic.hh:166
EntityType::Geometry GeometryType
Definition: dgelliptic.hh:181