From: bangerth Date: Wed, 31 Aug 2011 01:19:35 +0000 (+0000) Subject: Implement getting and setting boundary indicators of triangulations in 1d as well. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7b82cf41586aed524a1fdc98224c769c4886a2c0;p=dealii-svn.git Implement getting and setting boundary indicators of triangulations in 1d as well. git-svn-id: https://svn.dealii.org/trunk@24224 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 5ca9275299..7e18a31931 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -283,6 +283,11 @@ and DoF handlers embedded in higher dimensional space have been added.

Specific improvements

    +
  1. New: Code like cell-@>face(1)-@>set_boundary_indicator(42); +now also works in 1d. +
    +(Wolfgang Bangerth, 2011/08/30) +
  2. Fixed: The TimerOutput::print_summary() function changed the precision of output on the stream it prints to, but didn't restore the previous value. This is now fixed. diff --git a/deal.II/include/deal.II/grid/tria.h b/deal.II/include/deal.II/grid/tria.h index ff8aa0c49a..6912b7d3f1 100644 --- a/deal.II/include/deal.II/grid/tria.h +++ b/deal.II/include/deal.II/grid/tria.h @@ -38,7 +38,7 @@ template class Boundary; template class StraightBoundary; template class TriaAccessor; -template class TriaAccessor; +template class TriaAccessor<0,1,spacedim>; namespace internal { @@ -667,7 +667,7 @@ namespace internal * refinement smoothes the mesh a bit. * * The function will make sure that vertices on restricted faces (hanging - * nodes) will result in the correct place, i.e. in the middle of the two + * nodes) will end up in the correct place, i.e. in the middle of the two * other vertices of the mother line, and the analogue in higher space * dimensions (vertices on the boundary are not corrected, so don't distort * boundary vertices in more than two space dimension, i.e. in dimensions @@ -945,7 +945,7 @@ namespace internal * *

    Material and boundary information

    * - * Each line, quad, etc stores one byte of information denoting the + * Each cell, face or edge stores information denoting the * material or the part of the boundary that an object * belongs to. The material of a cell may be used * during matrix generation in order to implement different @@ -986,16 +986,13 @@ namespace internal * are separated by a vertex (in 2D) or a line (in 3D) such that each boundary * line or quad has a unique boundary indicator. * - * Since in one dimension, no substructures of lower dimension exist to - * cells (of course apart from vertices, but these are handled - * in another way than the structures and substructures with dimension one - * or greater), there is no way to denote boundary indicators to boundary - * vertices (the endpoints). This is not a big thing, however, since you - * will normally not want to do a loop over all vertices, but rather work - * on the cells, which in turn have a possibility to find out whether they - * are at one of the endpoints. Only the handling of boundary values gets - * a bit more complicated, but this seems to be the price to be paid for - * the different handling of vertices from lines and quads. + * By default (unless otherwise specified during creation of a + * triangulation), all parts of the boundary have boundary indicator + * zero. As a historical wart, this isn't true for 1d meshes, however: For + * these, leftmost vertices have boundary indicator zero while rightmost + * vertices have boundary indicator one. In either case, the boundary + * indicator of a face can be changed using a call of the kind + * cell-@>face(1)-@>set_boundary_indicator(42);. * * * @@ -3744,6 +3741,31 @@ class Triangulation : public Subscriptor */ internal::Triangulation::NumberCache number_cache; + /** + * A map that relates the number of a + * boundary vertex to the boundary + * indicator. This field is only used in + * 1d. We have this field because we + * store boundary 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_boundary_id_map_1d; + /** * A map that correlates each refinement listener that has been added * through the outdated RefinementListener interface via @@ -3763,6 +3785,8 @@ class Triangulation : public Subscriptor // friends template friend class TriaAccessorBase; template friend class TriaAccessor; + friend class TriaAccessor<0, 1, spacedim>; + friend class CellAccessor; friend struct internal::TriaAccessor::Implementation; @@ -3874,6 +3898,9 @@ Triangulation::save (Archive & ar, ar & number_cache; ar & check_for_distorted_cells; + + if (dim == 1) + ar & vertex_to_boundary_id_map_1d; } @@ -3909,6 +3936,9 @@ Triangulation::load (Archive & ar, "same setting with regard to reporting distorted " "cell as the one previously stored.")); + if (dim == 1) + ar & vertex_to_boundary_id_map_1d; + // trigger the create signal to indicate // that new content has been imported into // the triangulation diff --git a/deal.II/include/deal.II/grid/tria_accessor.h b/deal.II/include/deal.II/grid/tria_accessor.h index d6564f7134..e5ba1977bc 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.h +++ b/deal.II/include/deal.II/grid/tria_accessor.h @@ -2020,17 +2020,40 @@ class TriaAccessor<0, 1, spacedim> int isotropic_child_index (const unsigned int i); /** - * @brief Do nothing and throw and error + * Set the boundary indicator. + * The same applies as for the + * boundary_indicator() + * function. + * + * @warning You should never set the + * boundary indicator of an interior face + * (a face not at the boundary of the + * domain), or set set the boundary + * indicator of an exterior face to 255 + * (this value is reserved for another + * purpose). Algorithms may not work or + * produce very confusing results if + * boundary cells have a boundary + * indicator of 255 or if interior cells + * have boundary indicators other than + * 255. Unfortunately, the current object + * has no means of finding out whether it + * really is at the boundary of the + * domain and so cannot determine whether + * the value you are trying to set makes + * sense under the current circumstances. */ - static void set_boundary_indicator (const unsigned char); /** - * @brief Do nothing and throw and error + * Since this object only represents a + * single vertex, call + * set_boundary_indicator with the same + * argument. */ - static - void set_all_boundary_indicators (const unsigned char); + void + set_all_boundary_indicators (const unsigned char); /** * 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 7ac0529a88..eb0879bbb8 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.templates.h +++ b/deal.II/include/deal.II/grid/tria_accessor.templates.h @@ -2209,9 +2209,15 @@ TriaAccessor<0, 1, spacedim>::boundary_indicator () const switch (vertex_kind) { case left_vertex: - return 0; case right_vertex: - return 1; + { + Assert (tria->vertex_to_boundary_id_map_1d->find (this->vertex_index()) + != tria->vertex_to_boundary_id_map_1d->end(), + ExcInternalError()); + + return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()]; + } + default: return 255; } @@ -2334,22 +2340,22 @@ int TriaAccessor<0, 1, spacedim>::isotropic_child_index (const unsigned int) template inline void -TriaAccessor<0, 1, spacedim>::set_boundary_indicator (const unsigned char) +TriaAccessor<0, 1, spacedim>::set_boundary_indicator (const unsigned char b) { - Assert(false, - ExcMessage ("The boundary indicator of vertices is determined by their " - "location within a triangulation and can not be set explicitly")); + Assert (tria->vertex_to_boundary_id_map_1d->find (this->vertex_index()) + != tria->vertex_to_boundary_id_map_1d->end(), + ExcInternalError()); + + (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b; } template inline -void TriaAccessor<0, 1, spacedim>::set_all_boundary_indicators (const unsigned char) +void TriaAccessor<0, 1, spacedim>::set_all_boundary_indicators (const unsigned char b) { - Assert(false, - ExcMessage ("The boundary indicator of vertices is determined by their " - "location within a triangulation and can not be set explicitly")); + set_boundary_indicator (b); } diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index 15431b627e..705431cefe 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -1570,6 +1570,19 @@ namespace internal lines_at_vertex[line->vertex_index(vertex)][0]); line->set_neighbor (vertex, neighbor); } + + // finally set the + // vertex_to_boundary_id_map_1d + // map + triangulation.vertex_to_boundary_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; } @@ -8944,7 +8957,7 @@ namespace internal void prevent_distorted_boundary_cells (const Triangulation<1,spacedim> &); - + template static void @@ -8954,7 +8967,7 @@ namespace internal // one, we cannot perform // this check yet. if(spacedim>dim) return; - + for (typename Triangulation::cell_iterator cell=triangulation.begin(); cell!=triangulation.end(); ++cell) if (cell->at_boundary() && @@ -9378,13 +9391,18 @@ Triangulation (const MeshSmoothing smooth_grid, smooth_grid(smooth_grid), faces(NULL), anisotropic_refinement(false), - check_for_distorted_cells(check_for_distorted_cells) + check_for_distorted_cells(check_for_distorted_cells), + vertex_to_boundary_id_map_1d (0) { // set default boundary for all // possible components for (unsigned int i=0;i<255;++i) boundary[i] = &straight_boundary; - + + if (dim == 1) + vertex_to_boundary_id_map_1d + = new std::map(); + // connect the any_change signal to the other signals signals.create.connect (signals.any_change); signals.post_refinement.connect (signals.any_change); @@ -9400,7 +9418,8 @@ Triangulation (const Triangulation & other) // is an error! : Subscriptor(), - check_for_distorted_cells(other.check_for_distorted_cells) + check_for_distorted_cells(other.check_for_distorted_cells), + vertex_to_boundary_id_map_1d (0) { Assert (false, ExcInternalError()); } @@ -9414,6 +9433,14 @@ Triangulation::~Triangulation () delete levels[i]; levels.clear (); delete faces; + + // the vertex_to_boundary_id_map_1d field + // should be unused except in 1d + Assert ((dim == 1) + || + (vertex_to_boundary_id_map_1d == 0), + ExcInternalError()); + delete vertex_to_boundary_id_map_1d; } @@ -9472,24 +9499,40 @@ template std::vector Triangulation::get_boundary_indicators () const { + // in 1d, we store a map of all used + // boundary indicators. use it for our + // purposes + if (dim == 1) + { + std::vector boundary_indicators; + for (std::map::const_iterator + p = vertex_to_boundary_id_map_1d->begin(); + p != vertex_to_boundary_id_map_1d->end(); + ++p) + boundary_indicators.push_back (p->second); + + return boundary_indicators; + } + else + { + std::vector bi_exists(255, false); + active_cell_iterator cell=begin_active(); + for (; cell!=end(); ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->at_boundary(face)) + bi_exists[cell->face(face)->boundary_indicator()]=true; - std::vector bi_exists(255, false); - active_cell_iterator cell=begin_active(); - for (; cell!=end(); ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->at_boundary(face)) - bi_exists[cell->face(face)->boundary_indicator()]=true; - - const unsigned int n_bi= - std::count(bi_exists.begin(), bi_exists.end(), true); + const unsigned int n_bi= + std::count(bi_exists.begin(), bi_exists.end(), true); - std::vector boundary_indicators(n_bi); - unsigned int bi_counter=0; - for (unsigned int i=0; i boundary_indicators(n_bi); + unsigned int bi_counter=0; + for (unsigned int i=0; i &old_tria) number_cache = old_tria.number_cache; + if (dim == 1) + { + delete vertex_to_boundary_id_map_1d; + vertex_to_boundary_id_map_1d + = (new std::map + (*old_tria.vertex_to_boundary_id_map_1d)); + } + // inform RefinementListeners of old_tria of // the copy operation old_tria.signals.copy (*this); @@ -9594,7 +9645,7 @@ create_triangulation (const std::vector > &v, } compute_number_cache (*this, levels.size(), number_cache); - + // now verify that there are indeed // no distorted cells. as per the // documentation of this class, we @@ -12437,7 +12488,7 @@ Triangulation::execute_coarsening_and_refinement () // Inform RefinementListeners // about beginning of refinement. signals.pre_refinement(); - + execute_coarsening(); const DistortedCellList @@ -14507,8 +14558,8 @@ template void Triangulation::add_refinement_listener (RefinementListener &listener) const { - // in this compatibility mode with the old-style refinement listeners, an - // external class presents itself as one that may or may not have + // in this compatibility mode with the old-style refinement listeners, an + // external class presents itself as one that may or may not have // overloaded all of the functions that the RefinementListener // class has. consequently, we need to connect each of its functions // to the relevant signals. for those functions that haven't been @@ -14516,25 +14567,25 @@ Triangulation::add_refinement_listener (RefinementListener &liste // to the function in the RefinementListener base class which simply // does nothing std::vector connections; - - connections.push_back + + connections.push_back (signals.create.connect (std_cxx1x::bind (&RefinementListener::create_notification, std_cxx1x::ref(listener), std_cxx1x::cref(*this)))); - connections.push_back + connections.push_back (signals.copy.connect (std_cxx1x::bind (&RefinementListener::copy_notification, std_cxx1x::ref(listener), std_cxx1x::cref(*this), std_cxx1x::_1))); - connections.push_back + connections.push_back (signals.pre_refinement.connect (std_cxx1x::bind (&RefinementListener::pre_refinement_notification, std_cxx1x::ref(listener), std_cxx1x::cref(*this)))); - connections.push_back + connections.push_back (signals.post_refinement.connect (std_cxx1x::bind (&RefinementListener::post_refinement_notification, std_cxx1x::ref(listener), std_cxx1x::cref(*this)))); - + // now push the set of connections into the map refinement_listener_map.insert (std::make_pair(&listener, connections)); } @@ -14548,14 +14599,14 @@ Triangulation::remove_refinement_listener (RefinementListener &li Assert (refinement_listener_map.find (&listener) != refinement_listener_map.end(), ExcMessage("You try to remove a refinement listener that does " "not appear to have been added previously.")); - + // get the element of the map, and terminate these // connections. then erase the element from the list std::vector connections = refinement_listener_map.find(&listener)->second; for (unsigned int i=0; iface() in 1d as well. here: +// setting boundary indicators + + +#include "../tests.h" +#include +#include +#include +#include + +#include + + +template +void test () +{ + Triangulation<1,spacedim> tria; + GridGenerator::hyper_cube (tria); + + tria.begin_active()->face(0)->set_boundary_indicator(2); + tria.begin_active()->face(1)->set_boundary_indicator(4); + + std::vector + boundary_ids = tria.get_boundary_indicators (); + + for (unsigned int i=0; i (); + test<2> (); + + return 0; +} diff --git a/tests/deal.II/vertex_as_face_10/cmp/generic b/tests/deal.II/vertex_as_face_10/cmp/generic new file mode 100644 index 0000000000..a9927d4493 --- /dev/null +++ b/tests/deal.II/vertex_as_face_10/cmp/generic @@ -0,0 +1,5 @@ + +DEAL::2 +DEAL::4 +DEAL::2 +DEAL::4 diff --git a/tests/deal.II/vertex_as_face_11.cc b/tests/deal.II/vertex_as_face_11.cc new file mode 100644 index 0000000000..534b552b88 --- /dev/null +++ b/tests/deal.II/vertex_as_face_11.cc @@ -0,0 +1,52 @@ +//---------------------------- vertex_as_face_11.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2010, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- vertex_as_face_11.cc --------------------------- + +// verify that we can do things like cell->face() in 1d as well. here: +// getting boundary indicators + + +#include "../tests.h" +#include +#include +#include +#include + +#include + + +template +void test () +{ + Triangulation<1,spacedim> tria; + GridGenerator::hyper_cube (tria); + + tria.begin_active()->face(0)->set_boundary_indicator(2); + tria.begin_active()->face(1)->set_boundary_indicator(4); + + deallog << (int)tria.begin_active()->face(0)->boundary_indicator() << std::endl; + deallog << (int)tria.begin_active()->face(1)->boundary_indicator() << std::endl; +} + + + +int main () +{ + std::ofstream logfile("vertex_as_face_11/output"); + deallog.attach(logfile); + deallog.depth_console(0); + + test<1> (); + test<2> (); + + return 0; +} diff --git a/tests/deal.II/vertex_as_face_11/cmp/generic b/tests/deal.II/vertex_as_face_11/cmp/generic new file mode 100644 index 0000000000..a9927d4493 --- /dev/null +++ b/tests/deal.II/vertex_as_face_11/cmp/generic @@ -0,0 +1,5 @@ + +DEAL::2 +DEAL::4 +DEAL::2 +DEAL::4