adds DGFunction::dotGlobal and DGFunction::dotLocal

src/hyteg/dgfunctionspace/DGEdge.hpp

@@ -554,6 +554,39 @@ inline void add( const uint_t&                                               lev
    }
 }
 
+template < typename ValueType >
+inline ValueType dot( const uint_t&                                                 level,
+                      Edge&                                                         edge,
+                      const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& lhsMemoryId,
+                      const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& rhsMemoryId )
+{
+   walberla::math::KahanAccumulator< ValueType > scalarProduct;
+
+   ValueType* lhsPtr = edge.getData( lhsMemoryId )->getPointer( level );
+   ValueType* rhsPtr = edge.getData( rhsMemoryId )->getPointer( level );
+
+   size_t rowsize = levelinfo::num_microvertices_per_edge( level );
+
+   // gray south cells
+   for ( uint_t i = 1; i < rowsize - 2; ++i )
+   {
+      uint_t index = facedof::macroedge::indexFaceFromVertex( level, i, stencilDirection::CELL_GRAY_SE );
+      scalarProduct += lhsPtr[index] * rhsPtr[index];
+   }
+
+   if ( edge.getNumNeighborFaces() == 2 )
+   {
+      // gray north cells
+      for ( uint_t i = 1; i < rowsize - 2; ++i )
+      {
+         uint_t index = facedof::macroedge::indexFaceFromVertex( level, i, stencilDirection::CELL_GRAY_NE );
+         scalarProduct += lhsPtr[index] * rhsPtr[index];
+      }
+   }
+
+   return scalarProduct.get();
+}
+
 template < typename ValueType >
 inline void swap( const uint_t&                                               level,
                   Edge&                                                       edge,

src/hyteg/dgfunctionspace/DGFace.hpp

@@ -694,6 +694,46 @@ inline void add( const uint_t&                                               lev
    }
 }
 
+template < typename ValueType >
+inline ValueType dot( const uint_t&                                                 level,
+                      Face&                                                         face,
+                      const PrimitiveDataID< FunctionMemory< ValueType >, Face >& lhsMemoryId,
+                      const PrimitiveDataID< FunctionMemory< ValueType >, Face >& rhsMemoryId )
+{
+   walberla::math::KahanAccumulator< ValueType > scalarProduct;
+
+   ValueType* lhsPtr = face.getData( lhsMemoryId )->getPointer( level );
+   ValueType* rhsPtr = face.getData( rhsMemoryId )->getPointer( level );
+
+   size_t rowsize       = levelinfo::num_microvertices_per_edge( level );
+   size_t inner_rowsize = rowsize;
+
+   // gray cells
+   for ( size_t j = 1; j < rowsize - 2; ++j )
+   {
+      for ( size_t i = 1; i < inner_rowsize - 3; ++i )
+      {
+         auto cellIndex = facedof::macroface::indexFaceFromGrayFace( level, i, j, stencilDirection::CELL_GRAY_C );
+         scalarProduct += lhsPtr[cellIndex] * rhsPtr[cellIndex];
+      }
+      --inner_rowsize;
+   }
+
+   // blue cells
+   inner_rowsize = rowsize;
+   for ( size_t j = 0; j < rowsize - 2; ++j )
+   {
+      for ( size_t i = 0; i < inner_rowsize - 2; ++i )
+      {
+         auto cellIndex = facedof::macroface::indexFaceFromBlueFace( level, i, j, stencilDirection::CELL_BLUE_C );
+         scalarProduct += lhsPtr[cellIndex] * rhsPtr[cellIndex];
+      }
+      --inner_rowsize;
+   }
+
+   return scalarProduct.get();
+}
+
 template < typename ValueType >
 inline void swap( const uint_t&                                               level,
                   Face&                                                       face,

src/hyteg/dgfunctionspace/DGFunction.hpp

@@ -168,15 +168,9 @@ class DGFunction final : public Function< DGFunction< ValueType > >
 
    void setBoundaryCondition( BoundaryCondition bc ){WALBERLA_ABORT( "DGFunction::setBoundaryCondition not implemented!" )}
 
-   ValueType dotLocal( const DGFunction< ValueType >& secondOp, uint_t level, DoFType flag ) const
-   {
-      WALBERLA_ABORT( "DGFunction::dotLocal not implemented!" )
-   }
+   inline ValueType dotLocal( const DGFunction< ValueType >& secondOp, uint_t level, DoFType flag ) const;
 
-   ValueType dotGlobal( const DGFunction< ValueType >& secondOp, uint_t level, DoFType flag ) const
-   {
-      WALBERLA_ABORT( "DGFunction::dotGlobal not implemented!" )
-   }
+   inline ValueType dotGlobal( const DGFunction< ValueType >& secondOp, uint_t level, DoFType flag ) const;
 
    inline void swap( const DGFunction< ValueType >& other, const uint_t& level, const DoFType& flag = All ) const;
 
@@ -685,4 +679,56 @@ inline void DGFunction< ValueType >::swap( const DGFunction< ValueType >& other,
    this->stopTiming( "Swap" );
 }
 
+template < typename ValueType >
+ValueType DGFunction< ValueType >::dotLocal( const DGFunction< ValueType >& secondOp, uint_t level, DoFType flag ) const
+{
+   this->startTiming( "Dot (local)" );
+
+   walberla::math::KahanAccumulator< ValueType > scalarProduct;
+
+   for ( const auto& it : this->getStorage()->getVertices() )
+   {
+      Vertex& vertex = *it.second;
+
+      if ( testFlag( boundaryCondition_.getBoundaryType( vertex.getMeshBoundaryFlag() ), flag ) )
+      {
+         scalarProduct += DGVertex::dot( level, vertex, vertexDataID_, secondOp.vertexDataID_ );
+      }
+   }
+
+   for ( const auto& it : this->getStorage()->getEdges() )
+   {
+      Edge& edge = *it.second;
+
+      if ( testFlag( boundaryCondition_.getBoundaryType( edge.getMeshBoundaryFlag() ), flag ) )
+      {
+         scalarProduct += DGEdge::dot< ValueType >( level, edge, edgeDataID_, secondOp.edgeDataID_ );
+      }
+   }
+
+   for ( const auto& it : this->getStorage()->getFaces() )
+   {
+      Face& face = *it.second;
+
+      if ( testFlag( boundaryCondition_.getBoundaryType( face.getMeshBoundaryFlag() ), flag ) )
+      {
+         scalarProduct += DGFace::dot< ValueType >( level, face, faceDataID_, secondOp.faceDataID_ );
+      }
+   }
+
+   this->stopTiming( "Dot (local)" );
+
+   return scalarProduct.get();
+}
+
+template < typename ValueType >
+ValueType DGFunction< ValueType >::dotGlobal( const DGFunction< ValueType >& secondOp, uint_t level, DoFType flag ) const
+{
+   ValueType scalarProduct = dotLocal( secondOp, level, flag );
+   this->startTiming( "Dot (reduce)" );
+   walberla::mpi::allReduceInplace( scalarProduct, walberla::mpi::SUM, walberla::mpi::MPIManager::instance()->comm() );
+   this->stopTiming( "Dot (reduce)" );
+   return scalarProduct;
+}
+
 } // namespace hyteg

