Initial draft of mesh partitioning capability that allows big
serial meshes to be partitioned and saved in parallel format using one processor. This capability allows MFEM to work around the current bottleneck which requires every MPI rank to load the big serial mesh before it can be partitioned. This new capability is added to the `mesh-explorer` miniapp with the new menu option `D` and is based on two new classes: * `MeshPartitioner`, constructed from a serial mesh and any partitioning array. Once constructed, it can extract any part of the mesh consisting of the elements with ids `elem_id` such that `partitioning[elem_id] == part_id`. The extracted mesh part is given in the form of a `MeshPart` object. * `MeshPart`, which is currently construced by a `MeshPartitioner`. In the future, it can be created by other methods to facilitate other capabilities such as parallel mesh re-partitioing. Once constructed, the `MeshPart` can be saved to a file using the text-based parallel MFEM format. Support for other formats can be added as well. Another capability that can be added is the MPI communication of `MeshPart` objects between different ranks which can be used, for example, for implementing a `ParMesh` constructor that needs the serial mesh only on one processor. Current limitations: * Non-conforming and NURBS meshes are not supported. * Meshes with nodes (e.g. high-order or periodic meshes) are not supported. Small extension: if `Mesh::SetCurvature` is called with `order = 0` then the method will remove the nodal `GridFunction` and use the vertices array instead. This "curvature removal" capability can be used from the `mesh-explorer` miniapp with the `c` menu option, by specifying 0 when prompted to enter order. Temporary change: `ex1p.cpp` is modified to disregard the (serial) mesh specified with the option `-m`/`--mesh`. Instead, it loads the parallel mesh `../miniapps/meshing/mesh-explorer.mesh.<rank>` produced by the `D` menu option of the `mesh-explorer` miniapp.
This commit is contained in:
Veselin Dobrev
+1
-1
@@ -234,7 +234,7 @@ miniapps/meshing/mobius-strip.mesh
|
||||
miniapps/meshing/klein-bottle.mesh
|
||||
miniapps/meshing/toroid-*.mesh
|
||||
miniapps/meshing/twist-*.mesh
|
||||
miniapps/meshing/mesh-explorer.mesh
|
||||
miniapps/meshing/mesh-explorer.mesh*
|
||||
miniapps/meshing/partitioning.txt
|
||||
miniapps/meshing/mesh-explorer-visit*
|
||||
miniapps/meshing/mesh-explorer-paraview/
|
||||
|
||||
+8
-2
@@ -147,10 +147,16 @@ int main(int argc, char *argv[])
|
||||
// 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
|
||||
// this mesh further in parallel to increase the resolution. Once the
|
||||
// parallel mesh is defined, the serial mesh can be deleted.
|
||||
ParMesh pmesh(MPI_COMM_WORLD, mesh);
|
||||
// ParMesh pmesh(MPI_COMM_WORLD, mesh);
|
||||
mesh.Clear();
|
||||
ifstream mesh_ifs(
|
||||
MakeParFilename("../miniapps/meshing/mesh-explorer.mesh.",
|
||||
myid));
|
||||
ParMesh pmesh(MPI_COMM_WORLD, mesh_ifs, /* refine: */ false);
|
||||
dim = pmesh.Dimension();
|
||||
pmesh.PrintInfo(cout);
|
||||
{
|
||||
int par_ref_levels = 2;
|
||||
int par_ref_levels = 0;
|
||||
for (int l = 0; l < par_ref_levels; l++)
|
||||
{
|
||||
pmesh.UniformRefinement();
|
||||
|
||||
+12
-8
@@ -74,10 +74,14 @@ public:
|
||||
inline Array(int asize, MemoryType mt)
|
||||
: size(asize) { asize > 0 ? data.New(asize, mt) : data.Reset(mt); }
|
||||
|
||||
/** @brief Creates array using an externally allocated pointer @a data_ to
|
||||
@a asize elements. The data pointer will not be deleted by Array. */
|
||||
inline Array(T *data_, int asize)
|
||||
{ data.Wrap(data_, asize, false); size = asize; }
|
||||
/** @brief Creates array using an externally allocated host pointer @a data_
|
||||
to @a asize elements. If @a own_data is true, the array takes ownership
|
||||
of the pointer.
|
||||
|
||||
When @a own_data is true, the pointer @a data_ must be allocated with
|
||||
MemoryType given by MemoryManager::GetHostMemoryType(). */
|
||||
inline Array(T *data_, int asize, bool own_data = false)
|
||||
{ data.Wrap(data_, asize, own_data); size = asize; }
|
||||
|
||||
/// Copy constructor: deep copy from @a src
|
||||
/** This method supports source arrays using any MemoryType. */
|
||||
@@ -205,7 +209,7 @@ public:
|
||||
inline void Copy(Array ©) const;
|
||||
|
||||
/// Make this Array a reference to a pointer.
|
||||
inline void MakeRef(T *, int);
|
||||
inline void MakeRef(T *data_, int size_, bool own_data = false);
|
||||
|
||||
/// Make this Array a reference to 'master'.
|
||||
inline void MakeRef(const Array &master);
|
||||
@@ -868,11 +872,11 @@ inline void Array<T>::Copy(Array ©) const
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline void Array<T>::MakeRef(T *p, int s)
|
||||
inline void Array<T>::MakeRef(T *data_, int size_, bool own_data)
|
||||
{
|
||||
data.Delete();
|
||||
data.Wrap(p, s, false);
|
||||
size = s;
|
||||
data.Wrap(data_, size_, own_data);
|
||||
size = size_;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
|
||||
@@ -275,7 +275,7 @@ void GroupTopology::Save(ostream &os) const
|
||||
os << "\ncommunication_groups\n";
|
||||
os << "number_of_groups " << NGroups() << "\n\n";
|
||||
|
||||
os << "# number of entities in each group, followed by group ids in group\n";
|
||||
os << "# number of entities in each group, followed by ranks in group\n";
|
||||
for (int group_id = 0; group_id < NGroups(); ++group_id)
|
||||
{
|
||||
int group_size = GetGroupSize(group_id);
|
||||
|
||||
@@ -208,6 +208,7 @@ void Transpose (const Table &A, Table &At, int ncols_A_ = -1);
|
||||
Table * Transpose (const Table &A);
|
||||
|
||||
/// Transpose an Array<int>
|
||||
/** @note The column (TYPE II) indices in each row of @a At will be sorted. */
|
||||
void Transpose(const Array<int> &A, Table &At, int ncols_A_ = -1);
|
||||
|
||||
/// C = A * B (as boolean matrices)
|
||||
|
||||
+741
@@ -30,6 +30,7 @@
|
||||
#include <cstring>
|
||||
#include <ctime>
|
||||
#include <functional>
|
||||
#include <map>
|
||||
#include <unordered_map>
|
||||
#include <unordered_set>
|
||||
|
||||
@@ -6198,6 +6199,12 @@ const FiniteElementSpace *Mesh::GetNodalFESpace() const
|
||||
|
||||
void Mesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
|
||||
{
|
||||
if (order <= 0)
|
||||
{
|
||||
delete Nodes;
|
||||
Nodes = nullptr;
|
||||
return;
|
||||
}
|
||||
space_dim = (space_dim == -1) ? spaceDim : space_dim;
|
||||
FiniteElementCollection* nfec;
|
||||
if (discont)
|
||||
@@ -13267,6 +13274,740 @@ void Mesh::GetGeometricParametersFromJacobian(const DenseMatrix &J,
|
||||
}
|
||||
|
||||
|
||||
void MeshPart::Print(std::ostream &out) const
|
||||
{
|
||||
out << "MFEM mesh v1.2\n";
|
||||
|
||||
// optional
|
||||
out <<
|
||||
"\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
|
||||
"# POINT = 0\n"
|
||||
"# SEGMENT = 1\n"
|
||||
"# TRIANGLE = 2\n"
|
||||
"# SQUARE = 3\n"
|
||||
"# TETRAHEDRON = 4\n"
|
||||
"# CUBE = 5\n"
|
||||
"# PRISM = 6\n"
|
||||
"# PYRAMID = 7\n"
|
||||
"#\n";
|
||||
|
||||
const int dim = dimension;
|
||||
out << "\ndimension\n" << dim;
|
||||
|
||||
out << "\n\nelements\n" << num_elements << '\n';
|
||||
if (element_map.Size() == 0)
|
||||
{
|
||||
int elem_id = 0;
|
||||
for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
|
||||
{
|
||||
const int nv = Geometry::NumVerts[g];
|
||||
const Array<int> &entity_vert = entity_to_vertex[g];
|
||||
for (int i = 0; i < entity_vert.Size(); i += nv)
|
||||
{
|
||||
out << attributes[elem_id++] << ' ' << g;
|
||||
for (int j = 0; j < nv; j++)
|
||||
{
|
||||
out << ' ' << entity_vert[i+j];
|
||||
}
|
||||
out << '\n';
|
||||
}
|
||||
}
|
||||
MFEM_ASSERT(elem_id == num_elements, "invalid MeshPart state");
|
||||
}
|
||||
else
|
||||
{
|
||||
MFEM_ASSERT(element_map.Size() == num_elements, "invalid MeshPart state");
|
||||
// Construct 'nat_to_bytype' -- the inverse of 'element_map'.
|
||||
Array<int> nat_to_bytype(num_elements);
|
||||
MFEM_DEBUG_DO(nat_to_bytype = -1);
|
||||
for (int bytype_elem_id = 0; bytype_elem_id < num_elements;
|
||||
bytype_elem_id++)
|
||||
{
|
||||
const int nat_elem_id = element_map[bytype_elem_id];
|
||||
MFEM_ASSERT(nat_to_bytype[nat_elem_id] == -1,
|
||||
"invalid 'element_map' array");
|
||||
nat_to_bytype[nat_elem_id] = bytype_elem_id;
|
||||
}
|
||||
int geom_offsets[Geometry::NumGeom+1];
|
||||
int geom_offset = 0;
|
||||
for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
|
||||
{
|
||||
geom_offsets[g] = geom_offset;
|
||||
geom_offset += entity_to_vertex[g].Size()/Geometry::NumVerts[g];
|
||||
}
|
||||
geom_offsets[Geometry::DimStart[dim+1]] = geom_offset;
|
||||
MFEM_ASSERT(geom_offset == num_elements, "invalid MeshPart state");
|
||||
for (int nat_elem_id = 0; nat_elem_id < num_elements; nat_elem_id++)
|
||||
{
|
||||
const int bytype_elem_id = nat_to_bytype[nat_elem_id];
|
||||
// Find the 'geom' that corresponds to 'bytype_elem_id'
|
||||
int geom = Geometry::DimStart[dim];
|
||||
while (geom_offsets[geom+1] <= bytype_elem_id) { geom++; }
|
||||
MFEM_ASSERT(geom < Geometry::NumGeom, "internal error");
|
||||
MFEM_ASSERT(Geometry::Dimension[geom] == dim, "internal error");
|
||||
const int nv = Geometry::NumVerts[geom];
|
||||
const int geom_elem_id = bytype_elem_id - geom_offsets[geom];
|
||||
const int *v = &entity_to_vertex[geom][nv*geom_elem_id];
|
||||
// Print the element
|
||||
out << attributes[nat_elem_id] << ' ' << geom;
|
||||
for (int i = 0; i < nv; i++)
|
||||
{
|
||||
out << ' ' << v[i];
|
||||
}
|
||||
out << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
out << "\nboundary\n" << num_bdr_elements << '\n';
|
||||
// TODO: it should be possible to handle the printing of the elements and the
|
||||
// boundary with one function.
|
||||
if (boundary_map.Size() == 0)
|
||||
{
|
||||
int bdr_elem_id = 0;
|
||||
for (int g = Geometry::DimStart[dim-1]; g < Geometry::DimStart[dim]; g++)
|
||||
{
|
||||
const int nv = Geometry::NumVerts[g];
|
||||
const Array<int> &entity_vert = entity_to_vertex[g];
|
||||
for (int i = 0; i < entity_vert.Size(); i += nv)
|
||||
{
|
||||
out << bdr_attributes[bdr_elem_id++] << ' ' << g;
|
||||
for (int j = 0; j < nv; j++)
|
||||
{
|
||||
out << ' ' << entity_vert[i+j];
|
||||
}
|
||||
out << '\n';
|
||||
}
|
||||
}
|
||||
MFEM_ASSERT(bdr_elem_id == num_bdr_elements, "invalid MeshPart state");
|
||||
}
|
||||
else
|
||||
{
|
||||
MFEM_ASSERT(boundary_map.Size() == num_bdr_elements,
|
||||
"invalid MeshPart state");
|
||||
// Construct 'nat_to_bytype' -- the inverse of 'boundary_map'.
|
||||
Array<int> nat_to_bytype(num_bdr_elements);
|
||||
MFEM_DEBUG_DO(nat_to_bytype = -1);
|
||||
for (int bytype_bdr_id = 0; bytype_bdr_id < num_bdr_elements;
|
||||
bytype_bdr_id++)
|
||||
{
|
||||
const int nat_bdr_id = boundary_map[bytype_bdr_id];
|
||||
MFEM_ASSERT(nat_to_bytype[nat_bdr_id] == -1,
|
||||
"invalid 'boundary_map' array");
|
||||
nat_to_bytype[nat_bdr_id] = bytype_bdr_id;
|
||||
}
|
||||
int geom_offsets[Geometry::NumGeom+1];
|
||||
int geom_offset = 0;
|
||||
for (int g = Geometry::DimStart[dim-1]; g < Geometry::DimStart[dim]; g++)
|
||||
{
|
||||
geom_offsets[g] = geom_offset;
|
||||
geom_offset += entity_to_vertex[g].Size()/Geometry::NumVerts[g];
|
||||
}
|
||||
geom_offsets[Geometry::DimStart[dim]] = geom_offset;
|
||||
MFEM_ASSERT(geom_offset == num_bdr_elements, "invalid MeshPart state");
|
||||
for (int nat_bdr_id = 0; nat_bdr_id < num_bdr_elements; nat_bdr_id++)
|
||||
{
|
||||
const int bytype_bdr_id = nat_to_bytype[nat_bdr_id];
|
||||
// Find the 'geom' that corresponds to 'bytype_bdr_id'
|
||||
int geom = Geometry::DimStart[dim-1];
|
||||
while (geom_offsets[geom+1] <= bytype_bdr_id) { geom++; }
|
||||
MFEM_ASSERT(geom < Geometry::NumGeom, "internal error");
|
||||
MFEM_ASSERT(Geometry::Dimension[geom] == dim-1, "internal error");
|
||||
const int nv = Geometry::NumVerts[geom];
|
||||
const int geom_bdr_id = bytype_bdr_id - geom_offsets[geom];
|
||||
const int *v = &entity_to_vertex[geom][nv*geom_bdr_id];
|
||||
// Print the boundary element
|
||||
out << bdr_attributes[nat_bdr_id] << ' ' << geom;
|
||||
for (int i = 0; i < nv; i++)
|
||||
{
|
||||
out << ' ' << v[i];
|
||||
}
|
||||
out << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
out << "\nvertices\n" << num_vertices << '\n';
|
||||
// if (Nodes == NULL)
|
||||
{
|
||||
const int sdim = space_dimension;
|
||||
out << sdim << '\n';
|
||||
for (int i = 0; i < num_vertices; i++)
|
||||
{
|
||||
out << vertex_coordinates[i*sdim];
|
||||
for (int d = 1; d < sdim; d++)
|
||||
{
|
||||
out << ' ' << vertex_coordinates[i*sdim+d];
|
||||
}
|
||||
out << '\n';
|
||||
}
|
||||
}
|
||||
// else
|
||||
// {
|
||||
// out << "\nnodes\n";
|
||||
// Nodes->Save(out);
|
||||
// }
|
||||
|
||||
out << "\nmfem_serial_mesh_end\n";
|
||||
|
||||
// Start: GroupTopology::Save
|
||||
const int num_groups = my_groups.Size();
|
||||
out << "\ncommunication_groups\n";
|
||||
out << "number_of_groups " << num_groups << "\n\n";
|
||||
|
||||
out << "# number of entities in each group, followed by ranks in group\n";
|
||||
for (int group_id = 0; group_id < num_groups; ++group_id)
|
||||
{
|
||||
const int group_size = my_groups.RowSize(group_id);
|
||||
const int *group_ptr = my_groups.GetRow(group_id);
|
||||
out << group_size;
|
||||
for (int group_member_index = 0; group_member_index < group_size;
|
||||
++group_member_index)
|
||||
{
|
||||
out << ' ' << group_ptr[group_member_index];
|
||||
}
|
||||
out << '\n';
|
||||
}
|
||||
// End: GroupTopology::Save
|
||||
|
||||
const Table &g2v = group__shared_entity_to_vertex[Geometry::POINT];
|
||||
const Table &g2ev = group__shared_entity_to_vertex[Geometry::SEGMENT];
|
||||
const Table &g2tv = group__shared_entity_to_vertex[Geometry::TRIANGLE];
|
||||
const Table &g2qv = group__shared_entity_to_vertex[Geometry::SQUARE];
|
||||
|
||||
MFEM_VERIFY(g2v.RowSize(0) == 0, "internal erroor");
|
||||
out << "\ntotal_shared_vertices " << g2v.Size_of_connections() << '\n';
|
||||
if (dimension >= 2)
|
||||
{
|
||||
MFEM_VERIFY(g2ev.RowSize(0) == 0, "internal erroor");
|
||||
out << "total_shared_edges " << g2ev.Size_of_connections()/2 << '\n';
|
||||
}
|
||||
if (dimension >= 3)
|
||||
{
|
||||
MFEM_VERIFY(g2tv.RowSize(0) == 0, "internal erroor");
|
||||
MFEM_VERIFY(g2qv.RowSize(0) == 0, "internal erroor");
|
||||
const int total_shared_faces =
|
||||
g2tv.Size_of_connections()/3 + g2qv.Size_of_connections()/4;
|
||||
out << "total_shared_faces " << total_shared_faces << '\n';
|
||||
}
|
||||
out << "\n# group 0 has no shared entities\n";
|
||||
for (int gr = 1; gr < num_groups; gr++)
|
||||
{
|
||||
{
|
||||
const int nv = g2v.RowSize(gr);
|
||||
const int *sv = g2v.GetRow(gr);
|
||||
out << "\n# group " << gr << "\nshared_vertices " << nv << '\n';
|
||||
for (int i = 0; i < nv; i++)
|
||||
{
|
||||
out << sv[i] << '\n';
|
||||
}
|
||||
}
|
||||
if (dimension >= 2)
|
||||
{
|
||||
const int ne = g2ev.RowSize(gr)/2;
|
||||
const int *se = g2ev.GetRow(gr);
|
||||
out << "\nshared_edges " << ne << '\n';
|
||||
for (int i = 0; i < ne; i++)
|
||||
{
|
||||
const int *v = se + 2*i;
|
||||
out << v[0] << ' ' << v[1] << '\n';
|
||||
}
|
||||
}
|
||||
if (dimension >= 3)
|
||||
{
|
||||
const int nt = g2tv.RowSize(gr)/3;
|
||||
const int *st = g2tv.GetRow(gr);
|
||||
const int nq = g2qv.RowSize(gr)/4;
|
||||
const int *sq = g2qv.GetRow(gr);
|
||||
out << "\nshared_faces " << nt+nq << '\n';
|
||||
for (int i = 0; i < nt; i++)
|
||||
{
|
||||
out << Geometry::TRIANGLE;
|
||||
const int *v = st + 3*i;
|
||||
for (int j = 0; j < 3; j++) { out << ' ' << v[j]; }
|
||||
out << '\n';
|
||||
}
|
||||
for (int i = 0; i < nq; i++)
|
||||
{
|
||||
out << Geometry::SQUARE;
|
||||
const int *v = sq + 4*i;
|
||||
for (int j = 0; j < 4; j++) { out << ' ' << v[j]; }
|
||||
out << '\n';
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Write out section end tag for mesh.
|
||||
out << "\nmfem_mesh_end" << endl;
|
||||
}
|
||||
|
||||
|
||||
MeshPartitioner::MeshPartitioner(Mesh &mesh_,
|
||||
int num_parts_,
|
||||
int *partitioning_,
|
||||
int part_method)
|
||||
: mesh(mesh_),
|
||||
partitioning(partitioning_),
|
||||
own_partitioning(false)
|
||||
{
|
||||
if (partitioning == nullptr)
|
||||
{
|
||||
partitioning = mesh.GeneratePartitioning(num_parts_, part_method);
|
||||
own_partitioning = true;
|
||||
}
|
||||
|
||||
Transpose(Array<int>(partitioning, mesh.GetNE()),
|
||||
part_to_element, num_parts_);
|
||||
// Note: the element ids in each row of 'part_to_element' are sorted.
|
||||
|
||||
const int dim = mesh.Dimension();
|
||||
if (dim >= 2)
|
||||
{
|
||||
Transpose(mesh.ElementToEdgeTable(), edge_to_element, mesh.GetNEdges());
|
||||
}
|
||||
|
||||
Array<int> boundary_to_part(mesh.GetNBE());
|
||||
// Same logic as in ParMesh::BuildLocalBoundary
|
||||
if (dim >= 3)
|
||||
{
|
||||
for (int i = 0; i < boundary_to_part.Size(); i++)
|
||||
{
|
||||
int face, o, el1, el2;
|
||||
mesh.GetBdrElementFace(i, &face, &o);
|
||||
mesh.GetFaceElements(face, &el1, &el2);
|
||||
boundary_to_part[i] =
|
||||
partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2];
|
||||
}
|
||||
}
|
||||
else if (dim == 2)
|
||||
{
|
||||
for (int i = 0; i < boundary_to_part.Size(); i++)
|
||||
{
|
||||
int edge = mesh.GetBdrElementEdgeIndex(i);
|
||||
int el1 = edge_to_element.GetRow(edge)[0];
|
||||
boundary_to_part[i] = partitioning[el1];
|
||||
}
|
||||
}
|
||||
else if (dim == 1)
|
||||
{
|
||||
for (int i = 0; i < boundary_to_part.Size(); i++)
|
||||
{
|
||||
int vert = mesh.GetBdrElementEdgeIndex(i);
|
||||
int el1, el2;
|
||||
mesh.GetFaceElements(vert, &el1, &el2);
|
||||
boundary_to_part[i] = partitioning[el1];
|
||||
}
|
||||
}
|
||||
Transpose(boundary_to_part, part_to_boundary, num_parts_);
|
||||
// Note: the boundary element ids in each row of 'part_to_boundary' are
|
||||
// sorted.
|
||||
boundary_to_part.DeleteAll();
|
||||
|
||||
Table *vert_element = mesh.GetVertexToElementTable(); // we must delete this
|
||||
vertex_to_element.Swap(*vert_element);
|
||||
delete vert_element;
|
||||
}
|
||||
|
||||
void MeshPartitioner::ExtractPart(int part_id, MeshPart &mesh_part) const
|
||||
{
|
||||
const int num_parts = part_to_element.Size();
|
||||
|
||||
MFEM_VERIFY(0 <= part_id && part_id < num_parts,
|
||||
"invalid part_id = " << part_id
|
||||
<< ", num_parts = " << num_parts);
|
||||
|
||||
const int dim = mesh.Dimension();
|
||||
const int sdim = mesh.SpaceDimension();
|
||||
const int num_elems = part_to_element.RowSize(part_id);
|
||||
const int *elem_list = part_to_element.GetRow(part_id);
|
||||
const int num_bdr_elems = part_to_boundary.RowSize(part_id);
|
||||
const int *bdr_elem_list = part_to_boundary.GetRow(part_id);
|
||||
|
||||
// Initialize 'mesh_part'
|
||||
mesh_part.dimension = dim;
|
||||
mesh_part.space_dimension = sdim;
|
||||
mesh_part.num_vertices = 0;
|
||||
mesh_part.num_elements = num_elems;
|
||||
mesh_part.num_bdr_elements = num_bdr_elems;
|
||||
for (int g = 0; g < Geometry::NumGeom; g++)
|
||||
{
|
||||
mesh_part.entity_to_vertex[g].SetSize(0); // can reuse Array allocation
|
||||
}
|
||||
mesh_part.element_map.SetSize(0); // 0 or 'num_elements', if needed
|
||||
mesh_part.boundary_map.SetSize(0); // 0 or 'num_bdr_elements', if needed
|
||||
mesh_part.attributes.SetSize(num_elems);
|
||||
mesh_part.bdr_attributes.SetSize(num_bdr_elems);
|
||||
mesh_part.vertex_coordinates.SetSize(0);
|
||||
|
||||
mesh_part.num_parts = num_parts;
|
||||
mesh_part.my_part_id = part_id;
|
||||
mesh_part.my_groups.Clear();
|
||||
for (int g = 0; g < Geometry::NumGeom; g++)
|
||||
{
|
||||
mesh_part.group__shared_entity_to_vertex[g].Clear();
|
||||
}
|
||||
|
||||
// Initialize:
|
||||
// - 'mesh_part.entity_to_vertex' for the elements (boundary elements are
|
||||
// set later); vertex ids are global at this point - they will be mapped to
|
||||
// local ids later
|
||||
// - 'mesh_part.attributes'
|
||||
int geom_marker = 0, num_geom = 0;
|
||||
for (int i = 0; i < num_elems; i++)
|
||||
{
|
||||
const Element *elem = mesh.GetElement(elem_list[i]);
|
||||
const int geom = elem->GetGeometryType();
|
||||
const int nv = Geometry::NumVerts[geom];
|
||||
const int *v = elem->GetVertices();
|
||||
MFEM_VERIFY(numeric_limits<int>::max() - nv >=
|
||||
mesh_part.entity_to_vertex[geom].Size(),
|
||||
"overflow in 'entity_to_vertex[geom]', geom: "
|
||||
<< Geometry::Name[geom]);
|
||||
mesh_part.entity_to_vertex[geom].Append(v, nv);
|
||||
mesh_part.attributes[i] = elem->GetAttribute();
|
||||
if ((geom_marker & (1 << geom)) == 0)
|
||||
{
|
||||
geom_marker |= (1 << geom);
|
||||
num_geom++;
|
||||
}
|
||||
}
|
||||
// Initialize 'mesh_part.element_map' if needed
|
||||
if (num_geom > 1)
|
||||
{
|
||||
int offsets[Geometry::NumGeom];
|
||||
int offset = 0;
|
||||
for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
|
||||
{
|
||||
offsets[g] = offset;
|
||||
offset += mesh_part.entity_to_vertex[g].Size()/Geometry::NumVerts[g];
|
||||
}
|
||||
mesh_part.element_map.SetSize(num_elems);
|
||||
for (int i = 0; i < num_elems; i++)
|
||||
{
|
||||
const int geom = mesh.GetElementGeometry(elem_list[i]);
|
||||
mesh_part.element_map[offsets[geom]++] = i;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize:
|
||||
// - 'mesh_part.entity_to_vertex' for the boundary elements; vertex ids are
|
||||
// global at this point - they will be mapped to local ids later
|
||||
// - 'mesh_part.bdr_attributes'
|
||||
geom_marker = 0; num_geom = 0;
|
||||
for (int i = 0; i < num_bdr_elems; i++)
|
||||
{
|
||||
const Element *bdr_elem = mesh.GetBdrElement(bdr_elem_list[i]);
|
||||
const int geom = bdr_elem->GetGeometryType();
|
||||
const int nv = Geometry::NumVerts[geom];
|
||||
const int *v = bdr_elem->GetVertices();
|
||||
MFEM_VERIFY(numeric_limits<int>::max() - nv >=
|
||||
mesh_part.entity_to_vertex[geom].Size(),
|
||||
"overflow in 'entity_to_vertex[geom]', geom: "
|
||||
<< Geometry::Name[geom]);
|
||||
mesh_part.entity_to_vertex[geom].Append(v, nv);
|
||||
mesh_part.bdr_attributes[i] = bdr_elem->GetAttribute();
|
||||
if ((geom_marker & (1 << geom)) == 0)
|
||||
{
|
||||
geom_marker |= (1 << geom);
|
||||
num_geom++;
|
||||
}
|
||||
}
|
||||
// Initialize 'mesh_part.boundary_map' if needed
|
||||
if (num_geom > 1)
|
||||
{
|
||||
int offsets[Geometry::NumGeom];
|
||||
int offset = 0;
|
||||
for (int g = Geometry::DimStart[dim-1]; g < Geometry::DimStart[dim]; g++)
|
||||
{
|
||||
offsets[g] = offset;
|
||||
offset += mesh_part.entity_to_vertex[g].Size()/Geometry::NumVerts[g];
|
||||
}
|
||||
mesh_part.boundary_map.SetSize(num_bdr_elems);
|
||||
for (int i = 0; i < num_bdr_elems; i++)
|
||||
{
|
||||
const int geom = mesh.GetBdrElementGeometry(bdr_elem_list[i]);
|
||||
mesh_part.boundary_map[offsets[geom]++] = i;
|
||||
}
|
||||
}
|
||||
|
||||
// Create the vertex id map, 'vertex_loc_to_glob', which maps local ids to
|
||||
// global ones; the map is sorted, preserving the global ordering.
|
||||
Array<int> vertex_loc_to_glob;
|
||||
{
|
||||
std::unordered_set<int> vertex_set;
|
||||
for (int i = 0; i < num_elems; i++)
|
||||
{
|
||||
const Element *elem = mesh.GetElement(elem_list[i]);
|
||||
const int geom = elem->GetGeometryType();
|
||||
const int nv = Geometry::NumVerts[geom];
|
||||
const int *v = elem->GetVertices();
|
||||
vertex_set.insert(v, v + nv);
|
||||
}
|
||||
vertex_loc_to_glob.SetSize(vertex_set.size());
|
||||
std::copy(vertex_set.begin(), vertex_set.end(), // src
|
||||
vertex_loc_to_glob.begin()); // dest
|
||||
}
|
||||
vertex_loc_to_glob.Sort();
|
||||
|
||||
// Initialize:
|
||||
// - 'mesh_part.num_vertices'
|
||||
// - 'mesh_part.vertex_coordinates', if needed
|
||||
mesh_part.num_vertices = vertex_loc_to_glob.Size();
|
||||
if (!mesh.GetNodes())
|
||||
{
|
||||
MFEM_VERIFY(numeric_limits<int>::max()/sdim >= vertex_loc_to_glob.Size(),
|
||||
"overflow in 'vertex_coordinates', num_vertices = "
|
||||
<< vertex_loc_to_glob.Size() << ", sdim = " << sdim);
|
||||
mesh_part.vertex_coordinates.SetSize(sdim*vertex_loc_to_glob.Size());
|
||||
for (int i = 0; i < vertex_loc_to_glob.Size(); i++)
|
||||
{
|
||||
const double *coord = mesh.GetVertex(vertex_loc_to_glob[i]);
|
||||
for (int d = 0; d < sdim; d++)
|
||||
{
|
||||
mesh_part.vertex_coordinates[i*sdim+d] = coord[d];
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
MFEM_ABORT("meshes with nodes are not implemented yet");
|
||||
}
|
||||
|
||||
// Update the vertex ids in the arrays 'mesh_part.entity_to_vertex' from
|
||||
// global to local.
|
||||
for (int g = 0; g < Geometry::NumGeom; g++)
|
||||
{
|
||||
Array<int> &vert_array = mesh_part.entity_to_vertex[g];
|
||||
for (int i = 0; i < vert_array.Size(); i++)
|
||||
{
|
||||
const int glob_id = vert_array[i];
|
||||
const int loc_id = vertex_loc_to_glob.FindSorted(glob_id);
|
||||
MFEM_ASSERT(loc_id >= 0, "internal error: global vertex id not found");
|
||||
vert_array[i] = loc_id;
|
||||
}
|
||||
}
|
||||
|
||||
// Begin constructing the "neighbor" groups, i.e. the groups that contain
|
||||
// 'part_id'.
|
||||
ListOfIntegerSets groups;
|
||||
{
|
||||
// the first group is the local one
|
||||
IntegerSet group;
|
||||
group.Recreate(1, &part_id);
|
||||
groups.Insert(group);
|
||||
}
|
||||
|
||||
// 'shared_faces' : shared face id -> (global_face_id, group_id)
|
||||
// Note: 'shared_faces' will be sorted by 'global_face_id'.
|
||||
Array<Pair<int,int>> shared_faces;
|
||||
|
||||
// Add "neighbor" groups defined by faces
|
||||
// Construct 'shared_faces'.
|
||||
if (dim >= 3)
|
||||
{
|
||||
std::unordered_set<int> face_set;
|
||||
// Construct 'face_set'
|
||||
const Table &elem_to_face = mesh.ElementToFaceTable();
|
||||
for (int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
|
||||
{
|
||||
const int glob_elem_id = elem_list[loc_elem_id];
|
||||
const int nfaces = elem_to_face.RowSize(glob_elem_id);
|
||||
const int *faces = elem_to_face.GetRow(glob_elem_id);
|
||||
face_set.insert(faces, faces + nfaces);
|
||||
}
|
||||
// Construct 'shared_faces'; add "neighbor" groups defined by faces.
|
||||
IntegerSet group;
|
||||
for (int glob_face_id : face_set)
|
||||
{
|
||||
int el[2];
|
||||
mesh.GetFaceElements(glob_face_id, &el[0], &el[1]);
|
||||
if (el[1] < 0) { continue; }
|
||||
el[0] = partitioning[el[0]];
|
||||
el[1] = partitioning[el[1]];
|
||||
MFEM_ASSERT(el[0] == part_id || el[1] == part_id, "internal error");
|
||||
if (el[0] != part_id || el[1] != part_id)
|
||||
{
|
||||
group.Recreate(2, el);
|
||||
const int group_id = groups.Insert(group);
|
||||
shared_faces.Append(Pair<int,int>(glob_face_id, group_id));
|
||||
}
|
||||
}
|
||||
shared_faces.Sort(); // sort the shared faces by 'glob_face_id'
|
||||
}
|
||||
|
||||
// 'shared_edges' : shared edge id -> (global_edge_id, group_id)
|
||||
// Note: 'shared_edges' will be sorted by 'global_edge_id'.
|
||||
Array<Pair<int,int>> shared_edges;
|
||||
|
||||
// Add "neighbor" groups defined by edges.
|
||||
// Construct 'shared_edges'.
|
||||
if (dim >= 2)
|
||||
{
|
||||
std::unordered_set<int> edge_set;
|
||||
// Construct 'edge_set'
|
||||
const Table &elem_to_edge = mesh.ElementToEdgeTable();
|
||||
for (int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
|
||||
{
|
||||
const int glob_elem_id = elem_list[loc_elem_id];
|
||||
const int nedges = elem_to_edge.RowSize(glob_elem_id);
|
||||
const int *edges = elem_to_edge.GetRow(glob_elem_id);
|
||||
edge_set.insert(edges, edges + nedges);
|
||||
}
|
||||
// Construct 'shared_edges'; add "neighbor" groups defined by edges.
|
||||
IntegerSet group;
|
||||
for (int glob_edge_id : edge_set)
|
||||
{
|
||||
const int nelem = edge_to_element.RowSize(glob_edge_id);
|
||||
const int *elem = edge_to_element.GetRow(glob_edge_id);
|
||||
Array<int> &gr = group; // reference to the 'group' internal Array
|
||||
gr.SetSize(nelem);
|
||||
for (int j = 0; j < nelem; j++)
|
||||
{
|
||||
gr[j] = partitioning[elem[j]];
|
||||
}
|
||||
gr.Sort();
|
||||
gr.Unique();
|
||||
MFEM_ASSERT(gr.FindSorted(part_id) >= 0, "internal error");
|
||||
if (group.Size() > 1)
|
||||
{
|
||||
const int group_id = groups.Insert(group);
|
||||
shared_edges.Append(Pair<int,int>(glob_edge_id, group_id));
|
||||
}
|
||||
}
|
||||
shared_edges.Sort(); // sort the shared edges by 'glob_edge_id'
|
||||
}
|
||||
|
||||
// 'shared_verts' : shared vertex id -> (global_vertex_id, group_id)
|
||||
// Note: 'shared_verts' will be sorted by 'global_vertex_id'.
|
||||
Array<Pair<int,int>> shared_verts;
|
||||
|
||||
// Add "neighbor" groups defined by vertices.
|
||||
// Construct 'shared_verts'.
|
||||
{
|
||||
IntegerSet group;
|
||||
for (int i = 0; i < vertex_loc_to_glob.Size(); i++)
|
||||
{
|
||||
// 'vertex_to_element' maps global vertex ids to global element ids
|
||||
const int glob_vertex_id = vertex_loc_to_glob[i];
|
||||
const int nelem = vertex_to_element.RowSize(glob_vertex_id);
|
||||
const int *elem = vertex_to_element.GetRow(glob_vertex_id);
|
||||
Array<int> &gr = group; // reference to the 'group' internal Array
|
||||
gr.SetSize(nelem);
|
||||
for (int j = 0; j < nelem; j++)
|
||||
{
|
||||
gr[j] = partitioning[elem[j]];
|
||||
}
|
||||
gr.Sort();
|
||||
gr.Unique();
|
||||
MFEM_ASSERT(gr.FindSorted(part_id) >= 0, "internal error");
|
||||
if (group.Size() > 1)
|
||||
{
|
||||
const int group_id = groups.Insert(group);
|
||||
shared_verts.Append(Pair<int,int>(glob_vertex_id, group_id));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Done constructing the "neighbor" groups in 'groups'.
|
||||
const int num_groups = groups.Size();
|
||||
|
||||
// Define 'mesh_part.my_groups'
|
||||
groups.AsTable(mesh_part.my_groups);
|
||||
|
||||
// Construct 'mesh_part.group__shared_entity_to_vertex[Geometry::POINT]'
|
||||
Table &group__shared_vertex_to_vertex =
|
||||
mesh_part.group__shared_entity_to_vertex[Geometry::POINT];
|
||||
group__shared_vertex_to_vertex.MakeI(num_groups);
|
||||
for (int sv = 0; sv < shared_verts.Size(); sv++)
|
||||
{
|
||||
const int group_id = shared_verts[sv].two;
|
||||
group__shared_vertex_to_vertex.AddAColumnInRow(group_id);
|
||||
}
|
||||
group__shared_vertex_to_vertex.MakeJ();
|
||||
for (int sv = 0; sv < shared_verts.Size(); sv++)
|
||||
{
|
||||
const int glob_vertex_id = shared_verts[sv].one;
|
||||
const int group_id = shared_verts[sv].two;
|
||||
const int loc_vertex_id = vertex_loc_to_glob.FindSorted(glob_vertex_id);
|
||||
MFEM_ASSERT(loc_vertex_id >= 0, "internal error");
|
||||
group__shared_vertex_to_vertex.AddConnection(group_id, loc_vertex_id);
|
||||
}
|
||||
group__shared_vertex_to_vertex.ShiftUpI();
|
||||
|
||||
// Construct 'mesh_part.group__shared_entity_to_vertex[Geometry::SEGMENT]'
|
||||
if (dim >= 2)
|
||||
{
|
||||
Table &group__shared_edge_to_vertex =
|
||||
mesh_part.group__shared_entity_to_vertex[Geometry::SEGMENT];
|
||||
group__shared_edge_to_vertex.MakeI(num_groups);
|
||||
for (int se = 0; se < shared_edges.Size(); se++)
|
||||
{
|
||||
const int group_id = shared_edges[se].two;
|
||||
group__shared_edge_to_vertex.AddColumnsInRow(group_id, 2);
|
||||
}
|
||||
group__shared_edge_to_vertex.MakeJ();
|
||||
const Table &edge_to_vertex = *mesh.GetEdgeVertexTable();
|
||||
for (int se = 0; se < shared_edges.Size(); se++)
|
||||
{
|
||||
const int glob_edge_id = shared_edges[se].one;
|
||||
const int group_id = shared_edges[se].two;
|
||||
const int *v = edge_to_vertex.GetRow(glob_edge_id);
|
||||
for (int i = 0; i < 2; i++)
|
||||
{
|
||||
const int loc_vertex_id = vertex_loc_to_glob.FindSorted(v[i]);
|
||||
MFEM_ASSERT(loc_vertex_id >= 0, "internal error");
|
||||
group__shared_edge_to_vertex.AddConnection(group_id, loc_vertex_id);
|
||||
}
|
||||
}
|
||||
group__shared_edge_to_vertex.ShiftUpI();
|
||||
}
|
||||
|
||||
// Construct 'mesh_part.group__shared_entity_to_vertex[Geometry::TRIANGLE]'
|
||||
// and 'mesh_part.group__shared_entity_to_vertex[Geometry::SQUARE]'.
|
||||
if (dim >= 3)
|
||||
{
|
||||
Table &group__shared_tria_to_vertex =
|
||||
mesh_part.group__shared_entity_to_vertex[Geometry::TRIANGLE];
|
||||
Table &group__shared_quad_to_vertex =
|
||||
mesh_part.group__shared_entity_to_vertex[Geometry::SQUARE];
|
||||
Array<int> vertex_ids;
|
||||
group__shared_tria_to_vertex.MakeI(num_groups);
|
||||
group__shared_quad_to_vertex.MakeI(num_groups);
|
||||
for (int sf = 0; sf < shared_faces.Size(); sf++)
|
||||
{
|
||||
const int glob_face_id = shared_faces[sf].one;
|
||||
const int group_id = shared_faces[sf].two;
|
||||
const int geom = mesh.GetFaceGeometry(glob_face_id);
|
||||
mesh_part.group__shared_entity_to_vertex[geom].
|
||||
AddColumnsInRow(group_id, Geometry::NumVerts[geom]);
|
||||
}
|
||||
group__shared_tria_to_vertex.MakeJ();
|
||||
group__shared_quad_to_vertex.MakeJ();
|
||||
for (int sf = 0; sf < shared_faces.Size(); sf++)
|
||||
{
|
||||
const int glob_face_id = shared_faces[sf].one;
|
||||
const int group_id = shared_faces[sf].two;
|
||||
const int geom = mesh.GetFaceGeometry(glob_face_id);
|
||||
mesh.GetFaceVertices(glob_face_id, vertex_ids);
|
||||
for (int i = 0; i < vertex_ids.Size(); i++)
|
||||
{
|
||||
const int glob_id = vertex_ids[i];
|
||||
const int loc_id = vertex_loc_to_glob.FindSorted(glob_id);
|
||||
MFEM_ASSERT(loc_id >= 0, "internal error");
|
||||
vertex_ids[i] = loc_id;
|
||||
}
|
||||
mesh_part.group__shared_entity_to_vertex[geom].
|
||||
AddConnections(group_id, vertex_ids, vertex_ids.Size());
|
||||
}
|
||||
group__shared_tria_to_vertex.ShiftUpI();
|
||||
group__shared_quad_to_vertex.ShiftUpI();
|
||||
}
|
||||
}
|
||||
|
||||
MeshPartitioner::~MeshPartitioner()
|
||||
{
|
||||
if (own_partitioning) { delete [] partitioning; }
|
||||
}
|
||||
|
||||
|
||||
GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
|
||||
int flags, MemoryType d_mt)
|
||||
{
|
||||
|
||||
+140
@@ -2436,6 +2436,146 @@ public:
|
||||
std::ostream &operator<<(std::ostream &out, const Mesh &mesh);
|
||||
|
||||
|
||||
// Class containing a minimal description of a part (a subset of the elements)
|
||||
// of a Mesh and its connectivity to other parts. The main purpose of this class
|
||||
// is to be communicated between MPI ranks for repartitioning purposes. It can
|
||||
// also be used to implement parallel mesh I/O functions with partitionings that
|
||||
// have number of parts different from the number of MPI tasks.
|
||||
//
|
||||
// Note: parts of NURBS or non-conforming meshes cannot be fully described by
|
||||
// this class alone.
|
||||
class MeshPart
|
||||
{
|
||||
public:
|
||||
// Reference space dimension of the elements
|
||||
int dimension;
|
||||
|
||||
// Dimension of the physical space into which the MeshPart is embedded.
|
||||
int space_dimension;
|
||||
|
||||
// Number of vertices
|
||||
int num_vertices;
|
||||
|
||||
// Number of elements with reference space dimension equal to 'dimension'.
|
||||
int num_elements;
|
||||
|
||||
// Number of boundary elements with reference space dimension equal to
|
||||
// 'dimension'-1.
|
||||
int num_bdr_elements;
|
||||
|
||||
// Each 'entity_to_vertex[geom]' describes the entities of Geometry::Type
|
||||
// 'geom' in terms of their vertices. The number of entities of type 'geom'
|
||||
// is:
|
||||
// num_entities[geom] = size('entity_to_vertex[geom]')/num_vertices[geom]
|
||||
// The number of all elements, 'num_elements', is:
|
||||
// 'num_elements' = sum_{dim[geom]=='dimension'} num_entities[geom]
|
||||
// and the number of all boundary elements, 'num_bdr_elements' is:
|
||||
// 'num_bdr_elements' = sum_{dim[geom]=='dimension'-1} num_entities[geom]
|
||||
// Note that 'entity_to_vertex' does NOT describe all "faces" in the mesh
|
||||
// part (i.e. all 'dimension'-1 entities) but only the boundary elements.
|
||||
Array<int> entity_to_vertex[Geometry::NumGeom];
|
||||
|
||||
// "By-type" element/boundary ordering: ordered by Geometry::Type and within
|
||||
// each Geometry::Type 'geom' ordered as in 'entity_to_vertex[geom]'.
|
||||
|
||||
// Optional re-ordering of the elements that will be used by (Par)Mesh
|
||||
// objects constructed from this MeshPart. This array maps "by-type" element
|
||||
// ids to desired/"natural" element ids:
|
||||
// "natural" element id = element_map["by-type" element id]
|
||||
// The size of the array is either 'num_elements' or 0 when no re-ordering is
|
||||
// needed (then "by-type" id == "natural" id).
|
||||
Array<int> element_map;
|
||||
|
||||
// Optional re-ordering for the boundary elements, similar to 'element_map'.
|
||||
Array<int> boundary_map;
|
||||
|
||||
// Element attributes. Ordered using the "natural" element ordering defined
|
||||
// by the array 'element_map'. The size of this array is 'num_elements'.
|
||||
Array<int> attributes;
|
||||
|
||||
// Boundary element attributes. Ordered using the "natural" boundary element
|
||||
// ordering defined by the array 'boundary_map'. The size of this array is
|
||||
// 'num_bdr_elements'.
|
||||
Array<int> bdr_attributes;
|
||||
|
||||
// Optional vertex coordinates. The size of the array is either
|
||||
// size = 'space_dimension' * 'num_vertices'
|
||||
// or 0 when the vertex coordinates are not used, i.e. when the MeshPart uses
|
||||
// a nodal GridFunction to describe its location in physical space. This
|
||||
// array uses Ordering::byVDIM: "X0,Y0,Z0, X1,Y1,Z1, ...".
|
||||
Array<double> vertex_coordinates;
|
||||
|
||||
// TODO
|
||||
// Add data members for the nodes.
|
||||
|
||||
// Connectivity to other MeshPart objects
|
||||
// --------------------------------------
|
||||
|
||||
// Total number of MeshParts
|
||||
int num_parts;
|
||||
|
||||
// Index of the part described by this MeshPart:
|
||||
// 0 <= 'my_part_id' < 'num_parts'
|
||||
int my_part_id;
|
||||
|
||||
// A group G is a subsets of the set { 0, 1, ..., 'num_parts'-1 } for which
|
||||
// there is a mesh entity E (of any dimension) in the global mesh such that
|
||||
// G is the set of the parts assigned to the elements adjacent to E. The
|
||||
// MeshPart describes only the "neighbor" groups, i.e. the groups that
|
||||
// contain 'my_part_id'. The Table 'my_groups' defines the "neighbor" groups
|
||||
// in terms of their part ids. In other words, it maps "neighbor" group ids
|
||||
// to a (sorted) list of part ids. In particular, the number of "neighbor"
|
||||
// groups is given by 'my_groups.Size()'. The "local" group { 'my_part_id' }
|
||||
// has index 0 in 'my_groups'.
|
||||
Table my_groups;
|
||||
|
||||
// Shared entities for this MeshPart are mesh entities of all dimensions less
|
||||
// than 'dimension' that are generated by the elements of this MeshPart and
|
||||
// at least one other MeshPart.
|
||||
//
|
||||
// The Table 'group__shared_entity_to_vertex[geom]' defines, for each group,
|
||||
// the shared entities of Geometry::Type 'geom'. Each row (corresponding to a
|
||||
// "neighbor" group, as defined by 'my_groups') in the Table defines the
|
||||
// shared entities in a way similar to the arrays 'entity_to_vertex[geom]'.
|
||||
// The "local" group (with index 0) does not have any shared entities, so the
|
||||
// 0-th row in the Table is always empty.
|
||||
//
|
||||
// IMPORTANT: the desciptions of the groups in this MeshPart must match their
|
||||
// descriptions in all neighboring MeshParts. This includes the ordering of
|
||||
// the shared entities within the group, as well as the vertex ordering of
|
||||
// each shared entity.
|
||||
Table group__shared_entity_to_vertex[Geometry::NumGeom];
|
||||
|
||||
// Write the MeshPart to a stream using the format "MFEM mesh v1.2".
|
||||
void Print(std::ostream &out) const;
|
||||
};
|
||||
|
||||
|
||||
// TODO: doocumentation
|
||||
class MeshPartitioner
|
||||
{
|
||||
protected:
|
||||
Mesh &mesh;
|
||||
int *partitioning;
|
||||
bool own_partitioning;
|
||||
Table part_to_element;
|
||||
Table part_to_boundary;
|
||||
Table edge_to_element;
|
||||
Table vertex_to_element;
|
||||
|
||||
public:
|
||||
// TODO: doocumentation
|
||||
MeshPartitioner(Mesh &mesh_, int num_parts_, int *partitioning_ = NULL,
|
||||
int part_method = 1);
|
||||
|
||||
// TODO: doocumentation
|
||||
void ExtractPart(int part_id, MeshPart &mesh_part) const;
|
||||
|
||||
// Destructor
|
||||
~MeshPartitioner();
|
||||
};
|
||||
|
||||
|
||||
/** @brief Structure for storing mesh geometric factors: coordinates, Jacobians,
|
||||
and determinants of the Jacobians. */
|
||||
/** Typically objects of this type are constructed and owned by objects of class
|
||||
|
||||
@@ -249,6 +249,8 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
|
||||
BuildSharedVertMapping(nsvert, vert_element, vert_global_local);
|
||||
delete vert_element;
|
||||
|
||||
// FIXME: the next two lines are already done above! Any reason to do them
|
||||
// again?
|
||||
SetMeshGen();
|
||||
meshgen = mesh.meshgen; // copy the global 'meshgen'
|
||||
}
|
||||
@@ -6305,6 +6307,7 @@ void ParMesh::ParPrint(ostream &os, const std::string &comments) const
|
||||
{
|
||||
os << "total_shared_faces " << sface_lface.Size() << '\n';
|
||||
}
|
||||
out << "\n# group 0 has no shared entities\n";
|
||||
for (int gr = 1; gr < GetNGroups(); gr++)
|
||||
{
|
||||
{
|
||||
|
||||
@@ -123,7 +123,7 @@ clean-build:
|
||||
rm -rf *.dSYM *.TVD.*breakpoints
|
||||
|
||||
clean-exec:
|
||||
@rm -f mobius-strip.mesh klein-bottle.mesh mesh-explorer.mesh
|
||||
@rm -f mobius-strip.mesh klein-bottle.mesh mesh-explorer.mesh*
|
||||
@rm -f toroid-*.mesh twist-*.mesh trimmer.mesh reflected.mesh
|
||||
@rm -f partitioning.txt shaper.mesh extruder.mesh
|
||||
@rm -f optimized* perturbed* polar-nc.mesh
|
||||
|
||||
@@ -308,6 +308,7 @@ int main (int argc, char *argv[])
|
||||
partitioning = 0;
|
||||
bdr_partitioning.SetSize(mesh->GetNBE());
|
||||
bdr_partitioning = 0;
|
||||
np = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
@@ -382,7 +383,8 @@ int main (int argc, char *argv[])
|
||||
"f) Find physical point in reference space\n"
|
||||
"p) Generate a partitioning\n"
|
||||
"o) Reorder elements\n"
|
||||
"S) Save in MFEM format\n"
|
||||
"S) Save in MFEM serial format\n"
|
||||
"D) Save in MFEM parallel format using the current partitioning\n"
|
||||
"V) Save in VTK format (only linear and quadratic meshes)\n"
|
||||
"D) Save as a DataCollection\n"
|
||||
"q) Quit\n"
|
||||
@@ -1037,7 +1039,7 @@ int main (int argc, char *argv[])
|
||||
partitioning.SetSize(mesh->GetNE());
|
||||
for (int i = 0; i < mesh->GetNE(); i++)
|
||||
{
|
||||
partitioning[i] = i * np / mesh->GetNE();
|
||||
partitioning[i] = (long long)i * np / mesh->GetNE();
|
||||
}
|
||||
recover_bdr_partitioning(mesh, partitioning, bdr_partitioning);
|
||||
}
|
||||
@@ -1250,6 +1252,25 @@ int main (int argc, char *argv[])
|
||||
cout << "New mesh file: " << omesh_file << endl;
|
||||
}
|
||||
|
||||
if (mk == 'D')
|
||||
{
|
||||
const char mesh_prefix[] = "mesh-explorer.mesh.";
|
||||
MeshPartitioner partitioner(*mesh, np, partitioning);
|
||||
MeshPart mesh_part;
|
||||
int precision;
|
||||
cout << "Enter desired precision: " << flush;
|
||||
cin >> precision;
|
||||
for (int i = 0; i < np; i++)
|
||||
{
|
||||
partitioner.ExtractPart(i, mesh_part);
|
||||
|
||||
ofstream omesh(MakeParFilename(mesh_prefix, i));
|
||||
omesh.precision(precision);
|
||||
mesh_part.Print(omesh);
|
||||
}
|
||||
cout << "New parallel mesh files: " << mesh_prefix << "<rank>" << endl;
|
||||
}
|
||||
|
||||
if (mk == 'V')
|
||||
{
|
||||
const char omesh_file[] = "mesh-explorer.vtk";
|
||||
|
||||
Reference in New Issue
Block a user