<commit_before>//---------------------------------------------------------------------------
//    $Id$
//    Version: $Name$
//
//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 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.
//
//---------------------------------------------------------------------------


#include <grid/tria_boundary_lib.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
#include <base/tensor.h>
#include <cmath>

DEAL_II_NAMESPACE_OPEN





template <int dim>
CylinderBoundary<dim>::CylinderBoundary (const double radius,
					 const unsigned int axis)
		:
		radius(radius),
		direction (get_axis_vector (axis)),
		point_on_axis (Point<dim>())
{}


template <int dim>
CylinderBoundary<dim>::CylinderBoundary (const double     radius,
					 const Point<dim> direction,
					 const Point<dim> point_on_axis)
		:
		radius(radius),
		direction (direction / direction.norm()),
		point_on_axis (point_on_axis)
{}


template <int dim>
Point<dim>
CylinderBoundary<dim>::get_axis_vector (const unsigned int axis)
{
  Assert (axis < dim, ExcIndexRange (axis, 0, dim));

  Point<dim> axis_vector;
  axis_vector[axis] = 1;
  return axis_vector;
}



template <int dim>
Point<dim>
CylinderBoundary<dim>::
get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
{
				   // compute a proposed new point
  const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);

				   // 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<dim> 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);
}


#if deal_II_dimension >=3

template<>
Point<3>
CylinderBoundary<3>::
get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
{
  const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);

				   // same algorithm as above
  const unsigned int dim = 3;
  const Point<dim> vector_from_axis = (middle-point_on_axis) -
				      ((middle-point_on_axis) * direction) * direction;
  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);
}

#endif


template <int dim>
Point<dim>
CylinderBoundary<dim>::
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
{
  Assert (false, ExcImpossibleInDim(dim));
  return Point<dim>();
}



template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_on_line (
  const typename Triangulation<dim>::line_iterator &line,
  std::vector<Point<dim> > &points) const
{
  if (points.size()==1)
    points[0]=get_new_point_on_line(line);
  else
    get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
}


template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_between_points (
  const Point<dim> &v0,
  const Point<dim> &v1,
  std::vector<Point<dim> > &points) const
{
  const unsigned int n=points.size();
  Assert(n>0, ExcInternalError());

				   // Do a simple linear interpolation
				   // followed by projection, using
				   // the same algorithm as above
  const Point<dim> ds = (v1-v0) / (n+1);

  for (unsigned int i=0; i<n; ++i)
    {
      const Point<dim> middle = v0 + (i+1)*ds;

      const Point<dim> vector_from_axis = (middle-point_on_axis) -
					  ((middle-point_on_axis) * direction) * direction;
      if (vector_from_axis.norm() <= 1e-10 * middle.norm())
	points[i] = middle;
      else
	points[i] = (vector_from_axis / vector_from_axis.norm() * radius +
		     ((middle-point_on_axis) * direction) * direction +
		     point_on_axis);
    }
}


#if deal_II_dimension == 3

template <>
void
CylinderBoundary<3>::get_intermediate_points_on_quad (
  const Triangulation<3>::quad_iterator &quad,
  std::vector<Point<3> > &points) const
{
  if (points.size()==1)
    points[0]=get_new_point_on_quad(quad);
  else
    {
      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
      Assert(points.size()==m*m, ExcInternalError());

      std::vector<Point<3> > lp0(m);
      std::vector<Point<3> > lp1(m);

      get_intermediate_points_on_line(quad->line(0), lp0);
      get_intermediate_points_on_line(quad->line(1), lp1);

      std::vector<Point<3> > lps(m);
      for (unsigned int i=0; i<m; ++i)
	{
	  get_intermediate_points_between_points(lp0[i], lp1[i], lps);

	  for (unsigned int j=0; j<m; ++j)
	    points[i*m+j]=lps[j];
	}
    }
}

#endif


template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_on_quad (
  const typename Triangulation<dim>::quad_iterator &,
  std::vector<Point<dim> > &) const
{
  Assert (false, ExcImpossibleInDim(dim));
}



#if deal_II_dimension == 1

