From 14ce1a731a00cc2e0d0facd8f4548f683759abbf Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 11 May 2018 17:31:07 +0200
Subject: [PATCH] Improved 1d mesh management
- Uses a new mesh interface (pretty crap, using dimension to deduce mesh kind)
- 1d mesh reader almost functionnal (missing boundaries)
- 1d mesh based scheme plugged in a really ugly way (proof of concept)
---
src/main.cpp | 229 +++++++++++++++++++++++++++-------------
src/mesh/GmshReader.cpp | 19 ++--
src/mesh/GmshReader.hpp | 11 +-
src/mesh/Mesh.hpp | 13 ++-
4 files changed, 181 insertions(+), 91 deletions(-)
diff --git a/src/main.cpp b/src/main.cpp
index fb2da6d32..8b76ec57f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -124,104 +124,183 @@ int main(int argc, char *argv[])
if (filename != "") {
std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
GmshReader gmsh_reader(filename);
+ std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
+
+ switch (p_mesh->meshDimension()) {
+ case 1: {
+ typedef Connectivity1D ConnectivityType;
+ typedef Mesh<ConnectivityType> MeshType;
+ typedef MeshData<MeshType> MeshDataType;
+ typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+ Kokkos::Timer timer;
+ timer.reset();
+ MeshDataType mesh_data(mesh);
+
+ std::vector<BoundaryConditionHandler> bc_list;
+ // { // quite dirty!
+ // for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
+ // ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
+ // unsigned int ref = nodes_boundary.first;
+ // TinyVector<2> normal(0,0);
+ // if ((ref == 5) or (ref == 6)) {
+ // normal = TinyVector<2>(0,1);
+ // } else {
+ // normal = TinyVector<2>(1,0);
+ // }
+ // const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second;
+ // std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0));
+ // for (size_t r=0; r<nodes_ids.extent(0); ++r) {
+ // node_boundary_vector[r] = nodes_ids[r];
+ // }
+ // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+ // = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal);
+ // std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+ // bc_list.push_back(BoundaryConditionHandler(bc));
+ // }
+ // }
+
+ UnknownsType unknowns(mesh_data);
+
+ unknowns.initializeSod();
+
+ AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
+
+ typedef TinyVector<MeshType::dimension> Rd;
+
+ const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
+
+ const double tmax=0.2;
+ double t=0;
+
+ int itermax=std::numeric_limits<int>::max();
+ int iteration=0;
+
+ Kokkos::View<double*> rhoj = unknowns.rhoj();
+ Kokkos::View<double*> ej = unknowns.ej();
+ Kokkos::View<double*> pj = unknowns.pj();
+ Kokkos::View<double*> gammaj = unknowns.gammaj();
+ Kokkos::View<double*> cj = unknowns.cj();
+
+ BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+
+ VTKWriter vtk_writer("mesh", 0.01);
+
+ while((t<tmax) and (iteration<itermax)) {
+ vtk_writer.write(mesh, t);
+ double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+ if (t+dt>tmax) {
+ dt=tmax-t;
+ }
+ acoustic_solver.computeNextStep(t,dt, unknowns);
- typedef Mesh<Connectivity2D> MeshType;
- typedef MeshData<MeshType> MeshDataType;
- typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+ block_eos.updatePandCFromRhoE();
+
+ t += dt;
+ ++iteration;
+ }
+ vtk_writer.write(mesh, t, true); // forces last output
- MeshType& mesh = gmsh_reader.mesh();
+ std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
+ << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
- Kokkos::Timer timer;
- timer.reset();
- MeshDataType mesh_data(mesh);
+ method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
- std::vector<BoundaryConditionHandler> bc_list;
- { // quite dirty!
- for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
- Connectivity2D::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
- unsigned int ref = nodes_boundary.first;
- TinyVector<2> normal(0,0);
- if ((ref == 5) or (ref == 6)) {
- normal = TinyVector<2>(0,1);
- } else {
- normal = TinyVector<2>(1,0);
- }
- const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second;
- std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0));
- for (size_t r=0; r<nodes_ids.extent(0); ++r) {
- node_boundary_vector[r] = nodes_ids[r];
+ break;
+ }
+ case 2: {
+ typedef Connectivity2D ConnectivityType;
+ typedef Mesh<ConnectivityType> MeshType;
+ typedef MeshData<MeshType> MeshDataType;
+ typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+ Kokkos::Timer timer;
+ timer.reset();
+ MeshDataType mesh_data(mesh);
+
+ std::vector<BoundaryConditionHandler> bc_list;
+ { // quite dirty!
+ for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
+ ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
+ unsigned int ref = nodes_boundary.first;
+ TinyVector<2> normal(0,0);
+ if ((ref == 5) or (ref == 6)) {
+ normal = TinyVector<2>(0,1);
+ } else {
+ normal = TinyVector<2>(1,0);
+ }
+ const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second;
+ std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0));
+ for (size_t r=0; r<nodes_ids.extent(0); ++r) {
+ node_boundary_vector[r] = nodes_ids[r];
+ }
+ SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+ = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal);
+ std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+ bc_list.push_back(BoundaryConditionHandler(bc));
}
- SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
- = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal);
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
}
- // PressureBoundaryCondition* pres_bc1
- // = new PressureBoundaryCondition(0.1,
- // std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
- // std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
- // bc_list.push_back(BoundaryConditionHandler(bc1));
- }
+ UnknownsType unknowns(mesh_data);
- UnknownsType unknowns(mesh_data);
-
- unknowns.initializeSod();
+ unknowns.initializeSod();
- AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
+ AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- typedef TinyVector<MeshType::dimension> Rd;
+ typedef TinyVector<MeshType::dimension> Rd;
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
- const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
+ const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
- const double tmax=0.2;
- double t=0;
+ const double tmax=0.2;
+ double t=0;
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
+ int itermax=std::numeric_limits<int>::max();
+ int iteration=0;
- Kokkos::View<double*> rhoj = unknowns.rhoj();
- Kokkos::View<double*> ej = unknowns.ej();
- Kokkos::View<double*> pj = unknowns.pj();
- Kokkos::View<double*> gammaj = unknowns.gammaj();
- Kokkos::View<double*> cj = unknowns.cj();
+ Kokkos::View<double*> rhoj = unknowns.rhoj();
+ Kokkos::View<double*> ej = unknowns.ej();
+ Kokkos::View<double*> pj = unknowns.pj();
+ Kokkos::View<double*> gammaj = unknowns.gammaj();
+ Kokkos::View<double*> cj = unknowns.cj();
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+ BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
- VTKWriter vtk_writer("mesh", 0.01);
+ VTKWriter vtk_writer("mesh", 0.01);
- while((t<tmax) and (iteration<itermax)) {
- vtk_writer.write(mesh, t);
- double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
- if (t+dt>tmax) {
- dt=tmax-t;
- }
- acoustic_solver.computeNextStep(t,dt, unknowns);
+ while((t<tmax) and (iteration<itermax)) {
+ vtk_writer.write(mesh, t);
+ double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+ if (t+dt>tmax) {
+ dt=tmax-t;
+ }
+ acoustic_solver.computeNextStep(t,dt, unknowns);
- block_eos.updatePandCFromRhoE();
+ block_eos.updatePandCFromRhoE();
- t += dt;
- ++iteration;
- }
- vtk_writer.write(mesh, t, true); // forces last output
-
- std::ofstream gnuplot("sol.gnu");
- for (size_t j=0; j<mesh.numberOfCells(); ++j) {
- for (int r=0; r<mesh.connectivity().cellNbNodes()[j]; ++r) {
- const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,r)];
- gnuplot << x[0] << ' ' << x[1] << '\n';
+ t += dt;
+ ++iteration;
}
- const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,0)];
- gnuplot << x[0] << ' ' << x[1] << "\n\n";
- }
- gnuplot.close();
+ vtk_writer.write(mesh, t, true); // forces last output
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+ std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
+ << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+ method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+ break;
+ }
+ case 3: {
+ break;
+ }
+ }
+ std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
+
} else {
// class for acoustic solver test
Kokkos::Timer timer;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 18e301b35..ba398607b 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -5,6 +5,7 @@
#include <set>
#include <rang.hpp>
+#include <Connectivity1D.hpp>
#include <Connectivity2D.hpp>
#include <Mesh.hpp>
#include <map>
@@ -805,7 +806,7 @@ GmshReader::__proceedData()
for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
cell_nb_nodes[j] = 4;
}
- Connectivity2D* p_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
+ std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_nb_nodes, cell_nodes));
Connectivity2D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -856,7 +857,7 @@ GmshReader::__proceedData()
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
}
- m_mesh = new MeshType(std::shared_ptr<Connectivity2D>(p_connectivity), xr);
+ m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
} else if ((dimension1_mask, elementNumber)>0) {
@@ -874,13 +875,17 @@ GmshReader::__proceedData()
cell_nb_nodes[j] = 2;
}
- // std::shared_ptr<Connectivity1D> connectivity(new Connectivity1D());
- // m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
- // Connectivity2D& connectivity = *m_connectivity;
+ std::shared_ptr<Connectivity1D> connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes));
+ typedef Mesh<Connectivity1D> MeshType;
+ typedef TinyVector<1, double> Rd;
+ Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+ for (size_t i=0; i<__vertices.extent(0); ++i) {
+ xr[i][0] = __vertices[i][0];
+ }
+
+ m_mesh = std::shared_ptr<IMesh>(new MeshType(connectivity, xr));
- std::cerr << "*** using a 1d mesh (Implementation in progress)\n";
- std::exit(0);
} else {
std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
std::exit(0);
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 1baddb1ee..fb41a5c35 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -184,17 +184,12 @@ private:
*/
void __readPhysicalNames2_2();
- Mesh<Connectivity2D>* m_mesh;
+ std::shared_ptr<IMesh> m_mesh;
public:
- // Connectivity2D& connectivity()
- // {
- // return *m_connectivity;
- // }
-
- Mesh<Connectivity2D>& mesh()
+ std::shared_ptr<IMesh> mesh() const
{
- return *m_mesh;
+ return m_mesh;
}
GmshReader(const std::string& filename);
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 40936ce84..2b05731f5 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -6,8 +6,14 @@
#include <memory>
+struct IMesh
+{
+ virtual const size_t meshDimension() const = 0;
+ ~IMesh() = default;
+};
+
template <typename ConnectivityType>
-class Mesh
+class Mesh final : public IMesh
{
public:
typedef ConnectivityType Connectivity;
@@ -20,6 +26,11 @@ private:
public:
+ const size_t meshDimension() const
+ {
+ return dimension;
+ }
+
const Connectivity& connectivity() const
{
return *m_connectivity;
--
GitLab