]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Backward compatible manifold_id branch
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Apr 2013 18:18:41 +0000 (18:18 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Apr 2013 18:18:41 +0000 (18:18 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@29311 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/base/types.h
deal.II/include/deal.II/grid/tria.h
deal.II/include/deal.II/grid/tria_accessor.h
deal.II/include/deal.II/grid/tria_accessor.templates.h
deal.II/include/deal.II/grid/tria_objects.h
deal.II/source/grid/tria.cc
deal.II/source/grid/tria.inst.in
deal.II/source/grid/tria_accessor.cc
deal.II/source/grid/tria_objects.cc

index 2a37e62c29b90cd115371dd10cb5191525de3bf1..ec54568f5e4f2bbd9fc7ad4558008d4fd32248df 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2009, 2012 by the deal.II authors
+//    Copyright (C) 2009, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -77,7 +77,19 @@ namespace types
    * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
   typedef unsigned char boundary_id;
+  
+  /**
+   * The type used to denote manifold indicators associated with every
+   * object of the mesh.
+   *
+   * There is a special value, numbers::flat_manifold_id
+   * that is used to indicate the standard cartesian manifold.
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  typedef unsigned char manifold_id;
 
+  
   /**
    * @deprecated Old name for the typedef above.
    */
@@ -164,6 +176,24 @@ namespace numbers
    */
   const types::boundary_id internal_face_boundary_id = static_cast<types::boundary_id>(-1);
 
+  /**
+   * Invalid manifold_id which we
+   * need in several places as a
+   * default value.  We assume that
+   * all valid manifold_ids lie in the
+   * range [0, invalid_maifold_id).
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  const types::manifold_id invalid_manifold_id = static_cast<types::manifold_id>(-1);
+  
+  /**
+   * A manifold_id we reserve for the default flat Cartesian manifold.
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  const types::manifold_id flat_manifold_id = static_cast<types::manifold_id>(-1);
+
   /**
    * A special id for an invalid
    * subdomain id. This value may not
index 1962bd3a044a60cc3eac6153ba9dc84e340142cc..ce1fa0c4a81c117b52b3b5dbd18e8e15d3a61908 100644 (file)
@@ -92,6 +92,11 @@ template <int dim, int spacedim> class MGDoFHandler;
 template <int structdim>
 struct CellData
 {
+                                    /** Default constructor. Set
+                                        members to meaningful default
+                                        values. */
+    CellData();
+
   /**
    * Indices of the vertices of
    * this cell.
@@ -115,6 +120,16 @@ struct CellData
     types::boundary_id boundary_id;
     types::material_id material_id;
   };
+                                    /**
+                                     * Manifold identificator of this
+                                     * object. This identificator
+                                     * should be used to identify the
+                                     * manifold to which this object
+                                     * belongs, and from which this
+                                     * object will collect information
+                                     * on how to add points upon
+                                     * refinement.*/
+  types::manifold_id manifold_id;
 };
 
 
@@ -1785,7 +1800,7 @@ public:
    *
    * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  void set_boundary (const types::boundary_id   number,
+  void set_boundary (const types::manifold_id   number,
                      const Boundary<dim,spacedim> &boundary_object);
 
   /**
@@ -1801,7 +1816,7 @@ public:
    *
    * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  void set_boundary (const types::boundary_id number);
+  void set_boundary (const types::manifold_id number);
 
   /**
    * Return a constant reference to
@@ -1814,7 +1829,7 @@ public:
    *
    * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  const Boundary<dim,spacedim> &get_boundary (const types::boundary_id number) const;
+  const Boundary<dim,spacedim> &get_boundary (const types::manifold_id number) const;
 
   /**
    * Returns a vector containing
@@ -1834,7 +1849,25 @@ public:
    * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
   std::vector<types::boundary_id> get_boundary_indicators() const;
-
+    
+  /**
+   * Returns a vector containing
+   * all manifold indicators
+   * assigned to the objects of
+   * this Triangulation. Note, that each
+   * manifold indicator is reported
+   * only once. The size of the
+   * return vector will represent
+   * the number of different
+   * indicators (which is greater
+   * or equal one).
+   *
+   * @ingroup manifold
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  std::vector<types::manifold_id> get_manifold_ids() const;
+    
   /**
    *  Copy a triangulation. This
    *  operation is not cheap, so
@@ -3544,7 +3577,7 @@ private:
    *  objects, which are not
    *  of type StraightBoundary.
    */
-  std::map<types::boundary_id, SmartPointer<const Boundary<dim, spacedim> , Triangulation<dim, spacedim> > >  boundary;
+  std::map<types::manifold_id, SmartPointer<const Boundary<dim, spacedim> , Triangulation<dim, spacedim> > >  boundary;
 
   /**
    * Flag indicating whether
@@ -3604,6 +3637,34 @@ private:
    * were not a pointer.
    */
   std::map<unsigned int, types::boundary_id> *vertex_to_boundary_id_map_1d;
+       
+
+       /**
+        * A map that relates the number of a
+        * boundary vertex to the manifold
+        * indicator. This field is only used in
+        * 1d. We have this field because we
+        * store manifold indicator information
+        * with faces in 2d and higher where we
+        * have space in the structures that
+        * store data for faces, but in 1d there
+        * is no such space for faces.
+        *
+        * The field is declared as a pointer for
+        * a rather mundane reason: all other
+        * fields of this class that can be
+        * modified by the TriaAccessor hierarchy
+        * are pointers, and so these accessor
+        * classes store a const pointer to the
+        * triangulation. We could no longer do
+        * so for TriaAccessor<0,1,spacedim> if
+        * this field (that can be modified by
+        * TriaAccessor::set_boundary_indicator)
+        * were not a pointer.
+        */
+       
+  std::map<unsigned int, types::manifold_id> *vertex_to_manifold_id_map_1d;
+
 
   /**
    * A map that correlates each refinement listener that has been added
index f4e2167f3f84a48f4bfadb12afc7a9194d5591b7..d8bb9d2833a0e14e0bbab76caa23fa7c8430f9fe 100644 (file)
@@ -1156,6 +1156,70 @@ public:
    * @}
    */
 
+  /**
+   *  @name Dealing with manifold indicators
+   */
+  /**
+   * @{
+   */
+
+  /**
+   * Manifold indicator of this
+   * object.
+   *
+   * If the return value is the special value
+   * numbers::flat_manifold_id, then this object is associated with a
+   * standard cartesian Manifold Description.
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  types::manifold_id manifold_id () const;
+
+  /**
+   * Set the manifold indicator.  The same applies as for the
+   * <tt>manifold_id()</tt> function.
+   *
+   * Note that it only sets the manifold object of the current object
+   * itself, not the indicators of the ones that bound it, nor of its
+   * children. For example, in 3d, if this function is called on a
+   * face, then the manifold indicator of the 4 edges that bound the
+   * face remain unchanged. If you want to set the manifold indicators
+   * of face, edges and all children at the same time, use the
+   * set_all_manifold_ids() function.
+   *
+   *
+   * @ingroup manifold
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  void set_manifold_id (const types::manifold_id) const;
+
+  /**
+   * Do as set_manifold_id()
+   * but also set the manifold
+   * indicators of the objects that
+   * bound the current object. For
+   * example, in 3d, if
+   * set_manifold_id() is
+   * called on a face, then the
+   * manifold indicator of the 4
+   * edges that bound the face
+   * remain unchanged. On the other
+   * hand, the manifold indicators
+   * of face and edges are all set
+   * at the same time using the
+   * current function.
+   *
+   * @ingroup manifold
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  void set_all_manifold_ids (const types::manifold_id) const;
+
+  /**
+   * @}
+   */    
+    
 
   /**
    * @name User data
@@ -2003,6 +2067,15 @@ public:
    */
   types::boundary_id boundary_indicator () const;
 
+  /**
+   * Manifold indicator of this
+   * object. 
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+  types::manifold_id manifold_id () const;
+
+
   /**
    *  @name Orientation of sub-objects
    */
@@ -2143,6 +2216,16 @@ public:
   void
   set_boundary_indicator (const types::boundary_id);
 
+                                    /**
+                                     * Set the manifold indicator of
+                                     * this vertex. Does nothing so
+                                     * far. This function is here
+                                     * only to allow dimension
+                                     * independent programming.
+                                     */
+  void
+  set_manifold_id (const types::manifold_id);
+
   /**
    * Since this object only represents a
    * single vertex, call
@@ -2158,6 +2241,22 @@ public:
   /**
    * @}
    */
+/**
+   * Since this object only represents a
+   * single vertex, call
+   * set_manifold_id with the same
+   * argument.
+   *
+   * @ingroup manifold
+   *
+   * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
+   */
+    
+    void
+    set_all_manifold_ids (const types::manifold_id);
+    /**
+   * @}
+   */
 
   /**
    * Return whether the vertex
index 5379e3a5af5ffe0f18f2777210bfca1941f30bf5..bb1b20ec80302e66be9c060008b02bf83a080a4d 100644 (file)
@@ -1918,11 +1918,79 @@ template <int structdim, int dim, int spacedim>
 const Boundary<dim,spacedim> &
 TriaAccessor<structdim, dim, spacedim>::get_boundary () const
 {
-  Assert (structdim<dim, ExcImpossibleInDim(dim));
+                                  // This assert should be removed in
+                                  // future releases...
+  Assert (structdim<spacedim, ExcImpossibleInDim(dim));
+  Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
+                                  // Get the default (manifold_id)
+  types::manifold_id mi = this->objects().manifold_id[this->present_index];
+
+                                  // If the manifold_id is set to
+                                  // invalid, then ask for a
+                                  // boundary_id
+  if(mi == numbers::invalid_manifold_id)
+    mi =  static_cast<types::manifold_id>
+         (structdim < dim ?
+          this->objects().boundary_or_material_id[this->present_index].boundary_id:
+          this->objects().boundary_or_material_id[this->present_index].material_id);
+  
+  return this->tria->get_boundary(mi);
+}
+
+
+template <int structdim, int dim, int spacedim>
+types::manifold_id
+TriaAccessor<structdim, dim, spacedim>::manifold_id () const
+{
+  Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
+
+  return this->objects().manifold_id[this->present_index];
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+void
+TriaAccessor<structdim, dim, spacedim>::
+set_manifold_id (const types::manifold_id manifold_ind) const
+{
   Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
 
-  return this->tria->get_boundary(this->objects()
-                                  .boundary_or_material_id[this->present_index].boundary_id);
+  this->objects().manifold_id[this->present_index] = manifold_ind;
+}
+
+
+template <int structdim, int dim, int spacedim>
+void
+TriaAccessor<structdim, dim, spacedim>::
+set_all_manifold_ids (const types::manifold_id manifold_ind) const
+{
+  set_manifold_id (manifold_ind);
+
+  if (this->has_children())
+    for (unsigned int c=0; c<this->n_children(); ++c)
+      this->child(c)->set_all_manifold_ids (manifold_ind);
+
+  switch (structdim)
+    {
+    case 1:
+         if(dim == 1)
+           {
+             (*this->tria->vertex_to_manifold_id_map_1d)
+               [vertex_index(0)] = manifold_ind;
+             (*this->tria->vertex_to_manifold_id_map_1d)
+               [vertex_index(1)] = manifold_ind;
+           }
+      break;
+
+    case 2:
+      // for quads also set manifold_id of bounding lines
+      for (unsigned int i=0; i<4; ++i)
+        this->line(i)->set_manifold_id (manifold_ind);
+      break;
+    default:
+      Assert (false, ExcNotImplemented());
+    }
 }
 
 
@@ -2294,6 +2362,18 @@ TriaAccessor<0, 1, spacedim>::boundary_indicator () const
 
 }
 
+template <int spacedim>
+inline
+types::manifold_id
+TriaAccessor<0, 1, spacedim>::manifold_id () const
+{
+  if( tria->vertex_to_manifold_id_map_1d->find (this->vertex_index())
+      != tria->vertex_to_boundary_id_map_1d->end())
+    return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()];
+  else
+    return numbers::invalid_manifold_id;
+}
+
 
 template <int spacedim>
 inline
@@ -2420,6 +2500,15 @@ TriaAccessor<0, 1, spacedim>::set_boundary_indicator (const types::boundary_id b
 }
 
 
+template <int spacedim>
+inline
+void
+TriaAccessor<0, 1, spacedim>::set_manifold_id (const types::manifold_id b)
+{
+  (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()] = b;
+}
+
+
 
 template <int spacedim>
 inline
@@ -2430,6 +2519,15 @@ void TriaAccessor<0, 1, spacedim>::set_all_boundary_indicators (const types::bou
 
 
 
+template <int spacedim>
+inline
+void TriaAccessor<0, 1, spacedim>::set_all_manifold_ids (const types::manifold_id b)
+{
+  set_manifold_id (b);
+}
+
+
+
 template <int spacedim>
 inline
 bool TriaAccessor<0, 1, spacedim>::used () const
index c90e306a628979de6e7311db09a6b01ec1190b8f..6c58fa78f6be0a756fa9f9146b483468fd46d04a 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -195,6 +195,14 @@ namespace internal
        */
       std::vector<BoundaryOrMaterialId> boundary_or_material_id;
 
+      /**
+       * Store manifold ids. This field
+       * stores the manifold id of each object, which
+       * is a number between 0 and
+       * numbers::invalid_manifold_id-1. 
+       */
+      std::vector<types::manifold_id> manifold_id;
+
       /**
        *  Assert that enough space
        *  is allocated to
@@ -879,6 +887,7 @@ namespace internal
       ar &used;
       ar &user_flags;
       ar &boundary_or_material_id;
+      ar &manifold_id;
       ar &next_free_single &next_free_pair &reverse_order_next_free_single;
       ar &user_data &user_data_type;
     }
index 4350049bccf5797eae5f9e33afe6816d3527f6d4..02d7e5d2bb54a969a1d8ed46a2498a8147ebd8db 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+template<int dim>
+CellData<dim>::CellData()
+{
+                                  // Set the face to be interior
+  boundary_id = numbers::internal_face_boundary_id;
+                                  // And the manifold to be invalid
+  manifold_id = numbers::invalid_manifold_id;
+}
+
+
 bool
 SubCellData::check_consistency (const unsigned int dim) const
 {
@@ -1483,6 +1493,7 @@ namespace internal
                                                   cells[cell].vertices[1]));
             next_free_line->set_used_flag ();
             next_free_line->set_material_id (cells[cell].material_id);
+           next_free_line->set_manifold_id (cells[cell].manifold_id);
             next_free_line->clear_user_data ();
             next_free_line->set_subdomain_id (0);
 
@@ -1570,16 +1581,24 @@ namespace internal
 
         // finally set the
         // vertex_to_boundary_id_map_1d
-        // map
+       // and vertex_to_manifold_id_map_1d
+        // maps
         triangulation.vertex_to_boundary_id_map_1d->clear();
+       triangulation.vertex_to_manifold_id_map_1d->clear();
         for (typename Triangulation<dim,spacedim>::active_cell_iterator
              cell = triangulation.begin_active();
              cell != triangulation.end(); ++cell)
           for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-            if (cell->at_boundary(f))
-              (*triangulation
-               .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()]
-                = f;
+           {
+             (*triangulation
+              .vertex_to_manifold_id_map_1d)[cell->face(f)->vertex_index()]
+               = numbers::invalid_manifold_id;
+             
+             if (cell->at_boundary(f))
+               (*triangulation
+                .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()]
+                 = f;
+           }
       }
 
 
@@ -1756,6 +1775,7 @@ namespace internal
 
               cell->set_used_flag ();
               cell->set_material_id (cells[c].material_id);
+             cell->set_manifold_id (cells[c].manifold_id);
               cell->clear_user_data ();
               cell->set_subdomain_id (0);
 
@@ -1788,6 +1808,7 @@ namespace internal
             else
               // interior line -> numbers::internal_face_boundary_id
               line->set_boundary_indicator (numbers::internal_face_boundary_id);
+           line->set_manifold_id(numbers::flat_manifold_id);
           }
 
         // set boundary indicators where
@@ -1831,6 +1852,7 @@ namespace internal
                          ExcInteriorLineCantBeBoundary());
 
             line->set_boundary_indicator (boundary_line->boundary_id);
+           line->set_manifold_id (boundary_line->manifold_id);
           }
 
 
@@ -2407,6 +2429,7 @@ namespace internal
 
               cell->set_used_flag ();
               cell->set_material_id (cells[c].material_id);
+             cell->set_manifold_id (cells[c].manifold_id);
               cell->clear_user_flag ();
               cell->clear_user_data ();
               cell->set_subdomain_id (0);
@@ -2508,6 +2531,10 @@ namespace internal
             else
               // interior quad -> numbers::internal_face_boundary_id
               quad->set_boundary_indicator (numbers::internal_face_boundary_id);
+                                            // Manifold ids are set
+                                            // independently of where
+                                            // they are
+           quad->set_manifold_id(numbers::flat_manifold_id);
           }
 
         /////////////////////////////////////////
@@ -2515,11 +2542,15 @@ namespace internal
         // the boundary and mark all others as
         // interior ones
         //
-        // for this: first mark all lines
-        // as interior
+        // for this: first mark all lines as interior. use this loop
+        // to also set all manifold ids of all lines
         for (typename Triangulation<dim,spacedim>::line_iterator
-             line=triangulation.begin_line(); line!=triangulation.end_line(); ++line)
-          line->set_boundary_indicator (numbers::internal_face_boundary_id);
+              line=triangulation.begin_line(); line!=triangulation.end_line(); ++line)
+         {
+           line->set_boundary_indicator (numbers::internal_face_boundary_id);
+           line->set_manifold_id(numbers::flat_manifold_id);
+         }
+       
         // next reset all lines bounding
         // boundary quads as on the
         // boundary also. note that since
@@ -2592,6 +2623,8 @@ namespace internal
                                        "if they carry the same boundary indicator."));
 
             line->set_boundary_indicator (boundary_line->boundary_id);
+                                            // Set manifold id if given
+           line->set_manifold_id(boundary_line->manifold_id);
           }
 
 
@@ -2724,6 +2757,7 @@ namespace internal
                                        "if they carry the same boundary indicator."));
 
             quad->set_boundary_indicator (boundary_quad->boundary_id);
+           quad->set_manifold_id (boundary_quad->manifold_id);
           }
 
 
@@ -3577,6 +3611,7 @@ namespace internal
                           switch_1->set_line_orientation(2, switch_2->line_orientation(2));
                           switch_1->set_line_orientation(3, switch_2->line_orientation(3));
                           switch_1->set_boundary_indicator(switch_2->boundary_indicator());
+                         switch_1->set_manifold_id(switch_2->manifold_id());
                           switch_1->set_user_index(switch_2->user_index());
                           if (switch_2->user_flag_set())
                             switch_1->set_user_flag();
@@ -3592,6 +3627,7 @@ namespace internal
                           switch_2->set_line_orientation(2, switch_1_line_orientations[2]);
                           switch_2->set_line_orientation(3, switch_1_line_orientations[3]);
                           switch_2->set_boundary_indicator(switch_1_boundary_indicator);
+                         switch_2->set_manifold_id(switch_1->manifold_id());
                           switch_2->set_user_index(switch_1_user_index);
                           if (switch_1_user_flag)
                             switch_2->set_user_flag();
@@ -4011,13 +4047,6 @@ namespace internal
                 cell->clear_user_flag();
 
 
-
-                // An assert to make sure that the
-                // static_cast in the next line has
-                // the chance to give reasonable results.
-                Assert(cell->material_id()<= std::numeric_limits<types::material_id>::max(),
-                       ExcIndexRange(cell->material_id(),0,std::numeric_limits<types::material_id>::max()));
-
                 // new vertex is
                 // placed on the
                 // surface according
@@ -4025,8 +4054,7 @@ namespace internal
                 // stored in the
                 // boundary class
                 triangulation.vertices[next_unused_vertex] =
-                  triangulation.get_boundary(static_cast<types::boundary_id>(cell->material_id()))
-                  .get_new_point_on_quad (cell);
+                 cell->get_boundary().get_new_point_on_quad (cell);
               }
           }
 
@@ -4121,6 +4149,7 @@ namespace internal
             new_lines[l]->clear_children();
             // interior line
             new_lines[l]->set_boundary_indicator(numbers::internal_face_boundary_id);
+           new_lines[l]->set_manifold_id(cell->manifold_id());
           }
 
         // Now add the four (two)
@@ -4241,6 +4270,7 @@ namespace internal
             // inherit material
             // properties
             subcells[i]->set_material_id (cell->material_id());
+           subcells[i]->set_manifold_id (cell->manifold_id());
             subcells[i]->set_subdomain_id (cell->subdomain_id());
 
             if (i%2==0)
@@ -4403,19 +4433,19 @@ namespace internal
                   Assert (next_unused_vertex < triangulation.vertices.size(),
                           ExcTooFewVerticesAllocated());
 
-                  // first insert new
-                  // vertex. if dim==spacedim
-                  // then simply use the
-                  // midpoint; otherwise we
-                  // have to ask the manifold
-                  // object
-                  if (dim == spacedim)
-                    triangulation.vertices[next_unused_vertex] =
-                      (cell->vertex(0) + cell->vertex(1)) / 2;
-                  else
-                    triangulation.vertices[next_unused_vertex] =
-                      triangulation.get_boundary(static_cast<types::boundary_id>(cell->material_id()))
-                      .get_new_point_on_line(cell);
+                   // first insert new
+                   // vertex. if dim==spacedim
+                   // then simply use the
+                   // midpoint; otherwise we
+                   // have to ask the manifold
+                   // object
+                 if (dim == spacedim)
+                   triangulation.vertices[next_unused_vertex] =
+                     (cell->vertex(0) + cell->vertex(1)) / 2;
+                 else
+                   triangulation.vertices[next_unused_vertex] =
+                     cell->get_boundary().get_new_point_on_line(cell);
+                 
                   triangulation.vertices_used[next_unused_vertex] = true;
 
                   // search for next two
@@ -4444,11 +4474,19 @@ namespace internal
                                     ::TriaObject<1> (cell->vertex_index(0),
                                                      next_unused_vertex));
                   first_child->set_material_id (cell->material_id());
+                 first_child->set_manifold_id (cell->manifold_id());
                   first_child->set_subdomain_id (cell->subdomain_id());
                   first_child->set_direction_flag (cell->direction_flag());
-
+                 
                   first_child->set_parent (cell->index ());
 
+                                                  // Set manifold id
+                                                  // of the right
+                                                  // face. Only do
+                                                  // this on the
+                                                  // first child.
+                 first_child->face(1)->set_manifold_id(cell->manifold_id());
+                 
                   // reset neighborship info (refer
                   // to
                   // internal::Triangulation::TriaLevel<0>
@@ -4504,6 +4542,7 @@ namespace internal
                                                      cell->vertex_index(1)));
                   second_child->set_neighbor (0, first_child);
                   second_child->set_material_id (cell->material_id());
+                 second_child->set_manifold_id (cell->manifold_id());
                   second_child->set_subdomain_id (cell->subdomain_id());
                   second_child->set_direction_flag (cell->direction_flag());
 
@@ -4813,8 +4852,15 @@ namespace internal
                           ExcTooFewVerticesAllocated());
                   triangulation.vertices_used[next_unused_vertex] = true;
 
-                  if (spacedim == dim)
-                    {
+                                                  // Ask the manifold
+                                                  // where to put the
+                                                  // new line
+                 triangulation.vertices[next_unused_vertex]
+                   = line->get_boundary().get_new_point_on_line(line);
+
+                 
+//                   if (spacedim == dim)
+//                     {
                       // for the case of a domain
                       // in an equal-dimensional
                       // space we only have to
@@ -4823,21 +4869,22 @@ namespace internal
                       // lines we can compute the
                       // midpoint as the mean of
                       // the two vertices:
-                      if (line->at_boundary())
-                        triangulation.vertices[next_unused_vertex]
-                          = triangulation.get_boundary(line->boundary_indicator())
-                            .get_new_point_on_line (line);
-                      else
-                        triangulation.vertices[next_unused_vertex]
-                          = (line->vertex(0) + line->vertex(1)) / 2;
-                    }
-                  else
+//                      if (line->at_boundary())
+//                         triangulation.vertices[next_unused_vertex]
+//                           = line->get_boundary().get_new_point_on_line(line);
+//                       triangulation.get_boundary(line->boundary_indicator())
+//                             .get_new_point_on_line (line);
+//                       else
+//                         triangulation.vertices[next_unused_vertex]
+//                          = (line->vertex(0) + line->vertex(1)) / 2;
+//                     }
+//                   else
                     // however, if spacedim>dim, we
                     // always have to ask the
                     // boundary object for its
-                    // answer
-                    triangulation.vertices[next_unused_vertex]
-                      = triangulation.get_boundary(line->user_index()).get_new_point_on_line (line);
+//                     // answer
+//                     triangulation.vertices[next_unused_vertex]
+//                       = triangulation.get_boundary(line->user_index()).get_new_point_on_line (line);
 
                   // now that we created
                   // the right point, make
@@ -4901,6 +4948,9 @@ namespace internal
                   children[0]->set_boundary_indicator (line->boundary_indicator());
                   children[1]->set_boundary_indicator (line->boundary_indicator());
 
+                  children[0]->set_manifold_id (line->manifold_id());
+                  children[1]->set_manifold_id (line->manifold_id());
+
                   // finally clear flag
                   // indicating the need
                   // for refinement
@@ -5366,12 +5416,15 @@ namespace internal
                           ExcTooFewVerticesAllocated());
                   triangulation.vertices_used[next_unused_vertex] = true;
 
-                  if (line->at_boundary())
-                    triangulation.vertices[next_unused_vertex]
-                      = triangulation.get_boundary(line->boundary_indicator()).get_new_point_on_line (line);
-                  else
-                    triangulation.vertices[next_unused_vertex]
-                      = (line->vertex(0) + line->vertex(1)) / 2;
+                 triangulation.vertices[next_unused_vertex]
+                   = line->get_boundary().get_new_point_on_line(line);
+
+//                   if (line->at_boundary())
+//                     triangulation.vertices[next_unused_vertex]
+//                       = triangulation.get_boundary(line->boundary_indicator()).get_new_point_on_line (line);
+//                   else
+//                     triangulation.vertices[next_unused_vertex]
+//                       = (line->vertex(0) + line->vertex(1)) / 2;
 
                   // now that we created
                   // the right point, make
@@ -5424,6 +5477,9 @@ namespace internal
 
                   children[0]->set_boundary_indicator (line->boundary_indicator());
                   children[1]->set_boundary_indicator (line->boundary_indicator());
+                 
+                 children[0]->set_manifold_id (line->manifold_id());
+                  children[1]->set_manifold_id (line->manifold_id());
 
                   // finally clear flag
                   // indicating the need
@@ -5547,6 +5603,7 @@ namespace internal
                     new_line->clear_user_data();
                     new_line->clear_children();
                     new_line->set_boundary_indicator(quad->boundary_indicator());
+                   new_line->set_manifold_id(quad->manifold_id());
 
                     // child 0 and 1 of a line are
                     // switched if the line
@@ -5611,6 +5668,7 @@ namespace internal
                         new_quads[i]->clear_user_data();
                         new_quads[i]->clear_children();
                         new_quads[i]->set_boundary_indicator (quad->boundary_indicator());
+                       new_quads[i]->set_manifold_id (quad->manifold_id());
                         // set all line orientations to
                         // true, change this after the
                         // loop, as we have to consider
@@ -5731,6 +5789,7 @@ namespace internal
                                 new_child[i]->set(internal::Triangulation::TriaObject<1>(old_child[i]->vertex_index(0),
                                                                                          old_child[i]->vertex_index(1)));
                                 new_child[i]->set_boundary_indicator(old_child[i]->boundary_indicator());
+                               new_child[i]->set_manifold_id(old_child[i]->manifold_id());
                                 new_child[i]->set_user_index(old_child[i]->user_index());
                                 if (old_child[i]->user_flag_set())
                                   new_child[i]->set_user_flag();
@@ -5845,6 +5904,7 @@ namespace internal
                             switch_1->set_line_orientation(2, switch_2->line_orientation(2));
                             switch_1->set_line_orientation(3, switch_2->line_orientation(3));
                             switch_1->set_boundary_indicator(switch_2->boundary_indicator());
+                           switch_1->set_manifold_id(switch_2->manifold_id());
                             switch_1->set_user_index(switch_2->user_index());
                             if (switch_2->user_flag_set())
                               switch_1->set_user_flag();
@@ -5867,6 +5927,7 @@ namespace internal
                             switch_2->set_line_orientation(2, switch_1_line_orientations[2]);
                             switch_2->set_line_orientation(3, switch_1_line_orientations[3]);
                             switch_2->set_boundary_indicator(switch_1_boundary_indicator);
+                           switch_2->set_manifold_id(switch_1->manifold_id());
                             switch_2->set_user_index(switch_1_user_index);
                             if (switch_1_user_flag)
                               switch_2->set_user_flag();
@@ -5966,7 +6027,8 @@ namespace internal
                             // care about boundary quads
                             // here
                             triangulation.vertices[next_unused_vertex]
-                              = (middle_line->vertex(0) + middle_line->vertex(1)) / 2;
+                             = middle_line->get_boundary().get_new_point_on_line(middle_line);
+//                              = (middle_line->vertex(0) + middle_line->vertex(1)) / 2;
                             triangulation.vertices_used[next_unused_vertex] = true;
 
                             // now search a slot for the two
@@ -6010,6 +6072,9 @@ namespace internal
 
                             children[0]->set_boundary_indicator (middle_line->boundary_indicator());
                             children[1]->set_boundary_indicator (middle_line->boundary_indicator());
+
+                            children[0]->set_manifold_id (middle_line->manifold_id());
+                            children[1]->set_manifold_id (middle_line->manifold_id());
                           }
                         // now remove the flag from the
                         // quad and go to the next
@@ -6029,9 +6094,10 @@ namespace internal
 
                     // set the middle vertex
                     // appropriately
-                    if (quad->at_boundary())
+                    if (quad->at_boundary() || quad->manifold_id()!= numbers::invalid_manifold_id)
                       triangulation.vertices[next_unused_vertex]
-                        = triangulation.get_boundary(quad->boundary_indicator()).get_new_point_on_quad (quad);
+                       = quad->get_boundary().get_new_point_on_quad(quad);
+//                        = triangulation.get_boundary(quad->boundary_indicator()).get_new_point_on_quad (quad);
                     else
                       // it might be that the
                       // quad itself is not
@@ -6078,14 +6144,14 @@ namespace internal
                       // over from the
                       // @p{MappingQ::set_laplace_on_vector}
                       // function
-                      triangulation.vertices[next_unused_vertex]
-                        = (quad->vertex(0) + quad->vertex(1) +
-                           quad->vertex(2) + quad->vertex(3) +
-                           3*(quad->line(0)->child(0)->vertex(1) +
-                              quad->line(1)->child(0)->vertex(1) +
-                              quad->line(2)->child(0)->vertex(1) +
-                              quad->line(3)->child(0)->vertex(1))   ) / 16;
-
+                     triangulation.vertices[next_unused_vertex]
+                       = (quad->vertex(0) + quad->vertex(1) +
+                          quad->vertex(2) + quad->vertex(3) +
+                          3*(quad->line(0)->child(0)->vertex(1) +
+                             quad->line(1)->child(0)->vertex(1) +
+                             quad->line(2)->child(0)->vertex(1) +
+                             quad->line(3)->child(0)->vertex(1))   ) / 16;
+                   
                     triangulation.vertices_used[next_unused_vertex] = true;
 
                     // now that we created
@@ -6159,6 +6225,7 @@ namespace internal
                         new_lines[i]->clear_user_data();
                         new_lines[i]->clear_children();
                         new_lines[i]->set_boundary_indicator(quad->boundary_indicator());
+                       new_lines[i]->set_manifold_id(quad->manifold_id());
                       }
 
                     // now for the
@@ -6264,6 +6331,7 @@ namespace internal
                         new_quads[i]->clear_user_data();
                         new_quads[i]->clear_children();
                         new_quads[i]->set_boundary_indicator (quad->boundary_indicator());
+                       new_quads[i]->set_manifold_id (quad->manifold_id());
                         // set all line orientations to
                         // true, change this after the
                         // loop, as we have to consider
@@ -6390,6 +6458,7 @@ namespace internal
                       new_lines[i]->clear_children();
                       // interior line
                       new_lines[i]->set_boundary_indicator(numbers::internal_face_boundary_id);
+                      new_lines[i]->set_manifold_id(numbers::flat_manifold_id);
                     }
 
                   // find some space for the newly to
@@ -6409,6 +6478,7 @@ namespace internal
                       new_quads[i]->clear_children();
                       // interior quad
                       new_quads[i]->set_boundary_indicator (numbers::internal_face_boundary_id);
+                     new_quads[i]->set_manifold_id (numbers::flat_manifold_id);
                       // set all line orientation
                       // flags to true by default,
                       // change this afterwards, if
@@ -6440,6 +6510,7 @@ namespace internal
                       // inherit material
                       // properties
                       new_hexes[i]->set_material_id (hex->material_id());
+                     new_hexes[i]->set_manifold_id (hex->manifold_id());
                       new_hexes[i]->set_subdomain_id (hex->subdomain_id());
 
                       if (i%2)
@@ -9043,8 +9114,9 @@ namespace internal
                         // boundary would be
                         // this one.
                         const Point<spacedim> new_bound
-                          = triangulation.get_boundary(face->boundary_indicator())
-                            .get_new_point_on_face (face);
+                          = face->get_boundary().get_new_point_on_face(face);
+//                        triangulation.get_boundary(face->boundary_indicator())
+//                             .get_new_point_on_face (face);
                         // to check it,
                         // transform to the
                         // unit cell with
@@ -9434,11 +9506,16 @@ Triangulation (const MeshSmoothing smooth_grid,
   faces(NULL),
   anisotropic_refinement(false),
   check_for_distorted_cells(check_for_distorted_cells),
-  vertex_to_boundary_id_map_1d (0)
+  vertex_to_boundary_id_map_1d (0),
+  vertex_to_manifold_id_map_1d (0)
 {
   if (dim == 1)
-    vertex_to_boundary_id_map_1d
-      = new std::map<unsigned int, types::boundary_id>();
+    {
+      vertex_to_boundary_id_map_1d
+       = new std::map<unsigned int, types::boundary_id>();
+      vertex_to_manifold_id_map_1d
+       = new std::map<unsigned int, types::manifold_id>();
+    }
 
   // connect the any_change signal to the other signals
   signals.create.connect (signals.any_change);
@@ -9456,7 +9533,8 @@ Triangulation (const Triangulation<dim, spacedim> &other)
   :
   Subscriptor(),
   check_for_distorted_cells(other.check_for_distorted_cells),
-  vertex_to_boundary_id_map_1d (0)
+  vertex_to_boundary_id_map_1d (0),
+  vertex_to_manifold_id_map_1d (0)
 {
   Assert (false, ExcInternalError());
 }
@@ -9478,6 +9556,13 @@ Triangulation<dim, spacedim>::~Triangulation ()
           (vertex_to_boundary_id_map_1d == 0),
           ExcInternalError());
   delete vertex_to_boundary_id_map_1d;
+  // the vertex_to_manifold_id_map_1d field
+  // should be unused except in 1d
+  Assert ((dim == 1)
+          ||
+          (vertex_to_manifold_id_map_1d == 0),
+          ExcInternalError());
+  delete vertex_to_manifold_id_map_1d;
 }
 
 
@@ -9503,38 +9588,41 @@ Triangulation<dim, spacedim>::set_mesh_smoothing(const MeshSmoothing mesh_smooth
 
 template <int dim, int spacedim>
 void
-Triangulation<dim, spacedim>::set_boundary (const types::boundary_id number,
+Triangulation<dim, spacedim>::set_boundary (const types::manifold_id m_number,
                                             const Boundary<dim, spacedim> &boundary_object)
 {
+  types::boundary_id number = static_cast<types::boundary_id>(m_number);
   Assert(number < numbers::internal_face_boundary_id, ExcIndexRange(number,0,numbers::internal_face_boundary_id));
 
-  boundary[number] = &boundary_object;
+  boundary[m_number] = &boundary_object;
 }
 
 
 
 template <int dim, int spacedim>
 void
-Triangulation<dim, spacedim>::set_boundary (const types::boundary_id number)
+Triangulation<dim, spacedim>::set_boundary (const types::manifold_id m_number)
 {
+  types::boundary_id number = static_cast<types::boundary_id>(m_number);
   Assert(number < numbers::internal_face_boundary_id, ExcIndexRange(number,0,numbers::internal_face_boundary_id));
 
   //delete the entry located at number.
-  boundary.erase(number);
+  boundary.erase(m_number);
 }
 
 
 
 template <int dim, int spacedim>
 const Boundary<dim, spacedim> &
-Triangulation<dim, spacedim>::get_boundary (const types::boundary_id number) const
+Triangulation<dim, spacedim>::get_boundary (const types::manifold_id m_number) const
 {
-  Assert(number<numbers::internal_face_boundary_id, ExcIndexRange(number,0,numbers::internal_face_boundary_id));
+  types::boundary_id number = static_cast<types::boundary_id>(m_number);
+//  Assert(number<numbers::internal_face_boundary_id, ExcIndexRange(number,0,numbers::internal_face_boundary_id));
 
   //look, if there is a boundary stored at
   //boundary_id number.
-  typename std::map<types::boundary_id, SmartPointer<const Boundary<dim, spacedim>, Triangulation<dim, spacedim> > >::const_iterator it
-    = boundary.find(number);
+  typename std::map<types::manifold_id, SmartPointer<const Boundary<dim, spacedim>, Triangulation<dim, spacedim> > >::const_iterator it
+    = boundary.find(m_number);
 
   if (it != boundary.end())
     {
@@ -9593,6 +9681,34 @@ Triangulation<dim, spacedim>::get_boundary_indicators () const
     }
 }
 
+template <int dim, int spacedim>
+std::vector<types::manifold_id>
+Triangulation<dim, spacedim>::get_manifold_ids () const
+{
+//TODO: if manifold_id_t is a real integer type, this might become
+//a humongously large array...
+                                  // Just store all possible values
+  std::vector<bool> mi_exists(numbers::invalid_manifold_id+1, false);
+  active_cell_iterator cell=begin_active();
+  for (; cell!=end(); ++cell)
+    {
+      mi_exists[cell->manifold_id()] = true;
+      if(dim>1)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         mi_exists[cell->face(face)->manifold_id()]=true;
+    }
+  
+  const unsigned int n_mi=
+    std::count(mi_exists.begin(), mi_exists.end(), true)-1;
+
+  std::vector<types::manifold_id> manifold_indicators(n_mi);
+  unsigned int mi_counter=0;
+  for (unsigned int i=0; i<mi_exists.size()-1; ++i)
+    if (mi_exists[i]==true)
+      manifold_indicators[mi_counter++]=i;
+  
+  return manifold_indicators;
+}
 
 
 /*-----------------------------------------------------------------*/
@@ -9640,6 +9756,11 @@ copy_triangulation (const Triangulation<dim, spacedim> &old_tria)
       vertex_to_boundary_id_map_1d
         = (new std::map<unsigned int, types::boundary_id>
            (*old_tria.vertex_to_boundary_id_map_1d));
+
+      delete vertex_to_manifold_id_map_1d;
+      vertex_to_manifold_id_map_1d
+        = (new std::map<unsigned int, types::manifold_id>
+           (*old_tria.vertex_to_manifold_id_map_1d));
     }
 
   // inform those who are listening on old_tria of
index 16a72fdfbe3ad93960e6ff46a6c8b89c02756e22..1be216faebe69fbd23e48a7ead7b1cad973e4469 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2010 by the deal.II authors
+//    Copyright (C) 2010, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -16,5 +16,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
   {
 #if deal_II_dimension <= deal_II_space_dimension
     template class Triangulation<deal_II_dimension, deal_II_space_dimension>;
+#endif
+#if deal_II_dimension == deal_II_space_dimension
+    template class CellData<deal_II_dimension>;
 #endif
   }
index 0661e0369605779a4351f5700f60d04dd7bac827..f91265003266e540009f5d5b068df6a431074c80 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1096,6 +1096,29 @@ double TriaAccessor<3,3,3>::extent_in_direction(const unsigned int axis) const
 }
 
 
