Select Git revision
RefFaceList.hpp
GmshReader.cpp 68.66 KiB
#include <GmshReader.hpp>
#include <PastisMacros.hpp>
#include <iostream>
#include <fstream>
#include <set>
#include <rang.hpp>
#include <CellType.hpp>
#include <Connectivity.hpp>
#include <Mesh.hpp>
#include <MeshData.hpp>
#include <RefFaceList.hpp>
#include <Messenger.hpp>
#include <Partitioner.hpp>
#include <ArrayUtils.hpp>
#include <unordered_map>
#include <map>
#include <regex>
#include <iomanip>
class ErrorHandler
{
public:
enum Type {
asked, /**< execution request by the user*/
compilation, /**< syntax error in a language */
normal, /**< normal error due to a bad use of ff3d */
unexpected /**< Unexpected execution error */
};
private:
const std::string __filename; /**< The source file name where the error occured */
const size_t __lineNumber; /**< The line number where exception was raised */
const std::string __errorMessage; /**< The reporting message */
const Type __type; /**< the type of the error */
public:
/**
* Prints the error message
*
*/
virtual void writeErrorMessage() const;
/**
* The copy constructor
*
* @param e an handled error
*/
ErrorHandler(const ErrorHandler& e)
: __filename(e.__filename),
__lineNumber(e.__lineNumber),
__errorMessage(e.__errorMessage),
__type(e.__type)
{
;
}
/**
* Constructor
*
* @param filename the source file name
* @param lineNumber the source line
* @param errorMessage the reported message
* @param type the type of the error
*/ ErrorHandler(const std::string& filename,
const size_t& lineNumber,
const std::string& errorMessage,
const Type& type);
/**
* The destructor
*
*/
virtual ~ErrorHandler()
{
;
}
};
void ErrorHandler::writeErrorMessage() const
{
switch(__type) {
case asked: {
perr() << "\nremark: exit command explicitly called\n";
[[fallthrough]];
}
case normal: {
perr() << '\n' << __filename << ':' << __lineNumber
<< ":remark: emitted the following message\n";
perr() << "error: " << __errorMessage << '\n';
break;
}
case compilation: {
perr() << "\nline " << __lineNumber << ':' << __errorMessage << '\n';
break;
}
case unexpected: {
perr() << '\n' << __filename << ':' << __lineNumber << ":\n" << __errorMessage << '\n';
perr() << "\nUNEXPECTED ERROR: this should not occure, please report it\n";
perr() << "\nBUG REPORT: Please send bug reports to:\n"
<< " ff3d-dev@nongnu.org or freefem@ann.jussieu.fr\n"
<< "or better, use the Bug Tracking System:\n"
<< " http://savannah.nongnu.org/bugs/?group=ff3d\n";
break;
}
default: {
perr() << __filename << ':' << __lineNumber << ": " << __errorMessage << '\n';
perr() << __FILE__ << ':' << __LINE__ << ":remark: error type not implemented!\n";
}
}
}
ErrorHandler::
ErrorHandler(const std::string& filename,
const size_t& lineNumber,
const std::string& errorMessage,
const Type& type)
: __filename(filename),
__lineNumber(lineNumber),
__errorMessage(errorMessage),
__type(type)
{
;
}
template <size_t Dimension>
class ConnectivityFace;
template<>
class ConnectivityFace<1>
{
public:
friend struct Hash;
struct Hash {
size_t operator()(const ConnectivityFace& f) const;
};
};
template<>
class ConnectivityFace<2>
{
public:
friend struct Hash;
struct Hash
{
size_t operator()(const ConnectivityFace& f) const {
size_t hash = 0;
hash ^= std::hash<unsigned int>()(f.m_node0_id);
hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
return hash;
}
};
unsigned int m_node0_id;
unsigned int m_node1_id;
friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
{
os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
return os;
}
PASTIS_INLINE
bool operator==(const ConnectivityFace& f) const
{
return ((m_node0_id == f.m_node0_id) and
(m_node1_id == f.m_node1_id));
}
PASTIS_INLINE
bool operator<(const ConnectivityFace& f) const
{
return ((m_node0_id<f.m_node0_id) or
((m_node0_id == f.m_node0_id) and
(m_node1_id<f.m_node1_id)));
}
PASTIS_INLINE
ConnectivityFace& operator=(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace& operator=(ConnectivityFace&&) = default;
PASTIS_INLINE
ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
{
Assert(given_node_id_list.size()==2);
#warning rework this dirty constructor
const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]);
m_node0_id = min;
m_node1_id = max;
}
PASTIS_INLINE
ConnectivityFace(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace(ConnectivityFace&&) = default;
PASTIS_INLINE
~ConnectivityFace() = default;
};
template <>
class ConnectivityFace<3>
{
private:
friend class GmshReader;
friend struct Hash;
struct Hash
{
size_t operator()(const ConnectivityFace& f) const {
size_t hash = 0;
for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
}
return hash;
}
};
bool m_reversed;
std::vector<NodeId::base_type> m_node_id_list;
friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
{
for (auto id : f.m_node_id_list) {
os << id << ' ';
}
return os;
}
PASTIS_INLINE
const bool& reversed() const
{
return m_reversed;
}
PASTIS_INLINE
const std::vector<unsigned int>& nodeIdList() const
{
return m_node_id_list;
}
PASTIS_INLINE
std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
{
const auto min_id = std::min_element(node_list.begin(), node_list.end());
const int shift = std::distance(node_list.begin(), min_id);
std::vector<unsigned int> rotated_node_list(node_list.size());
if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
for (size_t i=0; i<node_list.size(); ++i) {
rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
m_reversed = true;
}
} else {
for (size_t i=0; i<node_list.size(); ++i) {
rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
}
}
return rotated_node_list;
}
PASTIS_INLINE
ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
: m_reversed(false),
m_node_id_list(_sort(given_node_id_list))
{
;
}
public: bool operator==(const ConnectivityFace& f) const
{
if (m_node_id_list.size() == f.nodeIdList().size()) {
for (size_t j=0; j<m_node_id_list.size(); ++j) {
if (m_node_id_list[j] != f.nodeIdList()[j]) {
return false;
}
}
return true;
}
return false;
}
PASTIS_INLINE
bool operator<(const ConnectivityFace& f) const
{
const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
for (size_t i=0; i<min_nb_nodes; ++i) {
if (m_node_id_list[i] < f.m_node_id_list[i]) return true;
if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
}
return m_node_id_list.size() < f.m_node_id_list.size();
}
PASTIS_INLINE
ConnectivityFace& operator=(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace& operator=(ConnectivityFace&&) = default;
PASTIS_INLINE
ConnectivityFace(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace(ConnectivityFace&&) = default;
PASTIS_INLINE
ConnectivityFace() = delete;
PASTIS_INLINE
~ConnectivityFace() = default;
};
template <int Dimension>
class MeshDispatcher
{
public:
using MeshType = Mesh<Connectivity<Dimension>>;
private:
const MeshType& m_mesh;
CellValue<const int> m_cell_new_owner;
NodeValue<const int> m_node_new_owner;
using CellListToSendByProc = std::vector<Array<const CellId>>;
const CellListToSendByProc m_cell_list_to_send_by_proc;
Array<int> m_nb_cell_to_send_by_proc;
Array<int> m_nb_cell_to_recv_by_proc;
CellValue<int> _getCellNewOwner()
{
CSRGraph mesh_graph = m_mesh.cellToCellGraph();
Partitioner P;
CellValue<int> cell_new_owner(m_mesh.connectivity());
cell_new_owner = P.partition(mesh_graph);
return cell_new_owner;
}
NodeValue<int> _getNodeNewOwner()
{
const auto& node_to_cell_matrix
= m_mesh.connectivity().nodeToCellMatrix();
pout() << __FILE__ << ':' << __LINE__ << ": use Min function\n";
#warning could use a better policy
NodeValue<int> node_new_owner(m_mesh.connectivity());
parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
const auto& node_to_cell = node_to_cell_matrix[r];
CellId Jmin = node_to_cell[0];
for (size_t j=1; j<node_to_cell.size(); ++j) {
const CellId J = node_to_cell[j];
if (J<Jmin) {
Jmin=J;
}
}
node_new_owner[r] = m_cell_new_owner[Jmin];
});
#warning Add missing synchronize
return node_new_owner;
}
const CellListToSendByProc
_buildCellListToSend() const
{
const auto& node_to_cell_matrix
= m_mesh.connectivity().nodeToCellMatrix();
const auto& cell_to_node_matrix
= m_mesh.connectivity().cellToNodeMatrix();
std::vector<std::vector<CellId>> cell_vector_to_send_by_proc(parallel::size());
Array<bool> send_to_rank(parallel::size());
for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
send_to_rank.fill(false);
const auto& cell_to_node = cell_to_node_matrix[j];
for (size_t R=0; R<cell_to_node.size(); ++R) {
const NodeId& r = cell_to_node[R];
const auto& node_to_cell = node_to_cell_matrix[r];
for (size_t K=0; K<node_to_cell.size(); ++K) {
const CellId& k = node_to_cell[K];
send_to_rank[m_cell_new_owner[k]] = true;
}
}
for (size_t k=0; k<send_to_rank.size(); ++k) {
if (send_to_rank[k]) {
cell_vector_to_send_by_proc[k].push_back(j);
}
}
}
std::vector<Array<const CellId>> cell_list_to_send_by_proc(parallel::size());
for (size_t i=0; i<parallel::size(); ++i) {
cell_list_to_send_by_proc[i] = convert_to_array(cell_vector_to_send_by_proc[i]);
}
return cell_list_to_send_by_proc;
}
Array<int> _buildNbCellToSend()
{
Array<int> nb_cell_to_send_by_proc(parallel::size());
for (size_t i=0; i<parallel::size(); ++i) {
nb_cell_to_send_by_proc[i] = m_cell_list_to_send_by_proc[i].size();
} return nb_cell_to_send_by_proc;
}
public:
PASTIS_INLINE
const CellValue<const int>& cellNewOwner() const
{
return m_cell_new_owner;
}
PASTIS_INLINE
const NodeValue<const int>& nodeNewOwner() const
{
return m_node_new_owner;
}
template<typename DataType>
std::vector<Array<std::remove_const_t<DataType>>>
exchange(const CellValue<DataType>& cell_value) const
{
using MutableDataType = std::remove_const_t<DataType>;
std::vector<Array<DataType>> cell_value_to_send_by_proc(parallel::size());
for (size_t i=0; i<parallel::size(); ++i) {
const Array<const CellId>& cell_list = m_cell_list_to_send_by_proc[i];
Array<MutableDataType> cell_value_list(cell_list.size());
parallel_for (cell_list.size(), PASTIS_LAMBDA(const CellId& cell_id) {
cell_value_list[cell_id] = cell_value[cell_list[cell_id]];
});
cell_value_to_send_by_proc[i] = cell_value_list;
}
std::vector<Array<MutableDataType>> recv_cell_value_by_proc(parallel::size());
for (size_t i=0; i<parallel::size(); ++i) {
recv_cell_value_by_proc[i] = Array<MutableDataType>(m_nb_cell_to_recv_by_proc[i]);
}
parallel::exchange(cell_value_to_send_by_proc, recv_cell_value_by_proc);
return recv_cell_value_by_proc;
}
[[deprecated]]
const CellListToSendByProc& cell_list_to_send_by_proc() const
{
return m_cell_list_to_send_by_proc;
}
MeshDispatcher(const MeshType& mesh)
: m_mesh(mesh),
m_cell_new_owner(_getCellNewOwner()),
m_node_new_owner(_getNodeNewOwner()),
m_cell_list_to_send_by_proc(_buildCellListToSend()),
m_nb_cell_to_send_by_proc(_buildNbCellToSend()),
m_nb_cell_to_recv_by_proc(parallel::allToAll(m_nb_cell_to_send_by_proc))
{
;
}
MeshDispatcher(const MeshDispatcher&) = delete;
~MeshDispatcher() = default;
};
template <int Dimension>
void GmshReader::_dispatch()
{
using ConnectivityType = Connectivity<Dimension>;
using Rd = TinyVector<Dimension>;
using MeshType = Mesh<ConnectivityType>;
if (not m_mesh) {
ConnectivityDescriptor descriptor;
std::shared_ptr connectivity = std::make_shared<ConnectivityType>(descriptor);
NodeValue<Rd> xr;
m_mesh = std::make_shared<MeshType>(connectivity, xr);
}
const MeshType& mesh = static_cast<const MeshType&>(*m_mesh);
MeshDispatcher<Dimension> dispatcher(mesh);
std::vector<Array<int>> recv_cell_number_by_proc
= dispatcher.exchange(mesh.connectivity().cellNumber());
std::vector<Array<CellType>> recv_cell_type_by_proc
= dispatcher.exchange(mesh.connectivity().cellType());
std::vector<Array<int>> recv_cell_new_owner_by_proc
= dispatcher.exchange(dispatcher.cellNewOwner());
CellValue<int> number_of_node_per_cell(mesh.connectivity());
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
parallel_for(mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
number_of_node_per_cell[j] = cell_to_node_matrix[j].size();
});
std::vector<Array<int>> recv_number_of_node_per_cell_by_proc
= dispatcher.exchange(number_of_node_per_cell);
const auto& cell_list_to_send_by_proc = dispatcher.cell_list_to_send_by_proc();
const NodeValue<const int>& node_number = mesh.connectivity().nodeNumber();
std::vector<Array<int>> cell_node_number_to_send_by_proc(parallel::size());
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
std::vector<int> node_number_by_cell_vector;
for (size_t j=0; j<cell_list_to_send_by_proc[i_rank].size(); ++j) {
const CellId& cell_id = cell_list_to_send_by_proc[i_rank][j];
const auto& cell_node_list = cell_to_node_matrix[cell_id];
for (size_t r=0; r<cell_node_list.size(); ++r) {
const NodeId& node_id = cell_node_list[r];
node_number_by_cell_vector.push_back(node_number[node_id]);
}
}
cell_node_number_to_send_by_proc[i_rank] = convert_to_array(node_number_by_cell_vector);
}
std::vector<Array<int>> recv_cell_node_number_by_proc(parallel::size());
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
recv_cell_node_number_by_proc[i_rank]
= Array<int>(Sum(recv_number_of_node_per_cell_by_proc[i_rank]));
}
parallel::exchange(cell_node_number_to_send_by_proc, recv_cell_node_number_by_proc);
const std::unordered_map<int, int> node_number_id_map
= [&] () {
std::unordered_map<int, int> node_number_id_map;
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
int cpt=0;
for (size_t i=0; i<recv_cell_node_number_by_proc[i_rank].size(); ++i) {
int node_number = recv_cell_node_number_by_proc[i_rank][i];
auto [iterator, inserted] = node_number_id_map.insert(std::make_pair(node_number, cpt));
if (inserted) cpt++;
}
}
return node_number_id_map;
} ();
ConnectivityDescriptor new_descriptor;
new_descriptor.node_number_vector.resize(node_number_id_map.size());
for (const auto& [number, id] : node_number_id_map) {
new_descriptor.node_number_vector[id] = number;
}
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
int l=0;
for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
std::vector<unsigned int> node_vector;
for (int k=0; k<recv_number_of_node_per_cell_by_proc[i_rank][i]; ++k) {
const auto& searched_node_id = node_number_id_map.find(recv_cell_node_number_by_proc[i_rank][l++]);
Assert(searched_node_id != node_number_id_map.end());
node_vector.push_back(searched_node_id->second);
}
new_descriptor.cell_by_node_vector.emplace_back(node_vector);
}
}
const std::unordered_map<int, int> cell_number_id_map
= [&] () {
std::unordered_map<int, int> cell_number_id_map;
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
int cpt=0;
for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
int cell_number = recv_cell_number_by_proc[i_rank][i];
auto [iterator, inserted] = cell_number_id_map.insert(std::make_pair(cell_number, cpt));
if (inserted) cpt++;
}
}
return cell_number_id_map;
} ();
new_descriptor.cell_number_vector.resize(cell_number_id_map.size());
for (const auto& [number, id] : cell_number_id_map) {
new_descriptor.cell_number_vector[id] = number;
}
new_descriptor.cell_type_vector.resize(cell_number_id_map.size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
int cell_number = recv_cell_number_by_proc[i_rank][i];
const auto& searched_cell_id = cell_number_id_map.find(cell_number);
Assert(searched_cell_id != cell_number_id_map.end());
new_descriptor.cell_type_vector[searched_cell_id->second] = recv_cell_type_by_proc[i_rank][i];
}
}
std::vector<Array<const NodeId>> send_node_id_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
Array<bool> tag(mesh.numberOfNodes());
tag.fill(false);
std::vector<NodeId> node_id_vector;
for (size_t j=0; j<cell_list_to_send_by_proc[i_rank].size(); ++j) {
const CellId& cell_id = cell_list_to_send_by_proc[i_rank][j];
const auto& cell_node_list = cell_to_node_matrix[cell_id];
for (size_t r=0; r<cell_node_list.size(); ++r) {
const NodeId& node_id = cell_node_list[r];
if (not tag[node_id]) {
node_id_vector.push_back(node_id);
tag[node_id] = true;
}
}
}
send_node_id_by_proc[i_rank] = convert_to_array(node_id_vector);
}
std::vector<Array<const int>> send_node_number_by_proc(parallel::size()); for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
Array<int> send_node_number(send_node_id_by_proc[i_rank].size());
const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
parallel_for(send_node_number.size(), PASTIS_LAMBDA(const size_t& j){
send_node_number[j] = node_number[send_node_id[j]];
});
send_node_number_by_proc[i_rank] = send_node_number;
}
Array<unsigned int> nb_node_to_send_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
nb_node_to_send_by_proc[i_rank] = send_node_id_by_proc[i_rank].size();
}
Array<const unsigned int> nb_node_to_recv_by_proc
= parallel::allToAll(nb_node_to_send_by_proc);
std::vector<Array<int>> recv_node_number_by_proc(parallel::size());
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
recv_node_number_by_proc[i_rank] = Array<int>(nb_node_to_recv_by_proc[i_rank]);
}
parallel::exchange(send_node_number_by_proc, recv_node_number_by_proc);
std::vector<Array<const NodeId>> recv_node_id_correspondance_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<nb_node_to_recv_by_proc.size(); ++i_rank) {
Array<NodeId> node_id_correspondace(nb_node_to_recv_by_proc[i_rank]);
for (size_t l=0; l<nb_node_to_recv_by_proc[i_rank]; ++l) {
const int& node_number = recv_node_number_by_proc[i_rank][l];
const auto& searched_node_id = node_number_id_map.find(node_number);
Assert(searched_node_id != node_number_id_map.end());
node_id_correspondace[l] = searched_node_id->second;
}
recv_node_id_correspondance_by_proc[i_rank] = node_id_correspondace;
}
using ConnectivityType = Connectivity<Dimension>;
new_descriptor.cell_owner_vector.resize(cell_number_id_map.size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
int cell_number = recv_cell_number_by_proc[i_rank][i];
const auto& searched_cell_id = cell_number_id_map.find(cell_number);
Assert(searched_cell_id != cell_number_id_map.end());
new_descriptor.cell_owner_vector[searched_cell_id->second] = recv_cell_new_owner_by_proc[i_rank][i];
}
}
const NodeValue<const int>& new_node_owner = dispatcher.nodeNewOwner();
std::vector<Array<const int>> send_node_owner_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
Array<int> send_node_owner(nb_node_to_send_by_proc[i_rank]);
const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
parallel_for(send_node_id.size(), PASTIS_LAMBDA(const size_t& r) {
const NodeId& node_id = send_node_id[r];
send_node_owner[r] = new_node_owner[node_id];
});
send_node_owner_by_proc[i_rank] = send_node_owner;
}
std::vector<Array<int>> recv_node_owner_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
recv_node_owner_by_proc[i_rank] = Array<int>(nb_node_to_recv_by_proc[i_rank]);
}
parallel::exchange(send_node_owner_by_proc, recv_node_owner_by_proc);
new_descriptor.node_owner_vector.resize(new_descriptor.node_number_vector.size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
for (size_t r=0; r<recv_node_owner_by_proc[i_rank].size(); ++r) {
const NodeId& node_id = recv_node_id_correspondance_by_proc[i_rank][r];
new_descriptor.node_owner_vector[node_id] = recv_node_owner_by_proc[i_rank][r]; }
}
std::shared_ptr p_connectivity
= std::make_shared<ConnectivityType>(new_descriptor);
const NodeValue<const Rd>& xr = mesh.xr();
std::vector<Array<const Rd>> send_node_coord_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
Array<Rd> send_node_coord(nb_node_to_send_by_proc[i_rank]);
const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
parallel_for(send_node_id.size(), PASTIS_LAMBDA(const size_t& r) {
const NodeId& node_id = send_node_id[r];
send_node_coord[r] = xr[node_id];
});
send_node_coord_by_proc[i_rank] = send_node_coord;
}
std::vector<Array<Rd>> recv_node_coord_by_proc(parallel::size());
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
recv_node_coord_by_proc[i_rank] = Array<Rd>(nb_node_to_recv_by_proc[i_rank]);
}
parallel::exchange(send_node_coord_by_proc, recv_node_coord_by_proc);
NodeValue<Rd> new_xr(*p_connectivity);
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
for (size_t r=0; r<recv_node_coord_by_proc[i_rank].size(); ++r) {
const NodeId& node_id = recv_node_id_correspondance_by_proc[i_rank][r];
new_xr[node_id] = recv_node_coord_by_proc[i_rank][r];
}
}
m_mesh = std::make_shared<MeshType>(p_connectivity, new_xr);
}
GmshReader::GmshReader(const std::string& filename)
: m_filename(filename)
{
if (parallel::rank() == 0) {
try {
m_fin.open(m_filename);
if (not m_fin) {
perr() << "cannot read file '" << m_filename << "'\n";
std::exit(0);
}
// gmsh 2.2 format keywords
__keywordList["$MeshFormat"] = MESHFORMAT;
__keywordList["$EndMeshFormat"] = ENDMESHFORMAT;
__keywordList["$Nodes"] = NODES;
__keywordList["$EndNodes"] = ENDNODES;
__keywordList["$Elements"] = ELEMENTS;
__keywordList["$EndElements"] = ENDELEMENTS;
__keywordList["$PhysicalNames"] = PHYSICALNAMES;
__keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;
__numberOfPrimitiveNodes.resize(16);
__numberOfPrimitiveNodes[ 0] = 2; // edge
__numberOfPrimitiveNodes[ 1] = 3; // triangle
__numberOfPrimitiveNodes[ 2] = 4; // quadrangle
__numberOfPrimitiveNodes[ 3] = 4; // Tetrahedron
__numberOfPrimitiveNodes[ 4] = 8; // Hexaredron
__numberOfPrimitiveNodes[ 5] = 6; // Prism
__numberOfPrimitiveNodes[ 6] = 5; // Pyramid
__numberOfPrimitiveNodes[ 7] = 3; // second order edge
__numberOfPrimitiveNodes[ 8] = 6; // second order triangle
__numberOfPrimitiveNodes[ 9] = 9; // second order quadrangle
__numberOfPrimitiveNodes[10] = 10; // second order tetrahedron
__numberOfPrimitiveNodes[11] = 27; // second order hexahedron __numberOfPrimitiveNodes[12] = 18; // second order prism
__numberOfPrimitiveNodes[13] = 14; // second order pyramid
__numberOfPrimitiveNodes[14] = 1; // point
__primitivesNames[0] = "edges";
__supportedPrimitives[0] = true;
__primitivesNames[1] = "triangles";
__supportedPrimitives[1] = true;
__primitivesNames[2] = "quadrangles";
__supportedPrimitives[2] = true;
__primitivesNames[3] = "tetrahedra";
__supportedPrimitives[3] = true;
__primitivesNames[4] = "hexahedra";
__supportedPrimitives[4] = true;
__primitivesNames[5] = "prisms";
__supportedPrimitives[5] = false;
__primitivesNames[6] = "pyramids";
__supportedPrimitives[6] = false;
__primitivesNames[7] = "second order edges";
__supportedPrimitives[7] = false;
__primitivesNames[8] = "second order triangles";
__supportedPrimitives[8] = false;
__primitivesNames[9] = "second order quadrangles";
__supportedPrimitives[9] = false;
__primitivesNames[10] = "second order tetrahedra";
__supportedPrimitives[10]= false;
__primitivesNames[11] = "second order hexahedra";
__supportedPrimitives[11]= false;
__primitivesNames[12] = "second order prisms";
__supportedPrimitives[12]= false;
__primitivesNames[13] = "second order pyramids";
__supportedPrimitives[13]= false;
__primitivesNames[14] = "point";
__supportedPrimitives[14]= true;
pout() << "Reading file '" << m_filename << "'\n";
// Getting vertices list
GmshReader::Keyword kw = this->__nextKeyword();
switch(kw.second) {
// case NOD: {
// this->__readGmsh1();
// break;
// }
case MESHFORMAT: {
double fileVersion = this->_getReal();
if (fileVersion != 2.2) {
throw ErrorHandler(__FILE__,__LINE__,
"Cannot read Gmsh format '"+std::to_string(fileVersion)+"'",
ErrorHandler::normal);
}
int fileType = this->_getInteger();
__binary = (fileType == 1);
if ((fileType < 0) or (fileType > 1)) {
throw ErrorHandler(__FILE__,__LINE__,
"Cannot read Gmsh file type '"+std::to_string(fileType)+"'",
ErrorHandler::normal);
}
int dataSize = this->_getInteger();
if (dataSize != sizeof(double)) {
throw ErrorHandler(__FILE__,__LINE__,
"Data size not supported '"+std::to_string(dataSize)+"'",
ErrorHandler::normal);
}
if (__binary) {
// int one=0;
// fseek(__ifh,1L,SEEK_CUR); // fread(reinterpret_cast<char*>(&one),sizeof(int),1,__ifh);
// if (one == 1) {
// __convertEndian = false;
// } else {
// invertEndianess(one);
// if (one == 1) {
// __convertEndian = true;
// } else {
// throw ErrorHandler(__FILE__,__LINE__,
// "Cannot determine endianess",
// ErrorHandler::normal);
// }
// }
// pout() << "- Binary format: ";
// #ifdef WORDS_BIGENDIAN
// if (not __convertEndian) {
// pout() << "big endian\n";
// } else {
// pout() << "little endian\n";
// }
// #else // WORDS_BIGENDIAN
// if (not __convertEndian) {
// pout() << "little endian\n";
// } else {
// pout() << "big endian\n";
// }
// #endif // WORDS_BIGENDIAN
}
kw = this->__nextKeyword();
if (kw.second != ENDMESHFORMAT) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': expecting $EndMeshFormat, '"+kw.first+"' was found",
ErrorHandler::normal);
}
this->__readGmshFormat2_2();
break;
}
default: {
throw ErrorHandler(__FILE__,__LINE__,
"cannot determine format version of '"+m_filename+"'",
ErrorHandler::normal);
}
}
this->__proceedData();
// this->__createMesh();
}
catch(const ErrorHandler& e) {
e.writeErrorMessage();
std::exit(0);
}
}
pout() << std::flush;
if (parallel::size() > 1) {
pout() << "Sequential mesh read! Need to be dispatched\n" << std::flush;
const int mesh_dimension
= [&]() {
int mesh_dimension = -1; // unknown mesh dimension
if (m_mesh) {
mesh_dimension = m_mesh->meshDimension();
}
Array<int> dimensions = parallel::allGather(mesh_dimension);
std::set<int> dimension_set;
for (size_t i=0; i<dimensions.size(); ++i) {
const int i_dimension = dimensions[i];
if (i_dimension != -1) {
dimension_set.insert(i_dimension);
}
}
if (dimension_set.size() != 1) {
std::cerr << "error dimensions of read mesh parts differ!\n";
std::exit(1);
}
return *begin(dimension_set);
}();
switch (mesh_dimension) {
case 1: {
this->_dispatch<1>();
break;
}
case 2: {
this->_dispatch<2>();
break;
}
case 3: {
this->_dispatch<3>();
break;
}
default: {
perr() << "unexpected dimension " << mesh_dimension << '\n';
std::exit(1);
}
}
}
}
void GmshReader::__readVertices()
{
const int numberOfVerices = this->_getInteger();
pout() << "- Number Of Vertices: " << numberOfVerices << '\n';
if (numberOfVerices<0) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+this->m_filename
+"': number of vertices is negative",
ErrorHandler::normal);
}
__verticesNumbers.resize(numberOfVerices);
__vertices = Array<R3>(numberOfVerices);
if (not __binary) {
for (int i=0; i<numberOfVerices; ++i) {
__verticesNumbers[i] = this->_getInteger();
const double x = this->_getReal();
const double y = this->_getReal();
const double z = this->_getReal();
__vertices[i] = TinyVector<3, double>(x,y,z);
}
} else {
// fseek(m_fin.file()__ifh,1L,SEEK_CUR);
// for (int i=0; i<numberOfVerices; ++i) {
// int number = 0;
// fread(reinterpret_cast<char*>(&number),sizeof(int),1,__ifh);
// __verticesNumbers[i] = number;
// TinyVector<3,double> x;
// fread(reinterpret_cast<char*>(&(x[0])),sizeof(double),3,__ifh);
// for (size_t j=0; j<3; ++j) {
// __vertices[i][j] = x[j];
// }
// }
// if (__convertEndian) {
// for (int i=0; i<numberOfVerices; ++i) {
// invertEndianess(__verticesNumbers[i]);
// for (size_t j=0; j<3; ++j) {
// invertEndianess(vertices[i][j]);
// }
// }
// }
}
}
// void
// GmshReader::__readElements1()
// {
// const int numberOfElements = this->_getInteger();
// pout() << "- Number Of Elements: " << numberOfElements << '\n';
// if (numberOfElements<0) {
// throw ErrorHandler(__FILE__,__LINE__,
// "reading file '"+m_filename
// +"': number of elements is negative",
// ErrorHandler::normal);
// }
// __elementType.resize(numberOfElements);
// __references.resize(numberOfElements);
// __elementVertices.resize(numberOfElements);
// for (int i=0; i<numberOfElements; ++i) {
// this->_getInteger(); // drops element number
// const int elementType = this->_getInteger()-1;
// if ((elementType < 0) or (elementType > 14)) {
// throw ErrorHandler(__FILE__,__LINE__,
// "reading file '"+m_filename
// +"': unknown element type '"+std::to_string(elementType)+"'",
// ErrorHandler::normal);
// }
// __elementType[i] = elementType;
// __references[i] = this->_getInteger(); // physical reference
// this->_getInteger(); // drops entity reference
// const int numberOfVertices = this->_getInteger();
// if (numberOfVertices != __numberOfPrimitiveNodes[elementType]) {
// std::stringstream errorMsg;
// errorMsg << "reading file '" <<m_filename
// << "':number of vertices (" << numberOfVertices
// << ") does not match expected ("
// << __numberOfPrimitiveNodes[elementType] << ")" << std::ends;
// throw ErrorHandler(__FILE__,__LINE__,
// errorMsg.str(),
// ErrorHandler::normal);
// }
// __elementVertices[i].resize(numberOfVertices);
// for (int j=0; j<numberOfVertices; ++j) {
// __elementVertices[i][j] = this->_getInteger();
// }
// }
// }
void
GmshReader::__readElements2_2()
{
const int numberOfElements = this->_getInteger();
pout() << "- Number Of Elements: " << numberOfElements << '\n';
if (numberOfElements<0) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': number of elements is negative", ErrorHandler::normal);
}
__elementNumber.resize(numberOfElements);
__elementType.resize(numberOfElements);
__references.resize(numberOfElements);
__elementVertices.resize(numberOfElements);
if (not __binary) {
for (int i=0; i<numberOfElements; ++i) {
__elementNumber[i] = this->_getInteger();
const int elementType = this->_getInteger()-1;
if ((elementType < 0) or (elementType > 14)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': unknown element type '"+std::to_string(elementType)+"'",
ErrorHandler::normal);
}
__elementType[i] = elementType;
const int numberOfTags = this->_getInteger();
__references[i] = this->_getInteger(); // physical reference
for (int tag=1; tag<numberOfTags; ++tag) {
this->_getInteger(); // drops remaining tags
}
const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
__elementVertices[i].resize(numberOfVertices);
for (int j=0; j<numberOfVertices; ++j) {
__elementVertices[i][j] = this->_getInteger();
}
}
} else {
// fseek(__ifh,1L,SEEK_CUR);
// int i=0;
// for (;i<numberOfElements;) {
// int elementType = 0;
// fread(reinterpret_cast<char*>(&elementType),sizeof(int),1,__ifh);
// if (__convertEndian) {
// invertEndianess(elementType);
// }
// --elementType;
// if ((elementType < 0) or (elementType > 14)) {
// throw ErrorHandler(__FILE__,__LINE__,
// "reading file '"+m_filename
// +"': unknown element type '"+std::to_string(elementType)+"'",
// ErrorHandler::normal);
// }
// int elementTypeNumber = 0;
// fread(reinterpret_cast<char*>(&elementTypeNumber),sizeof(int),1,__ifh);
// if (__convertEndian) {
// invertEndianess(elementTypeNumber);
// }
// int numberOfTags = 0;
// fread(reinterpret_cast<char*>(&numberOfTags),sizeof(int),1,__ifh);
// if (__convertEndian) {
// invertEndianess(numberOfTags);
// }
// for (int k=0; k<elementTypeNumber; ++k) {
// __elementType[i] = elementType;
// int elementNumber = 0;
// fread(reinterpret_cast<char*>(&elementNumber),sizeof(int),1,__ifh);
// int reference = 0;
// fread(reinterpret_cast<char*>(&reference),sizeof(int),1,__ifh);
// __references[i] = reference; // physical reference
// for (int tag=1; tag<numberOfTags; ++tag) {
// int t;
// fread(reinterpret_cast<char*>(&t),sizeof(int),1,__ifh);
// }
// const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
// __elementVertices[i].resize(numberOfVertices);
// fread(reinterpret_cast<char*>(&(__elementVertices[i][0])),sizeof(int),numberOfVertices,__ifh);
// ++i;
// }
// }
// if (__convertEndian) {
// for (size_t i=0; i<__references.size(); ++i) {
// invertEndianess(__references[i]);
// }
// for (size_t i=0; i<__elementVertices.size(); ++i) {
// for (size_t j=0; j<__elementVertices[i].size(); ++i) {
// invertEndianess(__elementVertices[i][j]);
// }
// }
// }
}
}
void
GmshReader::
__readPhysicalNames2_2()
{
const int number_of_names = this->_getInteger();
for (int i=0; i<number_of_names; ++i) {
const int physical_dimension = this->_getInteger();
const int physical_number = this->_getInteger();
std::string physical_name;
m_fin >> physical_name;
physical_name = std::regex_replace(physical_name, std::regex("(\")"), "");
PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
auto searched_physical_ref_id = m_physical_ref_map.find(physical_number);
if (searched_physical_ref_id != m_physical_ref_map.end()) {
perr() << "Physical reference id '" << physical_ref_id
<< "' already defined as '" << searched_physical_ref_id->second << "'!";
std::exit(0);
}
m_physical_ref_map[physical_number] = physical_ref_id;
}
}
void
GmshReader::__compute2DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
{
// const auto& cell_to_node_matrix
// = this->_getMatrix(ItemType::cell, ItemType::node);
using Face = ConnectivityFace<2>;
// In 2D faces are simply define
using CellFaceId = std::pair<CellId, unsigned short>;
std::map<Face, std::vector<CellFaceId>> face_cells_map;
for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
const auto& cell_nodes = descriptor.cell_by_node_vector[j];
for (unsigned short r=0; r<cell_nodes.size(); ++r) {
NodeId node0_id = cell_nodes[r];
NodeId node1_id = cell_nodes[(r+1)%cell_nodes.size()];
if (node1_id<node0_id) {
std::swap(node0_id, node1_id);
}
face_cells_map[Face({node0_id, node1_id})].push_back(std::make_pair(j, r));
} }
{
descriptor.face_to_node_vector.resize(face_cells_map.size());
int l=0;
for (const auto& face_info : face_cells_map) {
const Face& face = face_info.first;
descriptor.face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
++l;
}
}
{
descriptor.face_to_cell_vector.resize(face_cells_map.size());
int l=0;
for (const auto& face_cells_vector : face_cells_map) {
const auto& [face, cell_info_vector] = face_cells_vector;
for (const auto& cell_info : cell_info_vector) {
descriptor.face_to_cell_vector[l].push_back(cell_info.first);
}
++l;
}
}
{
descriptor.face_number_vector.resize(face_cells_map.size());
for (size_t l=0; l<face_cells_map.size(); ++l) {
descriptor.face_number_vector[l] = l;
}
perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
}
}
void
GmshReader::__compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
{
using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
using Face = ConnectivityFace<3>;
// const auto& cell_to_node_matrix
// = this->_getMatrix(ItemType::cell, ItemType::node);
Array<unsigned short> cell_nb_faces(descriptor.cell_by_node_vector.size());
std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
const auto& cell_nodes = descriptor.cell_by_node_vector[j];
switch (descriptor.cell_type_vector[j]) {
case CellType::Tetrahedron: {
cell_nb_faces[j] = 4;
// face 0
Face f0({cell_nodes[1],
cell_nodes[2],
cell_nodes[3]});
face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
// face 1
Face f1({cell_nodes[0],
cell_nodes[3],
cell_nodes[2]});
face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
// face 2
Face f2({cell_nodes[0],
cell_nodes[1],
cell_nodes[3]});
face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
// face 3
Face f3({cell_nodes[0], cell_nodes[2],
cell_nodes[1]});
face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
break;
}
case CellType::Hexahedron: {
// face 0
Face f0({cell_nodes[3],
cell_nodes[2],
cell_nodes[1],
cell_nodes[0]});
face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
// face 1
Face f1({cell_nodes[4],
cell_nodes[5],
cell_nodes[6],
cell_nodes[7]});
face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
// face 2
Face f2({cell_nodes[0],
cell_nodes[4],
cell_nodes[7],
cell_nodes[3]});
face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
// face 3
Face f3({cell_nodes[1],
cell_nodes[2],
cell_nodes[6],
cell_nodes[5]});
face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
// face 4
Face f4({cell_nodes[0],
cell_nodes[1],
cell_nodes[5],
cell_nodes[4]});
face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
// face 5
Face f5({cell_nodes[3],
cell_nodes[7],
cell_nodes[6],
cell_nodes[2]});
face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
cell_nb_faces[j] = 6;
break;
}
default: {
perr() << "unexpected cell type!\n";
std::exit(0);
}
}
}
{
descriptor.cell_to_face_vector.resize(descriptor.cell_by_node_vector.size());
for (CellId j=0; j<descriptor.cell_to_face_vector.size(); ++j) {
descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
}
FaceId l=0;
for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
} ++l;
}
}
{
descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_by_node_vector.size());
for (CellId j=0; j<descriptor.cell_face_is_reversed_vector.size(); ++j) {
descriptor.cell_face_is_reversed_vector[j].resize(cell_nb_faces[j]);
}
for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
}
}
}
{
descriptor.face_to_node_vector.resize(face_cells_map.size());
int l=0;
for (const auto& face_info : face_cells_map) {
const Face& face = face_info.first;
descriptor.face_to_node_vector[l] = face.nodeIdList();
++l;
}
}
{
descriptor.face_number_vector.resize(face_cells_map.size());
for (size_t l=0; l<face_cells_map.size(); ++l) {
descriptor.face_number_vector[l] = l;
}
perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
}
}
void
GmshReader::__proceedData()
{
TinyVector<15,int> elementNumber = zero;
std::vector<std::set<size_t> > elementReferences(15);
for (size_t i=0; i<__elementType.size(); ++i) {
const int elementType = __elementType[i];
elementNumber[elementType]++;
elementReferences[elementType].insert(__references[i]);
}
for (size_t i=0; i<elementNumber.dimension(); ++i) {
if (elementNumber[i]>0) {
pout() << " - Number of "
<< __primitivesNames[i]
<< ": " << elementNumber[i];
if (not(__supportedPrimitives[i])) {
pout() << " [" << rang::fg::yellow << "IGNORED" << rang::style::reset << "]";
} else {
pout() << " references: ";
for(std::set<size_t>::const_iterator iref
= elementReferences[i].begin();
iref != elementReferences[i].end(); ++iref) {
if (iref != elementReferences[i].begin()) pout() << ',';
pout() << *iref;
}
switch (i) {
// Supported entities
case 0: { // edges
__edges = Array<Edge>(elementNumber[i]); __edges_ref.resize(elementNumber[i]);
__edges_number.resize(elementNumber[i]);
break;
}
case 1: { // triangles
__triangles = Array<Triangle>(elementNumber[i]);
__triangles_ref.resize(elementNumber[i]);
__triangles_number.resize(elementNumber[i]);
break;
}
case 2: { // quadrangles
__quadrangles = Array<Quadrangle>(elementNumber[i]);
__quadrangles_ref.resize(elementNumber[i]);
__quadrangles_number.resize(elementNumber[i]);
break;
}
case 3: { // tetrahedra
__tetrahedra = Array<Tetrahedron>(elementNumber[i]);
__tetrahedra_ref.resize(elementNumber[i]);
__tetrahedra_number.resize(elementNumber[i]);
break;
}
case 4: {// hexahedra
__hexahedra = Array<Hexahedron>(elementNumber[i]);
__hexahedra_ref.resize(elementNumber[i]);
__hexahedra_number.resize(elementNumber[i]);
break;
}
case 14: {// point
__points = Array<Point>(elementNumber[i]);
__points_ref.resize(elementNumber[i]);
__points_number.resize(elementNumber[i]);
break;
}
// Unsupported entities
case 5: // prism
case 6: // pyramid
case 7: // second order edge
case 8: // second order triangle
case 9: // second order quadrangle
case 10: // second order tetrahedron
case 11: // second order hexahedron
case 12: // second order prism
case 13: // second order pyramid
default: {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': ff3d cannot read this kind of element",
ErrorHandler::normal);
}
}
}
pout() << '\n';
}
}
pout() << "- Building correspondance table\n";
int minNumber = std::numeric_limits<int>::max();
int maxNumber = std::numeric_limits<int>::min();
for (size_t i=0; i<__verticesNumbers.size(); ++i) {
const int& vertexNumber = __verticesNumbers[i];
minNumber = std::min(minNumber,vertexNumber);
maxNumber = std::max(maxNumber,vertexNumber);
}
if (minNumber<0) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': ff3d does not implement negative vertices number",
ErrorHandler::normal); }
#warning should use an unordered_map
// A value of -1 means that the vertex is unknown
__verticesCorrepondance.resize(maxNumber+1,-1);
for (size_t i=0; i<__verticesNumbers.size(); ++i) {
__verticesCorrepondance[__verticesNumbers[i]] = i;
}
// reset element number to count them while filling structures
elementNumber = zero;
// Element structures filling
for (size_t i=0; i<__elementType.size(); ++i) {
std::vector<int>& elementVertices = __elementVertices[i];
switch (__elementType[i]) {
// Supported entities
case 0: { // edge
int& edgeNumber = elementNumber[0];
const int a = __verticesCorrepondance[elementVertices[0]];
const int b = __verticesCorrepondance[elementVertices[1]];
if ((a<0) or (b<0)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': error reading element "+std::to_string(i)
+" [bad vertices definition]",
ErrorHandler::normal);
}
__edges[edgeNumber]
= Edge(a, b);
__edges_ref[edgeNumber] = __references[i];
edgeNumber++; // one more edge
break;
}
case 1: { // triangles
int& triangleNumber = elementNumber[1];
const int a = __verticesCorrepondance[elementVertices[0]];
const int b = __verticesCorrepondance[elementVertices[1]];
const int c = __verticesCorrepondance[elementVertices[2]];
if ((a<0) or (b<0) or (c<0)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': error reading element "+std::to_string(i)
+" [bad vertices definition]",
ErrorHandler::normal);
}
__triangles[triangleNumber]
= Triangle(a, b, c);
__triangles_ref[triangleNumber] = __references[i];
__triangles_number[triangleNumber] = __elementNumber[i];
triangleNumber++; // one more triangle
break;
}
case 2: { // quadrangle
int& quadrilateralNumber = elementNumber[2];
const int a = __verticesCorrepondance[elementVertices[0]];
const int b = __verticesCorrepondance[elementVertices[1]];
const int c = __verticesCorrepondance[elementVertices[2]];
const int d = __verticesCorrepondance[elementVertices[3]];
if ((a<0) or (b<0) or (c<0) or (d<0)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': error reading element "+std::to_string(i)
+" [bad vertices definition]",
ErrorHandler::normal);
}
__quadrangles[quadrilateralNumber] = Quadrangle(a,b,c,d); __quadrangles_ref[quadrilateralNumber] = __references[i];
__quadrangles_number[quadrilateralNumber] = __elementNumber[i];
quadrilateralNumber++; // one more quadrangle
break;
}
case 3: { // tetrahedra
int& tetrahedronNumber = elementNumber[3];
const int a = __verticesCorrepondance[elementVertices[0]];
const int b = __verticesCorrepondance[elementVertices[1]];
const int c = __verticesCorrepondance[elementVertices[2]];
const int d = __verticesCorrepondance[elementVertices[3]];
if ((a<0) or (b<0) or (c<0) or (d<0)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': error reading element "+std::to_string(i)
+" [bad vertices definition]",
ErrorHandler::normal);
}
__tetrahedra[tetrahedronNumber] = Tetrahedron(a, b, c, d);
__tetrahedra_ref[tetrahedronNumber] = __references[i];
__tetrahedra_number[tetrahedronNumber] = __elementNumber[i];
tetrahedronNumber++; // one more tetrahedron
break;
}
case 4: { // hexaredron
int& hexahedronNumber = elementNumber[4];
const int a = __verticesCorrepondance[elementVertices[0]];
const int b = __verticesCorrepondance[elementVertices[1]];
const int c = __verticesCorrepondance[elementVertices[2]];
const int d = __verticesCorrepondance[elementVertices[3]];
const int e = __verticesCorrepondance[elementVertices[4]];
const int f = __verticesCorrepondance[elementVertices[5]];
const int g = __verticesCorrepondance[elementVertices[6]];
const int h = __verticesCorrepondance[elementVertices[7]];
if ((a<0) or (b<0) or (c<0) or (d<0) or (e<0) or (f<0) or (g<0) or (h<0)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': error reading element "+std::to_string(i)
+" [bad vertices definition]",
ErrorHandler::normal);
}
__hexahedra[hexahedronNumber]
= Hexahedron(a, b, c, d,
e, f, g, h);
__hexahedra_ref[hexahedronNumber] = __references[i];
__hexahedra_number[hexahedronNumber] = __elementNumber[i];
hexahedronNumber++; // one more hexahedron
break;
}
// Unsupported entities
case 14:{// point
int& point_number = elementNumber[14];
const int a = __verticesCorrepondance[elementVertices[0]];
__points[point_number] = a;
__points_ref[point_number] = __references[i];
__points_number[point_number] = __elementNumber[i];
point_number++;
break;
}
case 5: // prism
case 6: // pyramid
case 7: // second order edge
case 8: // second order triangle
case 9: // second order quadrangle
case 10: // second order tetrahedron
case 11: // second order hexahedron
case 12: // second order prism
case 13:{// second order pyramid }
default: {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': ff3d cannot read this kind of element",
ErrorHandler::normal);
}
}
}
TinyVector<15,int> dimension0_mask = zero;
dimension0_mask[14]=1;
TinyVector<15,int> dimension1_mask = zero;
dimension1_mask[0]=1;
dimension1_mask[7]=1;
TinyVector<15,int> dimension2_mask = zero;
dimension2_mask[1]=1;
dimension2_mask[2]=1;
dimension2_mask[8]=1;
dimension2_mask[9]=1;
TinyVector<15,int> dimension3_mask = zero;
dimension3_mask[3]=1;
dimension3_mask[4]=1;
dimension3_mask[5]=1;
dimension3_mask[6]=1;
dimension3_mask[10]=1;
dimension3_mask[11]=1;
dimension3_mask[12]=1;
dimension3_mask[13]=1;
pout() << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
pout() << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
pout() << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
pout() << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
if ((dimension3_mask, elementNumber)>0) {
const size_t nb_cells = (dimension3_mask, elementNumber);
ConnectivityDescriptor descriptor;
descriptor.node_number_vector = __verticesNumbers;
descriptor.cell_type_vector.resize(nb_cells);
descriptor.cell_number_vector.resize(nb_cells);
descriptor.cell_by_node_vector.resize(nb_cells);
const size_t nb_tetrahedra = __tetrahedra.size();
for (size_t j=0; j<nb_tetrahedra; ++j) {
descriptor.cell_by_node_vector[j].resize(4);
for (int r=0; r<4; ++r) {
descriptor.cell_by_node_vector[j][r] = __tetrahedra[j][r];
}
descriptor.cell_type_vector[j] = CellType::Tetrahedron;
descriptor.cell_number_vector[j] = __tetrahedra_number[j];
}
const size_t nb_hexahedra = __hexahedra.size();
for (size_t j=0; j<nb_hexahedra; ++j) {
const size_t jh = nb_tetrahedra+j;
descriptor.cell_by_node_vector[jh].resize(8);
for (int r=0; r<8; ++r) {
descriptor.cell_by_node_vector[jh][r] = __hexahedra[j][r];
}
descriptor.cell_type_vector[jh] = CellType::Hexahedron;
descriptor.cell_number_vector[jh] = __hexahedra_number[j];
}
__compute3DCellFaceAndFaceNodeConnectivities(descriptor);
descriptor.cell_owner_vector.resize(nb_cells);
std::fill(descriptor.cell_owner_vector.begin(),
descriptor.cell_owner_vector.end(),
parallel::rank());
descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
std::fill(descriptor.face_owner_vector.begin(),
descriptor.face_owner_vector.end(),
parallel::rank());
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(),
descriptor.node_owner_vector.end(),
parallel::rank());
std::shared_ptr p_connectivity = std::make_shared<Connectivity3D>(descriptor);
Connectivity3D& connectivity = *p_connectivity;
using Face = ConnectivityFace<3>;
const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
= [&] {
std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
const auto& face_node_list = connectivity.faceToNodeMatrix();
for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
const auto& node_list = face_node_list[l];
std::vector<unsigned int> node_vector(node_list.size());
for (size_t r=0; r<node_list.size(); ++r) {
node_vector[r] = node_list[r];
}
face_to_id_map[Face(node_vector)] = l;
}
return face_to_id_map;
} ();
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int f=0; f<__triangles.size(); ++f) {
const unsigned int face_number
= [&]{
auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]}));
if (i == face_to_id_map.end()) {
std::cerr << "face not found!\n";
std::terminate();
}
return i->second;
}();
const unsigned int& ref = __triangles_ref[f];
ref_faces_map[ref].push_back(face_number);
}
for (unsigned int f=0; f<__quadrangles.size(); ++f) {
const unsigned int face_number
= [&]{
auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
__quadrangles[f][2], __quadrangles[f][3]}));
if (i == face_to_id_map.end()) {
std::cerr << "face not found!\n";
std::terminate();
}
return i->second;
}();
const unsigned int& ref = __quadrangles_ref[f];
ref_faces_map[ref].push_back(face_number);
}
for (const auto& ref_face_list : ref_faces_map) {
Array<FaceId> face_list(ref_face_list.second.size());
for (size_t j=0; j<ref_face_list.second.size(); ++j) {
face_list[j]=ref_face_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
}
using MeshType = Mesh<Connectivity3D>; using Rd = TinyVector<3, double>;
NodeValue<Rd> xr(connectivity);
for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
xr[i][2] = __vertices[i][2];
}
m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
} else if ((dimension2_mask, elementNumber)>0) {
const size_t nb_cells = (dimension2_mask, elementNumber);
ConnectivityDescriptor descriptor;
descriptor.node_number_vector = __verticesNumbers;
descriptor.cell_type_vector.resize(nb_cells);
descriptor.cell_number_vector.resize(nb_cells);
descriptor.cell_by_node_vector.resize(nb_cells);
const size_t nb_triangles = __triangles.size();
for (size_t j=0; j<nb_triangles; ++j) {
descriptor.cell_by_node_vector[j].resize(3);
for (int r=0; r<3; ++r) {
descriptor.cell_by_node_vector[j][r] = __triangles[j][r];
}
descriptor.cell_type_vector[j] = CellType::Triangle;
descriptor.cell_number_vector[j] = __triangles_number[j];
}
const size_t nb_quadrangles = __quadrangles.size();
for (size_t j=0; j<nb_quadrangles; ++j) {
const size_t jq = j+nb_triangles;
descriptor.cell_by_node_vector[jq].resize(4);
for (int r=0; r<4; ++r) {
descriptor.cell_by_node_vector[jq][r] = __quadrangles[j][r];
}
descriptor.cell_type_vector[jq] = CellType::Quadrangle;
descriptor.cell_number_vector[jq] = __quadrangles_number[j];
}
this->__compute2DCellFaceAndFaceNodeConnectivities(descriptor);
descriptor.cell_owner_vector.resize(nb_cells);
std::fill(descriptor.cell_owner_vector.begin(),
descriptor.cell_owner_vector.end(),
parallel::rank());
descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
std::fill(descriptor.face_owner_vector.begin(),
descriptor.face_owner_vector.end(),
parallel::rank());
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(),
descriptor.node_owner_vector.end(),
parallel::rank());
std::shared_ptr p_connectivity = std::make_shared<Connectivity2D>(descriptor);
Connectivity2D& connectivity = *p_connectivity;
using Face = ConnectivityFace<2>;
const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
= [&] {
std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
const auto& face_node_list = connectivity.faceToNodeMatrix();
for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
const auto& node_list = face_node_list[l];
face_to_id_map[Face({node_list[0], node_list[1]})] = l;
} return face_to_id_map;
} ();
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int e=0; e<__edges.size(); ++e) {
const unsigned int edge_number
= [&]{
auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}));
if (i == face_to_id_map.end()) {
std::cerr << "face " << __edges[e][0] << " not found!\n";
std::terminate();
}
return i->second;
}();
const unsigned int& ref = __edges_ref[e];
ref_faces_map[ref].push_back(edge_number);
}
for (const auto& ref_face_list : ref_faces_map) {
Array<FaceId> face_list(ref_face_list.second.size());
for (size_t j=0; j<ref_face_list.second.size(); ++j) {
face_list[j]=ref_face_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
}
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
for (unsigned int r=0; r<__points.size(); ++r) {
const unsigned int point_number = __points[r];
const unsigned int& ref = __points_ref[r];
ref_points_map[ref].push_back(point_number);
}
for (const auto& ref_point_list : ref_points_map) {
Array<NodeId> point_list(ref_point_list.second.size());
for (size_t j=0; j<ref_point_list.second.size(); ++j) {
point_list[j]=ref_point_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
}
using MeshType = Mesh<Connectivity2D>;
using Rd = TinyVector<2, double>;
NodeValue<Rd> xr(connectivity);
for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
}
m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
} else if ((dimension1_mask, elementNumber)>0) {
const size_t nb_cells = (dimension1_mask, elementNumber);
ConnectivityDescriptor descriptor;
descriptor.node_number_vector = __verticesNumbers;
descriptor.cell_type_vector.resize(nb_cells);
descriptor.cell_number_vector.resize(nb_cells);
descriptor.cell_by_node_vector.resize(nb_cells);
for (size_t j=0; j<nb_cells; ++j) {
descriptor.cell_by_node_vector[j].resize(2);
for (int r=0; r<2; ++r) {
descriptor.cell_by_node_vector[j][r] = __edges[j][r];
}
descriptor.cell_type_vector[j] = CellType::Line;
descriptor.cell_number_vector[j] = __edges_number[j]; }
descriptor.cell_owner_vector.resize(nb_cells);
std::fill(descriptor.cell_owner_vector.begin(),
descriptor.cell_owner_vector.end(),
parallel::rank());
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(),
descriptor.node_owner_vector.end(),
parallel::rank());
std::shared_ptr p_connectivity = std::make_shared<Connectivity1D>(descriptor);
Connectivity1D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
for (unsigned int r=0; r<__points.size(); ++r) {
const unsigned int point_number = __points[r];
const unsigned int& ref = __points_ref[r];
ref_points_map[ref].push_back(point_number);
}
for (const auto& ref_point_list : ref_points_map) {
Array<NodeId> point_list(ref_point_list.second.size());
for (size_t j=0; j<ref_point_list.second.size(); ++j) {
point_list[j]=ref_point_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
}
using MeshType = Mesh<Connectivity1D>;
using Rd = TinyVector<1, double>;
NodeValue<Rd> xr(connectivity);
for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
}
m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
} else {
perr() << "*** using a dimension 0 mesh is forbidden!\n";
std::exit(0);
}
}
GmshReader::Keyword
GmshReader::__nextKeyword()
{
GmshReader::Keyword kw;
std::string aKeyword;
m_fin >> aKeyword;
if (not m_fin) {
kw.second = EndOfFile;
return kw;
}
KeywordList::iterator i = __keywordList.find(aKeyword.c_str());
if (i != __keywordList.end()) {
kw.first = (*i).first;
kw.second = (*i).second;
return kw;
}
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename +"': unknown keyword '"+aKeyword+"'",
ErrorHandler::normal);
kw.first = aKeyword;
kw.second = Unknown;
return kw;
}
void GmshReader::
__readGmshFormat2_2()
{
pout() << "- Reading Gmsh format 2.2\n";
GmshReader::Keyword kw = std::make_pair("", Unknown);
while (kw.second != EndOfFile) {
kw = this->__nextKeyword();
switch (kw.second) {
case NODES: {
this->__readVertices();
if (this->__nextKeyword().second != ENDNODES) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': expecting $EndNodes, '"+kw.first+"' was found",
ErrorHandler::normal);
}
break;
}
case ELEMENTS: {
this->__readElements2_2();
kw = this->__nextKeyword();
if (kw.second != ENDELEMENTS) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': expecting $EndElements, '"+kw.first+"' was found",
ErrorHandler::normal);
}
break;
}
case PHYSICALNAMES: {
this->__readPhysicalNames2_2();
if (this->__nextKeyword().second != ENDPHYSICALNAMES) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': expecting $EndNodes, '"+kw.first+"' was found",
ErrorHandler::normal);
}
break;
}
case EndOfFile: {
break;
}
default: {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': unexpected '"+kw.first+"'",
ErrorHandler::normal);
}
}
}
}