Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Commits on Source (16)
Showing
with 150 additions and 971 deletions
......@@ -60,12 +60,11 @@ stages:
-DWITH_EXAMPLES=${WITH_EXAMPLES}
-DWITH_TOOLS=${WITH_TOOLS}
-DWITH_PYTHON=${WITH_PYTHON}
# - make
# - make test
# "install" implies the "all" target
# - make install
- ninja ${NINJAFLAGS}
# - make test
- ninja ${NINJAFLAGS} install
- ninja test
- ninja install
- popd
variables:
<<: *default_cmake_flags
......@@ -73,15 +72,23 @@ stages:
# Cuda builds are specified first because they take more time than host-only builds,
# which can be allocated on hosts whitout GPUs.
# Similarly, release builds are launched first to avoid the tail effect (they take
# significantly more time than debug builds).
cuda_base_Debug:
cuda_full_Release:
<<: *build_template
tags:
- openmp
- gpu
variables:
<<: *default_cmake_flags
WITH_OPENMP: "yes"
WITH_CUDA: "yes"
BUILD_TYPE: Debug
BUILD_TYPE: Release
WITH_BENCHMARKS: "yes"
WITH_EXAMPLES: "yes"
WITH_TOOLS: "yes"
WITH_PYTHON: "yes"
cuda_base_Release:
<<: *build_template
......@@ -92,19 +99,6 @@ cuda_base_Release:
WITH_CUDA: "yes"
BUILD_TYPE: Release
cuda_mpi_Debug:
<<: *build_template
tags:
- openmp
- gpu
- mpi
variables:
<<: *default_cmake_flags
WITH_OPENMP: "yes"
WITH_CUDA: "yes"
WITH_MPI: "yes"
BUILD_TYPE: Debug
cuda_mpi_Release:
<<: *build_template
tags:
......@@ -133,20 +127,28 @@ cuda_full_Debug:
WITH_TOOLS: "yes"
WITH_PYTHON: "yes"
cuda_full_Release:
cuda_base_Debug:
<<: *build_template
tags:
- gpu
variables:
<<: *default_cmake_flags
WITH_CUDA: "yes"
BUILD_TYPE: Debug
cuda_mpi_Debug:
<<: *build_template
tags:
- openmp
- gpu
- mpi
variables:
<<: *default_cmake_flags
WITH_OPENMP: "yes"
WITH_CUDA: "yes"
BUILD_TYPE: Release
WITH_BENCHMARKS: "yes"
WITH_EXAMPLES: "yes"
WITH_TOOLS: "yes"
WITH_PYTHON: "yes"
WITH_MPI: "yes"
BUILD_TYPE: Debug
default_base_Debug:
<<: *build_template
......
## How to configure git
It is important to [configure your git username and email address](
https://mmg-gitlab.fjfi.cvut.cz/gitlab/help/gitlab-basics/start-using-git.md#add-your-git-username-and-set-your-email),
since every git commit will use this information to identify you as the author:
git config --global user.name "John Doe"
git config --global user.email "john.doe@example.com"
In the TNL project, this username should be your real name (given name and
family name) and the email address should be the email address used in the
Gitlab profile. You should use the same configuration on all computers where
you make commits. If you have made some commits with a different email address
in the past, you can also add secondary email addresses to the Gitlab profile.
## How to write good commit messages
Begin with a short summary line a.k.a. message subject:
- Use up to 50 characters; this is the git official preference.
- Finish without a sentence-ending period.
Continue with a longer description a.k.a. message body:
- Add a blank line after the summary line, then write as much as you want.
- Use up to 72 characters per line for typical text for word wrap.
- Use as many characters as needed for atypical text, such as URLs, terminal
output, formatted messages, etc.
- Include any kind of notes, links, examples, etc. as you want.
See [5 Useful Tips For A Better Commit Message](
https://robots.thoughtbot.com/5-useful-tips-for-a-better-commit-message) and
[How to Write a Git Commit Message](https://chris.beams.io/posts/git-commit/)
for examples and reasoning.
## How to split commits
Extensive changes should be split into multiple commits so that only related
changes are made in each commit. All changes made in a commit should be
described in its commit message. If describing all changes would not result in
a good commit message, you should probably make multiple separate commits.
Multiple small commits are better than one big commit, because later they can be
easily [squashed](https://git-scm.com/book/en/v2/Git-Tools-Rewriting-History)
together, whereas splitting a big commit into logical parts is significantly
more difficult.
> __Tip:__ Use `git add -p` to stage only part of your working tree for the next
> commit. See [git add -p: The most powerful git feature you're not using yet](
https://johnkary.net/blog/git-add-p-the-most-powerful-git-feature-youre-not-using-yet/).
## Rebase-based workflow
The development of new features should follow the rebase-based workflow:
- Create a new _feature_ branch based on the main branch (`develop`).
- Make commits with your work in the feature branch.
- When there are other new commits in the `develop` branch, you should do a
[rebase](https://git-scm.com/book/en/v2/Git-Branching-Rebasing) to rewind your
commits on top of the current develop branch.
- If there are conflicts, you will need to resolve them manually. Hence, it
is a good practice to rebase as often as possible (generally as soon as
the changes appear in the `develop` branch) and to split commits into
logical parts.
- When your work is ready for review, you can open a [merge request](
https://mmg-gitlab.fjfi.cvut.cz/gitlab/tnl/tnl-dev/merge_requests/new).
- If the branch is not ready for merge yet, prepend `[WIP]` to the merge
request title to indicate _work in progress_ and to prevent a premature
merge.
- This is also a good time to squash small commits (e.g. typos, forgotten
changes or trivial corrections) with relevant bigger commits to make the
review easier.
- When reviewed, the feature branch can be merged into the develop branch.
The main advantages of this workflow are linear history, clear commits and
reduction of merge conflicts. See [A rebase-based workflow](
https://brokenco.de/2010/04/02/a-rebase-based-workflow.html) and
[Why is rebase-then-merge better than just merge](https://stackoverflow.com/a/457988)
([complement](https://stackoverflow.com/a/804178)) for reference.
\ No newline at end of file
ADD_SUBDIRECTORY( Config )
ADD_SUBDIRECTORY( Containers )
ADD_SUBDIRECTORY( DistributedContainers )
ADD_SUBDIRECTORY( Communicators )
ADD_SUBDIRECTORY( Debugging )
ADD_SUBDIRECTORY( Devices )
......
......@@ -38,7 +38,7 @@ List< T >::~List()
template< typename T >
String List< T >::getType()
{
return String( "List< " ) + TNL::getType< T >() + String( " >" );
return String( "Containers::List< " ) + TNL::getType< T >() + String( " >" );
}
template< typename T >
......
......@@ -23,7 +23,7 @@ Timer Cuda::smartPointersSynchronizationTimer;
String Cuda::getDeviceType()
{
return String( "Cuda" );
return String( "Devices::Cuda" );
}
int Cuda::getNumberOfBlocks( const int threads,
......
......@@ -71,7 +71,7 @@ class MIC
static String getDeviceType()
{
return String( "MIC" );
return String( "Devices::MIC" );
};
#ifdef HAVE_MIC
......
SET( headers DistributedArray.h
DistributedArray_impl.h
DistributedArrayView.h
DistributedArrayView_impl.h
DistributedMatrix.h
DistributedMatrix_impl.h
DistributedSpMV.h
DistributedVector.h
DistributedVector_impl.h
DistributedVectorView.h
DistributedVectorView_impl.h
Partitioner.h
Subrange.h
)
INSTALL( FILES ${headers} DESTINATION ${TNL_TARGET_INCLUDE_DIRECTORY}/DistributedContainers )
......@@ -521,7 +521,10 @@ String
DistributedMesh< Grid< Dimension, Real, Device, Index > >::
printProcessDistr() const
{
return convertToString(this->domainDecomposition[0])+String("-")+convertToString(this->domainDecomposition[1])+String("-")+convertToString(this->domainDecomposition[2]);
String res = convertToString(this->domainDecomposition[0]);
for(int i=1; i<Dimension; i++)
res=res+String("-")+convertToString(this->domainDecomposition[i]);
return res;
};
template< int Dimension, typename Real, typename Device, typename Index >
......
......@@ -16,6 +16,7 @@ SET( headers BICGStab.h
SOR_impl.h
TFQMR.h
TFQMR_impl.h
Traits.h
UmfpackWrapper.h
UmfpackWrapper_impl.h
)
......
......@@ -15,8 +15,7 @@
#include <TNL/Containers/VectorView.h>
#include <TNL/Pointers/SharedPointer.h>
#include <TNL/Config/ParameterContainer.h>
#include "../Traits.h"
#include <TNL/Solvers/Linear/Traits.h>
namespace TNL {
namespace Solvers {
......
......@@ -31,9 +31,6 @@ target_link_libraries (tnl-dicom-reader tnl ${DCMTK_LIBRARIES} )
ADD_EXECUTABLE(tnl-lattice-init tnl-lattice-init.cpp )
target_link_libraries (tnl-lattice-init tnl )
ADD_EXECUTABLE( tnl-functions-benchmark functions-benchmark.cpp )
target_link_libraries( tnl-functions-benchmark tnl )
IF( BUILD_CUDA )
CUDA_ADD_EXECUTABLE( tnl-cuda-arch tnl-cuda-arch.cu )
INSTALL( TARGETS tnl-cuda-arch
......@@ -55,9 +52,7 @@ INSTALL( TARGETS tnl-init
INSTALL( FILES ${PROJECT_TOOLS_PATH}/tnl-bindir
${PROJECT_TOOLS_PATH}/tnl-compile
${PROJECT_TOOLS_PATH}/tnl-link
tnl-time-series2png
tnl-err2eoc
tnl-eoc-test-log
tnl-log-to-html.py
DESTINATION bin
PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
......
/***************************************************************************
functions-benchmark.cpp - description
-------------------
begin : Jul 4, 2010
copyright : (C) 2010 by Tomas Oberhuber
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
/* See Copyright Notice in tnl/Copyright */
#include "functions-benchmark.h"
int main( int argc, char* argv[] )
{
const long int loops = 1 << 24;
std::cout << "Runnning benchmarks in single precision on CPU ... " << std::endl;
benchmarkAddition< float >( loops );
benchmarkMultiplication< float >( loops );
benchmarkDivision< float >( loops );
benchmarkSqrt< float >( loops );
benchmarkSin< float >( loops );
benchmarkExp< float >( loops );
benchmarkPow< float >( loops );
std::cout << "Runnning benchmarks in double precision on CPU ... " << std::endl;
benchmarkAddition< double >( loops );
benchmarkMultiplication< float >( loops );
benchmarkDivision< double >( loops );
benchmarkSqrt< double >( loops );
benchmarkSin< double >( loops );
benchmarkExp< double >( loops );
benchmarkPow< double >( loops );
return EXIT_SUCCESS;
}
/***************************************************************************
functions-benchmark.h - description
-------------------
begin : Jul 4, 2010
copyright : (C) 2010 by Tomas Oberhuber
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
/* See Copyright Notice in tnl/Copyright */
#ifndef FUNCTIONSBENCHMARK_H_
#define FUNCTIONSBENCHMARK_H_
#include <iostream>
#include <math.h>
#include <TNL/Timer.h>
using namespace TNL;
template< typename REAL > void benchmarkAddition( long int loops )
{
std::cout << "Benchmarking addition on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1 = 1.2;
REAL a2 = 1.2;
REAL a3 = 1.2;
REAL a4 = 1.2;
for( long int i = 0; i < loops; i ++ )
{
a1 += REAL( 0.1 );
a2 += REAL( 0.1 );
a3 += REAL( 0.1 );
a4 += REAL( 0.1 );
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
template< typename REAL > void benchmarkMultiplication( const long int loops )
{
std::cout << "Benchmarking multiplication on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1 = 1.0e9;
REAL a2 = 1.0e9;
REAL a3 = 1.0e9;
REAL a4 = 1.0e9;
for( long int i = 0; i < loops; i ++ )
{
{
a1 *= REAL( 0.99 );
a2 *= REAL( 0.99 );
a3 *= REAL( 0.99 );
a4 *= REAL( 0.99 );
if( a1 < REAL( 0.01 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
}
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 * a2 * a3 * a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
template< typename REAL > void benchmarkDivision( long int loops )
{
std::cout << "Benchmarking division on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1( 1.0e9 );
REAL a2( 1.0e9 );
REAL a3( 1.0e9 );
REAL a4( 1.0e9 );
for( long int i = 0; i < loops; i ++ )
{
a1 /= REAL( 1.1 );
a2 /= REAL( 1.1 );
a3 /= REAL( 1.1 );
a4 /= REAL( 1.1 );
if( a1 < REAL( 0.01 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 / a2 / a3 / a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops / 2 ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
template< typename REAL > void benchmarkSqrt( long int loops )
{
std::cout << "Benchmarking sqrt on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1( 1.0e9 );
REAL a2( 1.0e9 );
REAL a3( 1.0e9 );
REAL a4( 1.0e9 );
for( long int i = 0; i < loops; i ++ )
{
a1 = ::sqrt( a1 );
a2 = ::sqrt( a2 );
a3 = ::sqrt( a3 );
a4 = ::sqrt( a4 );
if( a1 < REAL( 100.0 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops / 2 ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
template< typename REAL > void benchmarkSin( long int loops )
{
std::cout << "Benchmarking sin on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1( 1.0e9 );
REAL a2( 1.0e9 );
REAL a3( 1.0e9 );
REAL a4( 1.0e9 );
for( long int i = 0; i < loops; i ++ )
{
a1 = ::sin( a1 );
a2 = ::sin( a2 );
a3 = ::sin( a3 );
a4 = ::sin( a4 );
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
template< typename REAL > void benchmarkExp( long int loops )
{
std::cout << "Benchmarking exp on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1( 1.1 );
REAL a2( 1.1 );
REAL a3( 1.1 );
REAL a4( 1.1 );
for( long int i = 0; i < loops; i ++ )
{
a1 = exp( a1 );
a2 = exp( a2 );
a3 = exp( a3 );
a4 = exp( a4 );
if( a1 > REAL( 1.0e9 ) ) a1 = a2 = a3 = a4 = REAL( 1.1 );
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
template< typename REAL > void benchmarkPow( long int loops )
{
std::cout << "Benchmarking pow on CPU ( " << loops << " loops ) ... " << std::flush;
Timer timer;
timer.start();
REAL a1( 1.0e9 );
REAL a2( 1.0e9 );
REAL a3( 1.0e9 );
REAL a4( 1.0e9 );
for( long int i = 0; i < loops; i ++ )
{
a1 = ::pow( a1, REAL( 0.9 ) );
a2 = ::pow( a2, REAL( 0.9 ) );
a3 = ::pow( a3, REAL( 0.9 ) );
a4 = ::pow( a4, REAL( 0.9 ) );
if( a1 < REAL( 1.0 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
}
double cpu_time = timer.getCPUTime();
std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
}
#endif /* FUNCTIONSBENCHMARK_H_ */
/***************************************************************************
polygonizer.cpp - description
-------------------
begin : Mon Feb 11 2002
copyright : (C) 2002 by Tomas Oberhuber
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
/* See Copyright Notice in tnl/Copyright */
#include <string.h>
#include "polygonizer.h"
#include "../solids/reciever.h"
#include "../solids/solid.h"
#include "hash.h"
POLYGONIZER polygonizer;
static int corner1[ 12 ] =
{ LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF };
static int corner2[ 12 ] =
{ LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF };
// these are fields of corners for edges in order
// LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF
static int left_face[ 12 ] =
{ B, L, L, F, R, T, N, R, N, B, T, F };
// face on left when going corner1 to corner2
static int right_face[ 12 ] =
{ L, T, N, L, B, R, R, F, B, F, N, T };
// face on right when going coner1 to corner2
//---------------------------------------------------------------------------
CUBE :: CUBE( int _i, int _j, int _k )
: i( _i ), j( _j ), k( _k )
{
bzero( corners, 8 * sizeof( CORNER* ) );
}
//---------------------------------------------------------------------------
POLYGONIZER :: POLYGONIZER()
{
edge_hash = new ( list< EDGE* >* )[ 2 * hash_size ];
bzero( edge_hash, 2 * hash_size * sizeof( list< EDGE* >* ) );
corner_hash = new ( list< CORNER* >* )[ 2 * hash_size ];
bzero( corner_hash, 2 * hash_size * sizeof( list< CORNER* >* ) );
center_hash = new ( list< CENTER* >* )[ 2 * hash_size ];
bzero( center_hash, 2 * hash_size * sizeof( list< CENTER* >* ) );
Make_Cube_Table();
}
//---------------------------------------------------------------------------
POLYGONIZER :: ~POLYGONIZER()
{
assert( edge_hash );
delete[] edge_hash;
assert( corner_hash );
delete[] corner_hash;
assert( center_hash );
delete[] center_hash;
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Parametric( const SOLID* sld, RECIEVER* reciever,
const unsigned int slices, const unsigned int stacks)
{
assert( slices && stacks );
const double x_step = 1.0 / slices;
const double y_step = 1.0 / stacks;
for( unsigned int i = 0; i < stacks; i ++ )
for( unsigned int j = 0; j < slices; j ++ )
{
double x = j * x_step;
double y = i * y_step;
VECTOR n1, n2, n3, n4;
VECTOR u1 = sld -> Surface_Point( x, y, &n1 );
VECTOR u2 = sld -> Surface_Point( x + x_step, y, &n2 );
VECTOR u3 = sld -> Surface_Point( x + x_step, y + y_step, &n3 );
VECTOR u4 = sld -> Surface_Point( x, y + y_step, &n4 );
VERTEX v1( u1, n1 );
VERTEX v2( u2, n2 );
VERTEX v3( u3, n3 );
VERTEX v4( u4, n4 );
/*TRIANGLE* t1 = new TRIANGLE( new VERTEX( v1 ),
new VERTEX( v4 ),
new VERTEX( v2 ) ) ;
TRIANGLE* t2 = new TRIANGLE( new VERTEX( v3 ),
new VERTEX( v2 ),
new VERTEX( v4 ) ) ;*/
if( ! reciever -> Insert_Triangle( v1, v4, v2 ) ) return 0;
if( ! reciever -> Insert_Triangle( v3, v2, v4 ) ) return 0;
}
return 1;
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Init( const VECTOR& pos, const SOLID* sld, const double& cb_sz )
{
cube_size = cb_sz;
solid = sld;
VECTOR in( pos ), out( pos );
if( ! Find_Point( 1, in, 1.0, solid ) ||
! Find_Point( 0, out, 1.0, solid ) ) return 0;
CUBE *c = NULL;
if( cube_stack. empty() ) // this is the initial cube
{
position = Compute_Surface_Point( in, out, solid -> Continuous_Function( in ) );
c = new CUBE( 0, 0, 0 );
}
else
{
VECTOR _pos = Compute_Surface_Point( in, out, solid -> Continuous_Function( in ) ) - position;
_pos *= 1.0 / cube_size;
int i( ( int ) _pos. x ),
j( ( int ) _pos. y ),
k( ( int ) _pos. z );
if( Set_Center( i, j, k ) ) return 1; //we have already found this cube
c = new CUBE( i, j, k );
}
for( int i = 0; i < 8; i ++ ) c -> corners[ i ] =
Set_Corner( ( i >> 2 ) & 1, ( i >> 1 ) & 1, i & 1 );
cube_stack. push( c );
return 1;
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Implicit( const SOLID* sld, RECIEVER* _reciever,
const VECTOR& ps, const VECTOR& cr,
const double& cb_size, int mode )
{
solid = sld;
reciever = _reciever;
cube_size = cb_size;
bound_ps = ps;
bound_cr = cr;
if( cube_stack. empty() && ! Init( solid -> Position(), solid, cb_size ) )
{
std::cerr << "Can not find initial points for polygonization" << std::endl;
return 0;
}
std::cout << "Starting polygonizer ... " << std::endl;
/*VECTOR in( sld -> Position() ), out( sld -> Position() );
if( ! Find_Point( 1, in, 1.0 ) ||
! Find_Point( 0, out, 1.0 ) )
{
std::cerr << "Can not find initial points for polygonization" << std::endl;
return 0;
}
std::cout << "Starting polygonizer ... " << std::endl;
position = Compute_Surface_Point( in, out, solid -> Continuous_Function( in ) );
CUBE *c = new CUBE( 0, 0, 0 );
for( int i = 0; i < 8; i ++ ) c -> corners[ i ] =
Set_Corner( ( i >> 2 ) & 1, ( i >> 1 ) & 1, i & 1 );
cube_stack. push( c );*/
int jkl;
CUBE* c;
while( ! cube_stack. empty() )
{
c = cube_stack. top();
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
cube_stack. pop();
jkl = cube_stack. size();
if( ! Triangulate_Cube( c ) )
{
Free_Memory();
return 0;
}
Test_Face( c -> i - 1, c -> j, c -> k, c, L, LBN, LBF, LTN, LTF );
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
Test_Face( c -> i + 1, c -> j, c -> k, c, R, RBN, RBF, RTN, RTF );
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
Test_Face( c -> i, c -> j - 1, c -> k, c, B, LBN, LBF, RBN, RBF );
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
Test_Face( c -> i, c -> j + 1, c -> k, c, T, LTN, LTF, RTN, RTF );
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
Test_Face( c -> i, c -> j, c -> k - 1, c, N, LBN, LTN, RBN, RTN );
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
Test_Face( c -> i, c -> j, c -> k + 1, c, F, LBF, LTF, RBF, RTF );
std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
}
std::cout << std::endl;
Free_Memory();
return 1;
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Triangulate_Cube( CUBE* cube )
{
#ifdef DBG_POLYGONIZER
std::cout << "Cube polygonize: " << cube -> i << " " << cube -> j << " " << cube -> k << std::endl;;
#endif
int index = 0;
// this is a index of cube in cube_table, it will be count by
// the signs of function in corners
for( int i = 0; i < 8; i ++ )
if( cube -> corners[ i ] -> value > 0.0 ) index += ( 1 << i );
for( size_t i = 0; i < cube_table[ index ]. size(); i ++ )
{
vector< int >& edges_vector = ( cube_table[ index ] )[ i ];
VERTEX *v1( 0 ), *v2( 0 ), *v3( 0 );
for( size_t j = 0; j < edges_vector. size(); j ++ )
{
CORNER* c1 = cube -> corners[ corner1[ edges_vector[ j ] ] ];
CORNER* c2 = cube -> corners[ corner2[ edges_vector[ j ] ] ];
VERTEX* c = Get_Vertex( c1, c2 );
if( j == 0 ) v1 = c;
if( j == 1 ) v2 = c;
if( j > 1 )
{
v3 = c;
if( ! reciever -> Insert_Triangle( *v1, *v3, *v2 ) ) return 0;
v2 = v3;
}
}
}
return 1;
}
//---------------------------------------------------------------------------
VERTEX* POLYGONIZER :: Get_Vertex( CORNER* c1, CORNER* c2 )
{
VERTEX* vertex;
if( Check_Edge( c1, c2, vertex ) ) return vertex;
// the vertex has been already computed
VECTOR p = Compute_Surface_Point( c1 -> position, c2 -> position, c1 -> value );
vertex = new VERTEX( p, solid -> Normal( p ) );
Set_Edge( c1, c2, vertex );
return vertex;
}
//---------------------------------------------------------------------------
VECTOR POLYGONIZER :: Compute_Surface_Point( const VECTOR& v1, const VECTOR& v2, const double& val )
{
VECTOR pos, neg, p;
if( val < 0.0 )
{
pos = v2; neg = v1;
}
else
{
pos = v1; neg = v2;
}
int i = 0;
while( i ++ < max_iter )
{
p = 0.5 * ( pos + neg );
if( solid -> Continuous_Function( p ) > 0 ) pos = p;
else neg = p;
}
return p;
}
//---------------------------------------------------------------------------
void POLYGONIZER :: Set_Edge( CORNER* c1, CORNER* c2, VERTEX* v )
{
int& i1 = c1 -> i;
int& j1 = c1 -> j;
int& k1 = c1 -> k;
int& i2 = c2 -> i;
int& j2 = c2 -> j;
int& k2 = c2 -> k;
if( i1 > i2 || ( i1 == i2 && ( j1 > j2 || ( j1 == j2 && k1 > k2 ) ) ) )
{
CORNER* c = c1; c1 = c2; c2 = c;
}
unsigned int index = hash_index( i1, j1, k1 ) + hash_index( i2, j2, k2 );
if( ! edge_hash[ index ] ) edge_hash[ index ] = new list< EDGE* >;
edge_hash[ index ] -> push_back( new EDGE( c1, c2, v ) );
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Check_Edge( CORNER* c1, CORNER* c2, VERTEX*& vertex )
{
int& i1 = c1 -> i;
int& j1 = c1 -> j;
int& k1 = c1 -> k;
int& i2 = c2 -> i;
int& j2 = c2 -> j;
int& k2 = c2 -> k;
if( i1 > i2 || ( i1 == i2 && ( j1 > j2 || ( j1 == j2 && k1 > k2 ) ) ) )
{
CORNER* c = c1; c1 = c2; c2 = c;
}
list< EDGE* >* el = edge_hash[ hash_index( i1, j1, k1 ) + hash_index( i2, j2, k2 ) ];
if( ! el ) return 0;
list< EDGE* > :: iterator it = el -> begin();
while( it != el -> end() )
{
EDGE* e = * it ++;
if( e -> corners[ 0 ] == c1 && e -> corners[ 1 ] == c2 )
{
vertex = e -> vertex;
return 1;
}
}
return 0;
}
//---------------------------------------------------------------------------
CORNER* POLYGONIZER :: Set_Corner( int i, int j, int k )
{
#ifdef DBG_POLYGONIZER
std::cout << "Set Corner index: " << i << " " << j << " " << k << std::endl;
#endif
int index = hash_index( i, j, k );
list< CORNER* >*& cl = corner_hash[ index ];
if( ! cl ) cl = new list< CORNER* >;
list< CORNER* > :: iterator it = cl -> begin();
CORNER* cr;
while( it != cl -> end() )
{
cr = * it ++;
if( ( cr -> i == i ) && ( cr -> j == j ) && ( cr -> k == k ) ) return cr;
}
VECTOR tmp = ( VECTOR( i, j, k ) - VECTOR( 0.5 ) );
VECTOR p = position + cube_size * tmp;
double val = solid -> Continuous_Function( p );
#ifdef DBG_POLYGONIZER
std::cout << "Set Corer: value " << val << std::endl;
#endif
cr = new CORNER( i, j, k, p, val );
cl -> push_back( cr );
return cr;
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Set_Center( int i, int j, int k )
{
#ifdef DBG_POLYGONIZER
std::cout << "Set Center: " << i << " " << j << " " << k << std::endl;
#endif
list< CENTER* >*& cl = center_hash[ hash_index( i, j, k ) ];
if( ! cl ) cl = new list< CENTER* >;
list< CENTER* > :: iterator it = cl -> begin();
while( it != cl -> end() )
{
CENTER* cr = * it ++;
if( cr -> i == i && cr -> j == j && cr -> k == k ) return 1;
}
cl -> push_back( new CENTER( i, j, k ) );
return 0;
}
//---------------------------------------------------------------------------
void POLYGONIZER :: Make_Cube_Table()
{
int done_edges[ 12 ], positive[ 8 ];
// we have 256 different possibilities to sign all 8 corners of a cube
// they are describe here in cube_table
for( int i = 0; i < 256; i ++ )
{
for( int e = 0; e < 12; e ++ ) done_edges[ e ] = 0;
for( int c = 0; c < 8; c ++ ) positive[ c ] = ( i >> c ) & 1;
// each bit in 'i' represent one corner of cube and sign of function value
// '1' means positive value, '0' means negative value
for( int e = 0; e < 12; e ++ )
{
if( ! done_edges[ e ] &&
// we didn't process this edge yet ...
( positive[ corner1[ e ] ] != positive[ corner2[ e ] ] ) )
// ... and corners of this edge have different sign
{
int start = e;
int edge = e;
int face = positive[ corner1[ e ] ] ? right_face[ e ] : left_face[ e ];
// get face that is to right of edge from positive to negatve corner
vector< int > temp_vector;
while( 1 )
{
edge = Next_Clockwise_Edge( edge, face );
done_edges[ edge ] = 1;
if( positive[ corner1[ edge ] ] != positive[ corner2[ edge ] ] )
{
temp_vector. push_back( edge );
if( edge == start ) break;
face = Other_Face( edge, face );
}
}
cube_table[ i ]. push_back( temp_vector );
}
}
}
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Next_Clockwise_Edge( int edge, int face )
{
switch( edge )
{
case LB: return ( face == L ) ? LF : BN;
case LT: return ( face == L ) ? LN : TF;
case LN: return ( face == L ) ? LB : TN;
case LF: return ( face == L ) ? LT : BF;
case RB: return ( face == R ) ? RN : BF;
case RT: return ( face == R ) ? RF : TN;
case RN: return ( face == R ) ? RT : BN;
case RF: return ( face == R ) ? RB : TF;
case BN: return ( face == B ) ? RB : LN;
case BF: return ( face == B ) ? LB : RF;
case TN: return ( face == T ) ? LT : RN;
case TF: return ( face == T ) ? RT : LF;
}
return 0; // this is just for avoiding compiler warning
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Other_Face( int edge, int face )
{
int other_face = left_face[ edge ];
return face == other_face ? right_face[ edge ] : other_face;
}
//---------------------------------------------------------------------------
void POLYGONIZER :: Test_Face( int i, int j, int k, CUBE* old, int face,
int c1, int c2, int c3, int c4 )
{
static int face_bit[ 6 ] = { 2, 2, 1, 1, 0, 0 };
int bit = face_bit[ face ];
int pos = old -> corners[ c1 ] -> value > 0.0 ? 1 : 0;
// test id no surface crossing, cube out of bounds, or already visited
if( ( old -> corners[ c2 ] -> value > 0.0 ) == pos &&
( old -> corners[ c3 ] -> value > 0.0 ) == pos &&
( old -> corners[ c4 ] -> value > 0.0 ) == pos ) return;
// test bounds
if( bound_ps < bound_cr )
{
VECTOR p( position + cube_size * VECTOR( ( double ) i, ( double ) j, ( double ) k ) );
if( ! ( p >= bound_ps && p <= bound_cr ) ) return;
}
else
if( abs( i ) > 50 || abs( j ) > 50 || abs( k ) > 50 ) return;
if( Set_Center( i, j, k ) ) return;
CUBE* new_cube = new CUBE( i, j, k );
// given face of old cube is the same as the opposite face of new_cube
new_cube -> corners[ flip_bit( c1, bit ) ] = old -> corners[ c1 ];
new_cube -> corners[ flip_bit( c2, bit ) ] = old -> corners[ c2 ];
new_cube -> corners[ flip_bit( c3, bit ) ] = old -> corners[ c3 ];
new_cube -> corners[ flip_bit( c4, bit ) ] = old -> corners[ c4 ];
#ifdef DBG_POLYGONIZER
std::cout << "Test Face indexes: " << i << " " << j << " " << k << std::endl;
#endif
for( int n = 0; n < 8; n ++ )
if( ! new_cube -> corners[ n ] )
new_cube -> corners[ n ] = Set_Corner( i + ( ( n >> 2 ) & 1 ),
j + ( ( n >> 1 ) & 1 ),
k + ( n & 1 ) );
cube_stack. push( new_cube );
}
//---------------------------------------------------------------------------
int POLYGONIZER :: Find_Point( int sign, VECTOR& point, double size, const SOLID* sld )
{
VECTOR init = point;
VECTOR ps( -0.5 );
VECTOR cr( 0.5 );
for( int i = 0; i < 10000; i ++ )
{
point = init + size * Random_Vector( ps, cr );
//cout << point << " - " << sld -> Continuous_Function( point ) << std::endl;
if( sign == ( sld -> Sign_Function( point ) == 1 ) )
return 1;
size *= 1.005;
}
return 0;
}
//---------------------------------------------------------------------------
void POLYGONIZER :: Free_Memory()
{
// clear hashes
for( size_t i = 0; i < 2 * hash_size; i ++ )
{
if( edge_hash[ i ] )
{
list< EDGE* > :: iterator it = edge_hash[ i ] -> begin();
while( it != edge_hash[ i ] -> end() ) delete * it ++;
delete edge_hash[ i ];
edge_hash[ i ] = NULL;
}
if( corner_hash[ i ] )
{
list< CORNER* > :: iterator it = corner_hash[ i ] -> begin();
while( it != corner_hash[ i ] -> end() ) delete * it ++;
delete corner_hash[ i ];
corner_hash[ i ] = NULL;
}
if( center_hash[ i ] )
{
list< CENTER* > :: iterator it = center_hash[ i ] -> begin();
while( it != center_hash[ i ] -> end() ) delete * it ++;
delete center_hash[ i ];
center_hash[ i ] = NULL;
}
}
}
/***************************************************************************
polygonizer.h - description
-------------------
begin : Mon Feb 11 2002
copyright : (C) 2002 by Tomas Oberhuber
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
/* See Copyright Notice in tnl/Copyright */
#ifndef POLYGONIZER_H
#define POLYGONIZER_H
/**
*@author Tomas Oberhuber
*
* this class is rewritten polygonizer by Jules Bloomenthal, Xerox PARC.
* Copyright of original code (c) Xerox Corporation, 1991.
* See 'implicit.eecs.wsu.edu'.
*
*/
#include <list.h>
#include <stack.h>
#include <vector.h>
#include "vctr.h"
#include "vertex.h"
//using namespace :: std;
class SOLID;
class RECIEVER;
enum { poly_parametric, poly_implicit_cube, poly_implicit_tetrahedron };
const int max_iter = 10;
inline int flip_bit( int i, int bit )
{ return i ^ 1 << bit; };
enum { L = 0, // left direction - x
R, // right direction + x
B, // bottom direction - y
T, // top direction + y
N, // near direction - z
F // far direction + z
};
enum { LB = 0, // left bottom edge
LT, // left top edge
LN, // left near edge
LF, // left far edge
RB, // right bottom edge
RT, // right top edge
RN, // right near edge
RF, // right far edge
BN, // bottom near edge
BF, // bottom far edge
TN, // top near edge
TF // top far edge
};
enum { LBN = 0, // left bottom near corner
LBF, // left bottom far corner
LTN, // left top near corner
LTF, // left top far corner
RBN, // right bottom near corner
RBF, // right bottom far corner
RTN, // right top near corner
RTF // right top far corner
};
struct CORNER
{
int i, j, k;
VECTOR position;
double value;
CORNER( int _i, int _j, int _k, const VECTOR& p, const double& val )
:i( _i ), j( _j ), k( _k ), position( p ), value( val ){};
};
struct EDGE
{
CORNER* corners[ 2 ];
VERTEX* vertex;
EDGE( CORNER* c1, CORNER* c2, VERTEX* v )
: vertex( v ){ corners[ 0 ] = c1; corners[ 1 ] = c2;};
~EDGE() { assert( vertex ); delete vertex; };
};
struct CUBE
{
int i, j, k;
// cube lattice
CORNER* corners[ 8 ];
CUBE( int, int, int );
};
struct CENTER
{
int i, j, k;
CENTER( int _i, int _j, int _k )
: i( _i ), j( _j ), k( _k ){};
};
class POLYGONIZER
{
public:
POLYGONIZER();
~POLYGONIZER();
//! Creates initial cubes for implicite polygonizer
/*! We need to call this method several times in case of
non-connected surfaces. It si at least once called automaticly
if it was not called before starting the process of tesselation.
*/
int Init( const VECTOR& pos, const SOLID* sld, const double& cb_sz );
//! Method for implicite surfaces
/*! start polygonizer with given SOLID, triangles are inserted into TRIANGLES
two const VECTOR&s define the bound for polygonizations
const double& is size of elementar cube
int is poligonization mode - polygonize cube or tetrahedrons
*/
int Implicit( const SOLID*, RECIEVER*, const VECTOR&, const VECTOR&, const double&, int );
int Parametric( const SOLID*, RECIEVER*, const unsigned int, const unsigned int );
protected:
//! This solid will be polygonized using method SOLID :: Solid_Function( const VECTOR& );
const SOLID* solid;
//! Reciver of emitted triangles
RECIEVER* reciever;
//! Position of the first cube
/*! Thus, it is also position of the origin for cubes indexing
*/
VECTOR position;
//! Bounds for polygonizer
/*! All coubes outside culled.
*/
VECTOR bound_ps, bound_cr;
//! Size of the cube and
double cube_size;
stack< CUBE* > cube_stack;
vector< vector< int > > cube_table[ 256 ];
list< EDGE* >** edge_hash;
// hash field of edges lists
list< CORNER* >** corner_hash;
// another hash for corners
list< CENTER* >** center_hash;
// yet another hash for cube centers
int Triangulate_Cube( CUBE* );
// triangulate the cube directly, without decomposition
VERTEX* Get_Vertex( CORNER*, CORNER* );
// return vertex for given edge using edge hash table ( Get_Edge method )
// both corners values are presumed of different sign
VECTOR Compute_Surface_Point( const VECTOR&, const VECTOR&, const double& );
// compute vertex on given edge using solid function of given solid
void Set_Edge( CORNER*, CORNER*, VERTEX* );
// insert edge to hash table
int Check_Edge( CORNER*, CORNER*, VERTEX*& );
// get vertex on edge from hash table
// return 1 if edge is in hash table
// 0 if edge is not in hash table
CORNER* Set_Corner( int, int, int );
// return corner with given lattice location
// set and cache its function value
int Set_Center( int, int, int );
// set ( i, j, k ) entry to center_hash table
// return 1 if already set; otherwise set and return 0
void Make_Cube_Table();
// create cube_table
int Next_Clockwise_Edge( int, int );
// return next clockwise edge from given edge around given face
int Other_Face( int, int );
// return face adjoining edge is not the given face
void Test_Face( int, int, int, CUBE*, int, int, int, int, int );
// given cube at lattice ( i, j, k ) and four corners of face,
// if surface crosses face, compute other four corners of adjacent cube
// and add new cube to cube stack
int Find_Point( int, VECTOR&, double, const SOLID* sld );
// find point with given value sign
void Free_Memory();
};
extern POLYGONIZER polygonizer;
#endif
......@@ -12,5 +12,5 @@ do
esac
done
echo -I@CMAKE_INSTALL_PREFIX@/@TNL_TARGET_INCLUDE_DIRECTORY@ ${CUDA_FLAGS} ${CXX_STD_FLAGS} ${DEBUG_FLAGS}
echo -I@CMAKE_INSTALL_PREFIX@/include ${CUDA_FLAGS} ${CXX_STD_FLAGS} ${DEBUG_FLAGS}
#!/usr/bin/env python
import sys, string, math
......@@ -20,31 +20,21 @@ template<> struct ConfigTagIndex< {problemBaseName}BuildConfigTag, short int >{{
template<> struct ConfigTagIndex< {problemBaseName}BuildConfigTag, long int >{{ enum {{ enabled = false }}; }};
/****
* With how many dimensions may have the problem to be solved...
*/
template< int Dimension > struct ConfigTagDimension< {problemBaseName}BuildConfigTag, Dimension >{{ enum {{ enabled = ( Dimension == 1 ) }}; }};
/****
* Use of Meshes::Grid is enabled for allowed dimensions and Real, Device and Index types.
* The mesh type will be resolved by the Solver by default.
* (The detailed mesh configuration is in TNL/Meshes/BuildConfigTags.h)
*/
template< int Dimension, typename Real, typename Device, typename Index >
struct ConfigTagMesh< {problemBaseName}BuildConfigTag, Meshes::Grid< Dimension, Real, Device, Index > >
{{ enum {{ enabled = ConfigTagDimension< {problemBaseName}BuildConfigTag, Dimension >::enabled &&
ConfigTagReal< {problemBaseName}BuildConfigTag, Real >::enabled &&
ConfigTagDevice< {problemBaseName}BuildConfigTag, Device >::enabled &&
ConfigTagIndex< {problemBaseName}BuildConfigTag, Index >::enabled }}; }};
template<> struct ConfigTagMeshResolve< {problemBaseName}BuildConfigTag >{{ enum {{ enabled = true }}; }};
/****
* Please, chose your preferred time discretization here.
* All time discretisations (explicit, semi-impicit and implicit ) are
* enabled by default.
*/
template<> struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, ExplicitTimeDiscretisationTag >{{ enum {{ enabled = true }}; }};
template<> struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, SemiImplicitTimeDiscretisationTag >{{ enum {{ enabled = true }}; }};
template<> struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, ImplicitTimeDiscretisationTag >{{ enum {{ enabled = false }}; }};
template< typename TimeDiscretisation > struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, TimeDiscretisation >{{ enum {{ enabled = true }}; }};
/****
* Only the Runge-Kutta-Merson solver is enabled by default.
* All explicit solvers are enabled by default
*/
template<> struct ConfigTagExplicitSolver< {problemBaseName}BuildConfigTag, ExplicitEulerSolverTag >{{ enum {{ enabled = false }}; }};
template< typename ExplicitSolver > struct ConfigTagExplicitSolver< {problemBaseName}BuildConfigTag, ExplicitSolver >{{ enum {{ enabled = true }}; }};
}} // namespace Solvers
}} // namespace TNL
\ No newline at end of file
}} // namespace TNL
......@@ -47,7 +47,8 @@ template< typename Real,
typename Index,
typename MeshType,
typename ConfigTag,
typename SolverStarter >
typename SolverStarter,
typename Communicator >
class {problemBaseName}Setter
{{
public:
......@@ -75,12 +76,12 @@ class {problemBaseName}Setter
if( boundaryConditionsType == "dirichlet" )
{{
typedef Operators::DirichletBoundaryConditions< MeshType, ConstantFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}}
typedef Operators::NeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}}
......@@ -88,12 +89,12 @@ class {problemBaseName}Setter
if( boundaryConditionsType == "dirichlet" )
{{
typedef Operators::DirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}}
typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}}
......
......@@ -34,7 +34,7 @@ operator()( const MeshFunction& u,
* The following example is the Laplace operator approximated
* by the Finite difference method.
*/
static_assert( MeshEntity::entityDimension == {meshDimension}, "Wrong mesh entity dimensions." );
static_assert( MeshEntity::getEntityDimension() == {meshDimension}, "Wrong mesh entity dimensions." );
static_assert( MeshFunction::getEntitiesDimension() == {meshDimension}, "Wrong preimage function" );
const typename MeshEntity::template NeighborEntities< {meshDimension} >& neighborEntities = entity.getNeighborEntities();
......@@ -84,7 +84,7 @@ setMatrixElements( const PreimageFunction& u,
Matrix& matrix,
Vector& b ) const
{{
static_assert( MeshEntity::entityDimension == {meshDimension}, "Wrong mesh entity dimensions." );
static_assert( MeshEntity::getEntityDimension() == {meshDimension}, "Wrong mesh entity dimensions." );
static_assert( PreimageFunction::getEntitiesDimension() == {meshDimension}, "Wrong preimage function" );
/****
......