16 #include <deal.II/grid/tria_boundary_lib.h> 17 #include <deal.II/grid/tria.h> 18 #include <deal.II/grid/tria_iterator.h> 19 #include <deal.II/grid/tria_accessor.h> 20 #include <deal.II/base/tensor.h> 25 DEAL_II_NAMESPACE_OPEN
29 template <
int dim,
int spacedim>
31 const unsigned int axis)
34 direction (get_axis_vector (axis)),
35 point_on_axis (
Point<spacedim>())
39 template <
int dim,
int spacedim>
45 direction (direction / direction.norm()),
46 point_on_axis (point_on_axis)
50 template <
int dim,
int spacedim>
57 axis_vector[axis] = 1;
63 template <
int dim,
int spacedim>
83 if (vector_from_axis.
norm() <= 1e-10 * middle.
norm())
101 const unsigned int spacedim = 3;
105 if (vector_from_axis.
norm() <= 1e-10 * middle.
norm())
121 const unsigned int spacedim = 3;
124 if (vector_from_axis.
norm() <= 1e-10 * middle.
norm())
133 template <
int dim,
int spacedim>
144 template <
int dim,
int spacedim>
147 const typename Triangulation<dim,spacedim>::line_iterator &line,
157 template <
int dim,
int spacedim>
164 const unsigned int n=
points.size();
171 for (
unsigned int i=0; i<n; ++i)
173 const double x = line_points[i+1][0];
178 if (vector_from_axis.
norm() <= 1e-10 * middle.
norm())
192 const Triangulation<3>::quad_iterator &quad,
199 unsigned int m=
static_cast<unsigned int> (std::sqrt(static_cast<double>(
points.size())));
202 std::vector<Point<3> > lp0(m);
203 std::vector<Point<3> > lp1(m);
208 std::vector<Point<3> > lps(m);
209 for (
unsigned int i=0; i<m; ++i)
213 for (
unsigned int j=0; j<m; ++j)
221 template <
int dim,
int spacedim>
224 const typename Triangulation<dim,spacedim>::quad_iterator &,
245 template <
int dim,
int spacedim>
251 for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
258 face_vertex_normals[v] = (vector_from_axis / vector_from_axis.
norm());
264 template <
int dim,
int spacedim>
276 const double radius_1,
291 for (
unsigned int i = 0; i < dim; ++i)
292 if ((
x_1 (i) -
x_0 (i)) != 0)
307 const unsigned int n =
points.size ();
314 for (
unsigned int i=0; i<n; ++i)
316 const double x = line_points[i+1][0];
322 const double c = (x_i -
x_0) * axis / (axis*axis);
340 const double c = (middle -
x_0) * axis / (axis*axis);
343 return middle_p +
get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
361 const double c = (middle -
x_0) * axis / (axis*axis);
364 return middle_p +
get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
406 unsigned int n =
static_cast<unsigned int> (std::sqrt (static_cast<double> (
points.size ())));
410 std::vector<Point<3> > points_line_0 (n);
411 std::vector<Point<3> > points_line_1 (n);
416 std::vector<Point<3> > points_line_segment (n);
418 for (
unsigned int i = 0; i < n; ++i)
422 points_line_segment);
424 for (
unsigned int j = 0; j < n; ++j)
425 points[i * n + j] = points_line_segment[j];
463 for (
unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_face; ++vertex)
467 const double c = (face->vertex (vertex) -
x_0) * axis / (axis*axis);
471 const Tensor<1,dim> axis_to_vertex = face->vertex (vertex) - vertex_p;
473 face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.
norm ();
480 template <
int dim,
int spacedim>
486 compute_radius_automatically(false)
491 template <
int dim,
int spacedim>
501 r = (line->vertex(0) -
center).norm();
506 middle *= r / std::sqrt(middle.
square());
534 template <
int dim,
int spacedim>
545 r = (quad->vertex(0) -
center).norm();
550 middle *= r / std::sqrt(middle.
square());
561 const Triangulation<1>::line_iterator &,
569 template <
int dim,
int spacedim>
572 const typename Triangulation<dim,spacedim>::line_iterator &line,
583 template <
int dim,
int spacedim>
589 const unsigned int n=
points.size();
594 const double length=(v1-v0).norm();
604 Assert((std::fabs(v0*v0-r*r)<eps*r*r)
606 (std::fabs(v1*v1-r*r)<eps*r*r),
607 ExcMessage(
"In computing the location of midpoints of an edge of a " 608 "HyperBallBoundary, one of the vertices of the edge " 609 "does not have the expected distance from the center " 610 "of the sphere. This may happen if the geometry has " 611 "been deformed, or if the HyperBallBoundary object has " 612 "been assigned to a part of the boundary that is not " 613 "in fact a sphere (e.g., the sides of a quarter shell " 614 "or one octant of a ball)."));
616 const double alpha=std::acos((v0*v1)/std::sqrt((v0*v0)*(v1*v1)));
619 const double h=pm.
norm();
624 const unsigned int m=n/2;
625 for (
unsigned int i=0; i<m ; ++i)
627 const double beta = alpha * (line_points[i+1][0]-0.5);
628 const double d = h*std::tan(beta);
639 for (
unsigned int i=0; i<n; ++i)
651 const Triangulation<3>::quad_iterator &quad,
658 unsigned int m=
static_cast<unsigned int> (std::sqrt(static_cast<double>(
points.size())));
661 std::vector<Point<3> > lp0(m);
662 std::vector<Point<3> > lp1(m);
667 std::vector<Point<3> > lps(m);
668 for (
unsigned int i=0; i<m; ++i)
672 for (
unsigned int j=0; j<m; ++j)
683 const Triangulation<2,3>::quad_iterator &quad,
690 unsigned int m=
static_cast<unsigned int> (std::sqrt(static_cast<double>(
points.size())));
693 std::vector<Point<3> > lp0(m);
694 std::vector<Point<3> > lp1(m);
699 std::vector<Point<3> > lps(m);
700 for (
unsigned int i=0; i<m; ++i)
704 for (
unsigned int j=0; j<m; ++j)
712 template <
int dim,
int spacedim>
715 const typename Triangulation<dim,spacedim>::quad_iterator &,
723 template <
int dim,
int spacedim>
730 return unnormalized_normal/unnormalized_normal.
norm();
755 template <
int dim,
int spacedim>
761 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
762 face_vertex_normals[vertex] = face->vertex(vertex)-
center;
767 template <
int dim,
int spacedim>
776 template <
int dim,
int spacedim>
804 const Point<dim> line_center = line->center();
805 const Point<dim> vertices[2] = { line->vertex(0), line->vertex(1) };
807 if ((line_center(0) == this->
center(0))
809 ((std::fabs(vertices[0].distance(this->
center)-this->
radius) >
812 (std::fabs(vertices[1].distance(this->
center)-this->
radius) >
837 const Point<dim> quad_center = quad->center();
838 if (quad_center(0) == this->
center(0))
854 const Point<dim> line_center = line->center();
855 if (line_center(0) == this->
center(0))
875 const Point<dim> quad_center = quad->center();
876 if (quad_center(0) == this->
center(0))
915 const Point<dim> quad_center = face->center();
916 if (quad_center(0) == this->
center(0))
943 const double inner_radius,
944 const double outer_radius)
947 inner_radius (inner_radius),
948 outer_radius (outer_radius)
951 Assert ((inner_radius >= 0) &&
952 (outer_radius > 0) &&
953 (outer_radius > inner_radius),
954 ExcMessage (
"Inner and outer radii must be specified explicitly in 3d."));
970 if ((line->vertex(0)(0) == this->
center(0))
972 (line->vertex(1)(0) == this->
center(0)))
973 return (line->vertex(0) + line->vertex(1))/2;
985 if (((line->vertex(0)(0) == this->
center(0))
987 (line->vertex(1)(0) == this->
center(0)))
989 !(((std::fabs (line->vertex(0).distance (this->
center)
992 (std::fabs (line->vertex(1).distance (this->
center)
995 ((std::fabs (line->vertex(0).distance (this->
center)
996 - outer_radius) < 1e-12 * outer_radius)
998 (std::fabs (line->vertex(1).distance (this->
center)
999 - outer_radius) < 1e-12 * outer_radius))))
1000 return (line->vertex(0) + line->vertex(1))/2;
1035 if ((quad->vertex(0)(0) == this->
center(0)) &&
1036 (quad->vertex(1)(0) == this->
center(0)) &&
1037 (quad->vertex(2)(0) == this->
center(0)) &&
1038 (quad->vertex(3)(0) == this->
center(0)))
1040 const Point<dim> quad_center = (quad->vertex(0) + quad->vertex(1) +
1041 quad->vertex(2) + quad->vertex(3) )/4;
1045 if (std::fabs (quad->line(0)->center().distance(this->center) -
1046 quad->line(1)->center().distance(this->center))
1047 < 1e-12 * outer_radius)
1050 const double needed_radius
1051 = quad->line(0)->center().distance(this->center);
1053 return (this->center +
1054 quad_center_offset/quad_center_offset.
norm() * needed_radius);
1056 else if (std::fabs (quad->line(2)->center().distance(this->center) -
1057 quad->line(3)->center().distance(this->center))
1058 < 1e-12 * outer_radius)
1061 const double needed_radius
1062 = quad->line(2)->center().distance(this->center);
1064 return (this->center +
1065 quad_center_offset/quad_center_offset.
norm() * needed_radius);
1090 if ((line->vertex(0)(0) == this->
center(0))
1092 (line->vertex(1)(0) == this->
center(0)))
1105 if (((line->vertex(0)(0) == this->
center(0))
1107 (line->vertex(1)(0) == this->
center(0)))
1109 !(((std::fabs (line->vertex(0).distance (this->
center)
1112 (std::fabs (line->vertex(1).distance (this->
center)
1115 ((std::fabs (line->vertex(0).distance (this->
center)
1116 - outer_radius) < 1e-12 * outer_radius)
1118 (std::fabs (line->vertex(1).distance (this->
center)
1119 - outer_radius) < 1e-12 * outer_radius))))
1146 const Point<dim> quad_center = quad->center();
1147 if (quad_center(0) == this->
center(0))
1185 if (face->center()(0) == this->
center(0))
1194 template <
int dim,
int spacedim>
1218 template <
int dim,
int spacedim>
1222 const double y)
const 1244 const double theta=surfP(0);
1245 const double phi=surfP(1);
1247 return Point<3> ((
R+r*std::cos(phi))*std::cos(theta),
1249 (
R+r*std::cos(phi))*std::sin(theta));
1258 const double phi=std::asin(std::abs(p(1))/r);
1259 const double Rr_2=p(0)*p(0)+p(2)*p(2);
1264 if (std::abs(p(0))<1.E-5)
1273 const double theta = std::atan(std::abs(p(2)/p(0)));
1295 for (
unsigned int i=0; i<2; i++)
1315 for (
unsigned int i=0; i<4; i++)
1323 for (
unsigned int i=0; i<2; i++)
1324 for (
unsigned int j=1; j<4; j++)
1332 for (
unsigned int i=0; i<4; i++)
1349 double theta=surfP[0];
1350 double phi=surfP[1];
1352 double f=
R+r*std::cos(phi);
1354 n[0]=r*std::cos(phi)*std::cos(theta)*f;
1355 n[1]=r*std::sin(phi)*f;
1356 n[2]=r*std::sin(theta)*std::cos(phi)*f;
1383 unsigned int npoints=
points.size();
1384 if (npoints==0)
return;
1388 for (
unsigned int i=0; i<2; i++)
1391 unsigned int offset[2];
1399 for (
unsigned int i=0; i<2; i++)
1400 for (
unsigned int j=1; j<2; j++)
1409 for (
unsigned int i=0; i<2; i++)
1410 for (
unsigned int j=0; j<2; j++)
1411 if (p[j](i)<1.E-12 )
1417 for (
unsigned int i=0; i<npoints; i++)
1419 const double x = line_points[i+1][0];
1420 target= (1-x)*p[0] + x*p[1];
1434 const unsigned int n=
points.size(),
1435 m=
static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
1441 for (
unsigned int i=0; i<4; i++)
1445 unsigned int offset[2];
1453 for (
unsigned int i=0; i<2; i++)
1454 for (
unsigned int j=1; j<4; j++)
1463 for (
unsigned int i=0; i<2; i++)
1464 for (
unsigned int j=0; j<4; j++)
1465 if (p[j](i)<1.E-12 )
1469 for (
unsigned int i=0; i<m; ++i)
1471 const double y=line_points[i+1][0];
1472 for (
unsigned int j=0; j<m; ++j)
1474 const double x=line_points[j+1][0];
1475 target=((1-x) * p[0] +
1493 for (
unsigned int i=0; i<GeometryInfo<2>::vertices_per_face; i++)
1500 #include "tria_boundary_lib.inst" 1502 DEAL_II_NAMESPACE_CLOSE
virtual void get_intermediate_points_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line, std::vector< Point< spacedim > > &points) const
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad, std::vector< Point< spacedim > > &points) const
virtual Point< dim > get_new_point_on_quad(const typename Triangulation< dim >::quad_iterator &quad) const
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Point< spacedim > get_surf_norm_from_sp(const Point< dim > &surfP) const
const Point< spacedim > point_on_axis
void get_intermediate_points_between_points(const Point< spacedim > &p0, const Point< spacedim > &p1, std::vector< Point< spacedim > > &points) const
virtual void get_intermediate_points_on_line(const typename Triangulation< dim >::line_iterator &line, std::vector< Point< dim > > &points) const
const Point< spacedim > direction
Point< spacedim > get_center() const
virtual void get_intermediate_points_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line, std::vector< Point< spacedim > > &points) const
double get_radius(const Point< dim > x) const
virtual Point< dim > get_new_point_on_quad(const typename Triangulation< dim >::quad_iterator &quad) const
numbers::NumberTraits< Number >::real_type norm() const
virtual void get_normals_at_vertices(const typename Triangulation< dim >::face_iterator &face, typename Boundary< dim >::FaceVertexNormals &face_vertex_normals) const
bool compute_radius_automatically
Point< dim > get_surf_coord(const Point< spacedim > &p) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Boundary< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Boundary< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const
double get_radius() const
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
virtual Point< dim > get_new_point_on_line(const typename Triangulation< dim >::line_iterator &line) const
static Point< spacedim > get_axis_vector(const unsigned int axis)
HalfHyperShellBoundary(const Point< dim > ¢er=Point< dim >(), const double inner_radius=-1, const double outer_radius=-1)
static::ExceptionBase & ExcMessage(std::string arg1)
Point< spacedim > get_surf_norm(const Point< spacedim > &p) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Boundary< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const
std::vector< std_cxx11::shared_ptr< QGaussLobatto< 1 > > > points
static::ExceptionBase & ExcImpossibleInDim(int arg1)
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
HalfHyperBallBoundary(const Point< dim > p=Point< dim >(), const double radius=1.0)
#define Assert(cond, exc)
Point< spacedim > get_real_coord(const Point< dim > &surfP) const
double get_correct_angle(const double angle, const double x, const double y) const
virtual Point< dim > get_new_point_on_quad(const typename Triangulation< dim >::quad_iterator &quad) const
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
const std::vector< Point< 1 > > & get_line_support_points(const unsigned int n_intermediate_points) const
virtual Point< dim > get_new_point_on_line(const typename Triangulation< dim >::line_iterator &line) const
void get_intermediate_points_between_points(const Point< dim > &p0, const Point< dim > &p1, std::vector< Point< dim > > &points) const
void get_intermediate_points_between_points(const Point< spacedim > &p0, const Point< spacedim > &p1, std::vector< Point< spacedim > > &points) const
ConeBoundary(const double radius_0, const double radius_1, const Point< dim > x_0, const Point< dim > x_1)
const Point< spacedim > center
static::ExceptionBase & ExcRadiusNotSet()
numbers::NumberTraits< Number >::real_type square() const
HyperShellBoundary(const Point< dim > ¢er=Point< dim >())
virtual void get_normals_at_vertices(const typename Triangulation< dim >::face_iterator &face, typename Boundary< dim >::FaceVertexNormals &face_vertex_normals) const
virtual Point< dim > get_new_point_on_line(const typename Triangulation< dim >::line_iterator &line) const
virtual void get_intermediate_points_on_line(const typename Triangulation< dim >::line_iterator &line, std::vector< Point< dim > > &points) const
virtual void get_intermediate_points_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line, std::vector< Point< spacedim > > &points) const
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad, std::vector< Point< spacedim > > &points) const
TorusBoundary(const double R, const double r)
static::ExceptionBase & ExcNotImplemented()
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad, std::vector< Point< spacedim > > &points) const
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad, std::vector< Point< spacedim > > &points) const
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim >::quad_iterator &quad, std::vector< Point< dim > > &points) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Boundary< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const
const double inner_radius
HyperBallBoundary(const Point< spacedim > p=Point< spacedim >(), const double radius=1.0)
CylinderBoundary(const double radius=1.0, const unsigned int axis=0)
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
virtual void get_intermediate_points_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line, std::vector< Point< spacedim > > &points) const
virtual void get_normals_at_vertices(const typename Triangulation< dim >::face_iterator &face, typename Boundary< dim >::FaceVertexNormals &face_vertex_normals) const
virtual void get_intermediate_points_on_line(const typename Triangulation< dim >::line_iterator &line, std::vector< Point< dim > > &points) const
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim >::quad_iterator &quad, std::vector< Point< dim > > &points) const
static::ExceptionBase & ExcInternalError()
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
double get_radius() const
virtual void get_intermediate_points_on_quad(const typename Triangulation< dim >::quad_iterator &quad, std::vector< Point< dim > > &points) const