template <>
void
CylinderBoundary<1>::
get_normals_at_vertices (const Triangulation<1>::face_iterator &,
			 Boundary<1>::FaceVertexNormals &) const
{
  Assert (false, ExcImpossibleInDim(1));
}

#endif


template <int dim>
void
CylinderBoundary<dim>::
get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
			 typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
{
  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
    {
      const Point<dim> vertex = face->vertex(v);

      const Point<dim> vector_from_axis = (vertex-point_on_axis) -
					  ((vertex-point_on_axis) * direction) * direction;

      face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm());
    }
}



template <int dim>
double
CylinderBoundary<dim>::get_radius () const
{
  return radius;
}


//======================================================================//

template<int dim>
ConeBoundary<dim>::ConeBoundary (const double radius_0,
				 const double radius_1,
				 const Point<dim> x_0,
				 const Point<dim> x_1)
		:
		radius_0 (radius_0),
		radius_1 (radius_1),
		x_0 (x_0),
		x_1 (x_1)
{}



template<int dim>
double ConeBoundary<dim>::get_radius (Point<dim> x) const
{
  for (unsigned int i = 0; i < dim; ++i)
    if ((x_1 (i) - x_0 (i)) != 0)
      return (radius_1 - radius_0) * (x (i) - x_0 (i)) / (x_1 (i) - x_0 (i)) + radius_0;

  return 0;
}



template<int dim>
void
ConeBoundary<dim>::
get_intermediate_points_between_points (const Point<dim> &p0,
					const Point<dim> &p1,
					std::vector<Point<dim> > &points) const
{
  const unsigned int n = points.size ();
  const Point<dim> axis = x_1 - x_0;
  const Point<dim> dx = (p1 - p0) / (n + 1);

  Assert (n > 0, ExcInternalError ());

  for (unsigned int i = 0; i < n; ++i)
    {
				       // Compute the current point.
      const Point<dim> x_i = p0 + (i + 1) * dx;
				       // To project this point on the
				       // boundary of the cone we first
				       // compute the orthogonal
				       // projection of this point onto
				       // the axis of the cone.
      const double c = (x_i - x_0) * axis / axis.square ();
      const Point<dim> x_ip = x_0 + c * axis;
				       // Compute the projection of
				       // the middle point on the
				       // boundary of the cone.
      points[i] = x_ip + get_radius (x_ip) *  (x_i - x_ip) / (x_i - x_ip).norm ();
    }
}

template<int dim>
Point<dim>
ConeBoundary<dim>::
get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
{
  const Point<dim> axis = x_1 - x_0;
				   // Compute the middle point of the line.
  const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
				   // To project it on the boundary of
				   // the cone we first compute the
				   // orthogonal projection of the
				   // middle point onto the axis of
				   // the cone.
  const double c = (middle - x_0) * axis / axis.square ();
  const Point<dim> middle_p = x_0 + c * axis;
				   // Compute the projection of the
				   // middle point on the boundary
				   // of the cone.
  return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
}



template<>
Point<3>
ConeBoundary<3>::get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
{
  const Point<3> axis = x_1 - x_0;
				   // Compute the middle point of the
				   // quad.
  const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
				   // Same algorithm as above: To
				   // project it on the boundary of
				   // the cone we first compute the
				   // orthogonal projection of the
				   // middle point onto the axis of
				   // the cone.
  const double c = (middle - x_0) * axis / axis.square ();
  const Point<3> middle_p = x_0 + c * axis;
				   // Compute the projection of the
				   // middle point on the boundary
				   // of the cone.
  return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
}



template<int dim>
Point<dim>
ConeBoundary<dim>::
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
{
  Assert (false, ExcImpossibleInDim (dim));

  return Point<dim>();
}



template<int dim>
void
ConeBoundary<dim>::
get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
				 std::vector<Point<dim> > &points) const
{
  if (points.size () == 1)
    points[0] = get_new_point_on_line (line);
  else
    get_intermediate_points_between_points (line->vertex (0), line->vertex (1), points);
}<commit_msg>remove reference to 3D in 2D<commit_after>//---------------------------------------------------------------------------
//    $Id$
//    Version: $Name$
//
//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 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.
//
//---------------------------------------------------------------------------


