dune-fem  2.8-git
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 
22 namespace 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 
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
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
Definition: bindguard.hh:11
double pow(const Double &a, const double b)
Definition: double.hh:881
Double abs(const Double &a)
Definition: double.hh:871
DofStoragePolicy
Definition: dofstorage.hh:16
static constexpr std::decay_t< T > sum(T a)
Definition: utility.hh:33
Lagrange discrete function space.
Definition: lagrange/space.hh:131
BaseType::DomainFieldType DomainFieldType
Definition: automaticdifferenceoperator.hh:93
BaseType::RangeFieldType RangeFieldType
Definition: automaticdifferenceoperator.hh:92
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
HdivProjection(const DiscreteFunctionSpaceType &space)
constructor taking space
Definition: hdivprojection.hh:183
const DiscreteFunctionSpaceType & space() const
return reference to space
Definition: hdivprojection.hh:200
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:18
Definition: combinedspace/combinedspace.hh:32
A vector valued function space.
Definition: functionspace.hh:60
Definition: discontinuousgalerkin/legendre.hh:142