#include <deal.II/base/polynomial.h>
#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx17/algorithm.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/fe/fe_pyramid_p.h>
+namespace
+{
+ // Compute the nearest point to @p on the line segment and the square of its
+ // distance to @p.
+ template <int dim>
+ std::pair<Point<dim>, double>
+ project_to_line(const Point<dim> &x0,
+ const Point<dim> &x1,
+ const Point<dim> &p)
+ {
+ Assert(x0 != x1, ExcInternalError());
+ // t is the convex combination coefficient (x = (1 - t) * x0 + t * x1)
+ // defining the position of the closest point on the line (not line segment)
+ // to p passing through x0 and x1. This formula is equivalent to the
+ // standard 'project a vector onto another vector', where each vector is
+ // shifted to start at x0.
+ const double t = ((x1 - x0) * (p - x0)) / ((x1 - x0).norm_square());
+
+ // Only consider points between x0 and x1
+ if (t <= 0)
+ return std::make_pair(x0, x0.distance_square(p));
+ else if (t <= 1)
+ {
+ const auto p2 = x0 + t * (x1 - x0);
+ return std::make_pair(p2, p2.distance_square(p));
+ }
+ else
+ return std::make_pair(x1, x1.distance_square(p));
+ }
+
+ // template base-case
+ template <int dim>
+ std::pair<Point<dim>, double>
+ project_to_quad(const std::array<Point<dim>, 3> & /*vertices*/,
+ const Point<dim> & /*p*/,
+ const ReferenceCell /*reference_cell*/)
+ {
+ Assert(false, ExcInternalError());
+ return std::make_pair(Point<dim>(),
+ std::numeric_limits<double>::signaling_NaN());
+ }
+
+ /**
+ * Compute the nearest point on a quad (in the deal.II sense: i.e., something
+ * with structdim = 2 and spacedim = 3) and the square of that point's
+ * distance to @p p. Here, the quad is described with three vertices: either
+ * the three vertices of a Triangle or the first three of a Quadrilateral (as
+ * the fourth one can be computed, in that case, from the first three).
+ *
+ * If the given point cannot be projected via a normal vector (i.e., if the
+ * line parallel to the normal vector intersecting @p p does not intersect the
+ * quad) then this function returns the origin and the largest double
+ * precision number. distance_to_line_square() is in charge of computing the
+ * distance to lines.
+ *
+ * @note This function is for Quadrilaterals and Triangles because it is
+ * intended to find the shortest distance to the face of a reference cell
+ * (which must be a Triangle or Quadrilateral) in 3d.
+ */
+ template <>
+ std::pair<Point<3>, double>
+ project_to_quad(const std::array<Point<3>, 3> &vertices,
+ const Point<3> & p,
+ const ReferenceCell face_reference_cell)
+ {
+ Assert(face_reference_cell == ReferenceCells::Triangle ||
+ face_reference_cell == ReferenceCells::Quadrilateral,
+ ExcNotImplemented());
+
+ // Make the problem slightly easier by shifting everything to avoid a point
+ // at the origin (this way we can invert the matrix of vertices). Use 2.0 so
+ // that the bottom left vertex of a Pyramid is now at x = 1.
+ std::array<Point<3>, 3> shifted_vertices = vertices;
+ const Tensor<1, 3> shift{{2.0, 2.0, 2.0}};
+ for (Point<3> &shifted_vertex : shifted_vertices)
+ shifted_vertex += shift;
+ const Point<3> shifted_p = p + shift;
+
+ // As we are projecting onto a face of a reference cell, the vectors
+ // describing its local coordinate system should be orthogonal. We don't
+ // know which of the three vectors computed from p are mutually orthogonal
+ // for triangles so that case requires an extra check.
+ Tensor<1, 3> e0;
+ Tensor<1, 3> e1;
+ const Point<3> vertex = shifted_vertices[0];
+ // Triangles are difficult because of two cases:
+ // 1. the top face of a Tetrahedron, which does not have a right angle
+ // 2. wedges and pyramids, whose faces do not lie on the reference simplex
+ //
+ // Deal with both by creating a locally orthogonal (but not necessarily
+ // orthonormal) coordinate system and testing if the projected point is in
+ // the triangle by expressing it as a convex combination of the vertices.
+ if (face_reference_cell == ReferenceCells::Triangle)
+ {
+ e0 = shifted_vertices[1] - shifted_vertices[0];
+ e1 = shifted_vertices[2] - shifted_vertices[0];
+ e1 -= (e0 * e1) * e0 / (e0.norm_square());
+ }
+ else
+ {
+ e0 = shifted_vertices[1] - shifted_vertices[0];
+ e1 = shifted_vertices[2] - shifted_vertices[0];
+ }
+ Assert(std::abs(e0 * e1) <= 1e-14, ExcInternalError());
+ // the quadrilaterals on pyramids and wedges don't necessarily have edge
+ // lengths of 1 so we cannot skip the denominator
+ const double c0 = e0 * (shifted_p - vertex) / e0.norm_square();
+ const double c1 = e1 * (shifted_p - vertex) / e1.norm_square();
+ const Point<3> projected_shifted_p = vertex + c0 * e0 + c1 * e1;
+
+ bool in_quad = false;
+ if (face_reference_cell == ReferenceCells::Triangle)
+ {
+ Tensor<2, 3> shifted_vertex_matrix;
+ for (unsigned int i = 0; i < 3; ++i)
+ shifted_vertex_matrix[i] = shifted_vertices[i];
+ const Tensor<1, 3> combination_coordinates =
+ invert(transpose(shifted_vertex_matrix)) * projected_shifted_p;
+ bool is_convex_combination = true;
+ for (unsigned int i = 0; i < 3; ++i)
+ is_convex_combination = is_convex_combination &&
+ (0.0 <= combination_coordinates[i]) &&
+ (combination_coordinates[i] <= 1.0);
+ in_quad = is_convex_combination;
+ }
+ else
+ in_quad = (0.0 <= c0 && c0 <= 1.0 && 0.0 <= c1 && c1 <= 1.0);
+
+ if (in_quad)
+ return std::make_pair(projected_shifted_p - shift,
+ shifted_p.distance_square(projected_shifted_p));
+ else
+ return std::make_pair(Point<3>(), std::numeric_limits<double>::max());
+ }
+} // namespace
+
+
+
+template <int dim>
+Point<dim>
+ReferenceCell::closest_point(const Point<dim> &p) const
+{
+ AssertDimension(dim, get_dimension());
+
+ // Handle simple cases first:
+ if (dim == 0)
+ return p;
+ if (contains_point(p, 0.0))
+ return p;
+ if (dim == 1)
+ return project_to_line(vertex<dim>(0), vertex<dim>(1), p).first;
+
+ // Find the closest vertex so that we only need to check adjacent faces and
+ // lines.
+ Point<dim> result;
+ unsigned int closest_vertex_no = 0;
+ double closest_vertex_distance_square = vertex<dim>(0).distance_square(p);
+ for (unsigned int i = 1; i < n_vertices(); ++i)
+ {
+ const double new_vertex_distance_square =
+ vertex<dim>(i).distance_square(p);
+ if (new_vertex_distance_square < closest_vertex_distance_square)
+ {
+ closest_vertex_no = i;
+ closest_vertex_distance_square = new_vertex_distance_square;
+ }
+ }
+
+ double min_distance_square = std::numeric_limits<double>::max();
+ if (dim == 2)
+ {
+ for (const unsigned int face_no :
+ faces_for_given_vertex(closest_vertex_no))
+ {
+ const Point<dim> v0 = vertex<dim>(line_to_cell_vertices(face_no, 0));
+ const Point<dim> v1 = vertex<dim>(line_to_cell_vertices(face_no, 1));
+
+ auto pair = project_to_line(v0, v1, p);
+ if (pair.second < min_distance_square)
+ {
+ result = pair.first;
+ min_distance_square = pair.second;
+ }
+ }
+ }
+ else
+ {
+ // Check faces and then lines.
+ //
+ // For reference cells with sloped faces (i.e., all 3D shapes except
+ // Hexahedra), we might be able to do a valid normal projection to a face
+ // with a different slope which is on the 'other side' of the reference
+ // cell. To catch that case we have to unconditionally check lines after
+ // checking faces.
+ //
+ // For pyramids the closest vertex might not be on the closest face: for
+ // example, the origin is closest to vertex 4 which is not on the bottom
+ // plane. Get around that by just checking all faces for pyramids.
+ const std::array<unsigned int, 5> all_pyramid_faces{{0, 1, 2, 3, 4}};
+ const auto &faces = *this == ReferenceCells::Pyramid ?
+ ArrayView<const unsigned int>(all_pyramid_faces) :
+ faces_for_given_vertex(closest_vertex_no);
+ for (const unsigned int face_no : faces)
+ {
+ auto face_cell = face_reference_cell(face_no);
+ // We only need the first three points since for quads the last point
+ // is redundant
+ std::array<Point<dim>, 3> vertices;
+ for (unsigned int vertex_no = 0; vertex_no < 3; ++vertex_no)
+ vertices[vertex_no] = vertex<dim>(face_to_cell_vertices(
+ face_no, vertex_no, default_combined_face_orientation()));
+
+ auto pair = project_to_quad(vertices, p, face_cell);
+ if (pair.second < min_distance_square)
+ {
+ result = pair.first;
+ min_distance_square = pair.second;
+ }
+ }
+
+ for (const unsigned int face_no :
+ faces_for_given_vertex(closest_vertex_no))
+ {
+ auto face_cell = face_reference_cell(face_no);
+ for (const unsigned int face_line_no : face_cell.line_indices())
+ {
+ const auto cell_line_no =
+ face_to_cell_lines(face_no,
+ face_line_no,
+ default_combined_face_orientation());
+ const auto v0 =
+ vertex<dim>(line_to_cell_vertices(cell_line_no, 0));
+ const auto v1 =
+ vertex<dim>(line_to_cell_vertices(cell_line_no, 1));
+ auto pair = project_to_line(v0, v1, p);
+ if (pair.second < min_distance_square)
+ {
+ result = pair.first;
+ min_distance_square = pair.second;
+ }
+ }
+ }
+ }
+
+ Assert(min_distance_square < std::numeric_limits<double>::max(),
+ ExcInternalError());
+
+ // If necessary, slightly adjust the computed point so that it is closer to
+ // being on the surface of the reference cell. Due to roundoff it is difficult
+ // to place points on sloped surfaces (e.g., for Pyramids) so this check isn't
+ // perfect but does improve the accuracy of the projected point.
+ if (!contains_point(result, 0.0))
+ {
+ constexpr unsigned int x_index = 0;
+ constexpr unsigned int y_index = (dim >= 2 ? 1 : 0);
+ constexpr unsigned int z_index = (dim >= 3 ? 2 : 0);
+ switch (this->kind)
+ {
+ case ReferenceCells::Vertex:
+ Assert(false, ExcInternalError());
+ break;
+ // the bounds for each dimension of a hypercube are mutually
+ // independent:
+ case ReferenceCells::Line:
+ case ReferenceCells::Quadrilateral:
+ case ReferenceCells::Hexahedron:
+ for (unsigned int d = 0; d < dim; ++d)
+ result[d] = std_cxx17::clamp(result[d], 0.0, 1.0);
+ // simplices can use the standard definition of a simplex:
+ break;
+ case ReferenceCells::Triangle:
+ result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0);
+ result[y_index] =
+ std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
+ break;
+ case ReferenceCells::Tetrahedron:
+ result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0);
+ result[y_index] =
+ std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
+ result[z_index] =
+ std_cxx17::clamp(result[z_index],
+ 0.0,
+ 1.0 - result[x_index] - result[y_index]);
+ break;
+ // wedges and pyramids are more ad-hoc:
+ case ReferenceCells::Wedge:
+ result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0);
+ result[y_index] =
+ std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
+ result[z_index] = std_cxx17::clamp(result[z_index], 0.0, 1.0);
+ break;
+ case ReferenceCells::Pyramid:
+ {
+ result[x_index] = std_cxx17::clamp(result[x_index], -1.0, 1.0);
+ result[y_index] = std_cxx17::clamp(result[y_index], -1.0, 1.0);
+ // It suffices to transform everything to the first quadrant to
+ // adjust z:
+ const auto x_abs = std::abs(result[x_index]);
+ const auto y_abs = std::abs(result[y_index]);
+
+ if (y_abs <= x_abs)
+ result[z_index] =
+ std_cxx17::clamp(result[z_index], 0.0, 1.0 - x_abs);
+ else
+ result[z_index] =
+ std_cxx17::clamp(result[z_index], 0.0, 1.0 - y_abs);
+ }
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ // We should be within 4 * eps of the cell by this point. The roundoff error
+ // comes from, e.g., computing (1 - x) + x when moving points onto the top of
+ // a Pyramid.
+ Assert(contains_point(result, 4.0 * std::numeric_limits<double>::epsilon()),
+ ExcInternalError());
+
+ return result;
+}
+
+
+
std::ostream &
operator<<(std::ostream &out, const ReferenceCell &reference_cell)
{
--- /dev/null
+#include <deal.II/grid/reference_cell.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const ReferenceCell &reference_cell,
+ const Point<dim> & p,
+ const bool print = true)
+{
+ const Point<dim> projected = reference_cell.closest_point(p);
+ const double distance_square = projected.distance_square(p);
+ // its difficult to place points exactly on sloped surfaces (e.g., for
+ // pyramids), but we should always be close
+ AssertThrow(reference_cell.contains_point(
+ projected, 2.0 * std::numeric_limits<double>::epsilon()),
+ ExcInternalError());
+ if (print)
+ deallog << "point = " << p << std::endl
+ << " d^2 = " << distance_square << std::endl;
+
+ // The distance we get should be at least as good as the distance to the
+ // nearest vertex
+ for (unsigned int vertex_no : reference_cell.vertex_indices())
+ AssertThrow(distance_square <=
+ p.distance_square(
+ reference_cell.template vertex<dim>(vertex_no)),
+ ExcInternalError());
+
+ // Also generate a lot of points in the unit cell - we should be closer than
+ // any generated point
+ const unsigned int n_points_1d = 100;
+ if (dim == 2)
+ {
+ for (unsigned int i = 0; i < n_points_1d; ++i)
+ for (unsigned int j = 0; j < n_points_1d; ++j)
+ {
+ const Point<dim> grid_point(i / double(n_points_1d),
+ j / double(n_points_1d));
+ if (reference_cell.contains_point(grid_point))
+ AssertThrow(distance_square <= grid_point.distance_square(p),
+ ExcInternalError());
+ }
+ }
+ else if (dim == 3)
+ {
+ for (unsigned int i = 0; i < n_points_1d; ++i)
+ for (unsigned int j = 0; j < n_points_1d; ++j)
+ for (unsigned int k = 0; j < n_points_1d; ++j)
+ {
+ // scale things so that we cover the bounding box of both pyramids
+ // and wedges
+ const Point<dim> grid_point(-1.0 + 2.0 * i / double(n_points_1d),
+ -1.0 + 2.0 * j / double(n_points_1d),
+ k / double(n_points_1d));
+ if (reference_cell.contains_point(grid_point))
+ AssertThrow(distance_square <= grid_point.distance_square(p),
+ ExcInternalError());
+ }
+ }
+}
+
+int
+main()
+{
+ initlog();
+
+ for (const auto &reference_cell : {ReferenceCells::Line})
+ {
+ deallog << "reference cell = " << reference_cell.to_string() << std::endl;
+ test(reference_cell, Point<1>(-1));
+ test(reference_cell, Point<1>(-0.5));
+ test(reference_cell, Point<1>(0.5));
+ test(reference_cell, Point<1>(1.0));
+ test(reference_cell, Point<1>(1.25));
+ }
+
+ for (const auto &reference_cell :
+ {ReferenceCells::Triangle, ReferenceCells::Quadrilateral})
+ {
+ deallog << "reference cell = " << reference_cell.to_string() << std::endl;
+ test(reference_cell, Point<2>(-0.5, 0.5));
+ test(reference_cell, Point<2>(-1, -1));
+ test(reference_cell, Point<2>(0.5, -0.5));
+ test(reference_cell, Point<2>(1.0, 0.0));
+ test(reference_cell, Point<2>(1.0, 1.0));
+ test(reference_cell, Point<2>(0.25, 0.25));
+
+ for (unsigned int i = 0; i < 1e3; ++i)
+ test(reference_cell,
+ Point<2>(random_value(-4.0, 4.0), random_value(-4.0, 4.0)),
+ false);
+ }
+
+ for (const auto &reference_cell : {ReferenceCells::Tetrahedron,
+ ReferenceCells::Pyramid,
+ ReferenceCells::Wedge,
+ ReferenceCells::Hexahedron})
+ {
+ deallog << "reference cell = " << reference_cell.to_string() << std::endl;
+ test(reference_cell, Point<3>(-0.5, 0.5, 0.5));
+ test(reference_cell, Point<3>(-1, -1, -1));
+ test(reference_cell, Point<3>(0.5, -0.5, -0.5));
+ test(reference_cell, Point<3>(1.0, 0.0, 0.0));
+ test(reference_cell, Point<3>(1.0, 1.0, 1.0));
+ test(reference_cell, Point<3>(-0.75, 0.5, 1.5));
+ test(reference_cell, Point<3>(2.0, 2.0, 1.0));
+ test(reference_cell, Point<3>(0.25, 0.25, 0.25));
+ test(reference_cell, Point<3>(1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0));
+ test(reference_cell, Point<3>(0.7, 0.9, 0.1));
+ test(reference_cell, Point<3>(0.08, 0.18, -0.07));
+ test(reference_cell, Point<3>(0.0, 0.0, 0.0));
+
+ for (unsigned int i = 0; i < 1e3; ++i)
+ test(reference_cell,
+ Point<3>(random_value(-4.0, 4.0),
+ random_value(-4.0, 4.0),
+ random_value(-4.0, 4.0)),
+ false);
+ }
+ deallog << "OK!" << std::endl;
+}
--- /dev/null
+
+DEAL::reference cell = Line
+DEAL::point = -1.00000
+DEAL:: d^2 = 1.00000
+DEAL::point = -0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = 0.500000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.25000
+DEAL:: d^2 = 0.0625000
+DEAL::reference cell = Tri
+DEAL::point = -0.500000 0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = -1.00000 -1.00000
+DEAL:: d^2 = 2.00000
+DEAL::point = 0.500000 -0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = 1.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000 1.00000
+DEAL:: d^2 = 0.500000
+DEAL::point = 0.250000 0.250000
+DEAL:: d^2 = 0.00000
+DEAL::reference cell = Quad
+DEAL::point = -0.500000 0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = -1.00000 -1.00000
+DEAL:: d^2 = 2.00000
+DEAL::point = 0.500000 -0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = 1.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000 1.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.250000 0.250000
+DEAL:: d^2 = 0.00000
+DEAL::reference cell = Tet
+DEAL::point = -0.500000 0.500000 0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = -1.00000 -1.00000 -1.00000
+DEAL:: d^2 = 3.00000
+DEAL::point = 0.500000 -0.500000 -0.500000
+DEAL:: d^2 = 0.500000
+DEAL::point = 1.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000 1.00000 1.00000
+DEAL:: d^2 = 1.33333
+DEAL::point = -0.750000 0.500000 1.50000
+DEAL:: d^2 = 1.06250
+DEAL::point = 2.00000 2.00000 1.00000
+DEAL:: d^2 = 5.50000
+DEAL::point = 0.250000 0.250000 0.250000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.333333 0.333333 0.333333
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.700000 0.900000 0.100000
+DEAL:: d^2 = 0.190000
+DEAL::point = 0.0800000 0.180000 -0.0700000
+DEAL:: d^2 = 0.00490000
+DEAL::point = 0.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::reference cell = Pyramid
+DEAL::point = -0.500000 0.500000 0.500000
+DEAL:: d^2 = 0.00000
+DEAL::point = -1.00000 -1.00000 -1.00000
+DEAL:: d^2 = 1.00000
+DEAL::point = 0.500000 -0.500000 -0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = 1.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000 1.00000 1.00000
+DEAL:: d^2 = 0.666667
+DEAL::point = -0.750000 0.500000 1.50000
+DEAL:: d^2 = 0.875000
+DEAL::point = 2.00000 2.00000 1.00000
+DEAL:: d^2 = 3.00000
+DEAL::point = 0.250000 0.250000 0.250000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.333333 0.333333 0.333333
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.700000 0.900000 0.100000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.0800000 0.180000 -0.0700000
+DEAL:: d^2 = 0.00490000
+DEAL::point = 0.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::reference cell = Wedge
+DEAL::point = -0.500000 0.500000 0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = -1.00000 -1.00000 -1.00000
+DEAL:: d^2 = 3.00000
+DEAL::point = 0.500000 -0.500000 -0.500000
+DEAL:: d^2 = 0.500000
+DEAL::point = 1.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000 1.00000 1.00000
+DEAL:: d^2 = 0.500000
+DEAL::point = -0.750000 0.500000 1.50000
+DEAL:: d^2 = 0.812500
+DEAL::point = 2.00000 2.00000 1.00000
+DEAL:: d^2 = 4.50000
+DEAL::point = 0.250000 0.250000 0.250000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.333333 0.333333 0.333333
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.700000 0.900000 0.100000
+DEAL:: d^2 = 0.180000
+DEAL::point = 0.0800000 0.180000 -0.0700000
+DEAL:: d^2 = 0.00490000
+DEAL::point = 0.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::reference cell = Hex
+DEAL::point = -0.500000 0.500000 0.500000
+DEAL:: d^2 = 0.250000
+DEAL::point = -1.00000 -1.00000 -1.00000
+DEAL:: d^2 = 3.00000
+DEAL::point = 0.500000 -0.500000 -0.500000
+DEAL:: d^2 = 0.500000
+DEAL::point = 1.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = 1.00000 1.00000 1.00000
+DEAL:: d^2 = 0.00000
+DEAL::point = -0.750000 0.500000 1.50000
+DEAL:: d^2 = 0.812500
+DEAL::point = 2.00000 2.00000 1.00000
+DEAL:: d^2 = 2.00000
+DEAL::point = 0.250000 0.250000 0.250000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.333333 0.333333 0.333333
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.700000 0.900000 0.100000
+DEAL:: d^2 = 0.00000
+DEAL::point = 0.0800000 0.180000 -0.0700000
+DEAL:: d^2 = 0.00490000
+DEAL::point = 0.00000 0.00000 0.00000
+DEAL:: d^2 = 0.00000
+DEAL::OK!
projects / dealii.git / commitdiff