+// Recursively set manifold ids on hex iterators.
+template <>
+void
+TriaAccessor<3,3,3>::
+set_all_manifold_ids (const types::manifold_id manifold_ind) const
+{
+  set_manifold_id (manifold_ind);
+
+  if (this->has_children())
+    for (unsigned int c=0; c<this->n_children(); ++c)
+      this->child(c)->set_all_manifold_ids (manifold_ind);
+
+                                  // for hexes also set manifold_id
+                                  // of bounding quads and lines
+  
+                                  // Six bonding quads
+  for(unsigned int i=0; i<6; ++i)
+    this->quad(i)->set_manifold_id(manifold_ind);
+                                  // Twelve bounding lines
+  for (unsigned int i=0; i<12; ++i)
+    this->line(i)->set_manifold_id (manifold_ind);
+}
+
 
 /*------------------------ Functions: CellAccessor<1> -----------------------*/
 
index 27e09b129dc35fd11a9ec7cbbd985998fdb78233..b8e092d5c55692021adf4437573022519e7310e3 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2006, 2007, 2010, 2011, 2012 by the deal.II authors
+//    Copyright (C) 2006, 2007, 2010, 2011, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -195,6 +195,13 @@ namespace internal
                                           new_size-boundary_or_material_id.size(),
                                           BoundaryOrMaterialId());
 
