From ee02868ac3194d80e9db86032d26688fef5a34b4 Mon Sep 17 00:00:00 2001
From: amola <amola@terminus.sissa.it>
Date: Tue, 5 Apr 2016 09:50:30 +0200
Subject: [PATCH] added new option to read indicators in a ucd file as manifold
 ids and not as boundary ids

---
 include/deal.II/grid/grid_in.h |  8 +++++++-
 source/grid/grid_in.cc         | 19 ++++++++++++++-----
 source/grid/tria.cc            | 31 +++++++++++++++++++++++++------
 3 files changed, 46 insertions(+), 12 deletions(-)

diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h
index 1177d07955..30bf551332 100644
--- a/include/deal.II/grid/grid_in.h
+++ b/include/deal.II/grid/grid_in.h
@@ -366,8 +366,14 @@ public:
 
   /**
    * Read grid data from an ucd file. Numerical data is ignored.
+   * It is not possible to use a ucd file to set both boundary_id and
+   * manifold_id for the same cell. Yet it is possible to use
+   * the flag apply_all_indicators_to_manifolds to decide if
+   * the indicators in the file refer to manifolds (flag set to true)
+   * or boundaries (flag set to false).
    */
-  void read_ucd (std::istream &in);
+  void read_ucd (std::istream                                  &in,
+                 const bool apply_all_indicators_to_manifolds=false);
 
   /**
    * Read grid data from an Abaqus file. Numerical and constitutive data is
diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index df169e1498..119c196546 100644
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -610,7 +610,8 @@ void GridIn<dim, spacedim>::read_unv(std::istream &in)
 
 
 template <int dim, int spacedim>
-void GridIn<dim, spacedim>::read_ucd (std::istream &in)
+void GridIn<dim, spacedim>::read_ucd (std::istream                            &in,
+                                      const bool apply_all_indicators_to_manifolds)
 {
   Assert (tria != 0, ExcNoTriangulationSelected());
   AssertThrow (in, ExcIO());
@@ -727,8 +728,12 @@ void GridIn<dim, spacedim>::read_ucd (std::istream &in)
           Assert(material_id < numbers::internal_face_boundary_id,
                  ExcIndexRange(material_id,0,numbers::internal_face_boundary_id));
 
-          subcelldata.boundary_lines.back().boundary_id
-            = static_cast<types::boundary_id>(material_id);
+          if (apply_all_indicators_to_manifolds)
+            subcelldata.boundary_lines.back().manifold_id
+              = static_cast<types::manifold_id>(material_id);
+          else
+            subcelldata.boundary_lines.back().boundary_id
+              = static_cast<types::boundary_id>(material_id);
 
           // transform from ucd to
           // consecutive numbering
@@ -764,8 +769,12 @@ void GridIn<dim, spacedim>::read_ucd (std::istream &in)
           Assert(material_id < numbers::internal_face_boundary_id,
                  ExcIndexRange(material_id,0,numbers::internal_face_boundary_id));
 
-          subcelldata.boundary_quads.back().boundary_id
-            = static_cast<types::boundary_id>(material_id);
+          if (apply_all_indicators_to_manifolds)
+            subcelldata.boundary_quads.back().manifold_id
+              = static_cast<types::manifold_id>(material_id);
+          else
+            subcelldata.boundary_quads.back().boundary_id
+              = static_cast<types::boundary_id>(material_id);
 
           // transform from ucd to
           // consecutive numbering
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 7bd94ee3f8..fce3fc8300 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -2068,13 +2068,32 @@ namespace internal
             // need to allow that
             if (boundary_line->boundary_id != numbers::internal_face_boundary_id)
               {
-                AssertThrow (! (line->boundary_id() == numbers::internal_face_boundary_id),
-                             ExcInteriorLineCantBeBoundary(line->vertex_index(0),
-                                                           line->vertex_index(1),
-                                                           boundary_line->boundary_id));
-                line->set_boundary_id (boundary_line->boundary_id);
+                if (line->boundary_id() == numbers::internal_face_boundary_id)
+                  {
+                    // if we are here, it means that we want to assign a boundary indicator
+                    // different from numbers::internal_face_boundary_id to an internal line.
+                    // As said, this would be not allowed, and an exception should be immediately
+                    // thrown. Still, there is the possibility that one only wants to specifiy a
+                    // manifold_id here. If that is the case (manifold_id !=  numbers::flat_manifold_id)
+                    // the operation is allowed. Otherwise, we really tried to specify a boundary_id
+                    // (and not a manifold_id) to an internal face. The exception must be thrown.
+                    if (boundary_line->manifold_id == numbers::flat_manifold_id)
+                      {
+                        // If we are here, this assertion will surely fail, for the aforementioned
+                        // reasons
+                        AssertThrow (! (line->boundary_id() == numbers::internal_face_boundary_id),
+                                     ExcInteriorLineCantBeBoundary(line->vertex_index(0),
+                                                                   line->vertex_index(1),
+                                                                   boundary_line->boundary_id));
+                      }
+                    else
+                      {
+                        line->set_manifold_id (boundary_line->manifold_id);
+                      }
+                  }
+                else
+                  line->set_boundary_id (boundary_line->boundary_id);
               }
-
             line->set_manifold_id (boundary_line->manifold_id);
           }
 
-- 
2.39.5