]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implementation of direction_flags for dim=spacedim-1 meshes
authorpauletti <pauletti@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 Dec 2010 20:30:57 +0000 (20:30 +0000)
committerpauletti <pauletti@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 Dec 2010 20:30:57 +0000 (20:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@22944 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/grid/tria.h
deal.II/include/deal.II/grid/tria_accessor.h
deal.II/include/deal.II/grid/tria_levels.h
deal.II/source/grid/tria.cc
deal.II/source/grid/tria_accessor.cc
deal.II/source/grid/tria_levels.cc

index 0c634137c480c269e2b96f8792d1fb64ee3179f6..57011a9c89926f9cfcb608e05758dc3c7db69bd9 100644 (file)
@@ -1777,6 +1777,12 @@ class Triangulation : public Subscriptor
       const std::vector<CellData<dim> > &cells,
       const SubCellData                 &subcelldata);
 
+/**
+   Revert or flip the direction_flags of a dim<spacedim triangulation.
+   For other dim==spacedim throws an exception.
+ */
+    void flip_all_direction_flags();
+    
                                     /**
                                      * Distort the grid by randomly
                                      * moving around all the vertices
@@ -3170,6 +3176,11 @@ class Triangulation : public Subscriptor
                    int,
                    << "You tried to do something on level " << arg1
                    << ", but this level is empty.");
+                                    /**
+                                     * Exception
+                                     * @ingroup Exceptions
+                                     */
+    DeclException0 (ExcNonOrientableTriangulation);
                                     /*@}*/
   protected:
                                     /**
index b1cce904c960c75f13dd02fe718d7e9e86ed84b2..94dbceda5c9d2a2660d7361f253b877b59572973 100644 (file)
@@ -2610,6 +2610,25 @@ class CellAccessor :  public TriaAccessor<dim,dim,spacedim>
                                      * object.
                                      */
     void recursively_set_subdomain_id (const types::subdomain_id_t new_subdomain_id) const;
+                                     /**
+                                     * @}
+                                     */
+                                     
+                                     /**
+                                     *  @name Dealing with codim 1 cell orientation
+                                     */
+                                    /**
+                                     * @{
+                                     */
+
+                                    /**
+                                     * Return the orientation of
+                                     * this cell.
+                                     *
+                                     */
+    bool direction_flag () const;
+
+                                    
 
                                     /**
                                      *  Return an iterator to the
@@ -2790,7 +2809,12 @@ class CellAccessor :  public TriaAccessor<dim,dim,spacedim>
     unsigned int neighbor_of_neighbor_internal (const unsigned int neighbor) const;
 
   private:
-
+                                    /**
+                                     * Set the orientation of this
+                                     * cell.
+                                     *
+                                     */
+    void set_direction_flag (const bool new_direction_flag) const;
                                     /**
                                      *  Copy operator. This is
                                      *  normally used in a context
@@ -2812,6 +2836,11 @@ class CellAccessor :  public TriaAccessor<dim,dim,spacedim>
                                      *  anyway.
                                      */
     void operator = (const CellAccessor<dim, spacedim> &);
+
+    template <int, int> friend class Triangulation;
+
+    friend class internal::Triangulation::Implementation;
+    
 };
 
 
index 6f2f0dbd06db4ec60450680c38cc965b37ce2f57..0977a4274b893b523948d56938003c6bd2ca435d 100644 (file)
@@ -145,7 +145,16 @@ namespace internal
                                           * index their parent has.
                                           */
         std::vector<int> parents;
-
+       
+                                        /**
+                                         * One bool per cell to indicate the
+                                         * direction of the normal
+                                         * true:  use orientation from vertex
+                                         * false: revert the orientation.
+                                         * It is used for codim==1 meshes
+                                         */
+        std::vector<bool> direction_flags;
+       
                                          /**
                                           *  Reserve enough space to accomodate
                                           *  @p total_cells cells on this level.
@@ -161,6 +170,7 @@ namespace internal
                                           *  on the dimensions, you have to pass
                                           *  that additionally.
                                           */
+
         void reserve_space (const unsigned int total_cells,
                             const unsigned int dimension);
 
@@ -221,6 +231,7 @@ namespace internal
         std::vector<std::pair<int,int> > neighbors;
         std::vector<int> parents;
         std::vector<types::subdomain_id_t> subdomain_ids;
+        std::vector<bool> direction_flags; //not use; only needed to allow compilation
         void reserve_space (const unsigned int total_cells,
                             const unsigned int dimension);
         void monitor_memory (const unsigned int true_dimension) const;
