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
{
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);
// 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;
+ }
}
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);
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
ExcInteriorLineCantBeBoundary());
line->set_boundary_indicator (boundary_line->boundary_id);
+ line->set_manifold_id (boundary_line->manifold_id);
}
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);
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);
}
/////////////////////////////////////////
// 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
"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);
}
"if they carry the same boundary indicator."));
quad->set_boundary_indicator (boundary_quad->boundary_id);
+ quad->set_manifold_id (boundary_quad->manifold_id);
}
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();
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();
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
// 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);
}
}
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)
// 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)
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
::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>
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());
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
// 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
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
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
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
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
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
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();
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();
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();
// 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
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
// 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
// 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
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
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
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
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
// 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)
// 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
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);
:
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());
}
(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;
}
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())
{
}
}
+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;
+}
/*-----------------------------------------------------------------*/
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