--- /dev/null
+New: Now opencascade also works for spacedim == 2.
+<br>
+(Luca Heltai, 2017/12/08)
* cases they can be used in two dimensions as well, and the third dimension
* will be set to zero.
*
- * @author Luca Heltai, Andrea Mola, 2011--2014.
+ * If you whish to use these tools when the dimension of the space is two, then
+ * make sure your cad files are actually flat and that all z coordinates are equal
+ * to zero, as otherwise you will get many exceptions.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2017.
*/
namespace OpenCASCADE
{
const double tolerance=1e-7);
/**
- * Creates a 3D smooth BSpline curve passing through the points in the
- * assigned vector, and store it in the returned TopoDS_Shape (which is of
- * type TopoDS_Edge). The points are reordered internally according to their
+ * Creates a smooth BSpline curve passing through the points in the assigned
+ * vector, and store it in the returned TopoDS_Shape (which is of type
+ * TopoDS_Edge). The points are reordered internally according to their
* scalar product with the direction, if direction is different from zero,
* otherwise they are used as passed. Notice that this function changes the
* input points if required by the algorithm.
* points. If the algorithm fails in generating such a curve, an exception
* is thrown.
*/
- TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
- const Tensor<1,3> &direction=Tensor<1,3>(),
+ template<int dim>
+ TopoDS_Edge interpolation_curve(std::vector<Point<dim> > &curve_points,
+ const Tensor<1, dim> &direction=Tensor<1,dim>(),
const bool closed=false,
const double tolerance=1e-7);
std::vector<TopoDS_Vertex> &vertices);
/**
- * Create a triangulation from a single face. This class extract the first u
- * and v parameter of the parametric surface making up this face, and
- * creates a Triangulation<2,3> containing a single coarse cell reflecting
+ * Create a triangulation from a single face. This class extracts the first u
+ * and v parameter of the parametric surface making up this face, and creates
+ * a Triangulation<2,spacedim> containing a single coarse cell reflecting
* this face. If the surface is not a trimmed surface, the vertices of this
* cell will coincide with the TopoDS_Vertex vertices of the original
* TopoDS_Face. This, however, is often not the case, and the user should be
* careful on how this mesh is used.
+ *
+ * If you call this function with a Triangulation<2,2>, make sure that the
+ * input face has all z coordinates set to zero, or you'll get an exception.
*/
+ template <int spacedim>
void create_triangulation(const TopoDS_Face &face,
- Triangulation<2,3> &tria);
+ Triangulation<2,spacedim> &tria);
/**
* the u coordinate and the v coordinate (which is different from zero only
* if the resulting shape is a face).
*/
- std::tuple<Point<3>, TopoDS_Shape, double, double>
+ template <int dim>
+ std::tuple<Point<dim>, TopoDS_Shape, double, double>
project_point_and_pull_back(const TopoDS_Shape &in_shape,
- const Point<3> &origin,
+ const Point<dim> &origin,
const double tolerance=1e-7);
/**
* are iterated, faces first, then edges, and the returned point is the
* closest one to the @p in_shape, regardless of its type.
*/
- Point<3> closest_point(const TopoDS_Shape &in_shape,
- const Point<3> &origin,
- const double tolerance=1e-7);
+ template <int dim>
+ Point<dim> closest_point(const TopoDS_Shape &in_shape,
+ const Point<dim> &origin,
+ const double tolerance=1e-7);
/**
* Given an elementary shape @p in_shape and the reference coordinates
* used as input to this function. If this is not the case, an Exception is
* thrown.
*/
- Point<3> push_forward(const TopoDS_Shape &in_shape,
- const double u,
- const double v);
+ template <int dim>
+ Point<dim> push_forward(const TopoDS_Shape &in_shape,
+ const double u,
+ const double v);
/**
*
* The optional @p tolerance parameter is used to compute distances.
*/
- Point<3> line_intersection(const TopoDS_Shape &in_shape,
- const Point<3> &origin,
- const Tensor<1,3> &direction,
- const double tolerance=1e-7);
+ template <int dim>
+ Point<dim> line_intersection(const TopoDS_Shape &in_shape,
+ const Point<dim> &origin,
+ const Tensor<1,dim> &direction,
+ const double tolerance=1e-7);
/**
- * Convert OpenCASCADE point into a Point<3>.
+ * Convert OpenCASCADE point into a Point<spacedim>.
+ *
+ * The tolerance argument is used to check if the non used components of the
+ * OpenCASCADE point are close to zero. If this is not the case, an assertion
+ * is thrown in debug mode.
*/
- Point<3> point(const gp_Pnt &p);
+ template <int spacedim>
+ Point<spacedim> point(const gp_Pnt &p, const double &tolerance=1e-10);
/**
* Convert Point<3> into OpenCASCADE point.
*/
- gp_Pnt point(const Point<3> &p);
+ template <int spacedim>
+ gp_Pnt point(const Point<spacedim> &p);
/**
* optional parameter is used as a relative tolerance when comparing
* objects.
*/
- bool point_compare(const Point<3> &p1,
- const Point<3> &p2,
- const Tensor<1,3> &direction = Tensor<1,3>(),
+ template <int dim>
+ bool point_compare(const Point<dim> &p1,
+ const Point<dim> &p2,
+ const Tensor<1,dim> &direction = Tensor<1,dim>(),
const double tolerance = 1e-10);
* Exception thrown when the point specified as argument does not lie
* between @p tolerance from the given TopoDS_Shape.
*/
+ template <int dim>
DeclException1 (ExcPointNotOnManifold,
- Point<3>,
+ Point<dim>,
<<"The point [ "<<arg1<<" ] is not on the manifold.");
/**
* Exception thrown when the point specified as argument cannot be projected
* to the manifold.
*/
+ template <int dim>
DeclException1 (ExcProjectionFailed,
- Point<3>,
+ Point<dim>,
<<"Projection of point [ "<< arg1
<< " ] failed.");
SET(_inst
boundary_lib.inst.in
+ utilities.inst.in
)
FILE(GLOB _header
Assert(closest_point(sh, surrounding_points[i], tolerance)
.distance(surrounding_points[i]) <
std::max(tolerance*surrounding_points[i].norm(), tolerance),
- ExcPointNotOnManifold(surrounding_points[i]));
+ ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
#endif
return closest_point(sh, candidate,tolerance);
}
Assert(closest_point(sh, surrounding_points[i],tolerance)
.distance(surrounding_points[i]) <
std::max(tolerance*surrounding_points[i].norm(), tolerance),
- ExcPointNotOnManifold(surrounding_points[i]));
+ ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
#endif
return line_intersection(sh, candidate, direction, tolerance);
}
Assert(closest_point(sh, surrounding_points[i], tolerance)
.distance(surrounding_points[i]) <
std::max(tolerance*surrounding_points[i].norm(), tolerance),
- ExcPointNotOnManifold(surrounding_points[i]));
+ ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
}
#endif
curve(curve_adaptor(sh)),
tolerance(tolerance),
length(shape_length(sh))
- {
- Assert(spacedim == 3, ExcNotImplemented());
- }
+ {}
template <int dim, int spacedim>
ShapeAnalysis_Curve curve_analysis;
gp_Pnt proj;
const double dist = curve_analysis.Project(curve->GetCurve(), point(space_point), tolerance, proj, t, true);
- Assert(dist < tolerance*length, ExcPointNotOnManifold(space_point));
+ Assert(dist < tolerance*length, ExcPointNotOnManifold<spacedim>(space_point));
(void)dist; // Silence compiler warning in Release mode.
return Point<1>(GCPnts_AbscissaPoint::Length(curve->GetCurve(),curve->GetCurve().FirstParameter(),t));
}
{
GCPnts_AbscissaPoint AP(curve->GetCurve(), chart_point[0], curve->GetCurve().FirstParameter());
gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
- return point(P);
+ return point<spacedim>(P);
}
template <int dim, int spacedim>
tolerance(tolerance)
{}
- template <>
- Point<2>
- NURBSPatchManifold<2, 3>::
- pull_back(const Point<3> &space_point) const
+ template <int dim, int spacedim> Point<2>
+ NURBSPatchManifold<dim, spacedim>::
+ pull_back(const Point<spacedim> &space_point) const
{
Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
return Point<2>(u,v);
}
- template <>
- Point<3>
- NURBSPatchManifold<2, 3>::
+ template <int dim, int spacedim>
+ Point<spacedim>
+ NURBSPatchManifold<dim, spacedim>::
push_forward(const Point<2> &chart_point) const
{
- return ::dealii::OpenCASCADE::push_forward(face, chart_point[0], chart_point[1]);
+ return ::dealii::OpenCASCADE::push_forward<spacedim>(face, chart_point[0], chart_point[1]);
}
- template <>
- DerivativeForm<1,2,3>
- NURBSPatchManifold<2, 3>::
+ template <int dim, int spacedim>
+ DerivativeForm<1,2,spacedim>
+ NURBSPatchManifold<dim,spacedim>::
push_forward_gradient(const Point<2> &chart_point) const
{
- DerivativeForm<1,2,3> DX;
+ DerivativeForm<1,2,spacedim> DX;
Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
gp_Pnt q;
DX[0][0] = Du.X();
DX[1][0] = Du.Y();
- DX[2][0] = Du.Z();
+ if (spacedim>2)
+ DX[2][0] = Du.Z();
+ else
+ Assert(std::abs(Du.Z()) < tolerance,
+ ExcMessage("Expecting derivative along Z to be zero! Bailing out."));
DX[0][1] = Dv.X();
DX[1][1] = Dv.Y();
- DX[2][1] = Dv.Z();
-
+ if (spacedim>2)
+ DX[2][1] = Dv.Z();
+ else
+ Assert(std::abs(Dv.Z()) < tolerance,
+ ExcMessage("Expecting derivative along Z to be zero! Bailing out."));
return DX;
}
- template <>
+ template <int dim, int spacedim>
std::tuple<double, double, double, double>
- NURBSPatchManifold<2, 3>::
+ NURBSPatchManifold<dim,spacedim>::
get_uv_bounds() const
{
Standard_Real umin, umax, vmin, vmax;
}
// Explicit instantiations
- template class NURBSPatchManifold<2, 3 >;
#include "boundary_lib.inst"
} // end namespace OpenCASCADE
template class DirectionalProjectionBoundary<deal_II_dimension, 3>;
template class NormalToMeshProjectionBoundary<deal_II_dimension, 3>;
template class ArclengthProjectionLineManifold<deal_II_dimension, 3>;
+ template class NURBSPatchManifold<deal_II_dimension, 3>;
+#if deal_II_dimension <= 2
+ template class DirectionalProjectionBoundary<deal_II_dimension, 2>;
+ template class ArclengthProjectionLineManifold<deal_II_dimension, 2>;
+ template class NURBSPatchManifold<deal_II_dimension, 2>;
+#endif
}
}
}
- gp_Pnt point(const Point<3> &p)
+ template <int spacedim>
+ gp_Pnt point(const Point<spacedim> &p)
{
- return gp_Pnt(p(0), p(1), p(2));
+ switch (spacedim)
+ {
+ case 1:
+ return gp_Pnt(p[0], 0, 0);
+ case 2:
+ return gp_Pnt(p[0], p[1], 0);
+ case 3:
+ return gp_Pnt(p[0], p[1], p[2]);
+ }
}
-
- Point<3> point(const gp_Pnt &p)
+ template <int spacedim>
+ Point<spacedim> point(const gp_Pnt &p, const double &tolerance)
{
- return Point<3>(p.X(), p.Y(), p.Z());
+ (void) tolerance;
+ switch (spacedim)
+ {
+ case 1:
+ Assert(std::abs(p.Y()) <= tolerance,
+ ExcMessage("Cannot convert OpenCASCADE point to 1d if p.Y() != 0."));
+ Assert(std::abs(p.Z()) <= tolerance,
+ ExcMessage("Cannot convert OpenCASCADE point to 1d if p.Z() != 0."));
+ return Point<spacedim>(p.X());
+ case 2:
+ Assert(std::abs(p.Z()) <= tolerance,
+ ExcMessage("Cannot convert OpenCASCADE point to 2d if p.Z() != 0."));
+ return Point<spacedim>(p.X(), p.Y());
+ case 3:
+ return Point<spacedim>(p.X(), p.Y(), p.Z());
+ }
}
- bool point_compare(const Point<3> &p1,
- const Point<3> &p2,
- const Tensor<1,3> &direction,
+ template<int dim>
+ bool point_compare(const Point<dim> &p1,
+ const Point<dim> &p2,
+ const Tensor<1,dim> &direction,
const double tolerance)
{
const double rel_tol=std::max(tolerance, std::max(p1.norm(), p2.norm())*tolerance);
if (direction.norm() > 0.0)
return (p1*direction < p2*direction-rel_tol);
else
- for (int d=2; d>=0; --d)
+ for (int d=dim; d>=0; --d)
if (p1[d] < p2[d]-rel_tol)
return true;
else if (p2[d] < p1[d]-rel_tol)
return out_shape;
}
-
- Point<3> line_intersection(const TopoDS_Shape &in_shape,
- const Point<3> &origin,
- const Tensor<1,3> &direction,
- const double tolerance)
+ template <int dim>
+ Point<dim> line_intersection(const TopoDS_Shape &in_shape,
+ const Point<dim> &origin,
+ const Tensor<1,dim> &direction,
+ const double tolerance)
{
// translating original Point<dim> to gp point
gp_Pnt P0 = point(origin);
- gp_Ax1 gpaxis(P0, gp_Dir(direction[0], direction[1], direction[2]));
+ gp_Ax1 gpaxis(P0, gp_Dir(direction[0], dim > 1 ? direction[1] : 0, dim>2 ? direction[2] : 0));
gp_Lin line(gpaxis);
// destination point
Assert(Inters.IsDone(), ExcMessage("Could not project point."));
double minDistance = 1e7;
- Point<3> result;
+ Point<dim> result;
for (int i=0; i<Inters.NbPnt(); ++i)
{
const double distance = point(origin).Distance(Inters.Pnt(i+1));
if (distance < minDistance)
{
minDistance = distance;
- result = point(Inters.Pnt(i+1));
+ result = point<dim>(Inters.Pnt(i+1));
}
}
return result;
}
- TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
- const Tensor<1,3> &direction,
+ template <int dim>
+ TopoDS_Edge interpolation_curve(std::vector<Point<dim> > &curve_points,
+ const Tensor<1,dim> &direction,
const bool closed,
const double tolerance)
{
if (direction*direction > 0)
{
std::sort(curve_points.begin(), curve_points.end(),
- [&](const Point<3> &p1, const Point<3> &p2)
+ [&](const Point<dim> &p1, const Point<dim> &p2)
{
return OpenCASCADE::point_compare(p1, p2, direction, tolerance);
});
unsigned int point_index = start_point_index;
// point_index and face_index always run together
- std::vector<Point<3> > pointlist;
+ std::vector<Point<spacedim> > pointlist;
do
{
visited_faces[face_index] = true;
auto current_point = vert_to_point[point_index];
- if (spacedim==2)
- pointlist.push_back(Point<3>(current_point[0],current_point[1],0));
- else
- pointlist.push_back(Point<3>(current_point[0],current_point[1],current_point[2]));
+ pointlist.push_back(current_point);
// Get next point
if (face_to_verts[face_index].first != point_index )
}
while (point_index != start_point_index);
- interpolation_curves.push_back(interpolation_curve(pointlist, Tensor<1,3>(), true));
+ interpolation_curves.push_back(interpolation_curve(pointlist, Tensor<1,spacedim>(), true));
finished = true;
for (const auto &f : visited_faces)
}
-
- std::tuple<Point<3>, TopoDS_Shape, double, double>
- project_point_and_pull_back(const TopoDS_Shape &in_shape,
- const Point<3> &origin,
- const double tolerance)
+ template<int dim>
+ std::tuple<Point<dim>, TopoDS_Shape, double, double> project_point_and_pull_back(const TopoDS_Shape &in_shape,
+ const Point<dim> &origin,
+ const double tolerance)
{
TopExp_Explorer exp;
gp_Pnt Pproj = point(origin);
SurfToProj->D0(proj_params.X(),proj_params.Y(),tmp_proj);
- double distance = point(tmp_proj).distance(origin);
+ double distance = point<dim>(tmp_proj).distance(origin);
if (distance < minDistance)
{
minDistance = distance;
}
Assert(counter > 0, ExcMessage("Could not find projection points."));
- return std::tuple<Point<3>, TopoDS_Shape, double, double>
- (point(Pproj),out_shape, u, v);
+ return std::tuple<Point<dim>, TopoDS_Shape, double, double>
+ (point<dim>(Pproj),out_shape, u, v);
}
- Point<3> closest_point(const TopoDS_Shape &in_shape,
- const Point<3> &origin,
- const double tolerance)
+ template<int dim>
+ Point<dim> closest_point(const TopoDS_Shape &in_shape,
+ const Point<dim> &origin,
+ const double tolerance)
{
- std::tuple<Point<3>, TopoDS_Shape, double, double>
+ std::tuple<Point<dim>, TopoDS_Shape, double, double>
ref = project_point_and_pull_back(in_shape, origin, tolerance);
return std::get<0>(ref);
}
return push_forward_and_differential_forms(face, u, v, tolerance);
}
- Point<3> push_forward(const TopoDS_Shape &in_shape,
- const double u,
- const double v)
+ template<int dim>
+ Point<dim> push_forward(const TopoDS_Shape &in_shape,
+ const double u,
+ const double v)
{
switch (in_shape.ShapeType())
{
case TopAbs_FACE:
{
BRepAdaptor_Surface surf(TopoDS::Face(in_shape));
- return point(surf.Value(u,v));
+ return point<dim>(surf.Value(u,v));
}
case TopAbs_EDGE:
{
BRepAdaptor_Curve curve(TopoDS::Edge(in_shape));
- return point(curve.Value(u));
+ return point<dim>(curve.Value(u));
}
default:
Assert(false, ExcUnsupportedShape());
}
- return Point<3>();
+ return Point<dim>();
}
std::tuple<Point<3>, Tensor<1,3>, double, double>
Max_Curvature *= -1;
}
- return std::tuple<Point<3>, Tensor<1,3>, double, double>(point(Value), normal, Min_Curvature, Max_Curvature);
+ return std::tuple<Point<3>, Tensor<1,3>, double, double>(point<3>(Value), normal, Min_Curvature, Max_Curvature);
}
+ template<int spacedim>
void create_triangulation(const TopoDS_Face &face,
- Triangulation<2,3> &tria)
+ Triangulation<2,spacedim> &tria)
{
BRepAdaptor_Surface surf(face);
const double u0 = surf.FirstUParameter();
const double v1 = surf.LastVParameter();
std::vector<CellData<2> > cells;
- std::vector<Point<3> > vertices;
+ std::vector<Point<spacedim> > vertices;
SubCellData t;
- vertices.push_back(point(surf.Value(u0,v0)));
- vertices.push_back(point(surf.Value(u1,v0)));
- vertices.push_back(point(surf.Value(u0,v1)));
- vertices.push_back(point(surf.Value(u1,v1)));
+ vertices.push_back(point<spacedim>(surf.Value(u0,v0)));
+ vertices.push_back(point<spacedim>(surf.Value(u1,v0)));
+ vertices.push_back(point<spacedim>(surf.Value(u0,v1)));
+ vertices.push_back(point<spacedim>(surf.Value(u1,v1)));
CellData<2> cell;
for (unsigned int i=0; i<4; ++i)
tria.create_triangulation(vertices, cells, t);
}
-
-
- // Explicit instantiations
- template
- std::vector<TopoDS_Edge> create_curves_from_triangulation_boundary
- (const Triangulation<2, 2> &triangulation,
- const Mapping<2, 2> &mapping);
-
-
-
- template
- std::vector<TopoDS_Edge> create_curves_from_triangulation_boundary
- (const Triangulation<2, 3> &triangulation,
- const Mapping<2, 3> &mapping);
-
-
+#include "utilities.inst"
} // end namespace
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : DIMENSIONS)
+{
+#if deal_II_dimension > 1
+ // Explicit instantiations for dim = 2 and 3
+ template
+ std::vector<TopoDS_Edge> create_curves_from_triangulation_boundary
+ (const Triangulation<2, deal_II_dimension> &triangulation,
+ const Mapping<2, deal_II_dimension> &mapping);
+
+ template bool point_compare(const Point<deal_II_dimension>& p1,
+ const Point<deal_II_dimension>& p2,
+ const Tensor<1,deal_II_dimension>& direction,
+ const double tolerance );
+
+ template Point<deal_II_dimension> point(const gp_Pnt &p, const double &tolerance);
+
+ template gp_Pnt point(const Point<deal_II_dimension> &p);
+
+ template
+ TopoDS_Edge interpolation_curve(std::vector<Point<deal_II_dimension> >& curve_points,
+ const Tensor<1, deal_II_dimension>& direction=Tensor<1,deal_II_dimension>(),
+ const bool closed,
+ const double tolerance);
+
+ template Point<deal_II_dimension> push_forward(const TopoDS_Shape &in_shape,
+ const double u,
+ const double v);
+
+ template Point<deal_II_dimension> line_intersection(const TopoDS_Shape &in_shape,
+ const Point<deal_II_dimension> &origin,
+ const Tensor<1,deal_II_dimension> &direction,
+ const double tolerance);
+
+ template void create_triangulation(const TopoDS_Face &face,
+ Triangulation<2,deal_II_dimension> &tria);
+
+ template std::tuple<Point<deal_II_dimension>, TopoDS_Shape, double, double>
+ project_point_and_pull_back(const TopoDS_Shape &in_shape,
+ const Point<deal_II_dimension>& origin,
+ const double tolerance);
+
+ template Point<deal_II_dimension> closest_point(const TopoDS_Shape &in_shape,
+ const Point<deal_II_dimension> &origin,
+ const double tolerance);
+
+#endif
+}
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 - 2015 by the deal.II authors
+//
+// This file is subject to LGPL 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 "../tests.h"
+
+#include <deal.II/opencascade/utilities.h>
+#include <deal.II/opencascade/boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <Standard_Stream.hxx>
+
+// Create a Triangulation, interpolate its boundary points to a smooth
+// BSpline, and use that as an arlength Boundary Descriptor.
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ std::ofstream logfile("output");
+
+ // Create a bspline passign through the points
+ std::vector<Point<2> > pts;
+ pts.push_back(Point<2>(0,0));
+ pts.push_back(Point<2>(0,1));
+ pts.push_back(Point<2>(1,1));
+ pts.push_back(Point<2>(1,0));
+
+ TopoDS_Edge edge = interpolation_curve(pts, Tensor<1,2>(), true);
+ ArclengthProjectionLineManifold<2,2> manifold(edge);
+
+ Triangulation<2> tria;
+ GridGenerator::hyper_cube(tria);
+
+ tria.set_all_manifold_ids_on_boundary(1);
+ tria.begin_active()->face(2)->set_manifold_id(numbers::invalid_manifold_id);
+ tria.set_manifold(1, manifold);
+
+ tria.refine_global(3);
+ GridOut gridout;
+ gridout.write_msh(tria, logfile);
+
+ return 0;
+}
+
--- /dev/null
+$NOD
+81
+1 0 0 0
+2 1 0 0
+3 0 1 0
+4 1 1 0
+5 0.5 0 0
+6 -0.251337 0.515857 0
+7 1.16406 0.505545 0
+8 0.504074 1.14672 0
+9 0.4896 0.521016 0
+10 0.25 0 0
+11 0.75 0 0
+12 -0.158882 0.242372 0
+13 -0.20277 0.796991 0
+14 1.11495 0.242077 0
+15 1.13535 0.7715 0
+16 0.242844 1.10455 0
+17 0.767092 1.12113 0
+18 0.4948 0.260508 0
+19 0.496837 0.833869 0
+20 0.119131 0.518436 0
+21 0.826831 0.51328 0
+22 0.180414 0.257274 0
+23 0.79253 0.255303 0
+24 0.174797 0.80468 0
+25 0.79798 0.801633 0
+26 0.125 0 0
+27 0.375 0 0
+28 0.625 0 0
+29 0.875 0 0
+30 -0.08483 0.11767 0
+31 -0.217128 0.375124 0
+32 -0.249712 0.660431 0
+33 -0.114709 0.911629 0
+34 1.06551 0.117217 0
+35 1.1481 0.372213 0
+36 1.16089 0.639761 0
+37 1.08298 0.894921 0
+38 0.117841 1.06055 0
+39 0.372227 1.13328 0
+40 0.636548 1.14359 0
+41 0.891132 1.07494 0
+42 0.4974 0.130254 0
+43 0.4922 0.390762 0
+44 0.493218 0.677443 0
+45 0.500456 0.990296 0
+46 -0.0661032 0.517147 0
+47 0.304365 0.519726 0
+48 0.658215 0.517148 0
+49 0.995446 0.509412 0
+50 0.215207 0.128637 0
+51 0.149773 0.387855 0
+52 0.0107661 0.249823 0
+53 0.337607 0.258891 0
+54 0.771265 0.127652 0
+55 0.80968 0.384292 0
+56 0.643665 0.257905 0
+57 0.953738 0.24869 0
+58 0.146964 0.661558 0
+59 0.208821 0.954615 0
+60 -0.0139862 0.800836 0
+61 0.335817 0.819275 0
+62 0.812405 0.657456 0
+63 0.782536 0.961381 0
+64 0.647409 0.817751 0
+65 0.966663 0.786567 0
+66 0.0672094 0.124472 0
+67 0.356303 0.129445 0
+68 -0.033678 0.381489 0
+69 0.320986 0.389308 0
+70 0.634332 0.128953 0
+71 0.915373 0.123867 0
+72 0.65094 0.387527 0
+73 0.978892 0.378252 0
+74 -0.0513737 0.660995 0
+75 0.320091 0.6695 0
+76 0.0516048 0.929231 0
+77 0.354022 0.976278 0
+78 0.652812 0.667449 0
+79 0.986646 0.648609 0
+80 0.641978 0.980671 0
+81 0.927966 0.926509 0
+$ENDNOD
+$ELM
+64
+1 3 0 0 4 1 26 66 30
+2 3 0 0 4 26 10 50 66
+3 3 0 0 4 30 66 52 12
+4 3 0 0 4 66 50 22 52
+5 3 0 0 4 10 27 67 50
+6 3 0 0 4 27 5 42 67
+7 3 0 0 4 50 67 53 22
+8 3 0 0 4 67 42 18 53
+9 3 0 0 4 12 52 68 31
+10 3 0 0 4 52 22 51 68
+11 3 0 0 4 31 68 46 6
+12 3 0 0 4 68 51 20 46
+13 3 0 0 4 22 53 69 51
+14 3 0 0 4 53 18 43 69
+15 3 0 0 4 51 69 47 20
+16 3 0 0 4 69 43 9 47
+17 3 0 0 4 5 28 70 42
+18 3 0 0 4 28 11 54 70
+19 3 0 0 4 42 70 56 18
+20 3 0 0 4 70 54 23 56
+21 3 0 0 4 11 29 71 54
+22 3 0 0 4 29 2 34 71
+23 3 0 0 4 54 71 57 23
+24 3 0 0 4 71 34 14 57
+25 3 0 0 4 18 56 72 43
+26 3 0 0 4 56 23 55 72
+27 3 0 0 4 43 72 48 9
+28 3 0 0 4 72 55 21 48
+29 3 0 0 4 23 57 73 55
+30 3 0 0 4 57 14 35 73
+31 3 0 0 4 55 73 49 21
+32 3 0 0 4 73 35 7 49
+33 3 0 0 4 6 46 74 32
+34 3 0 0 4 46 20 58 74
+35 3 0 0 4 32 74 60 13
+36 3 0 0 4 74 58 24 60
+37 3 0 0 4 20 47 75 58
+38 3 0 0 4 47 9 44 75
+39 3 0 0 4 58 75 61 24
+40 3 0 0 4 75 44 19 61
+41 3 0 0 4 13 60 76 33
+42 3 0 0 4 60 24 59 76
+43 3 0 0 4 33 76 38 3
+44 3 0 0 4 76 59 16 38
+45 3 0 0 4 24 61 77 59
+46 3 0 0 4 61 19 45 77
+47 3 0 0 4 59 77 39 16
+48 3 0 0 4 77 45 8 39
+49 3 0 0 4 9 48 78 44
+50 3 0 0 4 48 21 62 78
+51 3 0 0 4 44 78 64 19
+52 3 0 0 4 78 62 25 64
+53 3 0 0 4 21 49 79 62
+54 3 0 0 4 49 7 36 79
+55 3 0 0 4 62 79 65 25
+56 3 0 0 4 79 36 15 65
+57 3 0 0 4 19 64 80 45
+58 3 0 0 4 64 25 63 80
+59 3 0 0 4 45 80 40 8
+60 3 0 0 4 80 63 17 40
+61 3 0 0 4 25 65 81 63
+62 3 0 0 4 65 15 37 81
+63 3 0 0 4 63 81 41 17
+64 3 0 0 4 81 37 4 41
+$ENDELM