* exception of type ExcPureFunctionCalled.
*/
virtual
- void project_to_manifold (Point<spacedim> &candidate) const;
+ Point<spacedim> project_to_manifold (const Point<spacedim> candidate) const;
+
+ /**
+ * Given a vector of points, return the normals to the Manifold in
+ * those points. Notice that the Manifold may or may not be a
+ * codimension one manifold. If it is not, this function will throw
+ * an Error. Input arguments must be of the same size. The default
+ * implementation calls for each point the function
+ * normal_vector. Derived classes can overload that function.
+ */
+ virtual
+ void get_normals_at_points(std::vector<Point<spacedim> > &normals,
+ const std::vector<Point<spacedim> > &points) const;
+
+ /**
+ * Compute the normal to the manifold at the given point. The
+ * default implementation of this function trhows an
+ * Exception. Notice that this function only makes sense when the
+ * manifold is a codimension one manifold.
+ */
+ virtual
+ Point<spacedim> normal_vector(const Point<spacedim> point) const;
+
+
};
/**
* Let the new point be the average sum of surrounding vertices.
*
- * Refer to the general documentation of this class and the
- * documentation of the base class for more information.
+ * This particular implementation constructs the weighted average of
+ * the surrounding points, and then calss internally the function
+ * project_to_manifold. The reason why we do it this way, is to
+ * allow lazy programmers to implement only the project_to_manifold
+ * function for their own Manifold classes which are small (or
+ * trivial) perturbations of a flat manifold. For most simple
+ * geometries, it is possible to get reasonable results by deriving
+ * your own Manifold class from FlatManifold, and write a new
+ * interface only for the project_to_manifold function.
*/
virtual Point<spacedim>
get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const;
/**
- * Project to FlatManifold: do nothing.
+ * Project to FlatManifold: do nothing. Note however that this
+ * function can be overloaded by derived classes, which will then
+ * benefit from the logic behind the get_new_point class which are
+ * often very similar (if not identical) to the one implemented in
+ * this class.
*/
virtual
- void project_to_manifold (Point<spacedim> &candidate) const;
+ Point<spacedim> project_to_manifold (const Point<spacedim> candidate) const;
};
DEAL_II_NAMESPACE_OPEN
-template <int dim, int spacedim> class Boundary;
-template <int dim, int spacedim> class StraightBoundary;
template <int spacedim> class Manifold;
template <int spacedim> class FlatManifold;
typedef dealii::internal::Triangulation::Iterators<dim, spacedim> IteratorSelector;
public:
- /**
- * Default boundary object. This is used
- * for those boundaries for which no
- * boundary object has been explicitly
- * set using set_boundary().
- */
- static const StraightBoundary<dim,spacedim> straight_boundary;
-
/**
* Default manifold object. This is used
* for those objects for which no
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*/
void set_boundary (const types::manifold_id number,
- const Boundary<dim,spacedim> &boundary_object) DEAL_II_DEPRECATED;
+ const Manifold<spacedim> &boundary_object) DEAL_II_DEPRECATED;
/**
* Assign a Manifold object to a certain part of the
* triangulation. The Manifold object is used to find the location
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*/
- const Boundary<dim,spacedim> &get_boundary (const types::manifold_id number) const;
+ const Manifold<spacedim> &get_boundary (const types::manifold_id number) const;
/**
* Return a constant reference to a Manifold object used for this
* Collection of boundary objects. We store only objects, which are
* not of type StraightBoundary.
*/
- std::map<types::manifold_id, SmartPointer<const Boundary<dim, spacedim> , Triangulation<dim, spacedim> > > boundary;
+ std::map<types::manifold_id, SmartPointer<const Manifold<spacedim> , Triangulation<dim, spacedim> > > boundary;
/**
template <typename Accessor> class TriaIterator;
template <typename Accessor> class TriaActiveIterator;
-template <int dim, int spacedim> class Boundary;
template <int spacedim> class Manifold;
* Triangulation::get_boundary() function
* for the boundary object.
*/
- const Boundary<dim,spacedim> &get_boundary () const;
+ const Manifold<spacedim> &get_boundary () const;
/**
* Return a constant reference to a
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.templates.h>
-#include <deal.II/grid/manifold.h>
-#include <deal.II/grid/tria_boundary.h>
#include <deal.II/distributed/tria.h>
#include <cmath>
template <int structdim, int dim, int spacedim>
-const Boundary<dim,spacedim> &
+const Manifold<spacedim> &
TriaAccessor<structdim, dim, spacedim>::get_boundary () const
{
Assert (structdim<spacedim, ExcImpossibleInDim(dim));
* @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008
*/
template <int dim, int spacedim=dim>
- class Boundary : public Manifold<spacedim>
+ class Boundary : public FlatManifold<spacedim>
{
public:
virtual Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<spacedim>
+ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
/**
* Refer to the general documentation of
* this class and the documentation of the
Point<dim>
get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<dim>
+ get_new_point (const std::vector<Point<dim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
/**
* Refer to the general
* documentation of this class
Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<spacedim>
+ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
/**
* Refer to the general documentation of
* this class and the documentation of the
virtual Point<dim>
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<dim>
+ get_new_point (const std::vector<Point<dim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
/**
* Refer to the general
* documentation of this class
*/
virtual Point<dim>
get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<dim>
+ get_new_point (const std::vector<Point<dim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
/**
* Construct a new point on a quad.
virtual Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<spacedim>
+ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
/**
* Construct a new points on a line.
*/
// case in tests/deal.II/no_flux_11.
Point<dim> normal_vector
= (cell->face(face_no)->get_boundary()
- .normal_vector (cell->face(face_no),
- fe_values.quadrature_point(i)));
+ .normal_vector (fe_values.quadrature_point(i)));
if (normal_vector * fe_values.normal_vector(i) < 0)
normal_vector *= -1;
Assert (std::fabs(normal_vector.norm() - 1) < 1e-14,
{
// first get the normal vectors at the two vertices of this line
// from the boundary description
- const Boundary<dim> &boundary
- = line->get_triangulation().get_boundary(line->boundary_indicator());
-
- Boundary<dim>::FaceVertexNormals face_vertex_normals;
- boundary.get_normals_at_vertices (line, face_vertex_normals);
+ const Manifold<dim> &boundary
+ = line->get_boundary();
+
+ std::vector<Point<dim> > points(2);
+ std::vector<Point<dim> > face_vertex_normals(2);
+ points[0] = line->vertex(0);
+ points[1] = line->vertex(1);
+ boundary.get_normals_at_points (face_vertex_normals, points);
// then transform them into interpolation points for a cubic
// polynomial
fix_up_object (const Iterator &object,
const bool respect_manifold)
{
- const Boundary<Iterator::AccessorType::dimension,
- Iterator::AccessorType::space_dimension>
+ const Manifold<Iterator::AccessorType::space_dimension>
*manifold = (respect_manifold ?
&object->get_boundary() :
0);
/
eps);
else
- gradient[d]
- = ((objective_function (object,
- manifold->project_to_surface(object,
- object_mid_point + h))
- -
- objective_function (object,
- manifold->project_to_surface(object,
- object_mid_point - h)))
- /
- eps);
+ gradient[d]
+ = ((objective_function (object,
+ manifold->project_to_manifold(object_mid_point + h))
+ -
+ objective_function (object,
+ manifold->project_to_manifold(object_mid_point - h)))
+ /
+ eps);
}
// sometimes, the
gradient;
if (respect_manifold == true)
- object_mid_point = manifold->project_to_surface(object,
- object_mid_point);
+ object_mid_point = manifold->project_to_manifold(object_mid_point);
// compute current value of the
// objective function
+template <int spacedim>
+Point<spacedim>
+Manifold<spacedim>::project_to_manifold (const Point<spacedim> candidate) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return candidate;
+}
+
template <int spacedim>
void
-Manifold<spacedim>::project_to_manifold (Point<spacedim> &candidate) const
+Manifold<spacedim>::get_normals_at_points(std::vector<Point<spacedim> > &normals,
+ const std::vector<Point<spacedim> > &points) const
+{
+ AssertDimension(normals.size(), points.size());
+ for(unsigned int i=0; i<normals.size(); ++i)
+ normals[i] = normal_vector(points[i]);
+}
+
+template <int spacedim>
+Point<spacedim>
+Manifold<spacedim>::normal_vector(const Point<spacedim> point) const
{
Assert (false, ExcPureFunctionCalled());
+ return point;
}
template <int spacedim>
get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const
{
- Assert(surrounding_points.size() == weights.size(),
- ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+ AssertDimension(surrounding_points.size(), weights.size());
#ifdef DEBUG
double sum=0;
Point<spacedim> p;
for(unsigned int i=0; i<surrounding_points.size(); ++i)
p += surrounding_points[i]*weights[i];
- return p;
+ return project_to_manifold(p);
}
template <int spacedim>
-void
-FlatManifold<spacedim>::project_to_manifold (Point<spacedim> &candidate) const
-{}
+Point<spacedim>
+FlatManifold<spacedim>::project_to_manifold (const Point<spacedim> candidate) const
+{
+ return candidate;
+}
/* -------------------------- ManifoldChart --------------------- */
get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const
{
- Assert(surrounding_points.size() == weights.size(),
- ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+ AssertDimension(surrounding_points.size(), weights.size());
chart_points.resize(surrounding_points.size());
for(unsigned int i=0; i<surrounding_points.size(); ++i)
void
prevent_distorted_boundary_cells (Triangulation<dim,spacedim> &triangulation)
{
- // If the codimension is
- // one, we cannot perform
- // this check yet.
- if (spacedim>dim) return;
-
+ std::vector<double> wp;
+ std::vector<Point<spacedim> > sp;
for (typename Triangulation<dim,spacedim>::cell_iterator
cell=triangulation.begin(); cell!=triangulation.end(); ++cell)
if (cell->at_boundary() &&
// the new point on the
// boundary would be
// this one.
+ get_default_points_and_weights(sp,wp,face);
const Point<spacedim> new_bound
- = face->get_boundary().get_new_point_on_face(face);
+ = face->get_boundary().get_new_point(sp,wp);
// triangulation.get_boundary(face->boundary_indicator())
// .get_new_point_on_face (face);
// to check it,
}
-template <int dim, int spacedim>
-const StraightBoundary<dim,spacedim>
-Triangulation<dim, spacedim>::straight_boundary = StraightBoundary<dim,spacedim>();
-
-
template <int dim, int spacedim>
const FlatManifold<spacedim>
Triangulation<dim, spacedim>::flat_manifold = FlatManifold<spacedim>();
template <int dim, int spacedim>
void
Triangulation<dim, spacedim>::set_boundary (const types::manifold_id m_number,
- const Boundary<dim, spacedim> &boundary_object)
+ const Manifold<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));
+ Assert(number < numbers::internal_face_boundary_id,
+ ExcIndexRange(number,0,numbers::internal_face_boundary_id));
boundary[m_number] = &boundary_object;
}
template <int dim, int spacedim>
-const Boundary<dim, spacedim> &
+const Manifold<spacedim> &
Triangulation<dim, spacedim>::get_boundary (const types::manifold_id m_number) const
{
//look, if there is a boundary stored at
//boundary_id number.
- typename std::map<types::manifold_id, SmartPointer<const Boundary<dim, spacedim>, Triangulation<dim, spacedim> > >::const_iterator it
+ typename std::map<types::manifold_id,
+ SmartPointer<const Manifold<spacedim>, Triangulation<dim, spacedim> > >::const_iterator it
= boundary.find(m_number);
if (it != boundary.end())
//if we have not found an entry
//connected with number, we return
//straight_boundary
- return straight_boundary;
+ return flat_manifold;
}
}
faces = new internal::Triangulation::TriaFaces<dim>(*old_tria.faces);
- typename std::map<types::boundary_id, SmartPointer<const Boundary<dim, spacedim> , Triangulation<dim, spacedim> > >::const_iterator
+ typename std::map<types::boundary_id,
+ SmartPointer<const Manifold<spacedim> , Triangulation<dim, spacedim> > >::const_iterator
bdry_iterator = old_tria.boundary.begin();
for (; bdry_iterator != old_tria.boundary.end() ; bdry_iterator++)
boundary[bdry_iterator->first] = bdry_iterator->second;
compute_radius_automatically(false)
{}
+template <int dim, int spacedim>
+Point<spacedim>
+HyperBallBoundary<dim,spacedim>::get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const
+{
+ Assert(surrounding_points.size() == weights.size(),
+ ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+
+ double radius = 0;
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ radius += weights[i]*(surrounding_points[i]-center).norm();
+
+ Point<spacedim> p = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+ p = p-center;
+ p = p/p.norm()*radius+center;
+ return p;
+}
template <int dim, int spacedim>
projects / dealii-svn.git / commitdiff