Commit b73dc746 authored by Kirill Terekhov's avatar Kirill Terekhov

Huge update

Huge update involving update of underlying data representation and
elements management. Different interaction with iterators. New
ElementSet representation.  Some annotation for documentation is added.
New data structures. Examples require update. (check
examples/CMakeLists.txt)
parent b214f48f
......@@ -20,7 +20,9 @@ set(SOURCE solver.cpp
element.cpp
mesh_parallel.cpp
modify.cpp
autodiff.cpp)
earray.cpp
comparator.cpp
autodiff.cpp)
set(HEADER inmost.h
inmost_options_cmake.h
inmost_common.h
......
mesh + !partitioner + !domain + solvers + +/-nonlinear solvers + autodiff + !system handler
:
0) () , data()
1) () erase array , resize'
2) () ReferenceArray
3) () char, short, - , -
4) ()
5) ( )
6) ()
7) () chunk_bulk_array - chunk_array - ,
8) ()
8.1)
10) ()
10.-1) () Element getNodes,getEdges,getFaces,getCells
10.0) (, ) array,
10.1) () adjacent ( ID , )
10.1.0) (, , 10.6) - array ( reference_array) dynarray ( )
10.1.1) ()
10.2) () ElementSet - , LowConn , ,
HighConn , 1- - , - , - ( INMOST_experimental)
tag_set_size, , LowConn
?
,
10.2.0) () ,
:
- , , , , insert
- array, , swap array
,
10.2.1) () - , , , 10.1
- ,
10.2.2) () ApplyModification
10.2.3) () VTK
10.2) () ID
10.3) () Cell, Face, Edge, Node, ElementSet - , ID
10.4) () Cell,Face,Edge,Node,ElementSet ID
10.5) () reference_array
10.6) () createCell, CreateFace.. - TieElement new, 10.3
11) () Markers, HighConn, LowConn - TagManager, , ,
26) ()
26.1) ( ) radix ( INMOST_experimental)
25*n 8- , 16
17*n 11- , 131
27) () Read* pmf,
28) () std::vector array
29) () std::vector array
31) () TagManager & TagManager::assign(Mesh * m, TagManager const & other)
33) () ElementNum, ElementType 6
34) () Storage ,
36) ( ) , Storage Mesh
50) () container.hpp enumerator size_type
51) () container.hpp ( iterator_category?)
41) ( ) HighConn, LowConn ,
44) ( ) PackElementsData tag
46) ( ) . ,
14) ReferenceArray?
16) parallel_storage
17) ?
18)
19) ( , )
20)
21) /- (kd examples/OldDrawGrid)
22) ResolveShared
23) ,
24)
25) ,
30) - sort ReorderApply , chunk_array
ReorderApply radix-?
32) , , Cell, Face, Edge Node
35) Bridge*
37) container.hpp,
38) GatherBoundaryFaces/GatherInteriorFaces
39) CastRay ,
40) ?
42) Inside discr_common.cpp
43) RestoreCellNodes?
45) ?
. ( )
47) , , ?
49) -
,
ResolveShared
ModificationResolve
52) , radix big-endian
53) SortByGlobalID, globalid
54) ( unit tests)
7.MSPP: : , ,
8.MSPP:
9.MSPP: visual studio
10.MSPP: ( )
11.
12. (, )
12. MSPP: Mesh
13. MSPP: Redistribute!!! ExchangeGhost
13. INMOST / MMTK = Mathematical Modelling ToolKit / NSTK
last_created_element -
1) - ,
( ),
-
,
-
2) small_hash -> dynarray
3)
4)
5) parent()
6)
7) , []
8)
9) ,
( )
10) ,
11) LLVM asmjit
12) openmp
13) opencl , openmp
14)
1) -
1)
2) ILUC2
3) /,
4) E,F
5) EU^-1 L^-1F
6)
7) (openmp)
8) (openmp) EU^-1, L^-1F
9) nicole spillane
10) ,
11) amls
12)
solver:
read mc64
read nested dissection
try to postpone ReorderEF to the end of run
5.MSPP:
6.MSPP:
6.0 bicgs(l)
6.0 OrderInfo -> ASM_preconditioner?
6.0.1 OrderInfo.PrepareMatrix(A,levels) -> ASM_preconditioner(A,levels) : Method
6.0.3 OrderInfo.PrepareVector(vector) -> Matrix.PrepareVector(vector)
6.0.4 OrderInfo.Update(vector) -> Matrix.Update(vector)
6.0.5 orderInfo.Accumulate(vector) -> Matrix.Accumulate(vector)
6.1 (, MPI_Waitall->MPI_Waitsome)
6.2 Update(vector) -> UpdateBegin, UpdateEnd, , (INIT,PREC,MATVEC)
6.3 (?) ilu2 ( superlu?)
6.3.1
info->PrepareMatrix(A,0)
6.3.2
6.4 tau
6.5
6.6 superlu
6.7 - ..
//FOR 2D
// Node -> Vertex
// Edge -> Vertex
......
......@@ -35,7 +35,7 @@ namespace INMOST
hi = (char *) base + width * (num - 1);
recurse:
size = (hi - lo) / width + 1;
size = static_cast<unsigned>((hi - lo) / width) + 1;
if (size <= CUTOFF)
{
......
This diff is collapsed.
This diff is collapsed.
#include "inmost_mesh.h"
#if defined(USE_MESH)
namespace INMOST
{
int Mesh::CentroidComparator::Compare(const real * a, const real * b)
{
real e = m->GetEpsilon();
for(integer i = 0; i < m->GetDimensions(); i++)
if( fabs(a[i]-b[i]) > e )
{
if( a[i] < b[i] )
return -1;
else
return 1;
}
return 0;
}
bool Mesh::CentroidComparator::operator () (HandleType a, HandleType b)
{
if( a == InvalidHandle() || b == InvalidHandle() ) return a > b;
real ca[3] = {0,0,0}, cb[3] = {0,0,0};
m->GetGeometricData(a,CENTROID,ca);
m->GetGeometricData(b,CENTROID,ca);
return Compare(ca,cb) <= 0;
}
bool Mesh::CentroidComparator::operator () (HandleType a, const real * cb)
{
if( a == InvalidHandle() ) return true;
real ca[3] = {0,0,0};
m->GetGeometricData(a,CENTROID,ca);
return Compare(ca,cb) <= 0;
}
int Mesh::IerarhyComparator::CompareNodes(HandleType a, HandleType b)
{
real_array ca = m->RealArrayDF(a,m->CoordsTag());
real_array cb = m->RealArrayDF(b,m->CoordsTag());
real e = m->GetEpsilon();
for(integer i = 0; i < m->GetDimensions(); ++i)
if( fabs(ca[i]-cb[i]) > e )
{
if( ca[i] > cb[i] )
return 1;
else
return -1;
}
return 0;
}
int Mesh::IerarhyComparator::CompareElements(HandleType a, HandleType b)
{
integer ia = GetHandleElementNum(a);
integer ib = GetHandleElementNum(b);
if( ia == ib )
{
if( ia == 0 ) return CompareNodes(a,b);
Element::adj_type const & adja = m->LowConn(a);
Element::adj_type const & adjb = m->LowConn(b);
if( adja.size() == adjb.size() )
{
for(Element::adj_type::size_type k = 0; k < adja.size(); ++k)
{
int ret = CompareElements(adja[k],adjb[k]);
if( ret != 0 ) return ret;
}
return 0;
}
else return static_cast<int>(adja.size()) - static_cast<int>(adjb.size());
}
else return ia - ib;
}
void Mesh::SortByGlobalID(HandleType * h, enumerator n)
{
if( sizeof(integer) == 4 && n > 64 )
{
Element::adj_type::size_type i, num = static_cast<Element::adj_type::size_type>(n);
array<HandleType> tmp(num);
HandleType * ptra = h;
HandleType * ptrb = tmp.data();
const unsigned int kHist = 256;
unsigned int b0[kHist * 4];
unsigned int *b1 = b0 + kHist;
unsigned int *b2 = b0 + kHist*2;
unsigned int *b3 = b0 + kHist*3;
memset(b0,0,sizeof(unsigned int)*kHist*4);
for (i = 0; i < num; i++)
{
unsigned char * rec = reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]));
++b0[rec[0]];
++b1[rec[1]];
++b2[rec[2]];
++b3[rec[3]];
}
unsigned int sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0, tsum;
for (i = 0; i < kHist; i++)
{
tsum = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = tsum;
tsum = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = tsum;
tsum = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = tsum;
tsum = b3[i] + sum3; b3[i] = sum3 - 1; sum3 = tsum;
}
for (i = 0; i < num; i++) ptrb[++b0[reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]))[0]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b1[reinterpret_cast<unsigned char *>(&GlobalID(ptrb[i]))[1]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b2[reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]))[2]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b3[reinterpret_cast<unsigned char *>(&GlobalID(ptrb[i]))[3]]] = ptrb[i];
}
else if( sizeof(integer) == 8 && n > 128 )
{
Element::adj_type::size_type i, num = static_cast<Element::adj_type::size_type>(n);
array<HandleType> tmp(num);
HandleType * ptra = h;
HandleType * ptrb = tmp.data();
const unsigned int kHist = 256;
unsigned int b0[kHist * 8];
unsigned int *b1 = b0 + kHist;
unsigned int *b2 = b0 + kHist*2;
unsigned int *b3 = b0 + kHist*3;
unsigned int *b4 = b0 + kHist*4;
unsigned int *b5 = b0 + kHist*5;
unsigned int *b6 = b0 + kHist*6;
unsigned int *b7 = b0 + kHist*7;
memset(b0,0,sizeof(unsigned int)*kHist*8);
for (i = 0; i < num; i++)
{
unsigned char * rec = reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]));
++b0[rec[0]];
++b1[rec[1]];
++b2[rec[2]];
++b3[rec[3]];
++b4[rec[4]];
++b5[rec[5]];
++b6[rec[6]];
++b7[rec[7]];
}
unsigned int sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0, tsum;
unsigned int sum4 = 0, sum5 = 0, sum6 = 0, sum7 = 0;
for (i = 0; i < kHist; i++)
{
tsum = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = tsum;
tsum = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = tsum;
tsum = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = tsum;
tsum = b3[i] + sum3; b3[i] = sum3 - 1; sum3 = tsum;
tsum = b4[i] + sum4; b4[i] = sum4 - 1; sum4 = tsum;
tsum = b5[i] + sum5; b5[i] = sum5 - 1; sum5 = tsum;
tsum = b6[i] + sum6; b6[i] = sum6 - 1; sum6 = tsum;
tsum = b7[i] + sum7; b7[i] = sum7 - 1; sum7 = tsum;
}
for (i = 0; i < num; i++) ptrb[++b0[reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]))[0]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b1[reinterpret_cast<unsigned char *>(&GlobalID(ptrb[i]))[1]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b2[reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]))[2]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b3[reinterpret_cast<unsigned char *>(&GlobalID(ptrb[i]))[3]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b4[reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]))[4]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b5[reinterpret_cast<unsigned char *>(&GlobalID(ptrb[i]))[5]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b6[reinterpret_cast<unsigned char *>(&GlobalID(ptra[i]))[6]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b7[reinterpret_cast<unsigned char *>(&GlobalID(ptrb[i]))[7]]] = ptrb[i];
}
else std::sort(h,h+n,Mesh::GlobalIDComparator(this));
}
void Mesh::SortHandles(HandleType * h, enumerator n)
{
if( sizeof(HandleType) == 4 && n > 32 )
{
Element::adj_type::size_type i, num = static_cast<Element::adj_type::size_type>(n);
array<HandleType> tmp(num);
HandleType * ptra = h;
HandleType * ptrb = tmp.data();
//make so that invalid handles appear at he end
for (i = 0; i < num; i++) if( ptra[i] == InvalidHandle() ) ptra[i] = ENUMUNDEF;
const unsigned int kHist = 256;
unsigned int b0[kHist * 4];
unsigned int *b1 = b0 + kHist;
unsigned int *b2 = b0 + kHist*2;
unsigned int *b3 = b0 + kHist*3;
unsigned char * l1 = reinterpret_cast<unsigned char *>(ptra);
unsigned char * l2 = reinterpret_cast<unsigned char *>(ptrb);
memset(b0,0,sizeof(unsigned int)*kHist*4);
for (i = 0; i < num; i++)
{
++b0[l1[4*i+0]];
++b1[l1[4*i+1]];
++b2[l1[4*i+2]];
++b3[l1[4*i+3]];
}
unsigned int sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0, tsum;
for (i = 0; i < kHist; i++)
{
tsum = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = tsum;
tsum = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = tsum;
tsum = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = tsum;
tsum = b3[i] + sum3; b3[i] = sum3 - 1; sum3 = tsum;
}
for (i = 0; i < num; i++) ptrb[++b0[l1[4*i+0]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b1[l2[4*i+1]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b2[l1[4*i+2]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b3[l2[4*i+3]]] = ptrb[i];
for (i = 0; i < num; i++) if( ptra[i] == ENUMUNDEF ) ptra[i] = InvalidHandle();
}
else if( sizeof(HandleType) == 8 && n > 64 )
{
Element::adj_type::size_type i, num = static_cast<Element::adj_type::size_type>(n);
array<HandleType> tmp(num);
HandleType * ptra = h;
HandleType * ptrb = tmp.data();
//make so that invalid handles appear at he end
for (i = 0; i < num; i++) if( ptra[i] == InvalidHandle() ) ptra[i] = ENUMUNDEF;
const unsigned int kHist = 256;
unsigned int b0[kHist * 8];
unsigned int *b1 = b0 + kHist;
unsigned int *b2 = b0 + kHist*2;
unsigned int *b3 = b0 + kHist*3;
unsigned int *b4 = b0 + kHist*4;
unsigned int *b5 = b0 + kHist*5;
unsigned int *b6 = b0 + kHist*6;
unsigned int *b7 = b0 + kHist*7;
unsigned char * l1 = reinterpret_cast<unsigned char *>(ptra);
unsigned char * l2 = reinterpret_cast<unsigned char *>(ptrb);
memset(b0,0,sizeof(unsigned int)*kHist*8);
for (i = 0; i < num; i++)
{
++b0[l1[8*i+0]];
++b1[l1[8*i+1]];
++b2[l1[8*i+2]];
++b3[l1[8*i+3]];
++b4[l1[8*i+4]];
++b5[l1[8*i+5]];
++b6[l1[8*i+6]];
++b7[l1[8*i+7]];
}
unsigned int sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0, tsum;
unsigned int sum4 = 0, sum5 = 0, sum6 = 0, sum7 = 0;
for (i = 0; i < kHist; i++)
{
tsum = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = tsum;
tsum = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = tsum;
tsum = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = tsum;
tsum = b3[i] + sum3; b3[i] = sum3 - 1; sum3 = tsum;
tsum = b4[i] + sum4; b4[i] = sum4 - 1; sum4 = tsum;
tsum = b5[i] + sum5; b5[i] = sum5 - 1; sum5 = tsum;
tsum = b6[i] + sum6; b6[i] = sum6 - 1; sum6 = tsum;
tsum = b7[i] + sum7; b7[i] = sum7 - 1; sum7 = tsum;
}
for (i = 0; i < num; i++) ptrb[++b0[l1[8*i+0]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b1[l2[8*i+1]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b2[l1[8*i+2]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b3[l2[8*i+3]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b4[l1[8*i+4]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b5[l2[8*i+5]]] = ptrb[i];
for (i = 0; i < num; i++) ptrb[++b6[l1[8*i+6]]] = ptra[i];
for (i = 0; i < num; i++) ptra[++b7[l2[8*i+7]]] = ptrb[i];
for (i = 0; i < num; i++) if( ptra[i] == ENUMUNDEF ) ptra[i] = InvalidHandle();
}
else std::sort(h,h+n);
}
}
#endif
\ No newline at end of file
This diff is collapsed.
#include "inmost_mesh.h"
#if defined(USE_MESH)
namespace INMOST
{
template<typename StorageType>
void ElementArray<StorageType>::Unite(const HandleType * h, INMOST_DATA_ENUM_TYPE num)
{
assert(GetMeshLink() == other.GetMeshLink());
if( empty() ) container.insert(container.end(),h,h+num);
else
{
Mesh * m = GetMeshLink();
size_type s = size();
MarkerType mrk = m->CreateMarker();
SetMarker(mrk);
for(INMOST_DATA_ENUM_TYPE it = 0; it < num; it++)
if( !m->GetMarker(h[it],mrk) )
{
container.push_back(h[it]);
m->SetMarker(h[it],mrk);
}
RemMarker(mrk);
m->ReleaseMarker(mrk);
}
}
template<typename StorageType>
void ElementArray<StorageType>::Subtract(const HandleType * h, INMOST_DATA_ENUM_TYPE num)
{
assert(GetMeshLink() == other.GetMeshLink());
if( !empty() )
{
Mesh * mesh = GetMeshLink();
MarkerType mrk = mesh->CreateMarker();
//other.SetMarker(mrk);
mesh->SetMarkerArray(h,num,mrk);
{
size_type m = 0, n = 0;
while( m < size() )
{
if( !mesh->GetMarker(container[m],mrk) )
container[n++] = container[m];
m++;
}
container.resize(n);
}
//other.RemMarker(mrk);
mesh->RemMarkerArray(h,num,mrk);
mesh->ReleaseMarker(mrk);
}
}
template <typename StorageType>
void ElementArray<StorageType>::Intersect(const HandleType * h, INMOST_DATA_ENUM_TYPE num)
{
if( !empty() )
{
Mesh * mesh = GetMeshLink();
MarkerType mrk = mesh->CreateMarker();
//other.SetMarker(mrk);
mesh->SetMarkerArray(h,num,mrk);
{
size_type m = 0, n = 0;
while( m < size() )
{
if( mesh->GetMarker(container[m],mrk) )
container[n++] = container[m];
m++;
}
container.resize(n);
}
//other.RemMarker(mrk);
mesh->RemMarkerArray(h,num,mrk);
mesh->ReleaseMarker(mrk);
}
}
template <typename StorageType>
void ElementArray<StorageType>::SetMarker(MarkerType m) const
{
Mesh * mesh = GetMeshLink();
for(size_type it = 0; it < size(); it++) mesh->SetMarker(container[it],m);
}
template <typename StorageType>
void ElementArray<StorageType>::RemMarker(MarkerType m) const
{
Mesh * mesh = GetMeshLink();
for(size_type it = 0; it < size(); it++) mesh->RemMarker(container[it],m);
}
//all possible templates
template class ElementArray<Element>;
template class ElementArray<Node>;
template class ElementArray<Edge>;
template class ElementArray<Face>;
template class ElementArray<Cell>;
template class ElementArray<Storage>;
}
#endif
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
add_subdirectory(DrawGrid)
#add_subdirectory(DrawGrid)
add_subdirectory(OldDrawGrid)
add_subdirectory(DrawMatrix)
add_subdirectory(MatSolve)
add_subdirectory(GridGen)
add_subdirectory(FVDiscr)
add_subdirectory(OctreeCutcell)
add_subdirectory(Solver)
#add_subdirectory(DrawMatrix)
#add_subdirectory(MatSolve)
#add_subdirectory(GridGen)
#add_subdirectory(FVDiscr)
#add_subdirectory(OctreeCutcell)
#add_subdirectory(Solver)
......@@ -20,6 +20,7 @@ Mesh * ParallelCubeGenerator(INMOST_MPI_Comm comm, int nx, int ny, int nz)
m->SetCommunicator(comm); // Set the MPI communicator, usually MPI_COMM_WORLD
#if defined(USE_MPI)
MPI_Comm_set_errhandler(comm,MPI_ERRORS_RETURN);
//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
#endif
......@@ -125,6 +126,7 @@ Mesh * ParallelCubePrismGenerator(INMOST_MPI_Comm comm, int nx, int ny, int nz)
m->SetCommunicator(comm); // Set the MPI communicator, usually MPI_COMM_WORLD
#if defined(USE_MPI)
MPI_Comm_set_errhandler(comm,MPI_ERRORS_RETURN);
//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
#endif
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -263,7 +263,7 @@ namespace INMOST
//return data[i];
//for dynarray or array
for(unsigned it = 0; it < data.size(); ++it)
for(Entries::size_type it = 0; it < data.size(); ++it)
if( data[it].first == i ) return data[it].second;
entry new_entry;
new_entry.first = i;
......@@ -279,7 +279,7 @@ namespace INMOST
//for sparse data type
//return data[i];
for (unsigned it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
for (Entries::size_type it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
//you should not come here
assert(false);
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment