Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77953179
mesh_partition_scotch.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 17, 10:57
Size
14 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 10:57 (2 d)
Engine
blob
Format
Raw Data
Handle
19817973
Attached To
rAKA akantu
mesh_partition_scotch.cc
View Options
/**
* @file mesh_partition_scotch.cc
*
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Jan 22 2016
*
* @brief implementation of the MeshPartitionScotch class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh_partition_scotch.hh"
#include "aka_common.hh"
#include "aka_random_generator.hh"
#include "aka_static_if.hh"
#include "mesh_utils.hh"
/* -------------------------------------------------------------------------- */
#include <cstdio>
#include <fstream>
/* -------------------------------------------------------------------------- */
#if !defined(AKANTU_USE_PTSCOTCH)
#ifndef AKANTU_SCOTCH_NO_EXTERN
extern
"C"
{
#endif
// AKANTU_SCOTCH_NO_EXTERN
#include <scotch.h>
#ifndef AKANTU_SCOTCH_NO_EXTERN
}
#endif
// AKANTU_SCOTCH_NO_EXTERN
#else
// AKANTU_USE_PTSCOTCH
#include <ptscotch.h>
#endif
// AKANTU_USE_PTSCOTCH
namespace
akantu
{
namespace
{
constexpr
int
scotch_version
=
int
(
SCOTCH_VERSION
);
}
/* -------------------------------------------------------------------------- */
MeshPartitionScotch
::
MeshPartitionScotch
(
const
Mesh
&
mesh
,
UInt
spatial_dimension
,
const
ID
&
id
,
const
MemoryID
&
memory_id
)
:
MeshPartition
(
mesh
,
spatial_dimension
,
id
,
memory_id
)
{
AKANTU_DEBUG_IN
();
// check if the akantu types and Scotch one are consistent
static_assert
(
sizeof
(
Int
)
==
sizeof
(
SCOTCH_Num
),
"The integer type of Akantu does not match the one from Scotch"
);
// static_if(scotch_version >= 6)
// .then([](auto && y) {
// SCOTCH_randomSeed(y);
// })
// .else_([](auto && y) {
srandom
(
RandomGenerator
<
UInt
>::
seed
());
//y);
// })
// (std::forward<UInt>(RandomGenerator<UInt>::seed()));
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
static
SCOTCH_Mesh
*
createMesh
(
const
Mesh
&
mesh
)
{
AKANTU_DEBUG_IN
();
UInt
spatial_dimension
=
mesh
.
getSpatialDimension
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
UInt
total_nb_element
=
0
;
UInt
nb_edge
=
0
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
);
Mesh
::
type_iterator
end
=
mesh
.
lastType
(
spatial_dimension
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
UInt
nb_element
=
mesh
.
getNbElement
(
type
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
total_nb_element
+=
nb_element
;
nb_edge
+=
nb_element
*
nb_nodes_per_element
;
}
SCOTCH_Num
vnodbas
=
0
;
SCOTCH_Num
vnodnbr
=
nb_nodes
;
SCOTCH_Num
velmbas
=
vnodnbr
;
SCOTCH_Num
velmnbr
=
total_nb_element
;
SCOTCH_Num
*
verttab
=
new
SCOTCH_Num
[
vnodnbr
+
velmnbr
+
1
];
SCOTCH_Num
*
vendtab
=
verttab
+
1
;
SCOTCH_Num
*
velotab
=
NULL
;
SCOTCH_Num
*
vnlotab
=
NULL
;
SCOTCH_Num
*
vlbltab
=
NULL
;
memset
(
verttab
,
0
,
(
vnodnbr
+
velmnbr
+
1
)
*
sizeof
(
SCOTCH_Num
));
it
=
mesh
.
firstType
(
spatial_dimension
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
if
(
Mesh
::
getSpatialDimension
(
type
)
!=
spatial_dimension
)
continue
;
UInt
nb_element
=
mesh
.
getNbElement
(
type
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
,
_not_ghost
);
/// count number of occurrence of each node
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
)
{
UInt
*
conn_val
=
connectivity
.
storage
()
+
el
*
nb_nodes_per_element
;
for
(
UInt
n
=
0
;
n
<
nb_nodes_per_element
;
++
n
)
{
verttab
[
*
(
conn_val
++
)]
++
;
}
}
}
/// convert the occurrence array in a csr one
for
(
UInt
i
=
1
;
i
<
nb_nodes
;
++
i
)
verttab
[
i
]
+=
verttab
[
i
-
1
];
for
(
UInt
i
=
nb_nodes
;
i
>
0
;
--
i
)
verttab
[
i
]
=
verttab
[
i
-
1
];
verttab
[
0
]
=
0
;
/// rearrange element to get the node-element list
SCOTCH_Num
edgenbr
=
verttab
[
vnodnbr
]
+
nb_edge
;
SCOTCH_Num
*
edgetab
=
new
SCOTCH_Num
[
edgenbr
];
UInt
linearized_el
=
0
;
it
=
mesh
.
firstType
(
spatial_dimension
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
UInt
nb_element
=
mesh
.
getNbElement
(
type
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
,
_not_ghost
);
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
,
++
linearized_el
)
{
UInt
*
conn_val
=
connectivity
.
storage
()
+
el
*
nb_nodes_per_element
;
for
(
UInt
n
=
0
;
n
<
nb_nodes_per_element
;
++
n
)
edgetab
[
verttab
[
*
(
conn_val
++
)]
++
]
=
linearized_el
+
velmbas
;
}
}
for
(
UInt
i
=
nb_nodes
;
i
>
0
;
--
i
)
verttab
[
i
]
=
verttab
[
i
-
1
];
verttab
[
0
]
=
0
;
SCOTCH_Num
*
verttab_tmp
=
verttab
+
vnodnbr
+
1
;
SCOTCH_Num
*
edgetab_tmp
=
edgetab
+
verttab
[
vnodnbr
];
it
=
mesh
.
firstType
(
spatial_dimension
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
UInt
nb_element
=
mesh
.
getNbElement
(
type
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
,
_not_ghost
);
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
)
{
*
verttab_tmp
=
*
(
verttab_tmp
-
1
)
+
nb_nodes_per_element
;
verttab_tmp
++
;
UInt
*
conn
=
connectivity
.
storage
()
+
el
*
nb_nodes_per_element
;
for
(
UInt
i
=
0
;
i
<
nb_nodes_per_element
;
++
i
)
{
*
(
edgetab_tmp
++
)
=
*
(
conn
++
)
+
vnodbas
;
}
}
}
SCOTCH_Mesh
*
meshptr
=
new
SCOTCH_Mesh
;
SCOTCH_meshInit
(
meshptr
);
SCOTCH_meshBuild
(
meshptr
,
velmbas
,
vnodbas
,
velmnbr
,
vnodnbr
,
verttab
,
vendtab
,
velotab
,
vnlotab
,
vlbltab
,
edgenbr
,
edgetab
);
/// Check the mesh
AKANTU_DEBUG_ASSERT
(
SCOTCH_meshCheck
(
meshptr
)
==
0
,
"Scotch mesh is not consistent"
);
#ifndef AKANTU_NDEBUG
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
{
/// save initial graph
FILE
*
fmesh
=
fopen
(
"ScotchMesh.msh"
,
"w"
);
SCOTCH_meshSave
(
meshptr
,
fmesh
);
fclose
(
fmesh
);
/// write geometry file
std
::
ofstream
fgeominit
;
fgeominit
.
open
(
"ScotchMesh.xyz"
);
fgeominit
<<
spatial_dimension
<<
std
::
endl
<<
nb_nodes
<<
std
::
endl
;
const
Array
<
Real
>
&
nodes
=
mesh
.
getNodes
();
Real
*
nodes_val
=
nodes
.
storage
();
for
(
UInt
i
=
0
;
i
<
nb_nodes
;
++
i
)
{
fgeominit
<<
i
<<
" "
;
for
(
UInt
s
=
0
;
s
<
spatial_dimension
;
++
s
)
fgeominit
<<
*
(
nodes_val
++
)
<<
" "
;
fgeominit
<<
std
::
endl
;
;
}
fgeominit
.
close
();
}
#endif
AKANTU_DEBUG_OUT
();
return
meshptr
;
}
/* -------------------------------------------------------------------------- */
static
void
destroyMesh
(
SCOTCH_Mesh
*
meshptr
)
{
AKANTU_DEBUG_IN
();
SCOTCH_Num
velmbas
,
vnodbas
,
vnodnbr
,
velmnbr
,
*
verttab
,
*
vendtab
,
*
velotab
,
*
vnlotab
,
*
vlbltab
,
edgenbr
,
*
edgetab
,
degrptr
;
SCOTCH_meshData
(
meshptr
,
&
velmbas
,
&
vnodbas
,
&
velmnbr
,
&
vnodnbr
,
&
verttab
,
&
vendtab
,
&
velotab
,
&
vnlotab
,
&
vlbltab
,
&
edgenbr
,
&
edgetab
,
&
degrptr
);
delete
[]
verttab
;
delete
[]
edgetab
;
SCOTCH_meshExit
(
meshptr
);
delete
meshptr
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
MeshPartitionScotch
::
partitionate
(
UInt
nb_part
,
const
EdgeLoadFunctor
&
edge_load_func
,
const
Array
<
UInt
>
&
pairs
)
{
AKANTU_DEBUG_IN
();
nb_partitions
=
nb_part
;
tweakConnectivity
(
pairs
);
AKANTU_DEBUG_INFO
(
"Partitioning the mesh "
<<
mesh
.
getID
()
<<
" in "
<<
nb_part
<<
" parts."
);
Array
<
Int
>
dxadj
;
Array
<
Int
>
dadjncy
;
Array
<
Int
>
edge_loads
;
buildDualGraph
(
dxadj
,
dadjncy
,
edge_loads
,
edge_load_func
);
/// variables that will hold our structures in scotch format
SCOTCH_Graph
scotch_graph
;
SCOTCH_Strat
scotch_strat
;
/// description number and arrays for struct mesh for scotch
SCOTCH_Num
baseval
=
0
;
// base numbering for element and
// nodes (0 -> C , 1 -> fortran)
SCOTCH_Num
vertnbr
=
dxadj
.
size
()
-
1
;
// number of vertexes
SCOTCH_Num
*
parttab
;
// array of partitions
SCOTCH_Num
edgenbr
=
dxadj
(
vertnbr
);
// twice the number of "edges"
//(an "edge" bounds two nodes)
SCOTCH_Num
*
verttab
=
dxadj
.
storage
();
// array of start indices in edgetab
SCOTCH_Num
*
vendtab
=
NULL
;
// array of after-last indices in edgetab
SCOTCH_Num
*
velotab
=
NULL
;
// integer load associated with
// every vertex ( optional )
SCOTCH_Num
*
edlotab
=
edge_loads
.
storage
();
// integer load associated with
// every edge ( optional )
SCOTCH_Num
*
edgetab
=
dadjncy
.
storage
();
// adjacency array of every vertex
SCOTCH_Num
*
vlbltab
=
NULL
;
// vertex label array (optional)
/// Allocate space for Scotch arrays
parttab
=
new
SCOTCH_Num
[
vertnbr
];
/// Initialize the strategy structure
SCOTCH_stratInit
(
&
scotch_strat
);
/// Initialize the graph structure
SCOTCH_graphInit
(
&
scotch_graph
);
/// Build the graph from the adjacency arrays
SCOTCH_graphBuild
(
&
scotch_graph
,
baseval
,
vertnbr
,
verttab
,
vendtab
,
velotab
,
vlbltab
,
edgenbr
,
edgetab
,
edlotab
);
#ifndef AKANTU_NDEBUG
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
{
/// save initial graph
FILE
*
fgraphinit
=
fopen
(
"GraphIniFile.grf"
,
"w"
);
SCOTCH_graphSave
(
&
scotch_graph
,
fgraphinit
);
fclose
(
fgraphinit
);
/// write geometry file
std
::
ofstream
fgeominit
;
fgeominit
.
open
(
"GeomIniFile.xyz"
);
fgeominit
<<
spatial_dimension
<<
std
::
endl
<<
vertnbr
<<
std
::
endl
;
const
Array
<
Real
>
&
nodes
=
mesh
.
getNodes
();
Mesh
::
type_iterator
f_it
=
mesh
.
firstType
(
spatial_dimension
,
_not_ghost
,
_ek_not_defined
);
Mesh
::
type_iterator
f_end
=
mesh
.
lastType
(
spatial_dimension
,
_not_ghost
,
_ek_not_defined
);
Array
<
Real
>::
const_vector_iterator
nodes_it
=
nodes
.
begin
(
spatial_dimension
);
UInt
out_linerized_el
=
0
;
for
(;
f_it
!=
f_end
;
++
f_it
)
{
ElementType
type
=
*
f_it
;
UInt
nb_element
=
mesh
.
getNbElement
(
*
f_it
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
);
Vector
<
Real
>
mid
(
spatial_dimension
);
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
)
{
mid
.
set
(
0.
);
for
(
UInt
n
=
0
;
n
<
nb_nodes_per_element
;
++
n
)
{
UInt
node
=
connectivity
.
storage
()[
nb_nodes_per_element
*
el
+
n
];
mid
+=
Vector
<
Real
>
(
nodes_it
[
node
]);
}
mid
/=
nb_nodes_per_element
;
fgeominit
<<
out_linerized_el
++
<<
" "
;
for
(
UInt
s
=
0
;
s
<
spatial_dimension
;
++
s
)
fgeominit
<<
mid
[
s
]
<<
" "
;
fgeominit
<<
std
::
endl
;
;
}
}
fgeominit
.
close
();
}
#endif
/// Check the graph
AKANTU_DEBUG_ASSERT
(
SCOTCH_graphCheck
(
&
scotch_graph
)
==
0
,
"Graph to partition is not consistent"
);
/// Partition the mesh
SCOTCH_graphPart
(
&
scotch_graph
,
nb_part
,
&
scotch_strat
,
parttab
);
/// Check the graph
AKANTU_DEBUG_ASSERT
(
SCOTCH_graphCheck
(
&
scotch_graph
)
==
0
,
"Partitioned graph is not consistent"
);
#ifndef AKANTU_NDEBUG
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
{
/// save the partitioned graph
FILE
*
fgraph
=
fopen
(
"GraphFile.grf"
,
"w"
);
SCOTCH_graphSave
(
&
scotch_graph
,
fgraph
);
fclose
(
fgraph
);
/// save the partition map
std
::
ofstream
fmap
;
fmap
.
open
(
"MapFile.map"
);
fmap
<<
vertnbr
<<
std
::
endl
;
for
(
Int
i
=
0
;
i
<
vertnbr
;
i
++
)
fmap
<<
i
<<
" "
<<
parttab
[
i
]
<<
std
::
endl
;
fmap
.
close
();
}
#endif
/// free the scotch data structures
SCOTCH_stratExit
(
&
scotch_strat
);
SCOTCH_graphFree
(
&
scotch_graph
);
SCOTCH_graphExit
(
&
scotch_graph
);
fillPartitionInformation
(
mesh
,
parttab
);
delete
[]
parttab
;
restoreConnectivity
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
MeshPartitionScotch
::
reorder
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_INFO
(
"Reordering the mesh "
<<
mesh
.
getID
());
SCOTCH_Mesh
*
scotch_mesh
=
createMesh
(
mesh
);
UInt
nb_nodes
=
mesh
.
getNbNodes
();
SCOTCH_Strat
scotch_strat
;
// SCOTCH_Ordering scotch_order;
SCOTCH_Num
*
permtab
=
new
SCOTCH_Num
[
nb_nodes
];
SCOTCH_Num
*
peritab
=
NULL
;
SCOTCH_Num
cblknbr
=
0
;
SCOTCH_Num
*
rangtab
=
NULL
;
SCOTCH_Num
*
treetab
=
NULL
;
/// Initialize the strategy structure
SCOTCH_stratInit
(
&
scotch_strat
);
SCOTCH_Graph
scotch_graph
;
SCOTCH_graphInit
(
&
scotch_graph
);
SCOTCH_meshGraph
(
scotch_mesh
,
&
scotch_graph
);
#ifndef AKANTU_NDEBUG
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
{
FILE
*
fgraphinit
=
fopen
(
"ScotchMesh.grf"
,
"w"
);
SCOTCH_graphSave
(
&
scotch_graph
,
fgraphinit
);
fclose
(
fgraphinit
);
}
#endif
/// Check the graph
// AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0,
// "Mesh to Graph is not consistent");
SCOTCH_graphOrder
(
&
scotch_graph
,
&
scotch_strat
,
permtab
,
peritab
,
&
cblknbr
,
rangtab
,
treetab
);
SCOTCH_graphExit
(
&
scotch_graph
);
SCOTCH_stratExit
(
&
scotch_strat
);
destroyMesh
(
scotch_mesh
);
/// Renumbering
UInt
spatial_dimension
=
mesh
.
getSpatialDimension
();
for
(
UInt
g
=
_not_ghost
;
g
<=
_ghost
;
++
g
)
{
GhostType
gt
=
(
GhostType
)
g
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
_all_dimensions
,
gt
);
Mesh
::
type_iterator
end
=
mesh
.
lastType
(
_all_dimensions
,
gt
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
UInt
nb_element
=
mesh
.
getNbElement
(
type
,
gt
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
,
gt
);
UInt
*
conn
=
connectivity
.
storage
();
for
(
UInt
el
=
0
;
el
<
nb_element
*
nb_nodes_per_element
;
++
el
,
++
conn
)
{
*
conn
=
permtab
[
*
conn
];
}
}
}
/// \todo think of a in-place way to do it
Real
*
new_coordinates
=
new
Real
[
spatial_dimension
*
nb_nodes
];
Real
*
old_coordinates
=
mesh
.
getNodes
().
storage
();
for
(
UInt
i
=
0
;
i
<
nb_nodes
;
++
i
)
{
memcpy
(
new_coordinates
+
permtab
[
i
]
*
spatial_dimension
,
old_coordinates
+
i
*
spatial_dimension
,
spatial_dimension
*
sizeof
(
Real
));
}
memcpy
(
old_coordinates
,
new_coordinates
,
nb_nodes
*
spatial_dimension
*
sizeof
(
Real
));
delete
[]
new_coordinates
;
delete
[]
permtab
;
AKANTU_DEBUG_OUT
();
}
}
// namespace akantu
Event Timeline
Log In to Comment