index 8cf4b63d0be3d9be52dba5f13737c05942a4dbfd..65025dee53e27ab1c23a1a10dcd96635aa205b9e 100644 (file)
@@ -12,6 +12,8 @@
 //---------------------------------------------------------------------------
 
 #include <base/memory_consumption.h>
+#include <base/table.h>
+
 #include <grid/tria.h>
 #include <grid/tria_levels.h>
 #include <grid/tria_faces.h>
 #include <grid/magic_numbers.h>
 #include <fe/mapping_q1.h>
 #include <lac/vector.h>
+#include <lac/full_matrix.h>
 
 #include <algorithm>
 #include <numeric>
 #include <map>
+#include <list>
 #include <cmath>
 #include <functional>
 
@@ -4177,7 +4181,7 @@ namespace internal
                                                 // properties
                subcells[i]->set_material_id (cell->material_id());
                subcells[i]->set_subdomain_id (cell->subdomain_id());
-
+       
                if (i%2==0)
                  subcells[i]->set_parent (cell->index ());
              }
@@ -4196,6 +4200,11 @@ namespace internal
                                             // refinement flag was
                                             // already cleared at the
                                             // beginning of this function
+
+           if (dim < spacedim)
+             for (unsigned int c=0; c<n_children; ++c)
+               cell->child(c)->set_direction_flag (cell->direction_flag());  
+         
          }
 
 
@@ -4374,6 +4383,8 @@ namespace internal
                                                         next_unused_vertex));
                      first_child->set_material_id (cell->material_id());
                      first_child->set_subdomain_id (cell->subdomain_id());
+                     first_child->set_direction_flag (cell->direction_flag());
+
                      first_child->set_parent (cell->index ());
 
                                                       // reset neighborship info (refer
@@ -4433,6 +4444,8 @@ namespace internal
                      second_child->set_neighbor (0, first_child);
                      second_child->set_material_id (cell->material_id());
                      second_child->set_subdomain_id (cell->subdomain_id());
+                     second_child->set_direction_flag (cell->direction_flag());
+
                      if (cell->neighbor(1).state() != IteratorState::valid)
                        second_child->set_neighbor (1, cell->neighbor(1));
                      else
@@ -9527,6 +9540,142 @@ create_triangulation (const std::vector<Point<spacedim> >    &v,
       AssertThrow (distorted_cells.distorted_cells.size() == 0,
                   distorted_cells);
     }
