From af86461c441f25f8a29300318d7449c76ef99cb4 Mon Sep 17 00:00:00 2001 From: heltai Date: Tue, 16 Apr 2013 18:18:41 +0000 Subject: [PATCH] Backward compatible manifold_id branch 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 | 32 +- deal.II/include/deal.II/grid/tria.h | 71 ++++- deal.II/include/deal.II/grid/tria_accessor.h | 99 +++++++ .../deal.II/grid/tria_accessor.templates.h | 104 ++++++- deal.II/include/deal.II/grid/tria_objects.h | 11 +- deal.II/source/grid/tria.cc | 275 +++++++++++++----- deal.II/source/grid/tria.inst.in | 5 +- deal.II/source/grid/tria_accessor.cc | 25 +- deal.II/source/grid/tria_objects.cc | 23 +- 9 files changed, 555 insertions(+), 90 deletions(-) diff --git a/deal.II/include/deal.II/base/types.h b/deal.II/include/deal.II/base/types.h index 2a37e62c29..ec54568f5e 100644 --- a/deal.II/include/deal.II/base/types.h +++ b/deal.II/include/deal.II/base/types.h @@ -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(-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(-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(-1); + /** * A special id for an invalid * subdomain id. This value may not diff --git a/deal.II/include/deal.II/grid/tria.h b/deal.II/include/deal.II/grid/tria.h index 1962bd3a04..ce1fa0c4a8 100644 --- a/deal.II/include/deal.II/grid/tria.h +++ b/deal.II/include/deal.II/grid/tria.h @@ -92,6 +92,11 @@ template class MGDoFHandler; template 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 &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 &get_boundary (const types::boundary_id number) const; + const Boundary &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 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 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 , Triangulation > > boundary; + std::map , Triangulation > > boundary; /** * Flag indicating whether @@ -3604,6 +3637,34 @@ private: * were not a pointer. */ std::map *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 *vertex_to_manifold_id_map_1d; + /** * A map that correlates each refinement listener that has been added diff --git a/deal.II/include/deal.II/grid/tria_accessor.h b/deal.II/include/deal.II/grid/tria_accessor.h index f4e2167f3f..d8bb9d2833 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.h +++ b/deal.II/include/deal.II/grid/tria_accessor.h @@ -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 + * manifold_id() 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 diff --git a/deal.II/include/deal.II/grid/tria_accessor.templates.h b/deal.II/include/deal.II/grid/tria_accessor.templates.h index 5379e3a5af..bb1b20ec80 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.templates.h +++ b/deal.II/include/deal.II/grid/tria_accessor.templates.h @@ -1918,11 +1918,79 @@ template const Boundary & TriaAccessor::get_boundary () const { - Assert (structdimused(), 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 + (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 +types::manifold_id +TriaAccessor::manifold_id () const +{ + Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + + return this->objects().manifold_id[this->present_index]; +} + + + +template +void +TriaAccessor:: +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 +void +TriaAccessor:: +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; cn_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 +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 inline @@ -2420,6 +2500,15 @@ TriaAccessor<0, 1, spacedim>::set_boundary_indicator (const types::boundary_id b } +template +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 inline @@ -2430,6 +2519,15 @@ void TriaAccessor<0, 1, spacedim>::set_all_boundary_indicators (const types::bou +template +inline +void TriaAccessor<0, 1, spacedim>::set_all_manifold_ids (const types::manifold_id b) +{ + set_manifold_id (b); +} + + + template inline bool TriaAccessor<0, 1, spacedim>::used () const diff --git a/deal.II/include/deal.II/grid/tria_objects.h b/deal.II/include/deal.II/grid/tria_objects.h index c90e306a62..6c58fa78f6 100644 --- a/deal.II/include/deal.II/grid/tria_objects.h +++ b/deal.II/include/deal.II/grid/tria_objects.h @@ -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 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 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; } diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index 4350049bcc..02d7e5d2bb 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -39,6 +39,16 @@ DEAL_II_NAMESPACE_OPEN +template +CellData::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::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) for (unsigned int f=0; f::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::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::max(), - ExcIndexRange(cell->material_id(),0,std::numeric_limits::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(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(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 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(); + { + vertex_to_boundary_id_map_1d + = new std::map(); + vertex_to_manifold_id_map_1d + = new std::map(); + } // connect the any_change signal to the other signals signals.create.connect (signals.any_change); @@ -9456,7 +9533,8 @@ Triangulation (const Triangulation &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::~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::set_mesh_smoothing(const MeshSmoothing mesh_smooth template void -Triangulation::set_boundary (const types::boundary_id number, +Triangulation::set_boundary (const types::manifold_id m_number, const Boundary &boundary_object) { + types::boundary_id number = static_cast(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 void -Triangulation::set_boundary (const types::boundary_id number) +Triangulation::set_boundary (const types::manifold_id m_number) { + types::boundary_id number = static_cast(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 const Boundary & -Triangulation::get_boundary (const types::boundary_id number) const +Triangulation::get_boundary (const types::manifold_id m_number) const { - Assert(number(m_number); +// Assert(number, Triangulation > >::const_iterator it - = boundary.find(number); + typename std::map, Triangulation > >::const_iterator it + = boundary.find(m_number); if (it != boundary.end()) { @@ -9593,6 +9681,34 @@ Triangulation::get_boundary_indicators () const } } +template +std::vector +Triangulation::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 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::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 manifold_indicators(n_mi); + unsigned int mi_counter=0; + for (unsigned int i=0; i &old_tria) vertex_to_boundary_id_map_1d = (new std::map (*old_tria.vertex_to_boundary_id_map_1d)); + + delete vertex_to_manifold_id_map_1d; + vertex_to_manifold_id_map_1d + = (new std::map + (*old_tria.vertex_to_manifold_id_map_1d)); } // inform those who are listening on old_tria of diff --git a/deal.II/source/grid/tria.inst.in b/deal.II/source/grid/tria.inst.in index 16a72fdfbe..1be216faeb 100644 --- a/deal.II/source/grid/tria.inst.in +++ b/deal.II/source/grid/tria.inst.in @@ -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; +#endif +#if deal_II_dimension == deal_II_space_dimension + template class CellData; #endif } diff --git a/deal.II/source/grid/tria_accessor.cc b/deal.II/source/grid/tria_accessor.cc index 0661e03696..f912650032 100644 --- a/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/source/grid/tria_accessor.cc @@ -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; cn_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> -----------------------*/ diff --git a/deal.II/source/grid/tria_objects.cc b/deal.II/source/grid/tria_objects.cc index 27e09b129d..b8e092d5c5 100644 --- a/deal.II/source/grid/tria_objects.cc +++ b/deal.II/source/grid/tria_objects.cc @@ -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)); } -- 2.39.5