src/hyteg/dgfunctionspace/DGVertex.hpp

@@ -130,6 +130,22 @@ inline void add( const uint_t&                                                 l
       dstPtr[i * 2] += scalar;
 }
 
+template < typename ValueType >
+inline ValueType dot( const uint_t&                                                 level,
+                      Vertex&                                                       vertex,
+                      const PrimitiveDataID< FunctionMemory< ValueType >, Vertex >& lhsMemoryId,
+                      const PrimitiveDataID< FunctionMemory< ValueType >, Vertex >& rhsMemoryId )
+{
+   walberla::math::KahanAccumulator< ValueType > scalarProduct;
+
+   ValueType* lhsPtr = vertex.getData( lhsMemoryId )->getPointer( level );
+   ValueType* rhsPtr = vertex.getData( rhsMemoryId )->getPointer( level );
+   for ( uint_t i = 0; i < vertex.getNumNeighborFaces(); ++i )
+      scalarProduct += rhsPtr[i * 2] * lhsPtr[i * 2];
+
+   return scalarProduct.get();
+}
+
 template< typename ValueType >
 inline void upwind(const uint_t & Level, Vertex &vertex,
                    const std::shared_ptr< PrimitiveStorage >& storage,

tests/hyteg/CMakeLists.txt

@@ -514,6 +514,9 @@ waLBerla_execute_test(NAME DGAddTest)
 waLBerla_compile_test(FILES DGSwapTest.cpp DEPENDS hyteg core)
 waLBerla_execute_test(NAME DGSwapTest)
 
+waLBerla_compile_test(FILES DGDotTest.cpp DEPENDS hyteg core)
+waLBerla_execute_test(NAME DGDotTest)
+
 waLBerla_compile_test(FILES vCycleTest.cpp DEPENDS hyteg core)
 waLBerla_execute_test(NAME vCycleTest1 COMMAND $<TARGET_FILE:vCycleTest> )
 

tests/hyteg/DGDotTest.cpp

+/*
+ * Copyright (c) 2021 Andreas Wagner.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+#include "core/debug/all.h"
+#include "core/mpi/all.h"
+
+#include "hyteg/dataexport/VTKOutput.hpp"
+#include "hyteg/dgfunctionspace/DGFunction.hpp"
+#include "hyteg/mesh/MeshInfo.hpp"
+#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
+
+using namespace hyteg;
+
+using walberla::real_t;
+
+int main( int argc, char** argv )
+{
+   walberla::mpi::Environment MPIenv( argc, argv );
+   walberla::MPIManager::instance()->useWorldComm();
+   walberla::debug::enterTestMode();
+
+   MeshInfo              meshInfo = MeshInfo::meshRectangle( Point2D( { 0, 0 } ), Point2D( { 1, 1 } ), MeshInfo::CRISS, 1, 1 );
+   SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+   const uint_t minLevel = 2;
+   const uint_t maxLevel = 4;
+
+   hyteg::DGFunction< real_t > x( "x", storage, minLevel, maxLevel );
+   hyteg::DGFunction< real_t > y( "y", storage, minLevel, maxLevel );
+
+   x.interpolate( 2, maxLevel, hyteg::All );
+   y.interpolate( 3, maxLevel, hyteg::All );
+
+   const real_t numFaces      = levelinfo::num_microfaces_per_face( maxLevel ) * 2.;
+   const real_t numInnerFaces = numFaces - ( 4 * levelinfo::num_microedges_per_edge( maxLevel ) - 2 );
+
+   WALBERLA_CHECK_FLOAT_EQUAL( x.dotGlobal( y, maxLevel, All ), numFaces * 6. );
+   WALBERLA_CHECK_FLOAT_EQUAL( x.dotGlobal( y, maxLevel, Inner ), numInnerFaces * 6. );
+
+   return 0;
+}