From: Wolfgang Bangerth Date: Sat, 1 Mar 2008 13:45:53 +0000 (+0000) Subject: Add David Neckel's compressible Euler flow tutorial program. X-Git-Tag: v8.0.0~9337 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a98ab64fda76ed21f5acc880510b4a3da01794b0;p=dealii.git Add David Neckel's compressible Euler flow tutorial program. git-svn-id: https://svn.dealii.org/trunk@15823 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-33/Makefile b/deal.II/examples/step-33/Makefile new file mode 100644 index 0000000000..2102793621 --- /dev/null +++ b/deal.II/examples/step-33/Makefile @@ -0,0 +1,154 @@ +# $Id: Makefile 14008 2006-10-17 04:05:33Z bangerth $ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = $(basename $(shell echo step-*.cc)) + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov + + + + +# +# +# Usually, you will not need to change something beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox new file mode 100644 index 0000000000..e720513c32 --- /dev/null +++ b/deal.II/examples/step-33/doc/intro.dox @@ -0,0 +1,152 @@ +

Introduction

+ +This program was written for fun by David Neckels (NCAR) while working +at Sandia (on the Wyoming Express bus to and from Corrales each day). +The main purpose was to better understand Euler flow. +The code solves the basic Euler equations of gas dynamics, by using a +fully implicit Newton iteration (inspired by Sandia's Aria code). The +code may be configured by an input deck to run different simulations +on different meshes, with differing boundary conditions. + +The program also uses the Trilinos linear solvers (Aztec/Amesos) and +an automatic differentiation package, Sacado. + +

Euler flow

+ +The equations for a compressible, inviscid gas (the Euler equations) are +a basic system of conservation laws, in spatial dimension $d$, +@f[ +\partial_t \mathbf{w} + \nabla \cdot \mathbf{F}(\mathbf{w}) = \mathbf{0}, +@f] +with $\mathbf{w}=(\rho,\rho v_1,\dots,\rho v_d, E)^{\top}$ for $v_i$ equal to the +flow velocity in spatial direction $i$, $\rho$ equal to the fluid density, and +$E$ the energy of the gas. The flux matrix (or system of flux functions) +are defined such that the entire system of equations are +@f{eqnarray*} + \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0 \\ + \partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s + \delta_{is} p)}{\partial x_s} &=& 0, i=1,\dots,d \\ + \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=& 0 \\ + p = (\gamma -1)(E-\frac{1}{2} \rho |v|^2). +@f} + +For air, $\gamma=1.4$. + +

Discretization

+ +We choose a finite element space $V_h$, and integrate our conservation law against +our (vector) test function $\mathbf{v} \in V_h$. We integrate by parts and approximate the +boundary flux with a numerical flux $\mathbf{H}$, +@f{eqnarray*} +\int_{\Omega} (\partial_t \mathbf{w}, \mathbf{v}) + (\nabla \cdot \mathbf{F}(\mathbf{w}), \mathbf{v}) & \approx & \\ +\int_{\Omega} (\partial_t \mathbf{w}, \mathbf{v}) + (\mathbf{F}(\mathbf{w}), \nabla \mathbf{v}) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{v}) + \int_{\partial \Omega} (\mathbf{H}(\mathbf{w}^+, \mathbf{w}^-, \mathbf{n}), \mathbf{v}^+), +@f} +where $+$ is the interior trace of a function, and $-$ represents the outer trace. +The diffusion term $h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{v})$ is introduced strictly for stability, + where $h$ is the mesh size and $\eta$ is a parameter prescribing how much diffusion to add. +Depending on the boundary condition, we prescribe the outer trace $\mathbf{w}^-$, + + +For more information, please consult Ralf Hartmann's thesis TODO:Ref. + + +Our full discretization is thus +@f{eqnarray*} +R(\mathbf{W}_{n+1}) = \\ +\int_{\Omega} (\frac{\mathbf{w}_{n+1} - \mathbf{w}_n}{\delta t}, \mathbf{v}) + \int_{\Omega} (\mathbf{F}(\tilde{\mathbf{w}}), \mathbf{v}) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{v}) + \int_{\partial \Omega} (\mathbf{H}(\tilde{\mathbf{w}^+}), \mathbf{w}^-((\tilde{\mathbf{w}^+})), \mathbf{n}), \mathbf{v}) & = & 0 +@f} +where $\tilde{\mathbf{w}} = \theta \mathbf{w}_{n+1} + (1-\theta) \mathbf{w}_n$ for $0 \leq \theta \leq 1$ and +$\mathbf{w}_i = \sum_k \mathbf{W}_i^k \mathbf{\phi}_k$. + +We choose the Lax-Friedrich's flux, +$\mathbf{H}(\mathbf{a},\mathbf{b},\mathbf{n}) = \frac{1}{2}(\mathbf{F}(\mathbf{a})\cdot \mathbf{n} + \mathbf{F}(\mathbf{b})\cdot \mathbf{n} + \alpha (\mathbf{a} - \mathbf{b}))$. + +We solve the nonlinear system by a Newton iteration, i.e. by iterating +@f{eqnarray*} +\partial R(\mathbf{W}^k) \delta \mathbf{W} & = & - R(\mathbf{W}^{k}) \\ +\mathbf{W}^{k+1} &=& \mathbf{W}^k + \delta \mathbf{W}, +@f} +until $|R(\mathbf{W}^k)|$ (the residual) is sufficiently small. + +

Auto-Differentiation

+ +Since computing the Jacobian $\partial R$ is a terrible beast, we use an automatic differentiation package, +Sacado, to do this. Sacado is a C++ type that supports basic arithmetic operators and functions +such as sqrt, sin, cos, pow, etc. One declares a collection Sacado type and then declares +some of this collection as degrees of freedom. These variables are used in an algorithm, and as the variables +are used, their sensitivities with respect to these degrees of freedom are continuously updated. +One can imagine that for the full Jacobian, this could be prohibitively expensive. However, we do +not use the Sacado type for the entire computation, but only element by element. The author has +used this approach side by side with a hand coded Jacobian for the Incompressible Navier-Stokes problem +and found the Sacado approach to be just as fast as using a hand coded Jacobian. Since using the +auto-differentiation requires only that one code the Residual $R(\mathbf{W})$, ensuring code correctness +and maintaining code becomes tremendously more simple. We show a very simple Sacado example: + +@code +#include + +typedef Sacado::Fad::DFad fad_double; + +main() { + + fad_double a,b,c; + + a.diff(0,2); // Set a to be dof 0, in a 2-dof system. + + b.diff(1,2); // Set b to be dof 1, in a 2-dof system. + + a = 1; b = 2; + + c = 2*a+cos(a*b); + + double *derivs = c.fastAccessDx(0); // Access derivatives + + std::cout << "dc/da = " << derivs[0] << ", dc/db=" << derivs[1] << std::endl; + +} +@endcode + +It should be noted that Sacado provides more auto-differentation capabilities than the small subset +that I use in this program. However, if one understands the example above, they will understand +how I use Sacado in this Euler flow program. + +

Trilinos solvers

+The program uses either the Aztec iterative solvers, or the Amesos direct solver, as provided by +the Trilinos package. This package is inherently designed to be used in a parallel program, however, +it may be used in serial just as easily, as is done here. The Epetra package is the basic +vector/matrix library upon which the solvers are built. This very powerful package can be used +to describe the parallel distribution of a vector, and to define sparse matrices that operate +on these vectors. Please view the commented code for more details on how these solvers are used +within the example. + +

Adaptivity

+The example uses an ad-hoc refinement indicator that shows some usefullness in shock-type problems, and +in the downhill flow example included. We refine according to the squared gradient of the density. +Hanging nodes are handled by using computing the numerical flux across cells that are of differing +refinement levels. In this way, the example combines the continuous and DG methodologies. + +Further, we enforce a maximum number of refinement levels to keep refinement under check. It is the +author's experience that adaptivity in a time dependent problem, refinement can easily lead the simulation to +a screeching halt if care is not taken. The amount of refinement is +limited in the example by letting the user specify the +maximum level of refinement that will be present anywhere in the mesh. In this way, refinement +tends not to slow the simulation to a halt. This, of course, is purely a heuristic strategy, and +if the author's advisor heard about it, the author would likely be exiled forever from the finite + element error estimation community. + +

Input Deck

+ +We use an input file deck to drive the simulation. In this way, we can alter the boundary conditions +and other important properties of the simulation without having to recompile. For more information on +the format, look at the results section, where we describe an example file in more detail. + + + + diff --git a/deal.II/examples/step-33/doc/results.dox b/deal.II/examples/step-33/doc/results.dox new file mode 100644 index 0000000000..4e4ad15b43 --- /dev/null +++ b/deal.II/examples/step-33/doc/results.dox @@ -0,0 +1,134 @@ +

Results

