From: Peter Munch <peterrmuench@gmail.com>
Date: Thu, 2 Jul 2020 08:19:03 +0000 (+0200)
Subject: Add Simplex::ScalarPolynomial
X-Git-Tag: v9.3.0-rc1~1323^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3fe5311e2274765e382ce961239ae332a808db93;p=dealii.git

Add Simplex::ScalarPolynomial
---

diff --git a/doc/news/changes/minor/20200704Munch b/doc/news/changes/minor/20200704Munch
new file mode 100644
index 0000000000..858dbcb751
--- /dev/null
+++ b/doc/news/changes/minor/20200704Munch
@@ -0,0 +1,4 @@
+New: The new class Simplex::ScalarPolynomial provides polynomials defined on 
+simplices.
+<br>
+(Peter Munch, 2020/07/02)
diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h
new file mode 100644
index 0000000000..a99ceb2c31
--- /dev/null
+++ b/include/deal.II/simplex/polynomials.h
@@ -0,0 +1,167 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_simplex_polynomials_h
+#define dealii_simplex_polynomials_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/scalar_polynomials_base.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * A namespace for functions and classes that provide support for simplex
+ * reference cell entities, i.e., triangles and tetrahedrons.
+ */
+namespace Simplex
+{
+  /**
+   * Polynomials defined on dim-dimensional simplex entities. This class is
+   * basis of Simplex::FE_P.
+   */
+  template <int dim>
+  class ScalarPolynomial : public ScalarPolynomialsBase<dim>
+  {
+    static_assert(dim == 2 || dim == 3, "Dimension not supported!");
+
+  public:
+    /**
+     * Make the dimension available to the outside.
+     */
+    static const unsigned int dimension = dim;
+
+    /*
+     * Constructor taking the polynomial @p degree as input.
+     *
+     * @note Currently, only linear (degree=1) and quadratic polynomials
+     *   (degree=2) are implemented.
+     */
+    ScalarPolynomial(const unsigned int degree);
+
+    /**
+     * @copydoc ScalarPolynomialsBase::evaluate()
+     *
+     * @note Currently, only the vectors @p values and @grads are filled.
+     */
+    void
+    evaluate(const Point<dim> &           unit_point,
+             std::vector<double> &        values,
+             std::vector<Tensor<1, dim>> &grads,
+             std::vector<Tensor<2, dim>> &grad_grads,
+             std::vector<Tensor<3, dim>> &third_derivatives,
+             std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_value()
+     */
+    double
+    compute_value(const unsigned int i, const Point<dim> &p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_derivative()
+     *
+     * @note Currently, only implemented for first derivative.
+     */
+    template <int order>
+    Tensor<order, dim>
+    compute_derivative(const unsigned int i, const Point<dim> &p) const;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
+     */
+    Tensor<1, dim>
+    compute_1st_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_2nd_derivative()
+     *
+     * @note Not implemented yet.
+     */
+    Tensor<2, dim>
+    compute_2nd_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
+     *
+     * @note Not implemented yet.
+     */
+    Tensor<3, dim>
+    compute_3rd_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
+     *
+     * @note Not implemented yet.
+     */
+    Tensor<4, dim>
+    compute_4th_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_grad()
+     *
+     * @note Not implemented yet.
+     */
+    Tensor<1, dim>
+    compute_grad(const unsigned int i, const Point<dim> &p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_grad_grad()
+     *
+     * @note Not implemented yet.
+     */
+    Tensor<2, dim>
+    compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::name()
+     */
+    std::string
+    name() const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::clone()
+     */
+    virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+    clone() const override;
+  };
+
+  // template functions
+  template <int dim>
+  template <int order>
+  Tensor<order, dim>
+  ScalarPolynomial<dim>::compute_derivative(const unsigned int i,
+                                            const Point<dim> & p) const
+  {
+    Tensor<order, dim> der;
+
+    AssertDimension(order, 1);
+    const auto grad = compute_grad(i, p);
+
+    for (unsigned int i = 0; i < dim; i++)
+      der[i] = grad[i];
+
+    return der;
+  }
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt
index cfc0183c55..f0c06e1d0d 100644
--- a/source/simplex/CMakeLists.txt
+++ b/source/simplex/CMakeLists.txt
@@ -16,6 +16,7 @@
 INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 
 SET(_unity_include_src
+  polynomials.cc
   quadrature_lib.cc
   )
 
diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc
new file mode 100644
index 0000000000..c1823cfdff
--- /dev/null
+++ b/source/simplex/polynomials.cc
@@ -0,0 +1,444 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/simplex/polynomials.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+  namespace
+  {
+    unsigned int
+    compute_n_polynomials(const unsigned int dim, const unsigned int degree)
+    {
+      if (dim == 2)
+        {
+          if (degree == 1)
+            return 3;
+          if (degree == 2)
+            return 6;
+        }
+      else if (dim == 3)
+        {
+          if (degree == 1)
+            return 4;
+          if (degree == 2)
+            return 10;
+        }
+
+      Assert(false, ExcNotImplemented());
+
+      return 0;
+    }
+  } // namespace
+
+
+
+  template <int dim>
+  ScalarPolynomial<dim>::ScalarPolynomial(const unsigned int degree)
+    : ScalarPolynomialsBase<dim>(degree, compute_n_polynomials(dim, degree))
+  {}
+
+
+
+  template <int dim>
+  double
+  ScalarPolynomial<dim>::compute_value(const unsigned int i,
+                                       const Point<dim> & p) const
+  {
+    if (dim == 2)
+      {
+        if (this->degree() == 1)
+          {
+            if (i == 0)
+              return p[0];
+            else if (i == 1)
+              return p[1];
+            else if (i == 2)
+              return 1.0 - p[0] - p[1];
+          }
+        else if (this->degree() == 2)
+          {
+            const double t1 = p[0];
+            const double t2 = p[1];
+            const double t3 = 1.0 - p[0] - p[1];
+
+            if (i == 0)
+              return t1 * (2.0 * t1 - 1.0);
+            else if (i == 1)
+              return t2 * (2.0 * t2 - 1.0);
+            else if (i == 2)
+              return t3 * (2.0 * t3 - 1.0);
+            else if (i == 3)
+              return 4.0 * t2 * t1;
+            else if (i == 4)
+              return 4.0 * t2 * t3;
+            else if (i == 5)
+              return 4.0 * t3 * t1;
+          }
+      }
+    else if (dim == 3)
+      {
+        if (this->degree() == 1)
+          {
+            if (i == 0)
+              return 1.0 - p[0] - p[1] - p[2];
+            else if (i == 1)
+              return p[0];
+            else if (i == 2)
+              return p[1];
+            else if (i == 3)
+              return p[2];
+          }
+        else if (this->degree() == 2)
+          {
+            const double r = p[0];
+            const double s = p[1];
+            const double t = p[2];
+            const double u = 1.0 - p[0] - p[1] - p[2];
+            if (i == 0)
+              return u * (2.0 * u - 1.0);
+            else if (i == 1)
+              return r * (2.0 * r - 1.0);
+            else if (i == 2)
+              return s * (2.0 * s - 1.0);
+            else if (i == 3)
+              return t * (2.0 * t - 1.0);
+            else if (i == 4)
+              return 4.0 * r * u;
+            else if (i == 5)
+              return 4.0 * r * s;
+            else if (i == 6)
+              return 4.0 * s * u;
+            else if (i == 7)
+              return 4.0 * t * u;
+            else if (i == 8)
+              return 4.0 * r * t;
+            else if (i == 9)
+              return 4.0 * s * t;
+          }
+      }
+
+    Assert(false, ExcNotImplemented());
+
+    return 0;
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  ScalarPolynomial<dim>::compute_grad(const unsigned int i,
+                                      const Point<dim> & p) const
+  {
+    Tensor<1, dim> grad;
+
+    if (dim == 2)
+      {
+        if (this->degree() == 1)
+          {
+            if (i == 0)
+              {
+                grad[0] = +1.0;
+                grad[1] = +0.0;
+              }
+            else if (i == 1)
+              {
+                grad[0] = +0.0;
+                grad[1] = +1.0;
+              }
+            else if (i == 2)
+              {
+                grad[0] = -1.0;
+                grad[1] = -1.0;
+              }
+            else
+              {
+                Assert(false, ExcNotImplemented());
+              }
+          }
+        else if (this->degree() == 2)
+          {
+            if (i == 2)
+              {
+                grad[0] = -3.0 + 4.0 * (p[0] + p[1]);
+                grad[1] = -3.0 + 4.0 * (p[0] + p[1]);
+              }
+            else if (i == 0)
+              {
+                grad[0] = 4.0 * p[0] - 1.0;
+                grad[1] = 0.0;
+              }
+            else if (i == 1)
+              {
+                grad[0] = 0.0;
+                grad[1] = 4.0 * p[1] - 1.0;
+              }
+            else if (i == 5)
+              {
+                grad[0] = 4.0 * (1.0 - 2.0 * p[0] - p[1]);
+                grad[1] = -4.0 * p[0];
+              }
+            else if (i == 3)
+              {
+                grad[0] = 4.0 * p[1];
+                grad[1] = 4.0 * p[0];
+              }
+            else if (i == 4)
+              {
+                grad[0] = -4.0 * p[1];
+                grad[1] = 4.0 * (1.0 - p[0] - 2.0 * p[1]);
+              }
+            else
+              {
+                Assert(false, ExcNotImplemented());
+              }
+          }
+        else
+          {
+            Assert(false, ExcNotImplemented());
+          }
+      }
+    else if (dim == 3)
+      {
+        if (this->degree() == 1)
+          {
+            if (i == 0)
+              {
+                grad[0] = -1.0;
+                grad[1] = -1.0;
+                grad[2] = -1.0;
+              }
+            else if (i == 1)
+              {
+                grad[0] = +1.0;
+                grad[1] = +0.0;
+                grad[2] = +0.0;
+              }
+            else if (i == 2)
+              {
+                grad[0] = +0.0;
+                grad[1] = +1.0;
+                grad[2] = +0.0;
+              }
+            else if (i == 3)
+              {
+                grad[0] = +0.0;
+                grad[1] = +0.0;
+                grad[2] = +1.0;
+              }
+          }
+        else if (this->degree() == 2)
+          {
+            const double r = p[0];
+            const double s = p[1];
+            const double t = p[2];
+            const double u = 1.0 - p[0] - p[1] - p[2];
+
+            if (i == 0)
+              {
+                grad[0] = -4.0 * u + 1.;
+                grad[1] = grad[0];
+                grad[2] = grad[0];
+              }
+            else if (i == 1)
+              {
+                grad[0] = +4.0 * r - 1.;
+                grad[1] = +0.0;
+                grad[2] = +0.0;
+              }
+            else if (i == 2)
+              {
+                grad[0] = +0.0;
+                grad[1] = +4.0 * s - 1.;
+                grad[2] = +0.0;
+              }
+            else if (i == 3)
+              {
+                grad[0] = +0.0;
+                grad[1] = +0.0;
+                grad[2] = +4.0 * t - 1.;
+              }
+            else if (i == 4)
+              {
+                grad[0] = +4.0 * (u - r);
+                grad[1] = -4.0 * r;
+                grad[2] = -4.0 * r;
+              }
+            else if (i == 5)
+              {
+                grad[0] = +4.0 * s;
+                grad[1] = +4.0 * r;
+                grad[2] = +0.0;
+              }
+            else if (i == 6)
+              {
+                grad[0] = -4.0 * s;
+                grad[1] = +4.0 * (u - s);
+                grad[2] = -4.0 * s;
+              }
+            else if (i == 7)
+              {
+                grad[0] = -4.0 * t;
+                grad[1] = -4.0 * t;
+                grad[2] = +4.0 * (u - t);
+              }
+            else if (i == 8)
+              {
+                grad[0] = +4.0 * t;
+                grad[1] = +0.0;
+                grad[2] = +4.0 * r;
+              }
+            else if (i == 9)
+              {
+                grad[0] = +0.0;
+                grad[1] = +4.0 * t;
+                grad[2] = +4.0 * s;
+              }
+          }
+        else
+          {
+            Assert(false, ExcNotImplemented());
+          }
+      }
+    else
+      {
+        Assert(false, ExcNotImplemented());
+      }
+
+    return grad;
+  }
+
+
+
+  template <int dim>
+  Tensor<2, dim>
+  ScalarPolynomial<dim>::compute_grad_grad(const unsigned int i,
+                                           const Point<dim> & p) const
+  {
+    (void)i;
+    (void)p;
+
+    Assert(false, ExcNotImplemented());
+    return Tensor<2, dim>();
+  }
+
+
+
+  template <int dim>
+  void
+  ScalarPolynomial<dim>::evaluate(
+    const Point<dim> &           unit_point,
+    std::vector<double> &        values,
+    std::vector<Tensor<1, dim>> &grads,
+    std::vector<Tensor<2, dim>> &grad_grads,
+    std::vector<Tensor<3, dim>> &third_derivatives,
+    std::vector<Tensor<4, dim>> &fourth_derivatives) const
+  {
+    (void)grads;
+    (void)grad_grads;
+    (void)third_derivatives;
+    (void)fourth_derivatives;
+
+    if (values.size() == this->n())
+      for (unsigned int i = 0; i < this->n(); i++)
+        values[i] = compute_value(i, unit_point);
+
+    if (grads.size() == this->n())
+      for (unsigned int i = 0; i < this->n(); i++)
+        grads[i] = compute_grad(i, unit_point);
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  ScalarPolynomial<dim>::compute_1st_derivative(const unsigned int i,
+                                                const Point<dim> & p) const
+  {
+    return compute_grad(i, p);
+  }
+
+
+
+  template <int dim>
+  Tensor<2, dim>
+  ScalarPolynomial<dim>::compute_2nd_derivative(const unsigned int i,
+                                                const Point<dim> & p) const
+  {
+    (void)i;
+    (void)p;
+
+    Assert(false, ExcNotImplemented());
+
+    return {};
+  }
+
+
+
+  template <int dim>
+  Tensor<3, dim>
+  ScalarPolynomial<dim>::compute_3rd_derivative(const unsigned int i,
+                                                const Point<dim> & p) const
+  {
+    (void)i;
+    (void)p;
+
+    Assert(false, ExcNotImplemented());
+
+    return {};
+  }
+
+
+
+  template <int dim>
+  Tensor<4, dim>
+  ScalarPolynomial<dim>::compute_4th_derivative(const unsigned int i,
+                                                const Point<dim> & p) const
+  {
+    (void)i;
+    (void)p;
+
+    Assert(false, ExcNotImplemented());
+
+    return {};
+  }
+
+
+
+  template <int dim>
+  std::string
+  ScalarPolynomial<dim>::name() const
+  {
+    return "Simplex";
+  }
+
+
+
+  template <int dim>
+  std::unique_ptr<ScalarPolynomialsBase<dim>>
+  ScalarPolynomial<dim>::clone() const
+  {
+    return std::make_unique<ScalarPolynomial<dim>>(*this);
+  }
+
+  template class ScalarPolynomial<2>;
+  template class ScalarPolynomial<3>;
+
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/simplex/polynomials_01.cc b/tests/simplex/polynomials_01.cc
new file mode 100644
index 0000000000..86c276ef9e
--- /dev/null
+++ b/tests/simplex/polynomials_01.cc
@@ -0,0 +1,86 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test Simplex::ScalarPolynomial on an the points of an arbitrary quadrature
+// rule.
+
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/simplex/polynomials.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+  Simplex::ScalarPolynomial<dim> poly(degree);
+  QSimplex<dim>                  quad(QGauss<dim>(degree + 1));
+
+  std::vector<double>         values(poly.n());
+  std::vector<Tensor<1, dim>> grads(poly.n());
+  std::vector<Tensor<2, dim>> grad_grads;
+  std::vector<Tensor<3, dim>> third_derivatives;
+  std::vector<Tensor<4, dim>> fourth_derivatives;
+
+  for (unsigned int i = 0; i < quad.size(); i++)
+    {
+      poly.evaluate(quad.point(i),
+                    values,
+                    grads,
+                    grad_grads,
+                    third_derivatives,
+                    fourth_derivatives);
+
+      for (auto v : values)
+        deallog << v << " ";
+      deallog << std::endl;
+
+      for (auto v : grads)
+        deallog << v << " ";
+      deallog << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("2d-1");
+    test<2>(1);
+    deallog.pop();
+  }
+  {
+    deallog.push("2d-2");
+    test<2>(2);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-1");
+    test<3>(1);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-2");
+    test<3>(2);
+    deallog.pop();
+  }
+}
diff --git a/tests/simplex/polynomials_01.output b/tests/simplex/polynomials_01.output
new file mode 100644
index 0000000000..fe28ac73ea
--- /dev/null
+++ b/tests/simplex/polynomials_01.output
@@ -0,0 +1,29 @@
+
+DEAL:2d-1::0.211325 0.211325 0.577350 
+DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 
+DEAL:2d-1::0.788675 0.211325 2.77556e-17 
+DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 
+DEAL:2d-1::0.211325 0.788675 0.00000 
+DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 
+DEAL:2d-2::-0.0872983 -0.0872983 0.425403 0.0508067 0.349193 0.349193 
+DEAL:2d-2::-0.549193 0.00000 0.00000 -0.549193 -2.09839 -2.09839 0.450807 0.450807 -0.450807 2.64758 2.64758 -0.450807 
+DEAL:2d-2::0.00000 -0.0872983 -0.0872983 0.225403 0.174597 0.774597 
+DEAL:2d-2::1.00000 0.00000 0.00000 -0.549193 -0.549193 -0.549193 0.450807 2.00000 -0.450807 1.09839 -0.450807 -2.00000 
+DEAL:2d-2::0.687298 -0.0872983 1.38778e-17 0.400000 -6.25620e-18 -4.92550e-17 
+DEAL:2d-2::2.54919 0.00000 0.00000 -0.549193 1.00000 1.00000 0.450807 3.54919 -0.450807 -0.450807 -3.54919 -3.54919 
+DEAL:2d-2::-0.0872983 0.00000 -0.0872983 0.225403 0.774597 0.174597 
+DEAL:2d-2::-0.549193 0.00000 0.00000 1.00000 -0.549193 -0.549193 2.00000 0.450807 -2.00000 -0.450807 1.09839 -0.450807 
+DEAL:2d-2::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:2d-2::1.00000 0.00000 0.00000 1.00000 1.00000 1.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -2.00000 
+DEAL:2d-2::-0.0872983 0.687298 0.00000 0.400000 0.00000 0.00000 
+DEAL:2d-2::-0.549193 0.00000 0.00000 2.54919 1.00000 1.00000 3.54919 0.450807 -3.54919 -3.54919 -0.450807 -0.450807 
+DEAL:3d-1::0.366025 0.211325 0.211325 0.211325 
+DEAL:3d-1::-1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 
+DEAL:3d-2::0.214315 -0.0872983 -0.0872983 -0.0872983 0.298387 0.0508067 0.298387 0.298387 0.0508067 0.0508067 
+DEAL:3d-2::-1.64758 -1.64758 -1.64758 -0.549193 0.00000 0.00000 0.00000 -0.549193 0.00000 0.00000 0.00000 -0.549193 2.19677 -0.450807 -0.450807 0.450807 0.450807 0.00000 -0.450807 2.19677 -0.450807 -0.450807 -0.450807 2.19677 0.450807 0.00000 0.450807 0.00000 0.450807 0.450807 
+DEAL:3d-2::-0.123790 0.00000 -0.0872983 -0.0872983 0.549193 0.225403 0.123790 0.123790 0.225403 0.0508067 
+DEAL:3d-2::-0.0983867 -0.0983867 -0.0983867 1.00000 0.00000 0.00000 0.00000 -0.549193 0.00000 0.00000 0.00000 -0.549193 -0.901613 -2.00000 -2.00000 0.450807 2.00000 0.00000 -0.450807 0.647580 -0.450807 -0.450807 -0.450807 0.647580 0.450807 0.00000 2.00000 0.00000 0.450807 0.450807 
+DEAL:3d-2::-0.123790 -0.0872983 0.00000 -0.0872983 0.123790 0.225403 0.549193 0.123790 0.0508067 0.225403 
+DEAL:3d-2::-0.0983867 -0.0983867 -0.0983867 -0.549193 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 -0.549193 0.647580 -0.450807 -0.450807 2.00000 0.450807 0.00000 -2.00000 -0.901613 -2.00000 -0.450807 -0.450807 0.647580 0.450807 0.00000 0.450807 0.00000 0.450807 2.00000 
+DEAL:3d-2::-0.123790 -0.0872983 -0.0872983 0.00000 0.123790 0.0508067 0.123790 0.549193 0.225403 0.225403 
+DEAL:3d-2::-0.0983867 -0.0983867 -0.0983867 -0.549193 0.00000 0.00000 0.00000 -0.549193 0.00000 0.00000 0.00000 1.00000 0.647580 -0.450807 -0.450807 0.450807 0.450807 0.00000 -0.450807 0.647580 -0.450807 -2.00000 -2.00000 -0.901613 2.00000 0.00000 0.450807 0.00000 2.00000 0.450807