Skip to content

GitLab

Commit 8c42b86a authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

added torus parametrization

parent 6eb075b8
Hide whitespace changes
Inline Side-by-side
dune/curvedsurfacegrid/gridfunction/spheregridfunction.hh
View file @ 8c42b86a
......@@ -43,6 +43,26 @@ namespace Dune
return out;
};
}
//! normal vector = x/|x|
template< class Domain >
Domain normal (const Domain& x) const
{
return x / x.two_norm();
}
//! mean curvature of the sphere = 1/R
template< class Domain >
T mean_curvature (const Domain& /*x*/) const
{
return T(1)/radius_;
}
//! surface area of the sphere = 4*pi*r^2
T area () const
{
return 4*M_PI*radius_*radius_;
}
};
//! construct a grid function representing a sphere parametrization
......
dune/curvedsurfacegrid/gridfunction/torusgridfunction.hh
View file @ 8c42b86a
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
#ifndef DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH
#include <type_traits>
......@@ -49,9 +49,9 @@ namespace Dune
T x0 = x[0]*x[0], x1 = x[1]*x[1], x2 = x0+x1, x3 = sqrt(x2), x4 = T(1)/x3, x5 = r1*x4;
T x6 = sqrt(x2*x2*x2), x7 = r1/x6, x8 = x0*x7, x9 = T(1)-x5, x10 = x8+x9;
T x11 = x[2]*x[2], x12 = x9*x9, x13 = x0*x12, x14 = x1*x12, x15 = x11+x13+x14;
T x16 = r2/sqrt(x15), x17 = x1*x7, x18 = sqrt(power(r2*(x10 + x17)/x15, 3));
T x16 = r2/sqrt(x15), x17 = x1*x7, x18 = r2*(x10 + x17)/sqrt(power(x15,3));
T x19 = sqrt(power(x2,5)), x20 = r1*x19, x21 = x2*x2, x22 = r1*x3;
T x23 = x21-r1*x6; x24 = r1-x3, x25 = x24*x24, x26 = x0*x25+x1*x25+x11*x2;
T x23 = x21-r1*x6, x24 = r1-x3, x25 = x24*x24, x26 = x0*x25+x1*x25+x11*x2;
T x27 = x26/x2, x28 = sqrt(x27), x29 = r2*x25*x28, x30 = x2*x2*x2*x29, x31 = x26*x26;
T x32 = sqrt(x27*x27*x27), x33 = -r1*r2*sqrt(power(x2,13))*x32 + r1*sqrt(power(x2,9))*x31;
T x34 = T(1)/x31, x35 = x[0]*x[1]*x34/power(x2,6), x36 = r2*x[2], x37 = x[0]*x36;
......@@ -65,6 +65,30 @@ namespace Dune
};
};
}
// phi = (sqrt(x**2 + y**2) - r2)**2 + z**2 - r1**2
// N = grad(phi) / |grad(phi)|
template< class Domain >
Domain normal (const Domain& x)
{
using std::sqrt;
T x0 = x[0]*x[0], x1 = x[1]*x[1], x2 = x0+x1, x3 = sqrt(x2), x4 = radius2_-x3;
T x5 = x4*x4, x6 = sqrt((x0*x5 + x1*x5 + x2*x[2]*x[2])/x2), x7 = x4/(x6*x3);
return {-x[0]*x7, -x[1]*x7, x[2]/x6};
}
template< class Domain >
auto mean_curvature (const Domain& x)
{
return 0.0;
}
//! surface area of the torus = 4*pi^2*r1*r2
T area () const
{
return 4*M_PI*M_PI*radius1_*radius2_;
}
};
//! construct a grid function representing a torus parametrization
......@@ -76,4 +100,4 @@ namespace Dune
} // end namespace Dune
#endif // DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
#endif // DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH
dune/curvedsurfacegrid/test/convergence.cc
View file @ 8c42b86a
......@@ -7,29 +7,30 @@
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/curvedsurfacegrid/curvedsurfacegrid.hh>
#include <dune/curvedsurfacegrid/gridfunction/analyticgridfunction.hh>
#include <dune/curvedsurfacegrid/gridfunction/spheregridfunction.hh>
#include <dune/curvedsurfacegrid/gridfunction/torusgridfunction.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/localfunctions/lagrange/pk.hh>
//#include <dune/grid/albertagrid.hh>
#include <dune/foamgrid/foamgrid.hh>
#include <dune/functions/common/differentiablefunctionfromcallables.hh>
//#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
#define STR(s) STR_HELPER(s)
#define STR_HELPER(s) #s
#define TORUS
const int order = 4;
const int quad_order = order+10;
const int num_levels = 5;
// Test Hausdorf-distance
template <class Grid, class Projection>
typename Grid::ctype inf_error(Grid const& grid, Projection const& projection)
template <class Grid>
typename Grid::ctype inf_error(Grid const& grid)
{
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
typename Grid::ctype dist = 0;
auto projection = grid.gridFunction().f();
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
......@@ -46,11 +47,12 @@ typename Grid::ctype inf_error(Grid const& grid, Projection const& projection)
}
// Test integrated Hausdorf-distance
template <class Grid, class Projection>
typename Grid::ctype L2_error(Grid const& grid, Projection const& projection)
template <class Grid>
typename Grid::ctype L2_error(Grid const& grid)
{
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
typename Grid::ctype dist = 0;
auto projection = grid.gridFunction().f();
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
......@@ -67,11 +69,13 @@ typename Grid::ctype L2_error(Grid const& grid, Projection const& projection)
}
// Test integrated error in normal vectors
template <class Grid, class Projection>
typename Grid::ctype normal_error(Grid const& grid, Projection const& projection)
template <class Grid>
typename Grid::ctype normal_error(Grid const& grid)
{
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
typename Grid::ctype dist = 0;
auto projection = grid.gridFunction().f();
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
......@@ -90,8 +94,8 @@ typename Grid::ctype normal_error(Grid const& grid, Projection const& projection
}
// Test integrated error in mean curvature
template <class Grid, class Projection>
typename Grid::ctype curvature_error(Grid const& grid, Projection const& projection)
template <class Grid>
typename Grid::ctype curvature_error(Grid const& grid)
{
using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
......@@ -103,6 +107,8 @@ typename Grid::ctype curvature_error(Grid const& grid, Projection const& project
std::vector<GlobalCoordinate> gradients;
std::vector<GlobalCoordinate> normals;
auto projection = grid.gridFunction().f();
typename Grid::ctype dist = 0;
for (auto const& element : elements(grid.leafGridView()))
{
......@@ -135,7 +141,8 @@ typename Grid::ctype curvature_error(Grid const& grid, Projection const& project
for (std::size_t i = 0; i < gradients.size(); ++i)
H += normals[i].dot(gradients[i]);
dist += std::abs(H - 2) * qp.weight() * geometry.integrationElement(qp.position());
auto H_sphere = projection.mean_curvature(geometry.global(qp.position()));
dist += std::abs(H/2 - H_sphere) * qp.weight() * geometry.integrationElement(qp.position());
}
}
......@@ -161,38 +168,27 @@ int main(int argc, char** argv)
// using HostGrid = AlbertaGrid<2,3>;
using HostGrid = FoamGrid<2,3>;
#if defined(SPHERE)
std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( STR(DUNE_GRID_PATH) "sphere_coarse.msh");
using Signature = FieldVector<double,3>(FieldVector<double,3>);
auto sphere = Functions::makeDifferentiableFunctionFromCallables(
Functions::SignatureTag<Signature>{},
[](auto const& x) {
return x / x.two_norm();
},
[](auto const& x) {
double nrm = x.two_norm();
FieldMatrix<double,3,3> out;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
out[i][j] = ((i == j ? 1 : 0) - (x[i]/nrm) * (x[j]/nrm)) / nrm;
return out;
}
);
//auto sphereGridFct = Functions::makeAnalyticGridViewFunction(sphere, hostGrid->leafGridView());
auto sphereGridFct = analyticGridFunction<HostGrid>(sphere);
using Grid = CurvedSurfaceGrid<decltype(sphereGridFct)>;
auto sphereGridFct = sphereGridFunction<HostGrid>(1.0);
using Grid = CurvedSurfaceGrid<decltype(sphereGridFct),4>;
Grid grid(*hostGrid, sphereGridFct);
#elif defined(TORUS)
std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( STR(DUNE_GRID_PATH) "torus.msh");
auto torusGridFct = torusGridFunction<HostGrid>(2.0, 1.0);
using Grid = CurvedSurfaceGrid<decltype(torusGridFct),4>;
Grid grid(*hostGrid, torusGridFct);
#endif
std::vector<typename Grid::ctype> inf_errors, L2_errors, normal_errors, curvature_errors, edge_lengths;
for (int i = 0; i < num_levels; ++i) {
if (i > 0)
grid.globalRefine(1);
inf_errors.push_back( inf_error(grid, sphere) );
L2_errors.push_back( L2_error(grid, sphere) );
normal_errors.push_back( normal_error(grid, sphere) );
curvature_errors.push_back( curvature_error(grid, sphere) );
inf_errors.push_back( inf_error(grid) );
L2_errors.push_back( L2_error(grid) );
normal_errors.push_back( normal_error(grid) );
curvature_errors.push_back( curvature_error(grid) );
edge_lengths.push_back( edge_length(grid) );
}
......
example/torus.geo 0 → 100644
View file @ 8c42b86a
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {2, 0, 0, 1.0};
//+
Point(3) = {-2, 0, 0, 1.0};
//+
Point(4) = {0, 2, 0, 1.0};
//+
Point(5) = {0, -2, 0, 1.0};
//+
Point(6) = {1, 0, 0, 1.0};
//+
Point(7) = {3, 0, 0, 1.0};
//+
Point(8) = {-1, 0, 0, 1.0};
//+
Point(9) = {-3, 0, 0, 1.0};
//+
Point(10) = {0, 1, 0, 1.0};
//+
Point(11) = {0, 3, 0, 1.0};
//+
Point(12) = {0, -1, 0, 1.0};
//+
Point(13) = {0, -3, 0, 1.0};
//+
Point(14) = {2, 0, 1, 1.0};
//+
Point(15) = {2, 0, -1, 1.0};
//+
Point(16) = {-2, 0, 1, 1.0};
//+
Point(17) = {-2, 0, -1, 1.0};
//+
Point(18) = {0, 2, 1, 1.0};
//+
Point(19) = {0, 2, -1, 1.0};
//+
Point(20) = {0, -2, 1, 1.0};
//+
Point(21) = {0, -2, -1, 1.0};
//+
Circle(1) = {7, 2, 14};
//+
Circle(2) = {14, 2, 6};
//+
Circle(3) = {6, 2, 15};
//+
Circle(4) = {15, 2, 7};
//+
Circle(5) = {9, 3, 16};
//+
Circle(6) = {16, 3, 8};
//+
Circle(7) = {8, 3, 17};
//+
Circle(8) = {17, 3, 9};
//+
Circle(9) = {11, 4, 18};
//+
Circle(10) = {18, 4, 10};
//+
Circle(11) = {10, 4, 19};
//+
Circle(12) = {19, 4, 11};
//+
Circle(13) = {13, 5, 20};
//+
Circle(14) = {20, 5, 12};
//+
Circle(15) = {12, 5, 21};
//+
Circle(16) = {21, 5, 13};
//+
Circle(17) = {12, 1, 6};
//+
Circle(18) = {6, 1, 10};
//+
Circle(19) = {10, 1, 8};
//+
Circle(20) = {8, 1, 12};
//+
Circle(21) = {7, 1, 11};
//+
Circle(22) = {11, 1, 9};
//+
Circle(23) = {9, 1, 13};
//+
Circle(24) = {13, 1, 7};
//+
Point(22) = {0, 0, 1, 1.0};
//+
Point(23) = {0, 0, -1, 1.0};
//+
Circle(25) = {14, 22, 18};
//+
Circle(26) = {18, 22, 16};
//+
Circle(27) = {16, 22, 20};
//+
Circle(28) = {20, 22, 14};
//+
Circle(29) = {15, 23, 19};
//+
Circle(30) = {19, 23, 17};
//+
Circle(31) = {17, 23, 21};
//+
Circle(32) = {21, 23, 15};
//+
Line Loop(1) = {25, -9, -21, 1};
//+
Surface(1) = {1};
//+
Line Loop(2) = {26, -5, -22, 9};
//+
Surface(2) = {2};
//+
Line Loop(3) = {27, -13, -23, 5};
//+
Surface(3) = {3};
//+
Line Loop(4) = {28, -1, -24, 13};
//+
Surface(4) = {4};
//+
Line Loop(5) = {21, -12, -29, 4};
//+
Surface(5) = {5};
//+
Line Loop(6) = {4, -24, -16, 32};
//+
Surface(6) = {6};
//+
Line Loop(7) = {16, -23, -8, 31};
//+
Surface(7) = {7};
//+
Line Loop(8) = {8, -22, -12, 30};
//+
Surface(8) = {8};
//+
Line Loop(9) = {30, -7, -19, 11};
//+
Surface(9) = {9};
//+
Line Loop(10) = {11, -29, -3, 18};
//+
Surface(10) = {10};
//+
Line Loop(11) = {3, -32, -15, 17};
//+
Surface(11) = {11};
//+
Line Loop(12) = {15, -31, -7, 20};
//+
Surface(12) = {12};
//+
Line Loop(13) = {20, -14, -27, 6};
//+
Surface(13) = {13};
//+
Line Loop(14) = {14, 17, -2, -28};
//+
Surface(14) = {14};
//+
Line Loop(15) = {2, 18, -10, -25};
//+
Surface(15) = {15};
//+
Line Loop(16) = {10, 19, -6, -26};
//+
Surface(16) = {16};
example/torus.msh 0 → 100644
View file @ 8c42b86a
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
155
1 0 0 0
2 2 0 0
3 -2 0 0
4 0 2 0
5 0 -2 0
6 1 0 0
7 3 0 0
8 -1 0 0
9 -3 0 0
10 0 1 0
11 0 3 0
12 0 -1 0
13 0 -3 0
14 2 0 1
15 2 0 -1
16 -2 0 1
17 -2 0 -1
18 0 2 1
19 0 2 -1
20 0 -2 1
21 0 -2 -1
22 0 0 1
23 0 0 -1
24 2.707106781187844 0 0.7071067811852514
25 1.292893218814767 0 0.7071067811878623
26 1.292893218812007 0 -0.7071067811851022
27 2.707106781185225 0 -0.7071067811878706
28 -2.707106781187844 0 0.7071067811852514
29 -1.292893218814767 0 0.7071067811878623
30 -1.292893218812007 0 -0.7071067811851022
31 -2.707106781185225 0 -0.7071067811878706
32 0 2.707106781187844 0.7071067811852514
33 0 1.292893218814767 0.7071067811878623
34 0 1.292893218812007 -0.7071067811851022
35 0 2.707106781185225 -0.7071067811878706
36 0 -2.707106781187844 0.7071067811852514
37 0 -1.292893218814767 0.7071067811878623
38 0 -1.292893218812007 -0.7071067811851022
39 0 -2.707106781185225 -0.7071067811878706
40 0.707106781185206 -0.7071067811878892 0
41 0.7071067811878892 0.707106781185206 0
42 -0.707106781185206 0.7071067811878892 0
43 -0.7071067811878892 -0.707106781185206 0
44 2.853169548886248 0.9270509831224205 0
45 2.427050983127515 1.763355756873741 0
46 1.763355756881216 2.427050983122084 0
47 0.9270509831266447 2.853169548884875 0
48 -0.9270509831224205 2.853169548886248 0
49 -1.763355756873741 2.427050983127515 0
50 -2.427050983122084 1.763355756881216 0
51 -2.853169548884875 0.9270509831266447 0
52 -2.853169548886248 -0.9270509831224205 0
53 -2.427050983127515 -1.763355756873741 0
54 -1.763355756881216 -2.427050983122084 0
55 -0.9270509831266447 -2.853169548884875 0
56 0.9270509831224205 -2.853169548886248 0
57 1.763355756873741 -2.427050983127515 0
58 2.427050983122084 -1.763355756881216 0
59 2.853169548884875 -0.9270509831266447 0
60 1.84775906502338 0.765366864728233 1
61 1.414213562375778 1.414213562370412 1
62 0.765366864731742 1.847759065021926 1
63 -0.765366864728233 1.84775906502338 1
64 -1.414213562370412 1.414213562375778 1
65 -1.847759065021926 0.765366864731742 1
66 -1.84775906502338 -0.765366864728233 1
67 -1.414213562375778 -1.414213562370412 1
68 -0.765366864731742 -1.847759065021926 1
69 0.765366864728233 -1.84775906502338 1
70 1.414213562370412 -1.414213562375778 1
71 1.847759065021926 -0.765366864731742 1
72 1.84775906502338 0.765366864728233 -1
73 1.414213562375778 1.414213562370412 -1
74 0.765366864731742 1.847759065021926 -1
75 -0.765366864728233 1.84775906502338 -1
76 -1.414213562370412 1.414213562375778 -1
77 -1.847759065021926 0.765366864731742 -1
78 -1.84775906502338 -0.765366864728233 -1
79 -1.414213562375778 -1.414213562370412 -1
80 -0.765366864731742 -1.847759065021926 -1
81 0.765366864728233 -1.84775906502338 -1
82 1.414213562370412 -1.414213562375778 -1
83 1.847759065021926 -0.765366864731742 -1
84 2.256761916668028 1.393519944959817 0.7071067811865475
85 1.334041972098247 2.294044207478844 0.7071067811865475
86 1.93263800118319 1.953283852568877 0.5877852522924731
87 2.481399437256458 0.6301002387975739 0.809016994374732
88 0.6172986465386927 2.485134578607818 0.8090169943747321
89 -1.393519944959817 2.256761916668028 0.7071067811865475
90 -2.294044207478844 1.334041972098247 0.7071067811865475
91 -1.953283852568877 1.93263800118319 0.5877852522924731
92 -0.6301002387975739 2.481399437256458 0.809016994374732
93 -2.485134578607818 0.6172986465386927 0.8090169943747321
94 -2.256761916668028 -1.393519944959817 0.7071067811865475
95 -1.334041972098247 -2.294044207478844 0.7071067811865475
96 -1.93263800118319 -1.953283852568877 0.5877852522924731
97 -2.481399437256458 -0.6301002387975739 0.809016994374732
98 -0.6172986465386927 -2.485134578607818 0.8090169943747321
99 1.393519944959817 -2.256761916668028 0.7071067811865475
100 2.294044207478844 -1.334041972098247 0.7071067811865475
101 1.953283852568877 -1.93263800118319 0.5877852522924731
102 0.6301002387975739 -2.481399437256458 0.809016994374732
103 2.485134578607818 -0.6172986465386927 0.8090169943747321
104 2.256761916668028 1.393519944959817 -0.7071067811865477
105 1.334041972098247 2.294044207478844 -0.7071067811865475
106 1.93263800118319 1.953283852568877 -0.5877852522924731
107 2.481399437255894 0.6301002387974353 -0.8090169943751674
108 0.6172986465385571 2.485134578607252 -0.8090169943751674
109 1.393519944959817 -2.256761916668028 -0.7071067811865475
110 2.294044207478843 -1.334041972098247 -0.7071067811865477
111 1.953283852568877 -1.932638001183191 -0.587785252292473
112 0.6301002387974354 -2.481399437255894 -0.8090169943751676
113 2.485134578607252 -0.617298646538557 -0.8090169943751676
114 -2.256761916668028 -1.393519944959817 -0.7071067811865475
115 -1.334041972098247 -2.294044207478843 -0.7071067811865477
116 -1.932638001183191 -1.953283852568877 -0.587785252292473
117 -2.481399437255894 -0.6301002387974354 -0.8090169943751676
118 -0.617298646538557 -2.485134578607252 -0.8090169943751676
119 -1.393519944959817 2.256761916668028 -0.7071067811865475
120 -2.294044207478843 1.334041972098247 -0.7071067811865477
121 -1.953283852568877 1.932638001183191 -0.587785252292473
122 -0.6301002387974354 2.481399437255894 -0.8090169943751676
123 -2.485134578607252 0.617298646538557 -0.8090169943751676
124 -0.7429126824586613 1.179366135379631 -0.7518398074786405
125 -1.18231399064295 0.7535058031604537 -0.7589872812467737
126 -1.050969392999293 0.3788098644879105 -0.4002490705072542
127 -0.3769070614562796 1.050252288341556 -0.397748474526425
128 1.179366135379631 0.7429126824586614 -0.7518398074786404
129 0.7535058031604539 1.182313990642949 -0.7589872812467736
130 1.050252288341556 0.3769070614562796 -0.397748474526425
131 0.3788098644879104 1.050969392999293 -0.4002490705072542
132 0.7429126824586614 -1.179366135379631 -0.7518398074786404
133 1.182313990642949 -0.7535058031604539 -0.7589872812467736
134 1.050969392999293 -0.3788098644879104 -0.4002490705072542
135 0.3769070614562796 -1.050252288341556 -0.397748474526425
136 -1.179366135379631 -0.7429126824586614 -0.7518398074786404
137 -0.7535058031604539 -1.182313990642949 -0.7589872812467736
138 -1.050252288341556 -0.3769070614562796 -0.397748474526425
139 -0.3788098644879104 -1.050969392999293 -0.4002490705072542
140 -1.179366135380231 -0.7429126824590293 0.7518398074792839
141 -0.7535058031608894 -1.182313990643652 0.7589872812475151
142 -1.050252288342061 -0.3769070614564729 0.3977484745275441
143 -0.3788098644881118 -1.050969392999817 0.4002490705084094
144 0.7429126824590291 -1.179366135380232 0.751839807479284
145 1.182313990643652 -0.7535058031608894 0.7589872812475149
146 1.050969392999817 -0.3788098644881119 0.4002490705084094
147 0.376907061456473 -1.050252288342061 0.3977484745275444
148 1.179366135380232 0.7429126824590291 0.751839807479284
149 0.7535058031608894 1.182313990643652 0.7589872812475149
150 1.050252288342061 0.376907061456473 0.3977484745275444
151 0.3788098644881119 1.050969392999817 0.4002490705084094
152 -0.7429126824590291 1.179366135380232 0.751839807479284
153 -1.182313990643652 0.7535058031608894 0.7589872812475149
154 -1.050969392999817 0.3788098644881119 0.4002490705084094
155 -0.376907061456473 1.050252288342061 0.3977484745275444
$EndNodes
$Elements
411
1 15 2 0 1 1
2 15 2 0 2 2
3 15 2 0 3 3
4 15 2 0 4 4
5 15 2 0 5 5
6 15 2 0 6 6
7 15 2 0 7 7
8 15 2 0 8 8
9 15 2 0 9 9
10 15 2 0 10 10
11 15 2 0 11 11
12 15 2 0 12 12
13 15 2 0 13 13
14 15 2 0 14 14
15 15 2 0 15 15
16 15 2 0 16 16
17 15 2 0 17 17
18 15 2 0 18 18
19 15 2 0 19 19
20 15 2 0 20 20
21 15 2 0 21 21
22 15 2 0 22 22
23 15 2 0 23 23
24 1 2 0 1 7 24
25 1 2 0 1 24 14
26 1 2 0 2 14 25
27 1 2 0 2 25 6
28 1 2 0 3 6 26
29 1 2 0 3 26 15
30 1 2 0 4 15 27
31 1 2 0 4 27 7
32 1 2 0 5 9 28
33 1 2 0 5 28 16
34 1 2 0 6 16 29
35 1 2 0 6 29 8
36 1 2 0 7 8 30
37 1 2 0 7 30 17
38 1 2 0 8 17 31
39 1 2 0 8 31 9
40 1 2 0 9 11 32
41 1 2 0 9 32 18
42 1 2 0 10 18 33
43 1 2 0 10 33 10
44 1 2 0 11 10 34
45 1 2 0 11 34 19
46 1 2 0 12 19 35
47 1 2 0 12 35 11
48 1 2 0 13 13 36
49 1 2 0 13 36 20
50 1 2 0 14 20 37
51 1 2 0 14 37 12
52 1 2 0 15 12 38
53 1 2 0 15 38 21
54 1 2 0 16 21 39
55 1 2 0 16 39 13
56 1 2 0 17 12 40
57 1 2 0 17 40 6
58 1 2 0 18 6 41
59 1 2 0 18 41 10