<ol>
+ <li> Fixed: In 1d, GridIn::read_msh() ignored boundary indicators
+ associated with vertices. This is now fixed.
+ <br>
+ (Jan Strebel, Wolfgang Bangerth, 2015/08/25)
+ </li>
+
<li> Improved: The interface and documentation for periodic boundary
conditions have been restructured. A
@ref GlossPeriodicConstraints "glossary entry" has been written.
DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+ /**
+ * In 1d, boundary indicators are associated with vertices, but this is not
+ * currently passed through the SubcellData structure. This function sets
+ * boundary indicators on vertices after the triangulation has already been
+ * created.
+ *
+ * TODO: This this properly via SubcellData
+ */
+ template <int spacedim>
+ void
+ assign_1d_boundary_indicators (const std::map<unsigned int, types::boundary_id> &boundary_ids,
+ Triangulation<1,spacedim> &triangulation)
+ {
+ if (boundary_ids.size() > 0)
+ for (typename Triangulation<1,spacedim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<1>::faces_per_cell; ++f)
+ if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end())
+ {
+ AssertThrow (cell->at_boundary(f),
+ ExcMessage ("You are trying to prescribe boundary ids on the face "
+ "of a 1d cell (i.e., on a vertex), but this face is not actually at "
+ "the boundary of the mesh. This is not allowed."));
+ cell->face(f)->set_boundary_id (boundary_ids.find(cell->vertex_index(f))->second);
+ }
+ }
+
+
+ template <int dim, int spacedim>
+ void
+ assign_1d_boundary_indicators (const std::map<unsigned int, types::boundary_id> &,
+ Triangulation<dim,spacedim> &)
+ {
+ // we shouldn't get here since boundary ids are not assigned to
+ // vertices except in 1d
+ Assert (dim != 1, ExcInternalError());
+ }
+}
+
template <int dim, int spacedim>
GridIn<dim, spacedim>::GridIn () :
tria(0, typeid(*this).name()), default_format(ucd)
in >> n_cells;
- // set up array of cells
- std::vector<CellData<dim> > cells;
- SubCellData subcelldata;
+ // set up array of cells and subcells (faces). In 1d, there is currently no
+ // standard way in deal.II to pass boundary indicators attached to individual
+ // vertices, so do this by hand via the boundary_ids_1d array
+ std::vector<CellData<dim> > cells;
+ SubCellData subcelldata;
+ std::map<unsigned int, types::boundary_id> boundary_ids_1d;
for (unsigned int cell=0; cell<n_cells; ++cell)
{
case 2:
{
- // read the tags; ignore
- // all but the first one
+ // read the tags; ignore all but the first one which we will
+ // interpret as the material_id (for cells) or boundary_id
+ // (for faces)
unsigned int n_tags;
in >> n_tags;
if (n_tags > 0)
subcelldata.boundary_quads.back().vertices[i]));
subcelldata.boundary_quads.back().vertices[i] =
numbers::invalid_unsigned_int;
- };
+ }
}
else if (cell_type == 15)
{
- // Ignore vertices
- // but read the
- // number of nodes
- // given
+ // read the indices of nodes given
+ unsigned int node_index;
switch (gmsh_file_format)
{
case 1:
{
for (unsigned int i=0; i<nod_num; ++i)
- in >> dummy;
+ in >> node_index;
break;
}
case 2:
{
- in >> dummy;
+ in >> node_index;
break;
}
}
+
+ // we only care about boundary indicators assigned to individual
+ // vertices in 1d (because otherwise the vertices are not faces)
+ if (dim == 1)
+ boundary_ids_1d[vertex_indices[node_index]] = material_id;
}
else
// cannot read this, so throw
AssertThrow (false, ExcGmshUnsupportedGeometry(cell_type));
}
- };
+ }
// Assert we reached the end of the block
in >> line;
GridReordering<dim,spacedim>::invert_all_cells_of_negative_grid (vertices, cells);
GridReordering<dim,spacedim>::reorder_cells (cells);
tria->create_triangulation_compatibility (vertices, cells, subcelldata);
+
+ // in 1d, we also have to attach boundary ids to vertices, which does not
+ // currently work through the call above
+ if (dim == 1)
+ assign_1d_boundary_indicators (boundary_ids_1d, *tria);
}
// ---------------------------------------------------------------------
-// read a file in the MSH format used by the GMSH program. this
-// particular mesh has type-15 cells (nodes) with more than one
-// associated vertex. we failed to read the vertices
+// in 1d, we have to read vertex information to set boundary
+// indicators
+//
+// test case by Jan Strebel
+
#include "../tests.h"
#include <deal.II/dofs/dof_handler.h>
std::ofstream logfile("output");
-template<int dim>
-void check_file (const std::string name,
- typename GridIn<dim>::Format format)
+void check_file ()
{
- Triangulation<dim> tria;
- GridIn<dim> gi;
+ Triangulation<1> tria;
+ GridIn<1> gi;
gi.attach_triangulation (tria);
- gi.read(name, format);
- deallog << '\t' << tria.n_vertices()
- << '\t' << tria.n_cells()
- << std::endl;
+ std::ifstream in (SOURCE_DIR "/../grid/grids/grid_in_msh_02.msh");
+ gi.read_msh(in);
- GridOut grid_out;
- grid_out.write_gnuplot (tria, deallog.get_file_stream());
-}
-
-void filename_resolution()
-{
- check_file<2> (std::string(SOURCE_DIR "/grid_in_msh_02/mesh"), GridIn<2>::msh);
+ for (Triangulation<1>::active_cell_iterator cell = tria.begin_active(); cell != tria.end(); ++cell)
+ {
+ for (unsigned int face = 0; face < 2; ++face)
+ {
+ if (cell->at_boundary(face))
+ deallog << "vertex " << cell->face_index(face)
+ << " has boundary indicator " << (int)cell->face(face)->boundary_indicator()
+ << std::endl;
+ }
+ }
}
deallog.depth_console(0);
deallog.threshold_double(1.e-10);
- filename_resolution();
+ check_file ();
}
-DEAL:: 28 18
-2.0e+03 0.0 0 0
-2.0e+03 6.7e+02 0 0
-1.4e+03 5.0e+02 0 0
-1.5e+03 0.0 0 0
-2.0e+03 0.0 0 0
-
-
-1.5e+03 0.0 0 0
-1.4e+03 5.0e+02 0 0
-9.6e+02 4.0e+02 0 0
-1.0e+03 0.0 0 0
-1.5e+03 0.0 0 0
-
-
-1.0e+03 0.0 0 0
-9.6e+02 4.0e+02 0 0
-5.0e+02 3.3e+02 0 0
-5.0e+02 0.0 0 0
-1.0e+03 0.0 0 0
-
-
-2.0e+03 6.7e+02 0 0
-2.0e+03 1.3e+03 0 0
-1.4e+03 9.6e+02 0 0
-1.4e+03 5.0e+02 0 0
-2.0e+03 6.7e+02 0 0
-
-
-1.4e+03 5.0e+02 0 0
-1.4e+03 9.6e+02 0 0
-9.0e+02 7.8e+02 0 0
-9.6e+02 4.0e+02 0 0
-1.4e+03 5.0e+02 0 0
-
-
-9.6e+02 4.0e+02 0 0
-9.0e+02 7.8e+02 0 0
-5.0e+02 6.7e+02 0 0
-5.0e+02 3.3e+02 0 0
-9.6e+02 4.0e+02 0 0
-
-
-2.0e+03 1.3e+03 0 0
-2.0e+03 2.0e+03 0 0
-1.2e+03 1.3e+03 0 0
-1.4e+03 9.6e+02 0 0
-2.0e+03 1.3e+03 0 0
-
-
-1.4e+03 9.6e+02 0 0
-1.2e+03 1.3e+03 0 0
-7.9e+02 1.1e+03 0 0
-9.0e+02 7.8e+02 0 0
-1.4e+03 9.6e+02 0 0
-
-
-9.0e+02 7.8e+02 0 0
-7.9e+02 1.1e+03 0 0
-5.0e+02 1.0e+03 0 0
-5.0e+02 6.7e+02 0 0
-9.0e+02 7.8e+02 0 0
-
-
-2.0e+03 2.0e+03 0 0
-1.3e+03 2.0e+03 0 0
-8.4e+02 1.6e+03 0 0
-1.2e+03 1.3e+03 0 0
-2.0e+03 2.0e+03 0 0
-
-
-1.2e+03 1.3e+03 0 0
-8.4e+02 1.6e+03 0 0
-6.0e+02 1.4e+03 0 0
-7.9e+02 1.1e+03 0 0
-1.2e+03 1.3e+03 0 0
-
-
-7.9e+02 1.1e+03 0 0
-6.0e+02 1.4e+03 0 0
-4.3e+02 1.2e+03 0 0
-5.0e+02 1.0e+03 0 0
-7.9e+02 1.1e+03 0 0
-
-
-1.3e+03 2.0e+03 0 0
-6.7e+02 2.0e+03 0 0
-4.4e+02 1.8e+03 0 0
-8.4e+02 1.6e+03 0 0
-1.3e+03 2.0e+03 0 0
-
-
-8.4e+02 1.6e+03 0 0
-4.4e+02 1.8e+03 0 0
-3.3e+02 1.6e+03 0 0
-6.0e+02 1.4e+03 0 0
-8.4e+02 1.6e+03 0 0
-
-
-6.0e+02 1.4e+03 0 0
-3.3e+02 1.6e+03 0 0
-2.5e+02 1.4e+03 0 0
-4.3e+02 1.2e+03 0 0
-6.0e+02 1.4e+03 0 0
-
-
-6.7e+02 2.0e+03 0 0
-0.0 2.0e+03 0 0
-0.0 1.8e+03 0 0
-4.4e+02 1.8e+03 0 0
-6.7e+02 2.0e+03 0 0
-
-
-4.4e+02 1.8e+03 0 0
-0.0 1.8e+03 0 0
-0.0 1.7e+03 0 0
-3.3e+02 1.6e+03 0 0
-4.4e+02 1.8e+03 0 0
-
-
-3.3e+02 1.6e+03 0 0
-0.0 1.7e+03 0 0
-0.0 1.5e+03 0 0
-2.5e+02 1.4e+03 0 0
-3.3e+02 1.6e+03 0 0
-
-
+DEAL::vertex 0 has boundary indicator 10
+DEAL::vertex 1 has boundary indicator 20
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check whether we can read in with the gmsh format
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_in.h>
+
+#include <fstream>
+#include <cmath>
+
+
+template <int dim>
+void gmsh_grid (const char *name)
+{
+ Triangulation<dim> tria;
+ GridIn<dim> grid_in;
+ grid_in.attach_triangulation (tria);
+ std::ifstream input_file(name);
+ grid_in.read_msh(input_file);
+
+ deallog << " " << tria.n_active_cells() << " active cells" << std::endl;
+
+ int hash = 0;
+ int index = 0;
+ for (typename Triangulation<dim>::active_cell_iterator c=tria.begin_active();
+ c!=tria.end(); ++c, ++index)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ hash += (index * i * c->vertex_index(i)) % (tria.n_active_cells()+1);
+ deallog << " hash=" << hash << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ gmsh_grid<2> (SOURCE_DIR "/grids/grid_in_msh_01.2d.msh");
+ gmsh_grid<2> (SOURCE_DIR "/grids/grid_in_msh_01.2da.msh");
+ gmsh_grid<3> (SOURCE_DIR "/grids/grid_in_msh_01.3d.msh");
+ gmsh_grid<3> (SOURCE_DIR "/grids/grid_in_msh_01.3da.msh");
+ gmsh_grid<3> (SOURCE_DIR "/grids/grid_in_msh_01.3d_neg.msh");
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+
+ return 0;
+}
--- /dev/null
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+11
+1 0 0 0
+2 1 0 0
+3 0.1 0 0
+4 0.2 0 0
+5 0.3 0 0
+6 0.4 0 0
+7 0.5 0 0
+8 0.6 0 0
+9 0.7 0 0
+10 0.8 0 0
+11 0.9 0 0
+$EndNodes
+$Elements
+12
+1 15 2 10 1 1
+2 15 2 20 2 2
+3 1 2 100 1 1 3
+4 1 2 100 1 3 4
+5 1 2 100 1 4 5
+6 1 2 100 1 5 6
+7 1 2 100 1 6 7
+8 1 2 100 1 7 8
+9 1 2 100 1 8 9
+10 1 2 100 1 9 10
+11 1 2 100 1 10 11
+12 1 2 100 1 11 2
+$EndElements
projects / dealii.git / commitdiff