+          manifold_id.reserve (new_size);
+          manifold_id.insert (manifold_id.end(),
+                             new_size-manifold_id.size(),
+                             numbers::flat_manifold_id);
+
+
+         
           user_data.reserve (new_size);
           user_data.insert (user_data.end(),
                             new_size-user_data.size(),
@@ -275,6 +282,12 @@ namespace internal
                                           new_size-boundary_or_material_id.size(),
                                           BoundaryOrMaterialId());
 
+          manifold_id.reserve (new_size);
+          manifold_id.insert (manifold_id.end(),
+                             new_size-manifold_id.size(),
+                             numbers::flat_manifold_id);
+         
+
           user_data.reserve (new_size);
           user_data.insert (user_data.end(),
                             new_size-user_data.size(),
@@ -389,6 +402,8 @@ namespace internal
               ExcMemoryInexact (cells.size(), children.size()));
       Assert (cells.size() == boundary_or_material_id.size(),
               ExcMemoryInexact (cells.size(), boundary_or_material_id.size()));
+      Assert (cells.size() == manifold_id.size(),
+              ExcMemoryInexact (cells.size(), manifold_id.size()));
       Assert (cells.size() == user_data.size(),
               ExcMemoryInexact (cells.size(), user_data.size()));
     }
@@ -408,6 +423,8 @@ namespace internal
               ExcMemoryInexact (cells.size(), refinement_cases.size()));
       Assert (cells.size() == boundary_or_material_id.size(),
               ExcMemoryInexact (cells.size(), boundary_or_material_id.size()));
