]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Default manifold is set to boundary_id when that one is set and manifold_id is invali...
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 22 Aug 2013 03:05:52 +0000 (03:05 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 22 Aug 2013 03:05:52 +0000 (03:05 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30403 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-1/step-1.cc
deal.II/include/deal.II/grid/manifold_lib.h
deal.II/include/deal.II/grid/tria_accessor.templates.h
deal.II/source/grid/manifold_lib.cc
deal.II/source/grid/manifold_lib.inst.in

index 42def746f606e087cff426dd6cba3135fe984bd0..724dd45e95f0afd702e6ba03da5007adc0b09780 100644 (file)
@@ -105,16 +105,7 @@ void second_grid ()
 
    Triangulation<2>::active_cell_iterator
       cell = triangulation.begin_active(),
-      endc = triangulation.end();
-   for(; cell != endc; ++cell)
-     cell->set_all_manifold_ids(0);
-   
-       
-     // for(unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
-     //   if(cell->face(f)->at_boundary())
-     //         cell->face(f)->set_manifold_id(0);
-   
-       
+      endc = triangulation.end();       
   
   // By default, the triangulation assumes that all boundaries are straight
   // and given by the cells of the coarse grid (which we just created). It
index 97d203e8fcfd6e03da95e8d50101bec974fcb9fa..e426dfa67f4d9570b5cfe4e7751ed3f2f9625962 100644 (file)
@@ -25,190 +25,143 @@ DEAL_II_NAMESPACE_OPEN
 
 
 /**
- * Manifold description for a polar space coordinate system. Given any
- * two points in space, tranform them in polar coordinate centered at
- * center, compute intermediate points in polar coordinate systems,
- * and returns its space representation.
+ * Manifold description for a polar space coordinate system. Given a
+ * set of points in space, this class computes their weighted average
+ * and the weighted average of their distance from the center. The
+ * average of the points is then shifted in the radial direction, to
+ * match the weighted average of the distances.
+ *
+ * You can use this Manifold object to describe any sphere, circle,
+ * hypersphere or hyperdisc in two or three dimensions, both as a
+ * co-dimension one manifold descriptor or as co-dimension zero
+ * manifold descriptor.
+ *
+ * @ingroup manifold
+ *
+ * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007, Luca Heltai,
+ * 2013
  */
 template <int spacedim>
 class PolarManifold : public FlatManifold<spacedim> 
 {
   public:
+                                    /**
+                                     * The Constructor takes the
+                                     * center of the polar
+                                     * coordinates system.
+                                     */
     PolarManifold(const Point<spacedim> center = Point<spacedim>());
 
+                                    /**
+                                     * Given a set of points in
+                                     * space, this function computes
+                                     * their weighted average and the
+                                     * weighted average of their
+                                     * distance from the center. The
+                                     * average of the points is then
+                                     * shifted in the radial
+                                     * direction, to match the
+                                     * weighted average of the
+                                     * distances.
+                                     */
     virtual Point<spacedim>
     get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
                   const std::vector<double> &weights) const;
 
     
   private:
-    Point<spacedim> center;
+                                    /**
+                                     * The center of the polar
+                                     * coordinate system.
+                                     */
+    const Point<spacedim> center;
+};
+
+/**
+ * Manifold description for a cylindrical space coordinate system. In
+ * three dimensions, points are projected on a cylindrical tube along
+ * the <tt>x-</tt>, <tt>y-</tt> or <tt>z</tt>-axis (when using the
+ * first constructor of this class), or an arbitrarily oriented
+ * cylinder described by the direction of its axis and a point located
+ * on the axis. Similar to PolarManifold, new points are projected by
+ * taking the weighetd averages of the old points and adjusting the
+ * radius from the axis.
+ *
+ * This class was developed to be used in conjunction with the @p
+ * cylinder function of GridGenerator. Its use is discussed in detail
+ * in the results section of step-49.
+ *
+ * @ingroup manifold
+ *
+ * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007, Luca Heltai,
+ * 2013
+ */
+
+template <int spacedim>
+class CylinderManifold : public FlatManifold<spacedim> 
+{
+  public:
+                                    /**
+                                     * Constructor. Using default
+                                     * values for the constructor
+                                     * arguments yields a cylindrical
+                                     * manifold tube along the x-axis
+                                     * (<tt>axis=0</tt>). Choose
+                                     * <tt>axis=1</tt> or
+                                     * <tt>axis=2</tt> for a
+                                     * cylindrical manifold along the
+                                     * y- or z-axis, respectively.
+                                     */
+    CylinderManifold(const int axis=0);
+    
+                                    /**
+                                     * Constructor. If constructed
+                                     * with this constructor, the
+                                     * boundary described is a
+                                     * cylinder with an axis that
+                                     * points in direction #direction
+                                     * and goes through the given
+                                     * #point_on_axis. The direction
+                                     * may be arbitrarily scaled, and
+                                     * the given point may be any
+                                     * point on the axis.
+                                     */
+    CylinderManifold(const Point<spacedim> &direction,
+                    const Point<spacedim> &point_on_axis);
+    
+                                    /**
+                                     * 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;
+
+    
+  private:
+                                    /**
+                                     * Direction of the cylindrical
+                                     * coordinate system.
+                                     */
+    const Point<spacedim> direction;
+    
+                                    /**
+                                     *  One point on the axis of the
+                                     *  cylindrical coordinate systems
+                                     */
+    const Point<spacedim> point_on_axis;
+
+                                    /**
+                                     * Returns the axis vector
+                                     * associated with the given axis
+                                     * number.
+                                     */    
+    static Point<spacedim> get_axis_vector (const unsigned int axis);
 };
 
 
-//  
-//  /**
-//   * Boundary object for the hull of a cylinder.  In three dimensions,
-//   * points are projected on a circular tube along the <tt>x-</tt>,
-//   * <tt>y-</tt> or <tt>z</tt>-axis (when using the first constructor of
-//   * this class), or an arbitrarily oriented cylinder described by the
-//   * direction of its axis and a point located on the axis. The radius
-//   * of the tube can be given independently. Similar to
-//   * HyperBallBoundary, new points are projected by dividing the
-//   * straight line between the old two points and adjusting the radius
-//   * from the axis.
-//   *
-//   * This class was developed to be used in conjunction with the
-//   * @p cylinder function of GridGenerator. It should be used for
-//   * the hull of the cylinder only (boundary indicator 0). Its use is
-//   * discussed in detail in the results section of step-49.
-//   *
-//   * This class is derived from StraightBoundary rather than from
-//   * Boundary, which would seem natural, since this way we can use the
-//   * StraightBoundary::in_between() function.
-//   *
-//   * @ingroup boundary
-//   *
-//   * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007
-//   */
-//  template <int dim, int spacedim = dim>
-//  class CylinderBoundary : public StraightBoundary<dim,spacedim>
-//  {
-//  public:
-//    /**
-//     * Constructor. Using default values for the constructor arguments yields a
-//     * circular tube along the x-axis (<tt>axis=0</tt>). Choose <tt>axis=1</tt>
-//     * or <tt>axis=2</tt> for a tube along the y- or z-axis, respectively.
-//     */
-//    CylinderBoundary (const double radius = 1.0,
-//                      const unsigned int axis = 0);
-//  
-//    /**
-//     * Constructor. If constructed
-//     * with this constructor, the
-//     * boundary described is a
-//     * cylinder with an axis that
-//     * points in direction #direction
-//     * and goes through the given
-//     * #point_on_axis. The direction
-//     * may be arbitrarily scaled, and
-//     * the given point may be any
-//     * point on the axis.
-//     */
-//    CylinderBoundary (const double           radius,
-//                      const Point<spacedim> &direction,
-//                      const Point<spacedim> &point_on_axis);
-//  
-//    /**
-//     * Refer to the general documentation of
-//     * this class and the documentation of the
-//     * base class.
-//     */
-//    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_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.
-//     *
-//     * Calls
-//     * @p get_intermediate_points_between_points.
-//     */
-//    virtual void
-//    get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
-//                                     std::vector<Point<spacedim> > &points) const;
-//  
-//    /**
-//     * Refer to the general
-//     * documentation of this class
-//     * and the documentation of the
-//     * base class.
-//     *
-//     * Only implemented for <tt>dim=3</tt>
-//     * and for <tt>points.size()==1</tt>.
-//     */
-//    virtual void
-//    get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
-//                                     std::vector<Point<spacedim> > &points) const;
-//  
-//    /**
-//     * Compute the normals to the
-//     * boundary at the vertices of
-//     * the given face.
-//     *
-//     * Refer to the general
-//     * documentation of this class
-//     * and the documentation of the
-//     * base class.
-//     */
-//    virtual void
-//    get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-//                             typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
-//  
-//    /**
-//     * Return the radius of the cylinder.
-//     */
-//    double get_radius () const;
-//  
-//    /**
-//     * Exception. Thrown by the
-//     * @p get_radius if the
-//     * @p compute_radius_automatically,
-//     * see below, flag is set true.
-//     */
-//    DeclException0 (ExcRadiusNotSet);
-//  
-//  
-//  protected:
-//    /**
-//     * Radius of the cylinder.
-//     */
-//    const double radius;
-//  
-//    /**
-//     * The direction vector of the axis.
-//     */
-//    const Point<spacedim> direction;
-//  
-//    /**
-//     * An arbitrary point on the axis.
-//     */
-//    const Point<spacedim> point_on_axis;
-//  
-//  private:
-//  
-//    /**
-//     * Called by
-//     * @p get_intermediate_points_on_line
-//     * and by
-//     * @p get_intermediate_points_on_quad.
-//     *
-//     * Refer to the general
-//     * documentation of
-//     * @p get_intermediate_points_on_line
-//     * in the documentation of the
-//     * base class.
-//     */
-//    void get_intermediate_points_between_points (const Point<spacedim> &p0, const Point<spacedim> &p1,
-//                                                 std::vector<Point<spacedim> > &points) const;
-//  
-//    /**
-//     * Given a number for the axis,
-//     * return a vector that denotes
-//     * this direction.
-//     */
-//    static Point<spacedim> get_axis_vector (const unsigned int axis);
-//  };
 //  
 //  
 //  
index e7d0d573256ed4c4cc35e29ebadcd3b75c4d742a..14fade0c7b3c9a01b9db50e3e41cb1d3765817e4 100644 (file)
@@ -1939,8 +1939,6 @@ template <int structdim, int dim, int spacedim>
 const Boundary<dim,spacedim> &
 TriaAccessor<structdim, dim, spacedim>::get_boundary () const
 {
-                                  // This assert should be removed in
-                                  // future releases...
   Assert (structdim<spacedim, ExcImpossibleInDim(dim));
   Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
                                   // Get the default (manifold_id)
@@ -1968,7 +1966,17 @@ TriaAccessor<structdim, dim, spacedim>::get_manifold () const
   Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
                                   // Get the default (manifold_id)
   types::manifold_id mi = this->objects().manifold_id[this->present_index];
+  
+                                  // In case this is not valid, check
+                                  // the boundary id, after having
+                                  // casted it to a manifold id
+  if(mi == numbers::invalid_manifold_id)
+    mi =  static_cast<types::manifold_id>
+         (structdim < dim ?
+          this->objects().boundary_or_material_id[this->present_index].boundary_id:
+          this->objects().boundary_or_material_id[this->present_index].material_id);
 
+                                  // Now recover the manifold itself
   return this->tria->get_manifold(mi);
 }
 
index 5df357c341b376fe28880d96370f51d44b78c1bb..a7b5ba4a3258b1df8aadc3c4eb61e95c3b16f1b8 100644 (file)
@@ -49,8 +49,85 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
   return p; 
 }
 