+ +We run the problem with the mesh slide.inp and the following input deck: +@verbatim +# Listing of Parameters +# --------------------- + +# The input grid +set mesh=slide.inp + +# Stabilization parameter +set diffusion power = 2.0 + +# Scaled value for gravity. Positive means gravity points down. +set gravity = 1.0 + +# Boundary conditions +# We may specify boundary conditions for up to MAX_BD boundaries. +# Your .inp file should have these boundaries designated. +# +subsection boundary_1 + set no penetration = true # reflective boundary condition +end + +subsection boundary_2 +# outflow boundary + #set w_3 = pressure + #set w_3 value = 1.5 - y +end + +subsection boundary_3 + set no penetration = true # reflective + #set w_3 = pressure + #set w_3 value = 1.0 +end +subsection boundary_4 + set no penetration = true #reflective +end + +# Initial Conditions +# We set the initial conditions of the conservative variables. These lines +# are passed to the expression parsing function. You should use x,y,z for +# the coordinate variables. +# +subsection initial condition + set w_0 value = 0 + set w_1 value = 0 + set w_2 value = 10*(x<-0.7)*(y> 0.3)*(y< 0.45) + (1-(x<-0.7)*(y> 0.3)*(y< 0.45))*1.0 + set w_3 value = (1.5-(1.0*1.0*y))/0.4 +end + +# Time stepping control +# +subsection time stepping + set final time = 10.0 # simulation end time + set time step = 0.02 # simulation time step +end + +subsection linear solver + set output = quiet + set method = gmres + set ilut fill = 1.5 + set ilut drop tolerance = 1e-6 + set ilut absolute tolerance = 1e-6 + set ilut relative tolerance = 1.0 +end + +# Output frequency. +# You may wish to set this > time step if you dont want output at every step +subsection output + set step = 0.01 +end + +# Refinement control +subsection refinement + set refinement = shock # none only other option + set shock value = 1.5 + set shock levels = 1 # how many levels of refinement to allow +end + +# Flux parameters +subsection flux + set stab = mesh + #set stab value = 1.0 +end +@endverbatim + +When we run the program, we get the following kind of output: +@code +T=3.14, Number of active cells: 2617 + Number of degrees of freedom: 11184 +NonLin Res: Lin Iter Lin Res +______________________________________ +* 1.684e-02 0007 1.84e-13 +* 4.414e-05 0006 3.82e-15 +* 1.750e-09 0006 1.56e-19 +* 6.831e-16 0000 1.56e-19 +max_df:124 +T=3.16, Number of active cells: 2626 + Number of degrees of freedom: 11220 +NonLin Res: Lin Iter Lin Res +______________________________________ +* 1.673e-02 0007 9.64e-14 +* 4.230e-05 0006 3.83e-15 +* 8.748e-10 0006 7.67e-20 +* 6.804e-16 0000 7.67e-20 +max_df:124 +T=3.18, Number of active cells: 2644 + Number of degrees of freedom: 11296 +NonLin Res: Lin Iter Lin Res +______________________________________ +* 1.674e-02 0007 9.29e-14 +* 3.306e-05 0007 8.07e-17 +* 4.660e-10 0006 4.05e-20 +* 6.898e-16 0000 4.05e-20 +max_df:128 +T=3.2, Number of active cells: 2647 + Number of degrees of freedom: 11312 +NonLin Res: Lin Iter Lin Res +______________________________________ +... +@endcode + +This output reports the progress of the Newton iterations and the time stepping. + +@image html step-33.slide.gif + +As we see, the heavy mass of fluid is drawn down the slope by gravity, where it collides +with the ski lodge and is flung into the air! Hopefully everyone escapes! + +@image html step-33.slide_adapt.gif + +The adaptivity follows and preceeds the flow pattern, based on this heuristic refinement scheme. + diff --git a/deal.II/examples/step-33/slide.inp b/deal.II/examples/step-33/slide.inp new file mode 100644 index 0000000000..66684bc8b6 --- /dev/null +++ b/deal.II/examples/step-33/slide.inp @@ -0,0 +1,3560 @@ +1731 1827 0 0 0 +195 -0.964976 0.485465 0.000000 +133 -0.970149 0.500000 0.000000 +67 -1.000000 0.500000 0.000000 +180 -0.978723 0.478723 0.000000 +196 -0.945882 0.477585 0.000000 +132 -0.940299 0.500000 0.000000 +179 -0.957447 0.457447 0.000000 +197 -0.919731 0.466897 0.000000 +131 -0.910448 0.500000 0.000000 +178 -0.936170 0.436170 0.000000 +181 0.150000 -0.500000 0.000000 +187 0.150000 -0.471429 0.000000 +198 0.119014 -0.470800 0.000000 +185 0.120000 -0.500000 0.000000 +188 0.150000 -0.442857 0.000000 +199 0.117550 -0.441479 0.000000 +189 0.150000 -0.414286 0.000000 +200 0.117242 -0.412206 0.000000 +190 0.150000 -0.385714 0.000000 +201 0.118129 -0.383064 0.000000 +191 0.150000 -0.357143 0.000000 +202 0.120210 -0.354337 0.000000 +192 0.150000 -0.328571 0.000000 +203 0.121793 -0.325840 0.000000 +186 0.150000 -0.300000 0.000000 +204 0.122717 -0.296869 0.000000 +206 0.153091 -0.267708 0.000000 +205 0.123759 -0.266752 0.000000 +193 0.183333 -0.300000 0.000000 +207 0.184560 -0.267950 0.000000 +194 0.216667 -0.300000 0.000000 +208 0.216606 -0.268111 0.000000 +1 0.250000 -0.300000 0.000000 +209 0.248350 -0.268410 0.000000 +211 0.279437 -0.298917 0.000000 +210 0.278782 -0.268564 0.000000 +3 0.250000 -0.328571 0.000000 +212 0.279638 -0.328111 0.000000 +4 0.250000 -0.357143 0.000000 +213 0.279710 -0.356801 0.000000 +5 0.250000 -0.385714 0.000000 +214 0.279677 -0.385395 0.000000 +6 0.250000 -0.414286 0.000000 +215 0.279641 -0.413938 0.000000 +7 0.250000 -0.442857 0.000000 +216 0.279666 -0.442534 0.000000 +8 0.250000 -0.471429 0.000000 +217 0.279811 -0.471266 0.000000 +2 0.250000 -0.500000 0.000000 +10 0.280000 -0.500000 0.000000 +34 1.000000 0.500000 0.000000 +68 0.970149 0.500000 0.000000 +218 0.970096 0.469780 0.000000 +66 1.000000 0.469697 0.000000 +69 0.940299 0.500000 0.000000 +219 0.940231 0.469840 0.000000 +70 0.910448 0.500000 0.000000 +220 0.910371 0.469856 0.000000 +71 0.880597 0.500000 0.000000 +221 0.880527 0.469909 0.000000 +72 0.850746 0.500000 0.000000 +222 0.850684 0.469926 0.000000 +73 0.820896 0.500000 0.000000 +223 0.820842 0.469978 0.000000 +74 0.791045 0.500000 0.000000 +224 0.791001 0.470016 0.000000 +75 0.761194 0.500000 0.000000 +225 0.761159 0.470055 0.000000 +76 0.731343 0.500000 0.000000 +226 0.731316 0.470091 0.000000 +77 0.701493 0.500000 0.000000 +227 0.701472 0.470118 0.000000 +78 0.671642 0.500000 0.000000 +228 0.671627 0.470134 0.000000 +79 0.641791 0.500000 0.000000 +229 0.641780 0.470142 0.000000 +80 0.611940 0.500000 0.000000 +230 0.611934 0.470065 0.000000 +81 0.582090 0.500000 0.000000 +231 0.582089 0.470002 0.000000 +82 0.552239 0.500000 0.000000 +232 0.552248 0.469853 0.000000 +83 0.522388 0.500000 0.000000 +233 0.522413 0.469718 0.000000 +84 0.492537 0.500000 0.000000 +234 0.492577 0.469656 0.000000 +85 0.462687 0.500000 0.000000 +235 0.462743 0.469571 0.000000 +86 0.432836 0.500000 0.000000 +236 0.432895 0.469630 0.000000 +87 0.402985 0.500000 0.000000 +237 0.403046 0.469669 0.000000 +88 0.373134 0.500000 0.000000 +238 0.373191 0.469745 0.000000 +89 0.343284 0.500000 0.000000 +239 0.343335 0.469723 0.000000 +90 0.313433 0.500000 0.000000 +240 0.313478 0.469741 0.000000 +91 0.283582 0.500000 0.000000 +241 0.283622 0.469783 0.000000 +92 0.253731 0.500000 0.000000 +242 0.253768 0.469742 0.000000 +93 0.223881 0.500000 0.000000 +243 0.223918 0.469737 0.000000 +94 0.194030 0.500000 0.000000 +244 0.194064 0.469810 0.000000 +95 0.164179 0.500000 0.000000 +245 0.164210 0.469881 0.000000 +96 0.134328 0.500000 0.000000 +246 0.134355 0.469940 0.000000 +97 0.104478 0.500000 0.000000 +247 0.104500 0.469997 0.000000 +98 0.074627 0.500000 0.000000 +248 0.074645 0.470051 0.000000 +99 0.044776 0.500000 0.000000 +249 0.044790 0.470093 0.000000 +100 0.014925 0.500000 0.000000 +250 0.014935 0.470120 0.000000 +101 -0.014925 0.500000 0.000000 +251 -0.014918 0.470135 0.000000 +102 -0.044776 0.500000 0.000000 +252 -0.044771 0.470143 0.000000 +103 -0.074627 0.500000 0.000000 +253 -0.074623 0.470146 0.000000 +104 -0.104478 0.500000 0.000000 +254 -0.104496 0.469979 0.000000 +105 -0.134328 0.500000 0.000000 +255 -0.134359 0.469935 0.000000 +106 -0.164179 0.500000 0.000000 +256 -0.164225 0.469839 0.000000 +107 -0.194030 0.500000 0.000000 +257 -0.194083 0.469797 0.000000 +108 -0.223881 0.500000 0.000000 +258 -0.223940 0.469753 0.000000 +109 -0.253731 0.500000 0.000000 +259 -0.253798 0.469819 0.000000 +110 -0.283582 0.500000 0.000000 +260 -0.283647 0.470019 0.000000 +111 -0.313433 0.500000 0.000000 +261 -0.313491 0.470151 0.000000 +112 -0.343284 0.500000 0.000000 +262 -0.343348 0.470381 0.000000 +113 -0.373134 0.500000 0.000000 +263 -0.373195 0.470395 0.000000 +114 -0.402985 0.500000 0.000000 +264 -0.403049 0.470404 0.000000 +115 -0.432836 0.500000 0.000000 +265 -0.432937 0.470297 0.000000 +116 -0.462687 0.500000 0.000000 +266 -0.462848 0.469995 0.000000 +117 -0.492537 0.500000 0.000000 +267 -0.492807 0.469765 0.000000 +118 -0.522388 0.500000 0.000000 +268 -0.522760 0.469580 0.000000 +119 -0.552239 0.500000 0.000000 +269 -0.552730 0.469652 0.000000 +120 -0.582090 0.500000 0.000000 +270 -0.582662 0.469731 0.000000 +121 -0.611940 0.500000 0.000000 +271 -0.612594 0.469999 0.000000 +122 -0.641791 0.500000 0.000000 +272 -0.642647 0.470419 0.000000 +123 -0.671642 0.500000 0.000000 +273 -0.672638 0.470543 0.000000 +124 -0.701493 0.500000 0.000000 +274 -0.702968 0.469800 0.000000 +125 -0.731343 0.500000 0.000000 +275 -0.733232 0.468149 0.000000 +126 -0.761194 0.500000 0.000000 +276 -0.765399 0.466448 0.000000 +127 -0.791045 0.500000 0.000000 +277 -0.795702 0.469563 0.000000 +128 -0.820896 0.500000 0.000000 +278 -0.822797 0.472070 0.000000 +129 -0.850746 0.500000 0.000000 +279 -0.850723 0.467038 0.000000 +130 -0.880597 0.500000 0.000000 +303 -0.887028 0.454038 0.000000 +11 0.310000 -0.500000 0.000000 +280 0.309694 -0.471123 0.000000 +12 0.340000 -0.500000 0.000000 +281 0.339686 -0.471033 0.000000 +13 0.370000 -0.500000 0.000000 +282 0.369721 -0.471055 0.000000 +14 0.400000 -0.500000 0.000000 +283 0.399756 -0.470957 0.000000 +15 0.430000 -0.500000 0.000000 +284 0.429805 -0.470822 0.000000 +16 0.460000 -0.500000 0.000000 +285 0.459835 -0.470712 0.000000 +17 0.490000 -0.500000 0.000000 +286 0.489873 -0.470554 0.000000 +18 0.520000 -0.500000 0.000000 +287 0.519906 -0.470379 0.000000 +19 0.550000 -0.500000 0.000000 +288 0.549930 -0.470248 0.000000 +20 0.580000 -0.500000 0.000000 +289 0.579948 -0.470155 0.000000 +21 0.610000 -0.500000 0.000000 +290 0.609962 -0.470091 0.000000 +22 0.640000 -0.500000 0.000000 +291 0.639972 -0.470050 0.000000 +23 0.670000 -0.500000 0.000000 +292 0.669980 -0.470026 0.000000 +24 0.700000 -0.500000 0.000000 +293 0.699986 -0.470013 0.000000 +25 0.730000 -0.500000 0.000000 +294 0.729990 -0.470006 0.000000 +26 0.760000 -0.500000 0.000000 +295 0.759994 -0.470003 0.000000 +27 0.790000 -0.500000 0.000000 +296 0.789980 -0.469916 0.000000 +28 0.820000 -0.500000 0.000000 +297 0.819956 -0.469877 0.000000 +29 0.850000 -0.500000 0.000000 +298 0.849929 -0.469833 0.000000 +30 0.880000 -0.500000 0.000000 +299 0.879924 -0.469805 0.000000 +31 0.910000 -0.500000 0.000000 +300 0.909913 -0.469792 0.000000 +32 0.940000 -0.500000 0.000000 +301 0.939922 -0.469775 0.000000 +33 0.970000 -0.500000 0.000000 +302 0.969956 -0.469740 0.000000 +9 1.000000 -0.500000 0.000000 +35 1.000000 -0.469697 0.000000 +177 -0.914894 0.414894 0.000000 +176 -0.893617 0.393617 0.000000 +304 -0.867787 0.420844 0.000000 +175 -0.872340 0.372340 0.000000 +305 -0.847479 0.396796 0.000000 +174 -0.851064 0.351064 0.000000 +306 -0.826567 0.374510 0.000000 +173 -0.829787 0.329787 0.000000 +307 -0.806205 0.352141 0.000000 +172 -0.808511 0.308511 0.000000 +308 -0.786305 0.329477 0.000000 +171 -0.787234 0.287234 0.000000 +309 -0.766015 0.307496 0.000000 +170 -0.765957 0.265957 0.000000 +310 -0.745195 0.286345 0.000000 +169 -0.744681 0.244681 0.000000 +311 -0.723868 0.265704 0.000000 +168 -0.723404 0.223404 0.000000 +312 -0.702394 0.245107 0.000000 +167 -0.702128 0.202128 0.000000 +313 -0.680976 0.224300 0.000000 +166 -0.680851 0.180851 0.000000 +314 -0.659597 0.203254 0.000000 +165 -0.659574 0.159574 0.000000 +315 -0.638233 0.182030 0.000000 +164 -0.638298 0.138298 0.000000 +316 -0.616915 0.160630 0.000000 +163 -0.617021 0.117021 0.000000 +317 -0.595665 0.139192 0.000000 +162 -0.595745 0.095745 0.000000 +318 -0.574431 0.117764 0.000000 +161 -0.574468 0.074468 0.000000 +319 -0.553204 0.096309 0.000000 +160 -0.553191 0.053191 0.000000 +320 -0.531840 0.075023 0.000000 +159 -0.531915 0.031915 0.000000 +321 -0.510492 0.053733 0.000000 +158 -0.510638 0.010638 0.000000 +322 -0.489127 0.032461 0.000000 +157 -0.489362 -0.010638 0.000000 +323 -0.467791 0.011178 0.000000 +156 -0.468085 -0.031915 0.000000 +324 -0.446536 -0.010172 0.000000 +155 -0.446809 -0.053191 0.000000 +325 -0.425240 -0.031466 0.000000 +154 -0.425532 -0.074468 0.000000 +326 -0.404012 -0.052824 0.000000 +153 -0.404255 -0.095745 0.000000 +327 -0.382783 -0.074175 0.000000 +152 -0.382979 -0.117021 0.000000 +328 -0.361545 -0.095511 0.000000 +151 -0.361702 -0.138298 0.000000 +329 -0.340288 -0.116831 0.000000 +150 -0.340426 -0.159574 0.000000 +330 -0.319034 -0.138159 0.000000 +149 -0.319149 -0.180851 0.000000 +331 -0.297766 -0.159484 0.000000 +148 -0.297872 -0.202128 0.000000 +332 -0.276489 -0.180797 0.000000 +147 -0.276596 -0.223404 0.000000 +333 -0.255087 -0.202105 0.000000 +146 -0.255319 -0.244681 0.000000 +334 -0.233674 -0.223466 0.000000 +145 -0.234043 -0.265957 0.000000 +335 -0.212119 -0.244916 0.000000 +144 -0.212766 -0.287234 0.000000 +336 -0.190657 -0.266392 0.000000 +143 -0.191489 -0.308511 0.000000 +337 -0.169261 -0.287934 0.000000 +142 -0.170213 -0.329787 0.000000 +338 -0.147854 -0.309624 0.000000 +141 -0.148936 -0.351064 0.000000 +339 -0.126431 -0.331364 0.000000 +140 -0.127660 -0.372340 0.000000 +340 -0.105256 -0.353368 0.000000 +139 -0.106383 -0.393617 0.000000 +341 -0.083924 -0.375263 0.000000 +138 -0.085106 -0.414894 0.000000 +342 -0.062097 -0.396824 0.000000 +137 -0.063830 -0.436170 0.000000 +343 -0.040132 -0.418076 0.000000 +136 -0.042553 -0.457447 0.000000 +344 -0.019646 -0.439891 0.000000 +135 -0.021277 -0.478723 0.000000 +345 -0.002475 -0.461393 0.000000 +134 0.000000 -0.500000 0.000000 +346 0.009599 -0.478226 0.000000 +182 0.030000 -0.500000 0.000000 +347 0.029922 -0.474701 0.000000 +183 0.060000 -0.500000 0.000000 +348 0.057397 -0.471689 0.000000 +184 0.090000 -0.500000 0.000000 +349 0.087960 -0.470463 0.000000 +36 1.000000 -0.439394 0.000000 +350 0.969923 -0.439465 0.000000 +37 1.000000 -0.409091 0.000000 +351 0.969894 -0.409159 0.000000 +38 1.000000 -0.378788 0.000000 +352 0.969889 -0.378860 0.000000 +39 1.000000 -0.348485 0.000000 +353 0.969858 -0.348548 0.000000 +40 1.000000 -0.318182 0.000000 +354 0.969843 -0.318237 0.000000 +41 1.000000 -0.287879 0.000000 +355 0.969809 -0.287925 0.000000 +42 1.000000 -0.257576 0.000000 +356 0.969784 -0.257613 0.000000 +43 1.000000 -0.227273 0.000000 +357 0.969757 -0.227302 0.000000 +44 1.000000 -0.196970 0.000000 +358 0.969734 -0.196993 0.000000 +45 1.000000 -0.166667 0.000000 +359 0.969717 -0.166684 0.000000 +46 1.000000 -0.136364 0.000000 +360 0.969707 -0.136377 0.000000 +47 1.000000 -0.106061 0.000000 +361 0.969701 -0.106070 0.000000 +48 1.000000 -0.075758 0.000000 +362 0.969699 -0.075764 0.000000 +49 1.000000 -0.045455 0.000000 +363 0.969698 -0.045459 0.000000 +50 1.000000 -0.015152 0.000000 +364 0.969697 -0.015155 0.000000 +51 1.000000 0.015152 0.000000 +365 0.969697 0.015149 0.000000 +52 1.000000 0.045455 0.000000 +366 0.969697 0.045453 0.000000 +53 1.000000 0.075758 0.000000 +367 0.969697 0.075757 0.000000 +54 1.000000 0.106061 0.000000 +368 0.969697 0.106060 0.000000 +55 1.000000 0.136364 0.000000 +369 0.969697 0.136363 0.000000 +56 1.000000 0.166667 0.000000 +370 0.969697 0.166666 0.000000 +57 1.000000 0.196970 0.000000 +371 0.969697 0.196970 0.000000 +58 1.000000 0.227273 0.000000 +372 0.969697 0.227273 0.000000 +59 1.000000 0.257576 0.000000 +373 0.969795 0.257643 0.000000 +60 1.000000 0.287879 0.000000 +374 0.969874 0.287995 0.000000 +61 1.000000 0.318182 0.000000 +375 0.969958 0.318338 0.000000 +62 1.000000 0.348485 0.000000 +376 0.969973 0.348658 0.000000 +63 1.000000 0.378788 0.000000 +377 0.970000 0.378980 0.000000 +64 1.000000 0.409091 0.000000 +378 0.970010 0.409285 0.000000 +65 1.000000 0.439394 0.000000 +379 0.970062 0.439531 0.000000 +503 -0.841822 0.439272 0.000000 +380 -0.823383 0.418799 0.000000 +381 -0.801366 0.398075 0.000000 +382 -0.781899 0.373291 0.000000 +383 -0.763952 0.349541 0.000000 +384 -0.745019 0.327062 0.000000 +385 -0.724690 0.306382 0.000000 +386 -0.703265 0.286716 0.000000 +387 -0.681523 0.267067 0.000000 +388 -0.660018 0.246574 0.000000 +389 -0.638443 0.225771 0.000000 +390 -0.616833 0.204679 0.000000 +391 -0.595551 0.183139 0.000000 +392 -0.574409 0.161470 0.000000 +393 -0.553222 0.139792 0.000000 +394 -0.531865 0.118362 0.000000 +395 -0.510475 0.096979 0.000000 +396 -0.489009 0.075696 0.000000 +397 -0.467590 0.054392 0.000000 +398 -0.446270 0.033017 0.000000 +399 -0.424906 0.011735 0.000000 +400 -0.403621 -0.009619 0.000000 +401 -0.382400 -0.031023 0.000000 +402 -0.361173 -0.052410 0.000000 +403 -0.339922 -0.073800 0.000000 +404 -0.318715 -0.095237 0.000000 +405 -0.297453 -0.116646 0.000000 +406 -0.276156 -0.138018 0.000000 +407 -0.254723 -0.159496 0.000000 +408 -0.233335 -0.180914 0.000000 +409 -0.211596 -0.202409 0.000000 +410 -0.189882 -0.223912 0.000000 +411 -0.168328 -0.245541 0.000000 +412 -0.146829 -0.267432 0.000000 +413 -0.125347 -0.289413 0.000000 +414 -0.103926 -0.311617 0.000000 +415 -0.082894 -0.334299 0.000000 +416 -0.061472 -0.356952 0.000000 +417 -0.038982 -0.378982 0.000000 +418 -0.016521 -0.400436 0.000000 +419 0.006313 -0.423591 0.000000 +421 0.051537 -0.441884 0.000000 +420 0.021625 -0.449713 0.000000 +422 0.085122 -0.440903 0.000000 +423 0.082947 -0.409815 0.000000 +424 0.086501 -0.379229 0.000000 +425 0.090943 -0.350788 0.000000 +426 0.094089 -0.322730 0.000000 +427 0.095265 -0.294581 0.000000 +428 0.094987 -0.265825 0.000000 +430 0.123468 -0.236299 0.000000 +429 0.093679 -0.236645 0.000000 +431 0.153728 -0.236070 0.000000 +432 0.184970 -0.235916 0.000000 +433 0.216576 -0.236122 0.000000 +434 0.247804 -0.236880 0.000000 +435 0.278352 -0.237921 0.000000 +437 0.308849 -0.268892 0.000000 +436 0.308561 -0.238863 0.000000 +438 0.309174 -0.298672 0.000000 +439 0.309365 -0.327764 0.000000 +440 0.309409 -0.356452 0.000000 +441 0.309369 -0.385004 0.000000 +442 0.309363 -0.413591 0.000000 +443 0.309512 -0.442324 0.000000 +444 0.940150 0.439639 0.000000 +445 0.910293 0.439719 0.000000 +446 0.880439 0.439764 0.000000 +447 0.850606 0.439851 0.000000 +448 0.820773 0.439905 0.000000 +449 0.790942 0.439997 0.000000 +450 0.761110 0.440072 0.000000 +451 0.731277 0.440144 0.000000 +452 0.701441 0.440206 0.000000 +453 0.671603 0.440248 0.000000 +454 0.641769 0.440131 0.000000 +455 0.611935 0.440030 0.000000 +456 0.582116 0.439803 0.000000 +457 0.552311 0.439544 0.000000 +458 0.522496 0.439372 0.000000 +459 0.492687 0.439151 0.000000 +460 0.462850 0.439112 0.000000 +461 0.433008 0.439157 0.000000 +462 0.403152 0.439190 0.000000 +463 0.373296 0.439261 0.000000 +464 0.343430 0.439326 0.000000 +465 0.313568 0.439388 0.000000 +466 0.283709 0.439386 0.000000 +467 0.253856 0.439367 0.000000 +468 0.223995 0.439416 0.000000 +469 0.194134 0.439557 0.000000 +470 0.164274 0.439658 0.000000 +471 0.134410 0.439807 0.000000 +472 0.104545 0.439947 0.000000 +473 0.074682 0.440045 0.000000 +474 0.044821 0.440112 0.000000 +475 0.014962 0.440171 0.000000 +476 -0.014897 0.440222 0.000000 +477 -0.044754 0.440258 0.000000 +478 -0.074665 0.440098 0.000000 +479 -0.104552 0.439944 0.000000 +480 -0.134459 0.439734 0.000000 +481 -0.164338 0.439555 0.000000 +482 -0.194223 0.439404 0.000000 +483 -0.224113 0.439388 0.000000 +484 -0.253960 0.439636 0.000000 +485 -0.283803 0.439999 0.000000 +486 -0.313653 0.440461 0.000000 +487 -0.343481 0.440757 0.000000 +488 -0.373319 0.440824 0.000000 +489 -0.403191 0.440764 0.000000 +490 -0.433152 0.440474 0.000000 +491 -0.463150 0.439907 0.000000 +492 -0.493171 0.439333 0.000000 +493 -0.523226 0.438998 0.000000 +494 -0.553219 0.439109 0.000000 +495 -0.583244 0.439382 0.000000 +496 -0.613191 0.439846 0.000000 +497 -0.642956 0.440951 0.000000 +498 -0.672668 0.441524 0.000000 +499 -0.702061 0.440023 0.000000 +500 -0.733531 0.435747 0.000000 +501 -0.771806 0.426884 0.000000 +553 -0.803350 0.441135 0.000000 +502 -0.822712 0.450765 0.000000 +504 0.939858 -0.439508 0.000000 +505 0.939812 -0.409239 0.000000 +506 0.939773 -0.378933 0.000000 +507 0.939746 -0.348626 0.000000 +508 0.939690 -0.318305 0.000000 +509 0.939650 -0.287984 0.000000 +510 0.939589 -0.257662 0.000000 +511 0.939539 -0.227343 0.000000 +512 0.939492 -0.197025 0.000000 +513 0.939452 -0.166709 0.000000 +514 0.939425 -0.136395 0.000000 +515 0.939409 -0.106084 0.000000 +516 0.939401 -0.075774 0.000000 +517 0.939397 -0.045467 0.000000 +518 0.939395 -0.015160 0.000000 +519 0.939394 0.015146 0.000000 +520 0.939394 0.045451 0.000000 +521 0.939394 0.075755 0.000000 +522 0.939394 0.106059 0.000000 +523 0.939394 0.136362 0.000000 +524 0.939394 0.166666 0.000000 +525 0.939394 0.196969 0.000000 +526 0.939509 0.227355 0.000000 +527 0.939655 0.257739 0.000000 +528 0.939814 0.288143 0.000000 +529 0.939913 0.318502 0.000000 +530 0.939975 0.348849 0.000000 +531 0.940009 0.379160 0.000000 +532 0.940091 0.409414 0.000000 +533 0.339432 -0.442176 0.000000 +534 0.369451 -0.442035 0.000000 +535 0.399523 -0.441897 0.000000 +536 0.429589 -0.441684 0.000000 +537 0.459678 -0.441399 0.000000 +538 0.489748 -0.441062 0.000000 +539 0.519800 -0.440764 0.000000 +540 0.549844 -0.440522 0.000000 +541 0.579882 -0.440330 0.000000 +542 0.609913 -0.440196 0.000000 +543 0.639937 -0.440110 0.000000 +544 0.669955 -0.440059 0.000000 +545 0.699968 -0.440030 0.000000 +546 0.729978 -0.440014 0.000000 +547 0.759947 -0.439909 0.000000 +548 0.789937 -0.439803 0.000000 +549 0.819882 -0.439709 0.000000 +550 0.849843 -0.439644 0.000000 +551 0.879819 -0.439597 0.000000 +552 0.909827 -0.439560 0.000000 +590 0.040432 -0.405830 0.000000 +670 -0.755455 0.393182 0.000000 +554 -0.740317 0.367025 0.000000 +555 -0.724443 0.344933 0.000000 +556 -0.704864 0.325740 0.000000 +557 -0.683231 0.307949 0.000000 +558 -0.661196 0.289226 0.000000 +559 -0.639331 0.269137 0.000000 +560 -0.617409 0.248586 0.000000 +561 -0.595465 0.227639 0.000000 +562 -0.574245 0.205848 0.000000 +563 -0.553198 0.183785 0.000000 +564 -0.531973 0.162012 0.000000 +565 -0.510563 0.140460 0.000000 +566 -0.489049 0.119070 0.000000 +567 -0.467535 0.097696 0.000000 +568 -0.446086 0.076355 0.000000 +569 -0.424583 0.055151 0.000000 +570 -0.403151 0.033875 0.000000 +571 -0.381891 0.012449 0.000000 +572 -0.360645 -0.008987 0.000000 +573 -0.339410 -0.030497 0.000000 +574 -0.318218 -0.052043 0.000000 +575 -0.296914 -0.073549 0.000000 +576 -0.275621 -0.095064 0.000000 +577 -0.254107 -0.116660 0.000000 +578 -0.232698 -0.138235 0.000000 +579 -0.211016 -0.159840 0.000000 +580 -0.189307 -0.181488 0.000000 +581 -0.167499 -0.203012 0.000000 +582 -0.145658 -0.224717 0.000000 +583 -0.124325 -0.246942 0.000000 +584 -0.102865 -0.269240 0.000000 +585 -0.081284 -0.291521 0.000000 +586 -0.060814 -0.314997 0.000000 +587 -0.039868 -0.338366 0.000000 +588 -0.017866 -0.359895 0.000000 +589 0.007087 -0.381124 0.000000 +591 0.053625 -0.372025 0.000000 +592 0.062247 -0.344115 0.000000 +593 0.067251 -0.318078 0.000000 +594 0.068547 -0.292412 0.000000 +595 0.066596 -0.265244 0.000000 +596 0.064013 -0.237360 0.000000 +598 0.092058 -0.207144 0.000000 +597 0.061821 -0.208936 0.000000 +599 0.122455 -0.205662 0.000000 +600 0.153304 -0.204437 0.000000 +601 0.184778 -0.203622 0.000000 +602 0.216464 -0.203889 0.000000 +603 0.247666 -0.205271 0.000000 +604 0.278273 -0.207077 0.000000 +605 0.308444 -0.208860 0.000000 +607 0.338640 -0.239619 0.000000 +606 0.338501 -0.210215 0.000000 +608 0.338849 -0.269123 0.000000 +609 0.339046 -0.298501 0.000000 +610 0.339149 -0.327425 0.000000 +611 0.339145 -0.356028 0.000000 +612 0.339123 -0.384605 0.000000 +613 0.339265 -0.413346 0.000000 +614 0.910184 0.409525 0.000000 +615 0.880341 0.409627 0.000000 +616 0.850499 0.409699 0.000000 +617 0.820678 0.409821 0.000000 +618 0.790857 0.409913 0.000000 +619 0.761038 0.410043 0.000000 +620 0.731217 0.410148 0.000000 +621 0.701394 0.410247 0.000000 +622 0.671582 0.410190 0.000000 +623 0.641769 0.410064 0.000000 +624 0.611986 0.409820 0.000000 +625 0.582223 0.409499 0.000000 +626 0.552442 0.409166 0.000000 +627 0.522678 0.408817 0.000000 +628 0.492873 0.408592 0.000000 +629 0.463046 0.408514 0.000000 +630 0.433212 0.408428 0.000000 +631 0.403368 0.408486 0.000000 +632 0.373497 0.408591 0.000000 +633 0.343636 0.408704 0.000000 +634 0.313765 0.408808 0.000000 +635 0.283900 0.408804 0.000000 +636 0.254031 0.408855 0.000000 +637 0.224157 0.409018 0.000000 +638 0.194295 0.409133 0.000000 +639 0.164411 0.409379 0.000000 +640 0.134527 0.409617 0.000000 +641 0.104649 0.409788 0.000000 +642 0.074780 0.409914 0.000000 +643 0.044908 0.410046 0.000000 +644 0.015035 0.410160 0.000000 +645 -0.014837 0.410241 0.000000 +646 -0.044786 0.410157 0.000000 +647 -0.074715 0.410004 0.000000 +648 -0.104689 0.409724 0.000000 +649 -0.134629 0.409413 0.000000 +650 -0.164573 0.409102 0.000000 +651 -0.194543 0.408821 0.000000 +652 -0.224425 0.408905 0.000000 +653 -0.254313 0.409349 0.000000 +654 -0.284136 0.410098 0.000000 +655 -0.313927 0.410835 0.000000 +656 -0.343709 0.411217 0.000000 +657 -0.373516 0.411287 0.000000 +658 -0.403429 0.411111 0.000000 +659 -0.433418 0.410580 0.000000 +660 -0.463454 0.409742 0.000000 +661 -0.493588 0.408650 0.000000 +662 -0.523716 0.408015 0.000000 +663 -0.553867 0.408390 0.000000 +664 -0.583894 0.408825 0.000000 +665 -0.613736 0.409400 0.000000 +666 -0.642923 0.411762 0.000000 +667 -0.671337 0.413661 0.000000 +668 -0.699264 0.412614 0.000000 +669 -0.726520 0.405751 0.000000 +671 0.369212 -0.413104 0.000000 +672 0.399261 -0.412881 0.000000 +673 0.429386 -0.412541 0.000000 +674 0.459495 -0.412121 0.000000 +675 0.489590 -0.411630 0.000000 +676 0.519678 -0.411233 0.000000 +677 0.549757 -0.410839 0.000000 +678 0.579820 -0.410560 0.000000 +679 0.609867 -0.410349 0.000000 +680 0.639903 -0.410201 0.000000 +681 0.669930 -0.410109 0.000000 +682 0.699950 -0.410056 0.000000 +683 0.729910 -0.409917 0.000000 +684 0.759884 -0.409777 0.000000 +685 0.789837 -0.409624 0.000000 +686 0.819781 -0.409515 0.000000 +687 0.849736 -0.409429 0.000000 +688 0.879719 -0.409378 0.000000 +689 0.909756 -0.409298 0.000000 +690 0.909695 -0.379025 0.000000 +691 0.909632 -0.348708 0.000000 +692 0.909581 -0.318392 0.000000 +693 0.909498 -0.288059 0.000000 +694 0.909434 -0.257728 0.000000 +695 0.909350 -0.227398 0.000000 +696 0.909280 -0.197070 0.000000 +697 0.909216 -0.166745 0.000000 +698 0.909165 -0.136423 0.000000 +699 0.909131 -0.106105 0.000000 +700 0.909111 -0.075791 0.000000 +701 0.909100 -0.045478 0.000000 +702 0.909095 -0.015168 0.000000 +703 0.909093 0.015140 0.000000 +704 0.909072 0.045446 0.000000 +705 0.909067 0.075752 0.000000 +706 0.909073 0.106056 0.000000 +707 0.909080 0.136361 0.000000 +708 0.909086 0.166664 0.000000 +709 0.909215 0.197066 0.000000 +710 0.909385 0.227464 0.000000 +711 0.909604 0.257898 0.000000 +712 0.909779 0.288313 0.000000 +713 0.909917 0.318693 0.000000 +714 0.909986 0.349031 0.000000 +715 0.910101 0.379296 0.000000 +716 0.024590 -0.357790 0.000000 +717 0.036820 -0.334922 0.000000 +718 0.044018 -0.313436 0.000000 +719 0.043162 -0.290917 0.000000 +720 0.038764 -0.265690 0.000000 +721 0.034906 -0.239103 0.000000 +722 0.031729 -0.211118 0.000000 +724 0.059321 -0.180222 0.000000 +723 0.028990 -0.183087 0.000000 +725 0.089882 -0.177435 0.000000 +726 0.120857 -0.174840 0.000000 +727 0.152193 -0.172458 0.000000 +728 0.184221 -0.171049 0.000000 +729 0.216369 -0.171146 0.000000 +730 0.247776 -0.173286 0.000000 +731 0.278417 -0.176312 0.000000 +732 0.308664 -0.179103 0.000000 +733 0.338547 -0.181101 0.000000 +735 0.368461 -0.210821 0.000000 +734 0.368350 -0.182042 0.000000 +736 0.368647 -0.239838 0.000000 +737 0.368823 -0.269040 0.000000 +738 0.368963 -0.298178 0.000000 +739 0.368962 -0.326977 0.000000 +740 0.368943 -0.355559 0.000000 +741 0.369058 -0.384300 0.000000 +742 0.880206 0.379422 0.000000 +743 0.850374 0.379547 0.000000 +744 0.820545 0.379654 0.000000 +745 0.790741 0.379815 0.000000 +746 0.760936 0.379938 0.000000 +747 0.731131 0.380103 0.000000 +748 0.701348 0.380156 0.000000 +749 0.671569 0.380011 0.000000 +750 0.641834 0.379835 0.000000 +751 0.612110 0.379504 0.000000 +752 0.582403 0.379062 0.000000 +753 0.552705 0.378593 0.000000 +754 0.522968 0.378130 0.000000 +755 0.493191 0.377826 0.000000 +756 0.463421 0.377567 0.000000 +757 0.433611 0.377484 0.000000 +758 0.403752 0.377581 0.000000 +759 0.373895 0.377725 0.000000 +760 0.344005 0.377927 0.000000 +761 0.314132 0.378023 0.000000 +762 0.284270 0.378127 0.000000 +763 0.254375 0.378254 0.000000 +764 0.224504 0.378365 0.000000 +765 0.194586 0.378654 0.000000 +766 0.164677 0.379020 0.000000 +767 0.134778 0.379321 0.000000 +768 0.104899 0.379529 0.000000 +769 0.074998 0.379743 0.000000 +770 0.045096 0.379901 0.000000 +771 0.015198 0.380037 0.000000 +772 -0.014766 0.380029 0.000000 +773 -0.044758 0.379983 0.000000 +774 -0.074799 0.379773 0.000000 +775 -0.104855 0.379434 0.000000 +776 -0.134915 0.378975 0.000000 +777 -0.165034 0.378407 0.000000 +778 -0.195055 0.378036 0.000000 +779 -0.225064 0.378155 0.000000 +780 -0.254915 0.379009 0.000000 +781 -0.284684 0.380269 0.000000 +782 -0.314377 0.381326 0.000000 +783 -0.344063 0.381861 0.000000 +784 -0.373781 0.381860 0.000000 +785 -0.403636 0.381479 0.000000 +786 -0.433541 0.380809 0.000000 +787 -0.463618 0.379545 0.000000 +788 -0.493878 0.377655 0.000000 +789 -0.524338 0.376897 0.000000 +790 -0.554662 0.377591 0.000000 +791 -0.584521 0.378512 0.000000 +792 -0.614901 0.378305 0.000000 +793 -0.644852 0.383150 0.000000 +794 -0.671012 0.387415 0.000000 +795 -0.694363 0.388117 0.000000 +796 -0.716721 0.380258 0.000000 +797 0.399067 -0.383941 0.000000 +798 0.429159 -0.383498 0.000000 +799 0.459314 -0.382900 0.000000 +800 0.489457 -0.382364 0.000000 +801 0.519590 -0.381760 0.000000 +802 0.549686 -0.381311 0.000000 +803 0.579760 -0.380889 0.000000 +804 0.609819 -0.380588 0.000000 +805 0.639864 -0.380360 0.000000 +806 0.669900 -0.380203 0.000000 +807 0.699859 -0.379982 0.000000 +808 0.729822 -0.379783 0.000000 +809 0.759760 -0.379578 0.000000 +810 0.789702 -0.379418 0.000000 +811 0.819647 -0.379286 0.000000 +812 0.849617 -0.379205 0.000000 +813 0.879641 -0.379103 0.000000 +814 -0.704962 0.359795 0.000000 +815 -0.687415 0.344916 0.000000 +816 -0.665253 0.330013 0.000000 +817 -0.640678 0.312888 0.000000 +818 -0.618867 0.291654 0.000000 +819 -0.596620 0.271672 0.000000 +820 -0.574345 0.250779 0.000000 +821 -0.552849 0.228703 0.000000 +822 -0.532003 0.206108 0.000000 +823 -0.510760 0.184088 0.000000 +824 -0.489162 0.162574 0.000000 +825 -0.467617 0.141083 0.000000 +826 -0.446103 0.119735 0.000000 +827 -0.424459 0.098611 0.000000 +828 -0.402806 0.077505 0.000000 +829 -0.381337 0.056223 0.000000 +830 -0.359980 0.034804 0.000000 +831 -0.338720 0.013213 0.000000 +832 -0.317525 -0.008472 0.000000 +833 -0.296233 -0.030164 0.000000 +834 -0.274944 -0.051850 0.000000 +835 -0.253374 -0.073554 0.000000 +836 -0.231843 -0.095322 0.000000 +837 -0.210150 -0.117070 0.000000 +838 -0.188532 -0.138903 0.000000 +839 -0.167016 -0.160779 0.000000 +840 -0.144837 -0.182173 0.000000 +841 -0.122682 -0.203880 0.000000 +842 -0.101823 -0.226595 0.000000 +843 -0.081098 -0.248884 0.000000 +844 -0.058277 -0.270876 0.000000 +845 -0.037972 -0.296496 0.000000 +1723 -0.019092 -0.320042 0.000000 +1725 0.000735 -0.340540 0.000000 +848 0.879560 -0.348818 0.000000 +849 0.879472 -0.318487 0.000000 +850 0.879399 -0.288159 0.000000 +851 0.879291 -0.257815 0.000000 +852 0.879205 -0.227473 0.000000 +853 0.879100 -0.197133 0.000000 +854 0.879013 -0.166796 0.000000 +855 0.878936 -0.136464 0.000000 +856 0.878856 -0.106138 0.000000 +857 0.878795 -0.075817 0.000000 +858 0.878746 -0.045502 0.000000 +859 0.878725 -0.015190 0.000000 +860 0.878647 0.015119 0.000000 +861 0.878640 0.045426 0.000000 +862 0.878639 0.075733 0.000000 +863 0.878671 0.106040 0.000000 +864 0.878713 0.136347 0.000000 +865 0.878883 0.166766 0.000000 +866 0.879094 0.197181 0.000000 +867 0.879362 0.227644 0.000000 +868 0.879600 0.258094 0.000000 +869 0.879805 0.288518 0.000000 +870 0.879937 0.318888 0.000000 +871 0.880089 0.349176 0.000000 +846 0.020727 -0.293356 0.000000 +847 0.026838 -0.309682 0.000000 +872 0.012547 -0.270085 0.000000 +873 0.004882 -0.242473 0.000000 +874 0.002208 -0.213381 0.000000 +875 -0.000861 -0.185685 0.000000 +877 0.026244 -0.155284 0.000000 +876 -0.002631 -0.158212 0.000000 +878 0.056250 -0.151549 0.000000 +879 0.087018 -0.147720 0.000000 +880 0.118301 -0.143854 0.000000 +881 0.150274 -0.140350 0.000000 +882 0.183089 -0.137656 0.000000 +883 0.216064 -0.138164 0.000000 +884 0.248380 -0.141045 0.000000 +885 0.279257 -0.145764 0.000000 +886 0.309290 -0.149826 0.000000 +887 0.338917 -0.152404 0.000000 +888 0.368368 -0.153538 0.000000 +890 0.398029 -0.182007 0.000000 +889 0.397803 -0.153563 0.000000 +891 0.398367 -0.210697 0.000000 +892 0.398637 -0.239587 0.000000 +893 0.398861 -0.268621 0.000000 +894 0.398942 -0.297627 0.000000 +895 0.398900 -0.326384 0.000000 +896 0.398951 -0.355111 0.000000 +897 -0.691273 0.370756 0.000000 +898 -0.674493 0.364355 0.000000 +1040 -0.650656 0.354888 0.000000 +976 -0.617833 0.341638 0.000000 +899 -0.598837 0.312875 0.000000 +900 -0.576909 0.294643 0.000000 +901 -0.552428 0.275224 0.000000 +902 -0.530818 0.250993 0.000000 +903 -0.510792 0.227821 0.000000 +904 -0.489498 0.205790 0.000000 +905 -0.467802 0.184485 0.000000 +906 -0.446323 0.163103 0.000000 +907 -0.424713 0.141919 0.000000 +908 -0.402894 0.120971 0.000000 +909 -0.381057 0.100028 0.000000 +910 -0.359274 0.078918 0.000000 +911 -0.337845 0.057399 0.000000 +912 -0.316626 0.035561 0.000000 +913 -0.295418 0.013637 0.000000 +914 -0.274138 -0.008245 0.000000 +915 -0.252589 -0.030100 0.000000 +916 -0.230965 -0.052030 0.000000 +917 -0.209177 -0.074014 0.000000 +918 -0.187347 -0.095982 0.000000 +919 -0.165790 -0.118137 0.000000 +920 -0.144877 -0.140439 0.000000 +921 -0.122656 -0.161118 0.000000 +922 -0.098962 -0.182854 0.000000 +923 -0.079724 -0.207401 0.000000 +924 -0.061150 -0.228771 0.000000 +925 -0.031281 -0.247098 0.000000 +926 -0.013937 -0.279764 0.000000 +927 0.850209 0.349326 0.000000 +928 0.820397 0.349486 0.000000 +929 0.790584 0.349621 0.000000 +930 0.760799 0.349820 0.000000 +931 0.731036 0.349946 0.000000 +932 0.701303 0.349866 0.000000 +933 0.671610 0.349779 0.000000 +934 0.641932 0.349508 0.000000 +935 0.612330 0.349074 0.000000 +936 0.582719 0.348493 0.000000 +937 0.553104 0.347828 0.000000 +938 0.523441 0.347216 0.000000 +939 0.493783 0.346743 0.000000 +940 0.464050 0.346432 0.000000 +941 0.434236 0.346369 0.000000 +942 0.404391 0.346489 0.000000 +943 0.374528 0.346684 0.000000 +944 0.344649 0.346945 0.000000 +945 0.314797 0.347163 0.000000 +946 0.284881 0.347290 0.000000 +947 0.254985 0.347409 0.000000 +948 0.225047 0.347654 0.000000 +949 0.195107 0.348125 0.000000 +950 0.165173 0.348601 0.000000 +951 0.135276 0.348956 0.000000 +952 0.105340 0.349258 0.000000 +953 0.075403 0.349451 0.000000 +954 0.045475 0.349683 0.000000 +955 0.015506 0.349758 0.000000 +956 -0.014542 0.349888 0.000000 +957 -0.044633 0.349863 0.000000 +958 -0.074856 0.349598 0.000000 +959 -0.105110 0.349093 0.000000 +960 -0.135408 0.348361 0.000000 +961 -0.165677 0.347487 0.000000 +962 -0.195905 0.346884 0.000000 +963 -0.226027 0.347092 0.000000 +964 -0.255914 0.348502 0.000000 +965 -0.285596 0.350414 0.000000 +966 -0.315123 0.352026 0.000000 +967 -0.344543 0.352795 0.000000 +968 -0.374040 0.352571 0.000000 +969 -0.403574 0.351874 0.000000 +970 -0.433297 0.351311 0.000000 +971 -0.463227 0.349634 0.000000 +972 -0.493249 0.346733 0.000000 +973 -0.525130 0.344192 0.000000 +974 -0.556069 0.347386 0.000000 +975 -0.584638 0.349299 0.000000 +977 0.429047 -0.354564 0.000000 +978 0.459194 -0.353904 0.000000 +979 0.489384 -0.353164 0.000000 +980 0.519514 -0.352512 0.000000 +981 0.549620 -0.351866 0.000000 +982 0.579698 -0.351370 0.000000 +983 0.609762 -0.350919 0.000000 +984 0.639814 -0.350586 0.000000 +985 0.669775 -0.350217 0.000000 +986 0.699737 -0.349897 0.000000 +987 0.729661 -0.349597 0.000000 +988 0.759593 -0.349367 0.000000 +989 0.789532 -0.349176 0.000000 +990 0.819498 -0.349051 0.000000 +991 0.849515 -0.348917 0.000000 +992 0.849408 -0.318619 0.000000 +993 0.849295 -0.288272 0.000000 +994 0.849198 -0.257931 0.000000 +995 0.849065 -0.227574 0.000000 +996 0.848962 -0.197220 0.000000 +997 0.848836 -0.166870 0.000000 +998 0.848656 -0.136533 0.000000 +999 0.848506 -0.106206 0.000000 +1000 0.848367 -0.075894 0.000000 +1001 0.848296 -0.045579 0.000000 +1002 0.848134 -0.015275 0.000000 +1003 0.848083 0.015034 0.000000 +1004 0.848033 0.045346 0.000000 +1005 0.848096 0.075663 0.000000 +1006 0.848193 0.105984 0.000000 +1007 0.848437 0.136431 0.000000 +1008 0.848728 0.166875 0.000000 +1009 0.849076 0.197378 0.000000 +1010 0.849382 0.227866 0.000000 +1011 0.849654 0.258328 0.000000 +1012 0.849850 0.288735 0.000000 +1013 0.850055 0.319057 0.000000 +1014 -0.026262 -0.214203 0.000000 +1015 -0.029383 -0.187423 0.000000 +1016 -0.032824 -0.158817 0.000000 +1018 -0.003084 -0.132754 0.000000 +1017 -0.029552 -0.130056 0.000000 +1019 0.022484 -0.128711 0.000000 +1020 0.051752 -0.123251 0.000000 +1021 0.082401 -0.118191 0.000000 +1022 0.113973 -0.113014 0.000000 +1023 0.146621 -0.107888 0.000000 +1024 0.180362 -0.103934 0.000000 +1025 0.216332 -0.102481 0.000000 +1026 0.250760 -0.108923 0.000000 +1027 0.282088 -0.115954 0.000000 +1028 0.310848 -0.121366 0.000000 +1029 0.339795 -0.124483 0.000000 +1030 0.368570 -0.125488 0.000000 +1031 0.397434 -0.125259 0.000000 +1033 0.427150 -0.152702 0.000000 +1032 0.426416 -0.124184 0.000000 +1034 0.427780 -0.181299 0.000000 +1035 0.428281 -0.210062 0.000000 +1036 0.428674 -0.238927 0.000000 +1037 0.428954 -0.267922 0.000000 +1038 0.429022 -0.296882 0.000000 +1039 0.429009 -0.325704 0.000000 +1041 0.820193 0.319232 0.000000 +1042 0.790405 0.319425 0.000000 +1043 0.760627 0.319590 0.000000 +1044 0.730930 0.319648 0.000000 +1045 0.701264 0.319642 0.000000 +1046 0.671658 0.319493 0.000000 +1047 0.642145 0.319171 0.000000 +1048 0.612661 0.318582 0.000000 +1049 0.583227 0.317841 0.000000 +1050 0.553746 0.316951 0.000000 +1051 0.524252 0.316219 0.000000 +1052 0.494666 0.315548 0.000000 +1053 0.464954 0.315177 0.000000 +1054 0.435153 0.315100 0.000000 +1055 0.405377 0.315185 0.000000 +1056 0.375483 0.315504 0.000000 +1057 0.345638 0.315877 0.000000 +1058 0.315734 0.316161 0.000000 +1059 0.285817 0.316366 0.000000 +1060 0.255889 0.316515 0.000000 +1061 0.225925 0.316910 0.000000 +1062 0.195940 0.317535 0.000000 +1063 0.166024 0.318176 0.000000 +1064 0.136095 0.318655 0.000000 +1065 0.106116 0.318949 0.000000 +1066 0.076161 0.319226 0.000000 +1067 0.046155 0.319410 0.000000 +1068 0.016077 0.319689 0.000000 +1069 -0.014080 0.319906 0.000000 +1070 -0.044435 0.319837 0.000000 +1071 -0.074887 0.319446 0.000000 +1072 -0.105364 0.318757 0.000000 +1073 -0.135875 0.317703 0.000000 +1074 -0.166440 0.316418 0.000000 +1075 -0.197098 0.315318 0.000000 +1076 -0.227516 0.315492 0.000000 +1077 -0.257562 0.317644 0.000000 +1078 -0.287116 0.320624 0.000000 +1079 -0.316202 0.323175 0.000000 +1080 -0.345132 0.324256 0.000000 +1081 -0.374089 0.323465 0.000000 +1082 -0.403254 0.322262 0.000000 +1083 -0.432294 0.322314 0.000000 +1084 -0.460563 0.320631 0.000000 +1085 -0.490322 0.315176 0.000000 +1086 -0.527019 0.305092 0.000000 +1087 -0.560434 0.318601 0.000000 +1088 -0.581094 0.326893 0.000000 +1089 0.459215 -0.324955 0.000000 +1090 0.489364 -0.324151 0.000000 +1091 0.519476 -0.323335 0.000000 +1092 0.549546 -0.322613 0.000000 +1093 0.579618 -0.321924 0.000000 +1094 0.609675 -0.321348 0.000000 +1095 0.639635 -0.320779 0.000000 +1096 0.669600 -0.320274 0.000000 +1097 0.699520 -0.319810 0.000000 +1098 0.729446 -0.319441 0.000000 +1099 0.759383 -0.319141 0.000000 +1100 0.789349 -0.318940 0.000000 +1101 0.819365 -0.318752 0.000000 +1213 -0.506862 0.272229 0.000000 +1102 -0.488578 0.247543 0.000000 +1103 -0.468378 0.226625 0.000000 +1104 -0.446529 0.206460 0.000000 +1105 -0.425174 0.185182 0.000000 +1106 -0.403551 0.164089 0.000000 +1107 -0.381461 0.143421 0.000000 +1108 -0.359047 0.122895 0.000000 +1109 -0.336992 0.101927 0.000000 +1110 -0.315505 0.080197 0.000000 +1111 -0.294359 0.058038 0.000000 +1112 -0.273215 0.035830 0.000000 +1113 -0.251740 0.013783 0.000000 +1114 -0.230105 -0.008250 0.000000 +1115 -0.208280 -0.030424 0.000000 +1116 -0.186340 -0.052658 0.000000 +1117 -0.164096 -0.074999 0.000000 +1118 -0.142249 -0.098075 0.000000 +1119 -0.124600 -0.120870 0.000000 +1120 -0.101689 -0.138969 0.000000 +1121 -0.070069 -0.159405 0.000000 +1143 -0.056720 -0.191602 0.000000 +1122 -0.048349 -0.211338 0.000000 +1123 0.819230 -0.288441 0.000000 +1124 0.819089 -0.258075 0.000000 +1125 0.818971 -0.227720 0.000000 +1126 0.818800 -0.197352 0.000000 +1127 0.818537 -0.167025 0.000000 +1128 0.818297 -0.136701 0.000000 +1129 0.818075 -0.106416 0.000000 +1130 0.817914 -0.076099 0.000000 +1131 0.817666 -0.045814 0.000000 +1132 0.817482 -0.015504 0.000000 +1133 0.817346 0.014821 0.000000 +1134 0.817349 0.045170 0.000000 +1135 0.817459 0.075531 0.000000 +1136 0.817790 0.106031 0.000000 +1137 0.818215 0.136528 0.000000 +1138 0.818711 0.167083 0.000000 +1139 0.819128 0.197617 0.000000 +1140 0.819470 0.228123 0.000000 +1141 0.819728 0.258571 0.000000 +1142 0.819993 0.288931 0.000000 +1244 -0.059543 -0.125203 0.000000 +1145 -0.022658 -0.104274 0.000000 +1144 -0.050604 -0.097371 0.000000 +1146 -0.003277 -0.113375 0.000000 +1147 0.015569 -0.103577 0.000000 +1148 0.043900 -0.096090 0.000000 +1149 0.074890 -0.089260 0.000000 +1150 0.106866 -0.082653 0.000000 +1151 0.139199 -0.075782 0.000000 +1152 0.174640 -0.068119 0.000000 +1153 0.218123 -0.058538 0.000000 +1154 0.257570 -0.075952 0.000000 +1155 0.287428 -0.087324 0.000000 +1156 0.314488 -0.094635 0.000000 +1157 0.341036 -0.097765 0.000000 +1158 0.368607 -0.097995 0.000000 +1159 0.396531 -0.096956 0.000000 +1160 0.424895 -0.095330 0.000000 +1162 0.455499 -0.122538 0.000000 +1161 0.453743 -0.093457 0.000000 +1163 0.456665 -0.151337 0.000000 +1164 0.457569 -0.180155 0.000000 +1165 0.458233 -0.209057 0.000000 +1166 0.458740 -0.238043 0.000000 +1167 0.459073 -0.267080 0.000000 +1168 0.459174 -0.296028 0.000000 +1169 0.790159 0.289139 0.000000 +1170 0.760439 0.289302 0.000000 +1171 0.730789 0.289384 0.000000 +1172 0.701251 0.289345 0.000000 +1173 0.671804 0.289210 0.000000 +1174 0.642476 0.288781 0.000000 +1175 0.613248 0.288165 0.000000 +1176 0.584018 0.287164 0.000000 +1177 0.554753 0.286148 0.000000 +1178 0.525405 0.285084 0.000000 +1179 0.495870 0.284318 0.000000 +1180 0.466191 0.283786 0.000000 +1181 0.436505 0.283622 0.000000 +1182 0.406661 0.283753 0.000000 +1183 0.376803 0.284156 0.000000 +1184 0.346949 0.284678 0.000000 +1185 0.317067 0.285146 0.000000 +1186 0.287175 0.285413 0.000000 +1187 0.257218 0.285557 0.000000 +1188 0.227216 0.286008 0.000000 +1189 0.197271 0.286938 0.000000 +1190 0.167376 0.287870 0.000000 +1191 0.137356 0.288464 0.000000 +1192 0.107379 0.288752 0.000000 +1193 0.077319 0.288996 0.000000 +1194 0.047198 0.289325 0.000000 +1195 0.016953 0.289780 0.000000 +1196 -0.013452 0.290044 0.000000 +1197 -0.044033 0.289941 0.000000 +1198 -0.074662 0.289458 0.000000 +1199 -0.105330 0.288557 0.000000 +1200 -0.136186 0.287158 0.000000 +1201 -0.167320 0.285224 0.000000 +1202 -0.198610 0.283265 0.000000 +1203 -0.229792 0.283738 0.000000 +1204 -0.260249 0.286423 0.000000 +1205 -0.289549 0.291136 0.000000 +1206 -0.317900 0.295204 0.000000 +1207 -0.345821 0.296788 0.000000 +1208 -0.373271 0.294906 0.000000 +1209 -0.402699 0.292329 0.000000 +1210 -0.431306 0.294539 0.000000 +1211 -0.456225 0.294743 0.000000 +1212 -0.480543 0.286589 0.000000 +1214 0.489324 -0.295158 0.000000 +1215 0.519386 -0.294275 0.000000 +1216 0.549442 -0.293393 0.000000 +1217 0.579473 -0.292561 0.000000 +1218 0.609427 -0.291788 0.000000 +1219 0.639393 -0.291062 0.000000 +1220 0.669310 -0.290383 0.000000 +1221 0.699237 -0.289798 0.000000 +1222 0.729178 -0.289297 0.000000 +1223 0.759153 -0.288937 0.000000 +1224 0.789181 -0.288640 0.000000 +1225 -0.468048 0.263946 0.000000 +1226 -0.449064 0.247078 0.000000 +1227 -0.425672 0.229002 0.000000 +1228 -0.403972 0.207099 0.000000 +1229 -0.382701 0.186041 0.000000 +1230 -0.360028 0.165890 0.000000 +1231 -0.336966 0.146039 0.000000 +1232 -0.314349 0.125401 0.000000 +1233 -0.292863 0.103310 0.000000 +1234 -0.271955 0.080617 0.000000 +1235 -0.250761 0.058067 0.000000 +1236 -0.229184 0.035846 0.000000 +1237 -0.207348 0.013734 0.000000 +1238 -0.185635 -0.008675 0.000000 +1239 -0.163704 -0.031228 0.000000 +1240 -0.140931 -0.053778 0.000000 +1241 -0.113524 -0.078759 0.000000 +1242 -0.108485 -0.104935 0.000000 +1243 -0.087393 -0.114472 0.000000 +1245 0.789012 -0.258320 0.000000 +1246 0.788816 -0.227953 0.000000 +1247 0.788528 -0.197675 0.000000 +1248 0.788213 -0.167365 0.000000 +1249 0.787924 -0.137121 0.000000 +1250 0.787641 -0.106821 0.000000 +1251 0.787323 -0.076569 0.000000 +1252 0.787001 -0.046263 0.000000 +1253 0.786721 -0.015923 0.000000 +1254 0.786547 0.014478 0.000000 +1255 0.786533 0.044918 0.000000 +1256 0.786809 0.075509 0.000000 +1257 0.787410 0.106098 0.000000 +1258 0.788167 0.136724 0.000000 +1259 0.788811 0.167320 0.000000 +1260 0.789306 0.197873 0.000000 +1261 0.789646 0.228368 0.000000 +1262 0.789928 0.258780 0.000000 +1263 -0.079885 -0.088544 0.000000 +1264 -0.068853 -0.063035 0.000000 +1354 -0.096410 -0.050332 0.000000 +1265 -0.040079 -0.071151 0.000000 +1266 -0.003904 -0.078351 0.000000 +1267 0.033276 -0.068635 0.000000 +1268 0.065450 -0.060908 0.000000 +1269 0.097012 -0.053165 0.000000 +1270 0.129040 -0.044661 0.000000 +1271 0.161078 -0.033951 0.000000 +1272 0.194763 -0.019351 0.000000 +1274 0.297376 -0.062036 0.000000 +1273 0.273270 -0.047050 0.000000 +1730 0.318766 -0.071448 0.000000 +1275 0.342006 -0.072798 0.000000 +1276 0.367455 -0.071083 0.000000 +1277 0.394160 -0.068575 0.000000 +1278 0.422042 -0.066164 0.000000 +1279 0.450643 -0.063564 0.000000 +1281 0.482961 -0.091446 0.000000 +1280 0.480092 -0.061554 0.000000 +1282 0.484919 -0.120729 0.000000 +1283 0.486360 -0.149771 0.000000 +1284 0.487432 -0.178843 0.000000 +1285 0.488230 -0.207943 0.000000 +1286 0.488820 -0.237072 0.000000 +1287 0.489125 -0.266135 0.000000 +1288 0.247345 -0.025834 0.000000 +1289 0.224714 0.002070 0.000000 +1290 0.760232 0.258961 0.000000 +1291 0.730702 0.259023 0.000000 +1292 0.701273 0.259041 0.000000 +1293 0.672070 0.258884 0.000000 +1294 0.643041 0.258580 0.000000 +1295 0.614100 0.257855 0.000000 +1296 0.585158 0.256709 0.000000 +1297 0.556141 0.255290 0.000000 +1298 0.526914 0.253967 0.000000 +1299 0.497488 0.252938 0.000000 +1300 0.467901 0.252297 0.000000 +1301 0.438185 0.251985 0.000000 +1302 0.408350 0.252105 0.000000 +1303 0.378486 0.252635 0.000000 +1304 0.348689 0.253392 0.000000 +1305 0.318849 0.254180 0.000000 +1306 0.288993 0.254506 0.000000 +1307 0.259104 0.254448 0.000000 +1308 0.229125 0.254911 0.000000 +1309 0.199160 0.256428 0.000000 +1310 0.169139 0.257831 0.000000 +1311 0.139113 0.258524 0.000000 +1312 0.109056 0.258739 0.000000 +1313 0.078905 0.258871 0.000000 +1314 0.048551 0.259413 0.000000 +1315 0.018071 0.260052 0.000000 +1316 -0.012555 0.260381 0.000000 +1317 -0.043290 0.260233 0.000000 +1318 -0.074004 0.259633 0.000000 +1319 -0.104865 0.258556 0.000000 +1320 -0.136045 0.256823 0.000000 +1321 -0.167724 0.254059 0.000000 +1322 -0.199349 0.251178 0.000000 +1323 -0.232977 0.249469 0.000000 +1324 -0.265381 0.255203 0.000000 +1325 -0.294998 0.262571 0.000000 +1326 -0.321525 0.268904 0.000000 +1327 -0.345370 0.271940 0.000000 +1328 -0.370494 0.266937 0.000000 +1329 -0.403103 0.256984 0.000000 +1330 -0.434068 0.268653 0.000000 +1331 0.519208 -0.265168 0.000000 +1332 0.549223 -0.264174 0.000000 +1333 0.579169 -0.263217 0.000000 +1334 0.609116 -0.262295 0.000000 +1335 0.639018 -0.261389 0.000000 +1336 0.668940 -0.260559 0.000000 +1337 0.698881 -0.259811 0.000000 +1338 0.728862 -0.259217 0.000000 +1339 0.758922 -0.258684 0.000000 +1340 -0.452781 0.275781 0.000000 +1341 -0.382340 0.228180 0.000000 +1342 -0.361422 0.207003 0.000000 +1343 -0.338786 0.188092 0.000000 +1344 -0.314784 0.169671 0.000000 +1345 -0.291764 0.149035 0.000000 +1346 -0.269929 0.126678 0.000000 +1347 -0.249214 0.103062 0.000000 +1348 -0.228144 0.080017 0.000000 +1349 -0.206339 0.057748 0.000000 +1350 -0.184396 0.035681 0.000000 +1351 -0.163152 0.013135 0.000000 +1352 -0.142212 -0.009361 0.000000 +1353 -0.120583 -0.030889 0.000000 +1355 0.758631 -0.228471 0.000000 +1356 0.758281 -0.198203 0.000000 +1357 0.757927 -0.168007 0.000000 +1358 0.757568 -0.137774 0.000000 +1359 0.757184 -0.107553 0.000000 +1360 0.756766 -0.077283 0.000000 +1361 0.756298 -0.046951 0.000000 +1362 0.755911 -0.016494 0.000000 +1363 0.755578 0.014071 0.000000 +1364 0.755530 0.044794 0.000000 +1365 0.756009 0.075565 0.000000 +1366 0.757080 0.106305 0.000000 +1367 0.758211 0.136989 0.000000 +1368 0.759083 0.167591 0.000000 +1369 0.759636 0.198105 0.000000 +1370 0.760021 0.228532 0.000000 +1371 -0.360079 0.243309 0.000000 +1372 -0.341982 0.226476 0.000000 +1373 -0.319699 0.210923 0.000000 +1374 -0.294047 0.194220 0.000000 +1375 -0.267776 0.174727 0.000000 +1376 -0.245955 0.149439 0.000000 +1377 -0.226125 0.124759 0.000000 +1378 -0.205333 0.101062 0.000000 +1379 -0.183133 0.079096 0.000000 +1380 -0.161064 0.056981 0.000000 +1381 -0.140320 0.034408 0.000000 +1382 -0.120789 0.011646 0.000000 +1383 -0.101307 -0.009581 0.000000 +1384 -0.080345 -0.027959 0.000000 +1385 -0.059504 -0.040679 0.000000 +1386 -0.035185 -0.045016 0.000000 +1387 -0.006250 -0.045396 0.000000 +1388 0.025205 -0.040029 0.000000 +1389 0.055491 -0.032869 0.000000 +1390 0.085590 -0.024698 0.000000 +1391 0.115671 -0.015090 0.000000 +1392 0.145238 -0.002575 0.000000 +1393 0.175423 0.012389 0.000000 +1394 0.204047 0.030821 0.000000 +1396 0.252433 0.023671 0.000000 +1395 0.231976 0.050390 0.000000 +1397 0.272784 -0.001727 0.000000 +1398 0.292997 -0.024140 0.000000 +1399 -0.344301 0.253627 0.000000 +1476 -0.327649 0.245851 0.000000 +1475 -0.303967 0.235261 0.000000 +1474 -0.275606 0.222900 0.000000 +1473 -0.238652 0.207257 0.000000 +1400 -0.218940 0.171385 0.000000 +1401 -0.201026 0.143814 0.000000 +1402 -0.182515 0.120273 0.000000 +1403 -0.160481 0.099745 0.000000 +1404 -0.136284 0.078764 0.000000 +1405 -0.115924 0.053849 0.000000 +1406 -0.097288 0.029482 0.000000 +1407 -0.080305 0.007789 0.000000 +1408 -0.064615 -0.010938 0.000000 +1409 -0.053222 -0.024526 0.000000 +1410 -0.035546 -0.021961 0.000000 +1411 -0.011462 -0.018372 0.000000 +1412 0.015936 -0.012857 0.000000 +1413 0.044720 -0.005800 0.000000 +1414 0.073013 0.002618 0.000000 +1415 0.101137 0.013052 0.000000 +1416 0.129049 0.026082 0.000000 +1417 0.156365 0.041756 0.000000 +1418 0.184159 0.058506 0.000000 +1419 0.210898 0.077145 0.000000 +1421 0.258587 0.071409 0.000000 +1420 0.235912 0.097537 0.000000 +1422 0.279603 0.044767 0.000000 +1423 0.299087 0.018547 0.000000 +1424 0.315969 -0.006303 0.000000 +1731 0.310464 -0.042060 0.000000 +1425 0.330327 -0.029348 0.000000 +1428 0.363244 -0.044857 0.000000 +1427 0.340084 -0.050081 0.000000 +1429 0.389469 -0.040337 0.000000 +1430 0.417204 -0.036504 0.000000 +1431 0.446200 -0.033641 0.000000 +1432 0.475860 -0.030939 0.000000 +1434 0.509909 -0.059392 0.000000 +1433 0.506170 -0.028724 0.000000 +1435 0.512650 -0.089356 0.000000 +1436 0.514653 -0.118826 0.000000 +1437 0.516198 -0.148172 0.000000 +1438 0.517372 -0.177496 0.000000 +1439 0.518251 -0.206797 0.000000 +1440 0.518809 -0.236000 0.000000 +1441 0.730529 0.228662 0.000000 +1442 0.701323 0.228754 0.000000 +1443 0.672423 0.228800 0.000000 +1444 0.643741 0.228640 0.000000 +1445 0.615194 0.227932 0.000000 +1446 0.586654 0.226401 0.000000 +1447 0.557941 0.224612 0.000000 +1448 0.528929 0.223020 0.000000 +1449 0.499611 0.221647 0.000000 +1450 0.470054 0.220728 0.000000 +1451 0.440302 0.220268 0.000000 +1452 0.410353 0.220340 0.000000 +1453 0.380500 0.220966 0.000000 +1454 0.350669 0.222138 0.000000 +1455 0.320983 0.223466 0.000000 +1456 0.291443 0.223868 0.000000 +1457 0.261694 0.223241 0.000000 +1458 0.231496 0.223560 0.000000 +1459 0.201231 0.226194 0.000000 +1460 0.171137 0.228421 0.000000 +1461 0.141193 0.229123 0.000000 +1462 0.111175 0.228865 0.000000 +1463 0.080816 0.228761 0.000000 +1464 0.050208 0.229564 0.000000 +1465 0.019472 0.230514 0.000000 +1466 -0.011343 0.230954 0.000000 +1467 -0.042117 0.230757 0.000000 +1468 -0.072814 0.229991 0.000000 +1469 -0.103589 0.228826 0.000000 +1470 -0.134529 0.227009 0.000000 +1471 -0.165046 0.223412 0.000000 +1472 -0.197706 0.217242 0.000000 +1477 0.548816 -0.234928 0.000000 +1478 0.578757 -0.233868 0.000000 +1479 0.608646 -0.232776 0.000000 +1480 0.638544 -0.231722 0.000000 +1481 0.668475 -0.230745 0.000000 +1482 0.698449 -0.229883 0.000000 +1483 0.728502 -0.229106 0.000000 +1484 0.728097 -0.198990 0.000000 +1485 0.727706 -0.168870 0.000000 +1486 0.727268 -0.138716 0.000000 +1487 0.726781 -0.108502 0.000000 +1488 0.726187 -0.078232 0.000000 +1489 0.725568 -0.047739 0.000000 +1490 0.724892 -0.017056 0.000000 +1491 0.724335 0.013822 0.000000 +1492 0.724164 0.044888 0.000000 +1493 0.724932 0.075908 0.000000 +1494 0.726778 0.106759 0.000000 +1495 0.728488 0.137429 0.000000 +1496 0.729633 0.167949 0.000000 +1497 0.730244 0.198323 0.000000 +1624 -0.176223 0.159072 0.000000 +1565 -0.188636 0.185770 0.000000 +1498 -0.161949 0.137398 0.000000 +1499 -0.140110 0.121576 0.000000 +1500 -0.106545 0.104454 0.000000 +1501 -0.087526 0.070737 0.000000 +1502 -0.072361 0.044417 0.000000 +1503 -0.057196 0.020338 0.000000 +1504 -0.044389 -0.001210 0.000000 +1505 -0.020536 0.005897 0.000000 +1506 0.005896 0.012768 0.000000 +1507 0.032832 0.019947 0.000000 +1508 0.059602 0.028562 0.000000 +1509 0.085773 0.039645 0.000000 +1510 0.111268 0.053907 0.000000 +1511 0.136949 0.069672 0.000000 +1512 0.163708 0.085560 0.000000 +1513 0.189111 0.102843 0.000000 +1514 0.212292 0.123719 0.000000 +1516 0.261207 0.118056 0.000000 +1515 0.235876 0.152123 0.000000 +1517 0.283993 0.091978 0.000000 +1518 0.308470 0.065373 0.000000 +1519 0.327160 0.035468 0.000000 +1520 0.342353 0.007735 0.000000 +1522 0.381592 -0.011628 0.000000 +1521 0.354627 -0.019214 0.000000 +1523 0.410351 -0.006674 0.000000 +1524 0.440111 -0.003133 0.000000 +1525 0.470506 -0.000324 0.000000 +1526 0.501426 0.002030 0.000000 +1528 0.536865 -0.026396 0.000000 +1527 0.532740 0.004491 0.000000 +1529 0.540168 -0.057063 0.000000 +1530 0.542678 -0.087141 0.000000 +1531 0.544636 -0.116906 0.000000 +1532 0.546173 -0.146556 0.000000 +1533 0.547355 -0.176132 0.000000 +1534 0.548177 -0.205590 0.000000 +1535 0.701238 0.198603 0.000000 +1536 0.672622 0.198890 0.000000 +1537 0.644404 0.199062 0.000000 +1538 0.616485 0.198422 0.000000 +1539 0.588583 0.196609 0.000000 +1540 0.560344 0.194273 0.000000 +1541 0.531588 0.192198 0.000000 +1542 0.502343 0.190448 0.000000 +1543 0.472715 0.189284 0.000000 +1544 0.442689 0.188639 0.000000 +1545 0.412525 0.188524 0.000000 +1546 0.382371 0.189211 0.000000 +1547 0.352532 0.190926 0.000000 +1548 0.323383 0.193301 0.000000 +1549 0.294480 0.194119 0.000000 +1550 0.265376 0.192145 0.000000 +1551 0.233926 0.191285 0.000000 +1552 0.201869 0.196460 0.000000 +1553 0.172943 0.200159 0.000000 +1554 0.143553 0.200482 0.000000 +1555 0.113782 0.199148 0.000000 +1556 0.083231 0.198393 0.000000 +1557 0.052135 0.199830 0.000000 +1558 0.021092 0.201266 0.000000 +1559 -0.009922 0.201866 0.000000 +1560 -0.040700 0.201511 0.000000 +1561 -0.071228 0.200508 0.000000 +1562 -0.101438 0.199427 0.000000 +1563 -0.131223 0.198305 0.000000 +1564 -0.159535 0.194395 0.000000 +1566 0.578134 -0.204375 0.000000 +1567 0.608049 -0.203173 0.000000 +1568 0.637979 -0.201995 0.000000 +1569 0.667952 -0.200914 0.000000 +1570 0.698007 -0.199916 0.000000 +1571 0.697532 -0.169905 0.000000 +1572 0.697003 -0.139791 0.000000 +1573 0.696380 -0.109617 0.000000 +1574 0.695597 -0.079155 0.000000 +1575 0.694688 -0.048469 0.000000 +1576 0.693667 -0.017488 0.000000 +1577 0.692727 0.013854 0.000000 +1578 0.692318 0.045364 0.000000 +1579 0.693432 0.076853 0.000000 +1580 0.696657 0.107710 0.000000 +1581 0.699224 0.138222 0.000000 +1582 0.700642 0.168450 0.000000 +1583 0.371214 0.017358 0.000000 +1584 0.401698 0.024012 0.000000 +1585 0.432909 0.027812 0.000000 +1586 0.464397 0.030558 0.000000 +1587 0.495928 0.033147 0.000000 +1588 0.527675 0.035870 0.000000 +1590 0.564245 0.007571 0.000000 +1589 0.559782 0.038972 0.000000 +1591 0.567852 -0.023731 0.000000 +1592 0.570714 -0.054539 0.000000 +1593 0.572988 -0.084893 0.000000 +1594 0.574798 -0.114991 0.000000 +1595 0.576239 -0.144944 0.000000 +1596 0.577286 -0.174743 0.000000 +1597 0.672269 0.169162 0.000000 +1598 0.644670 0.169920 0.000000 +1599 0.617885 0.169958 0.000000 +1600 0.591267 0.167167 0.000000 +1601 0.563698 0.164256 0.000000 +1602 0.535150 0.161558 0.000000 +1603 0.505769 0.159594 0.000000 +1604 0.475777 0.158082 0.000000 +1605 0.445365 0.156971 0.000000 +1606 0.414493 0.156592 0.000000 +1607 0.383621 0.157125 0.000000 +1608 0.353035 0.160313 0.000000 +1609 0.324219 0.164732 0.000000 +1610 0.298724 0.166543 0.000000 +1611 0.271789 0.160670 0.000000 +1651 0.199126 0.167293 0.000000 +1612 0.173017 0.174362 0.000000 +1613 0.147220 0.172881 0.000000 +1614 0.118370 0.169386 0.000000 +1615 0.086388 0.167114 0.000000 +1616 0.053408 0.170284 0.000000 +1617 0.022505 0.172511 0.000000 +1618 -0.008623 0.173287 0.000000 +1619 -0.039412 0.172568 0.000000 +1620 -0.069931 0.170982 0.000000 +1621 -0.100028 0.170519 0.000000 +1622 -0.128149 0.170934 0.000000 +1623 -0.152764 0.169059 0.000000 +1625 0.607304 -0.173381 0.000000 +1626 0.637335 -0.172134 0.000000 +1627 0.667402 -0.170981 0.000000 +1628 0.666783 -0.140930 0.000000 +1629 0.665958 -0.110595 0.000000 +1630 0.664919 -0.080023 0.000000 +1631 0.663654 -0.049157 0.000000 +1632 0.662211 -0.017858 0.000000 +1633 0.660734 0.013978 0.000000 +1634 0.659723 0.046364 0.000000 +1635 0.661116 0.078961 0.000000 +1636 0.667246 0.110337 0.000000 +1637 0.670745 0.139548 0.000000 +1638 -0.147799 0.150716 0.000000 +1674 -0.128684 0.145691 0.000000 +1673 -0.100884 0.140613 0.000000 +1640 -0.044570 0.054253 0.000000 +1639 -0.056446 0.080725 0.000000 +1641 -0.032234 0.029928 0.000000 +1642 -0.006024 0.037236 0.000000 +1643 0.020359 0.044478 0.000000 +1644 0.045874 0.052989 0.000000 +1645 0.069878 0.064328 0.000000 +1646 0.092192 0.079845 0.000000 +1647 0.116360 0.097485 0.000000 +1648 0.143869 0.111631 0.000000 +1649 0.168436 0.126174 0.000000 +1650 0.188283 0.142853 0.000000 +1652 0.358550 0.049279 0.000000 +1653 0.392847 0.055605 0.000000 +1654 0.425498 0.059383 0.000000 +1655 0.457980 0.061988 0.000000 +1656 0.489964 0.064583 0.000000 +1657 0.521756 0.067277 0.000000 +1658 0.553712 0.070546 0.000000 +1660 0.592058 0.043198 0.000000 +1659 0.586126 0.075499 0.000000 +1661 0.596024 0.010863 0.000000 +1662 0.599068 -0.020945 0.000000 +1663 0.601524 -0.052102 0.000000 +1664 0.603505 -0.082800 0.000000 +1665 0.605107 -0.113223 0.000000 +1666 0.606324 -0.143417 0.000000 +1667 0.171493 0.155016 0.000000 +1696 0.153180 0.147834 0.000000 +1695 0.126680 0.139244 0.000000 +1694 0.090766 0.129758 0.000000 +1668 0.053205 0.140969 0.000000 +1669 0.022718 0.144596 0.000000 +1670 -0.007526 0.145267 0.000000 +1671 -0.037468 0.143511 0.000000 +1672 -0.068877 0.141356 0.000000 +1675 0.643744 0.141837 0.000000 +1676 0.619340 0.143192 0.000000 +1677 0.595381 0.139507 0.000000 +1678 0.568637 0.134721 0.000000 +1679 0.539884 0.131489 0.000000 +1680 0.510058 0.128996 0.000000 +1681 0.479577 0.126939 0.000000 +1682 0.448267 0.125455 0.000000 +1683 0.416247 0.124554 0.000000 +1684 0.383738 0.125001 0.000000 +1685 0.350239 0.128717 0.000000 +1686 0.321325 0.138246 0.000000 +1687 0.301817 0.146183 0.000000 +1688 0.285401 0.133759 0.000000 +1689 0.636514 -0.142039 0.000000 +1719 0.635429 -0.111764 0.000000 +1718 0.634140 -0.081198 0.000000 +1717 0.632516 -0.050293 0.000000 +1716 0.630582 -0.018910 0.000000 +1715 0.628349 0.013302 0.000000 +1714 0.625872 0.047113 0.000000 +1713 0.623071 0.085699 0.000000 +1690 0.639907 0.114926 0.000000 +1707 0.307163 0.113588 0.000000 +1706 0.342214 0.089227 0.000000 +1711 -0.017903 0.061881 0.000000 +1710 0.008821 0.068462 0.000000 +1691 0.033317 0.075470 0.000000 +1692 0.053991 0.085182 0.000000 +1693 0.071183 0.102698 0.000000 +1697 0.620610 0.123584 0.000000 +1698 0.602483 0.113360 0.000000 +1699 0.575576 0.106296 0.000000 +1700 0.546046 0.101572 0.000000 +1701 0.515442 0.098403 0.000000 +1702 0.484200 0.095912 0.000000 +1703 0.452238 0.093868 0.000000 +1704 0.419736 0.091956 0.000000 +1705 0.384875 0.090284 0.000000 +1708 0.047562 0.116186 0.000000 +1709 0.022481 0.118825 0.000000 +1727 -0.004797 0.117886 0.000000 +1729 -0.034224 0.114874 0.000000 +1712 -0.065866 0.110365 0.000000 +1721 0.041866 0.099182 0.000000 +1720 0.024047 0.096178 0.000000 +1722 0.000778 -0.304198 0.000000 +1724 0.015932 -0.322856 0.000000 +1726 -0.000448 0.092548 0.000000 +1728 -0.027125 0.087460 0.000000 +1426 0.323105 -0.054530 0.000000 +1 1 quad 195 133 67 180 +2 1 quad 196 132 133 195 +3 1 quad 196 195 180 179 +4 1 quad 132 196 197 131 +5 1 quad 196 179 178 197 +6 1 quad 181 187 198 185 +7 1 quad 187 188 199 198 +8 1 quad 188 189 200 199 +9 1 quad 189 190 201 200 +10 1 quad 190 191 202 201 +11 1 quad 191 192 203 202 +12 1 quad 192 186 204 203 +13 1 quad 186 206 205 204 +14 1 quad 186 193 207 206 +15 1 quad 193 194 208 207 +16 1 quad 194 1 209 208 +17 1 quad 1 211 210 209 +18 1 quad 1 3 212 211 +19 1 quad 3 4 213 212 +20 1 quad 4 5 214 213 +21 1 quad 5 6 215 214 +22 1 quad 6 7 216 215 +23 1 quad 7 8 217 216 +24 1 quad 8 2 10 217 +25 1 quad 34 68 218 66 +26 1 quad 68 69 219 218 +27 1 quad 69 70 220 219 +28 1 quad 70 71 221 220 +29 1 quad 71 72 222 221 +30 1 quad 72 73 223 222 +31 1 quad 73 74 224 223 +32 1 quad 74 75 225 224 +33 1 quad 75 76 226 225 +34 1 quad 76 77 227 226 +35 1 quad 77 78 228 227 +36 1 quad 78 79 229 228 +37 1 quad 79 80 230 229 +38 1 quad 80 81 231 230 +39 1 quad 81 82 232 231 +40 1 quad 82 83 233 232 +41 1 quad 83 84 234 233 +42 1 quad 84 85 235 234 +43 1 quad 85 86 236 235 +44 1 quad 86 87 237 236 +45 1 quad 87 88 238 237 +46 1 quad 88 89 239 238 +47 1 quad 89 90 240 239 +48 1 quad 90 91 241 240 +49 1 quad 91 92 242 241 +50 1 quad 92 93 243 242 +51 1 quad 93 94 244 243 +52 1 quad 94 95 245 244 +53 1 quad 95 96 246 245 +54 1 quad 96 97 247 246 +55 1 quad 97 98 248 247 +56 1 quad 98 99 249 248 +57 1 quad 99 100 250 249 +58 1 quad 100 101 251 250 +59 1 quad 101 102 252 251 +60 1 quad 102 103 253 252 +61 1 quad 103 104 254 253 +62 1 quad 104 105 255 254 +63 1 quad 105 106 256 255 +64 1 quad 106 107 257 256 +65 1 quad 107 108 258 257 +66 1 quad 108 109 259 258 +67 1 quad 109 110 260 259 +68 1 quad 110 111 261 260 +69 1 quad 111 112 262 261 +70 1 quad 112 113 263 262 +71 1 quad 113 114 264 263 +72 1 quad 114 115 265 264 +73 1 quad 115 116 266 265 +74 1 quad 116 117 267 266 +75 1 quad 117 118 268 267 +76 1 quad 118 119 269 268 +77 1 quad 119 120 270 269 +78 1 quad 120 121 271 270 +79 1 quad 121 122 272 271 +80 1 quad 122 123 273 272 +81 1 quad 123 124 274 273 +82 1 quad 124 125 275 274 +83 1 quad 125 126 276 275 +84 1 quad 126 127 277 276 +85 1 quad 127 128 278 277 +86 1 quad 128 129 279 278 +87 1 quad 129 130 303 279 +88 1 quad 130 131 197 303 +89 1 quad 10 11 280 217 +90 1 quad 11 12 281 280 +91 1 quad 12 13 282 281 +92 1 quad 13 14 283 282 +93 1 quad 14 15 284 283 +94 1 quad 15 16 285 284 +95 1 quad 16 17 286 285 +96 1 quad 17 18 287 286 +97 1 quad 18 19 288 287 +98 1 quad 19 20 289 288 +99 1 quad 20 21 290 289 +100 1 quad 21 22 291 290 +101 1 quad 22 23 292 291 +102 1 quad 23 24 293 292 +103 1 quad 24 25 294 293 +104 1 quad 25 26 295 294 +105 1 quad 26 27 296 295 +106 1 quad 27 28 297 296 +107 1 quad 28 29 298 297 +108 1 quad 29 30 299 298 +109 1 quad 30 31 300 299 +110 1 quad 31 32 301 300 +111 1 quad 32 33 302 301 +112 1 quad 33 9 35 302 +113 1 quad 178 177 303 197 +114 1 quad 177 176 304 303 +115 1 quad 176 175 305 304 +116 1 quad 175 174 306 305 +117 1 quad 174 173 307 306 +118 1 quad 173 172 308 307 +119 1 quad 172 171 309 308 +120 1 quad 171 170 310 309 +121 1 quad 170 169 311 310 +122 1 quad 169 168 312 311 +123 1 quad 168 167 313 312 +124 1 quad 167 166 314 313 +125 1 quad 166 165 315 314 +126 1 quad 165 164 316 315 +127 1 quad 164 163 317 316 +128 1 quad 163 162 318 317 +129 1 quad 162 161 319 318 +130 1 quad 161 160 320 319 +131 1 quad 160 159 321 320 +132 1 quad 159 158 322 321 +133 1 quad 158 157 323 322 +134 1 quad 157 156 324 323 +135 1 quad 156 155 325 324 +136 1 quad 155 154 326 325 +137 1 quad 154 153 327 326 +138 1 quad 153 152 328 327 +139 1 quad 152 151 329 328 +140 1 quad 151 150 330 329 +141 1 quad 150 149 331 330 +142 1 quad 149 148 332 331 +143 1 quad 148 147 333 332 +144 1 quad 147 146 334 333 +145 1 quad 146 145 335 334 +146 1 quad 145 144 336 335 +147 1 quad 144 143 337 336 +148 1 quad 143 142 338 337 +149 1 quad 142 141 339 338 +150 1 quad 141 140 340 339 +151 1 quad 140 139 341 340 +152 1 quad 139 138 342 341 +153 1 quad 138 137 343 342 +154 1 quad 137 136 344 343 +155 1 quad 136 135 345 344 +156 1 quad 135 134 346 345 +157 1 quad 134 182 347 346 +158 1 quad 182 183 348 347 +159 1 quad 183 184 349 348 +160 1 quad 184 185 198 349 +161 1 quad 35 36 350 302 +162 1 quad 36 37 351 350 +163 1 quad 37 38 352 351 +164 1 quad 38 39 353 352 +165 1 quad 39 40 354 353 +166 1 quad 40 41 355 354 +167 1 quad 41 42 356 355 +168 1 quad 42 43 357 356 +169 1 quad 43 44 358 357 +170 1 quad 44 45 359 358 +171 1 quad 45 46 360 359 +172 1 quad 46 47 361 360 +173 1 quad 47 48 362 361 +174 1 quad 48 49 363 362 +175 1 quad 49 50 364 363 +176 1 quad 50 51 365 364 +177 1 quad 51 52 366 365 +178 1 quad 52 53 367 366 +179 1 quad 53 54 368 367 +180 1 quad 54 55 369 368 +181 1 quad 55 56 370 369 +182 1 quad 56 57 371 370 +183 1 quad 57 58 372 371 +184 1 quad 58 59 373 372 +185 1 quad 59 60 374 373 +186 1 quad 60 61 375 374 +187 1 quad 61 62 376 375 +188 1 quad 62 63 377 376 +189 1 quad 63 64 378 377 +190 1 quad 64 65 379 378 +191 1 quad 65 66 218 379 +192 1 quad 303 304 503 279 +193 1 quad 304 305 380 503 +194 1 quad 305 306 381 380 +195 1 quad 306 307 382 381 +196 1 quad 307 308 383 382 +197 1 quad 308 309 384 383 +198 1 quad 309 310 385 384 +199 1 quad 310 311 386 385 +200 1 quad 311 312 387 386 +201 1 quad 312 313 388 387 +202 1 quad 313 314 389 388 +203 1 quad 314 315 390 389 +204 1 quad 315 316 391 390 +205 1 quad 316 317 392 391 +206 1 quad 317 318 393 392 +207 1 quad 318 319 394 393 +208 1 quad 319 320 395 394 +209 1 quad 320 321 396 395 +210 1 quad 321 322 397 396 +211 1 quad 322 323 398 397 +212 1 quad 323 324 399 398 +213 1 quad 324 325 400 399 +214 1 quad 325 326 401 400 +215 1 quad 326 327 402 401 +216 1 quad 327 328 403 402 +217 1 quad 328 329 404 403 +218 1 quad 329 330 405 404 +219 1 quad 330 331 406 405 +220 1 quad 331 332 407 406 +221 1 quad 332 333 408 407 +222 1 quad 333 334 409 408 +223 1 quad 334 335 410 409 +224 1 quad 335 336 411 410 +225 1 quad 336 337 412 411 +226 1 quad 337 338 413 412 +227 1 quad 338 339 414 413 +228 1 quad 339 340 415 414 +229 1 quad 340 341 416 415 +230 1 quad 341 342 417 416 +231 1 quad 342 343 418 417 +232 1 quad 343 344 419 418 +233 1 quad 347 348 421 420 +234 1 quad 348 349 422 421 +235 1 quad 349 198 199 422 +236 1 quad 199 200 423 422 +237 1 quad 200 201 424 423 +238 1 quad 201 202 425 424 +239 1 quad 202 203 426 425 +240 1 quad 203 204 427 426 +241 1 quad 204 205 428 427 +242 1 quad 205 430 429 428 +243 1 quad 205 206 431 430 +244 1 quad 206 207 432 431 +245 1 quad 207 208 433 432 +246 1 quad 208 209 434 433 +247 1 quad 209 210 435 434 +248 1 quad 210 437 436 435 +249 1 quad 210 211 438 437 +250 1 quad 211 212 439 438 +251 1 quad 212 213 440 439 +252 1 quad 213 214 441 440 +253 1 quad 214 215 442 441 +254 1 quad 215 216 443 442 +255 1 quad 216 217 280 443 +256 1 quad 218 219 444 379 +257 1 quad 219 220 445 444 +258 1 quad 220 221 446 445 +259 1 quad 221 222 447 446 +260 1 quad 222 223 448 447 +261 1 quad 223 224 449 448 +262 1 quad 224 225 450 449 +263 1 quad 225 226 451 450 +264 1 quad 226 227 452 451 +265 1 quad 227 228 453 452 +266 1 quad 228 229 454 453 +267 1 quad 229 230 455 454 +268 1 quad 230 231 456 455 +269 1 quad 231 232 457 456 +270 1 quad 232 233 458 457 +271 1 quad 233 234 459 458 +272 1 quad 234 235 460 459 +273 1 quad 235 236 461 460 +274 1 quad 236 237 462 461 +275 1 quad 237 238 463 462 +276 1 quad 238 239 464 463 +277 1 quad 239 240 465 464 +278 1 quad 240 241 466 465 +279 1 quad 241 242 467 466 +280 1 quad 242 243 468 467 +281 1 quad 243 244 469 468 +282 1 quad 244 245 470 469 +283 1 quad 245 246 471 470 +284 1 quad 246 247 472 471 +285 1 quad 247 248 473 472 +286 1 quad 248 249 474 473 +287 1 quad 249 250 475 474 +288 1 quad 250 251 476 475 +289 1 quad 251 252 477 476 +290 1 quad 252 253 478 477 +291 1 quad 253 254 479 478 +292 1 quad 254 255 480 479 +293 1 quad 255 256 481 480 +294 1 quad 256 257 482 481 +295 1 quad 257 258 483 482 +296 1 quad 258 259 484 483 +297 1 quad 259 260 485 484 +298 1 quad 260 261 486 485 +299 1 quad 261 262 487 486 +300 1 quad 262 263 488 487 +301 1 quad 263 264 489 488 +302 1 quad 264 265 490 489 +303 1 quad 265 266 491 490 +304 1 quad 266 267 492 491 +305 1 quad 267 268 493 492 +306 1 quad 268 269 494 493 +307 1 quad 269 270 495 494 +308 1 quad 270 271 496 495 +309 1 quad 271 272 497 496 +310 1 quad 272 273 498 497 +311 1 quad 273 274 499 498 +312 1 quad 274 275 500 499 +313 1 quad 275 276 501 500 +314 1 quad 276 277 553 501 +315 1 quad 277 278 502 553 +316 1 quad 278 279 503 502 +317 1 quad 302 350 504 301 +318 1 quad 350 351 505 504 +319 1 quad 351 352 506 505 +320 1 quad 352 353 507 506 +321 1 quad 353 354 508 507 +322 1 quad 354 355 509 508 +323 1 quad 355 356 510 509 +324 1 quad 356 357 511 510 +325 1 quad 357 358 512 511 +326 1 quad 358 359 513 512 +327 1 quad 359 360 514 513 +328 1 quad 360 361 515 514 +329 1 quad 361 362 516 515 +330 1 quad 362 363 517 516 +331 1 quad 363 364 518 517 +332 1 quad 364 365 519 518 +333 1 quad 365 366 520 519 +334 1 quad 366 367 521 520 +335 1 quad 367 368 522 521 +336 1 quad 368 369 523 522 +337 1 quad 369 370 524 523 +338 1 quad 370 371 525 524 +339 1 quad 371 372 526 525 +340 1 quad 372 373 527 526 +341 1 quad 373 374 528 527 +342 1 quad 374 375 529 528 +343 1 quad 375 376 530 529 +344 1 quad 376 377 531 530 +345 1 quad 377 378 532 531 +346 1 quad 378 379 444 532 +347 1 quad 280 281 533 443 +348 1 quad 281 282 534 533 +349 1 quad 282 283 535 534 +350 1 quad 283 284 536 535 +351 1 quad 284 285 537 536 +352 1 quad 285 286 538 537 +353 1 quad 286 287 539 538 +354 1 quad 287 288 540 539 +355 1 quad 288 289 541 540 +356 1 quad 289 290 542 541 +357 1 quad 290 291 543 542 +358 1 quad 291 292 544 543 +359 1 quad 292 293 545 544 +360 1 quad 293 294 546 545 +361 1 quad 294 295 547 546 +362 1 quad 295 296 548 547 +363 1 quad 296 297 549 548 +364 1 quad 297 298 550 549 +365 1 quad 298 299 551 550 +366 1 quad 299 300 552 551 +367 1 quad 300 301 504 552 +368 1 quad 420 421 590 419 +369 1 quad 421 422 423 590 +370 1 quad 503 380 553 502 +371 1 quad 380 381 501 553 +372 1 quad 381 382 670 501 +373 1 quad 382 383 554 670 +374 1 quad 383 384 555 554 +375 1 quad 384 385 556 555 +376 1 quad 385 386 557 556 +377 1 quad 386 387 558 557 +378 1 quad 387 388 559 558 +379 1 quad 388 389 560 559 +380 1 quad 389 390 561 560 +381 1 quad 390 391 562 561 +382 1 quad 391 392 563 562 +383 1 quad 392 393 564 563 +384 1 quad 393 394 565 564 +385 1 quad 394 395 566 565 +386 1 quad 395 396 567 566 +387 1 quad 396 397 568 567 +388 1 quad 397 398 569 568 +389 1 quad 398 399 570 569 +390 1 quad 399 400 571 570 +391 1 quad 400 401 572 571 +392 1 quad 401 402 573 572 +393 1 quad 402 403 574 573 +394 1 quad 403 404 575 574 +395 1 quad 404 405 576 575 +396 1 quad 405 406 577 576 +397 1 quad 406 407 578 577 +398 1 quad 407 408 579 578 +399 1 quad 408 409 580 579 +400 1 quad 409 410 581 580 +401 1 quad 410 411 582 581 +402 1 quad 411 412 583 582 +403 1 quad 412 413 584 583 +404 1 quad 413 414 585 584 +405 1 quad 414 415 586 585 +406 1 quad 415 416 587 586 +407 1 quad 416 417 588 587 +408 1 quad 417 418 589 588 +409 1 quad 418 419 590 589 +410 1 quad 423 424 591 590 +411 1 quad 424 425 592 591 +412 1 quad 425 426 593 592 +413 1 quad 426 427 594 593 +414 1 quad 427 428 595 594 +415 1 quad 428 429 596 595 +416 1 quad 429 598 597 596 +417 1 quad 429 430 599 598 +418 1 quad 430 431 600 599 +419 1 quad 431 432 601 600 +420 1 quad 432 433 602 601 +421 1 quad 433 434 603 602 +422 1 quad 434 435 604 603 +423 1 quad 435 436 605 604 +424 1 quad 436 607 606 605 +425 1 quad 436 437 608 607 +426 1 quad 437 438 609 608 +427 1 quad 438 439 610 609 +428 1 quad 439 440 611 610 +429 1 quad 440 441 612 611 +430 1 quad 441 442 613 612 +431 1 quad 442 443 533 613 +432 1 quad 444 445 614 532 +433 1 quad 445 446 615 614 +434 1 quad 446 447 616 615 +435 1 quad 447 448 617 616 +436 1 quad 448 449 618 617 +437 1 quad 449 450 619 618 +438 1 quad 450 451 620 619 +439 1 quad 451 452 621 620 +440 1 quad 452 453 622 621 +441 1 quad 453 454 623 622 +442 1 quad 454 455 624 623 +443 1 quad 455 456 625 624 +444 1 quad 456 457 626 625 +445 1 quad 457 458 627 626 +446 1 quad 458 459 628 627 +447 1 quad 459 460 629 628 +448 1 quad 460 461 630 629 +449 1 quad 461 462 631 630 +450 1 quad 462 463 632 631 +451 1 quad 463 464 633 632 +452 1 quad 464 465 634 633 +453 1 quad 465 466 635 634 +454 1 quad 466 467 636 635 +455 1 quad 467 468 637 636 +456 1 quad 468 469 638 637 +457 1 quad 469 470 639 638 +458 1 quad 470 471 640 639 +459 1 quad 471 472 641 640 +460 1 quad 472 473 642 641 +461 1 quad 473 474 643 642 +462 1 quad 474 475 644 643 +463 1 quad 475 476 645 644 +464 1 quad 476 477 646 645 +465 1 quad 477 478 647 646 +466 1 quad 478 479 648 647 +467 1 quad 479 480 649 648 +468 1 quad 480 481 650 649 +469 1 quad 481 482 651 650 +470 1 quad 482 483 652 651 +471 1 quad 483 484 653 652 +472 1 quad 484 485 654 653 +473 1 quad 485 486 655 654 +474 1 quad 486 487 656 655 +475 1 quad 487 488 657 656 +476 1 quad 488 489 658 657 +477 1 quad 489 490 659 658 +478 1 quad 490 491 660 659 +479 1 quad 491 492 661 660 +480 1 quad 492 493 662 661 +481 1 quad 493 494 663 662 +482 1 quad 494 495 664 663 +483 1 quad 495 496 665 664 +484 1 quad 496 497 666 665 +485 1 quad 497 498 667 666 +486 1 quad 498 499 668 667 +487 1 quad 499 500 669 668 +488 1 quad 500 501 670 669 +489 1 quad 533 534 671 613 +490 1 quad 534 535 672 671 +491 1 quad 535 536 673 672 +492 1 quad 536 537 674 673 +493 1 quad 537 538 675 674 +494 1 quad 538 539 676 675 +495 1 quad 539 540 677 676 +496 1 quad 540 541 678 677 +497 1 quad 541 542 679 678 +498 1 quad 542 543 680 679 +499 1 quad 543 544 681 680 +500 1 quad 544 545 682 681 +501 1 quad 545 546 683 682 +502 1 quad 546 547 684 683 +503 1 quad 547 548 685 684 +504 1 quad 548 549 686 685 +505 1 quad 549 550 687 686 +506 1 quad 550 551 688 687 +507 1 quad 551 552 689 688 +508 1 quad 552 504 505 689 +509 1 quad 505 506 690 689 +510 1 quad 506 507 691 690 +511 1 quad 507 508 692 691 +512 1 quad 508 509 693 692 +513 1 quad 509 510 694 693 +514 1 quad 510 511 695 694 +515 1 quad 511 512 696 695 +516 1 quad 512 513 697 696 +517 1 quad 513 514 698 697 +518 1 quad 514 515 699 698 +519 1 quad 515 516 700 699 +520 1 quad 516 517 701 700 +521 1 quad 517 518 702 701 +522 1 quad 518 519 703 702 +523 1 quad 519 520 704 703 +524 1 quad 520 521 705 704 +525 1 quad 521 522 706 705 +526 1 quad 522 523 707 706 +527 1 quad 523 524 708 707 +528 1 quad 524 525 709 708 +529 1 quad 525 526 710 709 +530 1 quad 526 527 711 710 +531 1 quad 527 528 712 711 +532 1 quad 528 529 713 712 +533 1 quad 529 530 714 713 +534 1 quad 530 531 715 714 +535 1 quad 531 532 614 715 +536 1 quad 590 591 716 589 +537 1 quad 591 592 717 716 +538 1 quad 592 593 718 717 +539 1 quad 593 594 719 718 +540 1 quad 594 595 720 719 +541 1 quad 595 596 721 720 +542 1 quad 596 597 722 721 +543 1 quad 597 724 723 722 +544 1 quad 597 598 725 724 +545 1 quad 598 599 726 725 +546 1 quad 599 600 727 726 +547 1 quad 600 601 728 727 +548 1 quad 601 602 729 728 +549 1 quad 602 603 730 729 +550 1 quad 603 604 731 730 +551 1 quad 604 605 732 731 +552 1 quad 605 606 733 732 +553 1 quad 606 735 734 733 +554 1 quad 606 607 736 735 +555 1 quad 607 608 737 736 +556 1 quad 608 609 738 737 +557 1 quad 609 610 739 738 +558 1 quad 610 611 740 739 +559 1 quad 611 612 741 740 +560 1 quad 612 613 671 741 +561 1 quad 614 615 742 715 +562 1 quad 615 616 743 742 +563 1 quad 616 617 744 743 +564 1 quad 617 618 745 744 +565 1 quad 618 619 746 745 +566 1 quad 619 620 747 746 +567 1 quad 620 621 748 747 +568 1 quad 621 622 749 748 +569 1 quad 622 623 750 749 +570 1 quad 623 624 751 750 +571 1 quad 624 625 752 751 +572 1 quad 625 626 753 752 +573 1 quad 626 627 754 753 +574 1 quad 627 628 755 754 +575 1 quad 628 629 756 755 +576 1 quad 629 630 757 756 +577 1 quad 630 631 758 757 +578 1 quad 631 632 759 758 +579 1 quad 632 633 760 759 +580 1 quad 633 634 761 760 +581 1 quad 634 635 762 761 +582 1 quad 635 636 763 762 +583 1 quad 636 637 764 763 +584 1 quad 637 638 765 764 +585 1 quad 638 639 766 765 +586 1 quad 639 640 767 766 +587 1 quad 640 641 768 767 +588 1 quad 641 642 769 768 +589 1 quad 642 643 770 769 +590 1 quad 643 644 771 770 +591 1 quad 644 645 772 771 +592 1 quad 645 646 773 772 +593 1 quad 646 647 774 773 +594 1 quad 647 648 775 774 +595 1 quad 648 649 776 775 +596 1 quad 649 650 777 776 +597 1 quad 650 651 778 777 +598 1 quad 651 652 779 778 +599 1 quad 652 653 780 779 +600 1 quad 653 654 781 780 +601 1 quad 654 655 782 781 +602 1 quad 655 656 783 782 +603 1 quad 656 657 784 783 +604 1 quad 657 658 785 784 +605 1 quad 658 659 786 785 +606 1 quad 659 660 787 786 +607 1 quad 660 661 788 787 +608 1 quad 661 662 789 788 +609 1 quad 662 663 790 789 +610 1 quad 663 664 791 790 +611 1 quad 664 665 792 791 +612 1 quad 665 666 793 792 +613 1 quad 666 667 794 793 +614 1 quad 667 668 795 794 +615 1 quad 668 669 796 795 +616 1 quad 669 670 554 796 +617 1 quad 671 672 797 741 +618 1 quad 672 673 798 797 +619 1 quad 673 674 799 798 +620 1 quad 674 675 800 799 +621 1 quad 675 676 801 800 +622 1 quad 676 677 802 801 +623 1 quad 677 678 803 802 +624 1 quad 678 679 804 803 +625 1 quad 679 680 805 804 +626 1 quad 680 681 806 805 +627 1 quad 681 682 807 806 +628 1 quad 682 683 808 807 +629 1 quad 683 684 809 808 +630 1 quad 684 685 810 809 +631 1 quad 685 686 811 810 +632 1 quad 686 687 812 811 +633 1 quad 687 688 813 812 +634 1 quad 688 689 690 813 +635 1 quad 554 555 814 796 +636 1 quad 555 556 815 814 +637 1 quad 556 557 816 815 +638 1 quad 557 558 817 816 +639 1 quad 558 559 818 817 +640 1 quad 559 560 819 818 +641 1 quad 560 561 820 819 +642 1 quad 561 562 821 820 +643 1 quad 562 563 822 821 +644 1 quad 563 564 823 822 +645 1 quad 564 565 824 823 +646 1 quad 565 566 825 824 +647 1 quad 566 567 826 825 +648 1 quad 567 568 827 826 +649 1 quad 568 569 828 827 +650 1 quad 569 570 829 828 +651 1 quad 570 571 830 829 +652 1 quad 571 572 831 830 +653 1 quad 572 573 832 831 +654 1 quad 573 574 833 832 +655 1 quad 574 575 834 833 +656 1 quad 575 576 835 834 +657 1 quad 576 577 836 835 +658 1 quad 577 578 837 836 +659 1 quad 578 579 838 837 +660 1 quad 579 580 839 838 +661 1 quad 580 581 840 839 +662 1 quad 581 582 841 840 +663 1 quad 582 583 842 841 +664 1 quad 583 584 843 842 +665 1 quad 584 585 844 843 +666 1 quad 585 586 845 844 +667 1 quad 586 587 1723 845 +668 1 quad 587 588 1725 1723 +669 1 quad 588 589 716 1725 +670 1 quad 690 691 848 813 +671 1 quad 691 692 849 848 +672 1 quad 692 693 850 849 +673 1 quad 693 694 851 850 +674 1 quad 694 695 852 851 +675 1 quad 695 696 853 852 +676 1 quad 696 697 854 853 +677 1 quad 697 698 855 854 +678 1 quad 698 699 856 855 +679 1 quad 699 700 857 856 +680 1 quad 700 701 858 857 +681 1 quad 701 702 859 858 +682 1 quad 702 703 860 859 +683 1 quad 703 704 861 860 +684 1 quad 704 705 862 861 +685 1 quad 705 706 863 862 +686 1 quad 706 707 864 863 +687 1 quad 707 708 865 864 +688 1 quad 708 709 866 865 +689 1 quad 709 710 867 866 +690 1 quad 710 711 868 867 +691 1 quad 711 712 869 868 +692 1 quad 712 713 870 869 +693 1 quad 713 714 871 870 +694 1 quad 714 715 742 871 +695 1 quad 718 719 846 847 +696 1 quad 719 720 872 846 +697 1 quad 720 721 873 872 +698 1 quad 721 722 874 873 +699 1 quad 722 723 875 874 +700 1 quad 723 877 876 875 +701 1 quad 723 724 878 877 +702 1 quad 724 725 879 878 +703 1 quad 725 726 880 879 +704 1 quad 726 727 881 880 +705 1 quad 727 728 882 881 +706 1 quad 728 729 883 882 +707 1 quad 729 730 884 883 +708 1 quad 730 731 885 884 +709 1 quad 731 732 886 885 +710 1 quad 732 733 887 886 +711 1 quad 733 734 888 887 +712 1 quad 734 890 889 888 +713 1 quad 734 735 891 890 +714 1 quad 735 736 892 891 +715 1 quad 736 737 893 892 +716 1 quad 737 738 894 893 +717 1 quad 738 739 895 894 +718 1 quad 739 740 896 895 +719 1 quad 740 741 797 896 +720 1 quad 796 814 897 795 +721 1 quad 814 815 898 897 +722 1 quad 815 816 1040 898 +723 1 quad 816 817 976 1040 +724 1 quad 817 818 899 976 +725 1 quad 818 819 900 899 +726 1 quad 819 820 901 900 +727 1 quad 820 821 902 901 +728 1 quad 821 822 903 902 +729 1 quad 822 823 904 903 +730 1 quad 823 824 905 904 +731 1 quad 824 825 906 905 +732 1 quad 825 826 907 906 +733 1 quad 826 827 908 907 +734 1 quad 827 828 909 908 +735 1 quad 828 829 910 909 +736 1 quad 829 830 911 910 +737 1 quad 830 831 912 911 +738 1 quad 831 832 913 912 +739 1 quad 832 833 914 913 +740 1 quad 833 834 915 914 +741 1 quad 834 835 916 915 +742 1 quad 835 836 917 916 +743 1 quad 836 837 918 917 +744 1 quad 837 838 919 918 +745 1 quad 838 839 920 919 +746 1 quad 839 840 921 920 +747 1 quad 840 841 922 921 +748 1 quad 841 842 923 922 +749 1 quad 842 843 924 923 +750 1 quad 843 844 925 924 +751 1 quad 844 845 926 925 +752 1 quad 742 743 927 871 +753 1 quad 743 744 928 927 +754 1 quad 744 745 929 928 +755 1 quad 745 746 930 929 +756 1 quad 746 747 931 930 +757 1 quad 747 748 932 931 +758 1 quad 748 749 933 932 +759 1 quad 749 750 934 933 +760 1 quad 750 751 935 934 +761 1 quad 751 752 936 935 +762 1 quad 752 753 937 936 +763 1 quad 753 754 938 937 +764 1 quad 754 755 939 938 +765 1 quad 755 756 940 939 +766 1 quad 756 757 941 940 +767 1 quad 757 758 942 941 +768 1 quad 758 759 943 942 +769 1 quad 759 760 944 943 +770 1 quad 760 761 945 944 +771 1 quad 761 762 946 945 +772 1 quad 762 763 947 946 +773 1 quad 763 764 948 947 +774 1 quad 764 765 949 948 +775 1 quad 765 766 950 949 +776 1 quad 766 767 951 950 +777 1 quad 767 768 952 951 +778 1 quad 768 769 953 952 +779 1 quad 769 770 954 953 +780 1 quad 770 771 955 954 +781 1 quad 771 772 956 955 +782 1 quad 772 773 957 956 +783 1 quad 773 774 958 957 +784 1 quad 774 775 959 958 +785 1 quad 775 776 960 959 +786 1 quad 776 777 961 960 +787 1 quad 777 778 962 961 +788 1 quad 778 779 963 962 +789 1 quad 779 780 964 963 +790 1 quad 780 781 965 964 +791 1 quad 781 782 966 965 +792 1 quad 782 783 967 966 +793 1 quad 783 784 968 967 +794 1 quad 784 785 969 968 +795 1 quad 785 786 970 969 +796 1 quad 786 787 971 970 +797 1 quad 787 788 972 971 +798 1 quad 788 789 973 972 +799 1 quad 789 790 974 973 +800 1 quad 790 791 975 974 +801 1 quad 791 792 976 975 +802 1 quad 797 798 977 896 +803 1 quad 798 799 978 977 +804 1 quad 799 800 979 978 +805 1 quad 800 801 980 979 +806 1 quad 801 802 981 980 +807 1 quad 802 803 982 981 +808 1 quad 803 804 983 982 +809 1 quad 804 805 984 983 +810 1 quad 805 806 985 984 +811 1 quad 806 807 986 985 +812 1 quad 807 808 987 986 +813 1 quad 808 809 988 987 +814 1 quad 809 810 989 988 +815 1 quad 810 811 990 989 +816 1 quad 811 812 991 990 +817 1 quad 812 813 848 991 +818 1 quad 848 849 992 991 +819 1 quad 849 850 993 992 +820 1 quad 850 851 994 993 +821 1 quad 851 852 995 994 +822 1 quad 852 853 996 995 +823 1 quad 853 854 997 996 +824 1 quad 854 855 998 997 +825 1 quad 855 856 999 998 +826 1 quad 856 857 1000 999 +827 1 quad 857 858 1001 1000 +828 1 quad 858 859 1002 1001 +829 1 quad 859 860 1003 1002 +830 1 quad 860 861 1004 1003 +831 1 quad 861 862 1005 1004 +832 1 quad 862 863 1006 1005 +833 1 quad 863 864 1007 1006 +834 1 quad 864 865 1008 1007 +835 1 quad 865 866 1009 1008 +836 1 quad 866 867 1010 1009 +837 1 quad 867 868 1011 1010 +838 1 quad 868 869 1012 1011 +839 1 quad 869 870 1013 1012 +840 1 quad 870 871 927 1013 +841 1 quad 872 873 925 926 +842 1 quad 873 874 1014 925 +843 1 quad 874 875 1015 1014 +844 1 quad 875 876 1016 1015 +845 1 quad 876 1018 1017 1016 +846 1 quad 876 877 1019 1018 +847 1 quad 877 878 1020 1019 +848 1 quad 878 879 1021 1020 +849 1 quad 879 880 1022 1021 +850 1 quad 880 881 1023 1022 +851 1 quad 881 882 1024 1023 +852 1 quad 882 883 1025 1024 +853 1 quad 883 884 1026 1025 +854 1 quad 884 885 1027 1026 +855 1 quad 885 886 1028 1027 +856 1 quad 886 887 1029 1028 +857 1 quad 887 888 1030 1029 +858 1 quad 888 889 1031 1030 +859 1 quad 889 1033 1032 1031 +860 1 quad 889 890 1034 1033 +861 1 quad 890 891 1035 1034 +862 1 quad 891 892 1036 1035 +863 1 quad 892 893 1037 1036 +864 1 quad 893 894 1038 1037 +865 1 quad 894 895 1039 1038 +866 1 quad 895 896 977 1039 +867 1 quad 795 897 898 794 +868 1 quad 792 793 1040 976 +869 1 quad 1040 793 794 898 +870 1 quad 927 928 1041 1013 +871 1 quad 928 929 1042 1041 +872 1 quad 929 930 1043 1042 +873 1 quad 930 931 1044 1043 +874 1 quad 931 932 1045 1044 +875 1 quad 932 933 1046 1045 +876 1 quad 933 934 1047 1046 +877 1 quad 934 935 1048 1047 +878 1 quad 935 936 1049 1048 +879 1 quad 936 937 1050 1049 +880 1 quad 937 938 1051 1050 +881 1 quad 938 939 1052 1051 +882 1 quad 939 940 1053 1052 +883 1 quad 940 941 1054 1053 +884 1 quad 941 942 1055 1054 +885 1 quad 942 943 1056 1055 +886 1 quad 943 944 1057 1056 +887 1 quad 944 945 1058 1057 +888 1 quad 945 946 1059 1058 +889 1 quad 946 947 1060 1059 +890 1 quad 947 948 1061 1060 +891 1 quad 948 949 1062 1061 +892 1 quad 949 950 1063 1062 +893 1 quad 950 951 1064 1063 +894 1 quad 951 952 1065 1064 +895 1 quad 952 953 1066 1065 +896 1 quad 953 954 1067 1066 +897 1 quad 954 955 1068 1067 +898 1 quad 955 956 1069 1068 +899 1 quad 956 957 1070 1069 +900 1 quad 957 958 1071 1070 +901 1 quad 958 959 1072 1071 +902 1 quad 959 960 1073 1072 +903 1 quad 960 961 1074 1073 +904 1 quad 961 962 1075 1074 +905 1 quad 962 963 1076 1075 +906 1 quad 963 964 1077 1076 +907 1 quad 964 965 1078 1077 +908 1 quad 965 966 1079 1078 +909 1 quad 966 967 1080 1079 +910 1 quad 967 968 1081 1080 +911 1 quad 968 969 1082 1081 +912 1 quad 969 970 1083 1082 +913 1 quad 970 971 1084 1083 +914 1 quad 971 972 1085 1084 +915 1 quad 972 973 1086 1085 +916 1 quad 973 974 1087 1086 +917 1 quad 974 975 1088 1087 +918 1 quad 975 976 899 1088 +919 1 quad 977 978 1089 1039 +920 1 quad 978 979 1090 1089 +921 1 quad 979 980 1091 1090 +922 1 quad 980 981 1092 1091 +923 1 quad 981 982 1093 1092 +924 1 quad 982 983 1094 1093 +925 1 quad 983 984 1095 1094 +926 1 quad 984 985 1096 1095 +927 1 quad 985 986 1097 1096 +928 1 quad 986 987 1098 1097 +929 1 quad 987 988 1099 1098 +930 1 quad 988 989 1100 1099 +931 1 quad 989 990 1101 1100 +932 1 quad 990 991 992 1101 +933 1 quad 899 900 1087 1088 +934 1 quad 900 901 1086 1087 +935 1 quad 901 902 1213 1086 +936 1 quad 902 903 1102 1213 +937 1 quad 903 904 1103 1102 +938 1 quad 904 905 1104 1103 +939 1 quad 905 906 1105 1104 +940 1 quad 906 907 1106 1105 +941 1 quad 907 908 1107 1106 +942 1 quad 908 909 1108 1107 +943 1 quad 909 910 1109 1108 +944 1 quad 910 911 1110 1109 +945 1 quad 911 912 1111 1110 +946 1 quad 912 913 1112 1111 +947 1 quad 913 914 1113 1112 +948 1 quad 914 915 1114 1113 +949 1 quad 915 916 1115 1114 +950 1 quad 916 917 1116 1115 +951 1 quad 917 918 1117 1116 +952 1 quad 918 919 1118 1117 +953 1 quad 919 920 1119 1118 +954 1 quad 920 921 1120 1119 +955 1 quad 921 922 1121 1120 +956 1 quad 922 923 1143 1121 +957 1 quad 923 924 1122 1143 +958 1 quad 924 925 1014 1122 +959 1 quad 992 993 1123 1101 +960 1 quad 993 994 1124 1123 +961 1 quad 994 995 1125 1124 +962 1 quad 995 996 1126 1125 +963 1 quad 996 997 1127 1126 +964 1 quad 997 998 1128 1127 +965 1 quad 998 999 1129 1128 +966 1 quad 999 1000 1130 1129 +967 1 quad 1000 1001 1131 1130 +968 1 quad 1001 1002 1132 1131 +969 1 quad 1002 1003 1133 1132 +970 1 quad 1003 1004 1134 1133 +971 1 quad 1004 1005 1135 1134 +972 1 quad 1005 1006 1136 1135 +973 1 quad 1006 1007 1137 1136 +974 1 quad 1007 1008 1138 1137 +975 1 quad 1008 1009 1139 1138 +976 1 quad 1009 1010 1140 1139 +977 1 quad 1010 1011 1141 1140 +978 1 quad 1011 1012 1142 1141 +979 1 quad 1012 1013 1041 1142 +980 1 quad 1014 1015 1143 1122 +981 1 quad 1015 1016 1121 1143 +982 1 quad 1016 1017 1244 1121 +983 1 quad 1017 1145 1144 1244 +984 1 quad 1017 1018 1146 1145 +985 1 quad 1018 1019 1147 1146 +986 1 quad 1019 1020 1148 1147 +987 1 quad 1020 1021 1149 1148 +988 1 quad 1021 1022 1150 1149 +989 1 quad 1022 1023 1151 1150 +990 1 quad 1023 1024 1152 1151 +991 1 quad 1024 1025 1153 1152 +992 1 quad 1025 1026 1154 1153 +993 1 quad 1026 1027 1155 1154 +994 1 quad 1027 1028 1156 1155 +995 1 quad 1028 1029 1157 1156 +996 1 quad 1029 1030 1158 1157 +997 1 quad 1030 1031 1159 1158 +998 1 quad 1031 1032 1160 1159 +999 1 quad 1032 1162 1161 1160 +1000 1 quad 1032 1033 1163 1162 +1001 1 quad 1033 1034 1164 1163 +1002 1 quad 1034 1035 1165 1164 +1003 1 quad 1035 1036 1166 1165 +1004 1 quad 1036 1037 1167 1166 +1005 1 quad 1037 1038 1168 1167 +1006 1 quad 1038 1039 1089 1168 +1007 1 quad 1041 1042 1169 1142 +1008 1 quad 1042 1043 1170 1169 +1009 1 quad 1043 1044 1171 1170 +1010 1 quad 1044 1045 1172 1171 +1011 1 quad 1045 1046 1173 1172 +1012 1 quad 1046 1047 1174 1173 +1013 1 quad 1047 1048 1175 1174 +1014 1 quad 1048 1049 1176 1175 +1015 1 quad 1049 1050 1177 1176 +1016 1 quad 1050 1051 1178 1177 +1017 1 quad 1051 1052 1179 1178 +1018 1 quad 1052 1053 1180 1179 +1019 1 quad 1053 1054 1181 1180 +1020 1 quad 1054 1055 1182 1181 +1021 1 quad 1055 1056 1183 1182 +1022 1 quad 1056 1057 1184 1183 +1023 1 quad 1057 1058 1185 1184 +1024 1 quad 1058 1059 1186 1185 +1025 1 quad 1059 1060 1187 1186 +1026 1 quad 1060 1061 1188 1187 +1027 1 quad 1061 1062 1189 1188 +1028 1 quad 1062 1063 1190 1189 +1029 1 quad 1063 1064 1191 1190 +1030 1 quad 1064 1065 1192 1191 +1031 1 quad 1065 1066 1193 1192 +1032 1 quad 1066 1067 1194 1193 +1033 1 quad 1067 1068 1195 1194 +1034 1 quad 1068 1069 1196 1195 +1035 1 quad 1069 1070 1197 1196 +1036 1 quad 1070 1071 1198 1197 +1037 1 quad 1071 1072 1199 1198 +1038 1 quad 1072 1073 1200 1199 +1039 1 quad 1073 1074 1201 1200 +1040 1 quad 1074 1075 1202 1201 +1041 1 quad 1075 1076 1203 1202 +1042 1 quad 1076 1077 1204 1203 +1043 1 quad 1077 1078 1205 1204 +1044 1 quad 1078 1079 1206 1205 +1045 1 quad 1079 1080 1207 1206 +1046 1 quad 1080 1081 1208 1207 +1047 1 quad 1081 1082 1209 1208 +1048 1 quad 1082 1083 1210 1209 +1049 1 quad 1083 1084 1211 1210 +1050 1 quad 1084 1085 1212 1211 +1051 1 quad 1085 1086 1213 1212 +1052 1 quad 1089 1090 1214 1168 +1053 1 quad 1090 1091 1215 1214 +1054 1 quad 1091 1092 1216 1215 +1055 1 quad 1092 1093 1217 1216 +1056 1 quad 1093 1094 1218 1217 +1057 1 quad 1094 1095 1219 1218 +1058 1 quad 1095 1096 1220 1219 +1059 1 quad 1096 1097 1221 1220 +1060 1 quad 1097 1098 1222 1221 +1061 1 quad 1098 1099 1223 1222 +1062 1 quad 1099 1100 1224 1223 +1063 1 quad 1100 1101 1123 1224 +1064 1 quad 1213 1102 1225 1212 +1065 1 quad 1102 1103 1226 1225 +1066 1 quad 1103 1104 1227 1226 +1067 1 quad 1104 1105 1228 1227 +1068 1 quad 1105 1106 1229 1228 +1069 1 quad 1106 1107 1230 1229 +1070 1 quad 1107 1108 1231 1230 +1071 1 quad 1108 1109 1232 1231 +1072 1 quad 1109 1110 1233 1232 +1073 1 quad 1110 1111 1234 1233 +1074 1 quad 1111 1112 1235 1234 +1075 1 quad 1112 1113 1236 1235 +1076 1 quad 1113 1114 1237 1236 +1077 1 quad 1114 1115 1238 1237 +1078 1 quad 1115 1116 1239 1238 +1079 1 quad 1116 1117 1240 1239 +1080 1 quad 1117 1118 1241 1240 +1081 1 quad 1118 1119 1242 1241 +1082 1 quad 1119 1120 1243 1242 +1083 1 quad 1120 1121 1244 1243 +1084 1 quad 1123 1124 1245 1224 +1085 1 quad 1124 1125 1246 1245 +1086 1 quad 1125 1126 1247 1246 +1087 1 quad 1126 1127 1248 1247 +1088 1 quad 1127 1128 1249 1248 +1089 1 quad 1128 1129 1250 1249 +1090 1 quad 1129 1130 1251 1250 +1091 1 quad 1130 1131 1252 1251 +1092 1 quad 1131 1132 1253 1252 +1093 1 quad 1132 1133 1254 1253 +1094 1 quad 1133 1134 1255 1254 +1095 1 quad 1134 1135 1256 1255 +1096 1 quad 1135 1136 1257 1256 +1097 1 quad 1136 1137 1258 1257 +1098 1 quad 1137 1138 1259 1258 +1099 1 quad 1138 1139 1260 1259 +1100 1 quad 1139 1140 1261 1260 +1101 1 quad 1140 1141 1262 1261 +1102 1 quad 1141 1142 1169 1262 +1103 1 quad 1244 1144 1263 1243 +1104 1 quad 1243 1263 1241 1242 +1105 1 quad 1263 1264 1354 1241 +1106 1 quad 1263 1144 1265 1264 +1107 1 quad 1144 1145 1266 1265 +1108 1 quad 1146 1147 1266 1145 +1109 1 quad 1147 1148 1267 1266 +1110 1 quad 1148 1149 1268 1267 +1111 1 quad 1149 1150 1269 1268 +1112 1 quad 1150 1151 1270 1269 +1113 1 quad 1151 1152 1271 1270 +1114 1 quad 1152 1153 1272 1271 +1115 1 quad 1154 1155 1274 1273 +1116 1 quad 1155 1156 1730 1274 +1117 1 quad 1156 1157 1275 1730 +1118 1 quad 1157 1158 1276 1275 +1119 1 quad 1158 1159 1277 1276 +1120 1 quad 1159 1160 1278 1277 +1121 1 quad 1160 1161 1279 1278 +1122 1 quad 1161 1281 1280 1279 +1123 1 quad 1161 1162 1282 1281 +1124 1 quad 1162 1163 1283 1282 +1125 1 quad 1163 1164 1284 1283 +1126 1 quad 1164 1165 1285 1284 +1127 1 quad 1165 1166 1286 1285 +1128 1 quad 1166 1167 1287 1286 +1129 1 quad 1167 1168 1214 1287 +1130 1 quad 1153 1154 1273 1288 +1131 1 quad 1272 1153 1288 1289 +1132 1 quad 1169 1170 1290 1262 +1133 1 quad 1170 1171 1291 1290 +1134 1 quad 1171 1172 1292 1291 +1135 1 quad 1172 1173 1293 1292 +1136 1 quad 1173 1174 1294 1293 +1137 1 quad 1174 1175 1295 1294 +1138 1 quad 1175 1176 1296 1295 +1139 1 quad 1176 1177 1297 1296 +1140 1 quad 1177 1178 1298 1297 +1141 1 quad 1178 1179 1299 1298 +1142 1 quad 1179 1180 1300 1299 +1143 1 quad 1180 1181 1301 1300 +1144 1 quad 1181 1182 1302 1301 +1145 1 quad 1182 1183 1303 1302 +1146 1 quad 1183 1184 1304 1303 +1147 1 quad 1184 1185 1305 1304 +1148 1 quad 1185 1186 1306 1305 +1149 1 quad 1186 1187 1307 1306 +1150 1 quad 1187 1188 1308 1307 +1151 1 quad 1188 1189 1309 1308 +1152 1 quad 1189 1190 1310 1309 +1153 1 quad 1190 1191 1311 1310 +1154 1 quad 1191 1192 1312 1311 +1155 1 quad 1192 1193 1313 1312 +1156 1 quad 1193 1194 1314 1313 +1157 1 quad 1194 1195 1315 1314 +1158 1 quad 1195 1196 1316 1315 +1159 1 quad 1196 1197 1317 1316 +1160 1 quad 1197 1198 1318 1317 +1161 1 quad 1198 1199 1319 1318 +1162 1 quad 1199 1200 1320 1319 +1163 1 quad 1200 1201 1321 1320 +1164 1 quad 1201 1202 1322 1321 +1165 1 quad 1202 1203 1323 1322 +1166 1 quad 1203 1204 1324 1323 +1167 1 quad 1204 1205 1325 1324 +1168 1 quad 1205 1206 1326 1325 +1169 1 quad 1206 1207 1327 1326 +1170 1 quad 1207 1208 1328 1327 +1171 1 quad 1208 1209 1329 1328 +1172 1 quad 1209 1210 1330 1329 +1173 1 quad 1214 1215 1331 1287 +1174 1 quad 1215 1216 1332 1331 +1175 1 quad 1216 1217 1333 1332 +1176 1 quad 1217 1218 1334 1333 +1177 1 quad 1218 1219 1335 1334 +1178 1 quad 1219 1220 1336 1335 +1179 1 quad 1220 1221 1337 1336 +1180 1 quad 1221 1222 1338 1337 +1181 1 quad 1222 1223 1339 1338 +1182 1 quad 1223 1224 1245 1339 +1183 1 quad 1212 1225 1340 1211 +1184 1 quad 1211 1340 1330 1210 +1185 1 quad 1225 1226 1330 1340 +1186 1 quad 1226 1227 1329 1330 +1187 1 quad 1227 1228 1341 1329 +1188 1 quad 1228 1229 1342 1341 +1189 1 quad 1229 1230 1343 1342 +1190 1 quad 1230 1231 1344 1343 +1191 1 quad 1231 1232 1345 1344 +1192 1 quad 1232 1233 1346 1345 +1193 1 quad 1233 1234 1347 1346 +1194 1 quad 1234 1235 1348 1347 +1195 1 quad 1235 1236 1349 1348 +1196 1 quad 1236 1237 1350 1349 +1197 1 quad 1237 1238 1351 1350 +1198 1 quad 1238 1239 1352 1351 +1199 1 quad 1239 1240 1353 1352 +1200 1 quad 1240 1241 1354 1353 +1201 1 quad 1245 1246 1355 1339 +1202 1 quad 1246 1247 1356 1355 +1203 1 quad 1247 1248 1357 1356 +1204 1 quad 1248 1249 1358 1357 +1205 1 quad 1249 1250 1359 1358 +1206 1 quad 1250 1251 1360 1359 +1207 1 quad 1251 1252 1361 1360 +1208 1 quad 1252 1253 1362 1361 +1209 1 quad 1253 1254 1363 1362 +1210 1 quad 1254 1255 1364 1363 +1211 1 quad 1255 1256 1365 1364 +1212 1 quad 1256 1257 1366 1365 +1213 1 quad 1257 1258 1367 1366 +1214 1 quad 1258 1259 1368 1367 +1215 1 quad 1259 1260 1369 1368 +1216 1 quad 1260 1261 1370 1369 +1217 1 quad 1261 1262 1290 1370 +1218 1 quad 1329 1341 1371 1328 +1219 1 quad 1341 1342 1372 1371 +1220 1 quad 1342 1343 1373 1372 +1221 1 quad 1343 1344 1374 1373 +1222 1 quad 1344 1345 1375 1374 +1223 1 quad 1345 1346 1376 1375 +1224 1 quad 1346 1347 1377 1376 +1225 1 quad 1347 1348 1378 1377 +1226 1 quad 1348 1349 1379 1378 +1227 1 quad 1349 1350 1380 1379 +1228 1 quad 1350 1351 1381 1380 +1229 1 quad 1351 1352 1382 1381 +1230 1 quad 1352 1353 1383 1382 +1231 1 quad 1353 1354 1384 1383 +1232 1 quad 1354 1264 1385 1384 +1233 1 quad 1264 1265 1386 1385 +1234 1 quad 1265 1266 1387 1386 +1235 1 quad 1266 1267 1388 1387 +1236 1 quad 1267 1268 1389 1388 +1237 1 quad 1268 1269 1390 1389 +1238 1 quad 1269 1270 1391 1390 +1239 1 quad 1270 1271 1392 1391 +1240 1 quad 1271 1272 1393 1392 +1241 1 quad 1272 1289 1394 1393 +1242 1 quad 1289 1396 1395 1394 +1243 1 quad 1289 1288 1397 1396 +1244 1 quad 1288 1273 1398 1397 +1245 1 quad 1328 1371 1399 1327 +1246 1 quad 1371 1372 1476 1399 +1247 1 quad 1372 1373 1475 1476 +1248 1 quad 1373 1374 1474 1475 +1249 1 quad 1374 1375 1473 1474 +1250 1 quad 1375 1376 1400 1473 +1251 1 quad 1376 1377 1401 1400 +1252 1 quad 1377 1378 1402 1401 +1253 1 quad 1378 1379 1403 1402 +1254 1 quad 1379 1380 1404 1403 +1255 1 quad 1380 1381 1405 1404 +1256 1 quad 1381 1382 1406 1405 +1257 1 quad 1382 1383 1407 1406 +1258 1 quad 1383 1384 1408 1407 +1259 1 quad 1384 1385 1409 1408 +1260 1 quad 1385 1386 1410 1409 +1261 1 quad 1386 1387 1411 1410 +1262 1 quad 1387 1388 1412 1411 +1263 1 quad 1388 1389 1413 1412 +1264 1 quad 1389 1390 1414 1413 +1265 1 quad 1390 1391 1415 1414 +1266 1 quad 1391 1392 1416 1415 +1267 1 quad 1392 1393 1417 1416 +1268 1 quad 1393 1394 1418 1417 +1269 1 quad 1394 1395 1419 1418 +1270 1 quad 1395 1421 1420 1419 +1271 1 quad 1395 1396 1422 1421 +1272 1 quad 1396 1397 1423 1422 +1273 1 quad 1397 1398 1424 1423 +1274 1 quad 1398 1731 1425 1424 +1275 1 quad 1398 1273 1274 1731 +1276 1 quad 1275 1276 1428 1427 +1277 1 quad 1276 1277 1429 1428 +1278 1 quad 1277 1278 1430 1429 +1279 1 quad 1278 1279 1431 1430 +1280 1 quad 1279 1280 1432 1431 +1281 1 quad 1280 1434 1433 1432 +1282 1 quad 1280 1281 1435 1434 +1283 1 quad 1281 1282 1436 1435 +1284 1 quad 1282 1283 1437 1436 +1285 1 quad 1283 1284 1438 1437 +1286 1 quad 1284 1285 1439 1438 +1287 1 quad 1285 1286 1440 1439 +1288 1 quad 1286 1287 1331 1440 +1289 1 quad 1290 1291 1441 1370 +1290 1 quad 1291 1292 1442 1441 +1291 1 quad 1292 1293 1443 1442 +1292 1 quad 1293 1294 1444 1443 +1293 1 quad 1294 1295 1445 1444 +1294 1 quad 1295 1296 1446 1445 +1295 1 quad 1296 1297 1447 1446 +1296 1 quad 1297 1298 1448 1447 +1297 1 quad 1298 1299 1449 1448 +1298 1 quad 1299 1300 1450 1449 +1299 1 quad 1300 1301 1451 1450 +1300 1 quad 1301 1302 1452 1451 +1301 1 quad 1302 1303 1453 1452 +1302 1 quad 1303 1304 1454 1453 +1303 1 quad 1304 1305 1455 1454 +1304 1 quad 1305 1306 1456 1455 +1305 1 quad 1306 1307 1457 1456 +1306 1 quad 1307 1308 1458 1457 +1307 1 quad 1308 1309 1459 1458 +1308 1 quad 1309 1310 1460 1459 +1309 1 quad 1310 1311 1461 1460 +1310 1 quad 1311 1312 1462 1461 +1311 1 quad 1312 1313 1463 1462 +1312 1 quad 1313 1314 1464 1463 +1313 1 quad 1314 1315 1465 1464 +1314 1 quad 1315 1316 1466 1465 +1315 1 quad 1316 1317 1467 1466 +1316 1 quad 1317 1318 1468 1467 +1317 1 quad 1318 1319 1469 1468 +1318 1 quad 1319 1320 1470 1469 +1319 1 quad 1320 1321 1471 1470 +1320 1 quad 1321 1322 1472 1471 +1321 1 quad 1322 1323 1473 1472 +1322 1 quad 1323 1324 1474 1473 +1323 1 quad 1324 1325 1475 1474 +1324 1 quad 1325 1326 1476 1475 +1325 1 quad 1476 1326 1327 1399 +1326 1 quad 1331 1332 1477 1440 +1327 1 quad 1332 1333 1478 1477 +1328 1 quad 1333 1334 1479 1478 +1329 1 quad 1334 1335 1480 1479 +1330 1 quad 1335 1336 1481 1480 +1331 1 quad 1336 1337 1482 1481 +1332 1 quad 1337 1338 1483 1482 +1333 1 quad 1338 1339 1355 1483 +1334 1 quad 1355 1356 1484 1483 +1335 1 quad 1356 1357 1485 1484 +1336 1 quad 1357 1358 1486 1485 +1337 1 quad 1358 1359 1487 1486 +1338 1 quad 1359 1360 1488 1487 +1339 1 quad 1360 1361 1489 1488 +1340 1 quad 1361 1362 1490 1489 +1341 1 quad 1362 1363 1491 1490 +1342 1 quad 1363 1364 1492 1491 +1343 1 quad 1364 1365 1493 1492 +1344 1 quad 1365 1366 1494 1493 +1345 1 quad 1366 1367 1495 1494 +1346 1 quad 1367 1368 1496 1495 +1347 1 quad 1368 1369 1497 1496 +1348 1 quad 1369 1370 1441 1497 +1349 1 quad 1400 1401 1624 1565 +1350 1 quad 1401 1402 1498 1624 +1351 1 quad 1402 1403 1499 1498 +1352 1 quad 1403 1404 1500 1499 +1353 1 quad 1404 1405 1501 1500 +1354 1 quad 1405 1406 1502 1501 +1355 1 quad 1406 1407 1503 1502 +1356 1 quad 1407 1408 1504 1503 +1357 1 quad 1410 1411 1505 1504 +1358 1 quad 1411 1412 1506 1505 +1359 1 quad 1412 1413 1507 1506 +1360 1 quad 1413 1414 1508 1507 +1361 1 quad 1414 1415 1509 1508 +1362 1 quad 1415 1416 1510 1509 +1363 1 quad 1416 1417 1511 1510 +1364 1 quad 1417 1418 1512 1511 +1365 1 quad 1418 1419 1513 1512 +1366 1 quad 1419 1420 1514 1513 +1367 1 quad 1420 1516 1515 1514 +1368 1 quad 1420 1421 1517 1516 +1369 1 quad 1421 1422 1518 1517 +1370 1 quad 1422 1423 1519 1518 +1371 1 quad 1423 1424 1520 1519 +1372 1 quad 1428 1429 1522 1521 +1373 1 quad 1429 1430 1523 1522 +1374 1 quad 1430 1431 1524 1523 +1375 1 quad 1431 1432 1525 1524 +1376 1 quad 1432 1433 1526 1525 +1377 1 quad 1433 1528 1527 1526 +1378 1 quad 1433 1434 1529 1528 +1379 1 quad 1434 1435 1530 1529 +1380 1 quad 1435 1436 1531 1530 +1381 1 quad 1436 1437 1532 1531 +1382 1 quad 1437 1438 1533 1532 +1383 1 quad 1438 1439 1534 1533 +1384 1 quad 1439 1440 1477 1534 +1385 1 quad 1441 1442 1535 1497 +1386 1 quad 1442 1443 1536 1535 +1387 1 quad 1443 1444 1537 1536 +1388 1 quad 1444 1445 1538 1537 +1389 1 quad 1445 1446 1539 1538 +1390 1 quad 1446 1447 1540 1539 +1391 1 quad 1447 1448 1541 1540 +1392 1 quad 1448 1449 1542 1541 +1393 1 quad 1449 1450 1543 1542 +1394 1 quad 1450 1451 1544 1543 +1395 1 quad 1451 1452 1545 1544 +1396 1 quad 1452 1453 1546 1545 +1397 1 quad 1453 1454 1547 1546 +1398 1 quad 1454 1455 1548 1547 +1399 1 quad 1455 1456 1549 1548 +1400 1 quad 1456 1457 1550 1549 +1401 1 quad 1457 1458 1551 1550 +1402 1 quad 1458 1459 1552 1551 +1403 1 quad 1459 1460 1553 1552 +1404 1 quad 1460 1461 1554 1553 +1405 1 quad 1461 1462 1555 1554 +1406 1 quad 1462 1463 1556 1555 +1407 1 quad 1463 1464 1557 1556 +1408 1 quad 1464 1465 1558 1557 +1409 1 quad 1465 1466 1559 1558 +1410 1 quad 1466 1467 1560 1559 +1411 1 quad 1467 1468 1561 1560 +1412 1 quad 1468 1469 1562 1561 +1413 1 quad 1469 1470 1563 1562 +1414 1 quad 1470 1471 1564 1563 +1415 1 quad 1471 1472 1565 1564 +1416 1 quad 1477 1478 1566 1534 +1417 1 quad 1478 1479 1567 1566 +1418 1 quad 1479 1480 1568 1567 +1419 1 quad 1480 1481 1569 1568 +1420 1 quad 1481 1482 1570 1569 +1421 1 quad 1482 1483 1484 1570 +1422 1 quad 1484 1485 1571 1570 +1423 1 quad 1485 1486 1572 1571 +1424 1 quad 1486 1487 1573 1572 +1425 1 quad 1487 1488 1574 1573 +1426 1 quad 1488 1489 1575 1574 +1427 1 quad 1489 1490 1576 1575 +1428 1 quad 1490 1491 1577 1576 +1429 1 quad 1491 1492 1578 1577 +1430 1 quad 1492 1493 1579 1578 +1431 1 quad 1493 1494 1580 1579 +1432 1 quad 1494 1495 1581 1580 +1433 1 quad 1495 1496 1582 1581 +1434 1 quad 1496 1497 1535 1582 +1435 1 quad 1521 1522 1583 1520 +1436 1 quad 1522 1523 1584 1583 +1437 1 quad 1523 1524 1585 1584 +1438 1 quad 1524 1525 1586 1585 +1439 1 quad 1525 1526 1587 1586 +1440 1 quad 1526 1527 1588 1587 +1441 1 quad 1527 1590 1589 1588 +1442 1 quad 1527 1528 1591 1590 +1443 1 quad 1528 1529 1592 1591 +1444 1 quad 1529 1530 1593 1592 +1445 1 quad 1530 1531 1594 1593 +1446 1 quad 1531 1532 1595 1594 +1447 1 quad 1532 1533 1596 1595 +1448 1 quad 1533 1534 1566 1596 +1449 1 quad 1535 1536 1597 1582 +1450 1 quad 1536 1537 1598 1597 +1451 1 quad 1537 1538 1599 1598 +1452 1 quad 1538 1539 1600 1599 +1453 1 quad 1539 1540 1601 1600 +1454 1 quad 1540 1541 1602 1601 +1455 1 quad 1541 1542 1603 1602 +1456 1 quad 1542 1543 1604 1603 +1457 1 quad 1543 1544 1605 1604 +1458 1 quad 1544 1545 1606 1605 +1459 1 quad 1545 1546 1607 1606 +1460 1 quad 1546 1547 1608 1607 +1461 1 quad 1547 1548 1609 1608 +1462 1 quad 1548 1549 1610 1609 +1463 1 quad 1549 1550 1611 1610 +1464 1 quad 1550 1551 1515 1611 +1465 1 quad 1551 1552 1651 1515 +1466 1 quad 1552 1553 1612 1651 +1467 1 quad 1553 1554 1613 1612 +1468 1 quad 1554 1555 1614 1613 +1469 1 quad 1555 1556 1615 1614 +1470 1 quad 1556 1557 1616 1615 +1471 1 quad 1557 1558 1617 1616 +1472 1 quad 1558 1559 1618 1617 +1473 1 quad 1559 1560 1619 1618 +1474 1 quad 1560 1561 1620 1619 +1475 1 quad 1561 1562 1621 1620 +1476 1 quad 1562 1563 1622 1621 +1477 1 quad 1563 1564 1623 1622 +1478 1 quad 1564 1565 1624 1623 +1479 1 quad 1566 1567 1625 1596 +1480 1 quad 1567 1568 1626 1625 +1481 1 quad 1568 1569 1627 1626 +1482 1 quad 1569 1570 1571 1627 +1483 1 quad 1571 1572 1628 1627 +1484 1 quad 1572 1573 1629 1628 +1485 1 quad 1573 1574 1630 1629 +1486 1 quad 1574 1575 1631 1630 +1487 1 quad 1575 1576 1632 1631 +1488 1 quad 1576 1577 1633 1632 +1489 1 quad 1577 1578 1634 1633 +1490 1 quad 1578 1579 1635 1634 +1491 1 quad 1579 1580 1636 1635 +1492 1 quad 1580 1581 1637 1636 +1493 1 quad 1581 1582 1597 1637 +1494 1 quad 1624 1498 1638 1623 +1495 1 quad 1498 1499 1674 1638 +1496 1 quad 1499 1500 1673 1674 +1497 1 quad 1501 1502 1640 1639 +1498 1 quad 1502 1503 1641 1640 +1499 1 quad 1505 1506 1642 1641 +1500 1 quad 1506 1507 1643 1642 +1501 1 quad 1507 1508 1644 1643 +1502 1 quad 1508 1509 1645 1644 +1503 1 quad 1509 1510 1646 1645 +1504 1 quad 1510 1511 1647 1646 +1505 1 quad 1511 1512 1648 1647 +1506 1 quad 1512 1513 1649 1648 +1507 1 quad 1513 1514 1650 1649 +1508 1 quad 1514 1515 1651 1650 +1509 1 quad 1520 1583 1652 1519 +1510 1 quad 1583 1584 1653 1652 +1511 1 quad 1584 1585 1654 1653 +1512 1 quad 1585 1586 1655 1654 +1513 1 quad 1586 1587 1656 1655 +1514 1 quad 1587 1588 1657 1656 +1515 1 quad 1588 1589 1658 1657 +1516 1 quad 1589 1660 1659 1658 +1517 1 quad 1589 1590 1661 1660 +1518 1 quad 1590 1591 1662 1661 +1519 1 quad 1591 1592 1663 1662 +1520 1 quad 1592 1593 1664 1663 +1521 1 quad 1593 1594 1665 1664 +1522 1 quad 1594 1595 1666 1665 +1523 1 quad 1595 1596 1625 1666 +1524 1 quad 1651 1612 1667 1650 +1525 1 quad 1612 1613 1696 1667 +1526 1 quad 1613 1614 1695 1696 +1527 1 quad 1614 1615 1694 1695 +1528 1 quad 1615 1616 1668 1694 +1529 1 quad 1616 1617 1669 1668 +1530 1 quad 1617 1618 1670 1669 +1531 1 quad 1618 1619 1671 1670 +1532 1 quad 1619 1620 1672 1671 +1533 1 quad 1620 1621 1673 1672 +1534 1 quad 1621 1622 1674 1673 +1535 1 quad 1674 1622 1623 1638 +1536 1 quad 1597 1598 1675 1637 +1537 1 quad 1598 1599 1676 1675 +1538 1 quad 1599 1600 1677 1676 +1539 1 quad 1600 1601 1678 1677 +1540 1 quad 1601 1602 1679 1678 +1541 1 quad 1602 1603 1680 1679 +1542 1 quad 1603 1604 1681 1680 +1543 1 quad 1604 1605 1682 1681 +1544 1 quad 1605 1606 1683 1682 +1545 1 quad 1606 1607 1684 1683 +1546 1 quad 1607 1608 1685 1684 +1547 1 quad 1608 1609 1686 1685 +1548 1 quad 1609 1610 1687 1686 +1549 1 quad 1610 1611 1688 1687 +1550 1 quad 1611 1515 1516 1688 +1551 1 quad 1625 1626 1689 1666 +1552 1 quad 1626 1627 1628 1689 +1553 1 quad 1628 1629 1719 1689 +1554 1 quad 1629 1630 1718 1719 +1555 1 quad 1630 1631 1717 1718 +1556 1 quad 1631 1632 1716 1717 +1557 1 quad 1632 1633 1715 1716 +1558 1 quad 1633 1634 1714 1715 +1559 1 quad 1634 1635 1713 1714 +1560 1 quad 1635 1636 1690 1713 +1561 1 quad 1636 1637 1675 1690 +1562 1 quad 1516 1517 1707 1688 +1563 1 quad 1517 1518 1706 1707 +1564 1 quad 1518 1519 1652 1706 +1565 1 quad 1641 1642 1711 1640 +1566 1 quad 1642 1643 1710 1711 +1567 1 quad 1643 1644 1691 1710 +1568 1 quad 1644 1645 1692 1691 +1569 1 quad 1645 1646 1693 1692 +1570 1 quad 1646 1647 1694 1693 +1571 1 quad 1647 1648 1695 1694 +1572 1 quad 1648 1649 1696 1695 +1573 1 quad 1696 1649 1650 1667 +1574 1 quad 1675 1676 1697 1690 +1575 1 quad 1676 1677 1698 1697 +1576 1 quad 1677 1678 1699 1698 +1577 1 quad 1678 1679 1700 1699 +1578 1 quad 1679 1680 1701 1700 +1579 1 quad 1680 1681 1702 1701 +1580 1 quad 1681 1682 1703 1702 +1581 1 quad 1682 1683 1704 1703 +1582 1 quad 1683 1684 1705 1704 +1583 1 quad 1684 1685 1706 1705 +1584 1 quad 1685 1686 1707 1706 +1585 1 quad 1688 1707 1686 1687 +1586 1 quad 1694 1668 1708 1693 +1587 1 quad 1668 1669 1709 1708 +1588 1 quad 1669 1670 1727 1709 +1589 1 quad 1670 1671 1729 1727 +1590 1 quad 1671 1672 1712 1729 +1591 1 quad 1652 1653 1705 1706 +1592 1 quad 1653 1654 1704 1705 +1593 1 quad 1654 1655 1703 1704 +1594 1 quad 1655 1656 1702 1703 +1595 1 quad 1656 1657 1701 1702 +1596 1 quad 1657 1658 1700 1701 +1597 1 quad 1658 1659 1699 1700 +1598 1 quad 1699 1659 1713 1698 +1599 1 quad 1659 1660 1714 1713 +1600 1 quad 1660 1661 1715 1714 +1601 1 quad 1661 1662 1716 1715 +1602 1 quad 1662 1663 1717 1716 +1603 1 quad 1663 1664 1718 1717 +1604 1 quad 1664 1665 1719 1718 +1605 1 quad 1719 1665 1666 1689 +1606 1 quad 1713 1690 1697 1698 +1607 1 quad 1691 1692 1721 1720 +1608 1 quad 1693 1708 1721 1692 +1609 1 quad 1708 1709 1720 1721 +1610 1 quad 420 419 344 345 +1611 1 quad 345 346 347 420 +1612 1 quad 872 926 1722 846 +1613 1 quad 926 845 1723 1722 +1614 1 quad 1722 1723 1725 1724 +1615 1 quad 846 1722 1724 847 +1616 1 quad 847 1724 717 718 +1617 1 quad 1724 1725 716 717 +1618 1 quad 1504 1505 1641 1503 +1619 1 quad 1504 1408 1409 1410 +1620 1 quad 1521 1520 1424 1425 +1621 1 quad 1425 1427 1428 1521 +1622 1 quad 1473 1400 1565 1472 +1623 1 quad 1691 1720 1726 1710 +1624 1 quad 1720 1709 1727 1726 +1625 1 quad 1726 1727 1729 1728 +1626 1 quad 1710 1726 1728 1711 +1627 1 quad 1711 1728 1639 1640 +1628 1 quad 1728 1729 1712 1639 +1629 1 quad 1274 1730 1426 1731 +1630 1 quad 1730 1275 1427 1426 +1631 1 quad 1426 1427 1425 1731 +1632 1 quad 1712 1672 1673 1500 +1633 1 quad 1712 1500 1501 1639 +1633 1 line 67 180 +1634 1 line 180 179 +1635 1 line 179 178 +1636 1 line 178 177 +1637 1 line 177 176 +1638 1 line 176 175 +1639 1 line 175 174 +1640 1 line 174 173 +1641 1 line 173 172 +1642 1 line 172 171 +1643 1 line 171 170 +1644 1 line 170 169 +1645 1 line 169 168 +1646 1 line 168 167 +1647 1 line 167 166 +1648 1 line 166 165 +1649 1 line 165 164 +1650 1 line 164 163 +1651 1 line 163 162 +1652 1 line 162 161 +1653 1 line 161 160 +1654 1 line 160 159 +1655 1 line 159 158 +1656 1 line 158 157 +1657 1 line 157 156 +1658 1 line 156 155 +1659 1 line 155 154 +1660 1 line 154 153 +1661 1 line 153 152 +1662 1 line 152 151 +1663 1 line 151 150 +1664 1 line 150 149 +1665 1 line 149 148 +1666 1 line 148 147 +1667 1 line 147 146 +1668 1 line 146 145 +1669 1 line 145 144 +1670 1 line 144 143 +1671 1 line 143 142 +1672 1 line 142 141 +1673 1 line 141 140 +1674 1 line 140 139 +1675 1 line 139 138 +1676 1 line 138 137 +1677 1 line 137 136 +1678 1 line 136 135 +1679 1 line 135 134 +1680 2 line 66 34 +1681 2 line 9 35 +1682 2 line 35 36 +1683 2 line 36 37 +1684 2 line 37 38 +1685 2 line 38 39 +1686 2 line 39 40 +1687 2 line 40 41 +1688 2 line 41 42 +1689 2 line 42 43 +1690 2 line 43 44 +1691 2 line 44 45 +1692 2 line 45 46 +1693 2 line 46 47 +1694 2 line 47 48 +1695 2 line 48 49 +1696 2 line 49 50 +1697 2 line 50 51 +1698 2 line 51 52 +1699 2 line 52 53 +1700 2 line 53 54 +1701 2 line 54 55 +1702 2 line 55 56 +1703 2 line 56 57 +1704 2 line 57 58 +1705 2 line 58 59 +1706 2 line 59 60 +1707 2 line 60 61 +1708 2 line 61 62 +1709 2 line 62 63 +1710 2 line 63 64 +1711 2 line 64 65 +1712 2 line 65 66 +1713 3 line 133 67 +1714 3 line 132 133 +1715 3 line 131 132 +1716 3 line 34 68 +1717 3 line 68 69 +1718 3 line 69 70 +1719 3 line 70 71 +1720 3 line 71 72 +1721 3 line 72 73 +1722 3 line 73 74 +1723 3 line 74 75 +1724 3 line 75 76 +1725 3 line 76 77 +1726 3 line 77 78 +1727 3 line 78 79 +1728 3 line 79 80 +1729 3 line 80 81 +1730 3 line 81 82 +1731 3 line 82 83 +1732 3 line 83 84 +1733 3 line 84 85 +1734 3 line 85 86 +1735 3 line 86 87 +1736 3 line 87 88 +1737 3 line 88 89 +1738 3 line 89 90 +1739 3 line 90 91 +1740 3 line 91 92 +1741 3 line 92 93 +1742 3 line 93 94 +1743 3 line 94 95 +1744 3 line 95 96 +1745 3 line 96 97 +1746 3 line 97 98 +1747 3 line 98 99 +1748 3 line 99 100 +1749 3 line 100 101 +1750 3 line 101 102 +1751 3 line 102 103 +1752 3 line 103 104 +1753 3 line 104 105 +1754 3 line 105 106 +1755 3 line 106 107 +1756 3 line 107 108 +1757 3 line 108 109 +1758 3 line 109 110 +1759 3 line 110 111 +1760 3 line 111 112 +1761 3 line 112 113 +1762 3 line 113 114 +1763 3 line 114 115 +1764 3 line 115 116 +1765 3 line 116 117 +1766 3 line 117 118 +1767 3 line 118 119 +1768 3 line 119 120 +1769 3 line 120 121 +1770 3 line 121 122 +1771 3 line 122 123 +1772 3 line 123 124 +1773 3 line 124 125 +1774 3 line 125 126 +1775 3 line 126 127 +1776 3 line 127 128 +1777 3 line 128 129 +1778 3 line 129 130 +1779 3 line 130 131 +1780 4 line 185 181 +1781 4 line 134 182 +1782 4 line 182 183 +1783 4 line 183 184 +1784 4 line 184 185 +1785 4 line 181 187 +1786 4 line 187 188 +1787 4 line 188 189 +1788 4 line 189 190 +1789 4 line 190 191 +1790 4 line 191 192 +1791 4 line 192 186 +1792 4 line 186 193 +1793 4 line 193 194 +1794 4 line 194 1 +1795 4 line 1 3 +1796 4 line 3 4 +1797 4 line 4 5 +1798 4 line 5 6 +1799 4 line 6 7 +1800 4 line 7 8 +1801 4 line 8 2 +1802 4 line 2 10 +1803 4 line 10 11 +1804 4 line 11 12 +1805 4 line 12 13 +1806 4 line 13 14 +1807 4 line 14 15 +1808 4 line 15 16 +1809 4 line 16 17 +1810 4 line 17 18 +1811 4 line 18 19 +1812 4 line 19 20 +1813 4 line 20 21 +1814 4 line 21 22 +1815 4 line 22 23 +1816 4 line 23 24 +1817 4 line 24 25 +1818 4 line 25 26 +1819 4 line 26 27 +1820 4 line 27 28 +1821 4 line 28 29 +1822 4 line 29 30 +1823 4 line 30 31 +1824 4 line 31 32 +1825 4 line 32 33 +1826 4 line 33 9 + diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc new file mode 100644 index 0000000000..bb2389e563 --- /dev/null +++ b/deal.II/examples/step-33/step-33.cc @@ -0,0 +1,2114 @@ +/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */ +/* Author: David Neckels, Boulder Colorado 2007 */ +/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */ +/* Version: Version-5-2-0 */ +/* */ +/* Copyright (C) 2001, 2002, 2003, 2004, 2005 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. */ + + // This program solves the Euler equations + // of gas dynamics for a given configuration + // file. It uses a standard Galerkin approach + // with weakly applied boundary conditions. + + // @sect3{Include files} + + // Aztecoo require mpi (even though we run on only + // one processor in this example). +#include + + // Here we have the necessary TRILINOS includes. + // + // Epetra is the basic trilinos vector/matrix library. +#include +#include +#include +#include +#include + // Teuchos is a Trilinos utility library that is used + // to set parameters within the Aztec solver library. +#include "Teuchos_ParameterList.hpp" + // Aztec is the iterative solver library. +#include +#include +#define HAVE_IFPACK_TEUCHOS +#include + + // Amesos is a direct solver package within Trilinos. +#include + // Sacado is the automatic differentiation package, which + // is used to find the jacobian for a fully implicit Newton + // iteration. +#include + + // A standard set of dealii includes. Nothing special to + // comment on here. +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + // And this again is C++: +#include +#include +#include + + // Introduce the dealii library into the current namespace. +using namespace dealii; + + // We define a shorter name for the automatic differentiation + // type. +typedef Sacado::Fad::DFad fad_double; +typedef unsigned int UInt; + // The Epetra library requires a 'communicator', which describes + // the layout of a parallel (or serial) set of processors. +Epetra_MpiComm *Comm; + + //@sect3{Flux function definition} + // Here we define the flux function for this system of conservation + // laws. Note: it would be terribly difficult to use this example + // to solve some other system of conservation laws. + // + // We define the number of components in the system. Euler's has + // one entry for momenta in each spatial direction, plus the energy + // and density components. +#define N_COMP (2 + DIMENSION) + // Define a handle to the density and energy indices. We have arrange + // the momenta to be first, then density, and, lastly, energy. +#define DENS_IDX DIMENSION +#define ENERGY_IDX (DIMENSION+1) + + // The gas constant. This value is representative of air. +const double GAMMA = 1.4; + // We define the flux functions as one large matrix. Each row of this + // matrix represents a scalar conservation law for the component in + // that row. We template the numerical type of the flux function + // so that we may use the automatic differentiation type here. + // The flux functions are defined in terms of the + // conserved variables $\rho w_0, \dots, \rho w_{d-1}, \rho, E$, + // so they do not look exactly like the Euler equations one is + // used to seeing. We evaluate the flux at a single quadrature + // point. +template +void Flux(std::vector > &flux, + const Point &point, + const std::vector &W) +{ + + // Pressure is a dependent variable: $p = + // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$. + number rho_normVsqr; + for (int d0 = 0; d0 < dim; d0++) rho_normVsqr += W[d0]*W[d0]; + // Since W are $\rho v$, we get a $\rho^2$ in the + // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$. + rho_normVsqr /= W[DENS_IDX]; + + number pressure = (GAMMA-1.0)*(W[ENERGY_IDX] - number(0.5)*(rho_normVsqr)); + + // We compute the momentum terms. We divide by the + // density here to get $v_i \rho v_j$ + for (int d = 0; d < dim; d++) { + for (int d1 = 0; d1 < dim; d1++) { + flux[d][d1] = W[d]*W[d1]/W[DENS_IDX]; + } + // The pressure contribution, along the diagonal: + flux[d][d] += pressure; + // Advection/conservation of density: + flux[DENS_IDX][d] = W[d]; + // And, lastly, conservation of energy. + flux[ENERGY_IDX][d] = W[d]/W[DENS_IDX]* + (W[ENERGY_IDX] + pressure); // energy + } +} + + // On the boundaries of the domain and across `hanging nodes` we use + // a numerical flux function to enforce boundary conditions. This routine + // is the basic Lax-Friedrich's flux with a stabilization parameter + // $\alpha$. +template +void LFNumFlux( + std::vector > &nflux, + const std::vector > &points, + const std::vector > &normals, + const std::vector > &Wplus, + const std::vector > &Wminus, + double alpha) +{ + int n_q_points = points.size(); + + // We evaluate the flux at each of the quadrature points. + for (int q = 0; q < n_q_points; q++) { + std::vector > iflux(N_COMP, + std::vector(dim, 0)); + std::vector > oflux(N_COMP, + std::vector(dim, 0)); + + Flux(iflux, points[q], Wplus[q]); + Flux(oflux, points[q], Wminus[q]); + + for (int di = 0; di < N_COMP; di++) { + nflux[q][di] = 0; + for (int d = 0; d < dim; d++) { + nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d); + } + nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]); + } + } + +} + + // @sect3{Initial and side condition parsing} + // For the initial condition we use the expression parser function + // object. +template +class InitialCondition : public FunctionParser +{ + public: + InitialCondition (); + + // This function should be called after parsing, but before using + // the object. It formalizes the expressions and initializes the + // function parser with the appropriate expressions. + void Init(); + + // During parsing we call this function as the initial condition + // for one of the $\mathbf{w}$ variables is encountered. + + void set_ic(int _row, std::string &expr) { + expressions[_row] = expr; + } + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + std::vector expressions; +}; + +template +InitialCondition::InitialCondition () : + FunctionParser (N_COMP), + expressions(N_COMP, "0.0") +{} + + // Here we set up x,y,z as the variables that one should use in the input + // deck to describe their initial condition. +template +void InitialCondition::Init() { + std::map constants; + constants["M_PI"] = M_PI; + std::string variables = (dim == 2 ? "x,y" : "x,y,z"); + + FunctionParser::initialize(variables, expressions, constants); + +} + +template +void InitialCondition::vector_value_list (const std::vector > &points, + std::vector > &value_list) const +{ + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); +} + + // As above, we use the expression function parser for boundary conditions. +template +class SideCondition : public FunctionParser +{ + public: + SideCondition (int ncomp); + ~SideCondition (); + + // As above. + void Init(); + // As above. + void set_coeff_row(int _row_n, std::string &expr); + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + std::vector expressions; +}; + +template +SideCondition::SideCondition (int ncomp) : + FunctionParser (ncomp), + expressions(ncomp, "0.0") +{ +} +template +void SideCondition::set_coeff_row (int _row_n, std::string &expr) +{ + expressions[_row_n] = expr; +} + +template +void SideCondition::Init() { + std::map constants; + constants["M_PI"] = M_PI; + std::string variables = (dim == 2 ? "x,y" : "x,y,z"); + + FunctionParser::initialize(variables, expressions, constants); + +} + +template +SideCondition::~SideCondition () +{ +} + +template +void SideCondition::vector_value_list (const std::vector > &points, + std::vector > &value_list) const +{ + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); +} + //@sect3{Conservation Law class} + // Here we define a Conservation Law class that helps group + // operations and data for our Euler equations into a manageable + // entity. Functions will be described as their definitions appear. +template +class ConsLaw +{ + public: + ConsLaw (); + ~ConsLaw (); + + void run (); + void declare_parameters(); + void load_parameters(const char *); + + private: + void build_fe(); + void setup_system (); + void initialize_system (); + void assemble_system (double &res_norm); + void solve (Vector &solution, int &, double &); + void refine_grid (); + void output_results (const unsigned int cycle) const; + void initialize(); + void zero_matrix(); + void estimate(); + void postprocess(); + void compute_predictor(); + + Triangulation triangulation; + const MappingQ1 mapping; + + + FESystem *fe_ptr; + + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + const QGauss quadrature; + const QGauss face_quadrature; + + // The actual solution to the Euler equation + Vector solution; + // The current value of the solution during the Newton iteration + Vector nlsolution; + // An estimate of the next time value; used for adaptivity and as a + // guess for the next Newton iteration. + Vector predictor; + // Values after post-processing (used to output the physical variables). + Vector ppsolution; + // The solution to the linear problem during the Newton iteration + Vector dsolution; + Vector right_hand_side; + + public: + + void assemble_cell_term(const FEValues& fe_v, + std::vector &dofs, + unsigned int cell_no + ); + + void assemble_face_term( + int face_no, + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + std::vector &dofs, + std::vector &dofs_neighbor, + int boundary = -1 + ); + + unsigned int get_n_components() const { return N_COMP;} + + private: + // T = current time, dT = time step, TF = final time. + double T, dT, TF; + double face_diameter; + double cell_diameter; + // An object to handle parsing the input deck. + ParameterHandler prm; + // Name of the mesh to read in. + string mesh; + InitialCondition ic; + + // Enums for the various supported boundary conditions. + typedef enum {INFLOW_BC = 1, OUTFLOW_BC=2, NO_PENETRATION_BC=3, PRESSURE_BC=4} bc_type; + + // For each boundary we store a map from boundary # to the type + // of boundary condition. If the boundary condition is prescribed, + // we store a pointer to a function object that will hold the expression + // for that boundary condition. + typedef typename std::map, Function*> > bdry_map_type; + bdry_map_type bdry_map; + + void add_boundary(unsigned int bd, std::vector& flags, Function *bf); + + // An object to store parameter information about the Aztec solver. + typedef struct { + int LIN_OUTPUT; + typedef enum { GMRES = 0, DIRECT = 1} solver_type; + solver_type SOLVER; + typedef enum { QUIET = 0, VERBOSE = 1 } output_type; + output_type OUTPUT; + // Linear residual tolerance. + double RES; + int MAX_ITERS; + // We use the ILUT preconditioner. This is similar + // to the ILU. FILL is the number of extra entries + // to add when forming the ILU decomposition. + double ILUT_FILL; + // When forming the preconditioner, for certain problems + // bad conditioning (or just bad luck) can cause the + // preconditioner to be very poorly conditioned. Hence + // it can help to add diagonal perturbations to the + // original matrix and form the preconditioner for this + // slightly better matrix. ATOL is an absolute perturbation + // that is added to the diagonal before forming the + // prec, and RTOL is a scaling factor $rtol >= 1$. + double ILUT_ATOL; + double ILUT_RTOL; + // The ILUT will drop any values that have magnitude less + // than this value. This is a way to + // manage the amount of memory used by this preconditioner. + double ILUT_DROP; + } solver_params_type; + + solver_params_type solver_params; + + // Some refinement parameters. + typedef struct { + typedef enum { NONE = 0, FIXED_NUMBER = 1, SHOCK = 2} refine_type; + double high_frac; + double low_frac; + refine_type refine; + double high_frac_sav; + double max_cells; + double shock_val; + double shock_levels; + } refinement_params_type; + + refinement_params_type refinement_params; + + // The user can set the stabilization parameter $\alpha$ + // in the Lax-Friedrich's flux. + typedef struct { + typedef enum {CONSTANT=1,MESH=2} LF_stab_type; + LF_stab_type LF_stab; + double LF_stab_value; + } flux_params_type; + + flux_params_type flux_params; + + bool is_stationary; + + // Power for the mesh stabilization term. + double diffusion_power; + double gravity; + // If true, we output the squared gradient of the + // density instead of density. Using this one can + // create shock plots. + bool schlieren_plot; + // How often to create an output file. + double output_step; + + Epetra_CrsMatrix *Matrix; + Epetra_Map *Map; + Vector indicator; + + // Crank-Nicolson value + const double theta; + +}; + + + // Asign a row of the conservation law a specified + // boundary type and (possibly) function. +template +void ConsLaw::add_boundary(unsigned int bd, + std::vector &flags, Function *bf) { + + std::pair, Function *> entry(flags, bf); + bdry_map[bd] = entry; +} + + + // Apply the initialial condition. Simultaneously + // initialize the non-linear solution. +template +void ConsLaw::initialize() { + VectorTools::interpolate(dof_handler, + ic, solution); + VectorTools::interpolate(dof_handler, + ic, nlsolution); + +} + + // @sect3{Assembly} + // @sect4{Function: assemble_cell_term} + // + // Assembles the cell term, adding minus the residual + // to the right hand side, and adding in the Jacobian + // contributions. +template +void ConsLaw::assemble_cell_term( + const FEValues &fe_v, + std::vector &dofs, + unsigned int cell_no + ) +{ + // The residual for each row (i) will be accumulating + // into this fad variable. At the end of the assembly + // for this row, we will query for the sensitivities + // to this variable and add them into the Jacobian. + fad_double F_i; + int dofs_per_cell = fe_v.dofs_per_cell; + int n_q_points = fe_v.n_quadrature_points; + + // We will define the dofs on this cell in these fad variables. + std::vector DOF(dofs_per_cell); + + // Values of the conservative variables at the quadrature points. + std::vector > W (n_q_points, + std::vector(get_n_components())); + + // Values at the last time step of the conservative variables. + // Note that these do not use fad variables, since they do + // not depend on the 'variables to be sought'=DOFS. + std::vector > Wl (n_q_points, + std::vector(get_n_components())); + + // Here we will hold the averaged values of the conservative + // variables that we will linearize around (cn=Crank Nicholson). + std::vector > Wcn (n_q_points, + std::vector(get_n_components())); + + // Gradients of the current variables. It is a + // bit of a shame that we have to compute these; we almost don't. + // The nice thing about a simple conservation law is that the + // the flux doesn't generally involve any gradients. We do + // need these, however, for the diffusion stabilization. + std::vector > > Wgrads (n_q_points, + std::vector >(get_n_components(), + std::vector(dim))); + + + const std::vector &JxW = fe_v.get_JxW_values (); + + + // Here is the magical point where we declare a subset + // of the fad variables as degrees of freedom. All + // calculations that reference these variables (either + // directly or indirectly) will accumulate sensitivies + // with respect to these dofs. + for (int in = 0; in < dofs_per_cell; in++) { + DOF[in] = nlsolution(dofs[in]); + DOF[in].diff(in, dofs_per_cell); + } + + // Here we compute the shape function values and gradients + // at the quadrature points. Ideally, we could call into + // something like get_function_values, get_function_grads, + // but since we don't want to make the entire solution vector + // fad types, only the local cell variables, we explicitly + // code this loop; + for (int q = 0; q < n_q_points; q++) { + for (int di = 0; di < get_n_components(); di++) { + W[q][di] = 0; + Wl[q][di] = 0; + Wcn[q][di] = 0; + for (int d = 0; d < dim; d++) { + Wgrads[q][di][d] = 0; + } + } + for (int sf = 0; sf < dofs_per_cell; sf++) { + int di = fe_v.get_fe().system_to_component_index(sf).first; + W[q][di] += + DOF[sf]*fe_v.shape_value_component(sf, q, di); + Wl[q][di] += + solution(dofs[sf])*fe_v.shape_value_component(sf, q, di); + Wcn[q][di] += + (theta*DOF[sf]+(1-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di); + + for (int d = 0; d < dim; d++) { + Wgrads[q][di][d] += DOF[sf]* + fe_v.shape_grad_component(sf, q, di)[d]; + } // for d + + } + + } // for q + + // Gather the flux values for all components at + // all of the quadrature points. This also + // computes the matrix of sensitivities. Perhaps + // this could be done in a better way, since this + // could be a rather large object, but for now it + // seems to work just fine. + std::vector > > flux(n_q_points, + std::vector >(get_n_components(), + std::vector(dim, 0))); + + for (unsigned int q=0; q < n_q_points; ++q) { + Flux(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]); + } + + // We now have all of the function values/grads/fluxes, + // so perform the assembly. We have an outer loop + // through the components of the system, and an + // inner loop over the quadrature points, where we + // accumulate contributions to the ith residual. + // + // We initialy sum all contributions of the residual + // in the positive sense, so that we don't need to + // negative the Jacobian entries. Then, when we sum + // into the right_hand_side vector, + // we negate this residual. + for (unsigned int i=0; iSumIntoGlobalValues(dofs[i], + dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + + // Add minus the residual to the right hand side. + right_hand_side(dofs[i]) -= F_i.val(); + + } // for i +} + // @sect4{Function: assemble_face_term} + // These are either + // boundary terms or terms across differing + // levels of refinement. In the first case, + // fe_v==fe_v_neighbor and dofs==dofs_neighbor. + // The int boundary < 0 if not at a boundary, + // otherwise it is the boundary indicator. +template +void ConsLaw::assemble_face_term( + int face_no, + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + std::vector &dofs, + std::vector &dofs_neighbor, + int boundary + ) +{ + fad_double F_i; + const unsigned int n_q_points = fe_v.n_quadrature_points; + const unsigned int dofs_per_cell = fe_v.get_fe().dofs_per_cell; + const unsigned int ndofs_per_cell = fe_v_neighbor.get_fe().dofs_per_cell; + Assert(dofs_per_cell == ndofs_per_cell, + ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell)); + + // As above, the fad degrees of freedom + std::vector DOF(dofs_per_cell+ndofs_per_cell); + + // The conservative variables for this cell, + // and for + std::vector > Wplus (n_q_points, + std::vector(get_n_components())); + std::vector > Wminus (n_q_points, + std::vector(get_n_components())); + + + const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector > &normals = fe_v.get_normal_vectors (); + + + // If we are at a boundary, then dofs_neighbor are + // the same as dofs, so we do not want to duplicate them. + // If there is a neighbor cell, then we want to include + // them. + int ndofs = (boundary < 0 ? dofs_per_cell + ndofs_per_cell : dofs_per_cell); + // Set the local DOFS. + for (int in = 0; in < dofs_per_cell; in++) { + DOF[in] = nlsolution(dofs[in]); + DOF[in].diff(in, ndofs); + } + // If present, set the neighbor dofs. + if (boundary < 0) + for (int in = 0; in < ndofs_per_cell; in++) { + DOF[in+dofs_per_cell] = nlsolution(dofs_neighbor[in]); + DOF[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs); + } + + // Set the values of the local conservative variables. + // Initialize all variables to zero. + for (int q = 0; q < n_q_points; q++) { + for (int di = 0; di < get_n_components(); di++) { + Wplus[q][di] = 0; + Wminus[q][di] = 0; + } + for (int sf = 0; sf < dofs_per_cell; sf++) { + int di = fe_v.get_fe().system_to_component_index(sf).first; + Wplus[q][di] += + (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di); + } + + + // If there is a cell across, then initialize + // the exterior trace as a function of the other + // cell degrees of freedom. + if (boundary < 0) { + for (int sf = 0; sf < ndofs_per_cell; sf++) { + int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first; + Wminus[q][di] += + (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))* + fe_v_neighbor.shape_value_component(sf, q, di); + } + } + } // for q + + // If this is a boundary, then the values of $W^-$ will + // be either functions of $W^+$, or they will be prescribed. + // This switch sets them appropriately. Since we are + // using fad variables here, sensitivities will be updated + // appropriately. These sensitivities would be tremendously + // difficult to manage without fad!!! + if (boundary >= 0) { + // Get the boundary descriptor. + typename bdry_map_type::iterator bme = bdry_map.find(boundary); + assert(bme != bdry_map.end()); + + // Evaluate the function object. This is a bit + // tricky; a given boundary might have both prescribed + // and implicit values. If a particular component is not + // prescribed, the values evaluate to zero and are + // ignored, below. + std::vector > bvals(n_q_points, Vector(N_COMP)); + bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals); + + // We loop the quadrature points, and we treat each + // component individualy. + for (int q = 0; q < n_q_points; q++) { + for (int di = 0; di < get_n_components(); di++) { + + // An inflow/dirichlet type of boundary condition + if (bme->second.first[di] == INFLOW_BC) { + Wminus[q][di] = bvals[q](di); + } else if (bme->second.first[di] == PRESSURE_BC) { + // A prescribed pressure boundary condition. This boundary + // condition is complicated by the fact that even though + // the pressure is prescribed, we really are setting + // the energy index here, which will depend on velocity + // and pressure. So even though this seems like a dirichlet + // type boundary condition, we get sensitivities of + // energy to velocity and density (unless these + // are also prescribed. + fad_double rho_vel_sqr = 0; + fad_double dens; + + dens = bme->second.first[DENS_IDX] == INFLOW_BC ? bvals[q](DENS_IDX) : + Wplus[q][DENS_IDX]; + + for (int d=0; d < dim; d++) { + if (bme->second.first[d] == INFLOW_BC) + rho_vel_sqr += bvals[q](d)*bvals[q](d); + else + rho_vel_sqr += Wplus[q][d]*Wplus[q][d]; + } + rho_vel_sqr /= dens; + // Finally set the energy value as determined by the + // prescribed pressure and the other variables. + Wminus[q][di] = bvals[q](di)/(GAMMA-1.0) + + 0.5*rho_vel_sqr; + + } else if (bme->second.first[di] == OUTFLOW_BC) { + // A free/outflow boundary, very simple. + Wminus[q][di] = Wplus[q][di]; + + } else { + // We must be at a no-penetration boundary. We + // prescribe the velocity (we are dealing with a + // particular component here so that the average + // of the velocities is orthogonal to the surface + // normal. This creates sensitivies of across + // the velocity components. + fad_double vdotn = 0; + for (int d = 0; d < dim; d++) { + vdotn += Wplus[q][d]*normals[q](d); + } + + Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di); + } + } + } // for q + } // b>= 0 + + // Determine the Lax-Friedrich's stability parameter, + // and evaluate the numerical flux function at the quadrature points + std::vector > nflux(n_q_points, std::vector(get_n_components(), 0)); + double alpha = 1; + + switch(flux_params.LF_stab) { + case flux_params_type::CONSTANT: + alpha = flux_params.LF_stab_value; + break; + case flux_params_type::MESH: + alpha = face_diameter/(2.0*dT); + break; + } + + LFNumFlux(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus, + alpha); + + // Now assemble the face term + for (unsigned int i=0; iSumIntoGlobalValues(dofs[i], + dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + if (boundary < 0) { + Matrix->SumIntoGlobalValues(dofs[i], + dofs_per_cell, &values[dofs_per_cell], reinterpret_cast(&dofs_neighbor[0])); + } + + // And add into the residual + right_hand_side(dofs[i]) -= F_i.val(); + } + +} + // @sect4{Assembling the whole system} + // Now we put all of the assembly pieces together + // in a routine that dispatches the correct + // piece for each cell/face. We keep track of + // the norm of the resdual for the Newton iteration. +template +void ConsLaw::assemble_system (double &res_norm) +{ + FESystem &fe = *fe_ptr; + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + + // We track the dofs on this cell and (if necessary) + // the adjacent cell. + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); + + // First we create the + // ``UpdateFlags'' for the + // ``FEValues'' and the + // ``FEFaceValues'' objects. + UpdateFlags update_flags = update_values + | update_gradients + | update_q_points + | update_JxW_values; + + // Note, that on faces we do not + // need gradients but we need + // normal vectors. + UpdateFlags face_update_flags = update_values + | update_q_points + | update_JxW_values + | update_normal_vectors; + + // On the neighboring cell we only + // need the shape values. Given a + // specific face, the quadrature + // points and `JxW values' are the + // same as for the current cells, + // the normal vectors are known to + // be the negative of the normal + // vectors of the current cell. + UpdateFlags neighbor_face_update_flags = update_values; + + // Then we create the ``FEValues'' + // object. Note, that since version + // 3.2.0 of deal.II the constructor + // of this class takes a + // ``Mapping'' object as first + // argument. Although the + // constructor without ``Mapping'' + // argument is still supported it + // is recommended to use the new + // constructor. This reduces the + // effect of `hidden magic' (the + // old constructor implicitely + // assumes a ``MappingQ1'' mapping) + // and makes it easier to change + // the mapping object later. + FEValues fe_v ( + mapping, fe, quadrature, update_flags); + + // Similarly we create the + // ``FEFaceValues'' and + // ``FESubfaceValues'' objects for + // both, the current and the + // neighboring cell. Within the + // following nested loop over all + // cells and all faces of the cell + // they will be reinited to the + // current cell and the face (and + // subface) number. + FEFaceValues fe_v_face ( + mapping, fe, face_quadrature, face_update_flags); + FESubfaceValues fe_v_subface ( + mapping, fe, face_quadrature, face_update_flags); + FEFaceValues fe_v_face_neighbor ( + mapping, fe, face_quadrature, neighbor_face_update_flags); + FESubfaceValues fe_v_subface_neighbor ( + mapping, fe, face_quadrature, neighbor_face_update_flags); + + // Furthermore we need some cell + // iterators. + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + // Now we start the loop over all + // active cells. + int fdofs_per_cell = fe_v.dofs_per_cell; + int fn_q_points = face_quadrature.n_quadrature_points; + + unsigned int cell_no = 0; + for (;cell!=endc; ++cell, ++cell_no) + { + + // Now we reinit the ``FEValues'' + // object for the current cell + fe_v.reinit (cell); + + // Collect the local dofs and + // asssemble the cell term. + cell->get_dof_indices (dofs); + + cell_diameter = cell->diameter(); + + assemble_cell_term(fe_v, + dofs, + cell_no); + + // We use the DG style loop through faces + // to determine if we need to apply a + // 'hanging node' flux calculation or a boundary + // computation. + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + // First we set the face + // iterator + typename DoFHandler::face_iterator face=cell->face(face_no); + face_diameter = face->diameter(); + + if (face->at_boundary()) + { + // We reinit the + // ``FEFaceValues'' + // object to the + // current face + fe_v_face.reinit (cell, face_no); + + // and assemble the + // corresponding face + // terms. We send the same + // fe_v and dofs as described + // in the assembly routine. + assemble_face_term( + face_no, fe_v_face, + fe_v_face, + dofs, + dofs, + face->boundary_indicator()); + } + else + { + // Now we are not on + // the boundary of the + // domain, therefore + // there must exist a + // neighboring cell. + typename DoFHandler::cell_iterator neighbor= + cell->neighbor(face_no);; + + if (face->has_children()) + { + // case I: This cell refined compared to neighbor + + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face_no); + + + // We loop over + // subfaces + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; + ++subface_no) + { + typename DoFHandler::active_cell_iterator + neighbor_child + = cell->neighbor_child_on_subface (face_no, subface_no); + + face_diameter = neighbor_child->diameter(); // working on subface + + Assert (neighbor_child->face(neighbor2) == face->child(subface_no), + ExcInternalError()); + Assert (!neighbor_child->has_children(), ExcInternalError()); + + fe_v_subface.reinit (cell, face_no, subface_no); + fe_v_face_neighbor.reinit (neighbor_child, neighbor2); + neighbor_child->get_dof_indices (dofs_neighbor); + + // Assemble as if we are working with + // a DG element. + assemble_face_term( + face_no, fe_v_subface, + fe_v_face_neighbor, + dofs, + dofs_neighbor); + + } + // End of ``if + // (face->has_children())'' + } + else + { + // We have no children, but + // the neighbor cell may be refine + // compared to use + neighbor->get_dof_indices (dofs_neighbor); + if (neighbor->level() != cell->level()) + { + // case II: This is refined compared to neighbor + Assert(neighbor->level() < cell->level(), ExcInternalError()); + const std::pair faceno_subfaceno= + cell->neighbor_of_coarser_neighbor(face_no); + const unsigned int neighbor_face_no=faceno_subfaceno.first, + neighbor_subface_no=faceno_subfaceno.second; + + Assert (neighbor->neighbor_child_on_subface (neighbor_face_no, + neighbor_subface_no) + == cell, + ExcInternalError()); + + // Reinit the + // appropriate + // ``FEFaceValues'' + // and assemble + // the face + // terms. + fe_v_face.reinit (cell, face_no); + fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no, + neighbor_subface_no); + + assemble_face_term( + face_no, fe_v_face, + fe_v_subface_neighbor, + dofs, + dofs_neighbor); + + } + + } + // End of ``face not at boundary'': + } + // End of loop over all faces: + } + + // End iteration through cells. + } + + // Notify Epetra that the matrix is done. + Matrix->FillComplete(); + + // Compute the nonlinear residual. + res_norm = right_hand_side.l2_norm(); + +} + + // Create a conservation law with some defaults. +template +ConsLaw::ConsLaw () + : + mapping (), + fe_ptr(NULL), + dof_handler (triangulation), + quadrature (2), + face_quadrature (2), + T(0), + dT(0.05), + TF(10), + Map(NULL), + Matrix(NULL), + is_stationary(false), + theta(0.5) +{} + + // At one time this example could work for both DG and + // continuous finite elements. The choice was made here. +template +void ConsLaw::build_fe() { + fe_ptr = new FESystem(FE_Q(1), N_COMP); +} + + // Bye bye Conservation law. +template +ConsLaw::~ConsLaw () +{ + dof_handler.clear (); + delete fe_ptr; +} + + // @sect3{Initialize System} + // Sizes all of the vectors and sets up the + // sparsity patter. This function is called at + // the very beginning of a simulation. The function + // setup_system repeats some of these + // chores and is called after adaptivity in leiu + // of this function. +template +void ConsLaw::initialize_system () +{ + // First we need to distribute the + // DoFs. + dof_handler.clear(); + dof_handler.distribute_dofs (*fe_ptr); + + // Size all of the fields. + solution.reinit (dof_handler.n_dofs()); + nlsolution.reinit (dof_handler.n_dofs()); + predictor.reinit (dof_handler.n_dofs()); + ppsolution.reinit (dof_handler.n_dofs()); + dsolution.reinit (dof_handler.n_dofs()); + right_hand_side.reinit (dof_handler.n_dofs()); + indicator.reinit(triangulation.n_active_cells()); +} + + // @sect3{Setup System} + // We call this function to build the sparsity + // and the matrix. +template +void ConsLaw::setup_system () +{ + + // The DoFs of a cell are coupled + // with all DoFs of all neighboring + // cells. Therefore the maximum + // number of matrix entries per row + // is needed when all neighbors of + // a cell are once more refined + // than the cell under + // consideration. + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + (GeometryInfo::faces_per_cell + *GeometryInfo::subfaces_per_face+1)*fe_ptr->dofs_per_cell); + + // Since the continuous sparsity pattern is + // a subset of the DG one, and since we need + // the DG terms for handling hanging nodes, we use + // the flux pattern. + DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern); + + sparsity_pattern.compress(); + + // Rebuild the map. In serial this doesn't do much, + // but is needed. In parallel, this would desribe + // the parallel dof layout. + if (Map) delete Map; + Map = new Epetra_Map(dof_handler.n_dofs(), 0, *Comm); + + // Epetra can build a more efficient matrix if + // one knows ahead of time the maxiumum number of + // columns in any row entry. We traverse the sparsity + // to discover this. + int cur_row = 0; + int cur_col = 0; + int max_df = -1; + for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); + s_i != sparsity_pattern.end(); s_i++) { + if (s_i->row() != cur_row) { + cur_col = 0; + cur_row = s_i->row(); + } + cur_col++; + if (cur_col >= max_df) max_df = cur_col; + } + + if (cur_col >= max_df) max_df = cur_col; + std::cout << "max_df:" << max_df << std::endl; + + // Now we build the matrix, using the constructor + // that optimizes with the max_df variable. + if (Matrix) delete Matrix; + Matrix = new Epetra_CrsMatrix(Copy, *Map, max_df+1, true); + + // We add the sparsity pattern to the matrix by + // inserting zeros. + std::vector vals(max_df, 0); + std::vector row_indices(max_df); + + cur_row = 0; + cur_col = 0; + for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); + s_i != sparsity_pattern.end(); s_i++) { + if (s_i->row() != cur_row) { + Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]); + cur_col = 0; + cur_row = s_i->row(); + } + row_indices[cur_col++] = s_i->column(); + } + // The last row. + Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]); + + // Epetra requires this function after building or + // filling a matrix. It typically does some parallel + // bookeeping; perhaps more. + Matrix->FillComplete(); + +} + + // @sect3{Solving the linear system} + // Actually solve the linear system, using either + // Aztec of Amesos. +template +void ConsLaw::solve (Vector &dsolution, int &niter, double &lin_residual) +{ + + // We must hand the solvers Epetra vectors. + // Luckily, they support the concept of a + // 'view', so we just send in a pointer to our + // dealii vectors. + Epetra_Vector x(View, *Map, dsolution.begin()); + Epetra_Vector b(View, *Map, right_hand_side.begin()); + + // The Direct option selects the Amesos solver. + if (solver_params.SOLVER == solver_params_type::DIRECT) { + + // Setup for solving with + // Amesos. + Epetra_LinearProblem prob; + prob.SetOperator(Matrix); + Amesos_BaseSolver *solver; + Amesos Factory; + + // Other solvers are available + // and may be selected by changing this + // string. + char *stype = "Amesos_Klu"; + + solver = Factory.Create(stype, prob); + + Assert (solver != NULL, ExcInternalError()); + + // There are two parts to the direct solve. + // As I understand, the symbolic part figures + // out the sparsity patterns, and then the + // numerical part actually performs Gaussian + // elimination or whatever the approach is. + if (solver_params.OUTPUT == solver_params_type::VERBOSE) + std::cout << "Starting Symbolic fact\n" << std::flush; + + solver->SymbolicFactorization(); + + if (solver_params.OUTPUT == solver_params_type::VERBOSE) + std::cout << "Starting Numeric fact\n" << std::flush; + + solver->NumericFactorization(); + + + // Define the linear problem by setting the + // right hand and left hand sides. + prob.SetRHS(&b); + prob.SetLHS(&x); + // And finally solve the problem. + if (solver_params.OUTPUT == solver_params_type::VERBOSE) + std::cout << "Starting solve\n" << std::flush; + solver->Solve(); + niter = 0; + lin_residual = 0; + + // We must free the solver that was created + // for us. + delete solver; + + } else if (solver_params.SOLVER == solver_params_type::GMRES) { + + // For the iterative solvers, we use Aztec. + AztecOO Solver; + + // Select the appropriate level of verbosity. + if (solver_params.OUTPUT == solver_params_type::QUIET) + Solver.SetAztecOption(AZ_output, AZ_none); + + if (solver_params.OUTPUT == solver_params_type::VERBOSE) + Solver.SetAztecOption(AZ_output, AZ_all); + + // Select gmres. Other solvers are available. + Solver.SetAztecOption(AZ_solver, AZ_gmres); + Solver.SetRHS(&b); + Solver.SetLHS(&x); + + // Set up the ILUT preconditioner. I do not know + // why, but we must pretend like we are in parallel + // using domain decomposition or the preconditioner + // refuses to activate. + Solver.SetAztecOption(AZ_precond, AZ_dom_decomp); + Solver.SetAztecOption(AZ_subdomain_solve, AZ_ilut); + Solver.SetAztecOption(AZ_overlap, 0); + Solver.SetAztecOption(AZ_reorder, 0); + + // ILUT parameters as described above. + Solver.SetAztecParam(AZ_drop, solver_params.ILUT_DROP); + Solver.SetAztecParam(AZ_ilut_fill, solver_params.ILUT_FILL); + Solver.SetAztecParam(AZ_athresh, solver_params.ILUT_ATOL); + Solver.SetAztecParam(AZ_rthresh, solver_params.ILUT_RTOL); + Solver.SetUserMatrix(Matrix); + + // Run the solver iteration. Collect the number + // of iterations and the residual. + Solver.Iterate(solver_params.MAX_ITERS, solver_params.RES); + niter = Solver.NumIters(); + lin_residual = Solver.TrueResidual(); + } +} + + // @sect3{Postprocessing and Output} + // Recover the physical variables from the conservative + // variables so that output will be (perhaps) more + // meaningfull. +template +void ConsLaw::postprocess() { + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + UpdateFlags update_flags = update_values + | update_gradients + | update_q_points + | update_JxW_values; + UpdateFlags update_flags1 = update_values + | update_gradients + | update_q_points + | update_JxW_values; + + QGauss quadrature_formula(4); + + const std::vector > &us = fe_ptr->base_element(0).get_unit_support_points(); + + + Quadrature unit_support(us); + + int n_q_points = quadrature_formula.n_quadrature_points; + int n_uq_points = unit_support.n_quadrature_points; + + FEValues fe_v ( + mapping, *fe_ptr, quadrature_formula, update_flags); + + FEValues fe_v_unit ( + mapping, *fe_ptr, unit_support, update_flags1); + + std::vector > U(n_uq_points, + Vector(get_n_components())); + std::vector > UU(n_q_points, + Vector(get_n_components())); + std::vector > > dU(n_uq_points, + std::vector >(get_n_components())); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + // Loop the cells + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) { + cell->get_dof_indices (dofs); + fe_v_unit.reinit(cell); + fe_v.reinit(cell); + + fe_v_unit.get_function_values(solution, U); + fe_v_unit.get_function_grads(solution, dU); + fe_v.get_function_values(solution, UU); + + const std::vector &JxW = fe_v.get_JxW_values (); + + for (int q = 0; q < fe_v.get_fe().base_element(0).n_dofs_per_cell(); q++) { + unsigned int didx = fe_v.get_fe().component_to_system_index(DENS_IDX, q); + unsigned int eidx = fe_v.get_fe().component_to_system_index(ENERGY_IDX, q); + double rho_normVsqr = 0; + for (int d = 0; d < dim; d++) { + unsigned int vidx = fe_v.get_fe().component_to_system_index(d, q); + ppsolution(dofs[vidx]) = solution(dofs[vidx])/solution(dofs[didx]); + rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]); + } + rho_normVsqr /= solution(dofs[didx]); + // Pressure + ppsolution(dofs[eidx]) = (GAMMA-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr); + + // Either output density or gradient squared of density, + // depending on what the user wants. + if (!schlieren_plot) { + ppsolution(dofs[didx]) = solution(dofs[didx]); + } else { + double ng = 0; + for (int i = 0; i < dim; i++) ng += dU[q][DENS_IDX][i]*dU[q][DENS_IDX][i]; + ng = std::sqrt(ng); + ppsolution(dofs[didx]) = ng; + } + } + + } // cell + +} + + // Loop and assign a value for refinement. We + // simply use the density squared, which selects + // shocks with some success. +template +void ConsLaw::estimate() { + + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + UpdateFlags update_flags = update_values + | update_gradients + | update_q_points + | update_JxW_values; + + QGauss quadrature_formula(1); + int n_q_points = quadrature_formula.n_quadrature_points; + + + FEValues fe_v ( + mapping, *fe_ptr, quadrature_formula, update_flags); + + std::vector > U(n_q_points, + Vector(get_n_components())); + std::vector > > dU(n_q_points, + std::vector >(get_n_components())); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) { + fe_v.reinit(cell); + + fe_v.get_function_values(predictor, U); + fe_v.get_function_grads(predictor, dU); + + indicator(cell_no) = 0; + for (int q = 0; q < n_q_points; q++) { + double ng = 0; + for (int d = 0; d < dim; d++) ng += dU[q][DENS_IDX][d]*dU[q][DENS_IDX][d]; + + indicator(cell_no) += std::log(1+std::sqrt(ng)); + + } + indicator(cell_no) /= n_q_points; + + } +} + +template +void ConsLaw::refine_grid () +{ + + SolutionTransfer soltrans(dof_handler); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + // Loop cells. If the indicator + // for the cell matches the refinement criterion, + // refine, else unrefine. The unrefinement has + // a slight hysterisis to avoid 'flashing' from refined + // to unrefined. + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) { + cell->clear_coarsen_flag(); + cell->clear_refine_flag(); + if (cell->level() < refinement_params.shock_levels && + std::fabs(indicator(cell_no)) > refinement_params.shock_val ) { + cell->set_refine_flag(); + } else { + if (cell->level() > 0 && + std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val) + cell->set_coarsen_flag(); + } + } + + // The following code prolongs the solution + // to the new grid and carries out the refinement. + std::vector > interp_in; + std::vector > interp_out; + + interp_in.push_back(solution); + interp_in.push_back(predictor); + + triangulation.prepare_coarsening_and_refinement(); + soltrans.prepare_for_coarsening_and_refinement(interp_in); + + triangulation.execute_coarsening_and_refinement (); + + dof_handler.clear(); + dof_handler.distribute_dofs (*fe_ptr); + + { + Vector new_solution(1); + Vector new_predictor(1); + + interp_out.push_back(new_solution); + interp_out.push_back(new_predictor); + interp_out[0].reinit(dof_handler.n_dofs()); + interp_out[1].reinit(dof_handler.n_dofs()); + } + + soltrans.interpolate(interp_in, interp_out); + + // Let the vector delete a very small vector + solution.reinit(1); + predictor.reinit(1); + solution.swap(interp_out[0]); + predictor.swap(interp_out[1]); + + // resize these vectors for the new grid. + nlsolution.reinit(dof_handler.n_dofs()); + ppsolution.reinit(dof_handler.n_dofs()); + nlsolution = solution; + dsolution.reinit (dof_handler.n_dofs()); + right_hand_side.reinit (dof_handler.n_dofs()); + + indicator.reinit(triangulation.n_active_cells()); + +} + +template +void ConsLaw::output_results (const unsigned int cycle) const +{ + char filename[512]; + std::sprintf(filename, "solution-%03d.vtk", cycle); + std::ofstream output (filename); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + std::vector solution_names; + + // Rename the output with the physical variable + // names. Send the post-processed values. + solution_names.clear(); + for (int i = 0; i < dim; i++) { + char buf[512]; + std::sprintf(buf, "v_%d", i); + solution_names.push_back (buf); + } + solution_names.push_back("density"); + solution_names.push_back("pressure"); + data_out.add_data_vector (ppsolution, solution_names); + + data_out.add_data_vector (indicator, "error"); + data_out.build_patches (); + data_out.write_vtk (output); + + output.close(); +} + + // @sect3{Parsing the Input Deck} + // Declare the parameters for the + // input deck. We assume a certain + // maximum number of boundaries and process + // any boundary the user supplies up to + // that maximum number. We + // leave a detailed explanation of these + // parameters to our description of the input + // sample file. +const UInt MAX_BD = 10; +template +void ConsLaw::declare_parameters() { + + // Global scope parameters/ + prm.declare_entry("mesh", "grid.inp", + Patterns::Anything(), + "intput file"); + + prm.declare_entry("diffusion power", "2.0", + Patterns::Double(), + "power of mesh size for diffusion"); + + prm.declare_entry("gravity", "0.0", + Patterns::Double(), + "gravity forcing"); + + // Time stepping block + prm.enter_subsection("time stepping"); + prm.declare_entry("time step", "0.1", + Patterns::Double(), + "simulation time step"); + prm.declare_entry("final time", "10.0", + Patterns::Double(), + "simulation end time"); + prm.leave_subsection(); + + + // Declare the boundary parameters + for (int b = 0; b < MAX_BD; b++) { + char bd[512]; + std::sprintf(bd, "boundary_%d", b); + prm.enter_subsection(bd); + prm.declare_entry("no penetration", "false", + Patterns::Selection("true|false"), + ""); + // declare a slot for each of the conservative + // variables. + for (int di = 0; di < N_COMP; di++) { + char var[512]; + std::sprintf(var, "w_%d", di); + prm.declare_entry(var, "outflow", + Patterns::Selection( + "inflow|outflow|pressure"), + ""); + + // for dirichlet, a function in x,y,z + std::sprintf(var, "w_%d value", di); + prm.declare_entry(var, "0.0", + Patterns::Anything(), + "expression in x,y,z"); + } + + prm.leave_subsection(); + } + + // Initial condition block. + prm.enter_subsection("initial condition"); + for (int di = 0; di < N_COMP; di++) { + char var[512]; + std::sprintf(var, "w_%d", di); + + // for dirichlet, a function in x,y,z + std::sprintf(var, "w_%d value", di); + prm.declare_entry(var, "0.0", + Patterns::Anything(), + "expression in x,y,z"); + } + prm.leave_subsection(); + + // The linear solver block. + prm.enter_subsection("linear solver"); + prm.declare_entry("output", "quiet", + Patterns::Selection( + "quiet|verbose"), + ""); + prm.declare_entry("method", "gmres", + Patterns::Selection( + "gmres|direct"), + ""); + prm.declare_entry("residual", "1e-10", + Patterns::Double(), + "linear solver residual"); + prm.declare_entry("max iters", "300", + Patterns::Double(), + "maximum solver iterations"); + prm.declare_entry("ilut fill", "2", + Patterns::Double(), + "ilut preconditioner fill"); + prm.declare_entry("ilut absolute tolerance", "1e-9", + Patterns::Double(), + "ilut preconditioner tolerance"); + prm.declare_entry("ilut relative tolerance", "1.1", + Patterns::Double(), + "rel tol"); + prm.declare_entry("ilut drop tolerance", "1e-10", + Patterns::Double(), + "ilut drop tol"); + prm.leave_subsection(); + + + // A refinement controller block. + prm.enter_subsection("refinement"); + prm.declare_entry("refinement", "none", + Patterns::Selection( + "none|fixed number|shock"), + ""); + prm.declare_entry("refinement fraction", "0.1", + Patterns::Double(), + "Fraction of high refinement"); + prm.declare_entry("unrefinement fraction", "0.1", + Patterns::Double(), + "Fraction of low unrefinement"); + prm.declare_entry("max elements", "1000000", + Patterns::Double(), + "maximum number of elements"); + prm.declare_entry("shock value", "4.0", + Patterns::Double(), + "value for shock indicator"); + prm.declare_entry("shock levels", "3.0", + Patterns::Double(), + "number of shock refinement levels"); + prm.leave_subsection(); + + // Output control. + prm.enter_subsection("output"); + prm.declare_entry("density", "standard", + Patterns::Selection( + "standard|schlieren"), + ""); + prm.declare_entry("step", "-1", + Patterns::Double(), + "output once per this period"); + prm.leave_subsection(); + + // Flux control + prm.enter_subsection("flux"); + prm.declare_entry("stab", "alpha", + Patterns::Selection( + "alpha|constant|mesh"), + ""); + prm.declare_entry("stab value", "1", + Patterns::Double(), + "alpha stabilization"); + prm.leave_subsection(); + + +} + + // Code to actually parse an input file. This function + // matches the declarations above. +template +void ConsLaw::load_parameters(const char *infile){ + + prm.read_input(infile); + + // The global parameters. + mesh = prm.get("mesh"); + + diffusion_power = prm.get_double("diffusion power"); + + gravity = prm.get_double("gravity"); + + // The time stepping. + prm.enter_subsection("time stepping"); + dT = prm.get_double("time step"); + std::cout << "dT=" << dT << std::endl; + if (dT == 0) { + is_stationary = true; + dT = 1.0; + TF = 1.0; + std::cout << "Stationary mode" << std::endl; + } + TF = prm.get_double("final time"); + std::cout << "TF=" << TF << std::endl; + prm.leave_subsection(); + + // The boundary info + for (int b = 0; b < MAX_BD; b++) { + std::vector flags(N_COMP, OUTFLOW_BC); + + // Define a parser for every boundary, though it may be + // unused. + SideCondition *sd = new SideCondition(N_COMP); + char bd[512]; + std::sprintf(bd, "boundary_%d", b); + prm.enter_subsection(bd); + + const std::string &nopen = prm.get("no penetration"); + + // Determine how each component is handled. + for (int di = 0; di < N_COMP; di++) { + char var[512]; + std::sprintf(var, "w_%d", di); + std::string btype = prm.get(var); + std::sprintf(var, "w_%d value", di); + std::string var_value = prm.get(var); + + if (di < dim && nopen == "true") { + flags[di] = NO_PENETRATION_BC; + } else if (btype == "inflow") { + flags[di] = INFLOW_BC; + sd->set_coeff_row(di, var_value); + } else if (btype == "pressure") { + flags[di] = PRESSURE_BC; + sd->set_coeff_row(di, var_value); + } + } + prm.leave_subsection(); + + // Add the boundary condition to the law. + sd->Init(); + add_boundary(b, flags, sd); + } + + // Initial conditions. + prm.enter_subsection("initial condition"); + for (int di = 0; di < N_COMP; di++) { + char var[512]; + + std::sprintf(var, "w_%d value", di); + std::string var_value = prm.get(var); + ic.set_ic(di, var_value); + } + ic.Init(); + prm.leave_subsection(); + + // The linear solver. + prm.enter_subsection("linear solver"); + const std::string &op = prm.get("output"); + if (op == "verbose") solver_params.OUTPUT = solver_params_type::VERBOSE; + if (op == "quiet") solver_params.OUTPUT = solver_params_type::QUIET; + const std::string &sv = prm.get("method"); + if (sv == "direct") { + solver_params.SOLVER = solver_params_type::DIRECT; + } else if (sv == "gmres") { + solver_params.SOLVER = solver_params_type::GMRES; + } + + solver_params.RES = prm.get_double("residual"); + solver_params.MAX_ITERS = (int) prm.get_double("max iters"); + solver_params.ILUT_FILL = prm.get_double("ilut fill"); + solver_params.ILUT_ATOL = prm.get_double("ilut absolute tolerance"); + solver_params.ILUT_RTOL = prm.get_double("ilut relative tolerance"); + solver_params.ILUT_DROP = prm.get_double("ilut drop tolerance"); + solver_params.RES = prm.get_double("residual"); + prm.leave_subsection(); + + + // And refiement. + prm.enter_subsection("refinement"); + const std::string &ref = prm.get("refinement"); + if (ref == "none") { + refinement_params.refine = refinement_params_type::NONE; + } else if (ref == "fixed number") { + refinement_params.refine = refinement_params_type::FIXED_NUMBER; + } else if (ref == "shock") { + refinement_params.refine = refinement_params_type::SHOCK; + } else + refinement_params.high_frac = prm.get_double("refinement fraction"); + refinement_params.high_frac_sav = refinement_params.high_frac; + refinement_params.low_frac = prm.get_double("unrefinement fraction"); + refinement_params.max_cells = prm.get_double("max elements"); + refinement_params.shock_val = prm.get_double("shock value"); + refinement_params.shock_levels = prm.get_double("shock levels"); + prm.leave_subsection(); + + // Output control. + prm.enter_subsection("output"); + const std::string &dens = prm.get("density"); + schlieren_plot = dens == "schlieren" ? true : false; + output_step = prm.get_double("step"); + prm.leave_subsection(); + + // Flux control. + prm.enter_subsection("flux"); + const std::string &stab = prm.get("stab"); + if (stab == "constant") { + flux_params.LF_stab = flux_params_type::CONSTANT; + } else if (stab == "mesh ") { + flux_params.LF_stab = flux_params_type::MESH; + } + flux_params.LF_stab_value = prm.get_double("stab value"); + prm.leave_subsection(); + + +} + +template +void ConsLaw::zero_matrix() { + Matrix->PutScalar(0); Matrix->FillComplete(); +} + + // We use a predictor to try and make adaptivity + // work better. The idea is to try and refine ahead + // of a front, rather than stepping into a coarse + // set of elements and smearing the solution. This + // simple time extrapolator does the job. +template +void ConsLaw::compute_predictor() { + predictor = nlsolution; + predictor.sadd(3/2.0, -1/2.0, solution); +} + + // @sect3{Run the simulation} Contains the initialization, + // the time loop, and the inner Newton iteration. +template +void ConsLaw::run () +{ + + // Open and load the mesh. + GridIn grid_in; + grid_in.attach_triangulation(triangulation); + std::cout << "Opening mesh <" << mesh << ">" << std::endl; + std::ifstream input_file(mesh.c_str(), std::ios::in); + + Assert (infile, + ExcFileNotOpen()); + + grid_in.read_ucd(input_file); + input_file.close(); + + build_fe(); + + unsigned int nstep = 0; + + // Initialize fields and matrices. + initialize_system (); + setup_system(); + initialize(); + predictor = solution; + + // Initial refinement. We apply the ic, + // estimate, refine, and repeat until + // happy. + if (refinement_params.refine != refinement_params_type::NONE) + for (int i = 0; i < refinement_params.shock_levels; i++) { + estimate(); + refine_grid(); + setup_system(); + initialize(); + predictor = solution; + } + postprocess(); + output_results (nstep); + + // Determine when we will output next. + double next_output = T + output_step; + + // @sect4{Main time stepping loop} + predictor = solution; + while(T < TF) + { + std::cout << "T=" << T << ", "; + + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + bool nonlin_done = false; + double res_norm; + int lin_iter; + + // Print some relevant information during the + // Newton iteration. + std::cout << "NonLin Res: Lin Iter Lin Res" << std::endl; + std::cout << "______________________________________" << std::endl; + + int max_nonlin = 7; + int nonlin_iter = 0; + double lin_res; + + // @sect5{Newton iteration} + nlsolution = predictor; + while (!nonlin_done) { + lin_iter = 0; + zero_matrix(); + right_hand_side = 0; + assemble_system (res_norm); + // Flash a star to the screen so one can + // know when the assembly has stopped and the linear + // solution is starting. + std::cout << "* " << std::flush; + + // Test against a (hardcoded) nonlinear tolderance. + // Do not solve the linear system at the last step + // (since it would be a waste). + + if (fabs(res_norm) < 1e-10) { + nonlin_done = true; + } else { + // Solve the linear system and update with the + // delta. + dsolution = 0; + solve (dsolution, lin_iter, lin_res); + nlsolution.add(1.0, dsolution); + } + + // Print the residuals. + std::printf("%-16.3e %04d %-5.2e\n", + res_norm, lin_iter, lin_res); + + nonlin_iter++; + } + + // Various post convergence tasks. + compute_predictor(); + + solution = nlsolution; + + + estimate(); + + postprocess(); + + T += dT; + + // Output if it is time. + if (output_step < 0) { + output_results (++nstep); + } else if (T >= next_output) { + output_results (++nstep); + next_output += output_step; + } + + // Refine, if refinement is selected. + if (refinement_params.refine != refinement_params_type::NONE) { + refine_grid(); + setup_system(); + } + } +} + + // The following ``main'' function is + // similar to previous examples and + // need not to be commented on. +int main (int argc, char *argv[]) +{ + + MPI_Init(&argc, &argv); + Comm = new Epetra_MpiComm(MPI_COMM_WORLD); + + if (argc != 2) { + std::cout << "Usage:" << argv[0] << " infile" << std::endl; + std::exit(1); + } + try + { + ConsLaw cons; + cons.declare_parameters(); + cons.load_parameters(argv[1]); + cons.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + + return 0; +} +