dune-fem 2.8.0
Loading...
Searching...
No Matches
hdivprojection.hh
Go to the documentation of this file.
1#warning "dune/fem/operator/projection/hdivprojection.hh is deprecated and will be moved to dune-fem-dg."
2
3#ifndef DUNE_FEM_HDIV_PROJECTION_HH
4#define DUNE_FEM_HDIV_PROJECTION_HH
5
6//- Dune includes
7#include <dune/common/deprecated.hh>
8#include <dune/geometry/referenceelements.hh>
9
10//- Dune-fem includes
17
18// make sure higher order Lagrange works (define USE_TWISTFREE_MAPPER)
21
22namespace Dune
23{
24
25 // External Forward Declarations
26 // -----------------------------
27
28 template< int dim, class CoordCont >
29 class YaspGrid;
30
31#ifdef ENABLE_UG
32 template< int dim >
33 class UGGrid;
34#endif // #ifdef ENABLE_UG
35
36 namespace Fem
37 {
38
60 template <class DiscreteFunctionType>
61 class HdivProjection : public SpaceOperatorInterface<DiscreteFunctionType>
62 {
63 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
64 typedef typename DiscreteFunctionType :: LocalFunctionType LocalFunctionType;
65 typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
66 typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
67 typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
68 typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
69 typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
70 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
71
72 typedef CachingQuadrature <GridPartType , 0> VolumeQuadratureType;
73 //typedef ElementQuadrature <GridPartType , 0> VolumeQuadratureType;
74
75 // face quadrature type
76 //typedef CachingQuadrature<GridPartType, 1> FaceQuadratureType;
78
79 typedef typename GridPartType :: GridType GridType;
80
81 enum { dimFaceRange = 1 };
82 enum { dimDomain = DiscreteFunctionSpaceType::dimDomain };
83 enum { dimFaceDomain = dimDomain - 1 };
84 enum { polOrdN = DiscreteFunctionSpaceType :: polynomialOrder };
85
88
89 enum { gradPolOrd = ((polOrdN - 1) < 0) ? 0 : (polOrdN - 1) };
90
91 template <class Space> struct BubbleM;
92
93 template <class FunctionSpaceImp,
94 class GridPartImp,
95 int polOrd,
96 template <class> class StorageImp,
97 template <class,class,int,template <class> class> class DiscreteFunctionSpaceImp>
98 struct BubbleM <DiscreteFunctionSpaceImp<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> >
99 {
100 enum { bubblePModifier = dimDomain - 1 };
101 };
102
103 template <class FunctionSpaceImp,
104 class GridPartImp,
105 int polOrd,
106 template <class> class StorageImp>
107 struct BubbleM < LegendreDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> >
108 {
109 enum { bubblePModifier = 1 };
110 };
111
112 template <class DiscreteFunctionSpaceImp,
113 int N,
114 DofStoragePolicy policy>
115 struct BubbleM <CombinedSpace<DiscreteFunctionSpaceImp,N,policy> >
116 {
117 enum { bubblePModifier = BubbleM< DiscreteFunctionSpaceImp > :: bubblePModifier };
118 };
119
120 typedef BubbleM <DiscreteFunctionSpaceType> BubbleMType;
121
122#ifdef USE_TWISTFREE_MAPPER
123 // modifier for bubble pol order
124 enum { bubblePModifier = BubbleMType :: bubblePModifier };
125
126 // we need polOrd + 1 in 2d and polOrd + 2 in 3d
127 enum { bubblePolOrd = (polOrdN + bubblePModifier) };
128#else
129#ifndef NDEBUG
130#warning "Hdiv-Projection only working for polOrd = 1 (enable higher order Lagrange with -DUSE_TWISTFREE_MAPPER)"
131#endif
132 // modifier for bubble pol order
133 enum { bubblePModifier = 1 };
134
135 // limit bubble polord otherwise no compilation possible
136 enum { bubblePolOrd = (polOrdN > 1) ? 2 : (polOrdN + 1) };
137#endif
138
139 template <class Space> struct Spaces;
140
141 template <class FunctionSpaceImp,
142 class GridPartImp,
143 int polOrd,
144 template <class> class StorageImp,
145 template <class,class,int,template <class> class> class DiscreteFunctionSpaceImp>
146 struct Spaces<DiscreteFunctionSpaceImp<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> >
147 {
148 //typedef LagrangeDiscreteFunctionSpace<ElSpaceType,GridPartImp,gradPolOrd,StorageImp> ElementGradientSpaceType;
149 //typedef LegendreDiscontinuousGalerkinSpace<ElSpaceType,GridPartImp,gradPolOrd,StorageImp> ElementGradientSpaceType;
150 typedef DiscreteFunctionSpaceImp<ElSpaceType,GridPartImp,gradPolOrd,StorageImp> ElementGradientSpaceType;
151 //typedef DiscontinuousGalerkinSpace<ElSpaceType,GridPartImp,gradPolOrd,StorageImp> ElementGradientSpaceType;
152 typedef DiscreteFunctionSpaceImp<FaceSpaceType,GridPartType,polOrdN,StorageImp> FaceDiscreteSpaceType;
153 //typedef LegendreDiscontinuousGalerkinSpace<FaceSpaceType,GridPartType,polOrdN,StorageImp> FaceDiscreteSpaceType;
154 //typedef DiscontinuousGalerkinSpace<FaceSpaceType,GridPartType,polOrdN,StorageImp> FaceDiscreteSpaceType;
156 };
157
158 template <class DiscreteFunctionSpaceImp,
159 int N,
160 DofStoragePolicy policy>
161 struct Spaces<CombinedSpace<DiscreteFunctionSpaceImp,N,policy> >
162 {
163 typedef typename DiscreteFunctionSpaceImp :: GridPartType GridPartImp;
164 typedef Spaces<DiscreteFunctionSpaceImp> AllSpacesType;
165 typedef typename AllSpacesType :: ElementGradientSpaceType ElementGradientSpaceType;
166 typedef typename AllSpacesType :: ElementDiscreteSpaceType ElementDiscreteSpaceType;
167 typedef typename AllSpacesType :: FaceDiscreteSpaceType FaceDiscreteSpaceType;
168 };
169
170 typedef Spaces<DiscreteFunctionSpaceType> AllSpacesType;
171 typedef typename AllSpacesType :: ElementGradientSpaceType ElementGradientSpaceType;
172 typedef typename AllSpacesType :: ElementDiscreteSpaceType ElementDiscreteSpaceType;
173 typedef typename AllSpacesType :: FaceDiscreteSpaceType FaceDiscreteSpaceType;
174
175 const DiscreteFunctionSpaceType& space_;
176 GridPartType & gridPart_;
177 const FaceDiscreteSpaceType faceSpace_;
178 const ElementDiscreteSpaceType elSpace_;
179 const ElementGradientSpaceType gradSpace_;
180
181 public:
183 HdivProjection(const DiscreteFunctionSpaceType& space) DUNE_DEPRECATED_MSG( "Use dune-fem-dg implementation." ) :
184 space_(space),
185 gridPart_(space.gridPart()),
186 faceSpace_( gridPart_ ),
187 elSpace_( gridPart_ ),
188 gradSpace_( gridPart_ )
189 {}
190
191 HdivProjection(const HdivProjection& org) DUNE_DEPRECATED_MSG( "Use dune-fem-dg implementation." ) :
192 space_(org.space_),
193 gridPart_( org.gridPart_),
194 faceSpace_( gridPart_ ),
195 elSpace_( gridPart_ ),
196 gradSpace_( gridPart_ )
197 {}
198
200 const DiscreteFunctionSpaceType& space() const
201 {
202 return space_;
203 }
204 void setTime(double) {
205 }
206 double timeStepEstimate() const {
207 return 0.;
208 }
209
211 static double normalJump(const DiscreteFunctionType &discFunc, const int polyOrder = -1 )
212 {
213 typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
214 typedef typename IntersectionIteratorType :: Intersection IntersectionType;
215 typedef typename GridType :: template Codim<0> :: Entity EntityType;
216 typedef typename GridType :: Traits :: LocalIdSet LocalIdSetType;
217
218 enum { dim = GridType::dimension };
219
220 const DiscreteFunctionSpaceType &space = discFunc.space();
221 const GridPartType &gridPart = space.gridPart();
222 const int polOrd = (polyOrder <0) ? (2 * space.order() + 2) : polyOrder;
223
224 typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
225
226 RangeType ret (0.0);
227 RangeType neighRet (0.0);
228
229 // define type of quadrature
230 typedef ElementQuadrature <GridPartType , 1> FaceQuadratureType;
231
232 double sum = 0.0;
233
234 const LocalIdSetType &idSet = gridPart.grid().localIdSet();
235
236 for(const auto & en : space)
237 {
238 const LocalFuncType lf = discFunc.localFunction(en);
239
240 double localValue = 0.0;
241
242 IntersectionIteratorType endnit = gridPart.iend(en);
243 for(IntersectionIteratorType nit = gridPart.ibegin(en);
244 nit != endnit; ++nit)
245 {
246 const IntersectionType& inter=*nit;
247 // only interior faces are considered
248 if(inter.neighbor() )
249 {
250 EntityType nb = inter.outside();
251
252 if(idSet.id( en ) < idSet.id( nb ))
253 {
254 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
255 FaceQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
256
257 const LocalFuncType neighLf = discFunc.localFunction(nb);
258
259 const int quadNop = faceQuadInner.nop();
260 for (int l = 0; l < quadNop ; ++l)
261 {
262 DomainType normal =
263 inter.unitOuterNormal(faceQuadInner.localPoint(l));
264
265 lf.evaluate(faceQuadInner[l], ret);
266 neighLf.evaluate(faceQuadOuter[l], neighRet);
267
268 ret -= neighRet;
269
270 double val = ret * normal;
271 val *= faceQuadInner.weight(l);
272
273 localValue += std::abs(val);
274 }
275 }
276 }
277 } // end of intersection iterator
278 sum += std::abs(localValue);
279 }
280
281 return sum;
282 }
283
284 private:
285
286 // only works for 2d right now
287 void curl(const DomainType & arg, DomainType & dest, const int d ) const
288 {
289 if( DomainType :: dimension == 2 )
290 {
291 dest[0] = arg[1];
292 dest[1] = -arg[0];
293
294 return ;
295 }
296 else if( DomainType :: dimension == 3 )
297 {
298 if( d == 0 )
299 {
300 dest[0] = arg[1];
301 dest[1] = -arg[2];
302 dest[2] = 0;
303 return ;
304 }
305 else if ( d == 1 )
306 {
307 dest[0] = 0;
308 dest[1] = arg[2];
309 dest[2] = -arg[0];
310 return ;
311 }
312 else
313 {
314 dest[0] = -arg[2];
315 dest[1] = 0;
316 dest[2] = arg[1];
317 return ;
318 }
319 }
320 else
321 {
322 assert( false );
323 abort();
324 }
325 }
326
327 template <class EntityType,
328 class LocalFunctionType,
329 class ArrayType,
330 class MatrixType,
331 class VectorType>
332 void bubblePart(const ElementDiscreteSpaceType& space,
333 EntityType & en,
334 const LocalFunctionType & uLF, const int startRow ,
335 ArrayType& uRets,
336 MatrixType & matrix, VectorType& rhs ) const
337 {
338
339 typedef typename ElementDiscreteSpaceType :: BaseFunctionSetType BaseFunctionSetType;
340 typedef typename ElementDiscreteSpaceType :: LagrangePointSetType LagrangePointSetType;
341
342 typedef typename ElementDiscreteSpaceType :: JacobianRangeType JacobianRangeType;
343 typedef typename ElementDiscreteSpaceType :: DomainType DomainType;
344
345 enum { dim = GridType::dimension };
346
347 const LagrangePointSetType& lagrangePointSet = space.lagrangePointSet( en );
348 const BaseFunctionSetType bSet = space.baseFunctionSet( en );
349 const int polOrd = 2 * space.order(); // + 2;
350
351 const int cols = uLF.numDofs();
352 assert( uRets.size() == (unsigned int)cols );
353
354 VolumeQuadratureType quad (en,polOrd);
355 DomainType result;
356 std::vector< JacobianRangeType > valTmpVec( bSet.size() );
357 DomainType bVal;
358 DomainType aVal;
359
360 // get geometry
361 typedef typename EntityType :: Geometry Geometry ;
362 const Geometry& geo = en.geometry();
363 // get geometry type
364 const GeometryType& type = geo.type();
365 const int bubbleOffset = (type.isSimplex()) ? 0 : baseFunctionOffset( 0 );
366
367 // type of jacobian inverse
368 enum { cdim = Geometry :: coorddimension };
369 enum { mydim = Geometry :: mydimension };
370 typedef typename Geometry::JacobianInverseTransposed JacobianInverseType;
371
372 // get number of dofs for codim 0 (skip first for)
373 const int enDofs = numberOfBubbles( lagrangePointSet.numDofs( 0 ), type ,
374 cols, startRow );
375 const int bubbleMod = bubbleModifier( mydim );
376
377 const int quadNop = quad.nop();
378 for (int l = 0; l < quadNop ; ++l)
379 {
380 // get jacobian inverse
381 const JacobianInverseType& inv = geo.jacobianInverseTransposed( quad.point(l) );
382
383 // get integration element
384 const double intel = quad.weight(l) *
385 geo.integrationElement(quad.point(l));
386
387 // evaluate u
388 uLF.evaluate(quad[l], result);
389
390 // evaluate base functions of u
391 uLF.baseFunctionSet().evaluateAll( quad[l], uRets );
392
393 // evaluate gradients
394 bSet.jacobianAll( quad[l], inv, valTmpVec );
395
396 // for all bubble functions
397 for( int i = 0 ; i<enDofs; i += bubbleMod )
398 {
399 // we might have other row
400 int row = startRow + i;
401
402 // map to lagrange base function number
403 const int localBaseFct = ((int) i/bubbleMod) + bubbleOffset;
404 // get dof number of 'localBaseFct' dof on codim 0 subentity 0
405 const int baseFct = lagrangePointSet.entityDofNumber( 0, 0, localBaseFct );
406
407 // evaluate gradient
408 JacobianRangeType& val = valTmpVec[ baseFct ];
409
410 //apply inverse
411 //inv.mv( valTmp[0], val[0] );
412
413 for(int d = 0; d<bubbleMod; ++d )
414 {
415 // apply curl
416 curl(val[0], aVal, d );
417
418 double r = aVal * result;
419 r *= intel;
420 rhs[row] += r;
421
422 // for cols make matrix
423 for(int j=0; j<cols; ++j)
424 {
425 double t = aVal * uRets[ j ];
426 t *= intel;
427 matrix[ row ][ j ] += t;
428 }
429
430 // increase row
431 ++row;
432 }
433 }
434 }
435 }
436
437 template <class EntityType, class LocalFunctionType, class ArrayType,
438 class MatrixType, class VectorType>
439 void gradientPart(const ElementGradientSpaceType & space,
440 EntityType & en,
441 const LocalFunctionType & uLF, const int startRow ,
442 ArrayType& uRets, MatrixType& matrix, VectorType& rhs ) const
443 {
444 typedef typename ElementGradientSpaceType :: BaseFunctionSetType BaseFunctionSetType;
445 const BaseFunctionSetType bSet = space.baseFunctionSet( en );
446 int polOrd = 2 * space.order() + 1;
447
448 const int localRows = gradientBaseFct( bSet );
449 const int cols = uLF.numDofs();
450
451 VolumeQuadratureType quad (en,polOrd);
452
453 RangeType result;
454 RangeType uPhi;
455
456 typedef typename ElementGradientSpaceType :: JacobianRangeType GradJacobianRangeType;
457 std::vector< GradJacobianRangeType > gradTmpVec( bSet.size() );
458 GradJacobianRangeType gradPhi;
459
460 typedef typename EntityType :: Geometry Geometry ;
461 const Geometry& geo = en.geometry();
462
463 enum { cdim = Geometry :: coorddimension };
464 enum { mydim = Geometry :: mydimension };
465 typedef typename Geometry::JacobianInverseTransposed JacobianInverseType;
466
467 const int quadNop = quad.nop();
468 for (int l = 0; l < quadNop ; ++l)
469 {
470 // get jacobian inverse
471 const JacobianInverseType& inv = geo.jacobianInverseTransposed( quad.point( l ) );
472
473 // get integration element
474 const double intel = quad.weight(l) *
475 geo.integrationElement( quad.point( l ) );
476
477 // evaluate uLF
478 uLF.evaluate(quad[l], result);
479
480 // evaluate base function on quadrature point
481 uLF.baseFunctionSet().evaluateAll( quad[l], uRets );
482
483 // evaluate gradient (skip first function because this function
484 // is constant and the gradient therefore 0 )
485 bSet.jacobianAll( quad[l], inv, gradTmpVec );
486
487 // apply jacobian Inverse
488 for(int i=0; i<localRows; ++i)
489 {
490 // we might have other row
491 const int row = startRow + i;
492
493 // evaluate gradient (skip first function because this function
494 // is constant and the gradient therefore 0 )
495 GradJacobianRangeType& gradPhi = gradTmpVec[ baseFunctionOffset( i ) ];
496
497 const double uDGVal = result * gradPhi[0];
498 rhs[row] += uDGVal * intel;
499
500 // for cols make matrix
501 for(int j=0; j<cols; ++j)
502 {
503 //uLF.baseFunctionSet().evaluate(j, quad[l], uPhi);
504 //const double val = uPhi * gradPhi[0];
505 const double val = uRets[ j ] * gradPhi[0];
506 matrix[row][j] += val * intel;
507 }
508 }
509 }
510 }
511
512 template <class FaceBSetType, class GridType>
513 struct GetSubBaseFunctionSet
514 {
515 template <class EntityType, class SpaceType>
516 static inline FaceBSetType faceBaseSet(const EntityType& en, const SpaceType& space)
517 {
518 return space.baseFunctionSet( (en.template subEntity<1> (0) )->type() );
519 }
520 };
521
522 template <class FaceBSetType, int dim, class CoordCont>
523 struct GetSubBaseFunctionSet< FaceBSetType, YaspGrid< dim, CoordCont > >
524 {
525 template <class EntityType, class SpaceType>
526 static inline FaceBSetType faceBaseSet(const EntityType& en, const SpaceType& space)
527 {
528 return space.baseFunctionSet( GeometryTypes::cube(dim-1) );
529 }
530 };
531
532#ifdef ENABLE_UG
533 template <class FaceBSetType, int dim>
534 struct GetSubBaseFunctionSet<FaceBSetType, UGGrid<dim> >
535 {
536 template <class EntityType, class SpaceType>
537 static inline FaceBSetType faceBaseSet(const EntityType& en, const SpaceType& space)
538 {
539 const GeometryType geoType (en.geometry().type().basicType(),dim-1);
540 return space.baseFunctionSet( geoType );
541 }
542 };
543#endif
544
545 enum { gradFuncOffset = 1 };
546 template <class GradBaseFunctionSet>
547 int gradientBaseFct(const GradBaseFunctionSet& gradSet) const
548 {
549 return (gradPolOrd <= 0) ? 0 : gradSet.size() - gradFuncOffset;
550 }
551
552 int baseFunctionOffset(const int i) const
553 {
554 return i + gradFuncOffset;
555 }
556
557 int numberOfBubbles( const int bubbles , const GeometryType& type,
558 const int allDofs, const int allOther ) const
559 {
560 /*
561 // for hexahedrons this is different
562 if( type.isHexahedron() )
563 {
564 return (bubbleModifier( type.dim() )) * (bubbles - gradFuncOffset) - 1;
565 }
566 else
567 */
568 {
569 // the rest is padded with bubble functions
570 const int numBubble = allDofs - allOther ;
571 return (numBubble > 0) ? numBubble : 0;
572 }
573 }
574
575 int bubbleModifier( const int dim ) const
576 {
577 // return 1 for 2d and 3 for 3d
578 return (dim - 2) * (dim - 1) + 1;
579 }
580
582 void project(const DiscreteFunctionType &uDG,
583 DiscreteFunctionType & velo ) const
584 {
585 typedef typename DiscreteFunctionType::Traits::DiscreteFunctionSpaceType FunctionSpaceType;
586 enum { localBlockSize = FunctionSpaceType::localBlockSize };
587
588 typedef typename FunctionSpaceType::GridType GridType;
589 typedef typename FunctionSpaceType::GridPartType GridPartType;
590 typedef typename FunctionSpaceType::IteratorType Iterator;
591
592 typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
593 typedef typename FunctionSpaceType::RangeType RangeType;
594
595 enum { dim = GridType::dimension };
596 typedef typename GridType :: ctype coordType;
597
598 const FunctionSpaceType& space = uDG.space();
599
600 // for polOrd 0 this is not working
601 // then just copy the function
602 if(space.order() < 1 )
603 {
604 velo.assign( uDG );
605 return ;
606 }
607
608 const int polOrd = 2 * space.order() + 2;
609
610 typedef typename FaceDiscreteSpaceType :: BaseFunctionSetType FaceBSetType ;
611
612 typedef typename ElementGradientSpaceType :: BaseFunctionSetType GradientBaseSetType ;
613
614 typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
615
616 Iterator start = space.begin();
617 // if empty grid, do nothing
618 if( start == space.end() ) return;
619
620 // only working for spaces with one element type
621 if( space.multipleGeometryTypes() )
622 {
623 if( space.indexSet().geomTypes(0).size() > 1)
624 {
625 assert( space.indexSet().geomTypes(0).size() == 1 );
626 DUNE_THROW(NotImplemented,"H-div projection not implemented for hybrid grids");
627 }
628 }
629
630 const GeometryType startType = start->type();
631
632 // only working for spaces with one element type
633 if( startType.isHexahedron() && space.order() > 1 )
634 {
635 assert( ! startType.isHexahedron() || space.order() <= 1 );
636 DUNE_THROW(NotImplemented,"H-div projection not implemented for p > 1 on hexas! ");
637 }
638
639 const int desiredOrder = space.order() + bubblePModifier;
640 // check Lagrange space present
641 if( bubblePolOrd != desiredOrder )
642 {
643 assert( bubblePolOrd == desiredOrder );
644 DUNE_THROW(NotImplemented,"H-div projection not working for "
645 << space.order() << " when LagrangeSpace of order "<< desiredOrder << " is missing");
646 }
647
648 // colums are dofs of searched function
649 LocalFuncType lf = uDG.localFunction(*start);
650 const int numDofs = lf.numDofs();
651 //std::cout << numDofs << " numDofs \n";
652
653 const FaceBSetType faceSet =
654 GetSubBaseFunctionSet<FaceBSetType,GridType>::faceBaseSet( *start , faceSpace_ );
655 // number of dofs on faces
656 const int numFaceDofs = faceSet.size();
657
658 const GradientBaseSetType gradSet = gradSpace_.baseFunctionSet(*start);
659 // in case of linear space the is zero
660 const int numGradDofs = gradientBaseFct( gradSet );
661
662 const Dune::ReferenceElement< coordType, dim > & refElem =
663 Dune::ReferenceElements< coordType, dim >::general( startType );
664
665 // get number of faces
666 const int overallFaceDofs = numFaceDofs * refElem.size(1);
667
668 // get all element dofs from Lagrange space
669 const int numBubbleDofs = (space.order() <= 1) ? 0 :
670 numberOfBubbles( elSpace_.lagrangePointSet( *start ).numDofs ( 0 ) , startType,
671 numDofs, overallFaceDofs + numGradDofs );
672
673 const int rows = (overallFaceDofs + numGradDofs + numBubbleDofs);
674
675 // number of columns
676 const int cols = numDofs;
677
678 DynamicArray< RangeFieldType > rets(numDofs);
679 DynamicArray< RangeType > uRets(numDofs);
680
681 typedef FieldMatrix<RangeFieldType,localBlockSize,localBlockSize> FieldMatrixType;
682
683 // resulting system matrix
684 FieldMatrixType inv;
685
686 typedef FieldVector<RangeFieldType,localBlockSize> VectorType;
687 VectorType fRhs(0.0);
688
689 assert( numDofs == localBlockSize );
690 if( numDofs != localBlockSize )
691 {
692 DUNE_THROW(InvalidStateException,"wrong sizes ");
693 }
694
695 // flag to say whether we have non-symetric or symetric situation
696 const bool nonSymetric = (cols != rows);
697
698 // matrix type
699 typedef Fem::DenseMatrix<double> MatrixType;
700 MatrixType matrix(rows,cols);
701 typedef typename MatrixType :: RowType RowType;
702 // matrix type
703 MatrixType fakeMatrix(cols,cols);
704 RowType rhs(rows,0.0);
705 RowType fakeRhs(numDofs,0.0);
706
707 // iterate over grid
708 for(const auto & en : space)
709 {
710 if( nonSymetric )
711 {
712 // reset values
713 matrix = 0.0;
714
715 // reset rhs
716 for(int i=0; i<rows; ++i)
717 {
718 rhs[i] = 0.0;
719 }
720
721 // fill non-symetric matrix
722 fillMatrix(gridPart_,en,uDG,
723 faceSpace_,
724 gradSpace_, overallFaceDofs,
725 elSpace_, rows - numBubbleDofs,
726 polOrd,numDofs,numFaceDofs,
727 rets,uRets, matrix,rhs);
728
729 // apply least square
730 matrix.multTransposed(rhs, fakeRhs);
731 fakeMatrix.multiply_AT_A(matrix);
732
733 // copy values
734 for(int i=0; i<numDofs; ++i)
735 {
736 fRhs[i] = fakeRhs[i];
737 for(int j=0; j<numDofs; ++j)
738 {
739 inv[i][j] = fakeMatrix[i][j];
740 }
741 }
742 }
743 else
744 {
745 // reset values
746 inv = 0.0;
747 // reset rhs
748 fRhs = 0.0;
749
750 assert( cols == rows );
751 // fill inv and fRhs directly
752 fillMatrix(gridPart_,en,uDG,
753 faceSpace_,
754 gradSpace_, rows - numBubbleDofs - numGradDofs,
755 elSpace_, rows - numBubbleDofs,
756 polOrd,numDofs,numFaceDofs,
757 rets, uRets, inv,fRhs);
758 }
759
760 // set new values to new velocity function
761 {
762 LocalFuncType veloLF = velo.localFunction( en );
763
764 // solve linear system
765 luSolve( inv, fRhs );
766 const VectorType& x = fRhs ;
767
768 for(int i=0; i<localBlockSize; ++i)
769 {
770 veloLF[ i ] = x[ i ];
771 }
772 }
773 }
774 }
775
776 template <class GridPartType,
777 class EntityType,
778 class ArrayType,
779 class Array2Type,
780 class MatrixType,
781 class VectorType>
782 void fillMatrix(const GridPartType& gridPart,
783 const EntityType& en,
784 const DiscreteFunctionType& uDG,
785 const FaceDiscreteSpaceType& faceSpace,
786 const ElementGradientSpaceType& gradSpace, const int startGradDofs,
787 const ElementDiscreteSpaceType& elSpace, const int startBubbleDofs,
788 const int polOrd, const int numDofs, const int numFaceDofs,
789 ArrayType& rets, Array2Type& uRets,
790 MatrixType& matrix, VectorType& rhs) const
791 {
792 typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
793 typedef typename IntersectionIteratorType::Intersection IntersectionType;
794
795 typedef typename FaceDiscreteSpaceType :: BaseFunctionSetType FaceBSetType ;
796 typedef typename FaceDiscreteSpaceType :: RangeType FaceRangeType;
797 FaceRangeType faceVal;
798
799 typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
800 RangeType ret (0.0);
801 RangeType neighRet (0.0);
802 RangeType uPhi (0.0);
803
804 typedef typename DiscreteFunctionType :: LocalFunctionType LocalFuncType ;
805 typedef typename DiscreteFunctionType :: DiscreteFunctionSpaceType
806 :: BaseFunctionSetType BaseFunctionSetType;
807
808 // get uDg local on entity
809 const LocalFuncType uLF = uDG.localFunction(en);
810
811 // get base functions set
812 const BaseFunctionSetType & bSet = uLF.baseFunctionSet();
813
814 // iterate over intersections
815 IntersectionIteratorType endnit = gridPart.iend(en);
816 for(IntersectionIteratorType nit = gridPart.ibegin(en);
817 nit != endnit; ++nit)
818 {
819 const IntersectionType& inter=*nit;
820 // get base function set of face
821 const FaceBSetType &faceSet = faceSpace.baseFunctionSet( inter.type() );
822
823 const int firstRow = inter.indexInInside() * numFaceDofs;
824
825 // only interior faces are considered
826 if(inter.neighbor())
827 {
828 // get neighbor entity
829 EntityType nb = inter.outside();
830
831 // get local function of neighbor
832 const LocalFuncType uNeighLf = uDG.localFunction(nb);
833
834 //typedef TwistUtility<GridType> TwistUtilityType;
835 // for conforming situations apply Quadrature given
836 //if( TwistUtilityType::conforming(gridPart.grid(),inter) )
837 if( inter.conforming() )
838 {
839 // create quadratures
840 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
841 FaceQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
842
843 applyLocalNeighbor(nit,faceQuadInner,faceQuadOuter,
844 bSet,faceSet, uLF, uNeighLf,
845 firstRow, numFaceDofs,
846 rets,
847 ret,neighRet,faceVal,
848 matrix, rhs);
849 }
850 else
851 {
852 // type of quadrature for non-conforming intersections
853 typedef typename FaceQuadratureType ::
854 NonConformingQuadratureType NonConformingQuadratureType;
855 // create quadratures
856 NonConformingQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
857 NonConformingQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
858
859 applyLocalNeighbor(nit,faceQuadInner,faceQuadOuter,
860 bSet,faceSet, uLF, uNeighLf,
861 firstRow, numFaceDofs,
862 rets,
863 ret,neighRet,faceVal,
864 matrix, rhs);
865
866 }
867 }
868
869 // only interior faces are considered
870 if(inter.boundary())
871 {
872 // create quadrature
873 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
874 const int quadNop = faceQuadInner.nop();
875
876 std::vector< RangeType > uPhiVec( numDofs );
877 std::vector< FaceRangeType > faceValVec( numFaceDofs );
878
879 for (int l = 0; l < quadNop ; ++l)
880 {
881 DomainType unitNormal =
882 inter.integrationOuterNormal(faceQuadInner.localPoint(l));
883
884 const double faceVol = unitNormal.two_norm();
885 unitNormal *= 1.0/faceVol;
886
887 // get integration weight
888 const double intel = faceVol * faceQuadInner.weight(l);
889
890 // evaluate function
891 uLF.evaluate(faceQuadInner[l], ret);
892
893 RangeFieldType val = ret * unitNormal;
894 val *= intel;
895
896 bSet.evaluateAll( faceQuadInner[l], uPhiVec );
897
898 // evaluate base functions
899 for(int i=0; i<numDofs; ++i)
900 {
901 rets[i] = uPhiVec[ i ] * unitNormal;
902 rets[i] *= intel;
903 }
904
905 faceSet.evaluateAll( faceQuadInner[ l ], faceValVec );
906
907 int row = firstRow;
908 for(int j=0; j<numFaceDofs; ++j, ++row)
909 {
910 FaceRangeType& faceVal = faceValVec[ j ];
911 rhs[row] += val*faceVal[0];
912
913 for(int i=0; i<numDofs; ++i)
914 {
915 matrix[row][i] += (faceVal[0] * rets[i]);
916 }
917 }
918 }
919 }
920 }
921
922 // add gradient part
923 if( gradPolOrd > 0 )
924 {
925 gradientPart(gradSpace, en, uLF, startGradDofs, uRets, matrix, rhs );
926 }
927
928 // add bubble part
929 if( bubblePolOrd - bubblePModifier > 1 )
930 {
931 bubblePart(elSpace, en, uLF, startBubbleDofs, uRets, matrix, rhs);
932 }
933
934 // printMatrix( matrix );
935 }
936
937 template <class MatrixType>
938 void printMatrix(const MatrixType& matrix) const
939 {
940 std::cout << "Print Matrix \n";
941 for(size_t row = 0; row < matrix.N(); ++row)
942 {
943 std::cout << row << ": ";
944 for(size_t col = 0; col< matrix.M(); ++col)
945 {
946 if( std::abs( matrix[row][col] ) < 1e-12 )
947 std::cout << "0 ";
948 else
949 std::cout << matrix[row][col] << " ";
950 }
951 std::cout << std::endl;
952 }
953 std::cout << "Finished print Matrix \n";
954 }
955
956 template <class IntersectionIteratorType,
957 class QuadratureType,
958 class BaseFunctionSetType,
959 class FaceBaseFunctionSetType,
960 class LocalFunctionType,
961 class ArrayType,
962 class RangeType,
963 class FaceRangeType,
964 class MatrixType,
965 class RHSType>
966 static void applyLocalNeighbor(const IntersectionIteratorType& nit,
967 const QuadratureType& faceQuadInner,
968 const QuadratureType& faceQuadOuter,
969 const BaseFunctionSetType& bSet,
970 const FaceBaseFunctionSetType& faceSet,
971 const LocalFunctionType& uLF,
972 const LocalFunctionType& uNeighLf,
973 const int firstRow,
974 const int numFaceDofs,
975 ArrayType& rets,
976 RangeType& ret, RangeType& neighRet,
977 FaceRangeType& faceVal,
978 MatrixType& matrix,
979 RHSType& rhs)
980 {
981 const typename IntersectionIteratorType::Intersection& inter = *nit;
982 const int quadNop = faceQuadInner.nop();
983 const int numDofs = uLF.numDofs();
984
985 std::vector< RangeType > phiVec( numDofs );
986 std::vector< FaceRangeType > facePhiVec( numFaceDofs );
987
988 for (int l = 0; l < quadNop ; ++l)
989 {
990 DomainType unitNormal =
991 inter.integrationOuterNormal(faceQuadInner.localPoint(l));
992
993 // get unit outer normal
994 const double faceVol = unitNormal.two_norm();
995 unitNormal *= 1.0/faceVol;
996
997 // integration weight
998 const double intel = faceVol * faceQuadInner.weight(l);
999
1000 // evaluate dg velocity
1001 uLF.evaluate(faceQuadInner[l], ret);
1002 uNeighLf.evaluate(faceQuadOuter[l], neighRet);
1003
1004 // take mean value
1005 ret += neighRet;
1006 ret *= 0.5;
1007
1008 RangeFieldType val = ret * unitNormal;
1009 val *= intel;
1010
1011 bSet.evaluateAll( faceQuadInner[l], phiVec );
1012
1013 // evaluate base functions
1014 for(int i=0; i<numDofs; ++i)
1015 {
1016 rets[i] = phiVec[ i ] * unitNormal;
1017 rets[i] *= intel;
1018 }
1019
1020 // evaluate all basis functions
1021 faceSet.evaluateAll( faceQuadInner[ l ], facePhiVec );
1022
1023 int row = firstRow;
1024 for(int j=0; j<numFaceDofs; ++j, ++row )
1025 {
1026 FaceRangeType& faceVal = facePhiVec[ j ];
1027
1028 rhs[row] += val * faceVal[0];
1029
1030 for(int i=0; i<numDofs; ++i)
1031 {
1032 matrix[row][i] += (faceVal[0] * rets[i]);
1033 }
1034 }
1035 }
1036 }
1037
1038 public:
1040 virtual void operator () (const DiscreteFunctionType &arg,
1041 DiscreteFunctionType& dest) const
1042 {
1043 // apply H-div projection
1044 project(arg,dest);
1045 }
1046
1047 template <class AdaptationType>
1048 static void estimator(const DiscreteFunctionType &velo,
1049 AdaptationType& adaptation)
1050 {
1051 typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
1052 typedef typename GridType :: Traits :: LocalIdSet LocalIdSetType;
1053
1054 enum { dim = GridType :: dimension };
1055 typedef typename GridType :: template Codim<0> :: Entity EntityType;
1056
1057 const DiscreteFunctionSpaceType& space = velo.space();
1058 GridPartType & gridPart = space.gridPart();
1059 int polOrd = space.order() + 1;
1060
1061 // get local id set
1062 const LocalIdSetType& localIdSet = gridPart.grid().localIdSet();
1063
1064 // define type of face quadrature
1066
1067 for(const auto & en : space)
1068 {
1069 const LocalFuncType uLF = velo.localFunction(en);
1070 const double enVol = en.geometry().volume();
1071
1072 typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
1073 typedef typename IntersectionIteratorType :: Intersection IntersectionType;
1074 IntersectionIteratorType endnit = gridPart.iend(en);
1075 for(IntersectionIteratorType nit = gridPart.ibegin(en);
1076 nit != endnit; ++nit)
1077 {
1078 const IntersectionType& inter=*nit;
1079 double enError = 0.0;
1080 // only interior faces are considered
1081 if(inter.neighbor())
1082 {
1083 EntityType nb = inter.outside();
1084 const double enVol_nbVol = 0.5 * (enVol + nb.geometry().volume());
1085
1086#if HAVE_MPI
1087 // get partition type
1088 const bool interiorEntity = (nb.partitionType() == InteriorEntity);
1089#else
1090 const bool interiorEntity = true;
1091#endif
1092 if( (localIdSet.id( en ) < localIdSet.id( nb ))
1093#if HAVE_MPI
1094 || ( ! interiorEntity )
1095#endif
1096 )
1097 {
1098 const LocalFuncType uNeighLf = velo.localFunction(nb);
1099
1100 //typedef TwistUtility<GridType> TwistUtilityType;
1101 // for conforming situations apply Quadrature given
1102 //if( TwistUtilityType::conforming(gridPart.grid(),inter) )
1103 if( inter.conforming() )
1104 {
1105 // create quadratures
1106 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
1107 FaceQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
1108
1109 applyLocalNeighEstimator(nit,nb,faceQuadInner,faceQuadOuter,
1110 uLF, uNeighLf, enVol_nbVol, interiorEntity,
1111 enError, adaptation);
1112 }
1113 else
1114 {
1115 // type of quadrature for non-conforming intersections
1116 typedef typename FaceQuadratureType ::
1117 NonConformingQuadratureType NonConformingQuadratureType;
1118 // create quadratures
1119 NonConformingQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
1120 NonConformingQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
1121
1122 applyLocalNeighEstimator(nit,nb,faceQuadInner,faceQuadOuter,
1123 uLF, uNeighLf, enVol_nbVol, interiorEntity,
1124 enError, adaptation);
1125 }
1126
1127 } // end enId < nbId
1128 } // end neighbor
1129
1130 if(enError > 0.0)
1131 {
1132 adaptation.addToLocalIndicator( en , enError );
1133 }
1134 } // end intersection iterator
1135 }
1136 }
1137
1138 private:
1139 template <class IntersectionIteratorType,
1140 class EntityType,
1141 class QuadratureType,
1142 class LocalFunctionType,
1143 class AdaptationType>
1144 static inline void applyLocalNeighEstimator(const IntersectionIteratorType& nit,
1145 const EntityType& nb,
1146 const QuadratureType& faceQuadInner,
1147 const QuadratureType& faceQuadOuter,
1148 const LocalFunctionType& uLF,
1149 const LocalFunctionType& uNeighLf,
1150 const double enVol_nbVol,
1151 const bool interiorEntity,
1152 double& enError,
1153 AdaptationType& adaptation)
1154 {
1155 const typename IntersectionIteratorType::Intersection& inter=*nit;
1156 enum { dim = GridType :: dimension };
1157 RangeType jump;
1158 RangeType neighRet;
1159 double nbError = 0.0;
1160
1161 const int quadNop = faceQuadInner.nop();
1162 for (int l = 0; l < quadNop ; ++l)
1163 {
1164 DomainType unitNormal =
1165 inter.integrationOuterNormal(faceQuadInner.localPoint(l));
1166
1167 double faceVol = unitNormal.two_norm();
1168 unitNormal *= 1.0/faceVol;
1169
1170 // in case of power != 0
1171 if(dim > 2)
1172 {
1173 const double power = (1.0/(dim-1));
1174 faceVol = pow(faceVol, power);
1175 }
1176
1177 // evaluate | (u_l * n_l) + (u_r * n_r) | = | (u_l - u_r) * n_l |
1178 uLF.evaluate(faceQuadInner[l], jump);
1179 uNeighLf.evaluate(faceQuadOuter[l], neighRet);
1180
1181 // get difference
1182 jump -= neighRet;
1183
1184 const double weight = (enVol_nbVol) * faceQuadInner.weight(l);
1185
1186 double error = weight * SQR(jump * unitNormal);
1187 error = std::abs( error );
1188
1189 enError += error;
1190 nbError += error;
1191 }
1192
1193 if( (nbError > 0.0)
1194#if HAVE_MPI
1195 // only add neihgbor on interior entities
1196 && (interiorEntity)
1197#endif
1198 )
1199 {
1200 adaptation.addToLocalIndicator( nb , nbError );
1201 }
1202 }
1203
1204 // LU decomposition of matrix (matrix and b are overwritten)
1205 //
1206 // param[inout] a Matrix that LU decomposition is calculated for
1207 // param[in] b right hand side
1208 // param[out] solution solution of linear system
1209 template <class MatrixType, class VectorType>
1210 void luSolve(MatrixType& a,
1211 VectorType& x) const
1212 {
1213 typedef typename VectorType :: field_type ctype;
1214 enum { n = VectorType :: dimension };
1215
1216 // make sure we got the right dimensions
1217 assert( a.N() == a.M() );
1218 assert( a.N() == n );
1219
1220 // pivot storage
1221 int p[ n-1 ];
1222
1223 for(int i=0; i<n-1; ++i)
1224 {
1225 // initialize
1226 p[i] = i;
1227
1228 // Pivot search
1229 ctype max_abs = 0;
1230 for(int k=i; k<n; ++k)
1231 {
1232 if ( std::abs(a[k][i]) > max_abs )
1233 {
1234 max_abs = fabs(a[k][i]);
1235 p[i] = k;
1236 }
1237 }
1238
1239 if( p[ i ] != i )
1240 {
1241 // toggle row i with row argmax=p[i]
1242 for(int j=0; j<n; ++j)
1243 {
1244 const ctype tmp = a[ p[i] ][j];
1245 a[ p[i] ][j] = a[i][j];
1246 a[i][j] = tmp;
1247 }
1248 }
1249
1250 // elimination
1251 for(int k=i+1; k<n; ++k)
1252 {
1253 const ctype lambda = a[k][i] / a[i][i];
1254 a[k][i] = lambda;
1255 for(int j=i+1; j<n; ++j) a[k][j] -= lambda * a[i][j];
1256 }
1257 }
1258
1259 // 1. x = Px_old, permutation with right hand side
1260 for(int i=0; i<n-1; ++i)
1261 {
1262 const ctype tmp = x[ i ];
1263 x[ i ] = x[ p[ i ] ];
1264 x[ p[ i ] ] = tmp;
1265 }
1266
1267 // 1. Lx = x_old, forward loesen
1268 for(int i=0; i<n; ++i)
1269 {
1270 ctype dot = 0;
1271 for(int j=0; j<i; ++j) dot += a[i][j] * x[j];
1272 x[i] -= dot;
1273 }
1274
1275 // 2. Ux = x_old, backward solve
1276 for(int i=n-1; i>=0; --i)
1277 {
1278 ctype dot = 0;
1279 for(int j=i+1; j<n; ++j) dot += a[i][j] * x[j];
1280 x[i] = (x[i] - dot) / a[i][i];
1281 }
1282 }
1283 };
1284
1285 } // namespace Fem
1286
1287} // namespace Dune
1288#endif // #ifndef DUNE_FEM_HDIV_PROJECTION_HH
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
double pow(const Double &a, const double b)
Definition: double.hh:881
DofStoragePolicy
Definition: dofstorage.hh:16
Lagrange discrete function space.
Definition: lagrange/space.hh:131
interface for time evolution operators
Definition: spaceoperatorif.hh:38
DiscreteFunctionSpaceType SpaceType
convenience typedef for space type
Definition: spaceoperatorif.hh:49
H-div Projection for discontinuous discrete functions. The projection is described in detail in:
Definition: hdivprojection.hh:62
void setTime(double)
set time for operators
Definition: hdivprojection.hh:204
double timeStepEstimate() const
estimate maximum time step
Definition: hdivprojection.hh:206
static double normalJump(const DiscreteFunctionType &discFunc, const int polyOrder=-1)
return sum of jumps of discrete function normal to intersection
Definition: hdivprojection.hh:211
HdivProjection(const HdivProjection &org)
Definition: hdivprojection.hh:191
static void estimator(const DiscreteFunctionType &velo, AdaptationType &adaptation)
Definition: hdivprojection.hh:1048
const DiscreteFunctionSpaceType & space() const
return reference to space
Definition: hdivprojection.hh:200
HdivProjection(const DiscreteFunctionSpaceType &space)
constructor taking space
Definition: hdivprojection.hh:183
virtual void operator()(const DiscreteFunctionType &arg, DiscreteFunctionType &dest) const
application operator projection arg to H-div space
Definition: hdivprojection.hh:1040
quadrature on the codim-0 reference element
Definition: elementquadrature.hh:58
Definition: combinedspace/combinedspace.hh:32
A vector valued function space.
Definition: functionspace.hh:60
Definition: discontinuousgalerkin/legendre.hh:142