From 227c0436b2b4bb870d6b8d9fb6828b7f435d6064 Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Thu, 22 Sep 2005 12:50:02 +0000 Subject: [PATCH] New netcdf GridIn::Format. New read_netcdf function. git-svn-id: https://svn.dealii.org/trunk@11497 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_in.h | 12 +- deal.II/deal.II/source/grid/grid_in.cc | 411 ++++++++++++++++++++++++- 2 files changed, 420 insertions(+), 3 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_in.h b/deal.II/deal.II/include/grid/grid_in.h index 0c9330f767..8f65ae5604 100644 --- a/deal.II/deal.II/include/grid/grid_in.h +++ b/deal.II/deal.II/include/grid/grid_in.h @@ -180,7 +180,9 @@ class GridIn /// Use read_xda() xda, /// Use read_msh() - msh + msh, + /// Use read_netcdf() + netcdf }; /** @@ -237,6 +239,14 @@ class GridIn */ void read_msh (std::istream &); + /** + * Read grid data from a netcdf + * file. This requires the + * library to be linked with the + * netcdf library. + */ + void read_netcdf (const std::string &); + /** * Returns the standard suffix * for a file in this format. diff --git a/deal.II/deal.II/source/grid/grid_in.cc b/deal.II/deal.II/source/grid/grid_in.cc index 7761707bce..3319d4d68f 100644 --- a/deal.II/deal.II/source/grid/grid_in.cc +++ b/deal.II/deal.II/source/grid/grid_in.cc @@ -21,6 +21,10 @@ #include #include +#ifdef DEAL_II_HAVE_NETCDF +#include +#endif + template GridIn::GridIn () : @@ -734,6 +738,388 @@ void GridIn::read_msh (std::istream &in) } +#if deal_II_dimension == 1 + +template <> +void GridIn<1>::read_netcdf (const std::string &filename) +{ + AssertThrow(false, ExcImpossibleInDim(1)); +} + +#endif + + + +#if deal_II_dimension == 2 + +template <> +void GridIn<2>::read_netcdf (const std::string &filename) +{ +#ifndef DEAL_II_HAVE_NETCDF + AssertThrow(false, ExcNeedsNetCDF()); +#else + const unsigned int dim=2; + const bool output=false; + Assert (tria != 0, ExcNoTriangulationSelected()); + // this function assumes the TAU + // grid format. + // + // This format stores 2d grids as + // 3d grids. In particular, a 2d + // grid of n_cells quadrilaterals + // in the y=0 plane is duplicated + // to y=1 to build n_cells + // hexaeders. The surface + // quadrilaterals of this 3d grid + // are marked with boundary + // marker. In the following we read + // in all data required, find the + // boundary marker associated with + // the plane y=0, and extract the + // corresponding 2d data to build a + // Triangulation<2>. + + // In the following, we assume that + // the 2d grid lies in the x-z + // plane (y=0). I.e. we choose: + // point[coord]=0, with coord=1 + const unsigned int coord=1; + // Also x-y-z (0-1-2) point + // coordinates will be transformed + // to x-y (x2d-y2d) coordinates. + // With coord=1 as above, we have + // x-z (0-2) -> (x2d-y2d) + const unsigned int x2d=0; + const unsigned int y2d=2; + // For the case, the 2d grid lies + // in x-y or y-z plane instead, the + // following code must be extended + // to find the right value for + // coord, and setting x2d and y2d + // accordingly. + + // First, open the file + NcFile nc (filename.c_str()); + AssertThrow(nc.is_valid(), ExcIO()); + + // then read n_cells + NcDim *elements_dim=nc.get_dim("no_of_elements"); + AssertThrow(elements_dim->is_valid(), ExcIO()); + const unsigned int n_cells=elements_dim->size(); + + // then we read + // int marker(no_of_markers) + NcDim *marker_dim=nc.get_dim("no_of_markers"); + AssertThrow(marker_dim->is_valid(), ExcIO()); + const unsigned int n_markers=marker_dim->size(); + + NcVar *marker_var=nc.get_var("marker"); + AssertThrow(marker_var->is_valid(), ExcIO()); + AssertThrow(marker_var->num_dims()==1, ExcIO()); + AssertThrow(static_cast( + marker_var->get_dim(0)->size())==n_markers, ExcIO()); + + vector marker(n_markers); + marker_var->get(marker.begin(), n_markers); + + if (output) + { + cout << "n_cell=" << n_cells << endl; + cout << "marker: "; + for (unsigned int i=0; iis_valid(), ExcIO()); + const unsigned int n_bquads=bquads_dim->size(); + + NcVar *bmarker_var=nc.get_var("boundarymarker_of_surfaces"); + AssertThrow(bmarker_var->is_valid(), ExcIO()); + AssertThrow(bmarker_var->num_dims()==1, ExcIO()); + AssertThrow(static_cast( + bmarker_var->get_dim(0)->size())==n_bquads, ExcIO()); + + vector bmarker(n_bquads); + bmarker_var->get(bmarker.begin(), n_bquads); + + // for each marker count the + // number of boundary quads + // which carry this marker + map n_bquads_per_bmarker; + for (unsigned int i=0; i::const_iterator + iter=n_bquads_per_bmarker.begin(); + for (; iter!=n_bquads_per_bmarker.end(); ++iter) + cout << " n_bquads_per_bmarker[" << iter->first + << "]=" << iter->second << endl; + } + + // next we read + // int points_of_surfacequadrilaterals( + // no_of_surfacequadrilaterals, + // points_per_surfacequadrilateral) + NcDim *quad_vertices_dim=nc.get_dim("points_per_surfacequadrilateral"); + AssertThrow(quad_vertices_dim->is_valid(), ExcIO()); + const unsigned int vertices_per_quad=quad_vertices_dim->size(); + AssertThrow(vertices_per_quad==GeometryInfo::vertices_per_cell, ExcIO()); + + NcVar *vertex_indices_var=nc.get_var("points_of_surfacequadrilaterals"); + AssertThrow(vertex_indices_var->is_valid(), ExcIO()); + AssertThrow(vertex_indices_var->num_dims()==2, ExcIO()); + AssertThrow(static_cast( + vertex_indices_var->get_dim(0)->size())==n_bquads, ExcIO()); + AssertThrow(static_cast( + vertex_indices_var->get_dim(1)->size())==vertices_per_quad, ExcIO()); + + vector vertex_indices(n_bquads*vertices_per_quad); + vertex_indices_var->get(vertex_indices.begin(), n_bquads, vertices_per_quad); + + for (unsigned int i=0; i=0, ExcInternalError()); + + if (output) + { + cout << "vertex_indices:" << endl; + for (unsigned int i=0, v=0; iis_valid(), ExcIO()); + const unsigned int n_vertices=vertices_dim->size(); + if (output) + cout << "n_vertices=" << n_vertices << endl; + + NcVar *points_xc=nc.get_var("points_xc"); + NcVar *points_yc=nc.get_var("points_yc"); + NcVar *points_zc=nc.get_var("points_zc"); + AssertThrow(points_xc->is_valid(), ExcIO()); + AssertThrow(points_yc->is_valid(), ExcIO()); + AssertThrow(points_zc->is_valid(), ExcIO()); + AssertThrow(points_xc->num_dims()==1, ExcIO()); + AssertThrow(points_yc->num_dims()==1, ExcIO()); + AssertThrow(points_zc->num_dims()==1, ExcIO()); + AssertThrow(points_yc->get_dim(0)->size()== + static_cast(n_vertices), ExcIO()); + AssertThrow(points_zc->get_dim(0)->size()== + static_cast(n_vertices), ExcIO()); + AssertThrow(points_xc->get_dim(0)->size()== + static_cast(n_vertices), ExcIO()); + vector > point_values(3, vector (n_vertices)); + points_xc->get(point_values[0].begin(), n_vertices); + points_yc->get(point_values[1].begin(), n_vertices); + points_zc->get(point_values[2].begin(), n_vertices); + + // and fill the vertices + std::vector > vertices (n_vertices); + for (unsigned int i=0; iis_valid(), ExcIO()); + const unsigned int n_cells=elements_dim->size(); + if (output) + cout << "n_cell=" << n_cells << endl; + // and n_hexes + NcDim *hexes_dim=nc.get_dim("no_of_hexaeders"); + AssertThrow(hexes_dim->is_valid(), ExcIO()); + const unsigned int n_hexes=hexes_dim->size(); + AssertThrow(n_hexes==n_cells, + ExcMessage("deal.II can handle purely hexaedral grids, only.")); + + // next we read + // int points_of_hexaeders( + // no_of_hexaeders, + // points_per_hexaeder) + NcDim *hex_vertices_dim=nc.get_dim("points_per_hexaeder"); + AssertThrow(hex_vertices_dim->is_valid(), ExcIO()); + const unsigned int vertices_per_hex=hex_vertices_dim->size(); + AssertThrow(vertices_per_hex==GeometryInfo::vertices_per_cell, ExcIO()); + + NcVar *vertex_indices_var=nc.get_var("points_of_hexaeders"); + AssertThrow(vertex_indices_var->is_valid(), ExcIO()); + AssertThrow(vertex_indices_var->num_dims()==2, ExcIO()); + AssertThrow(static_cast( + vertex_indices_var->get_dim(0)->size())==n_cells, ExcIO()); + AssertThrow(static_cast( + vertex_indices_var->get_dim(1)->size())==vertices_per_hex, ExcIO()); + + vector vertex_indices(n_cells*vertices_per_hex); + vertex_indices_var->get(vertex_indices.begin(), n_cells, vertices_per_hex); + + for (unsigned int i=0; i=0, ExcInternalError()); + + if (output) + { + cout << "vertex_indices:" << endl; + for (unsigned int cell=0, v=0; cellis_valid(), ExcIO()); + const unsigned int n_vertices=vertices_dim->size(); + if (output) + cout << "n_vertices=" << n_vertices << endl; + + NcVar *points_xc=nc.get_var("points_xc"); + NcVar *points_yc=nc.get_var("points_yc"); + NcVar *points_zc=nc.get_var("points_zc"); + AssertThrow(points_xc->is_valid(), ExcIO()); + AssertThrow(points_yc->is_valid(), ExcIO()); + AssertThrow(points_zc->is_valid(), ExcIO()); + AssertThrow(points_xc->num_dims()==1, ExcIO()); + AssertThrow(points_yc->num_dims()==1, ExcIO()); + AssertThrow(points_zc->num_dims()==1, ExcIO()); + AssertThrow(points_yc->get_dim(0)->size()== + static_cast(n_vertices), ExcIO()); + AssertThrow(points_zc->get_dim(0)->size()== + static_cast(n_vertices), ExcIO()); + AssertThrow(points_xc->get_dim(0)->size()== + static_cast(n_vertices), ExcIO()); + vector > point_values(3, vector (n_vertices)); + points_xc->get(point_values[0].begin(), n_vertices); + points_yc->get(point_values[1].begin(), n_vertices); + points_zc->get(point_values[2].begin(), n_vertices); + + // and fill the vertices + std::vector > vertices (n_vertices); + for (unsigned int i=0; i > cells(n_cells); + for (unsigned int cell=0; cell::reorder_cells (cells); + tria->create_triangulation (vertices, cells, subcelldata); +#endif +} + +#endif + template @@ -1022,7 +1408,10 @@ void GridIn::read (const std::string& filename, format = parse_format(ext); } } - read(in, format); + if (format == netcdf) + read_netcdf(filename); + else + read(in, format); } @@ -1051,6 +1440,16 @@ void GridIn::read (std::istream& in, read_xda (in); return; + case netcdf: + // there is no + // read_netcdf(istream &) + // function. Use the + // read_netcdf(const + // string &filename + // function instead + Assert(false, ExcNotImplemented()); + return; + case Default: break; } @@ -1073,6 +1472,8 @@ GridIn::default_suffix (const Format format) return ".inp"; case xda: return ".xda"; + case netcdf: + return ".nc"; default: Assert (false, ExcNotImplemented()); return ".unknown_format"; @@ -1100,6 +1501,12 @@ GridIn::parse_format (const std::string &format_name) if (format_name == "xda") return xda; + if (format_name == "netcdf") + return netcdf; + + if (format_name == "nc") + return netcdf; + AssertThrow (false, ExcInvalidState ()); // return something weird return Format(Default); @@ -1110,7 +1517,7 @@ GridIn::parse_format (const std::string &format_name) template std::string GridIn::get_format_names () { - return "dbmesh|msh|ucd|xda"; + return "dbmesh|msh|ucd|xda|netcdf"; } -- 2.39.5