From 969e79a4c2ca23e96c163a8d3fbaf1b507f5d023 Mon Sep 17 00:00:00 2001
From: Ondrej Szekely <ondra.szekely@gmail.com>
Date: Mon, 4 May 2015 01:09:11 +0200
Subject: [PATCH] MEAN CURVATURE FLOW - EXACT ALL DIMENSIONS, SEMI-IMPLICIT ALL
DIMENSIONS
---
examples/mean-curvature-flow/CMakeLists.txt | 19 +-
examples/mean-curvature-flow/Makefile | 6 +-
.../tnl-mean-curvature-flow-eoc.h | 20 +-
.../tnl-mean-curvature-flow.h | 12 +-
.../tnl-run-mean-curvature-flow-eoc-test | 12 +-
.../tnlBackwardFiniteDifference.h | 128 -----
.../tnlBackwardFiniteDifference_impl.h | 91 ----
.../tnlExactFlowDiffusion_impl.h | 59 ---
.../tnlForwardFiniteDifference.h | 128 -----
.../tnlForwardFiniteDifference_impl.h | 91 ----
.../tnlMeanCurvatureFlowDiffusion_impl.h | 104 ----
examples/mean-curvature-flow/tnlOperatorQ.h | 135 ------
.../mean-curvature-flow/tnlOperatorQ_impl.h | 93 ----
src/functions/tnlExpBumpFunction_impl.h | 8 +
src/functions/tnlSinBumpsFunction_impl.h | 8 +
src/functions/tnlSinWaveFunction_impl.h | 8 +
src/operators/CMakeLists.txt | 4 +-
src/operators/diffusion/CMakeLists.txt | 8 +-
.../CMakeLists.txt | 4 +
.../tnlOneSideDiffNonlinearOperator.h | 187 ++++++++
.../tnlOneSideDiffNonlinearOperator.h~ | 75 ++-
.../tnlOneSideDiffNonlinearOperator_impl.h | 320 +++++++++++++
.../tnlOneSideDiffNonlinearOperator_impl.h~ | 84 ++++
.../diffusion/tnlExactNonlinearDiffusion.h | 42 +-
.../tnlExactNonlinearDiffusion_impl.h | 96 ++++
.../diffusion/tnlNonlinearDiffusion.h | 187 ++++++++
.../diffusion/tnlNonlinearDiffusion_impl.h | 249 ++++++++++
src/operators/operator-Q/CMakeLists.txt | 6 +
src/operators/operator-Q/tnlExactOperatorQ.h | 90 ++++
.../operator-Q/tnlExactOperatorQ_impl.h | 116 +++++
.../tnlOneSideDiffOperatorQForGraph.h | 327 +++++++++++++
.../tnlOneSideDiffOperatorQForGraph_impl.h | 450 ++++++++++++++++++
.../operator-curvature/CMakeLists.txt | 4 +
.../tnlExactOperatorCurvature.h | 93 ++++
.../tnlExactOperatorCurvature_impl.h | 88 ++++
.../tnlMeanCurvatureFlowEocProblem.h | 14 +-
.../tnlMeanCurvatureFlowEocProblem_impl.h | 0
.../problems}/tnlMeanCurvatureFlowEocRhs.h | 0
.../problems}/tnlMeanCurvatureFlowProblem.h | 24 +-
.../tnlMeanCurvatureFlowProblem_impl.h | 67 +--
40 files changed, 2470 insertions(+), 987 deletions(-)
delete mode 100644 examples/mean-curvature-flow/tnlBackwardFiniteDifference.h
delete mode 100644 examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h
delete mode 100644 examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h
delete mode 100644 examples/mean-curvature-flow/tnlForwardFiniteDifference.h
delete mode 100644 examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h
delete mode 100644 examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h
delete mode 100644 examples/mean-curvature-flow/tnlOperatorQ.h
delete mode 100644 examples/mean-curvature-flow/tnlOperatorQ_impl.h
create mode 100755 src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt
create mode 100644 src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h
rename examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h => src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h~ (73%)
create mode 100644 src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h
create mode 100644 src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h~
rename examples/mean-curvature-flow/tnlExactFlowDiffusion.h => src/operators/diffusion/tnlExactNonlinearDiffusion.h (77%)
create mode 100644 src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
create mode 100644 src/operators/diffusion/tnlNonlinearDiffusion.h
create mode 100644 src/operators/diffusion/tnlNonlinearDiffusion_impl.h
create mode 100755 src/operators/operator-Q/CMakeLists.txt
create mode 100644 src/operators/operator-Q/tnlExactOperatorQ.h
create mode 100644 src/operators/operator-Q/tnlExactOperatorQ_impl.h
create mode 100644 src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h
create mode 100644 src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph_impl.h
create mode 100755 src/operators/operator-curvature/CMakeLists.txt
create mode 100644 src/operators/operator-curvature/tnlExactOperatorCurvature.h
create mode 100644 src/operators/operator-curvature/tnlExactOperatorCurvature_impl.h
rename {examples/mean-curvature-flow => src/problems}/tnlMeanCurvatureFlowEocProblem.h (66%)
rename {examples/mean-curvature-flow => src/problems}/tnlMeanCurvatureFlowEocProblem_impl.h (100%)
rename {examples/mean-curvature-flow => src/problems}/tnlMeanCurvatureFlowEocRhs.h (100%)
rename {examples/mean-curvature-flow => src/problems}/tnlMeanCurvatureFlowProblem.h (78%)
rename {examples/mean-curvature-flow => src/problems}/tnlMeanCurvatureFlowProblem_impl.h (80%)
diff --git a/examples/mean-curvature-flow/CMakeLists.txt b/examples/mean-curvature-flow/CMakeLists.txt
index 9497b2f7a0..88bfec479c 100755
--- a/examples/mean-curvature-flow/CMakeLists.txt
+++ b/examples/mean-curvature-flow/CMakeLists.txt
@@ -1,23 +1,8 @@
-INSTALL( FILES tnlBackwardFiniteDifference.h
- tnlBackwardFiniteDifference_impl.h
- tnlExactFlowDiffusion.h
- tnlExactFlowDiffusion_impl.h
- tnlForwardFiniteDifference.h
- tnlForwardFiniteDifference_impl.h
- tnl-mean-curvature-flow.cpp
+INSTALL( FILES tnl-mean-curvature-flow.cpp
tnl-mean-curvature-flow.h
- tnlMeanCurvatureFlowDiffusion.h
- tnlMeanCurvatureFlowDiffusion_impl.h
tnl-mean-curvature-flow-eoc.cpp
tnl-mean-curvature-flow-eoc.h
- tnlMeanCurvatureFlowEocProblem.h
- tnlMeanCurvatureFlowEocProblem_impl.h
- tnlMeanCurvatureFlowEocRhs.h
- tnlMeanCurvatureFlowEocRhs_impl.h
- tnlMeanCurvatureFlowProblem.h
- tnlMeanCurvatureFlowProblem_impl.h
- tnlOperatorQ.h
- tnlOperatorQ_impl.h
tnl-run-mean-curvature-flow
tnl-run-mean-curvature-flow-eoc-test
+ Makefile
DESTINATION share/tnl-${tnlVersion}/examples/mean-curvature-flow )
diff --git a/examples/mean-curvature-flow/Makefile b/examples/mean-curvature-flow/Makefile
index 9ca9bc675c..71873be3fb 100644
--- a/examples/mean-curvature-flow/Makefile
+++ b/examples/mean-curvature-flow/Makefile
@@ -6,7 +6,7 @@ TARGET = mean-curvature-flow-eoc
INSTALL_DIR = ${HOME}/local
CXX = g++
CUDA_CXX = nvcc
-CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -g -O0
+CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O0 -g
LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-dbg-0.1
SOURCES = tnl-mean-curvature-flow-eoc.cpp
@@ -23,11 +23,11 @@ dist: $(DIST)
install: $(TARGET)
cp $(TARGET) $(INSTALL_DIR)/bin
- cp $(INSTALL_DIR)/share/tnl-0.1/examples/mean-curvature-flow-eoc
+ cp $(INSTALL_DIR)/share/tnl-0.1/examples/mean-curvature-flow
uninstall: $(TARGET)
rm -f $(INSTALL_DIR)/bin/$(TARGET)
- rm -f $(INSTALL_DIR)/share/tnl-0.1/examples/mean-curvature-flow-eoc
+ rm -f $(INSTALL_DIR)/share/tnl-0.1/examples/mean-curvature-flow
$(TARGET): $(OBJECTS)
$(CXX) -o $(TARGET) $(OBJECTS) $(LD_FLAGS)
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
index 0b09c6059b..90ff6ca323 100644
--- a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
@@ -24,10 +24,14 @@
#include <functions/tnlTestFunction.h>
#include <operators/tnlAnalyticDirichletBoundaryConditions.h>
#include <operators/tnlAnalyticNeumannBoundaryConditions.h>
-#include "tnlMeanCurvatureFlowEocRhs.h"
-#include "tnlMeanCurvatureFlowEocProblem.h"
-#include "tnlExactFlowDiffusion.h"
-#include "tnlMeanCurvatureFlowDiffusion.h"
+#include <problems/tnlMeanCurvatureFlowEocRhs.h>
+#include <problems/tnlMeanCurvatureFlowEocProblem.h>
+#include <operators/diffusion/tnlExactNonlinearDiffusion.h>
+#include <operators/diffusion/tnlNonlinearDiffusion.h>
+#include <operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h>
+#include <operators/diffusion/tnlExactNonlinearDiffusion.h>
+#include <operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h>
+#include <operators/operator-Q/tnlExactOperatorQ.h>
//typedef tnlDefaultConfigTag BuildConfig;
typedef tnlFastBuildConfig BuildConfig;
@@ -63,12 +67,14 @@ class meanCurvatureFlowEocSetter
static bool run( const tnlParameterContainer& parameters )
{
enum { Dimensions = MeshType::Dimensions };
- typedef tnlMeanCurvatureFlowDiffusion< MeshType, Real, Index > ApproximateOperator;
- typedef tnlExactFlowDiffusion< Dimensions > ExactOperator;
+ typedef tnlOneSideDiffOperatorQForGraph<MeshType, Real, Index, 0> OperatorQ;
+ typedef tnlOneSideDiffNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;
+ typedef tnlNonlinearDiffusion< MeshType, NonlinearOperator, Real, Index > ApproximateOperator;
+ typedef tnlExactNonlinearDiffusion< tnlExactOperatorQ<Dimensions>, Dimensions > ExactOperator;
typedef tnlTestFunction< MeshType::Dimensions, Real, Device > TestFunction;
typedef tnlMeanCurvatureFlowEocRhs< ExactOperator, TestFunction > RightHandSide;
typedef tnlStaticVector < MeshType::Dimensions, Real > Vertex;
- typedef tnlAnalyticNeumannBoundaryConditions< MeshType, TestFunction, Real, Index > BoundaryConditions;
+ typedef tnlAnalyticDirichletBoundaryConditions< MeshType, TestFunction, Real, Index > BoundaryConditions;
typedef tnlMeanCurvatureFlowEocProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
SolverStarter solverStarter;
return solverStarter.template run< Solver >( parameters );
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow.h b/examples/mean-curvature-flow/tnl-mean-curvature-flow.h
index a363b26c46..00a6d6dc7c 100644
--- a/examples/mean-curvature-flow/tnl-mean-curvature-flow.h
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow.h
@@ -27,10 +27,10 @@
#include <operators/tnlAnalyticNeumannBoundaryConditions.h>
#include <operators/tnlNeumannBoundaryConditions.h>
#include <functions/tnlConstantFunction.h>
-#include <problems/tnlHeatEquationProblem.h>
-
-#include "tnlMeanCurvatureFlowProblem.h"
-#include "tnlMeanCurvatureFlowDiffusion.h"
+#include <problems/tnlMeanCurvatureFlowProblem.h>
+#include <operators/diffusion/tnlNonlinearDiffusion.h>
+#include <operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h>
+#include <operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h>
//typedef tnlDefaultConfigTag BuildConfig;
typedef tnlFastBuildConfig BuildConfig;
@@ -71,7 +71,9 @@ class meanCurvatureFlowSetter
static bool run( const tnlParameterContainer& parameters )
{
enum { Dimensions = MeshType::Dimensions };
- typedef tnlMeanCurvatureFlowDiffusion< MeshType, Real, Index > ApproximateOperator;
+ typedef tnlOneSideDiffOperatorQForGraph<MeshType, Real, Index, 0> OperatorQ;
+ typedef tnlOneSideDiffNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;
+ typedef tnlNonlinearDiffusion< MeshType, NonlinearOperator, Real, Index > ApproximateOperator;
typedef tnlConstantFunction< Dimensions, Real > RightHandSide;
typedef tnlStaticVector < MeshType::Dimensions, Real > Vertex;
diff --git a/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
index abf3827cd5..aa151c5a8a 100755
--- a/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
+++ b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
@@ -2,10 +2,10 @@
device="host"
dimensions="2D"
-sizes2D="16 32 64"
+sizes2D="16 32"
testFunctions="sin-bumps"
snapshotPeriod=0.1
-finalTime=1.0
+finalTime=0.3
timeDependence="cosine"
solverName="mean-curvature-flow-eoc"
#solverName="gdb --args tnl-heat-equation-eoc-test-dbg"
@@ -83,12 +83,12 @@ solve()
--mesh mesh.tnl \
--initial-condition exact-u-00000.tnl \
--time-discretisation ${timeDiscretisation} \
- --time-step 10 \
+ --time-step 1 \
--time-step-order 2 \
--discrete-solver ${discreteSolver} \
--merson-adaptivity 1.0e-7 \
--sor-omega 1.95 \
- --gmres-restarting 10 \
+ --gmres-restarting 20 \
--min-iterations 20 \
--convergence-residue 1.0e-12 \
--test-function ${testFunction}\
@@ -168,8 +168,8 @@ runTest()
echo "========================================================================="
echo "=== STARTING THE SOLVER ==="
echo "========================================================================="
- solve explicit merson
- #solve semi-implicit gmres
+ #solve explicit merson
+ solve semi-implicit gmres
mv computation-in-progress computation-done
fi
echo "========================================================================="
diff --git a/examples/mean-curvature-flow/tnlBackwardFiniteDifference.h b/examples/mean-curvature-flow/tnlBackwardFiniteDifference.h
deleted file mode 100644
index b685d9df8f..0000000000
--- a/examples/mean-curvature-flow/tnlBackwardFiniteDifference.h
+++ /dev/null
@@ -1,128 +0,0 @@
-#ifndef TNLBACKWARDFINITEDIFFERENCE_H
-#define TNLBACKWARDFINITEDIFFERENCE_H
-
-#include <core/vectors/tnlVector.h>
-#include <mesh/tnlGrid.h>
-
-template< typename Mesh,
- typename Real = typename Mesh::RealType,
- typename Index = typename Mesh::IndexType >
-class tnlBackwardFiniteDifference
-{
-
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlBackwardFiniteDifference< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const;
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const;
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlBackwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-};
-
-
-#include "tnlBackwardFiniteDifference_impl.h"
-
-
-#endif /* TNLBACKWARDFINITEDIFFERENCE_H */
diff --git a/examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h b/examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h
deleted file mode 100644
index fe3c3e0774..0000000000
--- a/examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h
+++ /dev/null
@@ -1,91 +0,0 @@
-
-#ifndef TNLBACKWARDFINITEDIFFERENCE_IMPL_H
-#define TNLBACKWARDFINITEDIFFERENCE_IMPL_H
-
-#include "tnlBackwardFiniteDifference.h"
-#include <mesh/tnlGrid.h>
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlBackwardFiniteDifference< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlBackwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlBackwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
-{
- return ( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ]
- + u[ cellIndex ]) * mesh.getHxInverse();
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
-{
- return ( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ]
- + u[ cellIndex ]) * mesh.getHyInverse();
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlBackwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlBackwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-#endif /* TNLBACKWARDFINITEDIFFERENCE_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h b/examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h
deleted file mode 100644
index 05d623831e..0000000000
--- a/examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h
+++ /dev/null
@@ -1,59 +0,0 @@
-/***************************************************************************
- tnlExactLinearDiffusion_impl.h - description
- -------------------
- begin : Aug 8, 2014
- copyright : (C) 2014 by Tomas Oberhuber
- email : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- * *
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation; either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- ***************************************************************************/
-
-#ifndef TNLEXACTFLOWDIFFUSION_IMPL_H_
-#define TNLEXACTFLOWDIFFUSION_IMPL_H_
-
-tnlString
-tnlExactFlowDiffusion< 1 >::
-getType()
-{
- return "tnlExactFlowDiffusion< 1 >";
-}
-
-tnlString
-tnlExactFlowDiffusion< 2 >::
-getType()
-{
- return "tnlExactFlowDiffusion< 2 >";
-}
-
-template< typename Function, typename Vertex, typename Real >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
-Real
-tnlExactFlowDiffusion< 2 >::
-getValue( const Function& function,
- const Vertex& v,
- const Real& time )
-{
- return (function.template getValue< 2, 0, 0, Vertex >( v, time )*(1.0+ function.template getValue< 0, 1, 0, Vertex >( v, time )*function.template getValue< 0, 1, 0, Vertex >( v, time ))
- - 2.0*function.template getValue< 0, 1, 0, Vertex >( v, time )*function.template getValue< 1, 0, 0, Vertex >( v, time )*function.template getValue< 1, 1, 0, Vertex >( v, time )
- + function.template getValue< 0, 2, 0, Vertex >( v, time )*(1.0+ function.template getValue< 1, 0, 0, Vertex >( v, time )*function.template getValue< 1, 0, 0, Vertex >( v, time )))
- / (1.0+function.template getValue< 1, 0, 0, Vertex >( v, time )*function.template getValue< 1, 0, 0, Vertex >( v, time )
- +function.template getValue< 0, 1, 0, Vertex >( v, time )*function.template getValue< 0, 1, 0, Vertex >( v, time ));
-}
-
-tnlString
-tnlExactFlowDiffusion< 3 >::
-getType()
-{
- return "tnlExactFlowDiffusion< 3 >";
-}
-
-#endif /* TNLEXACTFLOWDIFFUSION_IMPL_H_ */
diff --git a/examples/mean-curvature-flow/tnlForwardFiniteDifference.h b/examples/mean-curvature-flow/tnlForwardFiniteDifference.h
deleted file mode 100644
index f207879fff..0000000000
--- a/examples/mean-curvature-flow/tnlForwardFiniteDifference.h
+++ /dev/null
@@ -1,128 +0,0 @@
-#ifndef TNLFORWARDFINITEDIFFERENCE_H
-#define TNLFORWARDFINITEDIFFERENCE_H
-
-#include <core/vectors/tnlVector.h>
-#include <mesh/tnlGrid.h>
-
-template< typename Mesh,
- typename Real = typename Mesh::RealType,
- typename Index = typename Mesh::IndexType >
-class tnlForwardFiniteDifference
-{
-
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlForwardFiniteDifference< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const;
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const;
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlForwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-};
-
-
-#include "tnlForwardFiniteDifference_impl.h"
-
-
-#endif /* TNLFORWARDFINITEDIFFERENCE_H */
diff --git a/examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h b/examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h
deleted file mode 100644
index 419adb4786..0000000000
--- a/examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h
+++ /dev/null
@@ -1,91 +0,0 @@
-
-#ifndef TNLFORWARDFINITEDIFFERENCE_IMPL_H
-#define TNLFORWARDFINITEDIFFERENCE_IMPL_H
-
-#include "tnlForwardFiniteDifference.h"
-#include <mesh/tnlGrid.h>
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlForwardFiniteDifference< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlForwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlForwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValueX( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
-{
- return ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ]
- - u[ cellIndex ]) * mesh.getHxInverse();
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValueY( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
-{
- return ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ]
- - u[ cellIndex ]) * mesh.getHyInverse();
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlForwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlForwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-#endif /* TNLFORWARDFINITEDIFFERENCE_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h
deleted file mode 100644
index 061af10f40..0000000000
--- a/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h
+++ /dev/null
@@ -1,104 +0,0 @@
-
-#ifndef TNLMEANCURVATIVEFLOWDIFFUSION_IMPL_H
-#define TNLMEANCURVATIVEFLOWDIFFUSION_IMPL_H
-
-#include "tnlMeanCurvatureFlowDiffusion.h"
-#include "tnlForwardFiniteDifference.h"
-#include "tnlBackwardFiniteDifference.h"
-
-#include <mesh/tnlGrid.h>
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlMeanCurvatureFlowDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlForwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlForwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValue( const MeshType& mesh,
- const IndexType cellIndex,
- const CoordinatesType& coordinates,
- const Vector& u,
- const Real& time ) const
-{
- return operatorQ.getValueStriped(mesh,cellIndex,u)*((fDifference.getValueX(mesh,cellIndex,u)/operatorQ.getValue(mesh,cellIndex,u)
- -bDifference.getValueX(mesh,cellIndex,u)/operatorQ.getValue(mesh,mesh.template getCellNextToCell<-1,0>(cellIndex),u))
- *mesh.getHxInverse()+(fDifference.getValueY(mesh,cellIndex,u)/operatorQ.getValue(mesh,cellIndex,u)
- -bDifference.getValueY(mesh,cellIndex,u)/operatorQ.getValue(mesh,mesh.template getCellNextToCell<0,-1>(cellIndex),u))*mesh.getHyInverse());
-}
-//template< typename MeshReal,
-// typename Device,
-// typename MeshIndex,
-// typename Real,
-// typename Index >
-//void tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::setupDofs ( const MeshType& mesh )
-//{
-// dofs.setSize(mesh.getNumberOfCells());
-//}
-//
-//template< typename MeshReal,
-// typename Device,
-// typename MeshIndex,
-// typename Real,
-// typename Index >
-//template< typename Vector>
-//void
-//tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-//computeFirstGradient( const MeshType& mesh,
-// const RealType& time,
-// const Vector u)
-//{
-//
-//}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlForwardFiniteDifference< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-#endif /* TNLMEANCURVATIVEFLOWDIFFUSION_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlOperatorQ.h b/examples/mean-curvature-flow/tnlOperatorQ.h
deleted file mode 100644
index 5887515c60..0000000000
--- a/examples/mean-curvature-flow/tnlOperatorQ.h
+++ /dev/null
@@ -1,135 +0,0 @@
-#ifndef TNLOPERATORQ_H
-#define TNLOPERATORQ_H
-
-#include <core/vectors/tnlVector.h>
-#include <mesh/tnlGrid.h>
-#include "tnlBackwardFiniteDifference.h"
-#include "tnlForwardFiniteDifference.h"
-
-template< typename Mesh,
- typename Real = typename Mesh::RealType,
- typename Index = typename Mesh::IndexType >
-class tnlOperatorQ
-{
-
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlOperatorQ< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValue( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueStriped( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValue( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const;
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueStriped( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const;
-
- private:
-
- tnlForwardFiniteDifference<MeshType> fDifference;
- tnlBackwardFiniteDifference<MeshType> bDifference;
-};
-
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-class tnlOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
-{
- public:
-
- typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
- typedef typename MeshType::CoordinatesType CoordinatesType;
- typedef Real RealType;
- typedef Device DeviceType;
- typedef Index IndexType;
-
- static tnlString getType();
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValue( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-
- template< typename Vector >
-#ifdef HAVE_CUDA
- __device__ __host__
-#endif
- Real getValueStriped( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
- {}
-};
-
-
-#include "tnlOperatorQ_impl.h"
-
-
-#endif /* TNLOPERATORQ_H */
diff --git a/examples/mean-curvature-flow/tnlOperatorQ_impl.h b/examples/mean-curvature-flow/tnlOperatorQ_impl.h
deleted file mode 100644
index dd60e10a32..0000000000
--- a/examples/mean-curvature-flow/tnlOperatorQ_impl.h
+++ /dev/null
@@ -1,93 +0,0 @@
-
-#ifndef TTNLOPERATORQ_IMPL_H
-#define TTNLOPERATORQ_IMPL_H
-
-#include "tnlOperatorQ.h"
-#include <mesh/tnlGrid.h>
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlOperatorQ< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlOperatorQ< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValue( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
-{
- return sqrt( 1.0 + (fDifference.getValueX(mesh,cellIndex,u))*(fDifference.getValueX(mesh,cellIndex,u))
- + (fDifference.getValueY(mesh,cellIndex,u))*(fDifference.getValueY(mesh,cellIndex,u)) );
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-template< typename Vector >
-#ifdef HAVE_CUDA
-__device__ __host__
-#endif
-Real
-tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getValueStriped( const MeshType& mesh,
- const IndexType cellIndex,
- const Vector& u) const
-{
- return sqrt( 1.0 + 0.5*( (fDifference.getValueX(mesh,cellIndex,u))*(fDifference.getValueX(mesh,cellIndex,u))
- + (fDifference.getValueY(mesh,cellIndex,u))*(fDifference.getValueY(mesh,cellIndex,u))
- + (bDifference.getValueX(mesh,cellIndex,u))*(bDifference.getValueX(mesh,cellIndex,u))
- + (bDifference.getValueY(mesh,cellIndex,u))*(bDifference.getValueY(mesh,cellIndex,u)) ) );
-}
-
-template< typename MeshReal,
- typename Device,
- typename MeshIndex,
- typename Real,
- typename Index >
-tnlString
-tnlOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-getType()
-{
- return tnlString( "tnlOperatorQ< " ) +
- MeshType::getType() + ", " +
- ::getType< Real >() + ", " +
- ::getType< Index >() + " >";
-}
-
-#endif /* TTNLOPERATORQ_IMPL_H */
diff --git a/src/functions/tnlExpBumpFunction_impl.h b/src/functions/tnlExpBumpFunction_impl.h
index e701d9ee06..d00ebae93f 100644
--- a/src/functions/tnlExpBumpFunction_impl.h
+++ b/src/functions/tnlExpBumpFunction_impl.h
@@ -138,6 +138,8 @@ getValue( const Vertex& v,
return -2.0 * y / ( this->sigma * this->sigma ) * this->amplitude * exp( (-x * x - y * y)/ ( this->sigma * this->sigma ) );
if( XDiffOrder == 0 && YDiffOrder == 2 )
return -2.0 / ( this->sigma * this->sigma ) * this->amplitude * exp( (-x*x - y*y) / ( this->sigma * this->sigma ) ) + 4.0 * y * y / ( this->sigma * this->sigma * this->sigma * this->sigma ) * this->amplitude * exp( (-x*x - y*y) / ( this->sigma * this->sigma ) );
+ if( XDiffOrder == 1 && YDiffOrder == 1 )
+ return 4.0 * x * y / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( (-x * x - y * y)/ ( this->sigma * this->sigma ) );
return 0.0;
}
@@ -187,6 +189,12 @@ getValue( const Vertex& v,
return -2.0 * z / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 2 )
return -2.0 / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) ) + 4.0 * z * z / ( this->sigma * this->sigma * this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
+ if( XDiffOrder == 1 && YDiffOrder == 1 && ZDiffOrder == 0 )
+ return 4.0 * x * y / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
+ if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 1 )
+ return 4.0 * x * z / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
+ if( XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 1 )
+ return 4.0 * y * z / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
return 0.0;
}
diff --git a/src/functions/tnlSinBumpsFunction_impl.h b/src/functions/tnlSinBumpsFunction_impl.h
index 20a58d2293..0ccf9f8a09 100644
--- a/src/functions/tnlSinBumpsFunction_impl.h
+++ b/src/functions/tnlSinBumpsFunction_impl.h
@@ -152,6 +152,8 @@ getValue( const Vertex& v,
return 2.0 * M_PI / this->waveLength.y() * this->amplitude * cos( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * sin( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() );
if( XDiffOrder == 0 && YDiffOrder == 2 )
return -4.0 * M_PI * M_PI / ( this->waveLength.y() * this->waveLength.y() ) * this->amplitude * sin( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * sin( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() );
+ if( XDiffOrder == 1 && YDiffOrder == 1 )
+ return 4.0 * M_PI * M_PI / ( this->waveLength.x() * this->waveLength.y() ) * this->amplitude * cos( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * cos( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() );
return 0.0;
}
@@ -212,6 +214,12 @@ getValue( const Vertex& v,
return 2.0 * M_PI / this->waveLength.z() * this->amplitude * cos( this->phase.z() + 2.0 * M_PI * z / this->waveLength.z() ) * sin( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * sin( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() );
if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 2)
return -4.0 * M_PI * M_PI / ( this->waveLength.z() * this->waveLength.z() ) * this->amplitude * sin( this->phase.z() + 2.0 * M_PI * z / this->waveLength.z() ) * sin( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * sin( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() );
+ if( XDiffOrder == 1 && YDiffOrder == 1 && ZDiffOrder == 0)
+ return 4.0 * M_PI * M_PI / ( this->waveLength.x() * this->waveLength.y() ) * this->amplitude * cos( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() ) * cos( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * sin( this->phase.z() + 2.0 * M_PI * z / this->waveLength.z() );
+ if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 1)
+ return 4.0 * M_PI * M_PI / ( this->waveLength.x() * this->waveLength.z() ) * this->amplitude * cos( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() ) * sin( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * cos( this->phase.z() + 2.0 * M_PI * z / this->waveLength.z() );
+ if( XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 1)
+ return 4.0 * M_PI * M_PI / ( this->waveLength.y() * this->waveLength.z() ) * this->amplitude * sin( this->phase.x() + 2.0 * M_PI * x / this->waveLength.x() ) * cos( this->phase.y() + 2.0 * M_PI * y / this->waveLength.y() ) * cos( this->phase.z() + 2.0 * M_PI * z / this->waveLength.z() );
return 0.0;
}
diff --git a/src/functions/tnlSinWaveFunction_impl.h b/src/functions/tnlSinWaveFunction_impl.h
index 49c671af5d..50a3683701 100644
--- a/src/functions/tnlSinWaveFunction_impl.h
+++ b/src/functions/tnlSinWaveFunction_impl.h
@@ -140,6 +140,8 @@ getValue( const Vertex& v,
return 2.0 * M_PI * y / ( this->waveLength * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
if( XDiffOrder == 0 && YDiffOrder == 2 )
return 2.0 * M_PI * y * y / ( this->waveLength * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) - 4.0 * M_PI * M_PI * y * y / ( this->waveLength * this->waveLength * ( x * x + y * y ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
+ if( XDiffOrder == 1 && YDiffOrder == 1 )
+ return -4.0 * M_PI * M_PI * x * y / ( this->waveLength * this->waveLength * x * x + y * y ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
return 0.0;
}
@@ -175,6 +177,12 @@ getValue( const Vertex& v,
return 2.0 * M_PI * z / ( this->waveLength * sqrt( x * x + y * y + z * z ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 2 )
return 2.0 * M_PI * ( x * x + y * y ) / ( this->waveLength * sqrt( x * x + y * y + z * z ) * sqrt( x * x + y * y + z * z ) * sqrt( x * x + y * y + z * z ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength ) - 4.0 * M_PI * M_PI * z * z / ( this->waveLength * this->waveLength * ( x * x + y * y + z * z ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
+ if( XDiffOrder == 1 && YDiffOrder == 1 && ZDiffOrder == 0 )
+ return -4.0 * M_PI * M_PI * x * y / ( this->waveLength * this->waveLength * x * x + y * y + z * z ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
+ if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 1 )
+ return -4.0 * M_PI * M_PI * x * z / ( this->waveLength * this->waveLength * x * x + y * y + z * z ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
+ if( XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 1 )
+ return -4.0 * M_PI * M_PI * y * z / ( this->waveLength * this->waveLength * x * x + y * y + z * z ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y + z * z ) / this->waveLength );
return 0.0;
}
diff --git a/src/operators/CMakeLists.txt b/src/operators/CMakeLists.txt
index 8d6c37b489..00448fe209 100755
--- a/src/operators/CMakeLists.txt
+++ b/src/operators/CMakeLists.txt
@@ -1,6 +1,8 @@
ADD_SUBDIRECTORY( gradient )
ADD_SUBDIRECTORY( diffusion )
ADD_SUBDIRECTORY( euler )
+ADD_SUBDIRECTORY( operator-Q )
+ADD_SUBDIRECTORY( operator-curvature )
SET( headers tnlFiniteDifferences.h
tnlFiniteDifferences_impl.h
@@ -32,4 +34,4 @@ set( tnl_operators_SOURCES
${common_SOURCES}
PARENT_SCOPE )
-INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators )
\ No newline at end of file
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators )
diff --git a/src/operators/diffusion/CMakeLists.txt b/src/operators/diffusion/CMakeLists.txt
index cd9601b16a..d772eef4bb 100755
--- a/src/operators/diffusion/CMakeLists.txt
+++ b/src/operators/diffusion/CMakeLists.txt
@@ -1,7 +1,13 @@
+ADD_SUBDIRECTORY( nonlinear-diffusion-operators )
+
SET( headers tnlLinearDiffusion.h
tnlLinearDiffusion_impl.h
tnlExactLinearDiffusion.h
- tnlExactLinearDiffusion_impl.h )
+ tnlExactLinearDiffusion_impl.h
+ tnlNonLinearDiffusion.h
+ tnlNonLinearDiffusion_impl.h
+ tnlExactNonLinearDiffusion.h
+ tnlExactNonLinearDiffusion_impl.h )
SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/operators/diffusion )
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt b/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt
new file mode 100755
index 0000000000..0962ac61ee
--- /dev/null
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt
@@ -0,0 +1,4 @@
+SET( headers tnlOneSideDiffNonlinearOperator.h
+ tnlOneSideDiffNonlinearOperator_impl.h )
+
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/diffusion/nonlinear-diffusion-operators )
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h
new file mode 100644
index 0000000000..068b39f682
--- /dev/null
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h
@@ -0,0 +1,187 @@
+#ifndef TNLONESIDEDIFFNONLINEAROPERATOR_H
+#define TNLONESIDEDIFFNONLINEAROPERATOR_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename NonlinearDiffusionOperator,
+ typename OperatorQ,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class tnlOneSideDiffNonlinearOperator
+{
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+class tnlOneSideDiffNonlinearOperator< tnlGrid< 1,MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef OperatorQ OperatorQType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ OperatorQ operatorQ;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+class tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef OperatorQ OperatorQType;
+
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ OperatorQ operatorQ;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+class tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef OperatorQ OperatorQType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ OperatorQ operatorQ;
+};
+
+
+#include "tnlOneSideDiffNonlinearOperator_impl.h"
+
+
+#endif /* TNLONESIDEDIFFNONLINEAROPERATOR_H */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h~
similarity index 73%
rename from examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h
rename to src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h~
index 2800f91511..041307fa9a 100644
--- a/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h~
@@ -1,16 +1,17 @@
-#ifndef TNLMEANCURVATIVEFLOWDIFFUSION_H
-#define TNLMEANCURVATIVEFLOWDIFFUSION_H
+#ifndef TNLONESIDEDIFFNONLINEAROPERATOR_H
+#define TNLONESIDEDIFFNONLINEAROPERATOR_H
#include <core/vectors/tnlVector.h>
#include "tnlForwardFiniteDifference.h"
#include "tnlBackwardFiniteDifference.h"
-#include "tnlOperatorQ.h"
#include <mesh/tnlGrid.h>
template< typename Mesh,
+ typename NonlinearDiffusionOperator,
+ typename OperatorQ,
typename Real = typename Mesh::RealType,
typename Index = typename Mesh::IndexType >
-class tnlMeanCurvatureFlowDiffusion
+class tnlOneSideDiffNonlinearOperator
{
};
@@ -20,8 +21,9 @@ template< typename MeshReal,
typename Device,
typename MeshIndex,
typename Real,
- typename Index >
-class tnlMeanCurvatureFlowDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+ typename Index,
+ typename OperatorQ >
+class tnlOneSideDiffNonlinearOperator< tnlGrid< 1,MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
{
public:
@@ -30,6 +32,7 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, R
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
+ typedef OperatorQ OperatorQType;
static tnlString getType();
@@ -41,19 +44,8 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, R
const IndexType cellIndex,
const CoordinatesType& coordinates,
const Vector& u,
- const RealType& time) const
- {}
-
-// void setupDofs ( const MeshType& mesh )
-// {}
-//
-// template< typename Vector >
-// void computeFirstGradient( const MeshType& mesh,
-// const RealType& time,
-// const Vector u)
-// {}
+ const RealType& time) const;
-
#ifdef HAVE_CUDA
__device__ __host__
#endif
@@ -76,6 +68,10 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, R
MatrixRow& matrixRow ) const
{}
+
+ private:
+
+ OperatorQ operatorQ;
};
@@ -83,8 +79,9 @@ template< typename MeshReal,
typename Device,
typename MeshIndex,
typename Real,
- typename Index >
-class tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+ typename Index,
+ typename OperatorQ >
+class tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
{
public:
@@ -93,8 +90,9 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >,
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
-// typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+ typedej OperatorQ OperatorQType:
+
static tnlString getType();
template< typename Vector >
@@ -107,13 +105,6 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >,
const Vector& u,
const RealType& time) const;
-// void setupDofs ( const MeshType& mesh );
-//
-// template< typename Vector >
-// void computeFirstGradient( const MeshType& mesh,
-// const RealType& time,
-// const Vector u);
-
#ifdef HAVE_CUDA
__device__ __host__
#endif
@@ -135,12 +126,10 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >,
Vector& b,
MatrixRow& matrixRow ) const
{}
-
- protected:
-
- tnlForwardFiniteDifference<MeshType> fDifference;
- tnlBackwardFiniteDifference<MeshType> bDifference;
- tnlOperatorQ<MeshType> operatorQ;
+
+ private:
+
+ OperatorQ operatorQ;
};
@@ -148,8 +137,9 @@ template< typename MeshReal,
typename Device,
typename MeshIndex,
typename Real,
- typename Index >
-class tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+ typename Index,
+ typename OperatorQ >
+class tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
{
public:
@@ -158,6 +148,7 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >,
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
+ typedej OperatorQ OperatorQType;
static tnlString getType();
@@ -172,14 +163,6 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >,
const RealType& time) const
{}
-// void setupDofs ( const MeshType& mesh )
-// {}
-//
-// template< typename Vector >
-// void computeFirstGradient( const MeshType& mesh,
-// const RealType& time,
-// const Vector u )
-// {}
#ifdef HAVE_CUDA
__device__ __host__
#endif
@@ -205,7 +188,7 @@ class tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >,
};
-#include "tnlMeanCurvatureFlowDiffusion_impl.h"
+#include "tnlOneSideDiffNonlinearOperator_impl.h"
-#endif /* TNLMEANCURVATIVEFLOWDIFFUSION_H */
+#endif /* TNLONESIDEDIFFNONLINEAROPERATOR_H */
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h
new file mode 100644
index 0000000000..095f1d437a
--- /dev/null
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h
@@ -0,0 +1,320 @@
+
+#ifndef TNLONESIDEDIFFNONLINEAROPERATOR_IMPL_H
+#define TNLONESIDEDIFFNONLINEAROPERATOR_IMPL_H
+
+#include "tnlOneSideDiffNonlinearOperator.h"
+
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlOneSideDiffNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return operatorQ.getValueStriped(mesh,cellIndex, coordinates, u, time)*((( u[ mesh.template getCellNextToCell< 1 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHxInverse() / operatorQ.getValue(mesh,cellIndex, coordinates, u, time) - ( - u[ mesh.template getCellNextToCell< -1>( cellIndex ) ]
+ + u[ cellIndex ] ) * mesh.getHxInverse() / operatorQ.getValue(mesh,mesh.template getCellNextToCell<-1>(cellIndex), coordinates, u, time))
+ *mesh.getHxInverse());
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlOneSideDiffNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 3;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlOneSideDiffNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType aCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1 >( index ), coordinates, u, time );
+ const RealType bCoef = 1 + tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time ) + mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1 >( index ), coordinates, u, time ) );
+ const RealType cCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ matrixRow.setElement( 0, mesh.template getCellNextToCell< -1 >( index ), aCoef );
+ matrixRow.setElement( 1, index, bCoef );
+ matrixRow.setElement( 2, mesh.template getCellNextToCell< 1 >( index ), cCoef );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return operatorQ.getValueStriped(mesh,cellIndex, coordinates, u, time)*(((u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHxInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time)
+ -( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ]+ u[ cellIndex ]) * mesh.getHxInverse()/
+ operatorQ.getValue(mesh,mesh.template getCellNextToCell<-1,0>(cellIndex), coordinates, u, time))
+ *mesh.getHxInverse()+(( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHyInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time)
+ -( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHyInverse()
+ /operatorQ.getValue(mesh,mesh.template getCellNextToCell<0,-1>(cellIndex),coordinates, u, time))*mesh.getHyInverse());
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 5;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType aCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1 >( index ), coordinates, u, time );
+ const RealType bCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1,0 >( index ), coordinates, u, time );
+ const RealType cCoef = 1 + tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time ) + mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1 >( index ), coordinates, u, time )
+ + mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time ) +
+ mesh.getHxSquareInverse() / operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1,0 >( index ), coordinates, u, time ) );
+ const RealType dCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ const RealType eCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,-1 >( index ), aCoef );
+ matrixRow.setElement( 1, mesh.template getCellNextToCell< -1,0 >( index ), bCoef );
+ matrixRow.setElement( 2, index , cCoef );
+ matrixRow.setElement( 3, mesh.template getCellNextToCell< 1,0 >( index ), dCoef );
+ matrixRow.setElement( 4, mesh.template getCellNextToCell< 0,1 >( index ), eCoef );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return operatorQ.getValueStriped(mesh,cellIndex, coordinates, u, time)*(((u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHxInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time)
+ -( - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ]+ u[ cellIndex ]) * mesh.getHxInverse()/
+ operatorQ.getValue(mesh,mesh.template getCellNextToCell<-1,0,0>(cellIndex), coordinates, u, time))
+ *mesh.getHxInverse()+(( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHyInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time)
+ -( - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHyInverse()
+ /operatorQ.getValue(mesh,mesh.template getCellNextToCell<0,-1,0>(cellIndex),coordinates, u, time))*mesh.getHyInverse()
+ +(( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHzInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time)
+ -( - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHzInverse()
+ /operatorQ.getValue(mesh,mesh.template getCellNextToCell<0,0,-1>(cellIndex),coordinates, u, time))*mesh.getHzInverse());
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 7;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType aCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHzSquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,0,-1 >( index ), coordinates, u, time );
+ const RealType bCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1,0 >( index ), coordinates, u, time );
+ const RealType cCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1,0,0 >( index ), coordinates, u, time );
+ const RealType dCoef = 1 + tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time ) + mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1,0 >( index ), coordinates, u, time )
+ + mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time ) +
+ mesh.getHxSquareInverse() / operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1,0,0 >( index ), coordinates, u, time )
+ + mesh.getHzSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time ) +
+ mesh.getHzSquareInverse() / operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,0,-1 >( index ), coordinates, u, time ) );
+ const RealType eCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ const RealType fCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ const RealType gCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHzSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,0,-1 >( index ), aCoef );
+ matrixRow.setElement( 1, mesh.template getCellNextToCell< 0,-1,0 >( index ), bCoef );
+ matrixRow.setElement( 2, mesh.template getCellNextToCell< -1,0,0 >( index ), cCoef );
+ matrixRow.setElement( 3, index, dCoef );
+ matrixRow.setElement( 4, mesh.template getCellNextToCell< 1,0,0 >( index ), eCoef );
+ matrixRow.setElement( 5, mesh.template getCellNextToCell< 0,1,0 >( index ), fCoef );
+ matrixRow.setElement( 6, mesh.template getCellNextToCell< 0,0,1 >( index ), gCoef );
+}
+#endif /* TNLONESIDEDIFFNONLINEAROPERATOR_IMPL_H */
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h~ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h~
new file mode 100644
index 0000000000..1bbd523ee0
--- /dev/null
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h~
@@ -0,0 +1,84 @@
+
+#ifndef TNLONESIDEDIFFNONLINEAROPERATOR_IMPL_H
+#define TNLONESIDEDIFFNONLINEAROPERATOR_IMPL_H
+
+#include "tnlOneSideDiffNonlinearOperator.h"
+
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlOneSideDiffNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return operatorQ.getValueStriped(mesh,cellIndex,u)*((fDifference.getValueX(mesh,cellIndex,u)/operatorQ.getValue(mesh,cellIndex,u)
+ -bDifference.getValueX(mesh,cellIndex,u)/operatorQ.getValue(mesh,mesh.template getCellNextToCell<-1,0>(cellIndex),u))
+ *mesh.getHxInverse()+(fDifference.getValueY(mesh,cellIndex,u)/operatorQ.getValue(mesh,cellIndex,u)
+ -bDifference.getValueY(mesh,cellIndex,u)/operatorQ.getValue(mesh,mesh.template getCellNextToCell<0,-1>(cellIndex),u))*mesh.getHyInverse());
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+#endif /* TNLONESIDEDIFFNONLINEAROPERATOR_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlExactFlowDiffusion.h b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
similarity index 77%
rename from examples/mean-curvature-flow/tnlExactFlowDiffusion.h
rename to src/operators/diffusion/tnlExactNonlinearDiffusion.h
index 43d206f255..07da1dc96c 100644
--- a/examples/mean-curvature-flow/tnlExactFlowDiffusion.h
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
@@ -15,17 +15,17 @@
* *
***************************************************************************/
-#ifndef TNLEXACTFLOWDIFFUSION_H_
-#define TNLEXACTFLOWDIFFUSION_H_
+#ifndef TNLEXACTNONLINEARDIFFUSION_H_
+#define TNLEXACTNONLINEARDIFFUSION_H_
#include <functions/tnlFunctionType.h>
-template< int Dimensions >
-class tnlExactFlowDiffusion
+template< typename OperatorQ, int Dimensions >
+class tnlExactNonlinearDiffusion
{};
-template<>
-class tnlExactFlowDiffusion< 1 >
+template< typename OperatorQ >
+class tnlExactNonlinearDiffusion< OperatorQ, 1 >
{
public:
@@ -43,14 +43,11 @@ class tnlExactFlowDiffusion< 1 >
#endif
static Real getValue( const Function& function,
const Vertex& v,
- const Real& time = 0.0 )
- {
- return 0;
- }
+ const Real& time = 0.0 );
};
-template<>
-class tnlExactFlowDiffusion< 2 >
+template< typename OperatorQ >
+class tnlExactNonlinearDiffusion< OperatorQ, 2 >
{
public:
@@ -62,7 +59,7 @@ class tnlExactFlowDiffusion< 2 >
template< typename Function, typename Vertex, typename Real >
#else
template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
-#endif
+#endif
#ifdef HAVE_CUDA
__device__ __host__
#endif
@@ -71,8 +68,8 @@ class tnlExactFlowDiffusion< 2 >
const Real& time = 0.0 );
};
-template<>
-class tnlExactFlowDiffusion< 3 >
+template< typename OperatorQ >
+class tnlExactNonlinearDiffusion< OperatorQ, 3 >
{
public:
@@ -84,25 +81,22 @@ class tnlExactFlowDiffusion< 3 >
template< typename Function, typename Vertex, typename Real >
#else
template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
-#endif
+#endif
#ifdef HAVE_CUDA
__device__ __host__
#endif
static Real getValue( const Function& function,
const Vertex& v,
- const Real& time = 0.0 )
- {
- return 0;
- }
+ const Real& time = 0.0 );
};
-template< int Dimensions >
-class tnlFunctionType< tnlExactFlowDiffusion< Dimensions > >
+template< typename OperatorQ, int Dimensions >
+class tnlFunctionType< tnlExactNonlinearDiffusion< OperatorQ, Dimensions > >
{
public:
enum { Type = tnlAnalyticFunction };
};
-#include "tnlExactFlowDiffusion_impl.h"
+#include "tnlExactNonlinearDiffusion_impl.h"
-#endif /* TNLEXACTFLOWDIFFUSION_H_ */
+#endif /* TNLEXACTNONLINEARDIFFUSION_H_ */
diff --git a/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h b/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
new file mode 100644
index 0000000000..b1554d13d4
--- /dev/null
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
@@ -0,0 +1,96 @@
+/***************************************************************************
+ tnlExactLinearDiffusion_impl.h - description
+ -------------------
+ begin : Aug 8, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLEXACTNONLINEARDIFFUSION_IMPL_H_
+#define TNLEXACTNONLINEARDIFFUSION_IMPL_H_
+
+template< typename OperatorQ >
+tnlString
+tnlExactNonlinearDiffusion< OperatorQ, 1 >::
+getType()
+{
+ return "tnlExactNonlinearDiffusion< " + OperatorQ::getType() + ", 1 >";
+}
+
+template< typename OperatorQ >
+template< typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactNonlinearDiffusion< OperatorQ, 1 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time )
+{
+ return function.template getValue< 2, 0, 0, Vertex >( v, time ) - function.template getValue< 1, 0, 0, Vertex >( v, time )
+ * OperatorQ::template getValue<1, 0, 0>(function, v, time ) / OperatorQ::template getValue<0, 0, 0>(function, v, time );
+}
+
+template< typename OperatorQ >
+tnlString
+tnlExactNonlinearDiffusion< OperatorQ, 2 >::
+getType()
+{
+ return "tnlExactNonlinearDiffusion< " + OperatorQ::getType() + ", 2 >";
+}
+
+template< typename OperatorQ >
+template< typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactNonlinearDiffusion< OperatorQ, 2 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time )
+{
+ return function.template getValue< 2, 0, 0, Vertex >( v, time ) + function.template getValue< 0, 2, 0, Vertex >( v, time )
+ -( OperatorQ::template getValue<1, 0, 0> (function, v, time) * function.template getValue< 1, 0, 0, Vertex >( v, time )
+ + OperatorQ::template getValue<0, 1, 0> (function, v, time) * function.template getValue< 0, 1, 0, Vertex >( v, time ) )
+ / OperatorQ::template getValue<0, 0, 0> (function, v, time);
+}
+
+template< typename OperatorQ >
+tnlString
+tnlExactNonlinearDiffusion< OperatorQ, 3 >::
+getType()
+{
+ return "tnlExactNonlinearDiffusion< " + OperatorQ::getType() + ", 3 >";
+}
+
+template< typename OperatorQ >
+template< typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactNonlinearDiffusion< OperatorQ, 3 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time )
+{
+ return function.template getValue< 2, 0, 0, Vertex >( v, time ) + function.template getValue< 0, 2, 0, Vertex >( v, time )
+ + function.template getValue< 0, 0, 2, Vertex >( v, time )
+ -( OperatorQ::template getValue<1, 0, 0> (function, v, time) * function.template getValue< 1, 0, 0, Vertex >( v, time )
+ + OperatorQ::template getValue<0, 1, 0> (function, v, time) * function.template getValue< 0, 1, 0, Vertex >( v, time )
+ + OperatorQ::template getValue<0, 0, 1> (function, v, time) * function.template getValue< 0, 0, 1, Vertex >( v, time ) )
+ / OperatorQ::template getValue<0, 0, 0> (function, v, time);
+}
+
+#endif /* TNLEXACTNONLINEARDIFFUSION_IMPL_H_ */
diff --git a/src/operators/diffusion/tnlNonlinearDiffusion.h b/src/operators/diffusion/tnlNonlinearDiffusion.h
new file mode 100644
index 0000000000..5ee922340a
--- /dev/null
+++ b/src/operators/diffusion/tnlNonlinearDiffusion.h
@@ -0,0 +1,187 @@
+#ifndef TNLNONLINEARDIFFUSION_H
+#define TNLNONLINEARDIFFUSION_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename NonlinearDiffusionOperator,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class tnlNonlinearDiffusion
+{
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+class tnlNonlinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef NonlinearDiffusionOperator NonlinearDiffusionOperatorType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ NonlinearDiffusionOperator nonlinearDiffusionOperator;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+class tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef NonlinearDiffusionOperator NonlinearDiffusionOperatorType;
+
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ NonlinearDiffusionOperator nonlinearDiffusionOperator;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+class tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef NonlinearDiffusionOperator NonlinearDiffusionOperatorType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+
+ public:
+
+ NonlinearDiffusionOperator nonlinearDiffusionOperator;
+};
+
+
+#include "tnlNonlinearDiffusion_impl.h"
+
+
+#endif /* TNLNONLINEARDIFFUSION_H */
diff --git a/src/operators/diffusion/tnlNonlinearDiffusion_impl.h b/src/operators/diffusion/tnlNonlinearDiffusion_impl.h
new file mode 100644
index 0000000000..30d8fef981
--- /dev/null
+++ b/src/operators/diffusion/tnlNonlinearDiffusion_impl.h
@@ -0,0 +1,249 @@
+
+#ifndef TNLNONLINEARDIFFUSION_IMPL_H
+#define TNLNONLINEARDIFFUSION_IMPL_H
+
+#include "tnlNonlinearDiffusion.h"
+
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+tnlString
+tnlNonlinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getType()
+{
+ return tnlString( "tnlNonlinearDiffusion< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + "," +
+ NonlinearDiffusionOperator::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlNonlinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return nonlinearDiffusionOperator.getValue( mesh, cellIndex, coordinates, u, time );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlNonlinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return nonlinearDiffusionOperator.getLinearSystemRowLength( mesh, index, coordinates );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlNonlinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ nonlinearDiffusionOperator.updateLinearSystem( time, tau, mesh, index, coordinates, u, b, matrixRow );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+tnlString
+tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getType()
+{
+ return tnlString( "tnlNonlinearDiffusion< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + "," +
+ NonlinearDiffusionOperator::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return nonlinearDiffusionOperator.getValue( mesh, cellIndex, coordinates, u, time );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return nonlinearDiffusionOperator.getLinearSystemRowLength( mesh, index, coordinates );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ nonlinearDiffusionOperator.updateLinearSystem( time, tau, mesh, index, coordinates, u, b, matrixRow );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+tnlString
+tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getType()
+{
+ return tnlString( "tnlNonlinearDiffusion< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + "," +
+ NonlinearDiffusionOperator::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return nonlinearDiffusionOperator.getValue( mesh, cellIndex, coordinates, u, time );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return nonlinearDiffusionOperator.getLinearSystemRowLength( mesh, index, coordinates );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename NonlinearDiffusionOperator >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, NonlinearDiffusionOperator, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ nonlinearDiffusionOperator.updateLinearSystem( time, tau, mesh, index, coordinates, u, b, matrixRow );
+}
+
+#endif /* TNLNONLINEARDIFFUSION_IMPL_H */
diff --git a/src/operators/operator-Q/CMakeLists.txt b/src/operators/operator-Q/CMakeLists.txt
new file mode 100755
index 0000000000..3dcfef67cc
--- /dev/null
+++ b/src/operators/operator-Q/CMakeLists.txt
@@ -0,0 +1,6 @@
+SET( headers tnlExactOperatorQ.h
+ tnlExactOperatorQ_impl.h
+ tnlOneSideDiffOperatorQForGraph.h
+ tnlOneSideDiffOperatorQForGraph_impl.h )
+
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/operator-Q )
diff --git a/src/operators/operator-Q/tnlExactOperatorQ.h b/src/operators/operator-Q/tnlExactOperatorQ.h
new file mode 100644
index 0000000000..e9e731aa75
--- /dev/null
+++ b/src/operators/operator-Q/tnlExactOperatorQ.h
@@ -0,0 +1,90 @@
+#ifndef TNLEXACTOPERATORQ_H
+#define TNLEXACTOPERATORQ_H
+
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlSharedVector.h>
+#include <mesh/tnlGrid.h>
+#include <functions/tnlFunctionType.h>
+
+template< int Dimensions >
+class tnlExactOperatorQ
+{};
+
+template<>
+class tnlExactOperatorQ< 1 >
+{
+ public:
+
+ enum { Dimensions = 1 };
+
+ static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real >
+#else
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ static Real getValue( const Function& function,
+ const Vertex& v,
+ const Real& time = 0.0, const Real& eps = 1.0 );
+
+};
+
+template<>
+class tnlExactOperatorQ< 2 >
+{
+ public:
+
+ enum { Dimensions = 2 };
+
+ static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real >
+#else
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ static Real getValue( const Function& function,
+ const Vertex& v,
+ const Real& time = 0.0, const Real& eps = 1.0 );
+};
+
+template<>
+class tnlExactOperatorQ< 3 >
+{
+ public:
+
+ enum { Dimensions = 3 };
+
+ static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real >
+#else
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ static Real getValue( const Function& function,
+ const Vertex& v,
+ const Real& time = 0.0, const Real& eps = 1.0 );
+};
+
+template< int Dimensions >
+class tnlFunctionType< tnlExactOperatorQ< Dimensions > >
+{
+ public:
+ enum { Type = tnlAnalyticFunction };
+};
+
+#include <operators/operator-Q/tnlExactOperatorQ_impl.h>
+
+
+#endif /* TNLEXACTOPERATORQ_H */
diff --git a/src/operators/operator-Q/tnlExactOperatorQ_impl.h b/src/operators/operator-Q/tnlExactOperatorQ_impl.h
new file mode 100644
index 0000000000..2aaa77f2ec
--- /dev/null
+++ b/src/operators/operator-Q/tnlExactOperatorQ_impl.h
@@ -0,0 +1,116 @@
+/***************************************************************************
+ tnlExactLinearDiffusion_impl.h - description
+ -------------------
+ begin : Aug 8, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLEXACTOPERATORQ_IMPL_H_
+#define TNLEXACTOPERATORQ_IMPL_H_
+
+#include <operators/operator-Q/tnlExactOperatorQ.h>
+
+tnlString
+tnlExactOperatorQ< 1 >::
+getType()
+{
+ return "tnlExactOperatorQ< 1 >";
+}
+
+template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactOperatorQ< 1 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time, const Real& eps )
+{
+ if( YDiffOrder != 0 || ZDiffOrder != 0 )
+ return 0.0;
+ if (XDiffOrder == 0)
+ return sqrt(eps * eps + function.template getValue< 1, 0, 0, Vertex >( v, time ) * function.template getValue< 1, 0, 0, Vertex >( v, time ) );
+ if (XDiffOrder == 1)
+ return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 2, 0, 0 >( v, time ) ) / getValue( function, v, time, eps);
+ return 0;
+}
+
+tnlString
+tnlExactOperatorQ< 2 >::
+getType()
+{
+ return "tnlExactOperatorQ< 2 >";
+}
+
+template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactOperatorQ< 2 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time, const Real& eps )
+{
+ if( ZDiffOrder != 0 )
+ return 0.0;
+ if (XDiffOrder == 0 && YDiffOrder == 0 )
+ return sqrt(eps * eps + function.template getValue< 1, 0, 0, Vertex >( v, time ) * function.template getValue< 1, 0, 0, Vertex >( v, time )
+ + function.template getValue< 0, 1, 0, Vertex >( v, time ) * function.template getValue< 0, 1, 0, Vertex >( v, time ) );
+ if (XDiffOrder == 1 && YDiffOrder == 0 )
+ return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 2, 0, 0 >( v, time ) +
+ function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time )) / getValue( function, v, time, eps);
+ if (XDiffOrder == 0 && YDiffOrder == 1 )
+ return (function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 2, 0 >( v, time ) +
+ function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time )) / getValue( function, v, time, eps);
+ return 0;
+}
+
+tnlString
+tnlExactOperatorQ< 3 >::
+getType()
+{
+ return "tnlExactOperatorQ< 3 >";
+}
+
+template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactOperatorQ< 3 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time, const Real& eps )
+{
+ if ( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
+ return sqrt(eps * eps + function.template getValue< 1, 0, 0, Vertex >( v, time ) * function.template getValue< 1, 0, 0, Vertex >( v, time )
+ + function.template getValue< 0, 1, 0, Vertex >( v, time ) * function.template getValue< 0, 1, 0, Vertex >( v, time )
+ + function.template getValue< 0, 0, 1, Vertex >( v, time ) * function.template getValue< 0, 0, 1, Vertex >( v, time ) );
+ if (XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 0 )
+ return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 2, 0, 0 >( v, time ) +
+ function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time ) +
+ function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 1, 0, 1 >( v, time )) / getValue( function, v, time, eps);
+ if (XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 0 )
+ return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time ) +
+ function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 2, 0 >( v, time ) +
+ function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 0, 1, 1 >( v, time )) / getValue( function, v, time, eps);
+ if (XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 1 )
+ return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 0, 1 >( v, time ) +
+ function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 1, 1 >( v, time ) +
+ function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 0, 0, 2 >( v, time )) / getValue( function, v, time, eps);
+ return 0;
+}
+
+#endif /* TNLEXACTOPERATORQ_IMPL_H_ */
diff --git a/src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h b/src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h
new file mode 100644
index 0000000000..a7629b2ca9
--- /dev/null
+++ b/src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h
@@ -0,0 +1,327 @@
+#ifndef TNLONESIDEDIFFOPERATORQFORGRAPH_H
+#define TNLONESIDEDIFFOPERATORQFORGRAPH_H
+
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlSharedVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType,
+ int Precomputation = 0 >
+class tnlOneSideDiffOperatorQForGraph
+{
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlOneSideDiffOperatorQForGraph< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index, 0 >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ IndexType bind( Vector& u)
+ { return 0; }
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time )
+ {}
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time )const;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlOneSideDiffOperatorQForGraph< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index, 0 >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ IndexType bind( Vector& u)
+ { return 0; }
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time )
+ {}
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time )const;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlOneSideDiffOperatorQForGraph< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index, 0 >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ IndexType bind( Vector& u)
+ { return 0; }
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time )
+ {}
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlOneSideDiffOperatorQForGraph< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index, 1 >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ Index bind( Vector& u);
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time );
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ private:
+
+ tnlSharedVector< RealType, DeviceType, IndexType > u;
+ tnlVector< RealType, DeviceType, IndexType> q;
+ tnlVector< RealType, DeviceType, IndexType> qStriped;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlOneSideDiffOperatorQForGraph< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index, 1 >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlSharedVector< RealType, DeviceType, IndexType > DofVectorType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ Index bind( Vector& u);
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time );
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time )const;
+
+ private:
+
+ tnlSharedVector< RealType, DeviceType, IndexType > u;
+ tnlVector< RealType, DeviceType, IndexType> q;
+ tnlVector< RealType, DeviceType, IndexType> qStriped;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlOneSideDiffOperatorQForGraph< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index, 1 >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ Index bind( Vector& u);
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time );
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time )const;
+
+ private:
+
+ tnlSharedVector< RealType, DeviceType, IndexType > u;
+ tnlVector< RealType, DeviceType, IndexType> q;
+ tnlVector< RealType, DeviceType, IndexType> qStriped;
+};
+
+#include <operators/operator-Q/tnlOneSideDiffOperatorQForGraph_impl.h>
+
+
+#endif /* TNLONESIDEDIFFOPERATORQFORGRAPH_H */
diff --git a/src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph_impl.h b/src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph_impl.h
new file mode 100644
index 0000000000..925c82b403
--- /dev/null
+++ b/src/operators/operator-Q/tnlOneSideDiffOperatorQForGraph_impl.h
@@ -0,0 +1,450 @@
+
+#ifndef TNLONESIDEDIFFOPERATORQFORGRAPH_IMPL_H
+#define TNLONESIDEDIFFOPERATORQFORGRAPH_IMPL_H
+
+#include <operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h>
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffOperatorQForGraph< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 0 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffOperatorQForGraph< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 1 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+Index
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+bind( Vector& u )
+{
+ this->u.bind(u);
+ if(q.setSize(u.getSize()))
+ return 1;
+ if(qStriped.setSize(u.getSize()))
+ return 1;
+ q.setValue(0);
+ qStriped.setValue(0);
+ return 0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+void
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+update( const MeshType& mesh, const RealType& time )
+{
+ CoordinatesType dimensions = mesh.getDimensions();
+ CoordinatesType coordinates;
+ for( coordinates.x()=1; coordinates.x() < dimensions.x()-1; coordinates.x()++ )
+ {
+ q.setElement( coordinates.x(), getValue( mesh, mesh.getCellIndex(coordinates), coordinates, u, time ) );
+ }
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return sqrt( 1.0 + ( u[ mesh.template getCellNextToCell< 1 >( cellIndex ) ] - u[ cellIndex ]) *
+ ( u[ mesh.template getCellNextToCell< 1 >( cellIndex ) ] - u[ cellIndex ]) *
+ mesh.getHxInverse() * mesh.getHxInverse() );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return q.getElement(cellIndex);
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return sqrt( 1.0 + 0.5*( ( u[ mesh.template getCellNextToCell< 1 >( cellIndex ) ] - u[ cellIndex ]) *
+ ( u[ mesh.template getCellNextToCell< 1 >( cellIndex ) ] - u[ cellIndex ]) *
+ mesh.getHxInverse() * mesh.getHxInverse() + ( - u[ mesh.template getCellNextToCell< -1 >( cellIndex ) ] + u[ cellIndex ] )
+ * ( - u[ mesh.template getCellNextToCell< -1 >( cellIndex ) ] + u[ cellIndex ] ) * mesh.getHxInverse() * mesh.getHxInverse() ) );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return qStriped.getElement(cellIndex);
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffOperatorQForGraph< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 0 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffOperatorQForGraph< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 1 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+Index
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+bind( Vector& u)
+{
+ this->u.bind(u);
+ if(q.setSize(u.getSize()))
+ return 1;
+ if(qStriped.setSize(u.getSize()))
+ return 1;
+ q.setValue(0);
+ qStriped.setValue(0);
+ return 0;
+}
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+void
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+update( const MeshType& mesh, const RealType& time )
+{
+ CoordinatesType dimensions = mesh.getDimensions();
+ CoordinatesType coordinates;
+
+ for( coordinates.x()=1; coordinates.x() < dimensions.x()-1; coordinates.x()++ )
+ for( coordinates.y()=1; coordinates.y() < dimensions.y()-1; coordinates.y()++ )
+ {
+ q.setElement( mesh.getCellIndex(coordinates), getValue( mesh, mesh.getCellIndex(coordinates), coordinates, u, time ) );
+ }
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return sqrt( 1.0 + ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] ) *
+ ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHyInverse() * mesh.getHyInverse() + ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHxInverse() * mesh.getHxInverse() );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return q.getElement(cellIndex);
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return sqrt( 1.0 + 0.5*( ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] ) *
+ ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHyInverse() * mesh.getHyInverse() + ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHxInverse() * mesh.getHxInverse()
+ + ( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ] + u[ cellIndex ]) * ( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ] + u[ cellIndex ])
+ * mesh.getHxInverse() * mesh.getHxInverse()
+ + ( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] + u[ cellIndex ]) * ( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] + u[ cellIndex ]) *
+ mesh.getHyInverse() * mesh.getHyInverse() ) );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return qStriped.getElement(cellIndex);
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffOperatorQForGraph< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 0 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getType()
+{
+ return tnlString( "tnlOneSideDiffOperatorQForGraph< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 1 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+Index
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+bind( Vector& u)
+{
+ this->u.bind(u);
+ if(q.setSize(u.getSize()))
+ return 1;
+ if(qStriped.setSize(u.getSize()))
+ return 1;
+ q.setValue(0);
+ qStriped.setValue(0);
+ return 0;
+}
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+void
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+update( const MeshType& mesh, const RealType& time )
+{
+ CoordinatesType dimensions = mesh.getDimensions();
+ CoordinatesType coordinates;
+
+ for( coordinates.x()=1; coordinates.x() < dimensions.x()-1; coordinates.x()++ )
+ for( coordinates.y()=1; coordinates.y() < dimensions.y()-1; coordinates.y()++ )
+ for( coordinates.z()=1; coordinates.z() < dimensions.z()-1; coordinates.z()++ )
+ q.setElement( mesh.getCellIndex(coordinates), getValue( mesh, mesh.getCellIndex(coordinates), coordinates, u, time ) );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return sqrt( 1.0 + ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] ) *
+ ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHyInverse() * mesh.getHyInverse() + ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHxInverse() * mesh.getHxInverse()
+ + ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] ) * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHzInverse() * mesh.getHzInverse() );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOneSideDiffOperatorQForGraph< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValueStriped( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return sqrt( 1.0 + 0.5*( ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] ) *
+ ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHyInverse() * mesh.getHyInverse()
+ + ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] ) * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHxInverse() * mesh.getHxInverse() + ( - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] + u[ cellIndex ])
+ * ( - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHxInverse() * mesh.getHxInverse()
+ + ( - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] + u[ cellIndex ]) *
+ ( - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHyInverse() * mesh.getHyInverse()
+ + ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] ) * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHzInverse() * mesh.getHzInverse() + ( - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] + u[ cellIndex ]) *
+ ( - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHzInverse() * mesh.getHzInverse()
+ ) );
+}
+#endif /* TNLONESIDEDIFFOPERATORQFORGRAPH_IMPL_H */
diff --git a/src/operators/operator-curvature/CMakeLists.txt b/src/operators/operator-curvature/CMakeLists.txt
new file mode 100755
index 0000000000..7a877b91de
--- /dev/null
+++ b/src/operators/operator-curvature/CMakeLists.txt
@@ -0,0 +1,4 @@
+SET( headers tnlExactOperatorCurvature.h
+ tnlExactOperatorCurvature_impl.h )
+
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/operator-curvature )
diff --git a/src/operators/operator-curvature/tnlExactOperatorCurvature.h b/src/operators/operator-curvature/tnlExactOperatorCurvature.h
new file mode 100644
index 0000000000..634f93e613
--- /dev/null
+++ b/src/operators/operator-curvature/tnlExactOperatorCurvature.h
@@ -0,0 +1,93 @@
+#ifndef TNLEXACTOPERATORCURVATURE_H
+#define TNLEXACTOPERATORCURVATURE_H
+
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlSharedVector.h>
+#include <mesh/tnlGrid.h>
+#include <functions/tnlFunctionType.h>
+
+template< typename ExactOperatorQ, int Dimensions >
+class tnlExactOperatorCurvature
+{};
+
+template< typename ExactOperatorQ >
+class tnlExactOperatorCurvature< OperatorQ, 1 >
+{
+ public:
+
+ enum { Dimensions = 1 };
+
+ static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real >
+#else
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ static Real getValue( const Function& function,
+ const Vertex& v,
+ const Real& time = 0.0, const Real& eps = 1.0 );
+
+};
+
+template< typename ExactOperatorQ >
+class tnlExactOperatorCurvature< ExactOperatorQ, 2 >
+{
+ public:
+
+ enum { Dimensions = 2 };
+
+ static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real >
+#else
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ static Real getValue( const Function& function,
+ const Vertex& v,
+ const Real& time = 0.0, const Real& eps = 1.0 );
+};
+
+template< typename ExactOperatorQ >
+class tnlExactOperatorCurvature< ExactOperatorQ, 3 >
+{
+ public:
+
+ enum { Dimensions = 3 };
+
+ static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real >
+#else
+ template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ static Real getValue( const Function& function,
+ const Vertex& v,
+ const Real& time = 0.0, const Real& eps = 1.0 )
+ {
+ return 0;
+ }
+};
+
+template< typename ExactOperatorQ, int Dimensions >
+class tnlFunctionType< tnlExactOperatorCurvature< ExactOperatorQ, Dimensions > >
+{
+ public:
+ enum { Type = tnlAnalyticFunction };
+};
+
+#include <operators/operator-curvature/tnlExactOperatorCurvature_impl.h>
+
+
+#endif /* TNLEXACTOPERATORCURVATURE_H */
diff --git a/src/operators/operator-curvature/tnlExactOperatorCurvature_impl.h b/src/operators/operator-curvature/tnlExactOperatorCurvature_impl.h
new file mode 100644
index 0000000000..f09f59bd83
--- /dev/null
+++ b/src/operators/operator-curvature/tnlExactOperatorCurvature_impl.h
@@ -0,0 +1,88 @@
+/***************************************************************************
+ tnlExactLinearDiffusion_impl.h - description
+ -------------------
+ begin : Aug 8, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLEXACTOPERATORCURVATURE_IMPL_H_
+#define TNLEXACTOPERATORCURVATURE_IMPL_H_
+
+#include <operators/operator-curvature/tnlExactOperatorCurvature.h>
+
+template< typename ExactOperatorQ >
+tnlString
+tnlExactOperatorCurvature< ExactOperatorQ, 1 >::
+getType()
+{
+ return "tnlExactOperatorCurvature< " + ExactOperatorQ::getType() + ",1 >";
+}
+
+template< typename OperatorQ >
+template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactOperatorQ< 1 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time, const Real& eps )
+{
+ if( YDiffOrder != 0 || ZDiffOrder != 0 )
+ return 0.0;
+ if (XDiffOrder == 0)
+ return function.template getValue< 2, 0, 0, Vertex >( v, time )/ExactOperatorQ::template getValue< 0, 0, 0 >( function, v, time, eps ) -
+ ( function.template getValue< 1, 0, 0, Vertex >( v, time ) * ExactOperatorQ::template getValue< 1, 0, 0 >( function, v, time, eps ) )
+ / ( ExactOperatorQ::template getValue< 0, 0, 0 >( function, v, time, eps ) * ExactOperatorQ::template getValue< 0, 0, 0 >( function, v, time, eps ) );
+ return 0;
+}
+
+template< typename ExactOperatorQ >
+tnlString
+tnlExactOperatorCurvature< ExactOperatorQ, 2 >::
+getType()
+{
+ return "tnlExactOperatorCurvature< " + ExactOperatorQ::getType() + ",2 >";
+}
+
+template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Real
+tnlExactOperatorQ< 2 >::
+getValue( const Function& function,
+ const Vertex& v,
+ const Real& time, const Real& eps )
+{
+ if( ZDiffOrder != 0 )
+ return 0.0;
+ if (XDiffOrder == 0 && YDiffOrder == 0 )
+ return ( function.template getValue< 2, 0, 0, Vertex >( v, time ) * function.template getValue< 0, 2, 0, Vertex >( v, time ) )
+ /ExactOperatorQ::template getValue< 0, 0, 0 >( function, v, time, eps ) - ( function.template getValue< 1, 0, 0, Vertex >( v, time ) *
+ ExactOperatorQ::template getValue< 1, 0, 0 >( function, v, time, eps ) + function.template getValue< 0, 1, 0, Vertex >( v, time ) *
+ ExactOperatorQ::template getValue< 0, 1, 0 >( function, v, time, eps ) )
+ / ( ExactOperatorQ::template getValue< 0, 0, 0 >( function, v, time, eps ) * ExactOperatorQ::template getValue< 0, 0, 0 >( function, v, time, eps ) );
+ return 0;
+}
+
+template< typename ExactOperatorQ >
+tnlString
+tnlExactOperatorCurvature< ExactOperatorQ, 3 >::
+getType()
+{
+ return "tnlExactOperatorCurvature< " + ExactOperatorQ::getType() + ",3 >";
+}
+
+#endif /* TNLEXACTOPERATORCURVATURE_IMPL_H_ */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem.h b/src/problems/tnlMeanCurvatureFlowEocProblem.h
similarity index 66%
rename from examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem.h
rename to src/problems/tnlMeanCurvatureFlowEocProblem.h
index 5c94ec7100..8960d2306b 100644
--- a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem.h
+++ b/src/problems/tnlMeanCurvatureFlowEocProblem.h
@@ -18,15 +18,17 @@
#ifndef TNLMEANCURVATUREFLOWEOCPROBLEM_H_
#define TNLMEANCURVATUREFLOWEOCPROBLEM_H_
-#include "tnlMeanCurvatureFlowProblem.h"
-
+#include <problems/tnlMeanCurvatureFlowProblem.h>
+#include <operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h>
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
- typename DifferentialOperator = tnlMeanCurvatureFlowDiffusion< Mesh,
- typename BoundaryCondition::RealType > >
-class tnlMeanCurvatureFlowEocProblem : public tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >
+ typename DifferentialOperator = tnlNonlinearDiffusion< Mesh,
+ tnlOneSideDiffNonlinearOperator< Mesh, tnlOneSideDiffOperatorQForGraph<Mesh, typename BoundaryCondition::RealType,
+ typename BoundaryCondition::IndexType, 0>, typename BoundaryCondition::RealType, typename BoundaryCondition::IndexType >,
+ typename BoundaryCondition::RealType, typename BoundaryCondition::IndexType > >
+class tnlMeanCurvatureFlowEocProblem : public tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >
{
public:
@@ -35,6 +37,6 @@ class tnlMeanCurvatureFlowEocProblem : public tnlMeanCurvatureFlowProblem< Mesh,
bool setup( const tnlParameterContainer& parameters );
};
-#include "tnlMeanCurvatureFlowEocProblem_impl.h"
+#include <problems/tnlMeanCurvatureFlowEocProblem_impl.h>
#endif /* TNLMEANCURVATUREFLOWEOCPROBLEM_H_ */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem_impl.h b/src/problems/tnlMeanCurvatureFlowEocProblem_impl.h
similarity index 100%
rename from examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem_impl.h
rename to src/problems/tnlMeanCurvatureFlowEocProblem_impl.h
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocRhs.h b/src/problems/tnlMeanCurvatureFlowEocRhs.h
similarity index 100%
rename from examples/mean-curvature-flow/tnlMeanCurvatureFlowEocRhs.h
rename to src/problems/tnlMeanCurvatureFlowEocRhs.h
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem.h b/src/problems/tnlMeanCurvatureFlowProblem.h
similarity index 78%
rename from examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem.h
rename to src/problems/tnlMeanCurvatureFlowProblem.h
index f2686b2dc2..56904812fc 100644
--- a/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem.h
+++ b/src/problems/tnlMeanCurvatureFlowProblem.h
@@ -18,18 +18,23 @@
#ifndef TNLMEANCURVATUREFLOWPROBLEM_H_
#define TNLMEANCURVATUREFLOWPROBLEM_H_
-#include "tnlMeanCurvatureFlowDiffusion.h"
+#include <operators/diffusion/tnlNonlinearDiffusion.h>
+#include <operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h>
#include <problems/tnlPDEProblem.h>
+#include <operators/operator-Q/tnlOneSideDiffOperatorQForGraph.h>
+#include <matrices/tnlCSRMatrix.h>
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
- typename DifferentialOperator = tnlMeanCurvatureFlowDiffusion< Mesh,
- typename BoundaryCondition::RealType > >
+ typename DifferentialOperator = tnlNonlinearDiffusion< Mesh,
+ tnlOneSideDiffNonlinearOperator< Mesh, tnlOneSideDiffOperatorQForGraph<Mesh, typename Mesh::RealType,
+ typename Mesh::IndexType, 0>, typename Mesh::RealType, typename Mesh::IndexType >,
+ typename Mesh::RealType, typename Mesh::IndexType > >
class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
typename DifferentialOperator::RealType,
typename Mesh::DeviceType,
- typename DifferentialOperator::IndexType >
+ typename DifferentialOperator::IndexType >
{
public:
@@ -37,6 +42,7 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
typedef typename Mesh::DeviceType DeviceType;
typedef typename DifferentialOperator::IndexType IndexType;
typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+ typedef tnlCSRMatrix< RealType, DeviceType, IndexType> MatrixType;
using typename BaseType::MeshType;
using typename BaseType::DofVectorType;
@@ -55,10 +61,10 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
DofVectorType& dofs,
DofVectorType& auxDofs );
- template< typename MatrixType >
+ template< typename Matrix >
bool setupLinearSystem( const MeshType& mesh,
- MatrixType& matrix );
-
+ Matrix& matrix );
+
bool makeSnapshot( const RealType& time,
const IndexType& step,
const MeshType& mesh,
@@ -76,13 +82,13 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
DofVectorType& _u,
DofVectorType& _fu );
- template< typename MatrixType >
+ template< typename Matrix >
void assemblyLinearSystem( const RealType& time,
const RealType& tau,
const MeshType& mesh,
DofVectorType& dofs,
DofVectorType& auxDofs,
- MatrixType& matrix,
+ Matrix& matrix,
DofVectorType& rightHandSide );
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem_impl.h b/src/problems/tnlMeanCurvatureFlowProblem_impl.h
similarity index 80%
rename from examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem_impl.h
rename to src/problems/tnlMeanCurvatureFlowProblem_impl.h
index d03ff41504..8c4d76d143 100644
--- a/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem_impl.h
+++ b/src/problems/tnlMeanCurvatureFlowProblem_impl.h
@@ -99,6 +99,7 @@ bindDofs( const MeshType& mesh,
{
const IndexType dofs = mesh.getNumberOfCells();
this->solution.bind( dofVector.getData(), dofs );
+ differentialOperator.nonlinearDiffusionOperator.operatorQ.bind(solution);
// this->differentialOperator.setupDofs(mesh);
}
@@ -127,27 +128,26 @@ template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
- template< typename MatrixType >
+template< typename Matrix >
bool
tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setupLinearSystem( const MeshType& mesh,
- MatrixType& matrix )
+ Matrix& matrix )
{
-// const IndexType dofs = this->getDofs( mesh );
-// typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;
-// RowLengthsVectorType rowLengths;
-// if( ! rowLengths.setSize( dofs ) )
-// return false;
-// tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;
-// matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
-// differentialOperator,
-// boundaryCondition,
-// rowLengths );
-// matrix.setDimensions( dofs, dofs );
-// if( ! matrix.setRowLengths( rowLengths ) )
-// return false;
-// return true;
-// //return tnlMultidiagonalMatrixSetter< Mesh >::setupMatrix( mesh, matrix );
+ const IndexType dofs = this->getDofs( mesh );
+ typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;
+ RowLengthsVectorType rowLengths;
+ if( ! rowLengths.setSize( dofs ) )
+ return false;
+ tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;
+ matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
+ differentialOperator,
+ boundaryCondition,
+ rowLengths );
+ matrix.setDimensions( dofs, dofs );
+ if( ! matrix.setRowLengths( rowLengths ) )
+ return false;
+ return true;
}
template< typename Mesh,
@@ -198,6 +198,7 @@ getExplicitRHS( const RealType& time,
//cout << "u = " << u << endl;
this->bindDofs( mesh, u );
+ differentialOperator.nonlinearDiffusionOperator.operatorQ.update( mesh, time );
tnlExplicitUpdater< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
explicitUpdater.template update< Mesh::Dimensions >( time,
mesh,
@@ -217,7 +218,7 @@ template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
- template< typename MatrixType >
+template< typename Matrix >
void
tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
assemblyLinearSystem( const RealType& time,
@@ -225,23 +226,23 @@ assemblyLinearSystem( const RealType& time,
const MeshType& mesh,
DofVectorType& u,
DofVectorType& auxDofs,
- MatrixType& matrix,
+ Matrix& matrix,
DofVectorType& b )
{
-// tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
-// systemAssembler.template assembly< Mesh::Dimensions >( time,
-// tau,
-// mesh,
-// this->differentialOperator,
-// this->boundaryCondition,
-// this->rightHandSide,
-// u,
-// matrix,
-// b );
-// /*matrix.print( cout );
-// cout << endl << b << endl;
-// cout << endl << u << endl;
-// abort();*/
+ tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
+ systemAssembler.template assembly< Mesh::Dimensions >( time,
+ tau,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ matrix,
+ b );
+ /*matrix.print( cout );
+ cout << endl << b << endl;
+ cout << endl << u << endl;
+ abort();*/
}
#endif /* TNLMEANCURVATUREFLOWPROBLEM_IMPL_H_ */
--
GitLab