From e8ef00b8a113991613946570300e36456575031e Mon Sep 17 00:00:00 2001 From: AlePalu Date: Wed, 18 Dec 2024 09:52:54 +0100 Subject: [PATCH] factory methods for intervals, rectangles, squares and cubic meshes --- fdaPDE/geometry/interval.h | 4 ++ fdaPDE/geometry/triangulation.h | 79 +++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/fdaPDE/geometry/interval.h b/fdaPDE/geometry/interval.h index 21c5d1b..6b43e00 100644 --- a/fdaPDE/geometry/interval.h +++ b/fdaPDE/geometry/interval.h @@ -144,6 +144,10 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula mutable typename Base::CellType cell_; // used in case cell caching is off }; +// interval mesh factory methods +inline Triangulation<1, 1> make_interval(double a, double b, int n_nodes) { return Triangulation<1, 1>(a, b, n_nodes); } +inline Triangulation<1, 1> make_unit_interval(int n_nodes) { return make_interval(0.0, 1.0, n_nodes); } + } // namespace fdapde #endif // __INTERVAL_H__ diff --git a/fdaPDE/geometry/triangulation.h b/fdaPDE/geometry/triangulation.h index 55fa867..76b1f37 100644 --- a/fdaPDE/geometry/triangulation.h +++ b/fdaPDE/geometry/triangulation.h @@ -466,6 +466,41 @@ template class Triangulation<2, N> : public TriangulationBase<2, N, Tria mutable typename Base::CellType cell_; // used in case cell caching is off }; +// square mesh factory methods +inline Triangulation<2, 2> +make_rectangle(double a_x, double b_x, double a_y, double b_y, int nx, int ny) { + // allocate memory + Eigen::Matrix nodes(nx * ny, 2); + Eigen::Matrix cells((nx - 1) * (ny - 1) * 2, 3); // each subrectangle splitted in 2 + Eigen::Matrix boundary(nx * ny, 1); + // compute steps along axis + double h_x = (b_x - a_x) / (nx - 1); + double h_y = (b_y - a_y) / (ny - 1); + int cell_id = 0; + std::array v; + for (int j = 0; j < ny; ++j) { + for (int i = 0; i < nx; ++i) { + int node_id = j * nx + i; + nodes.row(node_id) << (a_x + i * h_x), (a_y + j * h_y); // node coordinates + if (i < nx - 1 && j < ny - 1) { + // compute vector of vertices of subrectangle + int p = j * nx + i; + v = {p, p + 1, p + nx, p + nx + 1}; + // compute vertices of each triangle in the subrectangle + cells.row(cell_id + 0) << v[0], v[1], v[2]; + cells.row(cell_id + 1) << v[1], v[2], v[3]; + cell_id = cell_id + 2; + } + boundary(node_id, 0) = (i == 0 || i == nx - 1 || j == 0 || j == ny - 1) ? 1 : 0; + } + } + return Triangulation<2, 2>(nodes, cells, boundary); +} +inline Triangulation<2, 2> make_square(double a, double b, int n_nodes) { + return make_rectangle(a, b, a, b, n_nodes, n_nodes); +} +inline Triangulation<2, 2> make_unit_square(int n_nodes) { return make_square(0.0, 1.0, n_nodes); } + // face-based storage template <> class Triangulation<3, 3> : public TriangulationBase<3, 3, Triangulation<3, 3>> { public: @@ -776,6 +811,50 @@ template <> class Triangulation<3, 3> : public TriangulationBase<3, 3, Triangula mutable typename Base::CellType cell_; // used in case cell caching is off }; +// cubic mesh factory methods (nx, ny, nz: number of nodes along x,y,z axis respectively) +inline Triangulation<3, 3> +make_parallelepiped(double a_x, double b_x, double a_y, double b_y, double a_z, double b_z, int nx, int ny, int nz) { + // allocate memory + Eigen::Matrix nodes(nx * ny * nz, 3); + Eigen::Matrix cells((nx - 1) * (ny - 1) * (nz - 1) * 6, 4); // each subcube spliited in 6 + Eigen::Matrix boundary(nx * ny * nz, 1); + // compute steps along axis + double h_x = (b_x - a_x) / (nx - 1); + double h_y = (b_y - a_y) / (ny - 1); + double h_z = (b_z - a_z) / (nz - 1); + int cell_id = 0; + std::array v; + for (int k = 0; k < nz; ++k) { + for (int j = 0; j < ny; ++j) { + for (int i = 0; i < nx; ++i) { + int node_id = k * nx * ny + j * nx + i; + nodes.row(node_id) << (a_x + i * h_x), (a_y + j * h_y), (a_z + k * h_z); + if (i < nx - 1 && j < ny - 1 && k < nz - 1) { + // build vector of vertices of sub-cube + int p = k * nx * ny + j * nx + i; + v = {p, p + 1, p + nx, p + nx + 1, + p + nx * ny, p + nx * ny + 1, p + nx * ny + nx, p + nx * ny + nx + 1}; + // split sub-cube in 6 tetrahedra (this guarantees a conforming triangulation) + cells.row(cell_id + 0) << v[6], v[2], v[1], v[0]; + cells.row(cell_id + 1) << v[6], v[4], v[1], v[0]; + cells.row(cell_id + 2) << v[6], v[5], v[7], v[1]; + cells.row(cell_id + 3) << v[6], v[3], v[7], v[1]; + cells.row(cell_id + 4) << v[6], v[3], v[2], v[1]; + cells.row(cell_id + 5) << v[6], v[5], v[4], v[1]; + cell_id = cell_id + 6; + } + boundary(node_id, 0) = + (i == 0 || i == nx - 1 || j == 0 || j == ny - 1 || k == 0 || k == nz - 1) ? 1 : 0; + } + } + } + return Triangulation<3, 3>(nodes, cells, boundary); +} +inline Triangulation<3, 3> make_cube(double a, double b, int n_nodes) { + return make_parallelepiped(a, b, a, b, a, b, n_nodes, n_nodes, n_nodes); +} +inline Triangulation<3, 3> make_unit_cube(int n_nodes) { return make_cube(0.0, 1.0, n_nodes); } + } // namespace fdapde #endif // __TRIANGULATION_H__