phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into hypre-runtime-compute-policy

Browse Source
This commit is contained in:
Veselin Dobrev
2024-05-02 17:15:11 -07:00
committed by GitHub
parent 4ff3271a71 a379d5e92a
commit 5660111b37
14 changed files with 1321 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.gitignore
+1 -1
View File
@@ -232,7 +232,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/
CHANGELOG
+4
View File
@@ -17,6 +17,10 @@ Meshing improvements
arbitrary integer factors is also enabled, e.g. in the mesh-explorer miniapp.
NURBS coarsening and knot removal are also introduced.
- Added the capability to partition (big) serial meshes in serial code, see the
new classes MeshPartitioner and MeshPart. This capability is also exposed as a
menu option in the mesh-explorer miniapp in miniapps/meshing.
Discretization improvements
---------------------------
- Introduced support for higher order non conformal Nedelec elements on
config/defaults.mk
+1 -1
View File
@@ -388,7 +388,7 @@ GINKGO_LIB = $(XLINKER)-rpath,$(GINKGO_LINK_LIB_DIR) -L$(GINKGO_LINK_LIB_DIR)\
# AmgX library configuration
AMGX_DIR = @MFEM_DIR@/../amgx
AMGX_OPT = -I$(AMGX_DIR)/include
AMGX_LIB = -lcusparse -lcusolver -lcublas -lnvToolsExt -L$(AMGX_DIR)/lib -lamgx
AMGX_LIB = -L$(AMGX_DIR)/lib -lamgx -lcusparse -lcusolver -lcublas -lnvToolsExt
# GnuTLS library configuration
GNUTLS_OPT =
data/compass.mesh
+1
View File
@@ -92,4 +92,5 @@ vertices
-0.70710678 -0.70710678
0 -1
0.70710678 -0.70710678
mfem_mesh_end
general/array.hpp
+27 -8
View File
@@ -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,14 @@ public:
inline void Copy(Array &copy) const;
/// Make this Array a reference to a pointer.
inline void MakeRef(T *, int);
/** When @a own_data is true, the pointer @a data_ must be allocated with
MemoryType given by MemoryManager::GetHostMemoryType(). */
inline void MakeRef(T *data_, int size_, bool own_data = false);
/// Make this Array a reference to a pointer.
/** When @a own_data is true, the pointer @a data_ must be allocated with
MemoryType given by @a mt. */
inline void MakeRef(T *data_, int size, MemoryType mt, bool own_data);
/// Make this Array a reference to 'master'.
inline void MakeRef(const Array &master);
@@ -868,11 +879,19 @@ inline void Array<T>::Copy(Array &copy) 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>
inline void Array<T>::MakeRef(T *data_, int size_, MemoryType mt, bool own_data)
{
data.Delete();
data.Wrap(data_, size_, mt, own_data);
size = size_;
}
template <class T>
general/communication.cpp
+1 -1
View File
@@ -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);
general/table.hpp
+6 -1
View File
@@ -207,7 +207,12 @@ template <> inline void Swap<Table>(Table &a, Table &b)
void Transpose (const Table &A, Table &At, int ncols_A_ = -1);
Table * Transpose (const Table &A);
/// Transpose an Array<int>
/// @brief Transpose an Array<int>.
///
/// The array @a A represents a table where each row @a i has exactly one
/// connection to the column (TYPE II) index specified by @a A[i].
///
/// @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)
mesh/mesh.cpp
+884 -3
View File
@@ -20,9 +20,10 @@
#include "../general/tic_toc.hpp"
#include "../general/gecko.hpp"
#include "../general/kdtree.hpp"
#include "../general/sets.hpp"
#include "../fem/quadinterpolator.hpp"
#include <iostream>
// headers already included by mesh.hpp: <iostream>, <array>, <map>, <memory>
#include <sstream>
#include <fstream>
#include <limits>
@@ -1338,7 +1339,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
return res;
}
std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
std::ostream &operator<<(std::ostream &os, const Mesh::FaceInformation& info)
{
os << "face topology=";
switch (info.topology)
@@ -6209,6 +6210,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)
@@ -11419,7 +11426,8 @@ void Mesh::Printer(std::ostream &os, std::string section_delimiter,
if (!section_delimiter.empty())
{
os << section_delimiter << endl; // only with formats v1.2 and above
os << '\n'
<< section_delimiter << endl; // only with formats v1.2 and above
}
}
@@ -13299,6 +13307,879 @@ void Mesh::GetGeometricParametersFromJacobian(const DenseMatrix &J,
}
MeshPart::EntityHelper::EntityHelper(
int dim_, const Array<int> (&entity_to_vertex_)[Geometry::NumGeom])
: dim(dim_),
entity_to_vertex(entity_to_vertex_)
{
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;
num_entities = geom_offset;
}
MeshPart::Entity MeshPart::EntityHelper::FindEntity(int bytype_entity_id)
{
// Find the 'geom' that corresponds to 'bytype_entity_id'
int geom = Geometry::DimStart[dim];
while (geom_offsets[geom+1] <= bytype_entity_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_entity_id - geom_offsets[geom];
const int *v = &entity_to_vertex[geom][nv*geom_elem_id];
return { geom, nv, v };
}
void MeshPart::Print(std::ostream &os) const
{
os << "MFEM mesh v1.2\n";
// optional
os <<
"\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;
os << "\ndimension\n" << dim;
os << "\n\nelements\n" << num_elements << '\n';
{
const bool have_element_map = (element_map.Size() == num_elements);
MFEM_ASSERT(have_element_map || element_map.Size() == 0,
"invalid MeshPart state");
EntityHelper elem_helper(dim, entity_to_vertex);
MFEM_ASSERT(elem_helper.num_entities == num_elements,
"invalid MeshPart state");
for (int nat_elem_id = 0; nat_elem_id < num_elements; nat_elem_id++)
{
const int bytype_elem_id = have_element_map ?
element_map[nat_elem_id] : nat_elem_id;
const Entity ent = elem_helper.FindEntity(bytype_elem_id);
// Print the element
os << attributes[nat_elem_id] << ' ' << ent.geom;
for (int i = 0; i < ent.num_verts; i++)
{
os << ' ' << ent.verts[i];
}
os << '\n';
}
}
os << "\nboundary\n" << num_bdr_elements << '\n';
{
const bool have_boundary_map = (boundary_map.Size() == num_bdr_elements);
MFEM_ASSERT(have_boundary_map || boundary_map.Size() == 0,
"invalid MeshPart state");
EntityHelper bdr_helper(dim-1, entity_to_vertex);
MFEM_ASSERT(bdr_helper.num_entities == 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 = have_boundary_map ?
boundary_map[nat_bdr_id] : nat_bdr_id;
const Entity ent = bdr_helper.FindEntity(bytype_bdr_id);
// Print the boundary element
os << bdr_attributes[nat_bdr_id] << ' ' << ent.geom;
for (int i = 0; i < ent.num_verts; i++)
{
os << ' ' << ent.verts[i];
}
os << '\n';
}
}
os << "\nvertices\n" << num_vertices << '\n';
if (!nodes)
{
const int sdim = space_dimension;
os << sdim << '\n';
for (int i = 0; i < num_vertices; i++)
{
os << vertex_coordinates[i*sdim];
for (int d = 1; d < sdim; d++)
{
os << ' ' << vertex_coordinates[i*sdim+d];
}
os << '\n';
}
}
else
{
os << "\nnodes\n";
nodes->Save(os);
}
os << "\nmfem_serial_mesh_end\n";
// Start: GroupTopology::Save
const int num_groups = my_groups.Size();
os << "\ncommunication_groups\n";
os << "number_of_groups " << num_groups << "\n\n";
os << "# 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);
os << group_size;
for (int group_member_index = 0; group_member_index < group_size;
++group_member_index)
{
os << ' ' << group_ptr[group_member_index];
}
os << '\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");
os << "\ntotal_shared_vertices " << g2v.Size_of_connections() << '\n';
if (dimension >= 2)
{
MFEM_VERIFY(g2ev.RowSize(0) == 0, "internal erroor");
os << "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;
os << "total_shared_faces " << total_shared_faces << '\n';
}
os << "\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);
os << "\n# group " << gr << "\nshared_vertices " << nv << '\n';
for (int i = 0; i < nv; i++)
{
os << sv[i] << '\n';
}
}
if (dimension >= 2)
{
const int ne = g2ev.RowSize(gr)/2;
const int *se = g2ev.GetRow(gr);
os << "\nshared_edges " << ne << '\n';
for (int i = 0; i < ne; i++)
{
const int *v = se + 2*i;
os << 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);
os << "\nshared_faces " << nt+nq << '\n';
for (int i = 0; i < nt; i++)
{
os << Geometry::TRIANGLE;
const int *v = st + 3*i;
for (int j = 0; j < 3; j++) { os << ' ' << v[j]; }
os << '\n';
}
for (int i = 0; i < nq; i++)
{
os << Geometry::SQUARE;
const int *v = sq + 4*i;
for (int j = 0; j < 4; j++) { os << ' ' << v[j]; }
os << '\n';
}
}
}
// Write out section end tag for mesh.
os << "\nmfem_mesh_end" << endl;
}
Mesh &MeshPart::GetMesh()
{
if (mesh) { return *mesh; }
mesh.reset(new Mesh(dimension,
num_vertices,
num_elements,
num_bdr_elements,
space_dimension));
// Add elements
{
const bool have_element_map = (element_map.Size() == num_elements);
MFEM_ASSERT(have_element_map || element_map.Size() == 0,
"invalid MeshPart state");
EntityHelper elem_helper(dimension, entity_to_vertex);
MFEM_ASSERT(elem_helper.num_entities == num_elements,
"invalid MeshPart state");
const bool have_tet_refine_flags = (tet_refine_flags.Size() > 0);
for (int nat_elem_id = 0; nat_elem_id < num_elements; nat_elem_id++)
{
const int bytype_elem_id = have_element_map ?
element_map[nat_elem_id] : nat_elem_id;
const Entity ent = elem_helper.FindEntity(bytype_elem_id);
Element *el = mesh->NewElement(ent.geom);
el->SetVertices(ent.verts);
el->SetAttribute(attributes[nat_elem_id]);
if (ent.geom == Geometry::TETRAHEDRON && have_tet_refine_flags)
{
constexpr int geom_tet = Geometry::TETRAHEDRON;
const int tet_id = (ent.verts - entity_to_vertex[geom_tet])/4;
const int ref_flag = tet_refine_flags[tet_id];
static_cast<Tetrahedron*>(el)->SetRefinementFlag(ref_flag);
}
mesh->AddElement(el);
}
}
// Add boundary elements
{
const bool have_boundary_map = (boundary_map.Size() == num_bdr_elements);
MFEM_ASSERT(have_boundary_map || boundary_map.Size() == 0,
"invalid MeshPart state");
EntityHelper bdr_helper(dimension-1, entity_to_vertex);
MFEM_ASSERT(bdr_helper.num_entities == 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 = have_boundary_map ?
boundary_map[nat_bdr_id] : nat_bdr_id;
const Entity ent = bdr_helper.FindEntity(bytype_bdr_id);
Element *bdr = mesh->NewElement(ent.geom);
bdr->SetVertices(ent.verts);
bdr->SetAttribute(bdr_attributes[nat_bdr_id]);
mesh->AddBdrElement(bdr);
}
}
// Add vertices
if (vertex_coordinates.Size() == space_dimension*num_vertices)
{
MFEM_ASSERT(!nodes, "invalid MeshPart state");
for (int vert_id = 0; vert_id < num_vertices; vert_id++)
{
mesh->AddVertex(vertex_coordinates + space_dimension*vert_id);
}
}
else
{
MFEM_ASSERT(vertex_coordinates.Size() == 0, "invalid MeshPart state");
for (int vert_id = 0; vert_id < num_vertices; vert_id++)
{
mesh->AddVertex(0., 0., 0.);
}
// 'mesh.Nodes' cannot be set here -- they can be set later, if needed
}
mesh->FinalizeTopology(/* generate_bdr: */ false);
return *mesh;
}
MeshPartitioner::MeshPartitioner(Mesh &mesh_,
int num_parts_,
int *partitioning_,
int part_method)
: mesh(mesh_)
{
if (partitioning_)
{
partitioning.MakeRef(partitioning_, mesh.GetNE(), false);
}
else
{
partitioning_ = mesh.GeneratePartitioning(num_parts_, part_method);
// Mesh::GeneratePartitioning always uses new[] to allocate the,
// partitioning, so we need to tell the memory manager to free it with
// delete[] (even if a different host memory type has been selected).
const MemoryType mt = MemoryType::HOST;
partitioning.MakeRef(partitioning_, mesh.GetNE(), mt, true);
}
Transpose(partitioning, 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.GetBdrElementFaceIndex(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.GetBdrElementFaceIndex(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); // sorted
const int num_bdr_elems = part_to_boundary.RowSize(part_id);
const int *bdr_elem_list = part_to_boundary.GetRow(part_id); // sorted
// 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.tet_refine_flags.SetSize(0);
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();
}
mesh_part.nodes.reset(nullptr);
mesh_part.nodal_fes.reset(nullptr);
mesh_part.mesh.reset(nullptr);
// 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'
// - 'mesh_part.tet_refine_flags' if needed
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 == Geometry::TETRAHEDRON)
{
// Create 'mesh_part.tet_refine_flags' but only if we find at least one
// non-zero flag in a tetrahedron.
const Tetrahedron *tet = static_cast<const Tetrahedron*>(elem);
const int ref_flag = tet->GetRefinementFlag();
if (mesh_part.tet_refine_flags.Size() == 0)
{
if (ref_flag)
{
// This is the first time we encounter non-zero 'ref_flag'
const int num_tets = mesh_part.entity_to_vertex[geom].Size()/nv;
mesh_part.tet_refine_flags.SetSize(num_tets, 0);
mesh_part.tet_refine_flags.Last() = ref_flag;
}
}
else
{
mesh_part.tet_refine_flags.Append(ref_flag);
}
}
if ((geom_marker & (1 << geom)) == 0)
{
geom_marker |= (1 << geom);
num_geom++;
}
}
MFEM_ASSERT(mesh_part.tet_refine_flags.Size() == 0 ||
mesh_part.tet_refine_flags.Size() ==
mesh_part.entity_to_vertex[Geometry::TETRAHEDRON].Size()/4,
"internal error");
// 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[i] = offsets[geom]++;
}
}
// 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[i] = offsets[geom]++;
}
}
// 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.num_vertices = vertex_loc_to_glob.Size();
// 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;
}
}
// Initialize one of 'mesh_part.vertex_coordinates' or 'mesh_part.nodes'
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 real_t *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
{
const GridFunction &glob_nodes = *mesh.GetNodes();
mesh_part.nodal_fes = ExtractFESpace(mesh_part, *glob_nodes.FESpace());
// Initialized 'mesh_part.mesh'.
// Note: the nodes of 'mesh_part.mesh' are not set.
mesh_part.nodes = ExtractGridFunction(mesh_part, glob_nodes,
*mesh_part.nodal_fes);
// Attach the 'mesh_part.nodes' to the 'mesh_part.mesh'.
mesh_part.mesh->NewNodes(*mesh_part.nodes, /* make_owner: */ false);
// Note: the vertices of 'mesh_part.mesh' are not set.
}
// 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);
// Rotate shared triangles that have an adjacent tetrahedron with a
// nonzero refinement flag.
// See also ParMesh::BuildSharedFaceElems.
if (geom == Geometry::TRIANGLE)
{
int glob_el_id[2];
mesh.GetFaceElements(glob_face_id, &glob_el_id[0], &glob_el_id[1]);
int side = 0;
const Element *el = mesh.GetElement(glob_el_id[0]);
const Tetrahedron *tet = nullptr;
if (el->GetGeometryType() == Geometry::TETRAHEDRON)
{
tet = static_cast<const Tetrahedron*>(el);
}
else
{
side = 1;
el = mesh.GetElement(glob_el_id[1]);
if (el->GetGeometryType() == Geometry::TETRAHEDRON)
{
tet = static_cast<const Tetrahedron*>(el);
}
}
if (tet && tet->GetRefinementFlag())
{
// mark the shared face for refinement by reorienting
// it according to the refinement flag in the tetrahedron
// to which this shared face belongs to.
int info[2];
mesh.GetFaceInfos(glob_face_id, &info[0], &info[1]);
tet->GetMarkedFace(info[side]/64, &vertex_ids[0]);
}
}
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();
}
}
std::unique_ptr<FiniteElementSpace>
MeshPartitioner::ExtractFESpace(MeshPart &mesh_part,
const FiniteElementSpace &global_fespace) const
{
mesh_part.GetMesh(); // initialize 'mesh_part.mesh'
// Note: the nodes of 'mesh_part.mesh' are not set by GetMesh() unless they
// were already constructed, e.g. by ExtractPart().
return std::unique_ptr<FiniteElementSpace>(
new FiniteElementSpace(mesh_part.mesh.get(),
global_fespace.FEColl(),
global_fespace.GetVDim(),
global_fespace.GetOrdering()));
}
std::unique_ptr<GridFunction>
MeshPartitioner::ExtractGridFunction(const MeshPart &mesh_part,
const GridFunction &global_gf,
FiniteElementSpace &local_fespace) const
{
std::unique_ptr<GridFunction> local_gf(new GridFunction(&local_fespace));
// Transfer data from 'global_gf' to 'local_gf'.
Array<int> gvdofs, lvdofs;
Vector loc_vals;
const int part_id = mesh_part.my_part_id;
const int num_elems = part_to_element.RowSize(part_id);
const int *elem_list = part_to_element.GetRow(part_id); // sorted
for (int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
{
const int glob_elem_id = elem_list[loc_elem_id];
auto glob_dt = global_gf.FESpace()->GetElementVDofs(glob_elem_id, gvdofs);
global_gf.GetSubVector(gvdofs, loc_vals);
if (glob_dt) { glob_dt->InvTransformPrimal(loc_vals); }
auto local_dt = local_fespace.GetElementVDofs(loc_elem_id, lvdofs);
if (local_dt) { local_dt->TransformPrimal(loc_vals); }
local_gf->SetSubVector(lvdofs, loc_vals);
}
return local_gf;
}
GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
int flags, MemoryType d_mt)
{
mesh/mesh.hpp
+349 -17
View File
@@ -30,6 +30,7 @@
#include <iostream>
#include <array>
#include <map>
#include <memory>
namespace mfem
{
@@ -75,8 +76,10 @@ protected:
visualization purpose in GLVis. */
mutable int nbInteriorFaces, nbBoundaryFaces;
int meshgen; // see MeshGenerator()
int mesh_geoms; // sum of (1 << geom) for all geom of all dimensions
// see MeshGenerator(); global in parallel
int meshgen;
// sum of (1 << geom) for all geom of all dimensions; local in parallel
int mesh_geoms;
// Counter for Mesh transformations: refinement, derefinement, rebalancing.
// Used for checking during Update operations on objects depending on the
@@ -307,11 +310,11 @@ protected:
void Destroy(); // Delete all owned data.
void ResetLazyData();
Element *ReadElementWithoutAttr(std::istream &);
static void PrintElementWithoutAttr(const Element *, std::ostream &);
Element *ReadElementWithoutAttr(std::istream &input);
static void PrintElementWithoutAttr(const Element *el, std::ostream &os);
Element *ReadElement(std::istream &);
static void PrintElement(const Element *, std::ostream &);
Element *ReadElement(std::istream &input);
static void PrintElement(const Element *el, std::ostream &os);
// Readers for different mesh formats, used in the Load() method.
// The implementations of these methods are in mesh_readers.cpp.
@@ -558,7 +561,7 @@ protected:
mfem v1.2 format with the given section_delimiter at the end.
If @a comments is non-empty, it will be printed after the first line of
the file, and each line should begin with '#'. */
void Printer(std::ostream &out = mfem::out,
void Printer(std::ostream &os = mfem::out,
std::string section_delimiter = "",
const std::string &comments = "") const;
@@ -2124,7 +2127,11 @@ public:
/// Set the curvature of the mesh nodes using the given polynomial degree.
/** Creates a nodal GridFunction if one doesn't already exist.
@param[in] order Polynomial degree of the nodal FE space.
@param[in] order Polynomial degree of the nodal FE space. If this
value is <= 0 then the method will remove the
nodal GridFunction and the Mesh will use the
vertices array instead; the other arguments are
ignored in this case.
@param[in] discont Whether to use a discontinuous or continuous
finite element space (continuous is default).
@param[in] space_dim The space dimension (optional).
@@ -2330,7 +2337,7 @@ public:
std::ostream &os, int elem_attr = 0) const;
void PrintElementsWithPartitioning (int *partitioning,
std::ostream &out,
std::ostream &os,
int interior_faces = 0);
/// Print set of disjoint surfaces:
@@ -2338,13 +2345,13 @@ public:
* If Aface_face(i,j) != 0, print face j as a boundary
* element with attribute i+1.
*/
void PrintSurfaces(const Table &Aface_face, std::ostream &out) const;
void PrintSurfaces(const Table &Aface_face, std::ostream &os) const;
/// Auxiliary method used by PrintCharacteristics().
/** It is also used in the `mesh-explorer` miniapp. */
static void PrintElementsByGeometry(int dim,
const Array<int> &num_elems_by_geom,
std::ostream &out);
std::ostream &os);
/** @brief Compute and print mesh characteristics such as number of vertices,
number of elements, number of boundary elements, minimal and maximal
@@ -2364,7 +2371,7 @@ public:
#ifdef MFEM_DEBUG
/// Output an NCMesh-compatible debug dump.
void DebugDump(std::ostream &out) const;
void DebugDump(std::ostream &os) const;
#endif
/// @}
@@ -2445,7 +2452,334 @@ public:
/** Overload operator<< for std::ostream and Mesh; valid also for the derived
class ParMesh */
std::ostream &operator<<(std::ostream &out, const Mesh &mesh);
std::ostream &operator<<(std::ostream &os, const Mesh &mesh);
/// @brief Print function for Mesh::FaceInformation.
std::ostream& operator<<(std::ostream &os, const Mesh::FaceInformation& info);
/** @brief 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 facilitate the partitioning of serial
meshes (in serial, i.e. on one processor) and save the parts in parallel
MFEM mesh format.
Another potential futrure purpose of this class could be to facilitate
exchange of MeshParts between MPI ranks for repartitioning purposes. It can
also potentially 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 with its current data members. Such extensions may be added
in the future.
*/
class MeshPart
{
protected:
struct Entity { int geom; int num_verts; const int *verts; };
struct EntityHelper
{
int dim, num_entities;
int geom_offsets[Geometry::NumGeom+1];
typedef const Array<int> entity_to_vertex_type[Geometry::NumGeom];
entity_to_vertex_type &entity_to_vertex;
EntityHelper(int dim_,
const Array<int> (&entity_to_vertex_)[Geometry::NumGeom]);
Entity FindEntity(int bytype_entity_id);
};
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;
/** @brief 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.
Also, note that lower dimesional entities ('dimension'-2 and lower) are
NOT described by the respective array, i.e. the array will be empty.
*/
Array<int> entity_to_vertex[Geometry::NumGeom];
/** @brief Store the refinement flags for tetraheral elements. If all tets
have zero refinement flags then this array is empty, i.e. has size 0. */
Array<int> tet_refine_flags;
/**
Terminology: "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 "natural" element
ids (used by the Mesh/ParMesh objects) to "by-type" element ids (see
above):
"by-type" element id = element_map["natural" 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<real_t> vertex_coordinates;
/**
Optional serial Mesh object constructed on demand using the method
GetMesh(). One use case for it is when one wants to construct FE spaces
and GridFunction%s on the MeshPart for saving or MPI communication.
*/
std::unique_ptr<Mesh> mesh;
/**
Nodal FE space defined on 'mesh' used by the GridFunction 'nodes'. Uses
the FE collection from the global nodal FE space.
*/
std::unique_ptr<FiniteElementSpace> nodal_fes;
/**
'nodes': pointer to a GridFunction describing the physical location of the
MeshPart. Used for describing high-order and periodic meshes. This
GridFunction is defined on the FE space 'nodal_fes' which, in turn, is
defined on the Mesh 'mesh'.
*/
std::unique_ptr<GridFunction> nodes;
/** @name Connectivity to other MeshPart objects */
///@{
/// Total number of MeshParts
int num_parts;
/** @brief Index of the part described by this MeshPart:
0 <= 'my_part_id' < 'num_parts' */
int my_part_id;
/**
A group G is a subset 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 (by the partitioning array) 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 descriptions 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];
///@}
/** @brief Write the MeshPart to a stream using the parallel format
"MFEM mesh v1.2". */
void Print(std::ostream &os) const;
/** @brief Construct a serial Mesh object from the MeshPart.
The nodes of 'mesh' are NOT initialized by this method, however, the
nodal FE space and nodal GridFunction can be created and then attached to
the 'mesh'. The Mesh is constructed only if 'mesh' is empty, otherwise
the method simply returns the object held by 'mesh'.
*/
Mesh &GetMesh();
};
/** @brief Class that allows serial meshes to be partitioned into MeshPart
objects, typically one MeshPart at a time, which can then be used to write
the local mesh in parallel MFEM mesh format.
Sample usage of this class: partition a serial mesh and save it in parallel
MFEM format:
\code
// The array 'partitioning' can be obtained e.g. from
// mesh->GeneratePartitioning():
void usage1(Mesh *mesh, int num_parts, int *partitioning)
{
MeshPartitioner partitioner(*mesh, num_parts, partitioning);
MeshPart mesh_part;
for (int i = 0; i < num_parts; i++)
{
partitioner.ExtractPart(i, mesh_part);
ofstream omesh(MakeParFilename("my-mesh.", i));
mesh_part.Print(omesh);
}
}
\endcode
This class can also be used to partition a mesh and GridFunction(s) and save
them in parallel:
\code
// The array 'partitioning' can be obtained e.g. from
// mesh->GeneratePartitioning():
void usage2(Mesh *mesh, int num_parts, int *partitioning,
GridFunction *gf)
{
MeshPartitioner partitioner(*mesh, num_parts, partitioning);
MeshPart mesh_part;
for (int i = 0; i < num_parts; i++)
{
partitioner.ExtractPart(i, mesh_part);
ofstream omesh(MakeParFilename("my-mesh.", i));
mesh_part.Print(omesh);
auto lfes = partitioner.ExtractFESpace(mesh_part, *gf->FESpace());
auto lgf = partitioner.ExtractGridFunction(mesh_part, *gf, *lfes);
ofstream ofield(MakeParFilename("my-field.", i));
lgf->Save(ofield);
}
}
\endcode
*/
class MeshPartitioner
{
protected:
Mesh &mesh;
Array<int> partitioning;
Table part_to_element;
Table part_to_boundary;
Table edge_to_element;
Table vertex_to_element;
public:
/** @brief Construct a MeshPartitioner.
@param[in] mesh_ Mesh to be partitioned into MeshPart%s.
@param[in] num_parts_ Number of parts to partition the mesh into.
@param[in] partitioning_ Partitioning array: for every element in the
mesh gives the partition it belongs to; if NULL,
partitioning will be generated internally by
calling Mesh::GeneratePartitioning().
@param[in] part_method Partitioning method to be used in the call to
Mesh::GeneratePartitioning() when the provided
input partitioning is NULL.
*/
MeshPartitioner(Mesh &mesh_, int num_parts_, int *partitioning_ = NULL,
int part_method = 1);
/** @brief Construct a MeshPart corresponding to the given @a part_id.
@param[in] part_id Partition index to extract; valid values are in
the range [0, num_parts).
@param[out] mesh_part Output MeshPart object; its contents is
overwritten, while potentially reusing existing
dynamic memory allocations.
*/
void ExtractPart(int part_id, MeshPart &mesh_part) const;
/** @brief Construct a local version of the given FiniteElementSpace
@a global_fespace corresponding to the given @a mesh_part.
@param[in,out] mesh_part MeshPart on which to construct the local
FiniteElementSpace; this object is
generally modified by this call since it
calls mesh_part.GetMesh() to ensure the
local mesh is constructed.
@param[in] global_fespace The global FiniteElementSpace that should
be restricted to the @a mesh_part.
@returns A FiniteElementSpace pointer stored in a unique_ptr. The
returned local FiniteElementSpace is built on the Mesh object
contained in @a mesh_part (MeshPart::mesh) and it reuses the
FiniteElementCollection of the @a global_fespace.
*/
std::unique_ptr<FiniteElementSpace>
ExtractFESpace(MeshPart &mesh_part,
const FiniteElementSpace &global_fespace) const;
/** @brief Construct a local version of the given GridFunction, @a global_gf,
corresponding to the given @a mesh_part. The respective data is copied
from @a global_gf to the returned local GridFunction.
@param[in] mesh_part MeshPart on which to construct the local
GridFunction.
@param[in] global_gf The global GridFunction that should be
restricted to the @a mesh_part.
@param[in,out] local_fespace The local FiniteElementSpace corresponding
to @a mesh_part, e.g. constructed by the
method ExtractFESpace().
@returns A GridFunction pointer stored in a unique_ptr. The returned
local GridFunction is initialized with data appropriately copied
from @a global_gf.
*/
std::unique_ptr<GridFunction>
ExtractGridFunction(const MeshPart &mesh_part,
const GridFunction &global_gf,
FiniteElementSpace &local_fespace) const;
};
/** @brief Structure for storing mesh geometric factors: coordinates, Jacobians,
@@ -2454,7 +2788,6 @@ std::ostream &operator<<(std::ostream &out, const Mesh &mesh);
Mesh. See Mesh::GetGeometricFactors(). */
class GeometricFactors
{
private:
void Compute(const GridFunction &nodes,
MemoryType d_mt = MemoryType::DEFAULT);
@@ -2502,6 +2835,7 @@ public:
Vector detJ;
};
/** @brief Structure for storing face geometric factors: coordinates, Jacobians,
determinants of the Jacobians, and normal vectors. */
/** Typically objects of this type are constructed and owned by objects of class
@@ -2556,6 +2890,7 @@ public:
Vector normal;
};
/// Class used to extrude the nodes of a mesh
class NodeExtrudeCoefficient : public VectorCoefficient
{
@@ -2587,9 +2922,6 @@ inline void ShiftRight(int &a, int &b, int &c)
a = c; c = b; b = t;
}
/// @brief Print function for Mesh::FaceInformation.
std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info);
}
#endif
mesh/pmesh.cpp
+7 -8
View File
@@ -256,9 +256,6 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
// build svert_lvert mapping
BuildSharedVertMapping(nsvert, vert_element, vert_global_local);
delete vert_element;
SetMeshGen();
meshgen = mesh.meshgen; // copy the global 'meshgen'
}
if (mesh.NURBSext)
@@ -1527,6 +1524,7 @@ ParMesh ParMesh::MakeSimplicial(ParMesh &orig_mesh)
void ParMesh::Finalize(bool refine, bool fix_orientation)
{
const int meshgen_save = meshgen; // Mesh::Finalize() may call SetMeshGen()
// 'mesh_geoms' is local, so there's no need to save and restore it.
Mesh::Finalize(refine, fix_orientation);
@@ -4807,7 +4805,7 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
if (NURBSext)
{
Printer(os, comments); // does not print shared boundary
Printer(os, "", comments); // does not print shared boundary
return;
}
@@ -4935,7 +4933,7 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
if (set_names)
{
os << "mfem_mesh_end\n";
os << "\nmfem_mesh_end" << endl;
}
}
@@ -5286,7 +5284,7 @@ void ParMesh::PrintAsSerial(std::ostream &os, const std::string &comments) const
Mesh serialmesh = GetSerialMesh(save_rank);
if (MyRank == save_rank)
{
serialmesh.Printer(os, comments);
serialmesh.Printer(os, "", comments);
}
MPI_Barrier(MyComm);
}
@@ -6325,11 +6323,11 @@ void ParMesh::ParPrint(ostream &os, const std::string &comments) const
if (Nonconforming())
{
// the NC mesh format works both in serial and in parallel
Printer(os, comments);
Printer(os, "", comments);
return;
}
// Write out serial mesh. Tell serial mesh to deliniate the end of it's
// Write out serial mesh. Tell serial mesh to delineate the end of its
// output with 'mfem_serial_mesh_end' instead of 'mfem_mesh_end', as we will
// be adding additional parallel mesh information.
Printer(os, "mfem_serial_mesh_end", comments);
@@ -6346,6 +6344,7 @@ void ParMesh::ParPrint(ostream &os, const std::string &comments) const
{
os << "total_shared_faces " << sface_lface.Size() << '\n';
}
os << "\n# group 0 has no shared entities\n";
for (int gr = 1; gr < GetNGroups(); gr++)
{
{
mesh/tetrahedron.cpp
+4 -3
View File
@@ -53,7 +53,7 @@ void Tetrahedron::Init(int ind1, int ind2, int ind3, int ind4, int attr,
}
void Tetrahedron::ParseRefinementFlag(int refinement_edges[2], int &type,
int &flag)
int &flag) const
{
int i, f = refinement_flag;
@@ -134,9 +134,10 @@ void Tetrahedron::CreateRefinementFlag(int refinement_edges[2], int type,
refinement_flag |= refinement_edges[0];
}
void Tetrahedron::GetMarkedFace(const int face, int *fv)
void Tetrahedron::GetMarkedFace(const int face, int *fv) const
{
int re[2], type, flag, *tv = this->indices;
int re[2], type, flag;
const int *tv = this->indices;
ParseRefinementFlag(re, type, flag);
switch (face)
{
mesh/tetrahedron.hpp
+4 -3
View File
@@ -58,12 +58,13 @@ public:
/// Return element's type.
Type GetType() const override { return Element::TETRAHEDRON; }
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag);
void ParseRefinementFlag(int refinement_edges[2], int &type,
int &flag) const;
void CreateRefinementFlag(int refinement_edges[2], int type, int flag = 0);
void GetMarkedFace(const int face, int *fv);
void GetMarkedFace(const int face, int *fv) const;
int GetRefinementFlag() { return refinement_flag; }
int GetRefinementFlag() const { return refinement_flag; }
void SetRefinementFlag(int rf) { refinement_flag = rf; }
miniapps/meshing/makefile
+1 -1
View File
@@ -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
miniapps/meshing/mesh-explorer.cpp
+31 -2
View File
@@ -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"
"T) 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,33 @@ int main (int argc, char *argv[])
cout << "New mesh file: " << omesh_file << endl;
}
if (mk == 'T')
{
string mesh_prefix("mesh-explorer.mesh."), line;
MeshPartitioner partitioner(*mesh, np, partitioning);
MeshPart mesh_part;
cout << "Enter mesh file prefix or press <enter> to use \""
<< mesh_prefix << "\": " << flush;
// extract and ignore all characters after 'T' up to and including the
// new line:
cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
getline(cin, line);
if (!line.empty()) { mesh_prefix = line; }
int precision;
cout << "Enter floating point output precision (num. digits): "
<< 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";