1#warning "dune/fem/operator/projection/hdivprojection.hh is deprecated and will be moved to dune-fem-dg."
3#ifndef DUNE_FEM_HDIV_PROJECTION_HH
4#define DUNE_FEM_HDIV_PROJECTION_HH
7#include <dune/common/deprecated.hh>
8#include <dune/geometry/referenceelements.hh>
28 template<
int dim,
class CoordCont >
60 template <
class DiscreteFunctionType>
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;
79 typedef typename GridPartType :: GridType GridType;
81 enum { dimFaceRange = 1 };
82 enum { dimDomain = DiscreteFunctionSpaceType::dimDomain };
83 enum { dimFaceDomain = dimDomain - 1 };
84 enum { polOrdN = DiscreteFunctionSpaceType :: polynomialOrder };
89 enum { gradPolOrd = ((polOrdN - 1) < 0) ? 0 : (polOrdN - 1) };
91 template <
class Space>
struct BubbleM;
93 template <
class FunctionSpaceImp,
96 template <
class>
class StorageImp,
97 template <
class,
class,
int,
template <
class>
class>
class DiscreteFunctionSpaceImp>
98 struct BubbleM <DiscreteFunctionSpaceImp<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> >
100 enum { bubblePModifier = dimDomain - 1 };
103 template <
class FunctionSpaceImp,
106 template <
class>
class StorageImp>
109 enum { bubblePModifier = 1 };
112 template <
class DiscreteFunctionSpaceImp,
115 struct BubbleM <
CombinedSpace<DiscreteFunctionSpaceImp,N,policy> >
117 enum { bubblePModifier = BubbleM< DiscreteFunctionSpaceImp > :: bubblePModifier };
120 typedef BubbleM <DiscreteFunctionSpaceType> BubbleMType;
122#ifdef USE_TWISTFREE_MAPPER
124 enum { bubblePModifier = BubbleMType :: bubblePModifier };
127 enum { bubblePolOrd = (polOrdN + bubblePModifier) };
130#warning "Hdiv-Projection only working for polOrd = 1 (enable higher order Lagrange with -DUSE_TWISTFREE_MAPPER)"
133 enum { bubblePModifier = 1 };
136 enum { bubblePolOrd = (polOrdN > 1) ? 2 : (polOrdN + 1) };
139 template <
class Space>
struct Spaces;
141 template <
class FunctionSpaceImp,
144 template <
class>
class StorageImp,
145 template <
class,
class,
int,
template <
class>
class>
class DiscreteFunctionSpaceImp>
146 struct Spaces<DiscreteFunctionSpaceImp<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> >
150 typedef DiscreteFunctionSpaceImp<ElSpaceType,GridPartImp,gradPolOrd,StorageImp> ElementGradientSpaceType;
152 typedef DiscreteFunctionSpaceImp<FaceSpaceType,GridPartType,polOrdN,StorageImp> FaceDiscreteSpaceType;
158 template <
class DiscreteFunctionSpaceImp,
161 struct Spaces<
CombinedSpace<DiscreteFunctionSpaceImp,N,policy> >
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;
170 typedef Spaces<DiscreteFunctionSpaceType> AllSpacesType;
171 typedef typename AllSpacesType :: ElementGradientSpaceType ElementGradientSpaceType;
172 typedef typename AllSpacesType :: ElementDiscreteSpaceType ElementDiscreteSpaceType;
173 typedef typename AllSpacesType :: FaceDiscreteSpaceType FaceDiscreteSpaceType;
175 const DiscreteFunctionSpaceType& space_;
176 GridPartType & gridPart_;
177 const FaceDiscreteSpaceType faceSpace_;
178 const ElementDiscreteSpaceType elSpace_;
179 const ElementGradientSpaceType gradSpace_;
183 HdivProjection(
const DiscreteFunctionSpaceType&
space) DUNE_DEPRECATED_MSG(
"Use dune-fem-dg implementation." ) :
185 gridPart_(
space.gridPart()),
186 faceSpace_( gridPart_ ),
187 elSpace_( gridPart_ ),
188 gradSpace_( gridPart_ )
193 gridPart_( org.gridPart_),
194 faceSpace_( gridPart_ ),
195 elSpace_( gridPart_ ),
196 gradSpace_( gridPart_ )
200 const DiscreteFunctionSpaceType&
space()
const
211 static double normalJump(
const DiscreteFunctionType &discFunc,
const int polyOrder = -1 )
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;
218 enum { dim = GridType::dimension };
220 const DiscreteFunctionSpaceType &
space = discFunc.space();
221 const GridPartType &gridPart =
space.gridPart();
222 const int polOrd = (polyOrder <0) ? (2 *
space.order() + 2) : polyOrder;
224 typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
227 RangeType neighRet (0.0);
234 const LocalIdSetType &idSet = gridPart.grid().localIdSet();
236 for(
const auto & en :
space)
238 const LocalFuncType lf = discFunc.localFunction(en);
240 double localValue = 0.0;
242 IntersectionIteratorType endnit = gridPart.iend(en);
243 for(IntersectionIteratorType nit = gridPart.ibegin(en);
244 nit != endnit; ++nit)
246 const IntersectionType& inter=*nit;
248 if(inter.neighbor() )
250 EntityType nb = inter.outside();
252 if(idSet.id( en ) < idSet.id( nb ))
254 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
255 FaceQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
257 const LocalFuncType neighLf = discFunc.localFunction(nb);
259 const int quadNop = faceQuadInner.nop();
260 for (
int l = 0; l < quadNop ; ++l)
263 inter.unitOuterNormal(faceQuadInner.localPoint(l));
265 lf.evaluate(faceQuadInner[l], ret);
266 neighLf.evaluate(faceQuadOuter[l], neighRet);
270 double val = ret * normal;
271 val *= faceQuadInner.weight(l);
287 void curl(
const DomainType & arg, DomainType & dest,
const int d )
const
289 if( DomainType :: dimension == 2 )
296 else if( DomainType :: dimension == 3 )
327 template <
class EntityType,
328 class LocalFunctionType,
332 void bubblePart(
const ElementDiscreteSpaceType&
space,
334 const LocalFunctionType & uLF,
const int startRow ,
336 MatrixType & matrix, VectorType& rhs )
const
339 typedef typename ElementDiscreteSpaceType :: BaseFunctionSetType BaseFunctionSetType;
340 typedef typename ElementDiscreteSpaceType :: LagrangePointSetType LagrangePointSetType;
342 typedef typename ElementDiscreteSpaceType :: JacobianRangeType JacobianRangeType;
343 typedef typename ElementDiscreteSpaceType :: DomainType DomainType;
345 enum { dim = GridType::dimension };
347 const LagrangePointSetType& lagrangePointSet =
space.lagrangePointSet( en );
348 const BaseFunctionSetType bSet =
space.baseFunctionSet( en );
349 const int polOrd = 2 *
space.order();
351 const int cols = uLF.numDofs();
352 assert( uRets.size() == (
unsigned int)cols );
354 VolumeQuadratureType quad (en,polOrd);
356 std::vector< JacobianRangeType > valTmpVec( bSet.size() );
361 typedef typename EntityType :: Geometry Geometry ;
362 const Geometry& geo = en.geometry();
364 const GeometryType& type = geo.type();
365 const int bubbleOffset = (type.isSimplex()) ? 0 : baseFunctionOffset( 0 );
368 enum { cdim = Geometry :: coorddimension };
369 enum { mydim = Geometry :: mydimension };
370 typedef typename Geometry::JacobianInverseTransposed JacobianInverseType;
373 const int enDofs = numberOfBubbles( lagrangePointSet.numDofs( 0 ), type ,
375 const int bubbleMod = bubbleModifier( mydim );
377 const int quadNop = quad.nop();
378 for (
int l = 0; l < quadNop ; ++l)
381 const JacobianInverseType& inv = geo.jacobianInverseTransposed( quad.point(l) );
384 const double intel = quad.weight(l) *
385 geo.integrationElement(quad.point(l));
388 uLF.evaluate(quad[l], result);
391 uLF.baseFunctionSet().evaluateAll( quad[l], uRets );
394 bSet.jacobianAll( quad[l], inv, valTmpVec );
397 for(
int i = 0 ; i<enDofs; i += bubbleMod )
400 int row = startRow + i;
403 const int localBaseFct = ((int) i/bubbleMod) + bubbleOffset;
405 const int baseFct = lagrangePointSet.entityDofNumber( 0, 0, localBaseFct );
408 JacobianRangeType& val = valTmpVec[ baseFct ];
413 for(
int d = 0; d<bubbleMod; ++d )
416 curl(val[0], aVal, d );
418 double r = aVal * result;
423 for(
int j=0; j<cols; ++j)
425 double t = aVal * uRets[ j ];
427 matrix[ row ][ j ] += t;
437 template <
class EntityType,
class LocalFunctionType,
class ArrayType,
438 class MatrixType,
class VectorType>
439 void gradientPart(
const ElementGradientSpaceType &
space,
441 const LocalFunctionType & uLF,
const int startRow ,
442 ArrayType& uRets, MatrixType& matrix, VectorType& rhs )
const
444 typedef typename ElementGradientSpaceType :: BaseFunctionSetType BaseFunctionSetType;
445 const BaseFunctionSetType bSet =
space.baseFunctionSet( en );
446 int polOrd = 2 *
space.order() + 1;
448 const int localRows = gradientBaseFct( bSet );
449 const int cols = uLF.numDofs();
451 VolumeQuadratureType quad (en,polOrd);
456 typedef typename ElementGradientSpaceType :: JacobianRangeType GradJacobianRangeType;
457 std::vector< GradJacobianRangeType > gradTmpVec( bSet.size() );
458 GradJacobianRangeType gradPhi;
460 typedef typename EntityType :: Geometry Geometry ;
461 const Geometry& geo = en.geometry();
463 enum { cdim = Geometry :: coorddimension };
464 enum { mydim = Geometry :: mydimension };
465 typedef typename Geometry::JacobianInverseTransposed JacobianInverseType;
467 const int quadNop = quad.nop();
468 for (
int l = 0; l < quadNop ; ++l)
471 const JacobianInverseType& inv = geo.jacobianInverseTransposed( quad.point( l ) );
474 const double intel = quad.weight(l) *
475 geo.integrationElement( quad.point( l ) );
478 uLF.evaluate(quad[l], result);
481 uLF.baseFunctionSet().evaluateAll( quad[l], uRets );
485 bSet.jacobianAll( quad[l], inv, gradTmpVec );
488 for(
int i=0; i<localRows; ++i)
491 const int row = startRow + i;
495 GradJacobianRangeType& gradPhi = gradTmpVec[ baseFunctionOffset( i ) ];
497 const double uDGVal = result * gradPhi[0];
498 rhs[row] += uDGVal * intel;
501 for(
int j=0; j<cols; ++j)
505 const double val = uRets[ j ] * gradPhi[0];
506 matrix[row][j] += val * intel;
512 template <
class FaceBSetType,
class Gr
idType>
513 struct GetSubBaseFunctionSet
515 template <
class EntityType,
class SpaceType>
516 static inline FaceBSetType faceBaseSet(
const EntityType& en,
const SpaceType&
space)
518 return space.baseFunctionSet( (en.template subEntity<1> (0) )->type() );
522 template <
class FaceBSetType,
int dim,
class CoordCont>
523 struct GetSubBaseFunctionSet< FaceBSetType, YaspGrid< dim, CoordCont > >
525 template <
class EntityType,
class SpaceType>
526 static inline FaceBSetType faceBaseSet(
const EntityType& en,
const SpaceType&
space)
528 return space.baseFunctionSet( GeometryTypes::cube(dim-1) );
533 template <
class FaceBSetType,
int dim>
534 struct GetSubBaseFunctionSet<FaceBSetType, UGGrid<dim> >
536 template <
class EntityType,
class SpaceType>
537 static inline FaceBSetType faceBaseSet(
const EntityType& en,
const SpaceType&
space)
539 const GeometryType geoType (en.geometry().type().basicType(),dim-1);
540 return space.baseFunctionSet( geoType );
545 enum { gradFuncOffset = 1 };
546 template <
class GradBaseFunctionSet>
547 int gradientBaseFct(
const GradBaseFunctionSet& gradSet)
const
549 return (gradPolOrd <= 0) ? 0 : gradSet.size() - gradFuncOffset;
552 int baseFunctionOffset(
const int i)
const
554 return i + gradFuncOffset;
557 int numberOfBubbles(
const int bubbles ,
const GeometryType& type,
558 const int allDofs,
const int allOther )
const
570 const int numBubble = allDofs - allOther ;
571 return (numBubble > 0) ? numBubble : 0;
575 int bubbleModifier(
const int dim )
const
578 return (dim - 2) * (dim - 1) + 1;
582 void project(
const DiscreteFunctionType &uDG,
583 DiscreteFunctionType & velo )
const
585 typedef typename DiscreteFunctionType::Traits::DiscreteFunctionSpaceType FunctionSpaceType;
586 enum { localBlockSize = FunctionSpaceType::localBlockSize };
588 typedef typename FunctionSpaceType::GridType GridType;
589 typedef typename FunctionSpaceType::GridPartType GridPartType;
590 typedef typename FunctionSpaceType::IteratorType Iterator;
592 typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
593 typedef typename FunctionSpaceType::RangeType RangeType;
595 enum { dim = GridType::dimension };
596 typedef typename GridType :: ctype coordType;
598 const FunctionSpaceType&
space = uDG.space();
602 if(
space.order() < 1 )
608 const int polOrd = 2 *
space.order() + 2;
610 typedef typename FaceDiscreteSpaceType :: BaseFunctionSetType FaceBSetType ;
612 typedef typename ElementGradientSpaceType :: BaseFunctionSetType GradientBaseSetType ;
614 typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
616 Iterator start =
space.begin();
618 if( start ==
space.end() )
return;
621 if(
space.multipleGeometryTypes() )
623 if(
space.indexSet().geomTypes(0).size() > 1)
625 assert(
space.indexSet().geomTypes(0).size() == 1 );
626 DUNE_THROW(NotImplemented,
"H-div projection not implemented for hybrid grids");
630 const GeometryType startType = start->type();
633 if( startType.isHexahedron() &&
space.order() > 1 )
635 assert( ! startType.isHexahedron() ||
space.order() <= 1 );
636 DUNE_THROW(NotImplemented,
"H-div projection not implemented for p > 1 on hexas! ");
639 const int desiredOrder =
space.order() + bubblePModifier;
641 if( bubblePolOrd != desiredOrder )
643 assert( bubblePolOrd == desiredOrder );
644 DUNE_THROW(NotImplemented,
"H-div projection not working for "
645 <<
space.order() <<
" when LagrangeSpace of order "<< desiredOrder <<
" is missing");
649 LocalFuncType lf = uDG.localFunction(*start);
650 const int numDofs = lf.numDofs();
653 const FaceBSetType faceSet =
654 GetSubBaseFunctionSet<FaceBSetType,GridType>::faceBaseSet( *start , faceSpace_ );
656 const int numFaceDofs = faceSet.size();
658 const GradientBaseSetType gradSet = gradSpace_.baseFunctionSet(*start);
660 const int numGradDofs = gradientBaseFct( gradSet );
662 const Dune::ReferenceElement< coordType, dim > & refElem =
663 Dune::ReferenceElements< coordType, dim >::general( startType );
666 const int overallFaceDofs = numFaceDofs * refElem.size(1);
669 const int numBubbleDofs = (
space.order() <= 1) ? 0 :
670 numberOfBubbles( elSpace_.lagrangePointSet( *start ).numDofs ( 0 ) , startType,
671 numDofs, overallFaceDofs + numGradDofs );
673 const int rows = (overallFaceDofs + numGradDofs + numBubbleDofs);
676 const int cols = numDofs;
678 DynamicArray< RangeFieldType > rets(numDofs);
679 DynamicArray< RangeType > uRets(numDofs);
681 typedef FieldMatrix<RangeFieldType,localBlockSize,localBlockSize> FieldMatrixType;
686 typedef FieldVector<RangeFieldType,localBlockSize> VectorType;
687 VectorType fRhs(0.0);
689 assert( numDofs == localBlockSize );
690 if( numDofs != localBlockSize )
692 DUNE_THROW(InvalidStateException,
"wrong sizes ");
696 const bool nonSymetric = (cols != rows);
699 typedef Fem::DenseMatrix<double> MatrixType;
700 MatrixType matrix(rows,cols);
701 typedef typename MatrixType :: RowType RowType;
703 MatrixType fakeMatrix(cols,cols);
704 RowType rhs(rows,0.0);
705 RowType fakeRhs(numDofs,0.0);
708 for(
const auto & en :
space)
716 for(
int i=0; i<rows; ++i)
722 fillMatrix(gridPart_,en,uDG,
724 gradSpace_, overallFaceDofs,
725 elSpace_, rows - numBubbleDofs,
726 polOrd,numDofs,numFaceDofs,
727 rets,uRets, matrix,rhs);
730 matrix.multTransposed(rhs, fakeRhs);
731 fakeMatrix.multiply_AT_A(matrix);
734 for(
int i=0; i<numDofs; ++i)
736 fRhs[i] = fakeRhs[i];
737 for(
int j=0; j<numDofs; ++j)
739 inv[i][j] = fakeMatrix[i][j];
750 assert( cols == rows );
752 fillMatrix(gridPart_,en,uDG,
754 gradSpace_, rows - numBubbleDofs - numGradDofs,
755 elSpace_, rows - numBubbleDofs,
756 polOrd,numDofs,numFaceDofs,
757 rets, uRets, inv,fRhs);
762 LocalFuncType veloLF = velo.localFunction( en );
765 luSolve( inv, fRhs );
766 const VectorType& x = fRhs ;
768 for(
int i=0; i<localBlockSize; ++i)
770 veloLF[ i ] = x[ i ];
776 template <
class GridPartType,
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
792 typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
793 typedef typename IntersectionIteratorType::Intersection IntersectionType;
795 typedef typename FaceDiscreteSpaceType :: BaseFunctionSetType FaceBSetType ;
796 typedef typename FaceDiscreteSpaceType :: RangeType FaceRangeType;
797 FaceRangeType faceVal;
799 typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
801 RangeType neighRet (0.0);
802 RangeType uPhi (0.0);
804 typedef typename DiscreteFunctionType :: LocalFunctionType LocalFuncType ;
805 typedef typename DiscreteFunctionType :: DiscreteFunctionSpaceType
806 :: BaseFunctionSetType BaseFunctionSetType;
809 const LocalFuncType uLF = uDG.localFunction(en);
812 const BaseFunctionSetType & bSet = uLF.baseFunctionSet();
815 IntersectionIteratorType endnit = gridPart.iend(en);
816 for(IntersectionIteratorType nit = gridPart.ibegin(en);
817 nit != endnit; ++nit)
819 const IntersectionType& inter=*nit;
821 const FaceBSetType &faceSet = faceSpace.baseFunctionSet( inter.type() );
823 const int firstRow = inter.indexInInside() * numFaceDofs;
829 EntityType nb = inter.outside();
832 const LocalFuncType uNeighLf = uDG.localFunction(nb);
837 if( inter.conforming() )
840 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
841 FaceQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
843 applyLocalNeighbor(nit,faceQuadInner,faceQuadOuter,
844 bSet,faceSet, uLF, uNeighLf,
845 firstRow, numFaceDofs,
847 ret,neighRet,faceVal,
853 typedef typename FaceQuadratureType ::
854 NonConformingQuadratureType NonConformingQuadratureType;
856 NonConformingQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
857 NonConformingQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
859 applyLocalNeighbor(nit,faceQuadInner,faceQuadOuter,
860 bSet,faceSet, uLF, uNeighLf,
861 firstRow, numFaceDofs,
863 ret,neighRet,faceVal,
873 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
874 const int quadNop = faceQuadInner.nop();
876 std::vector< RangeType > uPhiVec( numDofs );
877 std::vector< FaceRangeType > faceValVec( numFaceDofs );
879 for (
int l = 0; l < quadNop ; ++l)
881 DomainType unitNormal =
882 inter.integrationOuterNormal(faceQuadInner.localPoint(l));
884 const double faceVol = unitNormal.two_norm();
885 unitNormal *= 1.0/faceVol;
888 const double intel = faceVol * faceQuadInner.weight(l);
891 uLF.evaluate(faceQuadInner[l], ret);
893 RangeFieldType val = ret * unitNormal;
896 bSet.evaluateAll( faceQuadInner[l], uPhiVec );
899 for(
int i=0; i<numDofs; ++i)
901 rets[i] = uPhiVec[ i ] * unitNormal;
905 faceSet.evaluateAll( faceQuadInner[ l ], faceValVec );
908 for(
int j=0; j<numFaceDofs; ++j, ++row)
910 FaceRangeType& faceVal = faceValVec[ j ];
911 rhs[row] += val*faceVal[0];
913 for(
int i=0; i<numDofs; ++i)
915 matrix[row][i] += (faceVal[0] * rets[i]);
925 gradientPart(gradSpace, en, uLF, startGradDofs, uRets, matrix, rhs );
929 if( bubblePolOrd - bubblePModifier > 1 )
931 bubblePart(elSpace, en, uLF, startBubbleDofs, uRets, matrix, rhs);
937 template <
class MatrixType>
938 void printMatrix(
const MatrixType& matrix)
const
940 std::cout <<
"Print Matrix \n";
941 for(
size_t row = 0; row < matrix.N(); ++row)
943 std::cout << row <<
": ";
944 for(
size_t col = 0; col< matrix.M(); ++col)
946 if(
std::abs( matrix[row][col] ) < 1e-12 )
949 std::cout << matrix[row][col] <<
" ";
951 std::cout << std::endl;
953 std::cout <<
"Finished print Matrix \n";
956 template <
class IntersectionIteratorType,
957 class QuadratureType,
958 class BaseFunctionSetType,
959 class FaceBaseFunctionSetType,
960 class LocalFunctionType,
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,
974 const int numFaceDofs,
976 RangeType& ret, RangeType& neighRet,
977 FaceRangeType& faceVal,
981 const typename IntersectionIteratorType::Intersection& inter = *nit;
982 const int quadNop = faceQuadInner.nop();
983 const int numDofs = uLF.numDofs();
985 std::vector< RangeType > phiVec( numDofs );
986 std::vector< FaceRangeType > facePhiVec( numFaceDofs );
988 for (
int l = 0; l < quadNop ; ++l)
990 DomainType unitNormal =
991 inter.integrationOuterNormal(faceQuadInner.localPoint(l));
994 const double faceVol = unitNormal.two_norm();
995 unitNormal *= 1.0/faceVol;
998 const double intel = faceVol * faceQuadInner.weight(l);
1001 uLF.evaluate(faceQuadInner[l], ret);
1002 uNeighLf.evaluate(faceQuadOuter[l], neighRet);
1008 RangeFieldType val = ret * unitNormal;
1011 bSet.evaluateAll( faceQuadInner[l], phiVec );
1014 for(
int i=0; i<numDofs; ++i)
1016 rets[i] = phiVec[ i ] * unitNormal;
1021 faceSet.evaluateAll( faceQuadInner[ l ], facePhiVec );
1024 for(
int j=0; j<numFaceDofs; ++j, ++row )
1026 FaceRangeType& faceVal = facePhiVec[ j ];
1028 rhs[row] += val * faceVal[0];
1030 for(
int i=0; i<numDofs; ++i)
1032 matrix[row][i] += (faceVal[0] * rets[i]);
1041 DiscreteFunctionType& dest)
const
1047 template <
class AdaptationType>
1049 AdaptationType& adaptation)
1051 typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
1052 typedef typename GridType :: Traits :: LocalIdSet LocalIdSetType;
1054 enum { dim = GridType :: dimension };
1055 typedef typename GridType :: template Codim<0> :: Entity EntityType;
1057 const DiscreteFunctionSpaceType&
space = velo.space();
1058 GridPartType & gridPart =
space.gridPart();
1059 int polOrd =
space.order() + 1;
1062 const LocalIdSetType& localIdSet = gridPart.grid().localIdSet();
1067 for(
const auto & en :
space)
1069 const LocalFuncType uLF = velo.localFunction(en);
1070 const double enVol = en.geometry().volume();
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)
1078 const IntersectionType& inter=*nit;
1079 double enError = 0.0;
1081 if(inter.neighbor())
1083 EntityType nb = inter.outside();
1084 const double enVol_nbVol = 0.5 * (enVol + nb.geometry().volume());
1088 const bool interiorEntity = (nb.partitionType() == InteriorEntity);
1090 const bool interiorEntity =
true;
1092 if( (localIdSet.id( en ) < localIdSet.id( nb ))
1094 || ( ! interiorEntity )
1098 const LocalFuncType uNeighLf = velo.localFunction(nb);
1103 if( inter.conforming() )
1106 FaceQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
1107 FaceQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
1109 applyLocalNeighEstimator(nit,nb,faceQuadInner,faceQuadOuter,
1110 uLF, uNeighLf, enVol_nbVol, interiorEntity,
1111 enError, adaptation);
1116 typedef typename FaceQuadratureType ::
1117 NonConformingQuadratureType NonConformingQuadratureType;
1119 NonConformingQuadratureType faceQuadInner(gridPart, inter, polOrd, FaceQuadratureType::INSIDE);
1120 NonConformingQuadratureType faceQuadOuter(gridPart, inter, polOrd, FaceQuadratureType::OUTSIDE);
1122 applyLocalNeighEstimator(nit,nb,faceQuadInner,faceQuadOuter,
1123 uLF, uNeighLf, enVol_nbVol, interiorEntity,
1124 enError, adaptation);
1132 adaptation.addToLocalIndicator( en , enError );
1139 template <
class IntersectionIteratorType,
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,
1153 AdaptationType& adaptation)
1155 const typename IntersectionIteratorType::Intersection& inter=*nit;
1156 enum { dim = GridType :: dimension };
1159 double nbError = 0.0;
1161 const int quadNop = faceQuadInner.nop();
1162 for (
int l = 0; l < quadNop ; ++l)
1164 DomainType unitNormal =
1165 inter.integrationOuterNormal(faceQuadInner.localPoint(l));
1167 double faceVol = unitNormal.two_norm();
1168 unitNormal *= 1.0/faceVol;
1173 const double power = (1.0/(dim-1));
1174 faceVol =
pow(faceVol, power);
1178 uLF.evaluate(faceQuadInner[l], jump);
1179 uNeighLf.evaluate(faceQuadOuter[l], neighRet);
1184 const double weight = (enVol_nbVol) * faceQuadInner.weight(l);
1186 double error = weight * SQR(jump * unitNormal);
1200 adaptation.addToLocalIndicator( nb , nbError );
1209 template <
class MatrixType,
class VectorType>
1210 void luSolve(MatrixType& a,
1211 VectorType& x)
const
1213 typedef typename VectorType :: field_type ctype;
1214 enum { n = VectorType :: dimension };
1217 assert( a.N() == a.M() );
1218 assert( a.N() == n );
1223 for(
int i=0; i<n-1; ++i)
1230 for(
int k=i; k<n; ++k)
1234 max_abs = fabs(a[k][i]);
1242 for(
int j=0; j<n; ++j)
1244 const ctype tmp = a[ p[i] ][j];
1245 a[ p[i] ][j] = a[i][j];
1251 for(
int k=i+1; k<n; ++k)
1253 const ctype lambda = a[k][i] / a[i][i];
1255 for(
int j=i+1; j<n; ++j) a[k][j] -= lambda * a[i][j];
1260 for(
int i=0; i<n-1; ++i)
1262 const ctype tmp = x[ i ];
1263 x[ i ] = x[ p[ i ] ];
1268 for(
int i=0; i<n; ++i)
1271 for(
int j=0; j<i; ++j) dot += a[i][j] * x[j];
1276 for(
int i=n-1; i>=0; --i)
1279 for(
int j=i+1; j<n; ++j) dot += a[i][j] * x[j];
1280 x[i] = (x[i] - dot) / a[i][i];
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