#include <grid/tria_boundary_lib.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
#include <base/tensor.h>
#include <cmath>

DEAL_II_NAMESPACE_OPEN





template <int dim>
CylinderBoundary<dim>::CylinderBoundary (const double radius,
					 const unsigned int axis)
		:
		radius(radius),
		direction (get_axis_vector (axis)),
		point_on_axis (Point<dim>())
{}


template <int dim>
CylinderBoundary<dim>::CylinderBoundary (const double     radius,
					 const Point<dim> direction,
					 const Point<dim> point_on_axis)
		:
		radius(radius),
		direction (direction / direction.norm()),
		point_on_axis (point_on_axis)
{}


template <int dim>
Point<dim>
CylinderBoundary<dim>::get_axis_vector (const unsigned int axis)
{
  Assert (axis < dim, ExcIndexRange (axis, 0, dim));

  Point<dim> axis_vector;
  axis_vector[axis] = 1;
  return axis_vector;
}



template <int dim>
Point<dim>
CylinderBoundary<dim>::
get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
{
				   // compute a proposed new point
  const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);

				   // 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<dim> 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);
}


#if deal_II_dimension >=3

template<>
Point<3>
CylinderBoundary<3>::
get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
{
  const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);

				   // same algorithm as above
  const unsigned int dim = 3;
  const Point<dim> vector_from_axis = (middle-point_on_axis) -
				      ((middle-point_on_axis) * direction) * direction;
  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);
}

#endif


template <int dim>
Point<dim>
CylinderBoundary<dim>::
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
{
  Assert (false, ExcImpossibleInDim(dim));
  return Point<dim>();
}



template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_on_line (
  const typename Triangulation<dim>::line_iterator &line,
  std::vector<Point<dim> > &points) const
{
  if (points.size()==1)
    points[0]=get_new_point_on_line(line);
  else
    get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
}


template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_between_points (
  const Point<dim> &v0,
  const Point<dim> &v1,
  std::vector<Point<dim> > &points) const
{
  const unsigned int n=points.size();
  Assert(n>0, ExcInternalError());

				   // Do a simple linear interpolation
				   // followed by projection, using
				   // the same algorithm as above
  const Point<dim> ds = (v1-v0) / (n+1);

  for (unsigned int i=0; i<n; ++i)
    {
      const Point<dim> middle = v0 + (i+1)*ds;

      const Point<dim> vector_from_axis = (middle-point_on_axis) -
					  ((middle-point_on_axis) * direction) * direction;
      if (vector_from_axis.norm() <= 1e-10 * middle.norm())
	points[i] = middle;
      else
	points[i] = (vector_from_axis / vector_from_axis.norm() * radius +
		     ((middle-point_on_axis) * direction) * direction +
		     point_on_axis);
    }
}


#if deal_II_dimension == 3

template <>
void
CylinderBoundary<3>::get_intermediate_points_on_quad (
  const Triangulation<3>::quad_iterator &quad,
  std::vector<Point<3> > &points) const
{
  if (points.size()==1)
    points[0]=get_new_point_on_quad(quad);
  else
    {
      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
      Assert(points.size()==m*m, ExcInternalError());

      std::vector<Point<3> > lp0(m);
      std::vector<Point<3> > lp1(m);

      get_intermediate_points_on_line(quad->line(0), lp0);
      get_intermediate_points_on_line(quad->line(1), lp1);

      std::vector<Point<3> > lps(m);
      for (unsigned int i=0; i<m; ++i)
	{
	  get_intermediate_points_between_points(lp0[i], lp1[i], lps);

	  for (unsigned int j=0; j<m; ++j)
	    points[i*m+j]=lps[j];
	}
    }
}

#endif


template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_on_quad (
  const typename Triangulation<dim>::quad_iterator &,
  std::vector<Point<dim> > &) const
{
  Assert (false, ExcImpossibleInDim(dim));
}



#if deal_II_dimension == 1

template <>
void
CylinderBoundary<1>::
get_normals_at_vertices (const Triangulation<1>::face_iterator &,
			 Boundary<1>::FaceVertexNormals &) const
{
  Assert (false, ExcImpossibleInDim(1));
}

#endif


