+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2001 - 2013 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-// Author: Guido Kanschat
-
-#include <base/quadrature_lib.h>
-#include <lac/vector.h>
-#include <lac/full_matrix.h>
-#include <lac/full_matrix.templates.h>
-#include <lac/solver_richardson.h>
-#include <lac/vector_memory.h>
-#include <grid/tria.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <grid/grid_generator.h>
-#include <fe/fe_q.h>
-#include <fe/fe_dgq.h>
-#include <fe/fe_dgp.h>
-#include <fe/fe_system.h>
-#include <fe/mapping_cartesian.h>
-#include <fe/fe_values.h>
-#include <vector>
-#include <iomanip>
-#include <fstream>
-#include <strstream>
-#include <string>
-
-/*
- * Enter name of the finite element family here.
- */
-
-#define ELNAME "FE_DGP"
-#define PUSH(k) FE_DGP<dim> q ## k (k);\
- elements.push_back (ElPair(&q ## k, "dgp" # k))
-
-
-
-template <int dim>
-void compute_embedding (unsigned int degree,
- Triangulation<dim>& tr_coarse,
- Triangulation<dim>& tr_fine,
- const FiniteElement<dim>& fe_coarse,
- const FiniteElement<dim>& fe_fine,
- const char* name)
-{
- deallog.push(name);
-
- DoFHandler<dim> dof_coarse(tr_coarse);
- dof_coarse.distribute_dofs(fe_coarse);
- DoFHandler<dim> dof_fine(tr_fine);
- dof_fine.distribute_dofs(fe_fine);
-
- FullMatrix<long double> A(dof_fine.n_dofs());
- Vector<long double> f(dof_fine.n_dofs());
- Vector<long double> u(dof_fine.n_dofs());
- vector<vector<vector<double> > >
- result(2<<dim,
- vector<vector<double> >(fe_coarse.dofs_per_cell,
- vector<double>(fe_fine.dofs_per_cell, 0.)));
-
- DoFHandler<dim>::active_cell_iterator coarse_cell
- = dof_coarse.begin_active();
- vector<unsigned int> indices(fe_fine.dofs_per_cell);
-
- MappingCartesian<dim> mapping;
- QGauss<dim> q_fine(degree+1);
-
- FEValues<dim> fine (mapping, fe_fine, q_fine,
- update_q_points
- | update_JxW_values
- | update_values);
-
-
- DoFHandler<dim>::active_cell_iterator c;
-
- for (unsigned int coarse_no=0;coarse_no<fe_coarse.dofs_per_cell;
- ++coarse_no)
- {
- A.clear();
- f.clear();
- u.clear();
-
- for (c = dof_fine.begin_active();
- c != dof_fine.end();
- ++c)
- {
- c->get_dof_indices(indices);
- fine.reinit(c);
- // Build mass matrix and RHS
-
- Quadrature<dim> q_coarse (fine.get_quadrature_points(),
- fine.get_JxW_values());
- FEValues<dim> coarse (mapping, fe_coarse, q_coarse,
- update_values);
- coarse.reinit(coarse_cell);
-
- for (unsigned int k=0;k<fine.n_quadrature_points;++k)
- {
- double dx = fine.JxW(k);
- for (unsigned int i=0;i<fe_fine.dofs_per_cell;++i)
- {
- double v = fine.shape_value(i,k);
- f(indices[i]) += dx *
- v * coarse.shape_value(coarse_no, k);
- for (unsigned int j=0;j<fe_fine.dofs_per_cell;++j)
- A(indices[i],indices[j]) += dx*v*fine.shape_value(j,k);
- }
- }
- }
-// A.print_formatted(cout, 2, false, 4, "~", 9);
- FullMatrix<double> P(A.n());
- P.invert(A);
- SolverControl control (100, 1.e-24, true, false);
- PrimitiveVectorMemory<Vector<long double> > mem;
- SolverRichardson<Vector<long double> > solver(control, mem);
- solver.solve(A,u,f,P);
-
- unsigned int cell_i=0;
- for (c = dof_fine.begin_active();
- c != dof_fine.end();
- ++c, ++cell_i)
- {
- c->get_dof_indices(indices);
- for (unsigned int i=0;i<fe_fine.dofs_per_cell;++i)
- result[cell_i][coarse_no][i] = (fabs(u(indices[i])) > 1.e-16)
- ? (27.*u(indices[i])) : 0.;
- }
- }
-
- for (unsigned int cell=0;cell<tr_fine.n_active_cells();++cell)
- {
- cout << "static const double "
- << name << "_"
- << cell << "[] =" << endl << '{' << endl;
- for (unsigned int i=0;i<fe_fine.dofs_per_cell;++i)
- {
- for (unsigned int j=0;j<fe_coarse.dofs_per_cell;++j)
- {
- double s = result[cell][j][i];
- if (fabs(s) < 1.e-13)
- cout << " 0,";
- else
- cout << ' ' << setprecision(10) << s << "/27.,";
- }
- cout << endl;
- }
- cout << "};" << endl << endl;
- }
- deallog.pop();
-}
-
-
-template <int dim>
-void loop ()
-{
- char prefix[3];
- sprintf(prefix, "%dd", dim);
- deallog.push(prefix);
-
- cout << "namespace " << ELNAME << "_" << dim << "d\n{";
-
- Triangulation<dim> tr_coarse;
- Triangulation<dim> tr_fine;
- GridGenerator::hyper_cube (tr_coarse, 0, 1);
- GridGenerator::hyper_cube (tr_fine, 0, 1);
- tr_fine.refine_global(1);
-
- typedef pair<const FiniteElement<dim>*, const char*> ElPair ;
- vector <ElPair> elements(0);
-
- /*
- * List element degrees for
- * computation here.
- */
- PUSH(0);
- PUSH(1);
- PUSH(2);
- PUSH(3);
- PUSH(4);
- PUSH(5);
- PUSH(6);
-
- char* name = new char[100];
-
- // Embed all lower spaces into
- // higher or just the same degree
- // on different grids.
- bool same_degree_only = true;
-
- unsigned int n = elements.size();
- for (unsigned int i=0;i<n;++i)
- for (unsigned int j=((same_degree_only) ? i : 0);j<=i;++j)
- {
- ostrstream os (name, 99);
- os << elements[i].second << "_into_"
- << elements[j].second << "_refined" << ends;
-
- compute_embedding (i, tr_coarse, tr_fine,
- *elements[i].first,
- *elements[j].first,
- name);
- }
-
- delete [] name;
- cout << "};\n";
- deallog.pop();
-
-}
-
-int main ()
-{
- loop<1> ();
- loop<2> ();
- loop<3> ();
-}
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2001 - 2013 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-// Author: Guido Kanschat
-
-// A little program reading a grid *.inp and writing it to *.eps.
-// Some more functionality should be added some time.
-
-#include <grid/tria.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-
-#include <iostream>
-#include <fstream>
-#include <string>
-
-#include <unistd.h>
-
-using namespace std;
-
-template <int dim>
-void convert(const char* infile,
- const char* outfile,
- GridOut::OutputFormat oformat)
-{
- Triangulation<dim> tr;
- GridIn<dim> gin;
- gin.attach_triangulation(tr);
- gin.read(infile);
-
- GridOut gout;
- GridOutFlags::DX dx_flags(true, true, true, false, true);
- gout.set_flags(dx_flags);
-
- ofstream out(outfile);
- gout.write(tr, out, oformat);
-}
-
-
-int main(int argc, char** argv)
-{
- if (argc<4)
- {
- cerr << "Usage: " << argv[0]
- << " dim infile outfile" << endl;
- exit(1);
- }
-
- const unsigned int d = atoi(argv[1]);
-
- std::string outstring(argv[3]);
- GridOut::OutputFormat oformat = GridOut::eps;
-
- const unsigned int dotpos = outstring.find_last_of('.');
- if (dotpos < outstring.length())
- {
- std::string ext = outstring.substr(dotpos+1);
- if (ext == "inp")
- ext = "ucd";
-
- oformat = GridOut::parse_output_format(ext);
- }
-
- if (d == 2)
- convert<2>(argv[2], argv[3], oformat);
- else if (d == 3)
- convert<3>(argv[2], argv[3], oformat);
-}
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2003 - 2013 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-// Author: Guido Kanschat
-
-#include <base/quadrature_lib.h>
-#include <lac/vector.h>
-#include <lac/full_matrix.h>
-#include <lac/full_matrix.templates.h>
-#include <grid/tria.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <grid/grid_generator.h>
-#include <fe/fe_q.h>
-#include <fe/fe_dgq.h>
-#include <fe/fe_system.h>
-#include <fe/mapping_q1.h>
-#include <fe/fe_values.h>
-#include <vector>
-#include <iomanip>
-#include <fstream>
-#include <strstream>
-#include <string>
-
-
-template <int dim>
-inline bool in_cube (const Point<dim>& p)
-{
- if ((p(0) < -1e-10) || (p(0) > 1.000000001))
- return false;
- if ((dim>1) && ((p(1) < -1e-10) || (p(1) > 1.000000001)))
- return false;
- if ((dim>2) && ((p(2) < -1e-10) || (p(2) > 1.000000001)))
- return false;
- return true;
-}
-
-template <int dim>
-inline void
-move_points (vector<Point<dim> >& v, unsigned int i)
-{
- Point<dim> p;
- switch (dim)
- {
- case 1:
- if (i) p(0) = -1.;
- break;
- case 2:
- if ((i==1) || (i==2))
- p(0) = -1.;
- if ((i==2) || (i==3))
- p(1) = -1.;
- break;
- case 3:
- if ((i==1) || (i==2) || (i==5) || (i==6))
- p(0) = -1.;
- if ((i==4) || (i==5) || (i==6) || (i==7))
- p(1) = -1.;
- if ((i==2) || (i==3) || (i==6) || (i==7))
- p(2) = -1.;
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- for (unsigned int j=0;j<v.size();++j)
- v[j] += p;
-}
-
-
-template <int dim>
-void compute_interpolation (Triangulation<dim>& tr_coarse,
- Triangulation<dim>& tr_fine,
- const FiniteElement<dim>& fe_coarse,
- const FiniteElement<dim>& fe_fine,
- const char* name)
-{
- DoFHandler<dim> dof_coarse(tr_coarse);
- dof_coarse.distribute_dofs(fe_coarse);
- DoFHandler<dim> dof_fine(tr_fine);
- dof_fine.distribute_dofs(fe_fine);
-
- vector<vector<vector<double> > >
- result(2<<dim,
- vector<vector<double> >(fe_coarse.dofs_per_cell,
- vector<double>(fe_fine.dofs_per_cell, 0.)));
-
- vector<unsigned int> indices(fe_fine.dofs_per_cell);
-
- MappingQ1<dim> mapping;
-
- vector<Point<dim> > points;
- fe_coarse.get_unit_support_points (points);
- vector<Point<dim> > q_points(points.size());
- vector<double> dummy_weights(points.size(), 0.);
-
- // Coarse level cell of fine grid
- DoFHandler<dim>::cell_iterator father = dof_fine.begin (0);
-
- for (unsigned int child_no = 0; child_no < GeometryInfo<dim>::children_per_cell;
- ++child_no)
- {
- DoFHandler<dim>::active_cell_iterator c = father->child(child_no);
-
- c->get_dof_indices(indices);
- // Evaluate at support points
- // of father cell.
- for (unsigned int i=0;i<points.size();++i)
- q_points[i] = 2.*points[i];
- move_points (q_points, child_no);
-
- Quadrature<dim> q_fine (q_points, dummy_weights);
-
- FEValues<dim> fine (mapping, fe_fine, q_fine,
- update_values);
- fine.reinit(c);
- // Build mass matrix and RHS
-
- for (unsigned int k=0;k<fine.n_quadrature_points;++k)
- {
-// cerr << k << '\t' << q_points[k];
- if (in_cube(q_points[k]))
- {
- for (unsigned int i=0;i<fe_fine.dofs_per_cell;++i)
- {
- double v = fine.shape_value(i,k);
-// cerr << '[' << i << ' ' << v << ']';
- result[child_no][k][i] = v;
- }
- }
-// cerr << endl;
- }
- }
-
- for (unsigned int cell=0;cell<tr_fine.n_active_cells();++cell)
- {
- cout << "static double[] "
- << name << "_"
- << dim << "d_"
- << cell << " =" << endl << '{' << endl;
- for (unsigned int i=0;i<fe_coarse.dofs_per_cell;++i)
- {
- for (unsigned int j=0;j<fe_fine.dofs_per_cell;++j)
- {
- double a = result[cell][i][j];
- if ((a<1.e-14) && (a>-1.e-14))
- a = 0.;
-// cout << ' ' << setprecision(10) << a << ",";
- if (a != 0.)
- cout << "restriction[" << cell << "]("
- << i << ',' << j << ") = "
- << a << ';' << endl;
- }
- cout << endl;
- }
- cout << "};" << endl << endl;
- }
-}
-
-
-template <int dim>
-void compute_projection (Triangulation<dim>& tr_coarse,
- Triangulation<dim>& tr_fine,
- const FiniteElement<dim>& fe_coarse,
- const FiniteElement<dim>& fe_fine,
- const char* name)
-{
- DoFHandler<dim> dof_coarse(tr_coarse);
- dof_coarse.distribute_dofs(fe_coarse);
- DoFHandler<dim> dof_fine(tr_fine);
- dof_fine.distribute_dofs(fe_fine);
- DoFHandler<dim>::cell_iterator cell_coarse = dof_coarse.begin_active();
-
- vector<vector<vector<double> > >
- result(2<<dim,
- vector<vector<double> >(fe_coarse.dofs_per_cell,
- vector<double>(fe_fine.dofs_per_cell, 0.)));
-
- vector<unsigned int> indices(fe_fine.dofs_per_cell);
-
- MappingQ1<dim> mapping;
- QGauss<dim> quadrature (9);
- FullMatrix<double> mass (fe_coarse.dofs_per_cell);
-
- FEValues<dim> coarse (mapping, fe_coarse, quadrature,
- update_values
- | update_JxW_values);
- coarse.reinit (cell_coarse);
-
- // Build coarse level mass matrix
- for (unsigned int k=0;k<coarse.n_quadrature_points;++k)
- {
- for (unsigned int i=0;i<fe_coarse.dofs_per_cell;++i)
- for (unsigned int j=0;j<fe_coarse.dofs_per_cell;++j)
- mass(i,j) += coarse.JxW(k)
- * coarse.shape_value(i,k) * coarse.shape_value(j,k);
- }
-
- FullMatrix<double> inverse (fe_coarse.dofs_per_cell);
- inverse.invert(mass);
- Vector<double> rhs (fe_coarse.dofs_per_cell);
- Vector<double> u (fe_coarse.dofs_per_cell);
-
- // Coarse level cell of fine grid
- DoFHandler<dim>::cell_iterator father = dof_fine.begin (0);
-
- FEValues<dim> fine (mapping, fe_fine, quadrature,
- update_values
- | update_JxW_values
- | update_q_points);
-
- // loop over all fine shape functions
- for (unsigned int n=0;n<fe_fine.dofs_per_cell;++n)
- {
- rhs = 0.;
- DoFHandler<dim>::active_cell_iterator c = father->child(0);
-
- c->get_dof_indices(indices);
-
- fine.reinit(c);
-
- // Build a quadrature formula
- // for the coarse cell
- Quadrature<dim> q_coarse (fine.get_quadrature_points(),
- fine.get_JxW_values());
- FEValues<dim> coarse (mapping, fe_coarse, q_coarse,
- update_values);
- coarse.reinit(cell_coarse);
- // Build RHS
-
- for (unsigned int k=0;k<fine.n_quadrature_points;++k)
- {
- for (unsigned int i=0;i<fe_coarse.dofs_per_cell;++i)
- {
- double v = coarse.shape_value(i,k);
- double f = fine.shape_value(n,k);
- rhs(i) += fine.JxW(k) * v*f;
- }
- }
- inverse.vmult(u,rhs);
- for (unsigned int i=0;i<fe_coarse.dofs_per_cell;++i)
- result[0][i][n] = u(i);
- }
-
- for (unsigned int cell=0;cell<1;++cell)
- {
- cout << "static double "
- << name << "_"
- << dim << "d[] =" << endl << '{' << endl;
- for (unsigned int i=0;i<fe_coarse.dofs_per_cell;++i)
- {
- for (unsigned int j=0;j<fe_fine.dofs_per_cell;++j)
- {
- double a = result[cell][i][j];
- if ((a<1.e-14) && (a>-1.e-14)) a = 0.;
- cout << ' ' << setprecision(8) << a << ",";
- }
- cout << endl;
- }
- cout << "};" << endl << endl;
- }
-}
-
-#define PUSH_Q(k) FE_Q<dim> q ## k(k);\
- elements.push_back (ElPair(&q ## k, "dgq" # k))
-#define PUSH_DGQ(k) FE_DGQ<dim> dgq ## k(k);\
- dgelements.push_back (ElPair(&dgq ## k, "dgq" # k))
-
-template <int dim>
-void loop ()
-{
- Triangulation<dim> tr_coarse;
- Triangulation<dim> tr_fine;
- GridGenerator::hyper_cube (tr_coarse, 0, 1);
- GridGenerator::hyper_cube (tr_fine, 0, 1);
- tr_fine.refine_global(1);
-
- typedef pair<const FiniteElement<dim>*, const char*> ElPair ;
- vector <ElPair> elements(0);
- vector <ElPair> dgelements(0);
-
- PUSH_Q(1);
- PUSH_Q(2);
- PUSH_Q(3);
- PUSH_Q(4);
- PUSH_DGQ(0);
- PUSH_DGQ(1);
- PUSH_DGQ(2);
- PUSH_DGQ(3);
- PUSH_DGQ(4);
- PUSH_DGQ(5);
- PUSH_DGQ(6);
- PUSH_DGQ(7);
-
- char* name = new char[100];
-
- unsigned int n = elements.size();
- for (unsigned int i=0;i<n;++i)
- for (unsigned int j=i;j<=i;++j)
- {
-// ostrstream os (name, 99);
-// os << "interpolate " << elements[j].second << " refined onto "
-// << elements[i].second << ends;
-
-// compute_interpolation (tr_coarse, tr_fine,
-// *elements[i].first,
-// *elements[j].first,
-// name);
- }
-
- n = dgelements.size();
- for (unsigned int i=0;i<n;++i)
- for (unsigned int j=i;j<=i;++j)
- {
- ostrstream os (name, 99);
- os << "project_" << dgelements[j].second << "_refined_onto_"
- << dgelements[i].second << ends;
-
- compute_projection (tr_coarse, tr_fine,
- *dgelements[i].first,
- *dgelements[j].first,
- name);
- }
-
- delete [] name;
-}
-
-int main ()
-{
- loop<1> ();
- loop<2> ();
- loop<3> ();
-}
+++ /dev/null
-## ---------------------------------------------------------------------
-##
-## Copyright (C) 2001 - 2013 by the deal.II authors
-##
-## This file is part of the deal.II library.
-##
-## The deal.II library is free software; you can use it, redistribute
-## it, and/or modify it under the terms of the GNU Lesser General
-## Public License as published by the Free Software Foundation; either
-## version 2.1 of the License, or (at your option) any later version.
-## The full text of the license can be found in the file LICENSE at
-## the top level of the deal.II distribution.
-##
-## ---------------------------------------------------------------------
-
-#
-# Author: Ralf Hartmann, 2001
-#
-
-#
-# Maple script to compute the coefficients of the LagrangeEquidistant
-# basis functions of degree p. These are used as shape functions for
-# Qp elements. For higher p just change variable p in line 10.
-# Call
-# perl -p -e 's/ *t0 = (.*);\n/ $1/g;' lagrange_txt
-# to get a c-code ready to be copied into the source codes.
-#
-
- p := 10:
-
- n_functions := p+1:
-
- # first compute the support points
- support_points := array(0..n_functions-1):
- for i from 0 to n_functions-1 do
- support_points[i] := i/(n_functions-1):
- od;
-
- poly := array(0..n_functions-1):
-
- for i from 0 to n_functions-1 do
- # note that the interp function wants vectors indexed from
- # one and not from zero.
- values := array(1..n_functions):
- for j from 1 to n_functions do
- values[j] := 0:
- od:
- values[i+1] := 1:
-
- shifted_support_points := array (1..n_functions):
- for j from 1 to n_functions do
- shifted_support_points[j] := support_points[j-1]:
- od:
-
- poly[i] := interp (shifted_support_points, values, x):
- od:
-
- readlib(C):
- writeto(lagrange_output):
- printf(` case %d:\n {\n static const double x%d[%d]=\n {`, p,p,(p+1)*(p+1)):
- a := array(0..n_functions-1, 0..n_functions-1):
- b := array(0..n_functions-1):
- # a[i,j] is the jth coefficient of the ith base function.
- for i from 0 to n_functions-1 do
- for j from 0 to n_functions-1 do
- b[j] := coeff(poly[i], x, j):
- od:
- C(b[0]):
- for j from 1 to n_functions-1 do
- printf(`,`):
- C(b[j]):
- od:
- if (i<n_functions-1) then
- printf(`,`):
- fi:
- od:
- printf(`};\n x=&x%d[0];\n break;\n }\n`, p):
-
-
-
+++ /dev/null
-## ---------------------------------------------------------------------
-##
-## Copyright (C) 2002 - 2013 by the deal.II authors
-##
-## This file is part of the deal.II library.
-##
-## The deal.II library is free software; you can use it, redistribute
-## it, and/or modify it under the terms of the GNU Lesser General
-## Public License as published by the Free Software Foundation; either
-## version 2.1 of the License, or (at your option) any later version.
-## The full text of the license can be found in the file LICENSE at
-## the top level of the deal.II distribution.
-##
-## ---------------------------------------------------------------------
-
-s! -?\d\.\d+e-1[0123456789]/27.! 0.!g;
-s! 0/27.! 0.!g;
-s! (-?)27/27.! ${1}1.!g;
-s! (-?)13\.5/27.! ${1}.5!g;
-s! (-?)20\.25/27.! ${1}.75!g;
-s! (-?)6\.75/27.! ${1}.25!g;
-s! (-?)3\.375/27.! ${1}.125!g;
-s! (-?)10\.125/27.! ${1}.375!g;
-s! (-?)16\.875/27.! ${1}.625!g;
-s! (-?)23\.625/27.! ${1}.875!g;
-s! (-?)1\.6875/27.! ${1}1./16.!g;
-s! (-?)5\.0625/27.! ${1}3./16.!g;
-s! (-?)8\.4375/27.! ${1}5./16.!g;
-s! (-?)11\.8125/27.! ${1}7./16.!g;
-s! (-?)15\.1875/27.! ${1}9./16.!g;
-s! (-?)18\.5625/27.! ${1}11./16.!g;
-s! (-?)21\.9375/27.! ${1}13./16.!g;
-s! (-?)25\.3125/27.! ${1}15./16.!g;
-s! (-?)0\.84375/27.! ${1}1./32.!g;
-s! (-?)2\.53125/27.! ${1}3./32.!g;
-s! (-?)4\.21875/27.! ${1}5./32.!g;
-s! (-?)5\.90625/27.! ${1}7./32.!g;
-s! (-?)7\.59375/27.! ${1}9./32.!g;
-s! (-?)9\.28125/27.! ${1}11./32.!g;
-s! (-?)10\.96875/27.! ${1}13./32.!g;
-s! (-?)12\.65625/27.! ${1}15./32.!g;
-s! (-?)14\.34375/27.! ${1}17./32.!g;
-s! (-?)16\.03125/27.! ${1}19./32.!g;
-s! (-?)17\.71875/27.! ${1}21./32.!g;
-s! (-?)19\.40625/27.! ${1}23./32.!g;
-s! (-?)21\.09375/27.! ${1}25./32.!g;
-s! (-?)22\.78125/27.! ${1}27./32.!g;
-s! (-?)24\.46875/27.! ${1}29./32.!g;
-s! (-?)26\.15625/27.! ${1}31./32.!g;
-s! (-?)0\.421875/27.! ${1}1./64.!g;
-s! (-?)1\.265625/27.! ${1}3./64.!g;
-s! (-?)2\.109375/27.! ${1}5./64.!g;
-s! (-?)2\.953125/27.! ${1}7./64.!g;
-s! (-?)3\.796875/27.! ${1}9./64.!g;
-s! (-?)4\.640625/27.! ${1}11./64.!g;
-s! (-?)5\.484375/27.! ${1}13./64.!g;
-s! (-?)6\.328125/27.! ${1}15./64.!g;
-s! (-?)7\.171875/27.! ${1}17./64.!g;
-s! (-?)8\.015625/27.! ${1}19./64.!g;
-s! (-?)8\.859375/27.! ${1}21./64.!g;
-s! (-?)9\.703125/27.! ${1}23./64.!g;
-s! (-?)10\.546875/27.! ${1}25./64.!g;
-s! (-?)11\.390625/27.! ${1}27./64.!g;
-s! (-?)12\.234375/27.! ${1}29./64.!g;
-s! (-?)13\.078125/27.! ${1}31./64.!g;
-s! (-?)13\.921875/27.! ${1}33./64.!g;
-s! (-?)14\.765625/27.! ${1}35./64.!g;
-s! (-?)15\.609375/27.! ${1}37./64.!g;
-s! (-?)16\.453125/27.! ${1}39./64.!g;
-s! (-?)17\.296875/27.! ${1}41./64.!g;
-s! (-?)18\.140625/27.! ${1}43./64.!g;
-s! (-?)18\.984375/27.! ${1}45./64.!g;
-s! (-?)19\.828125/27.! ${1}47./64.!g;
-s! (-?)20\.671875/27.! ${1}49./64.!g;
-s! (-?)21\.515625/27.! ${1}51./64.!g;
-s! (-?)22\.359375/27.! ${1}53./64.!g;
-s! (-?)23\.203125/27.! ${1}55./64.!g;
-s! (-?)24\.046875/27.! ${1}57./64.!g;
-s! (-?)24\.890625/27.! ${1}59./64.!g;
-s! (-?)25\.734375/27.! ${1}61./64.!g;
-s! (-?)26\.578125/27.! ${1}63./64.!g;
-s! (-?)23\.3826859/27.! ${1}SQRT3/2.!g;
-s! (-?)11\.69134295/27.! ${1}SQRT3/4.!g;
-s! (-?)5\.845671476/27.! ${1}SQRT3/8.!g;
-s! (-?)2\.922835738/27.! ${1}SQRT3/16.!g;
-s! (-?)1\.461417869/27.! ${1}SQRT3/32.!g;
-s! (-?)0\.7307089344/27.! ${1}SQRT3/64.!g;
-s! (-?)4\.384253607/27.! ${1}SQRT3*3./32.!g;
-s! (-?)17\.53701443/27.! ${1}SQRT3*3./8.!g;