+// ============================================================
+// CylinderManifold
+// ============================================================
 
+template <int spacedim>
+CylinderManifold<spacedim>::CylinderManifold(const int axis):
+               direction (get_axis_vector (axis)),
+               point_on_axis (Point<spacedim>())
+{}
+
+template <int spacedim>
+CylinderManifold<spacedim>::CylinderManifold(const Point<spacedim> &direction,
+                                            const Point<spacedim> &point_on_axis):
+               direction (direction),
+               point_on_axis (point_on_axis)
+{}
+
+
+
+template <int spacedim>
+Point<spacedim>
+CylinderManifold<spacedim>::get_axis_vector (const unsigned int axis)
+{
+
+  Assert (axis < spacedim, ExcIndexRange (axis, 0, spacedim));
+  
+  Point<spacedim> axis_vector;
+
+  axis_vector[axis] = 1;
+
+  return axis_vector;
+}
+
+
+template <int spacedim>
+Point<spacedim>
+CylinderManifold<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()));
+
+                                  // compute a proposed new point  
+  Point<spacedim> middle = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+
+  double radius = 0;
+  Point<spacedim> on_plane;
+  
+  for(unsigned int i=0; i<surrounding_points.size(); ++i)
+    {
+      on_plane = surrounding_points[i]-point_on_axis;
+      on_plane = on_plane - (on_plane*direction) * direction;
+      radius += weights[i]*on_plane.norm();
+    }
+  
+                                  // we then have to project this
+                                  // point out to the given radius
+                                  // from the axis. to this end, we
+                                  // have to take into account the
+                                  // offset point_on_axis and the
+                                  // direction of the axis
+  const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
+                                          ((middle-point_on_axis) * direction) * direction;
+
+                                  // scale it to the desired length
+                                  // and put everything back
+                                  // together, unless we have a point
+                                  // on the axis
+  if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+    return middle;
+
+  else
+    return (vector_from_axis / vector_from_axis.norm() * radius +
+           ((middle-point_on_axis) * direction) * direction +
+           point_on_axis);
+}
 
+               
 // explicit instantiations
 #include "manifold_lib.inst"
 
index 2a5ddfc97891f406f21e2b93f63ca0944e0c77b2..1191c0d1fed1a548ffc76fc577cc06acc0e61f64 100644 (file)
 
 for (deal_II_dimension : DIMENSIONS)
   {
-//    template class PolarManifold<deal_II_dimension>;
+#if deal_II_dimension != 1
+    template class PolarManifold<deal_II_dimension>;
+    template class CylinderManifold<deal_II_dimension>;
+#endif
 //     template class ConeBoundary<deal_II_dimension>;
-//     template class HyperBallBoundary<deal_II_dimension>;
-//     template class HalfHyperBallBoundary<deal_II_dimension>;
-//     template class HyperShellBoundary<deal_II_dimension>;
-//     template class HalfHyperShellBoundary<deal_II_dimension>;
-// 
-// #if deal_II_dimension != 3
-//     template class HyperBallBoundary<deal_II_dimension,deal_II_dimension+1>;
-// #endif
-// #if deal_II_dimension == 3
-//    template class CylinderBoundary<deal_II_dimension-1,deal_II_dimension>;
-// #endif
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.