+
+
+/*
+    When the triangulation is a manifold (dim < spacedim), the normal field
+    provided from the map class depends on the order of the vertices.
+    It may happen that this normal field is discontinous.
+    The following code takes care that this is not the case by setting the
+    cell direction flag on those cell that produce the wrong orientation.
+
+    To determine if 2 neighbours have the same or opposite orientation
+    we use a table of truth.
+    Its entries are indexes by the local indeces of the common face.
+    For example if two elements share a face, and this face is
+    face 0 for element 0 and face 1 for element 1, then
+    table(0,1) will tell wether the orientation are the same (true) or
+    oppositte (false).
+
+    Even though there may be a combinatorial/graph theory argument to get
+    this table in any dimension, I tested by hand all the different possible
+    cases in 1D and 2D to generate the table.
+
+    Assuming that a surface respects the standard orientation for 2d meshes,
+    the tables of truth are symmetric and their true values are the following
+    1D curves:  (0,1)
+    2D surface: (0,1),(0,2),(1,3),(2,3)
+
+    We store this data using an n_faces x n_faces full matrix, which is actually
+    much bigger than the minimal data required, but it makes the code more readable.
+    
+  */
+  if (dim < spacedim) {
+
+    Table<2,bool> correct(GeometryInfo< dim >::faces_per_cell,
+                         GeometryInfo< dim >::faces_per_cell);
+    switch (dim) 
+      {
+       case 1:
+       {
+         bool values [][2] = {{false,true},
+                              {true,false} };
+         for (unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
+           for (unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
+             correct(i,j) = ( values[i][j]);
+         break;
+       }
+       case 2:
+       {
+         bool values [][4]= {{false,true ,true , false},
+                             {true ,false,false, true },
+                             {true ,false,false, true },
+                             {false,true ,true , false} };
+         for (unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
+           for (unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
+             correct(i,j) = ( values[i][j]);
+         break;
+       }
+       default:
+             Assert (false, ExcNotImplemented());
+      }
+  
+    
+    std::list<active_cell_iterator> this_round, next_round;
+    active_cell_iterator neighbor;
+
+    this_round.push_back (begin_active());
+    begin_active()->set_direction_flag (true);
+    begin_active()->set_user_flag ();
+
+    int round = 0, ncell=0;
+    
+    while (this_round.size() > 0) {
+      
+      for ( typename std::list<active_cell_iterator>::iterator cell = this_round.begin();
+           cell != this_round.end(); ++cell) {
+       
+       for (unsigned int i = 0; i < GeometryInfo< dim >::faces_per_cell; ++i) 
+         {
+           if ( !((*cell)->face(i)->at_boundary()) )
+             {
+               neighbor = (*cell)->neighbor(i);
+                               
+               unsigned int cf = (*cell)->face_index(i);
+               unsigned int j = 0;
+               while(neighbor->face_index(j) != cf)
+                 {++j;}
+                
+               if ( (correct(i,j) && !(*cell)->direction_flag())
+                    ||
+                    (!correct(i,j) && (*cell)->direction_flag()) )
+                 {
+                   if (neighbor->user_flag_set() == false)
+                     {
+                       neighbor->set_direction_flag (false);
+                       neighbor->set_user_flag ();
+                       next_round.push_back (neighbor);
+                     }
+                   else 
+                     Assert (neighbor->direction_flag() == false,
+                             ExcNonOrientableTriangulation());
+                   
+                 }
+             }
+           
+         }
+      }
+                
+//Before we quit let's check that if the triangulation is disconnected
+      if (next_round.size() == 0) {
+       for (active_cell_iterator cell = begin_active();
+            cell != end(); ++cell)
+         if (cell->user_flag_set() == false) {
+           next_round.push_back (cell);
+           cell->set_direction_flag (true);
+           cell->set_user_flag ();
+           break;
+         }
+      }
+      
+      this_round = next_round;
+      next_round.clear();
+    }
+  }
+}
+
+
+
+
+template <int dim, int spacedim>
+void
+Triangulation<dim,spacedim>::
+flip_all_direction_flags()
+{
+  AssertThrow (dim+1 == spacedim, ExcMessage ("Only works for dim == spacedim-1"));
+  for (active_cell_iterator cell = begin_active();
+       cell != end(); ++cell)
+    cell->set_direction_flag (!cell->direction_flag());
 }
 
 
index 4fd0c88c7ca757f7976f189fc49abea3d634147f..b5eacc9375f8ada7a560810f35359c767a76e141 100644 (file)
@@ -1262,6 +1262,35 @@ CellAccessor<dim, spacedim>::set_subdomain_id (const types::subdomain_id_t new_s
 }
 
 
+template <int dim, int spacedim>
+bool CellAccessor<dim, spacedim>::direction_flag () const
+{
+  if (dim==spacedim)
+    return true;
+  else
+    {
+      Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
+      return this->tria->levels[this->present_level]->direction_flags[this->present_index];
+    }
+  
+}
+
+
+
+template <int dim, int spacedim>
+void
+CellAccessor<dim, spacedim>::set_direction_flag (const bool new_direction_flag) const
+{
+  if (dim<spacedim)
+    {
+      Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
+      this->tria->levels[this->present_level]->direction_flags[this->present_index]
+       = new_direction_flag;
+    }
+}
+
+
+
 template <int dim, int spacedim>
 TriaIterator<CellAccessor<dim,spacedim> >
 CellAccessor<dim, spacedim>::parent () const
index b5f4766605171855b661fa56f3a4c9058f3d2658..7718c84f23b997775513d80cc5b71196931befef 100644 (file)
@@ -51,6 +51,11 @@ namespace internal
                                 total_cells - subdomain_ids.size(),
                                 0);
 
+         direction_flags.reserve (total_cells);
+          direction_flags.insert (direction_flags.end(),
+                                 total_cells - direction_flags.size(),
+                                 true);
+
           parents.reserve ((int) (total_cells + 1) / 2);
           parents.insert (parents.end (),
                           (total_cells + 1) / 2 - parents.size (),
@@ -91,6 +96,11 @@ namespace internal
              subdomain_ids.capacity() + DEAL_II_MIN_VECTOR_CAPACITY,
               ExcMemoryWasted ("subdomain_ids",
                                subdomain_ids.size(), subdomain_ids.capacity()));
+      Assert (direction_flags.size() <=
+             direction_flags.capacity() + DEAL_II_MIN_VECTOR_CAPACITY,
+              ExcMemoryWasted ("direction_flags",
+                               direction_flags.size(), direction_flags.capacity()));
+
       Assert (2*true_dimension*refine_flags.size() == neighbors.size(),
               ExcMemoryInexact (refine_flags.size(), neighbors.size()));
       Assert (2*true_dimension*coarsen_flags.size() == neighbors.size(),

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.