+      Assert (cells.size() == manifold_id.size(),
+              ExcMemoryInexact (cells.size(), manifold_id.size()));
       Assert (cells.size() == user_data.size(),
               ExcMemoryInexact (cells.size(), user_data.size()));
     }
@@ -424,6 +441,8 @@ namespace internal
               ExcMemoryInexact (cells.size(), children.size()));
       Assert (cells.size() == boundary_or_material_id.size(),
               ExcMemoryInexact (cells.size(), boundary_or_material_id.size()));
+      Assert (cells.size() == manifold_id.size(),
+              ExcMemoryInexact (cells.size(), manifold_id.size()));
       Assert (cells.size() == user_data.size(),
               ExcMemoryInexact (cells.size(), user_data.size()));
       Assert (cells.size() * GeometryInfo<3>::faces_per_cell
@@ -466,6 +485,7 @@ namespace internal
       used.clear();
       user_flags.clear();
       boundary_or_material_id.clear();
+      manifold_id.clear();
       user_data.clear();
       user_data_type = data_unknown;
     }
@@ -498,6 +518,7 @@ namespace internal
               MemoryConsumption::memory_consumption (used) +
               MemoryConsumption::memory_consumption (user_flags) +
               MemoryConsumption::memory_consumption (boundary_or_material_id) +
+             MemoryConsumption::memory_consumption (manifold_id) +
               MemoryConsumption::memory_consumption (refinement_cases) +
               user_data.capacity() * sizeof(UserData) + sizeof(user_data));
     }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.