-//---------------------------------------------------------------------------
-// $Id$
-// Version: $Name$
-//
-// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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.
-//
-//---------------------------------------------------------------------------
-
#include <base/logstream.h>
#include <base/utilities.h>
#include <base/quadrature.h>
template <int dim>
void
FE_Nedelec<dim>::get_subface_interpolation_matrix(
- const FiniteElement<dim>&, const unsigned int,
- FullMatrix<double>&) const
+ const FiniteElement<dim>&,
+ const unsigned int,
+ FullMatrix<double>&)
+const
{
Assert (false, ExcNotImplemented ());
}
// interpolation part.
for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
{
- const double weight = 2.0 * edge_quadrature.weight (q_point);
+ const double weight = edge_quadrature.weight (q_point);
for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
interpolation_matrix (0, dof)
const double tmp = 2.0 * quadrature_points[q_point] (0)
+ shift[subface];
const double weight
- = std::sqrt (2.0 * edge_quadrature.weight (q_point));
+ = std::sqrt (edge_quadrature.weight (q_point));
for (unsigned int i = 0; i < deg; ++i)
assembling_matrix (i, q_point)
= 2.0 * quadrature_points[q_point] (0)
+ shift[subface];
const double weight
- = 2.0 * edge_quadrature.weight (q_point)
+ = edge_quadrature.weight (q_point)
* (this->shape_value_component
(dof, Point<dim> (0.0,
quadrature_points[q_point] (0)),
++q_point)
{
const double
- weight = 2.0 * edge_quadrature.weight (q_point);
+ weight = edge_quadrature.weight (q_point);
for (unsigned int i = 0; i < 2; ++i)
for (unsigned int dof = 0; dof < this->dofs_per_face;
// face.
std::vector<std::vector<Point<dim> > >
edge_quadrature_points_full_dim
- (GeometryInfo<dim>::lines_per_face);
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<Point<dim> > (n_edge_points));
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face; ++line)
- edge_quadrature_points_full_dim.resize (n_edge_points);
-
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
0.5, 0.0);
}
- Vector<double> solution (deg);
- Vector<double> system_rhs (deg);
-
- for (unsigned int dof = 0; dof < this->dofs_per_face;
- ++dof)
- {
// Set up the system matrix.
// This can be used for all
// edges.
- for (unsigned int q_point = 0;
- q_point < n_edge_points; ++q_point)
- {
- const double tmp
- = 2.0 * edge_quadrature_points[q_point] (0);
- const double weight
- = std::sqrt (2.0 * edge_quadrature.weight
- (q_point));
-
- for (unsigned int i = 0; i < deg; ++i)
- assembling_matrix (i, q_point)
- = weight
- * lobatto_polynomials_grad[i + 1].value
- (tmp);
- }
+ for (unsigned int q_point = 0;
+ q_point < n_edge_points; ++q_point)
+ {
+ const double tmp
+ = 2.0 * edge_quadrature_points[q_point] (0);
+ const double weight
+ = std::sqrt (edge_quadrature.weight (q_point));
+
+ for (unsigned int i = 0; i < deg; ++i)
+ assembling_matrix (i, q_point)
+ = weight * lobatto_polynomials_grad[i + 1].value
+ (tmp);
+ }
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.invert (system_matrix);
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.invert (system_matrix);
+ Vector<double> solution (deg);
+ Vector<double> system_rhs (deg);
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
for (unsigned int line = 0;
line < GeometryInfo<dim>::lines_per_face;
++line)
q_point < n_edge_points; ++q_point)
{
const double right_hand_side_value
- = std::sqrt (2.0 * edge_quadrature.weight
+ = std::sqrt (edge_quadrature.weight
(q_point))
* (this->shape_value_component
(this->face_to_cell_index (dof, 4),
= solution (i);
}
- assembling_matrix.reinit (deg * this->degree,
- n_face_points);
- system_rhs.reinit (assembling_matrix.m ());
+ assembling_matrix.reinit (deg * this->degree,
+ n_face_points);
+
+ for (unsigned int q_point = 0; q_point < n_face_points;
+ ++q_point)
+ {
+ const Point<dim - 1> quadrature_point
+ (2.0 * face_quadrature_points[q_point] (0),
+ 2.0 * face_quadrature_points[q_point] (1));
+ const double weight
+ = std::sqrt (face_quadrature.weight (q_point));
+
+ for (unsigned int i = 0; i <= deg; ++i)
+ {
+ const double tmp
+ = weight * legendre_polynomials[i].value
+ (quadrature_point (0));
+
+ for (unsigned int j = 0; j < deg; ++j)
+ assembling_matrix (i * deg + j, q_point)
+ = tmp * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
+ }
+ }
+
+ system_matrix.reinit (assembling_matrix.m (),
+ assembling_matrix.m ());
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.reinit (system_matrix.m (),
+ system_matrix.m ());
+ system_matrix_inv.invert (system_matrix);
+ solution.reinit (system_matrix_inv.m ());
+ system_rhs.reinit (system_matrix_inv.m ());
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
+ {
system_rhs = 0;
// Now we project the remaining
quadrature_point, 1);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
+ *= face_quadrature.weight (q_point);
- const double weight
- = std::sqrt (4.0 * face_quadrature.weight
- (q_point));
-
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
- const double tmp1 = weight * L_i;
- const double tmp2
- = right_hand_side_value * L_i;
+ const double tmp
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
- {
- const double l_j
- = lobatto_polynomials[j + 2].value
- (quadrature_point (1));
-
- assembling_matrix (i * deg + j, q_point)
- = tmp1 * l_j;
- system_rhs (i * deg + j) += tmp2 * l_j;
- }
+ system_rhs (i * deg + j)
+ += tmp
+ * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
}
}
- system_matrix.reinit (assembling_matrix.m (),
- assembling_matrix.m ());
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
- system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix_inv.m ());
system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
quadrature_point, 0);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
+ *= face_quadrature.weight (q_point);
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
- const double tmp = right_hand_side_value
- * L_i;
+ const double tmp
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
-// const double weight
-// = 2.0 * edge_quadrature_x.weight (q_point);
+ const double weight
+ = edge_quadrature_x.weight (q_point);
for (unsigned int i = 0; i < 2; ++i)
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
{
interpolation_matrix (i * source_fe.degree, dof)
- += this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- Point<dim>
- (0.5 * (i + 1),
- edge_quadrature_y_points[q_point] (0), 0.0),
- 1);
+ += weight
+ * this->shape_value_component
+ (this->face_to_cell_index (dof, 4),
+ Point<dim>
+ (0.5 * (i + 1),
+ edge_quadrature_y_points[q_point] (0), 0.0),
+ 1);
interpolation_matrix
((i + 2) * source_fe.degree, dof)
- += this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- Point<dim>
- (edge_quadrature_x_points[q_point] (0),
- 0.5 * i, 0.0), 0);
+ += weight
+ * this->shape_value_component
+ (this->face_to_cell_index (dof, 4),
+ Point<dim>
+ (edge_quadrature_x_points[q_point] (0),
+ 0.5 * i, 0.0), 0);
}
}
// weights on the four
// edges of a face.
std::vector<std::vector<double> > edge_quadrature_points
- (GeometryInfo<dim>::lines_per_face);
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<double> (n_edge_points));
std::vector<std::vector<double> >
edge_quadrature_weights
- (GeometryInfo<dim>::lines_per_face);
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<double> (n_edge_points));
std::vector<std::vector<Point<dim> > >
edge_quadrature_points_full_dim
- (GeometryInfo<dim>::lines_per_face);
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face; ++line)
- {
- edge_quadrature_points.resize (n_edge_points);
- edge_quadrature_points_full_dim.resize
- (n_edge_points);
- edge_quadrature_weights.resize (n_edge_points);
- }
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<Point<dim> > (n_edge_points));
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
(edge_quadrature_x_points[q_point] (0), 0.5,
0.0);
edge_quadrature_weights[0][q_point]
- = std::sqrt (2.0 * edge_quadrature_y.weight
- (q_point));
+ = std::sqrt (edge_quadrature_y.weight (q_point));
edge_quadrature_weights[1][q_point]
= edge_quadrature_weights[0][q_point];
edge_quadrature_weights[2][q_point]
- = std::sqrt (2.0 * edge_quadrature_x.weight
- (q_point));
+ = std::sqrt (edge_quadrature_x.weight (q_point));
edge_quadrature_weights[3][q_point]
= edge_quadrature_weights[2][q_point];
}
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- for (unsigned int dof = 0; dof < this->dofs_per_face;
- ++dof)
- {
// Set up the system matrix.
// This can be used for all
// edges.
- for (unsigned int q_point = 0;
- q_point < n_edge_points; ++q_point)
- {
- const double tmp
- = 2.0 * edge_quadrature_y_points[q_point] (0);
- const double weight
- = std::sqrt (2.0 * edge_quadrature_y.weight
- (q_point));
-
- for (unsigned int i = 0; i < deg; ++i)
- assembling_matrix (i, q_point)
- = weight
- * lobatto_polynomials_grad[i + 1].value
- (tmp);
- }
+ for (unsigned int q_point = 0;
+ q_point < n_edge_points; ++q_point)
+ {
+ const double tmp
+ = 2.0 * edge_quadrature_y_points[q_point] (0);
+ const double weight
+ = std::sqrt (edge_quadrature_y.weight (q_point));
+
+ for (unsigned int i = 0; i < deg; ++i)
+ assembling_matrix (i, q_point)
+ = weight
+ * lobatto_polynomials_grad[i + 1].value (tmp);
+ }
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.invert (system_matrix);
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.invert (system_matrix);
+ Vector<double> system_rhs (system_matrix.m ());
+ Vector<double> solution (system_rhs.size ());
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
for (unsigned int line = 0;
line < GeometryInfo<dim - 1>::lines_per_cell;
++line)
// functions. First on the
// horizontal ones, then on
// the vertical ones.
- assembling_matrix.reinit (deg * this->degree,
- n_face_points);
- system_rhs.reinit (assembling_matrix.m ());
+ assembling_matrix.reinit (deg * this->degree,
+ n_face_points);
+
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ const Point<dim - 1> quadrature_point
+ (2.0 * face_quadrature_points[q_point] (0)
+ - 1.0,
+ 2.0 * face_quadrature_points[q_point] (1));
+ const double weight
+ = std::sqrt (face_quadrature.weight (q_point));
+
+ for (unsigned int i = 0; i <= deg; ++i)
+ {
+ const double tmp
+ = weight * legendre_polynomials[i].value
+ (quadrature_point (0));
+
+ for (unsigned int j = 0; j < deg; ++j)
+ assembling_matrix (i * deg + j, q_point)
+ = tmp * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
+ }
+ }
+
+ system_matrix.reinit (assembling_matrix.m (),
+ assembling_matrix.m ());
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.reinit (system_matrix.m (),
+ system_matrix.m ());
+ system_matrix_inv.invert (system_matrix);
+ solution.reinit (system_matrix_inv.m ());
+ system_rhs.reinit (assembling_matrix.m ());
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
+ {
system_rhs = 0;
for (unsigned int q_point = 0;
quadrature_point, 1);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
-
- const double weight
- = std::sqrt (4.0 * face_quadrature.weight
- (q_point));
+ *= face_quadrature.weight (q_point);
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
- const double tmp1 = weight * L_i;
- const double tmp2 = right_hand_side_value
- * L_i;
+ const double tmp
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
- {
- const double l_j
- = lobatto_polynomials[j + 2].value
+ system_rhs (i * deg + j)
+ = tmp
+ * lobatto_polynomials[j + 2].value
(quadrature_point (1));
-
- assembling_matrix (i * deg + j, q_point)
- = tmp1 * l_j;
- system_rhs (i * deg + j) += tmp2 * l_j;
- }
}
}
- system_matrix.reinit (assembling_matrix.m (),
- assembling_matrix.m ());
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
- system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix_inv.m ());
system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
quadrature_point, 0);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
+ *= face_quadrature.weight (q_point);
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
const double tmp
- = right_hand_side_value * L_i;
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
-// const double weight
-// = 2.0 * edge_quadrature_x.weight (q_point);
+ const double weight
+ = edge_quadrature_x.weight (q_point);
for (unsigned int i = 0; i < 2; ++i)
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
{
interpolation_matrix (i * source_fe.degree, dof)
- += this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- Point<dim>
- (0.5 * i,
- edge_quadrature_y_points[q_point] (0), 0.0),
- 1);
+ += weight
+ * this->shape_value_component
+ (this->face_to_cell_index (dof, 4),
+ Point<dim>
+ (0.5 * i,
+ edge_quadrature_y_points[q_point] (0), 0.0),
+ 1);
interpolation_matrix ((i + 2) * source_fe.degree,
dof)
- += this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- Point<dim>
- (edge_quadrature_x_points[q_point] (0),
- 0.5 * (i + 1), 0.0), 0);
+ += weight
+ * this->shape_value_component
+ (this->face_to_cell_index (dof, 4),
+ Point<dim>
+ (edge_quadrature_x_points[q_point] (0),
+ 0.5 * (i + 1), 0.0), 0);
}
}
// edges of a face.
std::vector<std::vector<double> >
edge_quadrature_points
- (GeometryInfo<dim>::lines_per_face);
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<double> (n_edge_points));
std::vector<std::vector<double> >
edge_quadrature_weights
- (GeometryInfo<dim>::lines_per_face);
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<double> (n_edge_points));
std::vector<std::vector<Point<dim> > >
edge_quadrature_points_full_dim
- (GeometryInfo<dim>::lines_per_face);
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face; ++line)
- {
- edge_quadrature_points.resize (n_edge_points);
- edge_quadrature_points_full_dim.resize
- (n_edge_points);
- edge_quadrature_weights.resize (n_edge_points);
- }
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<Point<dim> > (n_edge_points));
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
(edge_quadrature_x_points[q_point] (0), 1.0,
0.0);
edge_quadrature_weights[0][q_point]
- = std::sqrt (2.0 * edge_quadrature_y.weight
- (q_point));
+ = std::sqrt (edge_quadrature_y.weight (q_point));
edge_quadrature_weights[1][q_point]
= edge_quadrature_weights[0][q_point];
edge_quadrature_weights[2][q_point]
- = std::sqrt (2.0 * edge_quadrature_x.weight
- (q_point));
+ = std::sqrt (edge_quadrature_x.weight (q_point));
edge_quadrature_weights[3][q_point]
= edge_quadrature_weights[2][q_point];
}
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- for (unsigned int dof = 0; dof < this->dofs_per_face;
- ++dof)
- {
// Set up the system matrix.
// This can be used for all
// edges.
- for (unsigned int q_point = 0;
- q_point < n_edge_points; ++q_point)
- {
- const double weight
- = std::sqrt (2.0 * edge_quadrature_y.weight
- (q_point));
- const double tmp
- = 2.0 * edge_quadrature_y_points[q_point] (0)
- - 1.0;
-
- for (unsigned int i = 0; i < deg; ++i)
- assembling_matrix (i, q_point)
- = weight
- * lobatto_polynomials_grad[i + 1].value
- (tmp);
- }
+ for (unsigned int q_point = 0;
+ q_point < n_edge_points; ++q_point)
+ {
+ const double weight
+ = std::sqrt (edge_quadrature_y.weight
+ (q_point));
+ const double tmp
+ = 2.0 * edge_quadrature_y_points[q_point] (0)
+ - 1.0;
+
+ for (unsigned int i = 0; i < deg; ++i)
+ assembling_matrix (i, q_point)
+ = weight
+ * lobatto_polynomials_grad[i + 1].value
+ (tmp);
+ }
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.invert (system_matrix);
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.invert (system_matrix);
+ Vector<double> system_rhs (system_matrix.m ());
+ Vector<double> solution (system_rhs.size ());
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
for (unsigned int line = 0;
line < GeometryInfo<dim - 1>::lines_per_cell;
++line)
= solution (i);
}
- assembling_matrix.reinit (deg * this->degree,
- n_face_points);
- system_rhs.reinit (assembling_matrix.m ());
- system_rhs = 0;
+ assembling_matrix.reinit (deg * this->degree,
+ n_face_points);
+
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ const double weight
+ = std::sqrt (face_quadrature.weight (q_point));
+ const Point<dim - 1> quadrature_point
+ (2.0 * face_quadrature_points[q_point] (0),
+ 2.0 * face_quadrature_points[q_point] (1) - 1.0);
+
+ for (unsigned int i = 0; i <= deg; ++i)
+ {
+ const double tmp
+ = weight * legendre_polynomials[i].value
+ (quadrature_point (0));
+
+ for (unsigned int j = 0; j < deg; ++j)
+ assembling_matrix (i * deg + j, q_point)
+ = tmp * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
+ }
+ }
+
+ system_matrix.reinit (assembling_matrix.m (),
+ assembling_matrix.m ());
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.reinit (system_matrix.m (),
+ system_matrix.m ());
+ system_matrix_inv.invert (system_matrix);
+ solution.reinit (system_matrix_inv.m ());
+ system_rhs.reinit (assembling_matrix.m ());
+ system_rhs = 0;
// Now we project the remaining
// part on the face shape
// functions. First on the
// horizontal ones, then on
// the vertical ones.
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
+ {
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
quadrature_point, 1);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
+ *= face_quadrature.weight (q_point);
- const double weight
- = std::sqrt (4.0 * reference_face_quadrature.weight
- (q_point));
-
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
- const double tmp1 = weight * L_i;
- const double tmp2
- = right_hand_side_value * L_i;
+ const double tmp
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
- {
- const double l_j
- = lobatto_polynomials[j + 2].value
- (quadrature_point (1));
-
- assembling_matrix (i * deg + j, q_point)
- = tmp1 * l_j;
- system_rhs (i * deg + j) += tmp2 * l_j;
- }
+ system_rhs (i * deg + j)
+ += tmp
+ * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
}
}
- system_matrix.reinit (assembling_matrix.m (),
- assembling_matrix.m ());
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
- system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix_inv.m ());
system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
(i * source_fe.degree + j,
quadrature_point, 0);
- right_hand_side_value *= 4.0 * face_quadrature.weight (q_point);
+ right_hand_side_value *= face_quadrature.weight
+ (q_point);
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
- const double tmp = right_hand_side_value * L_i;
+ const double tmp
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
-// const double weight
-// = 2.0 * edge_quadrature.weight (q_point);
+ const double weight = edge_quadrature.weight (q_point);
for (unsigned int i = 0; i < 2; ++i)
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
{
interpolation_matrix (i * source_fe.degree, dof)
- += this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- Point<dim>
- (0.5 * (i + 1),
- edge_quadrature_points[q_point] (0), 0.0),
- 1);
+ += weight
+ * this->shape_value_component
+ (this->face_to_cell_index (dof, 4),
+ Point<dim>
+ (0.5 * (i + 1),
+ edge_quadrature_points[q_point] (0), 0.0),
+ 1);
interpolation_matrix ((i + 2) * source_fe.degree,
dof)
- += this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- Point<dim>
- (edge_quadrature_points[q_point] (0),
- 0.5 * (i + 1), 0.0), 0);
+ += weight
+ * this->shape_value_component
+ (this->face_to_cell_index (dof, 4),
+ Point<dim>
+ (edge_quadrature_points[q_point] (0),
+ 0.5 * (i + 1), 0.0), 0);
}
}
// face.
std::vector<std::vector<Point<dim> > >
edge_quadrature_points_full_dim
- (GeometryInfo<dim>::lines_per_face);
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face; ++line)
- edge_quadrature_points_full_dim.resize
- (n_edge_points);
+ (GeometryInfo<dim>::lines_per_face,
+ std::vector<Point<dim> > (n_edge_points));
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
1.0, 0.0);
}
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- for (unsigned int dof = 0; dof < this->dofs_per_face;
- ++dof)
- {
// Set up the system matrix.
// This can be used for all
// edges.
- for (unsigned int q_point = 0;
- q_point < n_edge_points; ++q_point)
- {
- const double tmp
- = 2.0 * edge_quadrature_points[q_point] (0)
- - 1.0;
- const double weight
- = std::sqrt (2.0 * edge_quadrature.weight
- (q_point));
-
- for (unsigned int i = 0; i < deg; ++i)
- assembling_matrix (i, q_point)
- = weight
- * lobatto_polynomials_grad[i + 1].value
- (tmp);
- }
+ for (unsigned int q_point = 0;
+ q_point < n_edge_points; ++q_point)
+ {
+ const double tmp
+ = 2.0 * edge_quadrature_points[q_point] (0)
+ - 1.0;
+ const double weight
+ = std::sqrt (edge_quadrature.weight (q_point));
+
+ for (unsigned int i = 0; i < deg; ++i)
+ assembling_matrix (i, q_point)
+ = weight
+ * lobatto_polynomials_grad[i + 1].value
+ (tmp);
+ }
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.invert (system_matrix);
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.invert (system_matrix);
+ Vector<double> system_rhs (system_matrix.m ());
+ Vector<double> solution (system_rhs.size ());
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
for (unsigned int line = 0;
line < GeometryInfo<dim - 1>::lines_per_cell;
++line)
q_point < n_edge_points; ++q_point)
{
const double right_hand_side_value
- = std::sqrt (2.0 * edge_quadrature.weight
- (q_point))
+ = std::sqrt (edge_quadrature.weight
+ (q_point))
* (this->shape_value_component
(this->face_to_cell_index (dof, 4),
edge_quadrature_points_full_dim[line][q_point],
= solution (i);
}
- assembling_matrix.reinit (deg * this->degree,
- n_face_points);
- system_rhs.reinit (assembling_matrix.m ());
- system_rhs = 0;
+ assembling_matrix.reinit (deg * this->degree,
+ n_face_points);
+
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ const double weight
+ = std::sqrt (face_quadrature.weight
+ (q_point));
+ const Point<dim - 1> quadrature_point
+ (2.0 * face_quadrature_points[q_point] (0) - 1.0,
+ 2.0 * face_quadrature_points[q_point] (1) - 1.0);
+
+ for (unsigned int i = 0; i <= deg; ++i)
+ {
+ const double tmp
+ = weight * legendre_polynomials[i].value
+ (quadrature_point (0));
+
+ for (unsigned int j = 0; j < deg; ++j)
+ assembling_matrix (i * deg + j, q_point)
+ = tmp * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
+ }
+ }
+
+ system_matrix.reinit (assembling_matrix.m (),
+ assembling_matrix.m ());
+ assembling_matrix.mTmult (system_matrix,
+ assembling_matrix);
+ system_matrix_inv.reinit (system_matrix.m (),
+ system_matrix.m ());
+ system_matrix_inv.invert (system_matrix);
+ solution.reinit (system_matrix.m ());
+ system_rhs.reinit (assembling_matrix.m ());
+ system_rhs = 0;
+ for (unsigned int dof = 0; dof < this->dofs_per_face;
+ ++dof)
+ {
// Now we project the remaining
// part on the face shape
// functions. First on the
quadrature_point, 1);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
-
- const double weight
- = std::sqrt (4.0 * face_quadrature.weight
- (q_point));
+ *= face_quadrature.weight (q_point);
for (unsigned int i = 0; i <= deg; ++i)
{
- const double L_i
- = legendre_polynomials[i].value
- (quadrature_point (0));
- const double tmp1 = weight * L_i;
- const double tmp2 = right_hand_side_value
- * L_i;
+ const double tmp
+ = right_hand_side_value
+ * legendre_polynomials[i].value
+ (quadrature_point (0));
for (unsigned int j = 0; j < deg; ++j)
- {
- const double l_j
- = lobatto_polynomials[j + 2].value
- (quadrature_point (1));
-
- assembling_matrix (i * deg + j, q_point)
- = tmp1 * l_j;
- system_rhs (i * deg + j) += tmp2 * l_j;
- }
+ system_rhs (i * deg + j)
+ += tmp
+ * lobatto_polynomials[j + 2].value
+ (quadrature_point (1));
}
}
- system_matrix.reinit (assembling_matrix.m (),
- assembling_matrix.m ());
- assembling_matrix.mTmult (system_matrix,
- assembling_matrix);
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
- system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix.m ());
system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
quadrature_point, 0);
right_hand_side_value
- *= 4.0 * face_quadrature.weight (q_point);
+ *= face_quadrature.weight (q_point);
for (unsigned int i = 0; i <= deg; ++i)
{
local_dofs[(i + 8) * this->degree] = 0.0;
}
- if (offset+1 < dim)
+ if (offset + 1 < dim)
{
for (unsigned int i = 0; i < 2; ++i)
for (unsigned int j = 0; j < 2; ++j)
system_rhs = 0;
if ((((line == 0) || (line == 1) || (line == 4) ||
- (line == 5)) && (offset+1 < dim)) ||
+ (line == 5)) && (offset + 1 < dim)) ||
(((line == 2) || (line == 3) || (line == 6) ||
(line == 7)) && (offset == 0)) || (line > 7))
{
{
case 0:
{
- if (offset+1 < dim)
+ if (offset + 1 < dim)
{
// Set up the right hand side
// for the horizontal shape
case 1:
{
- if (offset+1 < dim)
+ if (offset + 1 < dim)
{
// Set up the right hand side
// for the horizontal shape
case 4:
{
- if (offset+1 < dim)
+ if (offset + 1 < dim)
{
// Set up the right hand side
// for the horizontal shape
}
default:
- if (offset+1 < dim)
+ if (offset + 1 < dim)
{
// Set up the right hand side
// for the horizontal shape
system_rhs.reinit (system_matrix_inv.m ());
solution.reinit (system_matrix.m ());
- if (offset+1 < dim)
+ if (offset + 1 < dim)
{
if (offset == 0)
{
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points];
- for (int i = 2;i < (int) GeometryInfo<dim>::lines_per_face; ++i)
+ for (unsigned int i = 2;
+ i < GeometryInfo<dim>::lines_per_face; ++i)
for (unsigned int j = 0; j <= deg; ++j)
tmp -= local_dofs[edge_indices[face][i]
* this->degree + j]
+DEAL::Degree 0:
DEAL::2d, 4 CELLS
DEAL:: CELL 1.0
DEAL::, xq=0.00 0.00, f=[4.00 0]
DEAL::, xq=0.25 0.50 0.50, f=[36.00 32.00 30.00]
DEAL::, xq=0.50 0.50 0.50, f=[36.00 8.00 20.00]
DEAL::
+DEAL::Degree 1:
+DEAL::2d, 4 CELLS
+DEAL:: CELL 1.0
+DEAL::, xq=0.00 0.00, f=[-9.32 -3.46]
+DEAL::, xq=0.50 0.00, f=[25.32 -6.39]
+DEAL::, xq=0.00 0.50, f=[-12.25 3.46]
+DEAL::, xq=0.50 0.50, f=[36.25 14.39]
+DEAL::
+DEAL:: CELL 1.1
+DEAL::, xq=0.50 0.00, f=[-23.96 -6.39]
+DEAL::, xq=1.00 0.00, f=[79.96 -21.03]
+DEAL::, xq=0.50 0.50, f=[-26.89 14.39]
+DEAL::, xq=1.00 0.50, f=[90.89 69.03]
+DEAL::
+DEAL:: CELL 1.2
+DEAL::, xq=0.00 0.50, f=[-12.25 -35.67]
+DEAL::, xq=0.50 0.50, f=[36.25 -38.60]
+DEAL::, xq=0.00 1.00, f=[-41.53 123.67]
+DEAL::, xq=0.50 1.00, f=[145.53 134.60]
+DEAL::
+DEAL:: CELL 1.3
+DEAL::, xq=0.50 0.50, f=[-26.89 -38.60]
+DEAL::, xq=1.00 0.50, f=[90.89 -50.32]
+DEAL::, xq=0.50 1.00, f=[-53.24 134.60]
+DEAL::, xq=1.00 1.00, f=[189.24 178.32]
+DEAL::
+DEAL::2d, 7 CELLS
+DEAL:: CELL 1.1
+DEAL::, xq=0.50 0.00, f=[-9.32 -3.46]
+DEAL::, xq=1.00 0.00, f=[25.32 -6.39]
+DEAL::, xq=0.50 0.50, f=[-12.25 3.46]
+DEAL::, xq=1.00 0.50, f=[36.25 14.39]
+DEAL::
+DEAL:: CELL 1.2
+DEAL::, xq=0.00 0.50, f=[-26.89 -21.03]
+DEAL::, xq=0.50 0.50, f=[90.89 -23.96]
+DEAL::, xq=0.00 1.00, f=[-29.82 69.03]
+DEAL::, xq=0.50 1.00, f=[101.82 79.96]
+DEAL::
+DEAL:: CELL 1.3
+DEAL::, xq=0.50 0.50, f=[-12.25 -23.96]
+DEAL::, xq=1.00 0.50, f=[36.25 -38.60]
+DEAL::, xq=0.50 1.00, f=[-41.53 79.96]
+DEAL::, xq=1.00 1.00, f=[145.53 134.60]
+DEAL::
+DEAL:: CELL 2.0
+DEAL::, xq=0.00 0.00, f=[-112.34 -100.63]
+DEAL::, xq=0.25 0.00, f=[400.34 -106.49]
+DEAL::, xq=0.00 0.25, f=[-118.20 356.63]
+DEAL::, xq=0.25 0.25, f=[422.20 378.49]
+DEAL::
+DEAL:: CELL 2.1
+DEAL::, xq=0.25 0.00, f=[-141.63 -106.49]
+DEAL::, xq=0.50 0.00, f=[509.63 -5.20]
+DEAL::, xq=0.25 0.25, f=[-147.48 378.49]
+DEAL::, xq=0.50 0.25, f=[531.48 1.73]
+DEAL::
+DEAL:: CELL 2.2
+DEAL::, xq=0.00 0.25, f=[-118.20 -165.05]
+DEAL::, xq=0.25 0.25, f=[422.20 -170.91]
+DEAL::, xq=0.00 0.50, f=[-56.33 597.05]
+DEAL::, xq=0.25 0.50, f=[61.44 618.91]
+DEAL::
+DEAL:: CELL 2.3
+DEAL::, xq=0.25 0.25, f=[-147.48 -170.91]
+DEAL::, xq=0.50 0.25, f=[531.48 -1.73]
+DEAL::, xq=0.25 0.50, f=[2.56 618.91]
+DEAL::, xq=0.50 0.50, f=[120.33 5.20]
+DEAL::
+DEAL::3d, 8 CELLS
+DEAL:: CELL 1.0
+DEAL::, xq=0.00 0.00 0.00, f=[-9.32 -3.46 -26.89]
+DEAL::, xq=0.50 0.00 0.00, f=[25.32 -6.39 -29.82]
+DEAL::, xq=0.00 0.50 0.00, f=[-12.25 3.46 -32.75]
+DEAL::, xq=0.50 0.50 0.00, f=[36.25 14.39 -35.67]
+DEAL::, xq=0.00 0.00 0.50, f=[-21.03 -15.18 90.89]
+DEAL::, xq=0.50 0.00 0.50, f=[69.03 -18.11 101.82]
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 47.18 112.75]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 58.11 123.67]
+DEAL::
+DEAL:: CELL 1.1
+DEAL::, xq=0.50 0.00 0.00, f=[-85.45 -6.39 -29.82]
+DEAL::, xq=1.00 0.00 0.00, f=[309.45 -82.53 -100.09]
+DEAL::, xq=0.50 0.50 0.00, f=[-88.38 14.39 -35.67]
+DEAL::, xq=1.00 0.50 0.00, f=[320.38 298.53 -103.02]
+DEAL::, xq=0.50 0.00 0.50, f=[-94.24 -18.11 101.82]
+DEAL::, xq=1.00 0.00 0.50, f=[342.24 -91.31 364.09]
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 58.11 123.67]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 331.31 375.02]
+DEAL::
+DEAL:: CELL 1.2
+DEAL::, xq=0.00 0.50 0.00, f=[-12.25 -144.02 -32.75]
+DEAL::, xq=0.50 0.50 0.00, f=[36.25 -146.95 -35.67]
+DEAL::, xq=0.00 1.00 0.00, f=[-149.87 528.02 -161.59]
+DEAL::, xq=0.50 1.00 0.00, f=[549.87 538.95 -164.52]
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 -152.80 112.75]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 -155.73 123.67]
+DEAL::, xq=0.00 1.00 0.50, f=[-158.66 560.80 593.59]
+DEAL::, xq=0.50 1.00 0.50, f=[582.66 571.73 604.52]
+DEAL::
+DEAL:: CELL 1.3
+DEAL::, xq=0.50 0.50 0.00, f=[-88.38 -146.95 -35.67]
+DEAL::, xq=1.00 0.50 0.00, f=[320.38 -205.51 -103.02]
+DEAL::, xq=0.50 1.00 0.00, f=[-208.44 538.95 -164.52]
+DEAL::, xq=1.00 1.00 0.00, f=[768.44 757.51 -217.22]
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 -155.73 123.67]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 -211.37 375.02]
+DEAL::, xq=0.50 1.00 0.50, f=[-214.29 571.73 604.52]
+DEAL::, xq=1.00 1.00 0.50, f=[790.29 779.37 801.22]
+DEAL::
+DEAL:: CELL 1.4
+DEAL::, xq=0.00 0.00 0.50, f=[-21.03 -15.18 -264.07]
+DEAL::, xq=0.50 0.00 0.50, f=[69.03 -18.11 -267.00]
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 47.18 -269.93]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 58.11 -272.86]
+DEAL::, xq=0.00 0.00 1.00, f=[-258.22 -252.36 976.07]
+DEAL::, xq=0.50 0.00 1.00, f=[954.22 -255.29 987.00]
+DEAL::, xq=0.00 0.50 1.00, f=[-261.15 932.36 997.93]
+DEAL::, xq=0.50 0.50 1.00, f=[965.15 943.29 1008.86]
+DEAL::
+DEAL:: CELL 1.5
+DEAL::, xq=0.50 0.00 0.50, f=[-94.24 -18.11 -267.00]
+DEAL::, xq=1.00 0.00 0.50, f=[342.24 -91.31 -322.64]
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 58.11 -272.86]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 331.31 -325.57]
+DEAL::, xq=0.50 0.00 1.00, f=[-316.78 -255.29 987.00]
+DEAL::, xq=1.00 0.00 1.00, f=[1172.78 -313.85 1194.64]
+DEAL::, xq=0.50 0.50 1.00, f=[-319.71 943.29 1008.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1183.71 1161.85 1205.57]
+DEAL::
+DEAL:: CELL 1.6
+DEAL::, xq=0.00 0.50 0.50, f=[-23.96 -152.80 -269.93]
+DEAL::, xq=0.50 0.50 0.50, f=[79.96 -155.73 -272.86]
+DEAL::, xq=0.00 1.00 0.50, f=[-158.66 560.80 -369.49]
+DEAL::, xq=0.50 1.00 0.50, f=[582.66 571.73 -372.42]
+DEAL::, xq=0.00 0.50 1.00, f=[-261.15 -360.70 997.93]
+DEAL::, xq=0.50 0.50 1.00, f=[965.15 -363.63 1008.86]
+DEAL::, xq=0.00 1.00 1.00, f=[-366.56 1336.70 1369.49]
+DEAL::, xq=0.50 1.00 1.00, f=[1358.56 1347.63 1380.42]
+DEAL::
+DEAL:: CELL 1.7
+DEAL::, xq=0.50 0.50 0.50, f=[-97.17 -155.73 -272.86]
+DEAL::, xq=1.00 0.50 0.50, f=[353.17 -211.37 -325.57]
+DEAL::, xq=0.50 1.00 0.50, f=[-214.29 571.73 -372.42]
+DEAL::, xq=1.00 1.00 0.50, f=[790.29 779.37 -413.41]
+DEAL::, xq=0.50 0.50 1.00, f=[-319.71 -363.63 1008.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1183.71 -407.56 1205.57]
+DEAL::, xq=0.50 1.00 1.00, f=[-410.48 1347.63 1380.42]
+DEAL::, xq=1.00 1.00 1.00, f=[1522.48 1511.56 1533.41]
+DEAL::
+DEAL::3d, 15 CELLS
+DEAL:: CELL 1.1
+DEAL::, xq=0.50 0.00 0.00, f=[-9.32 -3.46 -26.89]
+DEAL::, xq=1.00 0.00 0.00, f=[25.32 -6.39 -29.82]
+DEAL::, xq=0.50 0.50 0.00, f=[-12.25 3.46 -32.75]
+DEAL::, xq=1.00 0.50 0.00, f=[36.25 14.39 -35.67]
+DEAL::, xq=0.50 0.00 0.50, f=[-21.03 -15.18 90.89]
+DEAL::, xq=1.00 0.00 0.50, f=[69.03 -18.11 101.82]
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 47.18 112.75]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 58.11 123.67]
+DEAL::
+DEAL:: CELL 1.2
+DEAL::, xq=0.00 0.50 0.00, f=[-88.38 -82.53 -105.95]
+DEAL::, xq=0.50 0.50 0.00, f=[320.38 -85.45 -32.75]
+DEAL::, xq=0.00 1.00 0.00, f=[-91.31 298.53 -108.88]
+DEAL::, xq=0.50 1.00 0.00, f=[331.31 309.45 -111.81]
+DEAL::, xq=0.00 0.50 0.50, f=[-100.09 -94.24 385.95]
+DEAL::, xq=0.50 0.50 0.50, f=[364.09 -97.17 112.75]
+DEAL::, xq=0.00 1.00 0.50, f=[-103.02 342.24 396.88]
+DEAL::, xq=0.50 1.00 0.50, f=[375.02 353.17 407.81]
+DEAL::
+DEAL:: CELL 1.3
+DEAL::, xq=0.50 0.50 0.00, f=[-12.25 -85.45 -32.75]
+DEAL::, xq=1.00 0.50 0.00, f=[36.25 -158.66 -35.67]
+DEAL::, xq=0.50 1.00 0.00, f=[-161.59 309.45 -111.81]
+DEAL::, xq=1.00 1.00 0.00, f=[593.59 582.66 -170.37]
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 -97.17 112.75]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 -164.52 123.67]
+DEAL::, xq=0.50 1.00 0.50, f=[-167.44 353.17 407.81]
+DEAL::, xq=1.00 1.00 0.50, f=[615.44 604.52 626.37]
+DEAL::
+DEAL:: CELL 1.4
+DEAL::, xq=0.00 0.00 0.50, f=[-208.44 -205.51 -223.08]
+DEAL::, xq=0.50 0.00 0.50, f=[768.44 -15.18 -226.01]
+DEAL::, xq=0.00 0.50 0.50, f=[-100.09 757.51 -228.94]
+DEAL::, xq=0.50 0.50 0.50, f=[364.09 47.18 -231.86]
+DEAL::, xq=0.00 0.00 1.00, f=[-217.22 -211.37 823.08]
+DEAL::, xq=0.50 0.00 1.00, f=[801.22 -214.29 834.01]
+DEAL::, xq=0.00 0.50 1.00, f=[-220.15 779.37 844.94]
+DEAL::, xq=0.50 0.50 1.00, f=[812.15 790.29 855.86]
+DEAL::
+DEAL:: CELL 1.5
+DEAL::, xq=0.50 0.00 0.50, f=[-21.03 -15.18 -226.01]
+DEAL::, xq=1.00 0.00 0.50, f=[69.03 -18.11 -287.50]
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 47.18 -231.86]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 58.11 -290.43]
+DEAL::, xq=0.50 0.00 1.00, f=[-281.64 -214.29 834.01]
+DEAL::, xq=1.00 0.00 1.00, f=[1041.64 -278.72 1063.50]
+DEAL::, xq=0.50 0.50 1.00, f=[-284.57 790.29 855.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1052.57 1030.72 1074.43]
+DEAL::
+DEAL:: CELL 1.6
+DEAL::, xq=0.00 0.50 0.50, f=[-100.09 -94.24 -228.94]
+DEAL::, xq=0.50 0.50 0.50, f=[364.09 -97.17 -231.86]
+DEAL::, xq=0.00 1.00 0.50, f=[-103.02 342.24 -334.35]
+DEAL::, xq=0.50 1.00 0.50, f=[375.02 353.17 -337.28]
+DEAL::, xq=0.00 0.50 1.00, f=[-220.15 -325.57 844.94]
+DEAL::, xq=0.50 0.50 1.00, f=[812.15 -328.49 855.86]
+DEAL::, xq=0.00 1.00 1.00, f=[-331.42 1205.57 1238.35]
+DEAL::, xq=0.50 1.00 1.00, f=[1227.42 1216.49 1249.28]
+DEAL::
+DEAL:: CELL 1.7
+DEAL::, xq=0.50 0.50 0.50, f=[-23.96 -97.17 -231.86]
+DEAL::, xq=1.00 0.50 0.50, f=[79.96 -164.52 -290.43]
+DEAL::, xq=0.50 1.00 0.50, f=[-167.44 353.17 -337.28]
+DEAL::, xq=1.00 1.00 0.50, f=[615.44 604.52 -378.27]
+DEAL::, xq=0.50 0.50 1.00, f=[-284.57 -328.49 855.86]
+DEAL::, xq=1.00 0.50 1.00, f=[1052.57 -372.42 1074.43]
+DEAL::, xq=0.50 1.00 1.00, f=[-375.35 1216.49 1249.28]
+DEAL::, xq=1.00 1.00 1.00, f=[1391.35 1380.42 1402.27]
+DEAL::
+DEAL:: CELL 2.0
+DEAL::, xq=0.00 0.00 0.00, f=[-826.83 -815.11 -861.96]
+DEAL::, xq=0.25 0.00 0.00, f=[3066.83 -820.97 -867.82]
+DEAL::, xq=0.00 0.25 0.00, f=[-832.68 3023.11 -873.68]
+DEAL::, xq=0.25 0.25 0.00, f=[3088.68 3044.97 -879.53]
+DEAL::, xq=0.00 0.00 0.25, f=[-850.25 -838.54 3197.96]
+DEAL::, xq=0.25 0.00 0.25, f=[3154.25 -844.39 3219.82]
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 3110.54 3241.68]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 3132.39 3263.53]
+DEAL::
+DEAL:: CELL 2.1
+DEAL::, xq=0.25 0.00 0.00, f=[-979.09 -820.97 -867.82]
+DEAL::, xq=0.50 0.00 0.00, f=[3635.09 -31.84 -115.22]
+DEAL::, xq=0.25 0.25 0.00, f=[-984.95 3044.97 -879.53]
+DEAL::, xq=0.50 0.25 0.00, f=[3656.95 28.38 -49.81]
+DEAL::, xq=0.25 0.00 0.25, f=[-996.66 -844.39 3219.82]
+DEAL::, xq=0.50 0.00 0.25, f=[3700.66 34.34 120.33]
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 3132.39 3263.53]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -39.73 51.46]
+DEAL::
+DEAL:: CELL 2.2
+DEAL::, xq=0.00 0.25 0.00, f=[-832.68 -1096.22 -873.68]
+DEAL::, xq=0.25 0.25 0.00, f=[3088.68 -1102.08 -879.53]
+DEAL::, xq=0.00 0.50 0.00, f=[-394.95 4072.22 -83.31]
+DEAL::, xq=0.25 0.50 0.00, f=[422.57 4094.08 134.29]
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 -1113.79 3241.68]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 -1119.65 3263.53]
+DEAL::, xq=0.00 0.50 0.25, f=[-169.82 4137.79 117.36]
+DEAL::, xq=0.25 0.50 0.25, f=[181.13 4159.65 -131.63]
+DEAL::
+DEAL:: CELL 2.3
+DEAL::, xq=0.25 0.25 0.00, f=[-984.95 -1102.08 -879.53]
+DEAL::, xq=0.50 0.25 0.00, f=[3656.95 43.56 -49.81]
+DEAL::, xq=0.25 0.50 0.00, f=[-190.57 4094.08 134.29]
+DEAL::, xq=0.50 0.50 0.00, f=[626.95 -40.09 -141.87]
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 -1119.65 3263.53]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -75.91 51.46]
+DEAL::, xq=0.25 0.50 0.25, f=[-82.09 4159.65 -131.63]
+DEAL::, xq=0.50 0.50 0.25, f=[268.87 30.16 149.12]
+DEAL::
+DEAL:: CELL 2.4
+DEAL::, xq=0.00 0.00 0.25, f=[-850.25 -838.54 -1336.33]
+DEAL::, xq=0.25 0.00 0.25, f=[3154.25 -844.39 -1342.19]
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 3110.54 -1348.05]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 3132.39 -1353.90]
+DEAL::, xq=0.00 0.00 0.50, f=[-380.71 -927.78 4968.33]
+DEAL::, xq=0.25 0.00 0.50, f=[452.27 16.47 4990.19]
+DEAL::, xq=0.00 0.25 0.50, f=[571.27 998.27 5012.05]
+DEAL::, xq=0.25 0.25 0.50, f=[-558.41 -20.16 5033.90]
+DEAL::
+DEAL:: CELL 2.5
+DEAL::, xq=0.25 0.00 0.25, f=[-996.66 -844.39 -1342.19]
+DEAL::, xq=0.50 0.00 0.25, f=[3700.66 34.34 -56.33]
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 3132.39 -1353.90]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -39.73 -24.49]
+DEAL::, xq=0.25 0.00 0.50, f=[560.05 16.47 4990.19]
+DEAL::, xq=0.50 0.00 0.50, f=[488.39 104.59 179.22]
+DEAL::, xq=0.25 0.25 0.50, f=[59.35 -20.16 5033.90]
+DEAL::, xq=0.50 0.25 0.50, f=[149.04 -103.77 76.78]
+DEAL::
+DEAL:: CELL 2.6
+DEAL::, xq=0.00 0.25 0.25, f=[-856.11 -1113.79 -1348.05]
+DEAL::, xq=0.25 0.25 0.25, f=[3176.11 -1119.65 -1353.90]
+DEAL::, xq=0.00 0.50 0.25, f=[-169.82 4137.79 386.73]
+DEAL::, xq=0.25 0.50 0.25, f=[181.13 4159.65 -131.78]
+DEAL::, xq=0.00 0.25 0.50, f=[571.27 -446.27 5012.05]
+DEAL::, xq=0.25 0.25 0.50, f=[-558.41 7.31 5033.90]
+DEAL::, xq=0.00 0.50 0.50, f=[-448.24 1479.78 139.22]
+DEAL::, xq=0.25 0.50 0.50, f=[480.14 -29.31 198.13]
+DEAL::
+DEAL:: CELL 2.7
+DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 -1119.65 -1353.90]
+DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -75.91 -24.49]
+DEAL::, xq=0.25 0.50 0.25, f=[-82.09 4159.65 -131.78]
+DEAL::, xq=0.50 0.50 0.25, f=[268.87 30.16 76.37]
+DEAL::, xq=0.25 0.25 0.50, f=[59.35 7.31 5033.90]
+DEAL::, xq=0.50 0.25 0.50, f=[149.04 31.59 76.78]
+DEAL::, xq=0.25 0.50 0.50, f=[248.05 -29.31 198.13]
+DEAL::, xq=0.50 0.50 0.50, f=[248.05 31.59 76.37]
+DEAL::