diff --git a/src/TNL/Containers/Algorithms/ArrayOperations.h b/src/TNL/Containers/Algorithms/ArrayOperations.h
index 4861cc460262ec74cad2a35400eb472cf2e551dc..e5d5c0deb04da19cf2c1c3c498ac6c2e3fbbde7e 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperations.h
+++ b/src/TNL/Containers/Algorithms/ArrayOperations.h
@@ -81,10 +81,12 @@ class ArrayOperations< Devices::Cuda >
static void freeMemory( Element* data );
template< typename Element >
+ __cuda_callable__
static void setMemoryElement( Element* data,
const Element& value );
template< typename Element >
+ __cuda_callable__
static Element getMemoryElement( const Element* data );
// TODO: does not make sense for CUDA - remove?
diff --git a/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h b/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h
index f7648dbbcd400136ed5812e14eb84f0352d546d9..429ed9b9549775e1f87cc8e070e5ec0b9e882a84 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h
+++ b/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h
@@ -60,24 +60,32 @@ freeMemory( Element* data )
}
template< typename Element >
-void
+__cuda_callable__ void
ArrayOperations< Devices::Cuda >::
setMemoryElement( Element* data,
const Element& value )
{
TNL_ASSERT_TRUE( data, "Attempted to set data through a nullptr." );
+#ifdef __CUDAARCH__
+ *data = value;
+#else
ArrayOperations< Devices::Cuda >::setMemory( data, value, 1 );
+#endif
}
template< typename Element >
-Element
+__cuda_callable__ Element
ArrayOperations< Devices::Cuda >::
getMemoryElement( const Element* data )
{
TNL_ASSERT_TRUE( data, "Attempted to get data through a nullptr." );
+#ifdef __CUDAARCH__
+ return *data;
+#else
Element result;
ArrayOperations< Devices::Host, Devices::Cuda >::copyMemory< Element, Element, int >( &result, data, 1 );
return result;
+#endif
}
template< typename Element, typename Index >
diff --git a/src/TNL/Containers/StaticVector.h b/src/TNL/Containers/StaticVector.h
index e1cf58ed7174cc28cdbe09a60d48de19353fddce..b6c89eefaf88e72ea1586aec469b43e95b1c1b45 100644
--- a/src/TNL/Containers/StaticVector.h
+++ b/src/TNL/Containers/StaticVector.h
@@ -57,7 +57,11 @@ class StaticVector : public StaticArray< Size, Real >
//! Multiplication with number
__cuda_callable__
StaticVector& operator *= ( const Real& c );
-
+
+ //! Division by number
+ __cuda_callable__
+ StaticVector& operator /= ( const Real& c );
+
//! Addition operator
__cuda_callable__
StaticVector operator + ( const StaticVector& u ) const;
@@ -155,6 +159,10 @@ class StaticVector< 1, Real > : public StaticArray< 1, Real >
//! Multiplication with number
__cuda_callable__
StaticVector& operator *= ( const Real& c );
+
+ //! Division by number
+ __cuda_callable__
+ StaticVector& operator /= ( const Real& c );
//! Addition operator
__cuda_callable__
@@ -257,6 +265,10 @@ class StaticVector< 2, Real > : public StaticArray< 2, Real >
__cuda_callable__
StaticVector& operator *= ( const Real& c );
+ //! Division by number
+ __cuda_callable__
+ StaticVector& operator /= ( const Real& c );
+
//! Adding operator
__cuda_callable__
StaticVector operator + ( const StaticVector& u ) const;
@@ -357,6 +369,11 @@ class StaticVector< 3, Real > : public StaticArray< 3, Real >
//! Multiplication with number
__cuda_callable__
StaticVector& operator *= ( const Real& c );
+
+ //! Division by number
+ __cuda_callable__
+ StaticVector& operator /= ( const Real& c );
+
//! Addition operator
__cuda_callable__
diff --git a/src/TNL/Containers/StaticVector1D_impl.h b/src/TNL/Containers/StaticVector1D_impl.h
index 3127a621910111ff1490ccdc94675c8e737d7fc8..f009b52afe4f7d630ad11c6324e5b5fe291dd71c 100644
--- a/src/TNL/Containers/StaticVector1D_impl.h
+++ b/src/TNL/Containers/StaticVector1D_impl.h
@@ -86,6 +86,14 @@ StaticVector< 1, Real >& StaticVector< 1, Real >::operator *= ( const Real& c )
return *this;
}
+template< typename Real >
+__cuda_callable__
+StaticVector< 1, Real >& StaticVector< 1, Real >::operator /= ( const Real& c )
+{
+ this->data[ 0 ] /= c;
+ return *this;
+}
+
template< typename Real >
__cuda_callable__
StaticVector< 1, Real > StaticVector< 1, Real >::operator + ( const StaticVector& u ) const
diff --git a/src/TNL/Containers/StaticVector2D_impl.h b/src/TNL/Containers/StaticVector2D_impl.h
index d66979bf7a3be3bed89929b3d84cdf5da61199d4..fdcc28aa51efcf1e3bf1da787d0957bf64a39db1 100644
--- a/src/TNL/Containers/StaticVector2D_impl.h
+++ b/src/TNL/Containers/StaticVector2D_impl.h
@@ -97,6 +97,16 @@ StaticVector< 2, Real >& StaticVector< 2, Real >::operator *= ( const Real& c )
return *this;
}
+template< typename Real >
+__cuda_callable__
+StaticVector< 2, Real >& StaticVector< 2, Real >::operator /= ( const Real& c )
+{
+ const RealType d = 1.0 /c;
+ this->data[ 0 ] *= d;
+ this->data[ 1 ] *= d;
+ return *this;
+}
+
template< typename Real >
__cuda_callable__
StaticVector< 2, Real > StaticVector< 2, Real >::operator + ( const StaticVector& u ) const
diff --git a/src/TNL/Containers/StaticVector3D_impl.h b/src/TNL/Containers/StaticVector3D_impl.h
index d01e82077afd4feba62657d769803a68e1fe8fa1..655ed87365b1b53cd45d4c33805108ccc155cba5 100644
--- a/src/TNL/Containers/StaticVector3D_impl.h
+++ b/src/TNL/Containers/StaticVector3D_impl.h
@@ -101,6 +101,17 @@ StaticVector< 3, Real >& StaticVector< 3, Real >::operator *= ( const Real& c )
return *this;
}
+template< typename Real >
+__cuda_callable__
+StaticVector< 3, Real >& StaticVector< 3, Real >::operator /= ( const Real& c )
+{
+ const RealType d = 1.0 / c;
+ this->data[ 0 ] *= d;
+ this->data[ 1 ] *= d;
+ this->data[ 2 ] *= d;
+ return *this;
+}
+
template< typename Real >
__cuda_callable__
StaticVector< 3, Real > StaticVector< 3, Real >::operator + ( const StaticVector& u ) const
diff --git a/src/TNL/Containers/StaticVector_impl.h b/src/TNL/Containers/StaticVector_impl.h
index 45e680a098f36aeeed127295f210dcad5d629d28..5f28e277af22fb6bdd6beab0a58d8a158359c38c 100644
--- a/src/TNL/Containers/StaticVector_impl.h
+++ b/src/TNL/Containers/StaticVector_impl.h
@@ -91,6 +91,16 @@ StaticVector< Size, Real >& StaticVector< Size, Real >::operator *= ( const Real
return *this;
}
+template< int Size, typename Real >
+__cuda_callable__
+StaticVector< Size, Real >& StaticVector< Size, Real >::operator /= ( const Real& c )
+{
+ const RealType d = 1.0 / c;
+ for( int i = 0; i < Size; i++ )
+ this->data[ i ] *= d;
+ return *this;
+}
+
template< int Size, typename Real >
__cuda_callable__
StaticVector< Size, Real > StaticVector< Size, Real >::operator + ( const StaticVector& u ) const
diff --git a/src/TNL/Functions/CMakeLists.txt b/src/TNL/Functions/CMakeLists.txt
index 92359a15bf92d3255f9f1d97495cd7efc45725d0..161858ac8700cd1470c29d1510951320ed0b4739 100644
--- a/src/TNL/Functions/CMakeLists.txt
+++ b/src/TNL/Functions/CMakeLists.txt
@@ -14,8 +14,12 @@ SET( headers Domain.h
TestFunction.h
TestFunction_impl.h
VectorField.h
+ VectorFieldEvaluator.h
+ VectorFieldEvaluator_impl.h
VectorFieldGnuplotWriter.h
VectorFieldGnuplotWriter_impl.h
+ VectorFieldVTKWriter.h
+ VectorFieldVTKWriter_impl.h
)
SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/TNL/Functions )
diff --git a/src/TNL/Functions/MeshFunction_impl.h b/src/TNL/Functions/MeshFunction_impl.h
index 437313830a76e8df576249fe0f43472261202540..4e340e3a524673ca50e6b16579607820f2c5418c 100644
--- a/src/TNL/Functions/MeshFunction_impl.h
+++ b/src/TNL/Functions/MeshFunction_impl.h
@@ -68,14 +68,12 @@ MeshFunction( const MeshPointer& meshPointer,
const IndexType& offset )
//: meshPointer( meshPointer )
{
-
+ TNL_ASSERT_GE( data.getSize(), meshPointer->template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
+ "The input vector is not large enough for binding to the mesh function." );
SetupSynchronizer(meshPointer->GetDistMesh());
this->meshPointer=meshPointer;
this->data.bind( data, offset, getMesh().template getEntitiesCount< typename Mesh::template EntityType< MeshEntityDimension > >() );
- TNL_ASSERT( this->data.getSize() == this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
- std::cerr << "this->data.getSize() = " << this->data.getSize() << std::endl
- << "this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() = " << this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() );
}
@@ -89,14 +87,13 @@ MeshFunction( const MeshPointer& meshPointer,
const IndexType& offset )
//: meshPointer( meshPointer )
{
+ TNL_ASSERT_GE( data->getSize(), offset + meshPointer->template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
+ "The input vector is not large enough for binding to the mesh function." );
SetupSynchronizer(meshPointer->GetDistMesh());
this->meshPointer=meshPointer;
this->data.bind( *data, offset, getMesh().template getEntitiesCount< typename Mesh::template EntityType< MeshEntityDimension > >() );
- TNL_ASSERT( this->data.getSize() == this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
- std::cerr << "this->data.getSize() = " << this->data.getSize() << std::endl
- << "this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() = " << this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() );
}
template< typename Mesh,
@@ -207,14 +204,17 @@ bind( const MeshPointer& meshPointer,
const Vector& data,
const IndexType& offset )
{
+ TNL_ASSERT_GE( data.getSize(), offset + meshPointer->template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
+ "The input vector is not large enough for binding to the mesh function." );
+
+ SetupSynchronizer(meshPointer->GetDistMesh());
+
+ this->meshPointer = meshPointer;
- SetupSynchronizer(meshPointer->GetDistMesh());
+
this->meshPointer=meshPointer;
this->data.bind( data, offset, getMesh().template getEntitiesCount< typename Mesh::template EntityType< MeshEntityDimension > >() );
- TNL_ASSERT( this->data.getSize() == this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
- std::cerr << "this->data.getSize() = " << this->data.getSize() << std::endl
- << "this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() = " << this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() );
}
template< typename Mesh,
@@ -227,14 +227,15 @@ bind( const MeshPointer& meshPointer,
const SharedPointer< Vector >& data,
const IndexType& offset )
{
+ TNL_ASSERT_GE( data->getSize(), offset + meshPointer->template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
+ "The input vector is not large enough for binding to the mesh function." );
- SetupSynchronizer(meshPointer->GetDistMesh());
+ SetupSynchronizer(meshPointer->GetDistMesh());
+
+ this->meshPointer = meshPointer;
this->meshPointer=meshPointer;
this->data.bind( *data, offset, getMesh().template getEntitiesCount< typename Mesh::template EntityType< MeshEntityDimension > >() );
- TNL_ASSERT( this->data.getSize() == this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
- std::cerr << "this->data.getSize() = " << this->data.getSize() << std::endl
- << "this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() = " << this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() );
}
@@ -246,13 +247,10 @@ MeshFunction< Mesh, MeshEntityDimension, Real >::
setMesh( const MeshPointer& meshPointer )
{
- SetupSynchronizer(meshPointer->GetDistMesh());
+ SetupSynchronizer(meshPointer->GetDistMesh());
this->meshPointer=meshPointer;
this->data.setSize( getMesh().template getEntitiesCount< typename Mesh::template EntityType< MeshEntityDimension > >() );
- TNL_ASSERT( this->data.getSize() == this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >(),
- std::cerr << "this->data.getSize() = " << this->data.getSize() << std::endl
- << "this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() = " << this->getMesh().template getEntitiesCount< typename MeshType::template EntityType< MeshEntityDimension > >() );
}
template< typename Mesh,
diff --git a/src/TNL/Functions/VectorField.h b/src/TNL/Functions/VectorField.h
index cbc143251962455b074d426f87928c2a25dd762a..b4c3af83abad6ecdb50e86702b99a60db8018820 100644
--- a/src/TNL/Functions/VectorField.h
+++ b/src/TNL/Functions/VectorField.h
@@ -15,6 +15,7 @@
#include <TNL/Config/ParameterContainer.h>
#include <TNL/Functions/MeshFunction.h>
#include <TNL/Functions/VectorFieldGnuplotWriter.h>
+#include <TNL/Functions/VectorFieldVTKWriter.h>
namespace TNL {
namespace Functions {
@@ -90,6 +91,11 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
typedef typename MeshType::GlobalIndexType IndexType;
typedef VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, RealType > > ThisType;
typedef Containers::StaticVector< Size, RealType > VectorType;
+
+ static constexpr int getEntitiesDimension() { return FunctionType::getEntitiesDimension(); }
+
+ static constexpr int getMeshDimension() { return MeshType::getMeshDimension(); }
+
static void configSetup( Config::ConfigDescription& config,
const String& prefix = "" )
@@ -146,6 +152,10 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
return this->vectorField[ 0 ].template getData< Device >().template getMesh< Device >();
}
+ const MeshPointer& getMeshPointer() const
+ {
+ return this->vectorField[ 0 ]->getMeshPointer();
+ }
bool setup( const MeshPointer& meshPointer,
const Config::ParameterContainer& parameters,
@@ -164,6 +174,42 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
{
return Size * FunctionType::getDofs( meshPointer );
}
+
+ void bind( ThisType& vectorField )
+ {
+ for( int i = 0; i < Size; i ++ )
+ {
+ this->vectorField[ i ]->bind( vectorField[ i ] );
+ }
+ };
+
+ template< typename Vector >
+ void bind( const MeshPointer& meshPointer,
+ const Vector& data,
+ IndexType offset = 0 )
+ {
+ TNL_ASSERT_GE( data.getSize(), offset + Size * this->vectorField[ 0 ]->getDofs( meshPointer ),
+ "Attempt to bind vector which is not large enough." );
+ for( int i = 0; i < Size; i ++ )
+ {
+ this->vectorField[ i ].bind( meshPointer, data, offset );
+ offset += this->vectorField[ i ]->getDofs();
+ }
+ }
+
+ template< typename Vector >
+ void bind( const MeshPointer& meshPointer,
+ const SharedPointer< Vector >& dataPtr,
+ IndexType offset = 0 )
+ {
+ TNL_ASSERT_GE( dataPtr->getSize(), offset + Size * this->vectorField[ 0 ]->getDofs( meshPointer ),
+ "Attempt to bind vector which is not large enough." );
+ for( int i = 0; i < Size; i ++ )
+ {
+ this->vectorField[ i ]->bind( meshPointer, dataPtr, offset );
+ offset += this->vectorField[ i ]->getDofs( meshPointer );
+ }
+ }
__cuda_callable__
const FunctionPointer& operator[]( int i ) const
@@ -177,13 +223,29 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
return this->vectorField[ i ];
}
+ __cuda_callable__
+ void setElement( const IndexType i, const VectorType& v )
+ {
+ for( int j = 0; j < Size; j++ )
+ ( *this )[ j ]->getData().setElement( i, v[ j ] );
+ }
+
+ __cuda_callable__
+ VectorType getElement( const IndexType i ) const
+ {
+ VectorType v;
+ for( int j = 0; j < Size; j++ )
+ v[ j ] = ( *this )[ j ]->getData().getElement( i );
+ return v;
+ }
+
__cuda_callable__
VectorType getVector( const IndexType index ) const
{
VectorType v;
for( int i = 0; i < Size; i++ )
- // FIXME: the dereferencing operator of FunctionPointer is not __cuda_callable__
- v[ i ] = ( *this->vectorField[ i ] )[ index ];
+ // FIXME: fix the dereferencing operator in smart pointers to be __cuda_callable__
+ v[ i ] = this->vectorField[ i ].template getData< Devices::Cuda >()[ index ];
return v;
}
@@ -193,8 +255,8 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
{
VectorType v;
for( int i = 0; i < Size; i++ )
- // FIXME: the dereferencing operator of FunctionPointer is not __cuda_callable__
- v[ i ] = ( *this->vectorField[ i ] )( meshEntity );
+ // FIXME: fix the dereferencing operator in smart pointers to be __cuda_callable__
+ v[ i ] = this->vectorField[ i ].template getData< Devices::Cuda >()( meshEntity );
return v;
}
@@ -240,7 +302,7 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
return false;
}
if( format == "vtk" )
- return false; //MeshFunctionVTKWriter< ThisType >::write( *this, file );
+ return VectorFieldVTKWriter< ThisType >::write( *this, file, scale );
else if( format == "gnuplot" )
return VectorFieldGnuplotWriter< ThisType >::write( *this, file, scale );
else {
diff --git a/src/TNL/Functions/VectorFieldEvaluator.h b/src/TNL/Functions/VectorFieldEvaluator.h
new file mode 100644
index 0000000000000000000000000000000000000000..c1f6c0caa4360fe3c4e0c9d82ad25caef6904f3c
--- /dev/null
+++ b/src/TNL/Functions/VectorFieldEvaluator.h
@@ -0,0 +1,173 @@
+/***************************************************************************
+ VectorFieldEvaluator.h - description
+ -------------------
+ begin : Dec 24, 2017
+ copyright : (C) 2017 by Tomas Oberhuber et al.
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Functions/MeshFunction.h>
+#include <TNL/Functions/OperatorFunction.h>
+#include <TNL/Functions/FunctionAdapter.h>
+
+namespace TNL {
+namespace Functions {
+
+template< typename OutVectorField,
+ typename InVectorField,
+ typename Real >
+class VectorFieldEvaluatorTraverserUserData
+{
+ public:
+ typedef InVectorField InVectorFieldType;
+
+ VectorFieldEvaluatorTraverserUserData( const InVectorField* inVectorField,
+ const Real& time,
+ OutVectorField* outVectorField,
+ const Real& outVectorFieldMultiplicator,
+ const Real& inVectorFieldMultiplicator )
+ : outVectorField( outVectorField ),
+ inVectorField( inVectorField ),
+ time( time ),
+ outVectorFieldMultiplicator( outVectorFieldMultiplicator ),
+ inVectorFieldMultiplicator( inVectorFieldMultiplicator )
+ {}
+
+ OutVectorField* outVectorField;
+ const InVectorField* inVectorField;
+ const Real time, outVectorFieldMultiplicator, inVectorFieldMultiplicator;
+
+};
+
+
+/***
+ * General vector field evaluator. As an input function any type implementing
+ * getValue( meshEntity, time ) may be substituted.
+ * Methods:
+ * evaluate() - evaluate the input function on its domain
+ * evaluateAllEntities() - evaluate the input function on ALL mesh entities
+ * evaluateInteriorEntities() - evaluate the input function only on the INTERIOR mesh entities
+ * evaluateBoundaryEntities() - evaluate the input function only on the BOUNDARY mesh entities
+ */
+template< typename OutVectorField,
+ typename InVectorField >
+class VectorFieldEvaluator
+{
+ static_assert( OutVectorField::getDomainDimension() == InVectorField::getDomainDimension(),
+ "Input and output vector field must have the same domain dimensions." );
+
+ public:
+ typedef typename OutVectorField::RealType RealType;
+ typedef typename OutVectorField::MeshType MeshType;
+ typedef typename MeshType::DeviceType DeviceType;
+ typedef Functions::VectorFieldEvaluatorTraverserUserData< OutVectorField, InVectorField, RealType > TraverserUserData;
+
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+ static void evaluate( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time = 0.0,
+ const RealType& outFunctionMultiplicator = 0.0,
+ const RealType& inFunctionMultiplicator = 1.0 );
+
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+ static void evaluateAllEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time = 0.0,
+ const RealType& outFunctionMultiplicator = 0.0,
+ const RealType& inFunctionMultiplicator = 1.0 );
+
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+ static void evaluateInteriorEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time = 0.0,
+ const RealType& outFunctionMultiplicator = 0.0,
+ const RealType& inFunctionMultiplicator = 1.0 );
+
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+ static void evaluateBoundaryEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time = 0.0,
+ const RealType& outFunctionMultiplicator = 0.0,
+ const RealType& inFunctionMultiplicator = 1.0 );
+
+ protected:
+
+ enum EntitiesType { all, boundary, interior };
+
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+ static void evaluateEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time,
+ const RealType& outFunctionMultiplicator,
+ const RealType& inFunctionMultiplicator,
+ EntitiesType entitiesType );
+
+
+};
+
+
+template< typename MeshType,
+ typename UserData >
+class VectorFieldEvaluatorAssignmentEntitiesProcessor
+{
+ public:
+
+ template< typename EntityType >
+ __cuda_callable__
+ static inline void processEntity( const MeshType& mesh,
+ UserData& userData,
+ const EntityType& entity )
+ {
+ userData.outVectorField->setElement(
+ entity.getIndex(),
+ userData.inVectorFieldMultiplicator * ( *userData.inVectorField )( entity, userData.time ) );
+
+ /*typedef FunctionAdapter< MeshType, typename UserData::InVectorFieldType > FunctionAdapter;
+ ( *userData.meshFunction )( entity ) =
+ userData.inFunctionMultiplicator *
+ FunctionAdapter::getValue( *userData.function, entity, userData.time );*/
+ /*cerr << "Idx = " << entity.getIndex()
+ << " Value = " << FunctionAdapter::getValue( *userData.function, entity, userData.time )
+ << " stored value = " << ( *userData.meshFunction )( entity )
+ << " multiplicators = " << std::endl;*/
+ }
+};
+
+template< typename MeshType,
+ typename UserData >
+class VectorFieldEvaluatorAdditionEntitiesProcessor
+{
+ public:
+
+ template< typename EntityType >
+ __cuda_callable__
+ static inline void processEntity( const MeshType& mesh,
+ UserData& userData,
+ const EntityType& entity )
+ {
+ const auto& i = entity.getIndex();
+ const auto v = userData.outVectorFieldMultiplicator * userData.outVectorField->getElement( i );
+ userData.outVectorField->setElement(
+ i,
+ v + userData.inVectorFieldMultiplicator * ( *userData.inVectorField )( entity, userData.time ) );
+
+ /*typedef FunctionAdapter< MeshType, typename UserData::InVectorFieldType > FunctionAdapter;
+ ( *userData.meshFunction )( entity ) =
+ userData.outFunctionMultiplicator * ( *userData.meshFunction )( entity ) +
+ userData.inFunctionMultiplicator *
+ FunctionAdapter::getValue( *userData.function, entity, userData.time );*/
+ /*cerr << "Idx = " << entity.getIndex()
+ << " Value = " << FunctionAdapter::getValue( *userData.function, entity, userData.time )
+ << " stored value = " << ( *userData.meshFunction )( entity )
+ << " multiplicators = " << std::endl;*/
+ }
+};
+
+} // namespace Functions
+} // namespace TNL
+
+#include <TNL/Functions/VectorFieldEvaluator_impl.h>
diff --git a/src/TNL/Functions/VectorFieldEvaluator_impl.h b/src/TNL/Functions/VectorFieldEvaluator_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..56668ef24f5d6a3d8163382a6b251b439cdd22f8
--- /dev/null
+++ b/src/TNL/Functions/VectorFieldEvaluator_impl.h
@@ -0,0 +1,173 @@
+/***************************************************************************
+ VectorFieldEvaluator.h - description
+ -------------------
+ begin : Dec 5, 2017
+ copyright : (C) 2017 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Functions/VectorFieldEvaluator.h>
+#include <TNL/Meshes/Traverser.h>
+
+namespace TNL {
+namespace Functions {
+
+template< typename OutVectorField,
+ typename InVectorField >
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+void
+VectorFieldEvaluator< OutVectorField, InVectorField >::
+evaluate( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time,
+ const RealType& outFunctionMultiplicator,
+ const RealType& inFunctionMultiplicator )
+{
+ static_assert( std::is_same< typename std::decay< typename OutVectorFieldPointer::ObjectType >::type, OutVectorField >::value, "expected a smart pointer" );
+ static_assert( std::is_same< typename std::decay< typename InVectorFieldPointer::ObjectType >::type, InVectorField >::value, "expected a smart pointer" );
+
+ switch( InVectorField::getDomainType() )
+ {
+ case NonspaceDomain:
+ case SpaceDomain:
+ case MeshDomain:
+ evaluateEntities( meshFunction, function, time, outFunctionMultiplicator, inFunctionMultiplicator, all );
+ break;
+ case MeshInteriorDomain:
+ evaluateEntities( meshFunction, function, time, outFunctionMultiplicator, inFunctionMultiplicator, interior );
+ break;
+ case MeshBoundaryDomain:
+ evaluateEntities( meshFunction, function, time, outFunctionMultiplicator, inFunctionMultiplicator, boundary );
+ break;
+ }
+}
+
+
+template< typename OutVectorField,
+ typename InVectorField >
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+void
+VectorFieldEvaluator< OutVectorField, InVectorField >::
+evaluateAllEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time,
+ const RealType& outFunctionMultiplicator,
+ const RealType& inFunctionMultiplicator )
+{
+ static_assert( std::is_same< typename std::decay< typename OutVectorFieldPointer::ObjectType >::type, OutVectorField >::value, "expected a smart pointer" );
+ static_assert( std::is_same< typename std::decay< typename InVectorFieldPointer::ObjectType >::type, InVectorField >::value, "expected a smart pointer" );
+
+ return evaluateEntities( meshFunction, function, time, outFunctionMultiplicator, inFunctionMultiplicator, all );
+}
+
+template< typename OutVectorField,
+ typename InVectorField >
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+void
+VectorFieldEvaluator< OutVectorField, InVectorField >::
+evaluateInteriorEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time,
+ const RealType& outFunctionMultiplicator,
+ const RealType& inFunctionMultiplicator )
+{
+ static_assert( std::is_same< typename std::decay< typename OutVectorFieldPointer::ObjectType >::type, OutVectorField >::value, "expected a smart pointer" );
+ static_assert( std::is_same< typename std::decay< typename InVectorFieldPointer::ObjectType >::type, InVectorField >::value, "expected a smart pointer" );
+
+ return evaluateEntities( meshFunction, function, time, outFunctionMultiplicator, inFunctionMultiplicator, interior );
+}
+
+template< typename OutVectorField,
+ typename InVectorField >
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+void
+VectorFieldEvaluator< OutVectorField, InVectorField >::
+evaluateBoundaryEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time,
+ const RealType& outFunctionMultiplicator,
+ const RealType& inFunctionMultiplicator )
+{
+ static_assert( std::is_same< typename std::decay< typename OutVectorFieldPointer::ObjectType >::type, OutVectorField >::value, "expected a smart pointer" );
+ static_assert( std::is_same< typename std::decay< typename InVectorFieldPointer::ObjectType >::type, InVectorField >::value, "expected a smart pointer" );
+
+ return evaluateEntities( meshFunction, function, time, outFunctionMultiplicator, inFunctionMultiplicator, boundary );
+}
+
+
+
+template< typename OutVectorField,
+ typename InVectorField >
+ template< typename OutVectorFieldPointer, typename InVectorFieldPointer >
+void
+VectorFieldEvaluator< OutVectorField, InVectorField >::
+evaluateEntities( OutVectorFieldPointer& meshFunction,
+ const InVectorFieldPointer& function,
+ const RealType& time,
+ const RealType& outFunctionMultiplicator,
+ const RealType& inFunctionMultiplicator,
+ EntitiesType entitiesType )
+{
+ static_assert( std::is_same< typename std::decay< typename OutVectorFieldPointer::ObjectType >::type, OutVectorField >::value, "expected a smart pointer" );
+ static_assert( std::is_same< typename std::decay< typename InVectorFieldPointer::ObjectType >::type, InVectorField >::value, "expected a smart pointer" );
+
+ typedef typename MeshType::template EntityType< OutVectorField::getEntitiesDimension() > MeshEntityType;
+ typedef Functions::VectorFieldEvaluatorAssignmentEntitiesProcessor< MeshType, TraverserUserData > AssignmentEntitiesProcessor;
+ typedef Functions::VectorFieldEvaluatorAdditionEntitiesProcessor< MeshType, TraverserUserData > AdditionEntitiesProcessor;
+ //typedef typename OutVectorField::MeshPointer OutMeshPointer;
+ typedef SharedPointer< TraverserUserData, DeviceType > TraverserUserDataPointer;
+
+ SharedPointer< TraverserUserData, DeviceType >
+ userData( &function.template getData< DeviceType >(),
+ time,
+ &meshFunction.template modifyData< DeviceType >(),
+ outFunctionMultiplicator,
+ inFunctionMultiplicator );
+ Meshes::Traverser< MeshType, MeshEntityType > meshTraverser;
+ switch( entitiesType )
+ {
+ case all:
+ if( outFunctionMultiplicator )
+ meshTraverser.template processAllEntities< TraverserUserData,
+ AdditionEntitiesProcessor >
+ ( meshFunction->getMeshPointer(),
+ userData );
+ else
+ meshTraverser.template processAllEntities< TraverserUserData,
+ AssignmentEntitiesProcessor >
+ ( meshFunction->getMeshPointer(),
+ userData );
+ break;
+ case interior:
+ if( outFunctionMultiplicator )
+ meshTraverser.template processInteriorEntities< TraverserUserData,
+ AdditionEntitiesProcessor >
+ ( meshFunction->getMeshPointer(),
+ userData );
+ else
+ meshTraverser.template processInteriorEntities< TraverserUserData,
+ AssignmentEntitiesProcessor >
+ ( meshFunction->getMeshPointer(),
+ userData );
+ break;
+ case boundary:
+ if( outFunctionMultiplicator )
+ meshTraverser.template processBoundaryEntities< TraverserUserData,
+ AdditionEntitiesProcessor >
+ ( meshFunction->getMeshPointer(),
+ userData );
+ else
+ meshTraverser.template processBoundaryEntities< TraverserUserData,
+ AssignmentEntitiesProcessor >
+ ( meshFunction->getMeshPointer(),
+ userData );
+ break;
+ }
+}
+
+} // namespace Functions
+} // namespace TNL
diff --git a/src/TNL/Functions/VectorFieldVTKWriter.h b/src/TNL/Functions/VectorFieldVTKWriter.h
new file mode 100644
index 0000000000000000000000000000000000000000..6d8b1a8535b25e076c83e688706f37490086c39d
--- /dev/null
+++ b/src/TNL/Functions/VectorFieldVTKWriter.h
@@ -0,0 +1,264 @@
+/***************************************************************************
+ VectorFieldVTKWriter.h - description
+ -------------------
+ begin : Jan 10, 2018
+ copyright : (C) 2018 by oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Meshes/Grid.h>
+
+namespace TNL {
+namespace Functions {
+
+template< int, typename > class VectorField;
+
+template< typename VectorField >
+class VectorFieldVTKWriter
+{
+ public:
+
+ static bool write( const VectorField& vectorField,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorField& vectorField,
+ std::ostream& str ){}
+
+};
+
+/***
+ * 1D grids cells
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+/***
+ * 1D grids vertices
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+
+/***
+ * 2D grids cells
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+/***
+ * 2D grids faces
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+/***
+ * 2D grids vertices
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+
+/***
+ * 3D grids cells
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 3, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+/***
+ * 3D grids faces
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+/***
+ * 3D grids edges
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+/***
+ * 3D grids vertices
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >
+{
+ public:
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
+ using VectorType = typename VectorFieldType::VectorType;
+
+ static bool write( const VectorFieldType& function,
+ std::ostream& str,
+ const double& scale );
+
+ static void writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str );
+
+};
+
+} // namespace Functions
+} // namespace TNL
+
+#include <TNL/Functions/VectorFieldVTKWriter_impl.h>
diff --git a/src/TNL/Functions/VectorFieldVTKWriter_impl.h b/src/TNL/Functions/VectorFieldVTKWriter_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..938227d22b57f61d2f6c4d5f4b7b13a9044d3aa8
--- /dev/null
+++ b/src/TNL/Functions/VectorFieldVTKWriter_impl.h
@@ -0,0 +1,881 @@
+/***************************************************************************
+ VectorFieldVTKWriter_impl.h - description
+ -------------------
+ begin : Jan 10, 2018
+ copyright : (C) 2018 by oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Functions/VectorFieldVTKWriter.h>
+#include <TNL/Functions/VectorField.h>
+
+namespace TNL {
+namespace Functions {
+
+template< typename VectorField >
+bool
+VectorFieldVTKWriter< VectorField >::
+write( const VectorField& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ std::cerr << "VTK writer for vector field defined on mesh type " << VectorField::MeshType::getType() << " is not (yet) implemented." << std::endl;
+ return false;
+}
+
+/****
+ * 1D grid, cells
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType origin = mesh.getOrigin().x();
+ const RealType spaceStep = mesh.getSpaceSteps().x();
+
+ str << "POINTS " << mesh.getDimensions().x() + 1 << " float" << std::endl;
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << origin + i * spaceStep << " 0 0" << std::endl;
+ }
+
+ str << std::endl << "CELLS " << mesh.getDimensions().x() << " " << mesh.getDimensions().x() * 3 << std::endl;
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "2 " << i << " " << i+1 << std::endl;
+ }
+
+ str << std::endl << "CELL_TYPES " << mesh.getDimensions().x() << std::endl;
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "3 " << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << mesh.getDimensions().x() << std::endl;
+ str << "VECTORS cellVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < mesh.template getEntitiesCount< typename MeshType::Cell >(); i++ )
+ {
+ typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 1D grid, vertices
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType origin = mesh.getOrigin().x();
+ const RealType spaceStep = mesh.getSpaceSteps().x();
+
+ str << "POINTS " << mesh.getDimensions().x() + 1 << " float" << std::endl;
+ for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
+ {
+ str << origin + i * spaceStep << " 0 0" << std::endl;
+ }
+
+ str << std::endl << "CELLS " << mesh.getDimensions().x() + 1 << " " << ( mesh.getDimensions().x() + 1 ) * 2 << std::endl;
+ for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
+ {
+ str << "1 " << i << std::endl;
+ }
+
+ str << std::endl << "CELL_TYPES " << mesh.getDimensions().x() + 1 << std::endl;
+ for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
+ {
+ str << "1 " << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << mesh.getDimensions().x() + 1 << std::endl;
+ str << "VECTORS VerticesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < mesh.template getEntitiesCount< typename MeshType::Vertex >(); i++ )
+ {
+ typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 2D grid, cells
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+ const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Cell >();
+
+ str << "POINTS " << verticesCount << " " << getType< RealType >() << std::endl;
+ for (int j = 0; j < mesh.getDimensions().y() + 1; j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " 0" << std::endl;
+ }
+ }
+
+ str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 5 << std::endl;
+ for (int j = 0; j < mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "4 " << j * ( mesh.getDimensions().x() + 1 ) + i << " " << j * ( mesh.getDimensions().x() + 1 )+ i + 1 <<
+ " " << (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " " << (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << mesh.getDimensions().x() * mesh.getDimensions().y() << std::endl;
+ for (int i = 0; i < mesh.getDimensions().x()*mesh.getDimensions().y(); i++)
+ {
+ str << "8 " << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+ str << "VECTORS cellVectorFieldValues " << getType< RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < entitiesCount; i++ )
+ {
+ typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ {
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 );
+ if( i < 2 )
+ str << " ";
+ }
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 2D grid, faces
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ typedef typename MeshType::template EntityType< 0 > Vertex;
+ typedef typename MeshType::template EntityType< 1 > Face;
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+ const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Face >();
+
+ str << "POINTS " << verticesCount << " float" << std::endl;
+ for (int j = 0; j < ( mesh.getDimensions().y() + 1); j++)
+ {
+ for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " 0" << std::endl;
+ }
+ }
+
+ str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 3 << std::endl;
+ for (int j = 0; j < mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
+ {
+ str << "2 " << j * ( mesh.getDimensions().x() + 1 ) + i << " " << (j+1) * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
+ }
+ }
+
+ for (int j = 0; j < (mesh.getDimensions().y()+1); j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "2 " << j * ( mesh.getDimensions().x() + 1 ) + i << " " <<j * ( mesh.getDimensions().x() + 1 ) + i + 1<< std::endl;
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
+ for (int i = 0; i < entitiesCount; i++)
+ {
+ str << "3" << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+ str << "VECTORS FaceslVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < entitiesCount; i++ )
+ {
+ typename MeshType::Face entity = mesh.template getEntity< typename MeshType::Face >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 2D grid, vertices
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ typedef typename MeshType::template EntityType< 0 > Vertex;
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+
+ str << "POINTS " << verticesCount << " float" << std::endl;
+ for (int j = 0; j < ( mesh.getDimensions().y() + 1); j++)
+ {
+ for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " 0" << std::endl;
+ }
+ }
+
+ str << std::endl << "CELLS " << verticesCount << " " << verticesCount * 2 << std::endl;
+ for (int j = 0; j < ( mesh.getDimensions().y() + 1 ); j++)
+ {
+ for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
+ {
+ str << "1 " << j * mesh.getDimensions().x() + i << std::endl;
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << verticesCount << std::endl;
+ for (int i = 0; i < verticesCount; i++)
+ {
+ str << "1" << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << verticesCount << std::endl;
+ str << "VECTORS VerticesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < verticesCount; i++ )
+ {
+ typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 3D grid, cells
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const RealType originZ = mesh.getOrigin().z();
+ const RealType spaceStepZ = mesh.getSpaceSteps().z();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+ const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Cell >();
+
+ str << "POINTS " << verticesCount << " float" << std::endl;
+ for (int k = 0; k <= mesh.getDimensions().y(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
+ originZ + k * spaceStepZ << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELLS " << entitiesCount << " " <<
+ entitiesCount * 9 << std::endl;
+ for (int k = 0; k < mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j < mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "8 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
+ for (int i = 0; i < entitiesCount; i++)
+ {
+ str << "11" << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+ str << "VECTORS cellVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < entitiesCount; i++ )
+ {
+ typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 3D grid, faces
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const RealType originZ = mesh.getOrigin().z();
+ const RealType spaceStepZ = mesh.getSpaceSteps().z();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+ const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Face >();
+
+ str << "POINTS " << verticesCount << " float" << std::endl;
+ for (int k = 0; k <= mesh.getDimensions().y(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
+ originZ + k * spaceStepZ << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 5 << std::endl;
+ for (int k = 0; k < mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j < mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << "4 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
+ }
+ }
+ }
+
+ for (int k = 0; k < mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "4 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
+ }
+ }
+ }
+
+ for (int k = 0; k <= mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j < mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "4 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1<< std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
+ for (int i = 0; i < entitiesCount; i++)
+ {
+ str << "8" << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+ str << "VECTORS facesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < entitiesCount; i++ )
+ {
+ typename MeshType::Face entity = mesh.template getEntity< typename MeshType::Face >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 3D grid, edges
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const RealType originZ = mesh.getOrigin().z();
+ const RealType spaceStepZ = mesh.getSpaceSteps().z();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+ const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Edge >();
+
+ str << "POINTS " << verticesCount << " float" << std::endl;
+ for (int k = 0; k <= mesh.getDimensions().y(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
+ originZ + k * spaceStepZ << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 3 << std::endl;
+ for (int k = 0; k <= mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i < mesh.getDimensions().x(); i++)
+ {
+ str << "3 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
+ }
+ }
+ }
+
+ for (int k = 0; k <= mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j < mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << "3 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
+ }
+ }
+ }
+
+ for (int k = 0; k < mesh.getDimensions().z(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << "3 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
+ << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
+ for (int i = 0; i < entitiesCount; i++)
+ {
+ str << "3" << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+ str << "VECTORS edgesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < entitiesCount; i++ )
+ {
+ typename MeshType::Edge entity = mesh.template getEntity< typename MeshType::Edge >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+/****
+ * 3D grid, vertices
+ */
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+void
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >::
+writeHeader( const VectorFieldType& vectorField,
+ std::ostream& str )
+{
+ const MeshType& mesh = vectorField.getMesh();
+ const typename MeshType::PointType& origin = mesh.getOrigin();
+ const typename MeshType::PointType& proportions = mesh.getProportions();
+ str << "# vtk DataFile Version 2.0" << std::endl;
+ str << "TNL DATA" << std::endl;
+ str << "ASCII" << std::endl;
+ str << "DATASET UNSTRUCTURED_GRID" << std::endl;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ int VectorFieldSize >
+bool
+VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >::
+write( const VectorFieldType& vectorField,
+ std::ostream& str,
+ const double& scale )
+{
+ writeHeader(vectorField, str);
+
+ const MeshType& mesh = vectorField.getMesh();
+ const RealType originX = mesh.getOrigin().x();
+ const RealType spaceStepX = mesh.getSpaceSteps().x();
+ const RealType originY = mesh.getOrigin().y();
+ const RealType spaceStepY = mesh.getSpaceSteps().y();
+ const RealType originZ = mesh.getOrigin().z();
+ const RealType spaceStepZ = mesh.getSpaceSteps().z();
+ const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
+
+ str << "POINTS " << verticesCount << " float" << std::endl;
+ for (int k = 0; k <= mesh.getDimensions().y(); k++)
+ {
+ for (int j = 0; j <= mesh.getDimensions().y(); j++)
+ {
+ for (int i = 0; i <= mesh.getDimensions().x(); i++)
+ {
+ str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
+ originZ + k * spaceStepZ << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELLS " << verticesCount << " " << verticesCount * 2 << std::endl;
+ for (int k = 0; k < ( mesh.getDimensions().z() + 1 ); k++)
+ {
+ for (int j = 0; j < ( mesh.getDimensions().y() + 1 ); j++)
+ {
+ for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
+ {
+ str << "1 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
+ }
+ }
+ }
+
+ str << std::endl << "CELL_TYPES " << verticesCount << std::endl;
+ for (int i = 0; i < verticesCount; i++)
+ {
+ str << "1" << std::endl;
+ }
+
+ str << std::endl << "CELL_DATA " << verticesCount << std::endl;
+ str << "VECTORS verticesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
+
+ for( MeshIndex i = 0; i < verticesCount; i++ )
+ {
+ typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
+ entity.refresh();
+ const VectorType v = vectorField.getElement( entity.getIndex() );
+ for( int i = 0; i < 3; i++ )
+ str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
+ str << std::endl;
+ }
+
+ return true;
+}
+
+} // namespace Functions
+} // namespace TNL
diff --git a/src/TNL/Math.h b/src/TNL/Math.h
index 0a40b248df9906936f86b9eb77339ddecfaeeb5a..1f754bcecaa88ffc973387b56aa9d692aaca344c 100644
--- a/src/TNL/Math.h
+++ b/src/TNL/Math.h
@@ -27,12 +27,17 @@ template< typename T1, typename T2, typename ResultType = typename std::common_t
__cuda_callable__ inline
ResultType min( const T1& a, const T2& b )
{
-#if defined(__CUDA_ARCH__)
+#if __cplusplus >= 201402L
+ // std::min is constexpr since C++14 so it can be reused directly
+ return std::min( (ResultType) a, (ResultType) b );
+#else
+ #if defined(__CUDA_ARCH__)
return ::min( (ResultType) a, (ResultType) b );
-#elif defined(__MIC__)
+ #elif defined(__MIC__)
return a < b ? a : b;
-#else
+ #else
return std::min( (ResultType) a, (ResultType) b );
+ #endif
#endif
}
@@ -46,12 +51,17 @@ template< typename T1, typename T2, typename ResultType = typename std::common_t
__cuda_callable__
ResultType max( const T1& a, const T2& b )
{
-#if defined(__CUDA_ARCH__)
+#if __cplusplus >= 201402L
+ // std::max is constexpr since C++14 so it can be reused directly
+ return std::max( (ResultType) a, (ResultType) b );
+#else
+ #if defined(__CUDA_ARCH__)
return ::max( (ResultType) a, (ResultType) b );
-#elif defined(__MIC__)
+ #elif defined(__MIC__)
return a > b ? a : b;
-#else
+ #else
return std::max( (ResultType) a, (ResultType) b );
+ #endif
#endif
}
diff --git a/src/TNL/ParallelFor.h b/src/TNL/ParallelFor.h
index ec66906a5e3c3de9930cc622e150dacccbe4d5f8..54a3a86412636fab14f0cff2ee84319a40d96a70 100644
--- a/src/TNL/ParallelFor.h
+++ b/src/TNL/ParallelFor.h
@@ -193,7 +193,7 @@ struct ParallelFor2D< Devices::Cuda >
dim3 blockSize;
if( sizeX >= sizeY * sizeY ) {
- blockSize.x = TNL::min( 256, sizeX );
+ blockSize.x = TNL::min( 128, sizeX );
blockSize.y = 1;
}
else if( sizeY >= sizeX * sizeX ) {
@@ -212,7 +212,7 @@ struct ParallelFor2D< Devices::Cuda >
gridCount.x = Devices::Cuda::getNumberOfGrids( sizeX );
gridCount.y = Devices::Cuda::getNumberOfGrids( sizeY );
- Devices::Cuda::synchronizeDevice();
+// Devices::Cuda::synchronizeDevice();
if( gridCount.x == 1 && gridCount.y == 1 )
ParallelFor2DKernel< false, false ><<< gridSize, blockSize >>>
@@ -281,7 +281,7 @@ struct ParallelFor3D< Devices::Cuda >
gridCount.y = Devices::Cuda::getNumberOfGrids( sizeY );
gridCount.z = Devices::Cuda::getNumberOfGrids( sizeZ );
- Devices::Cuda::synchronizeDevice();
+// Devices::Cuda::synchronizeDevice();
if( gridCount.x == 1 && gridCount.y == 1 && gridCount.z == 1 )
ParallelFor3DKernel< false, false, false ><<< gridSize, blockSize >>>
diff --git a/src/TNL/SharedPointer.h b/src/TNL/SharedPointer.h
index c9e71427518d373f26aec6ff57e386f5809d7365..2030e4898c196282b1310352405fc1a35218e69b 100644
--- a/src/TNL/SharedPointer.h
+++ b/src/TNL/SharedPointer.h
@@ -238,6 +238,11 @@ class SharedPointer< Object, Devices::Host > : public SmartPointer
{
this->free();
}
+
+ void swap( ThisType& ptr2 )
+ {
+ std::swap( this->pd, ptr2.pd );
+ }
~SharedPointer()
{
@@ -531,7 +536,13 @@ class SharedPointer< Object, Devices::Cuda > : public SmartPointer
void clear()
{
this->free();
- }
+ }
+
+ void swap( ThisType& ptr2 )
+ {
+ std::swap( this->pd, ptr2.pd );
+ std::swap( this->cuda_pointer, ptr2.cuda_pointer );
+ }
~SharedPointer()
{
@@ -854,6 +865,12 @@ class SharedPointer< Object, Devices::MIC > : public SmartPointer
this->free();
}
+ void swap( ThisType& ptr2 )
+ {
+ std::swap( this->pd, ptr2.pd );
+ std::swap( this->mic_pointer, ptr2.mic_pointer );
+ }
+
~SharedPointer()
{
this->free();
diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h
index 3180abb3c34641f7359710dd4e4d8fde44578e77..f29acee8fc08c380fa8b64e8237f95a92d4e71f7 100644
--- a/src/TNL/Solvers/IterativeSolverMonitor.h
+++ b/src/TNL/Solvers/IterativeSolverMonitor.h
@@ -35,6 +35,8 @@ public:
void setResidue( const RealType& residue );
void setVerbose( const Index& verbose );
+
+ void setNodesPerIteration( const IndexType& nodes );
virtual void refresh();
@@ -45,11 +47,13 @@ protected:
std::atomic_bool saved;
- RealType time, saved_time, timeStep, saved_timeStep, residue, saved_residue;
+ RealType time, saved_time, timeStep, saved_timeStep, residue, saved_residue, elapsed_time_before_refresh;
- IndexType iterations, saved_iterations;
+ IndexType iterations, saved_iterations, iterations_before_refresh;
IndexType verbose;
+
+ IndexType nodesPerIteration;
};
} // namespace Solvers
diff --git a/src/TNL/Solvers/IterativeSolverMonitor_impl.h b/src/TNL/Solvers/IterativeSolverMonitor_impl.h
index ef18320193dc83007fc87bd4e04e960f95c0ff0e..a94c033568a96e723fe4ba76c8b452f5621f5f9b 100644
--- a/src/TNL/Solvers/IterativeSolverMonitor_impl.h
+++ b/src/TNL/Solvers/IterativeSolverMonitor_impl.h
@@ -37,9 +37,12 @@ IterativeSolverMonitor< Real, Index > :: IterativeSolverMonitor()
saved_timeStep( 0.0 ),
residue( 0.0 ),
saved_residue( 0.0 ),
+ elapsed_time_before_refresh( 0.0 ),
iterations( 0 ),
saved_iterations( 0 ),
- verbose( 1 )
+ iterations_before_refresh( 0 ),
+ verbose( 1 ),
+ nodesPerIteration( 0 )
{
}
@@ -93,6 +96,12 @@ void IterativeSolverMonitor< Real, Index > :: setVerbose( const Index& verbose )
this->verbose = verbose;
}
+template< typename Real, typename Index>
+void IterativeSolverMonitor< Real, Index > :: setNodesPerIteration( const IndexType& nodes )
+{
+ this->nodesPerIteration = nodes;
+}
+
template< typename Real, typename Index>
void IterativeSolverMonitor< Real, Index > :: refresh()
{
@@ -133,7 +142,7 @@ void IterativeSolverMonitor< Real, Index > :: refresh()
if( (saved) ? saved_timeStep : timeStep > 0 ) {
// print_item( " TAU:" );
print_item( " TAU:", 0 );
- print_item( real_to_string( (saved) ? saved_timeStep : timeStep, 5 ), 8 );
+ print_item( real_to_string( (saved) ? saved_timeStep : timeStep, 5 ), 10 );
}
const std::string displayed_stage = (saved) ? saved_stage : stage;
@@ -159,6 +168,14 @@ void IterativeSolverMonitor< Real, Index > :: refresh()
print_item( real_to_string( (saved) ? saved_residue : residue, 5 ), 12 );
}
+ if( nodesPerIteration ) {
+ const RealType mlups = nodesPerIteration * (iterations - iterations_before_refresh) / (getElapsedTime() - elapsed_time_before_refresh) * 1e-6;
+ print_item( " MLUPS:", 0 );
+ print_item( real_to_string( mlups, 5 ), 7 );
+ }
+ iterations_before_refresh = iterations;
+ elapsed_time_before_refresh = getElapsedTime();
+
// return to the beginning of the line
std::cout << "\r" << std::flush;
}
diff --git a/src/TNL/Solvers/ODE/Euler_impl.h b/src/TNL/Solvers/ODE/Euler_impl.h
index 781379ecf3c27afdbb70224f772a52398ba85df5..92cb1e5c974b6a4b19c34b3d116f4fcb6d3ec568 100644
--- a/src/TNL/Solvers/ODE/Euler_impl.h
+++ b/src/TNL/Solvers/ODE/Euler_impl.h
@@ -136,7 +136,7 @@ bool Euler< Problem > :: solve( DofVectorPointer& u )
currentTau = this -> getStopTime() - time; //we don't want to keep such tau
else this -> tau = currentTau;
- this -> refreshSolverMonitor();
+ this->refreshSolverMonitor();
/****
* Check stop conditions.
@@ -158,8 +158,8 @@ bool Euler< Problem > :: solve( DofVectorPointer& u )
template< typename Problem >
void Euler< Problem > :: computeNewTimeLevel( DofVectorPointer& u,
- RealType tau,
- RealType& currentResidue )
+ RealType tau,
+ RealType& currentResidue )
{
RealType localResidue = RealType( 0.0 );
const IndexType size = k1->getSize();
diff --git a/src/TNL/Solvers/PDE/ExplicitTimeStepper.h b/src/TNL/Solvers/PDE/ExplicitTimeStepper.h
index a4ae5927f9eb5ec68e33ca0d41d0198db25c957d..32fbe5e933868216f08b1e648d4b320c53dd1128 100644
--- a/src/TNL/Solvers/PDE/ExplicitTimeStepper.h
+++ b/src/TNL/Solvers/PDE/ExplicitTimeStepper.h
@@ -27,78 +27,78 @@ class ExplicitTimeStepper
{
public:
- typedef Problem ProblemType;
- typedef OdeSolver< ExplicitTimeStepper< Problem, OdeSolver > > OdeSolverType;
- typedef typename Problem::RealType RealType;
- typedef typename Problem::DeviceType DeviceType;
- typedef typename Problem::IndexType IndexType;
- typedef typename Problem::MeshType MeshType;
- typedef SharedPointer< MeshType > MeshPointer;
- typedef typename ProblemType::DofVectorType DofVectorType;
- typedef typename ProblemType::MeshDependentDataType MeshDependentDataType;
- typedef SharedPointer< DofVectorType, DeviceType > DofVectorPointer;
- typedef SharedPointer< MeshDependentDataType, DeviceType > MeshDependentDataPointer;
- typedef IterativeSolverMonitor< RealType, IndexType > SolverMonitorType;
-
- static_assert( ProblemType::isTimeDependent(), "The problem is not time dependent." );
+ typedef Problem ProblemType;
+ typedef OdeSolver< ExplicitTimeStepper< Problem, OdeSolver > > OdeSolverType;
+ typedef typename Problem::RealType RealType;
+ typedef typename Problem::DeviceType DeviceType;
+ typedef typename Problem::IndexType IndexType;
+ typedef typename Problem::MeshType MeshType;
+ typedef SharedPointer< MeshType > MeshPointer;
+ typedef typename ProblemType::DofVectorType DofVectorType;
+ typedef typename ProblemType::MeshDependentDataType MeshDependentDataType;
+ typedef SharedPointer< DofVectorType, DeviceType > DofVectorPointer;
+ typedef SharedPointer< MeshDependentDataType, DeviceType > MeshDependentDataPointer;
+ typedef IterativeSolverMonitor< RealType, IndexType > SolverMonitorType;
- ExplicitTimeStepper();
+ static_assert( ProblemType::isTimeDependent(), "The problem is not time dependent." );
- static void configSetup( Config::ConfigDescription& config,
- const String& prefix = "" );
+ ExplicitTimeStepper();
- bool setup( const Config::ParameterContainer& parameters,
- const String& prefix = "" );
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" );
- bool init( const MeshPointer& meshPointer );
+ bool setup( const Config::ParameterContainer& parameters,
+ const String& prefix = "" );
- void setSolver( OdeSolverType& odeSolver );
+ bool init( const MeshPointer& meshPointer );
- void setSolverMonitor( SolverMonitorType& solverMonitor );
+ void setSolver( OdeSolverType& odeSolver );
- void setProblem( ProblemType& problem );
+ void setSolverMonitor( SolverMonitorType& solverMonitor );
- ProblemType* getProblem() const;
+ void setProblem( ProblemType& problem );
- bool setTimeStep( const RealType& tau );
+ ProblemType* getProblem() const;
- const RealType& getTimeStep() const;
+ bool setTimeStep( const RealType& tau );
- bool solve( const RealType& time,
- const RealType& stopTime,
- const MeshPointer& mesh,
- DofVectorPointer& dofVector,
- MeshDependentDataPointer& meshDependentData );
+ const RealType& getTimeStep() const;
- void getExplicitUpdate( const RealType& time,
- const RealType& tau,
- DofVectorPointer& _u,
- DofVectorPointer& _fu );
-
- bool writeEpilog( Logger& logger ) const;
+ bool solve( const RealType& time,
+ const RealType& stopTime,
+ const MeshPointer& mesh,
+ DofVectorPointer& dofVector,
+ MeshDependentDataPointer& meshDependentData );
+
+ void getExplicitUpdate( const RealType& time,
+ const RealType& tau,
+ DofVectorPointer& _u,
+ DofVectorPointer& _fu );
+
+ bool writeEpilog( Logger& logger ) const;
protected:
- OdeSolverType* odeSolver;
+ OdeSolverType* odeSolver;
+
+ SolverMonitorType* solverMonitor;
+
+ Problem* problem;
+
- SolverMonitorType* solverMonitor;
+ RealType timeStep;
- Problem* problem;
+ /****
+ * The pointers on the shared pointer is important here to avoid
+ * memory deallocation in the assignment operator in SharedPointer.
+ */
+ MeshDependentDataPointer* meshDependentData;
+ const MeshPointer* mesh;
- RealType timeStep;
+ Timer preIterateTimer, explicitUpdaterTimer, mainTimer, postIterateTimer;
- /****
- * The pointers on the shared pointer is important here to avoid
- * memory deallocation in the assignment operator in SharedPointer.
- */
- MeshDependentDataPointer* meshDependentData;
-
- const MeshPointer* mesh;
-
- Timer preIterateTimer, explicitUpdaterTimer, mainTimer, postIterateTimer;
-
- long long int allIterations;
+ long long int allIterations;
};
} // namespace PDE
diff --git a/src/TNL/Solvers/SolverStarter_impl.h b/src/TNL/Solvers/SolverStarter_impl.h
index f84699fb9596880a44de560d098025021fb9be25..7137c546e8de44733e24a35bf6abdae95a4b6648 100644
--- a/src/TNL/Solvers/SolverStarter_impl.h
+++ b/src/TNL/Solvers/SolverStarter_impl.h
@@ -135,7 +135,8 @@ class UserDefinedTimeDiscretisationSetter
return false;
}
SolverStarter< ConfigTag > solverStarter;
- return solverStarter.template runPDESolver< Problem, TimeStepper >( problem, parameters, timeStepper );
+ // TODO: Solve the set-up of the DiscreteSOlver type in some better way
+ return solverStarter.template runPDESolver< Problem, TimeStepper, typename Problem::DiscreteSolver >( problem, parameters );
}
};
diff --git a/src/Tools/tnl-view.cpp b/src/Tools/tnl-view.cpp
index 453586e65caf917d8952cda8180f5ac811e8aec7..25e38bc65d4f831379827118af9856a898e737de 100644
--- a/src/Tools/tnl-view.cpp
+++ b/src/Tools/tnl-view.cpp
@@ -26,7 +26,7 @@ namespace BuildConfigTags {
/****
* Turn off support for float and long double.
*/
-template<> struct GridRealTag< TNLViewBuildConfigTag, float > { enum { enabled = false }; };
+//template<> struct GridRealTag< TNLViewBuildConfigTag, float > { enum { enabled = false }; };
template<> struct GridRealTag< TNLViewBuildConfigTag, long double > { enum { enabled = false }; };
/****