From b85dc5c824cc511aab20da08ac5c5cb34e3a27ac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Thu, 14 Jul 2016 21:23:53 +0200
Subject: [PATCH] Refactoring the fast sweeping method.
---
.../tnl-direct-eikonal-solver.h | 1 +
.../tnlDirectEikonalMethodsBase.h | 85 ++++++++
.../tnlDirectEikonalMethodsBase_impl.h | 192 ++++++++++++++++++
.../hamilton-jacobi/tnlDirectEikonalProblem.h | 3 +
.../tnlDirectEikonalProblem_impl.h | 22 +-
.../hamilton-jacobi/tnlFastSweepingMethod.h | 79 ++++++-
.../tnlFastSweepingMethod1D_impl.h | 67 ++++++
.../tnlFastSweepingMethod2D_impl.h | 60 ++++++
.../tnlFastSweepingMethod3D_impl.h | 67 ++++++
src/core/tnlTypeInfo.h | 39 ++++
src/functions/tnlFunctions.h | 1 +
src/functions/tnlMeshFunction_impl.h | 1 +
src/functions/tnlParaboloid.h | 2 +-
src/functions/tnlParaboloidSDF.h | 2 +-
src/functions/tnlParaboloidSDF_impl.h | 14 +-
src/functions/tnlParaboloid_impl.h | 14 +-
src/functions/tnlTestFunction_impl.h | 2 +-
17 files changed, 625 insertions(+), 26 deletions(-)
create mode 100644 examples/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
create mode 100644 examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
create mode 100644 examples/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
create mode 100644 examples/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
create mode 100644 src/core/tnlTypeInfo.h
diff --git a/examples/hamilton-jacobi/tnl-direct-eikonal-solver.h b/examples/hamilton-jacobi/tnl-direct-eikonal-solver.h
index 144e4d540a..0e7a819fc9 100644
--- a/examples/hamilton-jacobi/tnl-direct-eikonal-solver.h
+++ b/examples/hamilton-jacobi/tnl-direct-eikonal-solver.h
@@ -32,6 +32,7 @@ class tnlDirectEikonalSolverConfig
static void configSetup( tnlConfigDescription& config )
{
config.addDelimiter( "Direct eikonal equation solver settings:" );
+ config.addRequiredEntry< tnlString >( "input-file", "Input file." );
};
};
diff --git a/examples/hamilton-jacobi/tnlDirectEikonalMethodsBase.h b/examples/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
new file mode 100644
index 0000000000..049798a911
--- /dev/null
+++ b/examples/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
@@ -0,0 +1,85 @@
+/*
+ * File: tnlDirectEikonalMethodsBase.h
+ * Author: oberhuber
+ *
+ * Created on July 14, 2016, 3:17 PM
+ */
+
+#pragma once
+
+#include <mesh/tnlGrid.h>
+#include <functions/tnlMeshFunction.h>
+
+template< typename Mesh >
+class tnlDirectEikonalMethodsBase
+{
+};
+
+template< typename Real,
+ typename Device,
+ typename Index >
+class tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >
+{
+ public:
+
+ typedef tnlGrid< 1, Real, Device, Index > MeshType;
+ typedef Real RealType;
+ typedef Device DevcieType;
+ typedef Index IndexType;
+
+ template< typename MeshFunction >
+ void initInterface( const MeshFunction& input,
+ MeshFunction& output );
+
+ template< typename MeshFunction, typename MeshEntity >
+ void updateCell( MeshFunction& u,
+ const MeshEntity& cell );
+
+};
+
+
+template< typename Real,
+ typename Device,
+ typename Index >
+class tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >
+{
+ public:
+
+ typedef tnlGrid< 2, Real, Device, Index > MeshType;
+ typedef Real RealType;
+ typedef Device DevcieType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+
+ template< typename MeshFunction >
+ void initInterface( const MeshFunction& input,
+ MeshFunction& output );
+
+ template< typename MeshFunction, typename MeshEntity >
+ void updateCell( MeshFunction& u,
+ const MeshEntity& cell );
+};
+
+template< typename Real,
+ typename Device,
+ typename Index >
+class tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >
+{
+ public:
+
+ typedef tnlGrid< 3, Real, Device, Index > MeshType;
+ typedef Real RealType;
+ typedef Device DevcieType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+
+ template< typename MeshFunction >
+ void initInterface( const MeshFunction& input,
+ MeshFunction& output );
+
+ template< typename MeshFunction, typename MeshEntity >
+ void updateCell( MeshFunction& u,
+ const MeshEntity& cell );
+};
+
+#include "tnlDirectEikonalMethodsBase_impl.h"
diff --git a/examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h b/examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
new file mode 100644
index 0000000000..52c9dc3d60
--- /dev/null
+++ b/examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
@@ -0,0 +1,192 @@
+/*
+ * File: tnlDirectEikonalMethodsBase_impl.h
+ * Author: oberhuber
+ *
+ * Created on July 14, 2016, 3:22 PM
+ */
+
+#pragma once
+
+#include <core/tnlTypeInfo.h>
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< typename MeshFunction >
+void
+tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::
+initInterface( const MeshFunction& input,
+ MeshFunction& output )
+{
+ const MeshType& mesh = input.getMesh();
+ typedef typename MeshType::Cell Cell;
+ Cell cell( mesh );
+ for( cell.getCoordinates().x() = 1;
+ cell.getCoordinates().x() < mesh.getDimensions().x() - 1;
+ cell.getCoordinates().x() ++ )
+ {
+ cell.refresh();
+ const auto& neighbours = cell.getNeighbourEntities();
+ //const IndexType& c = cell.getIndex();
+ const IndexType e = neighbours.template getEntityIndex< 1 >();
+ const IndexType w = neighbours.template getEntityIndex< -1 >();
+ const RealType& c = input( cell );
+ if( c * input[ e ] <= 0 || c * input[ w ] <= 0 )
+ output[ cell.getIndex() ] = c;
+ else output[ cell.getIndex() ] =
+ c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
+ -tnlTypeInfo< RealType >::getMaxValue();
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< typename MeshFunction, typename MeshEntity >
+void
+tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::
+updateCell( MeshFunction& u,
+ const MeshEntity& cell )
+{
+ const auto& neighbourEntities = cell.getNeighbourEntities< 2 >();
+ const MeshType& mesh = cell.getMesh();
+
+ const RealType& value = u( cell );
+ Real a,b, tmp;
+
+ if( cell.getCoordinates().x() == 0 )
+ a = u[ neighbourEntities.template getEntityIndex< 1, 0 >() ];
+ else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
+ a = u[ neighbourEntities.template getEntityIndex< -1, 0 >() ];
+ else
+ {
+ a = fabsMin( u[ neighbourEntities.template getEntityIndex< -1, 0 >() ],
+ u[ neighbourEntities.template getEntityIndex< 1, 0 >() ] );
+ }
+
+ if( cell.getCoordinates().y() == 0 )
+ b = u[ neighbourEntities.template getEntityIndex< 0, 1 >()];
+ else if( cell.getCoordinates().y() == mesh.getDimensions().y() - 1 )
+ b = u[ neighbourEntities.template getEntityIndex< 0, -1 >() ];
+ else
+ {
+ b = fabsMin( u[ neighbourEntities.template getEntityIndex< 0, -1 >() ],
+ u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] );
+ }
+
+ if( fabs( a - b ) >= h )
+ tmp = fabsMin( a, b ) + Sign( value ) * h;
+ else
+ tmp = 0.5 * (a + b + Sign(value)*sqrt(2.0 * h * h - (a - b) * (a - b) ) );
+
+ u[ cell.getIndex() ] = fabsMin( value, tmp );
+}
+
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< typename MeshFunction >
+void
+tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
+initInterface( const MeshFunction& input,
+ MeshFunction& output )
+{
+ const MeshType& mesh = input.getMesh();
+ typedef typename MeshType::Cell Cell;
+ Cell cell( mesh );
+ for( cell.getCoordinates().y() = 0;
+ cell.getCoordinates().y() < mesh.getDimensions().y();
+ cell.getCoordinates().y() ++ )
+ for( cell.getCoordinates().x() = 0;
+ cell.getCoordinates().x() < mesh.getDimensions().x();
+ cell.getCoordinates().x() ++ )
+ {
+ cell.refresh();
+ const RealType& c = input( cell );
+ if( ! cell.isBoundaryEntity() )
+ {
+ auto neighbours = cell.getNeighbourEntities();
+ const IndexType e = neighbours.template getEntityIndex< 1, 0 >();
+ const IndexType w = neighbours.template getEntityIndex< -1, 0 >();
+ const IndexType n = neighbours.template getEntityIndex< 0, 1 >();
+ const IndexType s = neighbours.template getEntityIndex< 0, -1 >();
+ if( c * input[ e ] <= 0 || c * input[ w ] <= 0 ||
+ c * input[ n ] <= 0 || c * input[ s ] <= 0 )
+ {
+ output[ cell.getIndex() ] = c;
+ continue;
+ }
+ }
+ output[ cell.getIndex() ] =
+ c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
+ -tnlTypeInfo< RealType >::getMaxValue();
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< typename MeshFunction, typename MeshEntity >
+void
+tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
+updateCell( MeshFunction& u,
+ const MeshEntity& cell )
+{
+
+}
+
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< typename MeshFunction >
+void
+tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::
+initInterface( const MeshFunction& input,
+ MeshFunction& output )
+{
+ const MeshType& mesh = input.getMesh();
+ typedef typename MeshType::Cell Cell;
+ Cell cell( mesh );
+ for( cell.getCoordinates().z() = 1;
+ cell.getCoordinates().z() < mesh.getDimensions().z() - 1;
+ cell.getCoordinates().z() ++ )
+ for( cell.getCoordinates().y() = 1;
+ cell.getCoordinates().y() < mesh.getDimensions().y() - 1;
+ cell.getCoordinates().y() ++ )
+ for( cell.getCoordinates().x() = 1;
+ cell.getCoordinates().x() < mesh.getDimensions().x() - 1;
+ cell.getCoordinates().x() ++ )
+ {
+ cell.refresh();
+ auto neighbours = cell.getNeighbourEntities();
+ //const IndexType& c = cell.getIndex();
+ const IndexType e = neighbours.template getEntityIndex< 1, 0, 0 >();
+ const IndexType w = neighbours.template getEntityIndex< -1, 0, 0 >();
+ const IndexType n = neighbours.template getEntityIndex< 0, 1, 0 >();
+ const IndexType s = neighbours.template getEntityIndex< 0, -1, 0 >();
+ const IndexType t = neighbours.template getEntityIndex< 0, 0, 1 >();
+ const IndexType b = neighbours.template getEntityIndex< 0, 0, -1 >();
+ const RealType& c = input( cell );
+ if( c * input[ e ] <= 0 || c * input[ w ] <= 0 ||
+ c * input[ n ] <= 0 || c * input[ s ] <= 0 ||
+ c * input[ t ] <= 0 || c * input[ b ] <= 0 )
+ output[ cell.getIndex() ] = c;
+ else output[ cell.getIndex() ] =
+ c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
+ -tnlTypeInfo< RealType >::getMaxValue();
+ }
+}
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< typename MeshFunction, typename MeshEntity >
+void
+tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::
+updateCell( MeshFunction& u,
+ const MeshEntity& cell )
+{
+
+}
diff --git a/examples/hamilton-jacobi/tnlDirectEikonalProblem.h b/examples/hamilton-jacobi/tnlDirectEikonalProblem.h
index 654c775ea9..abe8281426 100644
--- a/examples/hamilton-jacobi/tnlDirectEikonalProblem.h
+++ b/examples/hamilton-jacobi/tnlDirectEikonalProblem.h
@@ -35,6 +35,7 @@ class tnlDirectEikonalProblem
typedef Index IndexType;
typedef tnlMeshFunction< Mesh > MeshFunctionType;
typedef tnlPDEProblem< Mesh, TimeIndependentProblem, RealType, DeviceType, IndexType > BaseType;
+ typedef Anisotropy AnisotropyType;
using typename BaseType::MeshType;
using typename BaseType::DofVectorType;
@@ -71,6 +72,8 @@ class tnlDirectEikonalProblem
MeshFunctionType u;
MeshFunctionType initialData;
+
+ AnisotropyType anisotropy;
};
diff --git a/examples/hamilton-jacobi/tnlDirectEikonalProblem_impl.h b/examples/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
index 1a9f7bcf46..b73f553795 100644
--- a/examples/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
+++ b/examples/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
@@ -32,7 +32,7 @@ tnlString
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
getPrologHeader() const
{
-
+ return tnlString( "Direct eikonal solver" );
}
template< typename Mesh,
@@ -55,7 +55,7 @@ bool
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
writeEpilog( tnlLogger& logger )
{
-
+ return true;
}
template< typename Mesh,
@@ -66,7 +66,7 @@ bool
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
setup( const tnlParameterContainer& parameters )
{
-
+ return true;
}
template< typename Mesh,
@@ -77,7 +77,7 @@ Index
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
getDofs( const MeshType& mesh ) const
{
- return mesh.getEntitiesCount();
+ return mesh.template getEntitiesCount< typename MeshType::Cell >();
}
template< typename Mesh,
@@ -92,7 +92,12 @@ bindDofs( const MeshType& mesh,
this->u.bind( mesh, dofs );
}
+template< typename Mesh,
+ typename Anisotropy,
+ typename Real,
+ typename Index >
bool
+tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
setInitialData( const tnlParameterContainer& parameters,
const MeshType& mesh,
DofVectorType& dofs,
@@ -102,7 +107,7 @@ setInitialData( const tnlParameterContainer& parameters,
this->initialData.setMesh( mesh );
if( !this->initialData.boundLoad( inputFile ) )
return false;
-
+ return true;
}
@@ -113,7 +118,8 @@ template< typename Mesh,
bool
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
solve( const MeshType& mesh,
- DofVectorType& dosf )
+ DofVectorType& dofs )
{
-
-}
\ No newline at end of file
+ tnlFastSweepingMethod< MeshType, AnisotropyType > fsm;
+ fsm.solve( mesh, anisotropy, initialData );
+}
diff --git a/examples/hamilton-jacobi/tnlFastSweepingMethod.h b/examples/hamilton-jacobi/tnlFastSweepingMethod.h
index 0234111ea5..47ce2e32dd 100644
--- a/examples/hamilton-jacobi/tnlFastSweepingMethod.h
+++ b/examples/hamilton-jacobi/tnlFastSweepingMethod.h
@@ -15,18 +15,55 @@
#include <mesh/tnlGrid.h>
#include <functions/tnlConstantFunction.h>
+#include "tnlDirectEikonalMethodsBase.h"
template< typename Mesh,
- typename Anisotropy = tnlConstantFunction< Mesh::getMeshDimensions(), typename Mesh::RealType >
+ typename Anisotropy = tnlConstantFunction< Mesh::getMeshDimensions(), typename Mesh::RealType > >
class tnlFastSweepingMethod
{
};
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+class tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >
+ : public tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >
+{
+ static_assert( std::is_same< Device, tnlHost >::value, "The fast sweeping method works only on CPU." );
+
+ public:
+
+ typedef tnlGrid< 1, Real, Device, Index > MeshType;
+ typedef Real RealType;
+ typedef tnlHost DeviceType;
+ typedef Index IndexType;
+ typedef Anisotropy AnisotropyType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ typedef tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > > BaseType;
+
+ tnlFastSweepingMethod();
+
+ const IndexType& getMaxIterations() const;
+
+ void setMaxIterations( const IndexType& maxIterations );
+
+ void solve( const MeshType& mesh,
+ const AnisotropyType& anisotropy,
+ MeshFunctionType& u );
+
+
+ protected:
+
+ const IndexType maxIterations;
+};
+
template< typename Real,
typename Device,
typename Index,
typename Anisotropy >
class tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >
+ : public tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >
{
static_assert( std::is_same< Device, tnlHost >::value, "The fast sweeping method works only on CPU." );
@@ -38,6 +75,42 @@ class tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >
typedef Index IndexType;
typedef Anisotropy AnisotropyType;
typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ typedef tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > > BaseType;
+
+ tnlFastSweepingMethod();
+
+ const IndexType& getMaxIterations() const;
+
+ void setMaxIterations( const IndexType& maxIterations );
+
+ void solve( const MeshType& mesh,
+ const AnisotropyType& anisotropy,
+ MeshFunctionType& u );
+
+
+ protected:
+
+ const IndexType maxIterations;
+};
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+class tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >
+ : public tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >
+{
+ static_assert( std::is_same< Device, tnlHost >::value, "The fast sweeping method works only on CPU." );
+
+ public:
+
+ typedef tnlGrid< 3, Real, Device, Index > MeshType;
+ typedef Real RealType;
+ typedef tnlHost DeviceType;
+ typedef Index IndexType;
+ typedef Anisotropy AnisotropyType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ typedef tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > > BaseType;
tnlFastSweepingMethod();
@@ -55,4 +128,8 @@ class tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >
const IndexType maxIterations;
};
+
+
+#include "tnlFastSweepingMethod1D_impl.h"
#include "tnlFastSweepingMethod2D_impl.h"
+#include "tnlFastSweepingMethod3D_impl.h"
diff --git a/examples/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h b/examples/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
new file mode 100644
index 0000000000..d6605a954a
--- /dev/null
+++ b/examples/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
@@ -0,0 +1,67 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+/*
+ * File: tnlFastSweepingMethod2D_impl.h
+ * Author: oberhuber
+ *
+ * Created on July 14, 2016, 10:32 AM
+ */
+
+#pragma once
+
+#include "tnlFastSweepingMethod.h"
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
+tnlFastSweepingMethod()
+: maxIterations( 1 )
+{
+
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+const Index&
+tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
+getMaxIterations() const
+{
+
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+void
+tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
+setMaxIterations( const IndexType& maxIterations )
+{
+
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+void
+tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
+solve( const MeshType& mesh,
+ const AnisotropyType& anisotropy,
+ MeshFunctionType& u )
+{
+ MeshFunctionType aux;
+ aux.setMesh( mesh );
+ std::cout << "Initiating the interface cells ..." << std::endl;
+ BaseType::initInterface( u, aux );
+ aux.save( "aux.tnl" );
+}
+
diff --git a/examples/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h b/examples/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
index 1959be39fa..5ce32f5cba 100644
--- a/examples/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
+++ b/examples/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
@@ -21,6 +21,7 @@ template< typename Real,
typename Anisotropy >
tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >::
tnlFastSweepingMethod()
+: maxIterations( 1 )
{
}
@@ -57,6 +58,65 @@ solve( const MeshType& mesh,
const AnisotropyType& anisotropy,
MeshFunctionType& u )
{
+ MeshFunctionType aux;
+ aux.setMesh( mesh );
+ std::cout << "Initiating the interface cells ..." << std::endl;
+ BaseType::initInterface( u, aux );
+ aux.save( "aux.tnl" );
+
+ typename MeshType::Cell cell( mesh );
+
+ for( cell.getCoordinates().y() = 0;
+ cell.getCoordinates().y() < mesh.getDimensions().y();
+ cell.getCoordinates().y()++ )
+ {
+ for( cell.getCoordinates().x() = 0;
+ cell.getCoordinates().x() < mesh.getDimensions().x();
+ cell.getCoordinates().x()++ )
+ {
+ cell.refresh();
+ this->updateValue( aux, cell );
+ }
+ }
+
+ for( cell.getCoordinates().y() = 0;
+ cell.getCoordinates().y() < mesh.getDimensions().y();
+ cell.getCoordinates().y()++ )
+ {
+ for( cell.getCoordinates().x() = mesh.getDimensions().x() - 1;
+ cell.getCoordinates().x() >= 0 ;
+ cell.getCoordinates().x()-- )
+ {
+ cell.refresh();
+ this->updateValue( aux, cell );
+ }
+ }
+
+ for( cell.getCoordinates().y() = mesh.getDimensions().y() - 1;
+ cell.getCoordinates().y() >= 0 ;
+ cell.getCoordinates().y()-- )
+ {
+ for( cell.getCoordinates().x() = 0;
+ cell.getCoordinates().x() < mesh.getDimensions().x();
+ cell.getCoordinates().x()++ )
+ {
+ cell.refresh();
+ this->updateValue( aux, cell );
+ }
+ }
+
+ for( cell.getCoordinates().y() = mesh.getDimensions().y() - 1;
+ cell.getCoordinates().y() >= 0;
+ cell.getCoordinates().y()-- )
+ {
+ for( cell.getCoordinates().x() = mesh.getDimensions().x() - 1;
+ cell.getCoordinates().x() >= 0 ;
+ cell.getCoordinates().x()-- )
+ {
+ cell.refresh();
+ this->updateValue( aux, cell );
+ }
+ }
}
diff --git a/examples/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h b/examples/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
new file mode 100644
index 0000000000..d803fc5b86
--- /dev/null
+++ b/examples/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
@@ -0,0 +1,67 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+/*
+ * File: tnlFastSweepingMethod2D_impl.h
+ * Author: oberhuber
+ *
+ * Created on July 14, 2016, 10:32 AM
+ */
+
+#pragma once
+
+#include "tnlFastSweepingMethod.h"
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
+tnlFastSweepingMethod()
+: maxIterations( 1 )
+{
+
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+const Index&
+tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
+getMaxIterations() const
+{
+
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+void
+tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
+setMaxIterations( const IndexType& maxIterations )
+{
+
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename Anisotropy >
+void
+tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
+solve( const MeshType& mesh,
+ const AnisotropyType& anisotropy,
+ MeshFunctionType& u )
+{
+ MeshFunctionType aux;
+ aux.setMesh( mesh );
+ std::cout << "Initiating the interface cells ..." << std::endl;
+ BaseType::initInterface( u, aux );
+ aux.save( "aux.tnl" );
+}
+
diff --git a/src/core/tnlTypeInfo.h b/src/core/tnlTypeInfo.h
new file mode 100644
index 0000000000..5d5282f3c4
--- /dev/null
+++ b/src/core/tnlTypeInfo.h
@@ -0,0 +1,39 @@
+/*
+ * File: tnlTypeInfo.h
+ * Author: oberhuber
+ *
+ * Created on July 14, 2016, 3:46 PM
+ */
+
+#pragma once
+
+#include <limits>
+
+template< typename Type >
+class tnlTypeInfo
+{
+};
+
+template<>
+class tnlTypeInfo< double >
+{
+ public:
+
+ typedef double Type;
+
+ static __cuda_callable__
+ Type getMaxValue() { return DBL_MAX; };
+};
+
+template<>
+class tnlTypeInfo< float >
+{
+ public:
+
+ typedef float Type;
+
+ static __cuda_callable__
+ Type getMaxValue() { return FLT_MAX; };
+};
+
+
diff --git a/src/functions/tnlFunctions.h b/src/functions/tnlFunctions.h
index 9d6557216f..ea83b9a756 100644
--- a/src/functions/tnlFunctions.h
+++ b/src/functions/tnlFunctions.h
@@ -37,5 +37,6 @@ Real negativePart( const Real& arg)
return arg < 0.0 ? arg : 0.0;
}
+
#endif /* TNLFUNCTIONS_H */
diff --git a/src/functions/tnlMeshFunction_impl.h b/src/functions/tnlMeshFunction_impl.h
index 505a3ce816..fb066f83e2 100644
--- a/src/functions/tnlMeshFunction_impl.h
+++ b/src/functions/tnlMeshFunction_impl.h
@@ -166,6 +166,7 @@ tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
setMesh( const MeshType& mesh )
{
this->mesh = &mesh;
+ this->data.setSize( mesh.template getEntitiesCount< typename Mesh::template MeshEntity< MeshEntityDimensions > >() );
}
template< typename Mesh,
diff --git a/src/functions/tnlParaboloid.h b/src/functions/tnlParaboloid.h
index c29648f65d..45d2eb935c 100644
--- a/src/functions/tnlParaboloid.h
+++ b/src/functions/tnlParaboloid.h
@@ -54,7 +54,7 @@ class tnlParaboloidBase : public tnlDomain< dimensions, SpaceDomain >
protected:
- Real xCentre, yCentre, zCentre, coefficient, offset;
+ Real xCentre, yCentre, zCentre, coefficient, radius;
};
template< int Dimensions, typename Real >
diff --git a/src/functions/tnlParaboloidSDF.h b/src/functions/tnlParaboloidSDF.h
index f998448634..46d03a692d 100644
--- a/src/functions/tnlParaboloidSDF.h
+++ b/src/functions/tnlParaboloidSDF.h
@@ -54,7 +54,7 @@ class tnlParaboloidSDFBase : public tnlDomain< dimensions, SpaceDomain >
protected:
- Real xCentre, yCentre, zCentre, coefficient, offset;
+ Real xCentre, yCentre, zCentre, coefficient, radius;
};
template< int Dimensions, typename Real >
diff --git a/src/functions/tnlParaboloidSDF_impl.h b/src/functions/tnlParaboloidSDF_impl.h
index 63ebc4386b..12738687e5 100644
--- a/src/functions/tnlParaboloidSDF_impl.h
+++ b/src/functions/tnlParaboloidSDF_impl.h
@@ -22,7 +22,7 @@
template< int dimensions, typename Real >
tnlParaboloidSDFBase< dimensions, Real >::tnlParaboloidSDFBase()
: xCentre( 0 ), yCentre( 0 ), zCentre( 0 ),
- coefficient( 1 ), offset ( 0 )
+ coefficient( 1 ), radius ( 0 )
{
}
@@ -34,7 +34,7 @@ bool tnlParaboloidSDFBase< dimensions, Real >::setup( const tnlParameterContaine
this->yCentre = parameters.getParameter< double >( "y-centre" );
this->zCentre = parameters.getParameter< double >( "z-centre" );
this->coefficient = parameters.getParameter< double >( "coefficient" );
- this->offset = parameters.getParameter< double >( "offset" );
+ this->radius = parameters.getParameter< double >( "radius" );
return true;
}
@@ -89,13 +89,13 @@ Real tnlParaboloidSDFBase< dimensions, Real >::getCoefficient() const
template< int dimensions, typename Real >
void tnlParaboloidSDFBase< dimensions, Real >::setOffset( const Real& offset )
{
- this->offset = offset;
+ this->radius = offset;
}
template< int dimensions, typename Real >
Real tnlParaboloidSDFBase< dimensions, Real >::getOffset() const
{
- return this->offset;
+ return this->radius;
}
template< typename Real >
@@ -112,7 +112,7 @@ getPartialDerivative( const VertexType& v,
if( YDiffOrder != 0 || ZDiffOrder != 0 )
return 0.0;
if( XDiffOrder == 0 )
- return sqrt( ( x - this -> xCentre ) * ( x - this -> xCentre ) ) - this->offset;
+ return sqrt( ( x - this -> xCentre ) * ( x - this -> xCentre ) ) - this->radius;
if( XDiffOrder == 1 )
return 1.0;
return 0.0;
@@ -134,7 +134,7 @@ getPartialDerivative( const VertexType& v,
if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
{
return sqrt ( ( x - this -> xCentre ) * ( x - this -> xCentre )
- + ( y - this -> yCentre ) * ( y - this -> yCentre ) ) - this->offset;
+ + ( y - this -> yCentre ) * ( y - this -> yCentre ) ) - this->radius;
}
return 0.0;
}
@@ -156,7 +156,7 @@ getPartialDerivative( const VertexType& v,
{
return sqrt( ( x - this -> xCentre ) * ( x - this -> xCentre )
+ ( y - this -> yCentre ) * ( y - this -> yCentre )
- + ( z - this -> zCentre ) * ( z - this -> zCentre ) ) - this->offset;
+ + ( z - this -> zCentre ) * ( z - this -> zCentre ) ) - this->radius;
}
return 0.0;
}
diff --git a/src/functions/tnlParaboloid_impl.h b/src/functions/tnlParaboloid_impl.h
index ae4debbaff..f3b6572009 100644
--- a/src/functions/tnlParaboloid_impl.h
+++ b/src/functions/tnlParaboloid_impl.h
@@ -22,7 +22,7 @@
template< int dimensions, typename Real >
tnlParaboloidBase< dimensions, Real >::tnlParaboloidBase()
: xCentre( 0 ), yCentre( 0 ), zCentre( 0 ),
- coefficient( 1 ), offset ( 0 )
+ coefficient( 1 ), radius ( 0 )
{
}
@@ -34,7 +34,7 @@ bool tnlParaboloidBase< dimensions, Real >::setup( const tnlParameterContainer&
this->yCentre = parameters.getParameter< double >( "y-centre" );
this->zCentre = parameters.getParameter< double >( "z-centre" );
this->coefficient = parameters.getParameter< double >( "coefficient" );
- this->offset = parameters.getParameter< double >( "offset" );
+ this->radius = parameters.getParameter< double >( "radius" );
return true;
}
@@ -89,13 +89,13 @@ Real tnlParaboloidBase< dimensions, Real >::getCoefficient() const
template< int dimensions, typename Real >
void tnlParaboloidBase< dimensions, Real >::setOffset( const Real& offset )
{
- this->offset = offset;
+ this->radius = offset;
}
template< int dimensions, typename Real >
Real tnlParaboloidBase< dimensions, Real >::getOffset() const
{
- return this->offset;
+ return this->radius;
}
template< typename Real >
@@ -112,7 +112,7 @@ getPartialDerivative( const VertexType& v,
if( YDiffOrder != 0 || ZDiffOrder != 0 )
return 0.0;
if( XDiffOrder == 0 )
- return this->coefficient * ( ( x - this -> xCentre ) * ( x - this -> xCentre ) - this->offset*this->offset );
+ return this->coefficient * ( ( x - this -> xCentre ) * ( x - this -> xCentre ) - this->radius*this->radius );
if( XDiffOrder == 1 )
return 2.0 * this->coefficient * ( x - this -> xCentre );
return 0.0;
@@ -136,7 +136,7 @@ getPartialDerivative( const VertexType& v,
if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
{
return this->coefficient * ( ( x - this -> xCentre ) * ( x - this -> xCentre )
- + ( y - this -> yCentre ) * ( y - this -> yCentre ) - this->offset*this->offset );
+ + ( y - this -> yCentre ) * ( y - this -> yCentre ) - this->radius*this->radius );
}
if( XDiffOrder == 1 && YDiffOrder == 0)
return 2.0 * this->coefficient * ( x - this -> xCentre );
@@ -166,7 +166,7 @@ getPartialDerivative( const VertexType& v,
{
return this->coefficient * ( ( x - this -> xCentre ) * ( x - this -> xCentre )
+ ( y - this -> yCentre ) * ( y - this -> yCentre )
- + ( z - this -> zCentre ) * ( z - this -> zCentre ) - this->offset*this->offset );
+ + ( z - this -> zCentre ) * ( z - this -> zCentre ) - this->radius*this->radius );
}
if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 0)
return 2.0 * this->coefficient * ( x - this -> xCentre );
diff --git a/src/functions/tnlTestFunction_impl.h b/src/functions/tnlTestFunction_impl.h
index 4eaf0dfcf2..86483459b3 100644
--- a/src/functions/tnlTestFunction_impl.h
+++ b/src/functions/tnlTestFunction_impl.h
@@ -86,7 +86,7 @@ configSetup( tnlConfigDescription& config,
config.addEntry < double >( prefix + "waves-number-y", "Cut-off for the sine based test functions.", 0.0 );
config.addEntry < double >( prefix + "waves-number-z", "Cut-off for the sine based test functions.", 0.0 );
config.addEntry < double >( prefix + "sigma", "Sigma for the exp based test functions.", 1.0 );
- config.addEntry < double >( prefix + "offset", "Offset for paraboloids.", 1.0 );
+ config.addEntry < double >( prefix + "radius", "Radius for paraboloids.", 1.0 );
config.addEntry < double >( prefix + "coefficient", "Coefficient for paraboloids.", 1.0 );
config.addEntry < double >( prefix + "x-centre", "x-centre for paraboloids.", 0.0 );
config.addEntry < double >( prefix + "y-centre", "y-centre for paraboloids.", 0.0 );
--
GitLab