Select Git revision
gmsh_io.cpp
gmsh_io.cpp 17.81 KiB
#include <iostream>
#include <sstream>
#include <cassert>
#include "gmsh_io.h"
#include "entity.h"
#ifdef USE_MPI
#include <mpi/mpi.h>
#endif /* USE_MPI */
std::string
quadrature_name(int order)
{
std::stringstream ss;
ss << "Gauss" << order;
return ss.str();
}
std::string
basis_func_name(int order)
{
(void) order;
std::stringstream ss;
//ss << "Lagrange";
ss << "H1Legendre" << order;
return ss.str();
}
std::string
basis_grad_name(int order)
{
(void) order;
std::stringstream ss;
//ss << "GradLagrange";
ss << "GradH1Legendre" << order;
return ss.str();
}
model::model()
: geometric_order(1), approximation_order(1)
#ifdef USE_MPI
, num_partitions(1)
#endif /* USE_MPI */
{}
model::model(int pgo, int pao)
: geometric_order(pgo), approximation_order(pao)
#ifdef USE_MPI
, num_partitions(1)
#endif /* USE_MPI */
{
assert(geometric_order > 0);
assert(approximation_order > 0);
}
model::~model()
{
//gmsh::clear();
}
void
model::make_boundary_to_faces_map(void)
{
#ifdef USE_MPI
ASSERT_MPI_RANK_0;
#endif /* USE_MPI */
gvp_t gmsh_entities;
/* Create a table mapping the tag of a boundary to the list of
* its face keys. This must happen before import_gmsh_entities()
* because otherwise you get spurious 2D entities (and this in
* turn because the constructor of entity calls addDiscreteEntity()
* to get all the element faces) .*/
gm::getEntities(gmsh_entities, DIMENSION(2));
for (auto [dim, tag] : gmsh_entities)
{
#ifdef USE_MPI
if (num_partitions > 1)
{
/* If the model is partitioned and entity has no parents,
* it means that there are no elements in the entity on which
* we will compute. Skip it. */
int pdim, ptag;
gmsh::model::getParent(dim, tag, pdim, ptag);
if (ptag == -1)
continue;
/* This stores which partitions share a boundary. It will be
* used to determine who communicates with who. */
std::vector<int> parts;
gm::getPartitions(dim, tag, parts);
comm_map[tag] = parts;
if (pdim == dim)
surface_to_parent_map[tag] = ptag;
}
#endif /* USE_MPI */
std::vector<int> elemTypes;
gmm::getElementTypes(elemTypes, dim, tag);
assert(elemTypes.size() == 1);
/* Get all the face keys for a given boundary. */
for (auto& elemType : elemTypes)
{
element_key_factory ekf(DIMENSION(2), tag, elemType);
boundary_map[tag] = ekf.get_keys();
}
}
#ifdef USE_MPI
//for (auto& sp : surface_to_parent_map)
// std::cout << sp.first << " -> " << sp.second << std::endl;
assert( IFF(comm_map.size() == 0, num_partitions == 1) );
#endif /* USE_MPI */
}
#ifdef USE_MPI
void
model::make_partition_to_entities_map()
{
ASSERT_MPI_RANK_0;
gvp_t gmsh_entities;
gm::getEntities(gmsh_entities, DIMENSION(3));
for (auto [dim, tag] : gmsh_entities)
{
if (num_partitions > 1)
{
/* If the model is partitioned and entity has no parents,
* it means that there are no elements in the entity on which
* we will compute. Skip it. */
int pdim, ptag;
gmsh::model::getParent(dim, tag, pdim, ptag);
if (ptag == -1)
continue;
std::vector<int> parts;
gm::getPartitions(dim, tag, parts);
assert(parts.size() == 1);
partition_map[ parts[0] ].push_back(tag);
partition_inverse_map[tag] = parts[0];
}
}
for (auto& [pn, tags] : partition_map)
std::sort(tags.begin(), tags.end()); /* For lookups */
}
#endif /* USE_MPI */
void
model::update_connectivity(const entity& e, size_t entity_index)
{
#ifdef USE_MPI
//ASSERT_MPI_RANK_0;
#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
for (size_t iF = 0; iF < e.num_faces_per_elem(); iF++)
{
size_t findex = e.num_faces_per_elem() * iT + iF;
neighbour_descriptor cd(entity_index, iT, iF);
element_key fk( e.face(findex) );
conn.connect( fk, cd );
}
}
}
#ifdef USE_MPI
int
model::my_partition(void) const
{
/* Here we are assuming that GMSH enumerates mesh partitions
* sequentially and starting from 1. */
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
return rank+1;
}
#endif /* USE_MPI */
void
model::import_gmsh_entities_rank0(void)
{
#ifdef USE_MPI
ASSERT_MPI_RANK_0;
#endif /* USE_MPI */
gvp_t gmsh_entities;
gm::getEntities(gmsh_entities, DIMENSION(3));
size_t entity_index = 0;
size_t dof_base = 0;
size_t flux_base = 0;
size_t index_base = 0;
#ifdef USE_MPI
size_t global_dof_base = 0;
#endif /* USE_MPI */
for (auto [dim, tag] : gmsh_entities)
{
#ifdef USE_MPI
if (num_partitions > 1)
{ /* If the model is partitioned and entity has no parents,
* it means that there are no elements in the entity on which
* we will compute. Skip it. */
int pdim, ptag;
gmsh::model::getParent(dim, tag, pdim, ptag);
if (ptag == -1)
continue;
}
#endif /* USE_MPI */
std::vector<int> elemTypes;
gmm::getElementTypes(elemTypes, dim, tag);
assert(elemTypes.size() == 1);
if (elemTypes.size() != 1)
throw std::invalid_argument("Only one element type per entity is allowed");
for (auto& elemType : elemTypes)
{
entity e({.dim = dim, .tag = tag, .etype = elemType,
.gorder = geometric_order, .aorder = approximation_order});
e.sort_by_orientation();
#ifdef USE_MPI
if (num_partitions > 1)
{
/* If there are multiple partitions, decide to which one this
* entity belongs to. If the partition is remote, send entity
* to the remote process. The entity is saved also locally
* to have all the data for postprocessing later. */
auto mp = my_partition();
auto tag_partition = partition_inverse_map.at(tag);
bool entity_is_remote = ( tag_partition != mp );
if (entity_is_remote)
{
dof_range dr(e, global_dof_base);
global_dof_ranges.emplace(tag, std::move(dr));
global_dof_base += e.num_dofs();
assert(tag_partition > 0);
int rank = tag_partition - 1;
e.mpi_send(rank, MPI_COMM_WORLD);
std::cout << "Sending entity " << tag << " to " << rank << " " << e.num_dofs() << std::endl;
/* Save locally the entity */
remote_entities[tag_partition].push_back( std::move(e) );
continue;
}
}
std::cout << "Keeping entity " << tag << std::endl;
#endif /* USE_MPI */
for (size_t i = 0; i < e.num_cells(); i++)
{
const auto& pe = e.cell(i);
auto eitor = etag_to_entity_offset.find( pe.element_tag() );
if ( eitor != etag_to_entity_offset.end() )
throw std::logic_error("Duplicate tag");
etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index, i);
}
e.base(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
flux_base += e.num_fluxes();
index_base += e.num_cells();
#ifdef USE_MPI
dof_range dr(e, global_dof_base); global_dof_ranges.emplace(tag, std::move(dr));
global_dof_base += e.num_dofs();
#endif /* USE_MPI */
update_connectivity(e, entity_index);
entities.push_back( std::move(e) );
entity_index++;
}
}
#ifdef USE_MPI
/* Compute DOF offset in the remote entities saved locally. */
assert( IMPLIES(num_partitions == 1, remote_entities.size() == 0) );
for (auto& re : remote_entities)
{
entity_index = 0;
dof_base = 0;
flux_base = 0;
index_base = 0;
for (auto& e : re.second)
{
e.base(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
flux_base += e.num_fluxes();
index_base += e.num_cells();
entity_index++;
}
}
#endif
}
#ifdef USE_MPI
void
model::import_gmsh_entities_rankN(void)
{
int mp = my_partition();
size_t entity_index = 0;
size_t dof_base = 0;
size_t flux_base = 0;
size_t index_base = 0;
for (auto& rtag : partition_map.at(mp))
{
entity e;
e.mpi_recv(0, MPI_COMM_WORLD);
std::cout << "Receiving " << rtag << ", rank " << mp-1 << " " << e.num_dofs() << std::endl;
for (size_t i = 0; i < e.num_cells(); i++)
{
const auto& pe = e.cell(i);
auto eitor = etag_to_entity_offset.find( pe.element_tag() );
if ( eitor != etag_to_entity_offset.end() )
throw std::logic_error("Duplicate tag");
etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index, i);
}
e.base(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
flux_base += e.num_fluxes();
index_base += e.num_cells();
update_connectivity(e, entity_index);
entities.push_back( std::move(e) );
entity_index++;
}}
bool
model::is_interprocess_boundary(int tag)
{
if (num_partitions > 1)
return comm_map.at(tag).size() == 2;
return false;
}
#endif /* USE_MPI */
void
model::map_boundaries(void)
{
#ifdef USE_MPI
//ASSERT_MPI_RANK_0;
#endif /* USE_MPI */
/* Make a vector mapping element_key to entity tag */
using bfk_t = std::pair<element_key, int>;
std::vector<bfk_t> bfk;
for (auto& [tag, keys] : boundary_map)
for (auto& k : keys)
bfk.push_back( std::make_pair(k, tag) );
/* Sort it for lookups */
auto comp = [](const bfk_t& a, const bfk_t& b) -> bool {
return a.first < b.first;
};
std::sort(bfk.begin(), bfk.end(), comp);
bnd_descriptors.resize( num_faces() );
size_t fbase = 0;
/* For each entity */
for (auto& e : entities)
{
/* and for each face */
for (size_t iF = 0; iF < e.num_faces(); iF++)
{
/* get the element key */
element_key fk(e.face(iF));
auto lbcomp = [](const bfk_t& a, const element_key& b) -> bool {
return a.first < b;
};
/* and lookup it in the vector we created before. */
auto itor = std::lower_bound(bfk.begin(), bfk.end(), fk, lbcomp);
if ( (itor == bfk.end()) or not ((*itor).first == fk) )
continue;
/* If found */
boundary_descriptor bd;
/* determine if it is an interface or a boundary */
if (conn.num_neighbours(fk) == 1)
{
#ifdef USE_MPI
if ( is_interprocess_boundary( (*itor).second ) )
bd.type = face_type::INTERPROCESS_BOUNDARY;
else
#endif /* USE_MPI */
bd.type = face_type::BOUNDARY;
}
else
{
bd.type = face_type::INTERFACE;
}
/* and store the gmsh tag in the descriptor. */
bd.gmsh_entity = (*itor).second;
#ifdef USE_MPI if (num_partitions > 1)
{
auto sitor = surface_to_parent_map.find( (*itor).second );
if ( sitor != surface_to_parent_map.end() )
bd.parent_entity = surface_to_parent_map.at( (*itor).second );
}
#endif /* USE_MPI */
/* Finally, put the descriptor in the global array of faces. */
bnd_descriptors.at(fbase + iF) = bd;
}
fbase += e.num_faces();
}
#if 0
size_t normal = 0;
size_t interface = 0;
size_t boundary = 0;
size_t ipc_boundary = 0;
for (auto& bd : bnd_descriptors)
{
switch (bd.type)
{
case face_type::NONE: normal++; break;
case face_type::INTERFACE: interface++; break;
case face_type::BOUNDARY: boundary++; break;
#ifdef USE_MPI
case face_type::INTERPROCESS_BOUNDARY: ipc_boundary++; break;
#endif /* USE_MPI */
}
}
std::cout << bfk.size() << " " << bnd_descriptors.size() << std::endl;
std::cout << normal << " " << interface << " " << boundary << " " << ipc_boundary << std::endl;
#endif
}
void
model::populate_from_gmsh_rank0()
{
#ifdef USE_MPI
ASSERT_MPI_RANK_0;
#endif /* USE_MPI */
gmm::setOrder( geometric_order );
make_boundary_to_faces_map();
#ifdef USE_MPI
make_partition_to_entities_map();
/* At this point, the following members are populated:
* - geometric_order, approximation_order and num_partitions
* - boundary_map
* - conn_map
* - partition_map
* We broadcast them to the other processes
*/
MPI_Bcast(&geometric_order, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&approximation_order, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&num_partitions, 1, MPI_INT, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(boundary_map, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(comm_map, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(partition_map, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(surface_to_parent_map, 0, MPI_COMM_WORLD);
/*
for (auto& [tag, vec] : partition_map)
{
std::cout << tag << ": " << vec.size() << " -> ";
for (auto& v : vec) std::cout << v << " ";
std::cout << std::endl;
}
*/
#endif /* USE_MPI */
entities.clear();
import_gmsh_entities_rank0();
map_boundaries(); /* This must happen after import_gmsh_entities(). */
}
#ifdef USE_MPI
void
model::populate_from_gmsh_rankN()
{
MPI_Bcast(&geometric_order, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&approximation_order, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&num_partitions, 1, MPI_INT, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(boundary_map, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(comm_map, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(partition_map, 0, MPI_COMM_WORLD);
priv_MPI_Bcast(surface_to_parent_map, 0, MPI_COMM_WORLD);
/*
for (auto& [tag, vec] : partition_map)
{
std::cout << tag << ": " << vec.size() << " -> ";
for (auto& v : vec)
std::cout << v << " ";
std::cout << std::endl;
}
*/
entities.clear();
import_gmsh_entities_rankN();
map_boundaries();
}
#endif /* USE_MPI */
void
model::populate_from_gmsh()
{
#ifdef USE_MPI
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if (rank == 0)
populate_from_gmsh_rank0();
else
populate_from_gmsh_rankN();
#else /* USE_MPI */
populate_from_gmsh_rank0();
#endif /* USE_MPI */
}
void
model::build()
{
generate_mesh();
make_boundary_to_faces_map();
populate_from_gmsh_rank0();
}
void
model::build(double sf)
{
generate_mesh(sf);
make_boundary_to_faces_map();
populate_from_gmsh_rank0();
}
void
model::generate_mesh()
{
gmm::generate( DIMENSION(3) );
}
void
model::generate_mesh(double sf)
{
gmm::generate( DIMENSION(3) );
std::vector<double> tr = {sf, 0, 0, 0,
0, sf, 0, 0,
0, 0, sf, 0 };
gmm::affineTransform(tr);
}
#ifdef USE_MPI
void
model::partition_mesh(int parts)
{
gmm::partition(parts);
num_partitions = parts;
}
entity&
model::entity_global_lookup(int tag)
{
/* Look for tag between local entities */
for (auto& e : entities)
if (e.gmsh_tag() == tag)
return e;
/* Try between remote entities */
int partition = partition_inverse_map.at(tag);
for (auto& e : remote_entities[partition])
if (e.gmsh_tag() == tag)
return e;
throw std::invalid_argument("tag not found");
}
const entity&
model::entity_global_lookup(int tag) const
{
/* Look for tag between local entities */
for (const auto& e : entities)
if (e.gmsh_tag() == tag)
return e;
/* Try between remote entities */
int partition = partition_inverse_map.at(tag);
for (const auto& e : remote_entities.at(partition))
if (e.gmsh_tag() == tag)
return e;
throw std::invalid_argument("tag not found");
}
#endif /* USE_MPI */
std::vector<entity>::const_iterator
model::begin() const
{
return entities.begin();
}
std::vector<entity>::const_iterator
model::end() const
{ return entities.end();
}
std::vector<entity>::iterator
model::begin()
{
return entities.begin();
}
std::vector<entity>::iterator
model::end()
{
return entities.end();
}
entity&
model::entity_at(size_t pos)
{
return entities.at(pos);
}
const entity&
model::entity_at(size_t pos) const
{
return entities.at(pos);
}
size_t
model::num_dofs(void) const
{
size_t ret = 0;
for (const auto& e : entities)
ret += e.num_dofs();
return ret;
}
size_t
model::num_fluxes(void) const
{
size_t ret = 0;
for (const auto& e : entities)
ret += e.num_fluxes();
return ret;
}
size_t
model::num_cells(void) const
{
size_t ret = 0;
for (const auto& e : entities)
ret += e.num_cells();
return ret;
}
size_t
model::num_faces(void) const
{
size_t ret = 0;
for (const auto& e : entities)
ret += e.num_faces();
return ret;
}
size_t
model::num_entities(void) const
{ return entities.size();
}
model::entofs_pair
model::lookup_tag(size_t tag) const
{
return etag_to_entity_offset.at(tag);
}