diff --git a/CMakeLists.txt b/CMakeLists.txt
index d7c7318458823c05c9f385c1b5c50118534f6a76..93eb3fa005eda92b4884d61cfa8963fc237b5354 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -706,6 +706,12 @@ else()
message(" gnuplot: not found!")
endif()
+if (GMSH)
+ message(" gmsh: ${GMSH}")
+else()
+ message(" gmsh: not found!")
+endif()
+
if (PYGMENTIZE)
message(" pygmentize: ${PYGMENTIZE}")
else()
@@ -718,7 +724,7 @@ else()
message(" pdflatex: not found!")
endif()
-if (NOT EMACS OR NOT GNUPLOT_FOUND)
+if (NOT EMACS OR NOT GNUPLOT_FOUND OR NOT GMSH)
message(" ** Cannot build documentation: missing ")
elseif(NOT LATEX_PDFLATEX_FOUND OR NOT PYGMENTIZE)
message(" ** Cannot build pdf documentation: missing")
@@ -729,6 +735,9 @@ endif()
if (NOT GNUPLOT_FOUND)
message(" - gnuplot")
endif()
+if (NOT GMSH)
+ message(" - gmsh")
+endif()
if (NOT LATEX_PDFLATEX_FOUND)
message(" - pdflatex")
endif()
diff --git a/README.md b/README.md
index 62d6b43773cb0164a29e8a3716a818c5ec219c7d..36e460a851a53d40b4ab442460594091fa22b20f 100644
--- a/README.md
+++ b/README.md
@@ -86,7 +86,7 @@ apt install slepc-dev
### User documentation
-To build documentation one requires `emacs` and `gnuplot`,
+To build documentation one requires `emacs`, `gmsh` and `gnuplot`,
additionally since examples results are generated, the documentation
can only be produced after the compilation of `pugs` itself.
@@ -94,6 +94,10 @@ To install `emacs` on Debian-like systems
```shell
apt install emacs
```
+To install `gmsh` on Debian-like systems
+```shell
+apt install gmsh
+```
To install `gnuplot` one can either use
```shell
apt install gnuplot-nox
@@ -106,6 +110,7 @@ apt install gnuplot-x11
> When building the documentation for the first time, a local `emacs`
> configuration is generated. *This requires an internet connection.*
+
These packages are enough to build the html documentation. To build
the pdf documentation one requires a few more packages: `pdflatex`
(actually a fresh texlive installation is probably necessary) and `pygmentize`
diff --git a/cmake/PugsDoc.cmake b/cmake/PugsDoc.cmake
index 1f8701ccc20f6416b529908f3893d971e86c7452..b067a5cc828ec0c245e8d51c0bda46a96dcce702 100644
--- a/cmake/PugsDoc.cmake
+++ b/cmake/PugsDoc.cmake
@@ -12,10 +12,13 @@ find_program(PYGMENTIZE pygmentize)
# check for gnuplot
find_package(Gnuplot)
+# check for gmsh
+find_program(GMSH NAMES gmsh)
+
add_custom_target(userdoc)
add_custom_target(doc DEPENDS userdoc)
-if (EMACS AND GNUPLOT_FOUND)
+if (EMACS AND GNUPLOT_FOUND AND GMSH)
add_custom_command(
OUTPUT "${PUGS_BINARY_DIR}/doc"
@@ -129,6 +132,13 @@ else()
COMMAND ${CMAKE_COMMAND} -E cmake_echo_color --red --bold "gnuplot missing")
add_dependencies(userdoc userdoc-missing-gnuplot)
endif()
+
+ if (NOT GMSH)
+ add_custom_target(userdoc-missing-gmsh
+ COMMAND ${CMAKE_COMMAND} -E cmake_echo_color --no-newline "Cannot build documentation: "
+ COMMAND ${CMAKE_COMMAND} -E cmake_echo_color --red --bold "gmsh missing")
+ add_dependencies(userdoc userdoc-missing-gmsh)
+ endif()
endif()
add_dependencies(doc userdoc)
diff --git a/src/language/utils/OFStream.cpp b/src/language/utils/OFStream.cpp
index 39b7cc8ecedee3ff76e3d6d54790d25f90d01569..ca7c77fad8a69fd79ae2b512635acbb7f03277fd 100644
--- a/src/language/utils/OFStream.cpp
+++ b/src/language/utils/OFStream.cpp
@@ -1,10 +1,12 @@
#include <language/utils/OFStream.hpp>
+#include <utils/Filesystem.hpp>
#include <utils/Messenger.hpp>
OFStream::OFStream(const std::string& filename)
{
if (parallel::rank() == 0) {
+ createDirectoryIfNeeded(filename);
m_fstream.open(filename);
if (m_fstream.is_open()) {
m_ostream = &m_fstream;
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index b7de42715eff88510119b54507ed32323365c817..459a827c84ead76ed23acdf28cbd8525923694da 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -87,57 +87,108 @@ MeshFlatNodeBoundary<2>::_getNormal(const Mesh<Connectivity<2>>& mesh)
template <>
TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
+MeshFlatNodeBoundary<3>::_getFarestNode(const Mesh<Connectivity<3>>& mesh, const Rd& x0, const Rd& x1)
{
- using R3 = TinyVector<3, double>;
-
- std::array<R3, 6> bounds = this->_getBounds(mesh);
+ const NodeValue<const Rd>& xr = mesh.xr();
+ const auto node_number = mesh.connectivity().nodeNumber();
- const R3& xmin = bounds[0];
- const R3& ymin = bounds[1];
- const R3& zmin = bounds[2];
- const R3& xmax = bounds[3];
- const R3& ymax = bounds[4];
- const R3& zmax = bounds[5];
+ using NodeNumberType = std::remove_const_t<typename decltype(node_number)::data_type>;
- const R3 u = xmax - xmin;
- const R3 v = ymax - ymin;
- const R3 w = zmax - zmin;
+ Rd t = x1 - x0;
+ t *= 1. / l2Norm(t);
- const R3 uv = crossProduct(u, v);
- const double uv_l2 = dot(uv, uv);
+ double farest_distance = 0;
+ Rd farest_x = zero;
+ NodeNumberType farest_number = std::numeric_limits<NodeNumberType>::max();
- R3 normal = uv;
- double normal_l2 = uv_l2;
+ auto node_list = this->m_ref_node_list.list();
- const R3 uw = crossProduct(u, w);
- const double uw_l2 = dot(uw, uw);
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId& node_id = node_list[i_node];
+ const Rd& x = xr[node_id];
+ const double distance = l2Norm(crossProduct(t, x - x0));
- if (uw_l2 > uv_l2) {
- normal = uw;
- normal_l2 = uw_l2;
+ if ((distance > farest_distance) or ((distance == farest_distance) and (node_number[node_id] < farest_number))) {
+ farest_distance = distance;
+ farest_number = node_number[node_id];
+ farest_x = x;
+ }
}
- const R3 vw = crossProduct(v, w);
- const double vw_l2 = dot(vw, vw);
+ if (parallel::size()) {
+ Array<double> farest_distance_array = parallel::allGather(farest_distance);
+ Array<Rd> farest_x_array = parallel::allGather(farest_x);
+ Array<NodeNumberType> farest_number_array = parallel::allGather(farest_number);
- if (vw_l2 > normal_l2) {
- normal = vw;
- normal_l2 = vw_l2;
+ Assert(farest_distance_array.size() == farest_x_array.size());
+ Assert(farest_distance_array.size() == farest_number_array.size());
+
+ for (size_t i = 0; i < farest_distance_array.size(); ++i) {
+ if ((farest_distance_array[i] > farest_distance) or
+ ((farest_distance_array[i] == farest_distance) and (farest_number_array[i] < farest_number))) {
+ farest_distance = farest_distance_array[i];
+ farest_number = farest_number_array[i];
+ farest_x = farest_x_array[i];
+ }
+ }
}
- if (normal_l2 == 0) {
+ return farest_x;
+}
+
+template <>
+TinyVector<3, double>
+MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
+{
+ using R3 = TinyVector<3, double>;
+
+ std::array<R3, 2> diagonal = [](const std::array<R3, 6>& bounds) {
+ size_t max_i = 0;
+ size_t max_j = 0;
+ double max_length = 0;
+
+ for (size_t i = 0; i < bounds.size(); ++i) {
+ for (size_t j = i + 1; j < bounds.size(); ++j) {
+ double length = l2Norm(bounds[i] - bounds[j]);
+ if (length > max_length) {
+ max_i = i;
+ max_j = j;
+ max_length = length;
+ }
+ }
+ }
+
+ return std::array<R3, 2>{bounds[max_i], bounds[max_j]};
+ }(this->_getBounds(mesh));
+
+ const R3& x0 = diagonal[0];
+ const R3& x1 = diagonal[1];
+
+ if (x0 == x1) {
std::ostringstream ost;
ost << "invalid boundary \"" << rang::fgB::yellow << m_ref_node_list.refId() << rang::style::reset
<< "\": unable to compute normal";
throw NormalError(ost.str());
}
- const double length = sqrt(normal_l2);
+ const R3 x2 = this->_getFarestNode(mesh, x0, x1);
+
+ const R3 u = x1 - x0;
+ const R3 v = x2 - x0;
+
+ R3 normal = crossProduct(u, v);
+ const double normal_norm = l2Norm(normal);
+
+ if (normal_norm == 0) {
+ std::ostringstream ost;
+ ost << "invalid boundary \"" << rang::fgB::yellow << m_ref_node_list.refId() << rang::style::reset
+ << "\": unable to compute normal";
+ throw NormalError(ost.str());
+ }
- normal *= 1. / length;
+ normal *= (1. / normal_norm);
- this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
+ this->_checkBoundaryIsFlat(normal, 1. / 3. * (x0 + x1 + x2), normal_norm, mesh);
return normal;
}
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
index 8486f02c7377fe2987560d29f2ee7cc45c1b44cf..54a10ae1c1412005440cd9e6cc49973c017265e0 100644
--- a/src/mesh/MeshFlatNodeBoundary.hpp
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -13,21 +13,26 @@ class [[nodiscard]] MeshFlatNodeBoundary final
private:
const Rd m_outgoing_normal;
+ Rd _getFarestNode(const Mesh<Connectivity<Dimension>>& mesh, const Rd& x0, const Rd& x1);
+
Rd _getNormal(const Mesh<Connectivity<Dimension>>& mesh);
- void _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal, const TinyVector<Dimension, double>& origin,
- const double length, const Mesh<Connectivity<Dimension>>& mesh) const;
+ void _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+ const TinyVector<Dimension, double>& origin,
+ const double length,
+ const Mesh<Connectivity<Dimension>>& mesh) const;
Rd _getOutgoingNormal(const Mesh<Connectivity<Dimension>>& mesh);
public:
- const Rd& outgoingNormal() const
+ const Rd&
+ outgoingNormal() const
{
return m_outgoing_normal;
}
MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
- MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
+ MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
template <size_t MeshDimension>
friend MeshFlatNodeBoundary<MeshDimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
@@ -47,7 +52,7 @@ class [[nodiscard]] MeshFlatNodeBoundary final
public:
MeshFlatNodeBoundary() = default;
MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
- MeshFlatNodeBoundary(MeshFlatNodeBoundary &&) = default;
+ MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default;
~MeshFlatNodeBoundary() = default;
};
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 50969799eeb4250b28d21722e8476ab86a88d9eb..8ccf7767180b30524ba54f7c4f085c3ad35783d2 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -5,6 +5,7 @@
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
+#include <utils/Filesystem.hpp>
#include <utils/Messenger.hpp>
#include <utils/PugsTraits.hpp>
#include <utils/RevisionInfo.hpp>
@@ -243,8 +244,16 @@ GnuplotWriter::_write(const MeshType& mesh,
const OutputNamedItemDataSet& output_named_item_data_set,
std::optional<double> time) const
{
+ createDirectoryIfNeeded(_getFilename());
+
if (parallel::rank() == 0) {
std::ofstream fout{_getFilename()};
+ if (not fout) {
+ std::ostringstream error_msg;
+ error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilename() << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+
fout.precision(15);
fout.setf(std::ios_base::scientific);
@@ -286,6 +295,12 @@ GnuplotWriter::_write(const MeshType& mesh,
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
if (i_rank == parallel::rank()) {
std::ofstream fout(_getFilename(), std::ios_base::app);
+ if (not fout) {
+ std::ostringstream error_msg;
+ error_msg << "cannot open file \"" << rang::fgB::yellow << _getFilename() << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+
fout.precision(15);
fout.setf(std::ios_base::scientific);
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index a3156904843bdeb77be718de6a10614a3333afae..3e0740446e961505eb01e36eee5ba168e0f47dcc 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -5,6 +5,7 @@
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
+#include <utils/Filesystem.hpp>
#include <utils/Messenger.hpp>
#include <utils/PugsTraits.hpp>
#include <utils/RevisionInfo.hpp>
@@ -322,6 +323,12 @@ GnuplotWriter1D::_write(const MeshType& mesh,
if (parallel::rank() == 0) {
fout.open(_getFilename());
+ if (not fout) {
+ std::ostringstream error_msg;
+ error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilename() << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+
fout.precision(15);
fout.setf(std::ios_base::scientific);
fout << _getDateAndVersionComment();
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 0a40c7265e05a17d31e59bd41eeb10cb63935821..6ec6295e44b54bef9ea21c110120fdd873f71488 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -4,6 +4,7 @@
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
+#include <utils/Filesystem.hpp>
#include <utils/Messenger.hpp>
#include <utils/RevisionInfo.hpp>
#include <utils/Stringify.hpp>
@@ -364,8 +365,16 @@ VTKWriter::_write(const MeshType& mesh,
output_named_item_data_set.add(NamedItemData{"cell_number", mesh.connectivity().cellNumber()});
output_named_item_data_set.add(NamedItemData{"node_number", mesh.connectivity().nodeNumber()});
+ createDirectoryIfNeeded(m_base_filename);
+
if (parallel::rank() == 0) { // write PVTK file
std::ofstream fout(_getFilenamePVTU());
+ if (not fout) {
+ std::ostringstream error_msg;
+ error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilenamePVTU() << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+
fout << "<?xml version=\"1.0\"?>\n";
fout << _getDateAndVersionComment();
fout << "<VTKFile type=\"PUnstructuredGrid\">\n";
@@ -424,7 +433,7 @@ VTKWriter::_write(const MeshType& mesh,
fout << "</PCellData>\n";
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
- fout << "<Piece Source=\"" << _getFilenameVTU(i_rank) << "\"/>\n";
+ fout << "<Piece Source=" << std::filesystem::path{_getFilenameVTU(i_rank)}.filename() << "/>\n";
}
fout << "</PUnstructuredGrid>\n";
fout << "</VTKFile>\n";
@@ -435,6 +444,13 @@ VTKWriter::_write(const MeshType& mesh,
// write VTK files
std::ofstream fout(_getFilenameVTU(parallel::rank()));
+ if (not fout) {
+ std::ostringstream error_msg;
+ error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilenameVTU(parallel::rank()) << rang::fg::reset
+ << '"';
+ throw NormalError(error_msg.str());
+ }
+
fout << "<?xml version=\"1.0\"?>\n";
fout << _getDateAndVersionComment();
fout << "<VTKFile type=\"UnstructuredGrid\"";
@@ -632,7 +648,13 @@ VTKWriter::_write(const MeshType& mesh,
}
if (parallel::rank() == 0) { // write PVD file
- std::ofstream fout(m_base_filename + ".pvd");
+ const std::string pvd_filename = m_base_filename + ".pvd";
+ std::ofstream fout(pvd_filename);
+ if (not fout) {
+ std::ostringstream error_msg;
+ error_msg << "cannot create file \"" << rang::fgB::yellow << pvd_filename << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
fout << "<?xml version=\"1.0\"?>\n";
fout << _getDateAndVersionComment();
@@ -645,11 +667,13 @@ VTKWriter::_write(const MeshType& mesh,
sout << m_base_filename;
sout << '.' << std::setfill('0') << std::setw(4) << i_time << ".pvtu";
- fout << "<DataSet timestep=\"" << m_period_manager->savedTime(i_time) << "\" file=\"" << sout.str() << "\"/>\n";
+ fout << "<DataSet timestep=\"" << m_period_manager->savedTime(i_time)
+ << "\" file=" << std::filesystem::path{sout.str()}.filename() << "/>\n";
}
- fout << "<DataSet timestep=\"" << *time << "\" file=\"" << _getFilenamePVTU() << "\"/>\n";
+ fout << "<DataSet timestep=\"" << *time << "\" file=" << std::filesystem::path{_getFilenamePVTU()}.filename()
+ << "/>\n";
} else {
- fout << "<DataSet file=\"" << _getFilenamePVTU() << "\"/>\n";
+ fout << "<DataSet file=" << std::filesystem::path{_getFilenamePVTU()}.filename() << "/>\n";
}
fout << "</Collection>\n";
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 264ecef042956e6b77230e00f49cd321b50dcc6a..754c6e3ac6448a96e5e05c0a7f60967243b4798f 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -1,5 +1,8 @@
#include <scheme/FluxingAdvectionSolver.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/PrismTransformation.hpp>
#include <language/utils/EvaluateAtPoints.hpp>
#include <language/utils/InterpolateItemArray.hpp>
#include <language/utils/InterpolateItemValue.hpp>
@@ -322,6 +325,10 @@ template <>
FaceValue<double>
FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
{
+ // due to the z component of the jacobian determinant, degree 3
+ // polynomials must be exactly integrated
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(3));
+
const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
auto old_xr = m_old_mesh->xr();
@@ -329,31 +336,10 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
parallel_for(
m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
const auto& face_to_node = face_to_node_matrix[face_id];
- if (face_to_node.size() == 4) {
- const Rd& x0 = old_xr[face_to_node[0]];
- const Rd& x1 = old_xr[face_to_node[1]];
- const Rd& x2 = old_xr[face_to_node[2]];
- const Rd& x3 = old_xr[face_to_node[3]];
-
- const Rd& x4 = new_xr[face_to_node[0]];
- const Rd& x5 = new_xr[face_to_node[1]];
- const Rd& x6 = new_xr[face_to_node[2]];
- const Rd& x7 = new_xr[face_to_node[3]];
-
- const Rd& a1 = x6 - x1;
- const Rd& a2 = x6 - x3;
- const Rd& a3 = x6 - x4;
-
- const Rd& b1 = x7 - x0;
- const Rd& b2 = x5 - x0;
- const Rd& b3 = x2 - x0;
- TinyMatrix<3> M1(a1 + b1, a2, a3);
- TinyMatrix<3> M2(b1, a2 + b2, a3);
- TinyMatrix<3> M3(a1, b2, a3 + b3);
+ const size_t face_nb_node = face_to_node.size();
- algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
- } else if (face_to_node.size() == 3) {
+ if (face_nb_node == 3) {
const Rd& x0 = old_xr[face_to_node[0]];
const Rd& x1 = old_xr[face_to_node[1]];
const Rd& x2 = old_xr[face_to_node[2]];
@@ -362,21 +348,40 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
const Rd& x4 = new_xr[face_to_node[1]];
const Rd& x5 = new_xr[face_to_node[2]];
- const Rd& a1 = x5 - x1;
- const Rd& a2 = x5 - x2;
- const Rd& a3 = x5 - x3;
+ double volume = 0;
- const Rd& b1 = x5 - x0;
- const Rd& b2 = x4 - x0;
- const Rd& b3 = x2 - x0;
-
- TinyMatrix<3> M1(a1 + b1, a2, a3);
- TinyMatrix<3> M2(b1, a2 + b2, a3);
- TinyMatrix<3> M3(a1, b2, a3 + b3);
+ PrismTransformation T(x0, x1, x2, x3, x4, x5);
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ volume += gauss.weight(i) * T.jacobianDeterminant(gauss.point(i));
+ }
- algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
+ algebraic_fluxing_volume[face_id] = volume;
} else {
- throw NotImplementedError("Not implemented for non quad faces");
+ Rd xg_old = zero;
+ for (size_t i_node = 0; i_node < face_nb_node; ++i_node) {
+ xg_old += old_xr[face_to_node[i_node]];
+ }
+ xg_old *= 1. / face_nb_node;
+ Rd xg_new = zero;
+ for (size_t i_node = 0; i_node < face_nb_node; ++i_node) {
+ xg_new += new_xr[face_to_node[i_node]];
+ }
+ xg_new *= 1. / face_nb_node;
+
+ double volume = 0;
+ for (size_t i_node = 0; i_node < face_nb_node; ++i_node) {
+ PrismTransformation T(xg_old, //
+ old_xr[face_to_node[i_node]], //
+ old_xr[face_to_node[(i_node + 1) % face_nb_node]], //
+ xg_new, //
+ new_xr[face_to_node[i_node]], //
+ new_xr[face_to_node[(i_node + 1) % face_nb_node]]);
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ volume += gauss.weight(i) * T.jacobianDeterminant(gauss.point(i));
+ }
+ }
+
+ algebraic_fluxing_volume[face_id] = volume;
}
});
diff --git a/src/utils/Filesystem.hpp b/src/utils/Filesystem.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..63af52b2e432f1baf55d2cbfd895485a408f7e10
--- /dev/null
+++ b/src/utils/Filesystem.hpp
@@ -0,0 +1,23 @@
+#ifndef FILESYSTEM_HPP
+#define FILESYSTEM_HPP
+
+#include <utils/Exceptions.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <filesystem>
+
+PUGS_INLINE void
+createDirectoryIfNeeded(const std::string& filename)
+{
+ std::filesystem::path path = std::filesystem::path{filename}.parent_path();
+ if (not path.empty()) {
+ try {
+ std::filesystem::create_directories(path);
+ }
+ catch (std::filesystem::filesystem_error& e) {
+ throw NormalError(e.what());
+ }
+ }
+}
+
+#endif // FILESYSTEM_HPP
diff --git a/tests/test_MeshFlatNodeBoundary.cpp b/tests/test_MeshFlatNodeBoundary.cpp
index 9b9169f73c4306f9591ba0a9c0f5c1c622f7ef69..760f5564226f77fa53ec8de15f3f85552f02cbc4 100644
--- a/tests/test_MeshFlatNodeBoundary.cpp
+++ b/tests/test_MeshFlatNodeBoundary.cpp
@@ -1133,6 +1133,211 @@ TEST_CASE("MeshFlatNodeBoundary", "[mesh]")
}
}
+ SECTION("rotated diamond")
+ {
+ SECTION("2D")
+ {
+ static constexpr size_t Dimension = 2;
+
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ using R2 = TinyVector<2>;
+
+ auto T = [](const R2& x) -> R2 { return R2{x[0] + 0.1 * x[1], x[1] + 0.1 * x[0]}; };
+
+ SECTION("cartesian 2d")
+ {
+ std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian2DMesh();
+
+ const ConnectivityType& connectivity = p_mesh->connectivity();
+
+ auto xr = p_mesh->xr();
+
+ NodeValue<R2> rotated_xr{connectivity};
+
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { rotated_xr[node_id] = T(xr[node_id]); });
+
+ MeshType mesh{p_mesh->shared_connectivity(), rotated_xr};
+
+ {
+ const std::set<size_t> tag_set = {0, 1, 2, 3};
+
+ for (auto tag : tag_set) {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+ const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor);
+
+ auto node_list = get_node_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(node_boundary.nodeList(), node_list));
+
+ R2 normal = zero;
+
+ switch (tag) {
+ case 0: {
+ normal = 1. / std::sqrt(1.01) * R2{-1, 0.1};
+ break;
+ }
+ case 1: {
+ normal = 1. / std::sqrt(1.01) * R2{1, -0.1};
+ break;
+ }
+ case 2: {
+ normal = 1. / std::sqrt(1.01) * R2{0.1, -1};
+ break;
+ }
+ case 3: {
+ normal = 1. / std::sqrt(1.01) * R2{-0.1, 1};
+ break;
+ }
+ default: {
+ FAIL("unexpected tag number");
+ }
+ }
+ REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+
+ {
+ const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"};
+
+ for (const auto& name : name_set) {
+ NamedBoundaryDescriptor named_boundary_descriptor(name);
+ const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor);
+
+ auto node_list = get_node_list_from_name(name, connectivity);
+ REQUIRE(is_same(node_boundary.nodeList(), node_list));
+
+ R2 normal = zero;
+
+ if (name == "XMIN") {
+ normal = 1. / std::sqrt(1.01) * R2{-1, 0.1};
+ } else if (name == "XMAX") {
+ normal = 1. / std::sqrt(1.01) * R2{1, -0.1};
+ } else if (name == "YMIN") {
+ normal = 1. / std::sqrt(1.01) * R2{0.1, -1};
+ } else if (name == "YMAX") {
+ normal = 1. / std::sqrt(1.01) * R2{-0.1, 1};
+ } else {
+ FAIL("unexpected name: " + name);
+ }
+
+ REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ static constexpr size_t Dimension = 3;
+
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ using R3 = TinyVector<3>;
+
+ auto T = [](const R3& x) -> R3 {
+ return R3{x[0] + 0.1 * x[1] + 0.2 * x[2], x[1] + 0.1 * x[0] + 0.1 * x[2], x[2] + 0.1 * x[0]};
+ };
+
+ SECTION("cartesian 3d")
+ {
+ std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian3DMesh();
+
+ const ConnectivityType& connectivity = p_mesh->connectivity();
+
+ auto xr = p_mesh->xr();
+
+ NodeValue<R3> rotated_xr{connectivity};
+
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { rotated_xr[node_id] = T(xr[node_id]); });
+
+ MeshType mesh{p_mesh->shared_connectivity(), rotated_xr};
+
+ {
+ const std::set<size_t> tag_set = {0, 1, 2, 3, 4, 5};
+
+ for (auto tag : tag_set) {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+ const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor);
+
+ auto node_list = get_node_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(node_boundary.nodeList(), node_list));
+
+ R3 normal = zero;
+
+ switch (tag) {
+ case 0: {
+ normal = R3{-0.977717523265611, 0.0977717523265611, 0.185766329420466};
+ break;
+ }
+ case 1: {
+ normal = R3{0.977717523265611, -0.0977717523265612, -0.185766329420466};
+ break;
+ }
+ case 2: {
+ normal = R3{0.0911512175788074, -0.992535480302569, 0.0810233045144955};
+ break;
+ }
+ case 3: {
+ normal = R3{-0.0911512175788074, 0.992535480302569, -0.0810233045144955};
+ break;
+ }
+ case 4: {
+ normal = R3{0.100493631166705, -0.0100493631166705, -0.994886948550377};
+ break;
+ }
+ case 5: {
+ normal = R3{-0.100493631166705, 0.0100493631166705, 0.994886948550377};
+ break;
+ }
+ default: {
+ FAIL("unexpected tag number");
+ }
+ }
+
+ REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+
+ {
+ const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+
+ for (const auto& name : name_set) {
+ NamedBoundaryDescriptor named_boundary_descriptor(name);
+ const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor);
+
+ auto node_list = get_node_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(node_boundary.nodeList(), node_list));
+
+ R3 normal = zero;
+
+ if (name == "XMIN") {
+ normal = R3{-0.977717523265611, 0.0977717523265611, 0.185766329420466};
+ } else if (name == "XMAX") {
+ normal = R3{0.977717523265611, -0.0977717523265612, -0.185766329420466};
+ } else if (name == "YMIN") {
+ normal = R3{0.0911512175788074, -0.992535480302569, 0.0810233045144955};
+ } else if (name == "YMAX") {
+ normal = R3{-0.0911512175788074, 0.992535480302569, -0.0810233045144955};
+ } else if (name == "ZMIN") {
+ normal = R3{0.100493631166705, -0.0100493631166705, -0.994886948550377};
+ } else if (name == "ZMAX") {
+ normal = R3{-0.100493631166705, 0.0100493631166705, 0.994886948550377};
+ } else {
+ FAIL("unexpected name: " + name);
+ }
+
+ REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+ }
+
SECTION("curved mesh")
{
SECTION("2D")
diff --git a/tests/test_OFStream.cpp b/tests/test_OFStream.cpp
index fbb5818dbc70451ff212c820edacfa595963e961..6ca7d98df0a6f36bf30add2b891a659beba757c2 100644
--- a/tests/test_OFStream.cpp
+++ b/tests/test_OFStream.cpp
@@ -13,7 +13,8 @@ TEST_CASE("OFStream", "[language]")
{
SECTION("ofstream")
{
- const std::string basename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("ofstream_");
+ const std::string basename =
+ std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("ofstream_dir").append("ofstream_");
const std::string filename = basename + stringify(parallel::rank());
// Ensures that the file is closed after this line
@@ -47,13 +48,4 @@ TEST_CASE("OFStream", "[language]")
REQUIRE(not std::filesystem::exists(filename));
}
-
- SECTION("invalid filename")
- {
- if (parallel::rank() == 0) {
- const std::string filename = "badpath/invalidpath/ofstream";
-
- REQUIRE_THROWS_WITH(std::make_shared<OFStream>(filename), "error: cannot create file " + filename);
- }
- }
}