template <int dim>
void
CylinderBoundary<dim>::
get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
			 typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
{
  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
    {
      const Point<dim> vertex = face->vertex(v);

      const Point<dim> vector_from_axis = (vertex-point_on_axis) -
					  ((vertex-point_on_axis) * direction) * direction;

      face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm());
    }
}



template <int dim>
double
CylinderBoundary<dim>::get_radius () const
{
  return radius;
}


//======================================================================//

template<int dim>
ConeBoundary<dim>::ConeBoundary (const double radius_0,
				 const double radius_1,
				 const Point<dim> x_0,
				 const Point<dim> x_1)
		:
		radius_0 (radius_0),
		radius_1 (radius_1),
		x_0 (x_0),
		x_1 (x_1)
{}



template<int dim>
double ConeBoundary<dim>::get_radius (Point<dim> x) const
{
  for (unsigned int i = 0; i < dim; ++i)
    if ((x_1 (i) - x_0 (i)) != 0)
      return (radius_1 - radius_0) * (x (i) - x_0 (i)) / (x_1 (i) - x_0 (i)) + radius_0;

  return 0;
}



template<int dim>
void
ConeBoundary<dim>::
get_intermediate_points_between_points (const Point<dim> &p0,
					const Point<dim> &p1,
					std::vector<Point<dim> > &points) const
{
  const unsigned int n = points.size ();
  const Point<dim> axis = x_1 - x_0;
  const Point<dim> dx = (p1 - p0) / (n + 1);

  Assert (n > 0, ExcInternalError ());

  for (unsigned int i = 0; i < n; ++i)
    {
				       // Compute the current point.
      const Point<dim> x_i = p0 + (i + 1) * dx;
				       // To project this point on the
				       // boundary of the cone we first
				       // compute the orthogonal
				       // projection of this point onto
				       // the axis of the cone.
      const double c = (x_i - x_0) * axis / axis.square ();
      const Point<dim> x_ip = x_0 + c * axis;
				       // Compute the projection of
				       // the middle point on the
				       // boundary of the cone.
      points[i] = x_ip + get_radius (x_ip) *  (x_i - x_ip) / (x_i - x_ip).norm ();
    }
}

template<int dim>
Point<dim>
ConeBoundary<dim>::
get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
{
  const Point<dim> axis = x_1 - x_0;
				   // Compute the middle point of the line.
  const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
				   // To project it on the boundary of
				   // the cone we first compute the
				   // orthogonal projection of the
				   // middle point onto the axis of
				   // the cone.
  const double c = (middle - x_0) * axis / axis.square ();
  const Point<dim> middle_p = x_0 + c * axis;
				   // Compute the projection of the
				   // middle point on the boundary
				   // of the cone.
  return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
}


#if deal_II_dimension == 3

template<int dim>
Point<dim>
ConeBoundary<dim>::get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
{
  const Point<dim> axis = x_1 - x_0;
				   // Compute the middle point of the
				   // quad.
  const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
				   // Same algorithm as above: To
				   // project it on the boundary of
				   // the cone we first compute the
				   // orthogonal projection of the
				   // middle point onto the axis of
				   // the cone.
  const double c = (middle - x_0) * axis / axis.square ();
  const Point<dim> middle_p = x_0 + c * axis;
				   // Compute the projection of the
				   // middle point on the boundary
				   // of the cone.
  return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
}

#else

template<int dim>
Point<dim>
ConeBoundary<dim>::
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
{
  Assert (false, ExcImpossibleInDim (dim));

  return Point<dim>();
}

#endif

template<int dim>
void
ConeBoundary<dim>::
get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
				 std::vector<Point<dim> > &points) const
{
  if (points.size () == 1)
    points[0] = get_new_point_on_line (line);
  else
    get_intermediate_points_between_points (line->vertex (0), line->vertex (1), points);
}
{'fraction_non_alphanumeric': 0.1484949504853417, 'fraction_numerical': 0.014707324247475242, 'mean_word_length': 3.1844102564102563, 'pattern_counts': {'":': 0, '<': 271, '<?xml version=': 0, '>': 271, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5929020', 'n_tokens_mistral': 7631, 'n_tokens_neox': 6824, 'n_words': 1761}