Select Git revision
ExecutionStatManager.cpp
-
Stéphane Del Pino authored
First sequential implementation
Stéphane Del Pino authoredFirst sequential implementation
MeshBuilderBase.cpp 8.17 KiB
#include <mesh/MeshBuilderBase.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
#include <mesh/ConnectivityDispatcher.hpp>
#include <mesh/ConnectivityDispatcherVariant.hpp>
#include <mesh/ItemId.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshVariant.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <vector>
template <size_t Dimension>
void
MeshBuilderBase::_dispatch()
{
if (parallel::size() == 1) {
return;
}
using ConnectivityType = Connectivity<Dimension>;
using Rd = TinyVector<Dimension>;
using MeshType = Mesh<Dimension>;
if (not m_mesh) {
ConnectivityDescriptor descriptor;
std::shared_ptr connectivity = ConnectivityType::build(descriptor);
NodeValue<Rd> xr;
m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(connectivity, xr));
}
const MeshType& mesh = *(m_mesh->get<const MeshType>());
auto p_dispatcher = std::make_shared<const ConnectivityDispatcher<Dimension>>(mesh.connectivity());
m_connectivity_dispatcher = std::make_shared<ConnectivityDispatcherVariant>(p_dispatcher);
std::shared_ptr dispatched_connectivity = p_dispatcher->dispatchedConnectivity();
NodeValue<Rd> dispatched_xr = p_dispatcher->dispatch(mesh.xr());
m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(dispatched_connectivity, dispatched_xr));
}
template void MeshBuilderBase::_dispatch<1>();
template void MeshBuilderBase::_dispatch<2>();
template void MeshBuilderBase::_dispatch<3>();
template <size_t Dimension>
void
MeshBuilderBase::_checkMesh() const
{
using MeshType = Mesh<Dimension>;
using ConnectivityType = typename MeshType::Connectivity;
if (not m_mesh) {
throw UnexpectedError("mesh is not built yet");
}
const MeshType& mesh = *(m_mesh->get<const MeshType>());
const ConnectivityType& connectivity = mesh.connectivity();
if constexpr (Dimension > 2) { // check for duplicated edges
auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
std::vector<std::vector<EdgeId>> node_edges(mesh.numberOfNodes());
for (EdgeId edge_id = 0; edge_id < mesh.numberOfEdges(); ++edge_id) {
auto node_list = edge_to_node_matrix[edge_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
node_edges[node_list[i_node]].push_back(edge_id);
}
}
for (auto&& edge_list : node_edges) {
std::sort(edge_list.begin(), edge_list.end());
}
for (EdgeId edge_id = 0; edge_id < mesh.numberOfEdges(); ++edge_id) {
auto node_list = edge_to_node_matrix[edge_id];
std::vector<EdgeId> intersection = node_edges[node_list[0]];
for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
std::vector<EdgeId> local_intersection;
std::set_intersection(intersection.begin(), intersection.end(), //
node_edges[node_list[i_node]].begin(), node_edges[node_list[i_node]].end(),
std::back_inserter(local_intersection));
std::swap(local_intersection, intersection);
if (intersection.size() < 2) {
break;
}
}
if (intersection.size() > 1) {
std::ostringstream error_msg;
error_msg << "invalid mesh.\n\tFollowing edges\n";
for (EdgeId duplicated_edge_id : intersection) {
error_msg << "\t - id=" << duplicated_edge_id << " number=" << connectivity.edgeNumber()[duplicated_edge_id]
<< '\n';
}
error_msg << "\tare duplicated";
throw NormalError(error_msg.str());
}
}
}
if constexpr (Dimension > 1) { // check for duplicated faces
auto face_to_node_matrix = connectivity.faceToNodeMatrix();
std::vector<std::vector<FaceId>> node_faces(mesh.numberOfNodes());
for (FaceId face_id = 0; face_id < mesh.numberOfFaces(); ++face_id) {
auto node_list = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
node_faces[node_list[i_node]].push_back(face_id);
}
}
for (auto&& face_list : node_faces) {
std::sort(face_list.begin(), face_list.end());
}
for (FaceId face_id = 0; face_id < mesh.numberOfFaces(); ++face_id) {
auto node_list = face_to_node_matrix[face_id];
std::vector<FaceId> intersection = node_faces[node_list[0]];
for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
std::vector<FaceId> local_intersection;
std::set_intersection(intersection.begin(), intersection.end(), //
node_faces[node_list[i_node]].begin(), node_faces[node_list[i_node]].end(),
std::back_inserter(local_intersection));
std::swap(local_intersection, intersection);
if (intersection.size() < 2) {
break;
}
}
if (intersection.size() > 1) {
std::ostringstream error_msg;
error_msg << "invalid mesh.\n\tFollowing faces\n";
for (FaceId intersection_face_id : intersection) {
error_msg << "\t - id=" << intersection_face_id
<< " number=" << connectivity.faceNumber()[intersection_face_id] << '\n';
error_msg << "\t nodes:";
for (size_t i = 0; i < face_to_node_matrix[intersection_face_id].size(); ++i) {
error_msg << ' ' << face_to_node_matrix[intersection_face_id][i];
}
error_msg << '\n';
}
error_msg << "\tare duplicated";
throw NormalError(error_msg.str());
}
}
}
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
{ // check for duplicated cells
std::vector<std::vector<CellId>> node_cells(mesh.numberOfNodes());
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto node_list = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
node_cells[node_list[i_node]].push_back(cell_id);
}
}
for (auto&& cell_list : node_cells) {
std::sort(cell_list.begin(), cell_list.end());
}
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto node_list = cell_to_node_matrix[cell_id];
std::vector<CellId> intersection = node_cells[node_list[0]];
for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
std::vector<CellId> local_intersection;
std::set_intersection(intersection.begin(), intersection.end(), //
node_cells[node_list[i_node]].begin(), node_cells[node_list[i_node]].end(),
std::back_inserter(local_intersection));
std::swap(local_intersection, intersection);
if (intersection.size() < 2) {
break;
}
}
if (intersection.size() > 1) {
std::ostringstream error_msg;
error_msg << "invalid mesh.\n\tFollowing cells\n";
for (CellId duplicated_cell_id : intersection) {
error_msg << "\t - id=" << duplicated_cell_id << " number=" << connectivity.cellNumber()[duplicated_cell_id]
<< '\n';
}
error_msg << "\tare duplicated";
throw NormalError(error_msg.str());
}
}
}
const auto& Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
const auto& xr = mesh.xr();
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
double cell_volume = 0;
auto cell_nodes = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
cell_volume += dot(Cjr(cell_id, i_node), xr[cell_nodes[i_node]]);
}
if (cell_volume <= 0) {
std::ostringstream error_msg;
error_msg << "invalid mesh.\n\tThe following cell\n";
error_msg << "\t - id=" << cell_id << " number=" << connectivity.cellNumber()[cell_id] << '\n';
error_msg << "\thas non-positive volume: " << cell_volume / Dimension;
throw NormalError(error_msg.str());
}
}
}
template void MeshBuilderBase::_checkMesh<1>() const;
template void MeshBuilderBase::_checkMesh<2>() const;
template void MeshBuilderBase::_checkMesh<3>() const;