dune-fem  2.8-git
common/scalarproducts.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SCALARPRODURCTS_HH
2 #define DUNE_FEM_SCALARPRODURCTS_HH
3 
4 #include <iostream>
5 #include <memory>
6 #include <set>
7 #include <map>
8 #include <limits>
9 #include <algorithm>
10 
11 #include <dune/common/exceptions.hh>
12 #include <dune/common/genericiterator.hh>
13 #include <dune/common/ftraits.hh>
14 #include <dune/common/version.hh>
15 
16 #include <dune/grid/common/gridenums.hh>
17 #include <dune/grid/common/datahandleif.hh>
18 
19 #if HAVE_DUNE_ISTL
20 #include <dune/istl/scalarproducts.hh>
21 #endif
22 
30 
31 namespace Dune
32 {
33 
34  namespace Fem
35  {
36 
41 #if HAVE_DUNE_ISTL
42  template <class DofVector>
43  struct ISTLScalarProductSelector
44  {
45  typedef Dune::FieldVector< typename DofVector::FieldType, DofVector::blockSize > Block;
46  typedef Dune::BlockVector< Block > type;
47  };
48 
49  template <class Block>
50  struct ISTLScalarProductSelector< Dune::Fem::ISTLBlockVector< Block > >
51  : public Dune::ScalarProduct< typename Dune::Fem::ISTLBlockVector< Block > :: DofContainerType >
52  {
53 #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
55  enum { category=Dune::SolverCategory::sequential };
56 #endif // #if !DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
57 
58  typedef typename ISTLBlockVector< Block > :: DofContainerType type;
59 
60 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
61  Dune::SolverCategory::Category category () const override { return SolverCategory::sequential; }
62 #endif // #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
63  };
64 #endif
65 
68  template< class DiscreteFunction >
70 #if HAVE_DUNE_ISTL
71  : public ISTLScalarProductSelector< typename DiscreteFunction :: DofVectorType >
72 #endif
73  {
74  public:
75  typedef DiscreteFunction DiscreteFunctionType;
76 
78  typedef typename DiscreteFunctionType :: DiscreteFunctionSpaceType
80 
81  private:
83 
84  public:
86  typedef typename DiscreteFunctionSpaceType :: RangeFieldType RangeFieldType;
87 
89  typedef typename DiscreteFunctionSpaceType :: BlockMapperType MapperType;
90 
91  // type of communication manager object which does communication
93 
95  typedef typename Dune::FieldTraits< RangeFieldType >::real_type real_type;
96 
99  : space_( space )
100  {}
101 
103  {
104  return space_;
105  }
106 
108  template < class OtherDiscreteFunctionType >
109  RangeFieldType scalarProductDofs ( const DiscreteFunctionType &x, const OtherDiscreteFunctionType &y ) const
110  {
111  assert(x.space() == y.space());
112  assert(x.space() == space());
113  return dotProduct( x.dofVector(), y.dofVector() );
114  }
115 
117  {
118  return space().auxiliaryDofs();
119  }
120 
121  protected:
123  template < class DofVector, class OtherDofVector >
124  RangeFieldType dotProduct ( const DofVector &x, const OtherDofVector &y ) const
125  {
126  typedef typename DiscreteFunctionSpaceType::LocalBlockIndices LocalBlockIndices;
127 
128  RangeFieldType scp = 0;
129  for( const auto i : primaryDofs( space().auxiliaryDofs() ) )
130  Hybrid::forEach( LocalBlockIndices(), [ &x, &y, &scp, i ] ( auto &&j ) { scp += x[ i ][ j ] * y[ i ][ j ]; } );
131  return space().gridPart().comm().sum( scp );
132  }
133 
134 #if HAVE_DUNE_ISTL
135  protected:
136  typedef typename ISTLScalarProductSelector< typename DiscreteFunction :: DofVectorType > :: type BlockVectorType;
137 
138  public:
140  virtual field_type dot (const BlockVectorType& x,
141  const BlockVectorType& y) const
142  {
143  return dotProduct( x, y );
144  }
145 
147  virtual real_type norm( const BlockVectorType& x ) const
148  {
149  return std::abs( std::sqrt( dotProduct( x, x ) ) );
150  }
151 
153  void deleteNonInterior( BlockVectorType& x) const
154  {
155 #if HAVE_MPI
156  // case of ALUGrid and DGSpace or FVSpace
157  // BUG: We should not use the leafGridView to detect whether the grid has overlap!
158  const bool deleteGhostEntries = (space().gridPart().grid().leafGridView().overlapSize( 0 ) == 0) && !space().continuous();
159 
160  // only delete ghost entries
161  if( deleteGhostEntries )
162  {
163  const auto &auxiliaryDofs = space().auxiliaryDofs();
164 
165  // don't delete the last since this is the overall Size
166  const int auxiliarySize = auxiliaryDofs.size() - 1;
167  for(int auxiliary = 0; auxiliary<auxiliarySize; ++auxiliary)
168  x[ auxiliaryDofs[auxiliary] ] = 0;
169  }
170 #endif // #if HAVE_MPI
171  }
172 #endif // #if HAVE_DUNE_ISTL
174  };
175 
177 
178  } // end namespace Fem
179 
180 } // end namespace Dune
181 #endif // #ifndef DUNE_FEM_SCALARPRODURCTS_HH
const DiscreteFunctionSpaceType & space() const
Definition: common/scalarproducts.hh:102
const DiscreteFunctionSpaceType & space_
Definition: common/scalarproducts.hh:173
static PrimaryDofs< AuxiliaryDofs > primaryDofs(const AuxiliaryDofs &auxiliaryDofs)
Definition: auxiliarydofs.hh:341
RangeFieldType scalarProductDofs(const DiscreteFunctionType &x, const OtherDiscreteFunctionType &y) const
evaluate scalar product and omit auxiliary nodes
Definition: common/scalarproducts.hh:109
AuxiliaryDofs< typename DiscreteFunctionSpaceType::GridPartType, MapperType > AuxiliaryDofsType
Definition: common/scalarproducts.hh:92
DiscreteFunctionType ::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
type of the discrete function space
Definition: common/scalarproducts.hh:79
RangeFieldType field_type
Definition: common/scalarproducts.hh:94
DiscreteFunctionSpaceType ::RangeFieldType RangeFieldType
type of range field
Definition: common/scalarproducts.hh:86
const AuxiliaryDofsType & auxiliaryDofs() const
Definition: common/scalarproducts.hh:116
DiscreteFunction DiscreteFunctionType
Definition: common/scalarproducts.hh:75
ParallelScalarProduct(const DiscreteFunctionSpaceType &space)
constructor taking space
Definition: common/scalarproducts.hh:98
DiscreteFunctionSpaceType ::BlockMapperType MapperType
type of used mapper
Definition: common/scalarproducts.hh:89
RangeFieldType dotProduct(const DofVector &x, const OtherDofVector &y) const
evaluate scalar product on dofVector and omit auxiliary nodes
Definition: common/scalarproducts.hh:124
int size() const
return number of auxiliary dofs
Definition: auxiliarydofs.hh:122
Dune::FieldTraits< RangeFieldType >::real_type real_type
Definition: common/scalarproducts.hh:95
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
static double sqrt(const Double &v)
Definition: double.hh:886
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
Definition: common/scalarproducts.hh:73
In parallel computations the dofs of a discrete function are made up by all primary dofs....
Definition: auxiliarydofs.hh:47