Implementing vector field evaluator. (59f59bcd) · Commits · TNL / tnl-dev

src/TNL/Containers/Algorithms/ArrayOperations.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -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?

src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h

+10 −2

Original line number	Diff line number	Diff line
		@@ -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 >

src/TNL/Functions/VectorField.h

+61 −4

Original line number	Diff line number	Diff line
		@@ -91,6 +91,11 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimension, Real > >
		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 +151,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,
		@@ -165,6 +174,38 @@ 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 )
		{
		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 )
		{
		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 +218,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 +250,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;
		}

src/TNL/Functions/VectorFieldEvaluator.h

0 → 100644

+173 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		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>

src/TNL/Functions/VectorFieldEvaluator_impl.h

0 → 100644

+173 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		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