From: pauletti Date: Wed, 8 Dec 2010 20:30:57 +0000 (+0000) Subject: Implementation of direction_flags for dim=spacedim-1 meshes X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bd34964987d3196003baf437ab79f050269f3ea2;p=dealii-svn.git Implementation of direction_flags for dim=spacedim-1 meshes git-svn-id: https://svn.dealii.org/trunk@22944 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/grid/tria.h b/deal.II/include/deal.II/grid/tria.h index 0c634137c4..57011a9c89 100644 --- a/deal.II/include/deal.II/grid/tria.h +++ b/deal.II/include/deal.II/grid/tria.h @@ -1777,6 +1777,12 @@ class Triangulation : public Subscriptor const std::vector > &cells, const SubCellData &subcelldata); +/** + Revert or flip the direction_flags of a dim * 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 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 * anyway. */ void operator = (const CellAccessor &); + + template friend class Triangulation; + + friend class internal::Triangulation::Implementation; + }; diff --git a/deal.II/include/deal.II/grid/tria_levels.h b/deal.II/include/deal.II/grid/tria_levels.h index 6f2f0dbd06..0977a4274b 100644 --- a/deal.II/include/deal.II/grid/tria_levels.h +++ b/deal.II/include/deal.II/grid/tria_levels.h @@ -145,7 +145,16 @@ namespace internal * index their parent has. */ std::vector 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 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 > neighbors; std::vector parents; std::vector subdomain_ids; + std::vector 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; diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index 8cf4b63d0b..65025dee53 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -12,6 +12,8 @@ //--------------------------------------------------------------------------- #include +#include + #include #include #include @@ -23,10 +25,12 @@ #include #include #include +#include #include #include #include +#include #include #include @@ -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; cchild(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 > &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 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::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 +void +Triangulation:: +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()); } diff --git a/deal.II/source/grid/tria_accessor.cc b/deal.II/source/grid/tria_accessor.cc index 4fd0c88c7c..b5eacc9375 100644 --- a/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/source/grid/tria_accessor.cc @@ -1262,6 +1262,35 @@ CellAccessor::set_subdomain_id (const types::subdomain_id_t new_s } +template +bool CellAccessor::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 +void +CellAccessor::set_direction_flag (const bool new_direction_flag) const +{ + if (dimused(), TriaAccessorExceptions::ExcCellNotUsed()); + this->tria->levels[this->present_level]->direction_flags[this->present_index] + = new_direction_flag; + } +} + + + template TriaIterator > CellAccessor::parent () const diff --git a/deal.II/source/grid/tria_levels.cc b/deal.II/source/grid/tria_levels.cc index b5f4766605..7718c84f23 100644 --- a/deal.II/source/grid/tria_levels.cc +++ b/deal.II/source/grid/tria_levels.cc @@ -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(),