#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
void assemble_system (const bool reconstruct_trace = false);
void solve ();
void postprocess ();
- void refine_mesh ();
+ void refine_grid (const unsigned int cylce);
void output_results (const unsigned int cycle);
Triangulation<dim> triangulation;
const unsigned int n_q_points = quadrature.size();
std::vector<double> u_values(n_q_points);
std::vector<Tensor<1,dim> > u_gradients(n_q_points);
- FEValuesExtractors::Vector gradients(0);
- FEValuesExtractors::Scalar values(dim);
+ FEValuesExtractors::Vector fluxes(0);
+ FEValuesExtractors::Scalar scalar(dim);
FEValues<dim> fe_values_local(mapping, fe_local, quadrature, update_values);
FullMatrix<double> cell_matrix(fe_u_post.dofs_per_cell,
fe_u_post.dofs_per_cell);
fe_values.reinit(cell);
fe_values_local.reinit(cell_loc);
- fe_values_local[values].get_function_values(solution_local, u_values);
- fe_values_local[gradients].get_function_values(solution_local, u_gradients);
+ fe_values_local[scalar].get_function_values(solution_local, u_values);
+ fe_values_local[fluxes].get_function_values(solution_local, u_gradients);
for (unsigned int i=1; i<fe_u_post.dofs_per_cell; ++i)
{
for (unsigned int j=0; j<fe_u_post.dofs_per_cell; ++j)
-
template <int dim>
-void Step51<dim>::run ()
+void Step51<dim>::refine_grid (const unsigned int cycle)
{
const bool do_cube = true;
- if (!do_cube)
+ if (cycle == 0)
{
- GridGenerator::hyper_ball (triangulation);
- static const HyperBallBoundary<dim> boundary;
- triangulation.set_boundary(0, boundary);
- triangulation.refine_global(6-2*dim);
+ if (!do_cube)
+ {
+ GridGenerator::hyper_ball (triangulation);
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary(0, boundary);
+ triangulation.refine_global(6-2*dim);
+ }
+ else
+ GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1);
}
+ else
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ {
+ if (do_cube)
+ {
+ triangulation.clear();
+ GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+ triangulation.refine_global(3-dim+cycle/2);
+ }
+ else
+ triangulation.refine_global (1);
+ break;
+ }
+
+ case adaptive_refinement:
+ {
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ FEValuesExtractors::Scalar scalar(dim);
+ typename FunctionMap<dim>::type neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (dof_handler_local,
+ QGauss<dim-1>(3),
+ neumann_boundary,
+ solution_local,
+ estimated_error_per_cell,
+ fe_local.component_mask(scalar));
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.);
+
+ triangulation.execute_coarsening_and_refinement ();
+
+ break;
+ }
+
+ default:
+ {
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ }
+
+
+
+
+
+template <int dim>
+void Step51<dim>::run ()
+{
for (unsigned int cycle=0; cycle<10; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
-
- if (do_cube)
- {
- triangulation.clear();
- GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
- triangulation.refine_global(3-dim+cycle/2);
- }
- else triangulation.refine_global(1);
-
+
+ refine_grid (cycle);
setup_system ();
assemble_system (false);
solve ();
/**
- * A finite element, which is a tensor product polynomial on each face
- * and undefined in the interior of the cells. The basis functions on
- * the faces are from Polynomials::LagrangeEquidistant
+ * A finite element, which is a tensor product polynomial on each face and
+ * undefined in the interior of the cells. The basis functions on the faces
+ * are Lagrange polynomials based on the support points of the
+ * (dim-1)-dimensional Gauss--Lobatto quadrature rule. For element degree one
+ * and two, the polynomials hence correspond to the usual Lagrange polynomials
+ * on equidistant points.
*
* This finite element is the trace space of FE_RaviartThomas on the
* faces and serves in hybridized methods.
* element space, but all shape function values extracted will equal
* to zero.
*
- * @todo Polynomials::LagrangeEquidistant should be and will be
- * replaced by Polynomials::LagrangeGaussLobatto as soon as such a
- * polynomial set exists.
- *
* @ingroup fe
- * @author Guido Kanschat
- * @date 2009, 2011
+ * @author Guido Kanschat, Martin Kronbichler
+ * @date 2009, 2011, 2013
*/
template <int dim, int spacedim=dim>
class FE_FaceQ : public FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim>
*/
virtual std::string get_name () const;
+ /**
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ FullMatrix<double> &matrix) const;
+
+ /**
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ const unsigned int subface,
+ FullMatrix<double> &matrix) const;
+
/**
* Check for non-zero values on a face.
*
virtual bool has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const;
+ /**
+ * Return whether this element implements its hanging node constraints in
+ * the new way, which has to be used to make elements "hp compatible".
+ *
+ * For the FE_Q class the result is always true (independent of the degree
+ * of the element), as it implements the complete set of functions necessary
+ * for hp capability.
+ */
+ virtual bool hp_constraints_are_implemented () const;
+
+ /**
+ * Return whether this element dominates the one given as argument when they
+ * meet at a common face, whether it is the other way around, whether
+ * neither dominates, or if either could dominate.
+ *
+ * For a definition of domination, see FiniteElementBase::Domination and in
+ * particular the @ref hp_paper "hp paper".
+ */
+ virtual
+ FiniteElementDomination::Domination
+ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+
private:
/**
* Return vector with dofs per vertex, line, quad, hex.
/**
* Return the matrix interpolating from a face of of one element to the face
* of the neighboring element. The size of the matrix is then
- * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
- *
- * Derived elements will have to implement this function. They may only
- * provide interpolation matrices for certain source finite elements, for
- * example those from the same family. If they don't implement interpolation
- * from a given element, then they must throw an exception of type
- * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented.
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. The
+ * FE_Q element family only provides interpolation matrices for elements of
+ * the same type and FE_Nothing. For all other elements, an exception of
+ * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is
+ * thrown.
*/
virtual void
get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
/**
* Return the matrix interpolating from a face of of one element to the face
* of the neighboring element. The size of the matrix is then
- * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
- *
- * Derived elements will have to implement this function. They may only
- * provide interpolation matrices for certain source finite elements, for
- * example those from the same family. If they don't implement interpolation
- * from a given element, then they must throw an exception of type
- * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented.
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. The
+ * FE_Q element family only provides interpolation matrices for elements of
+ * the same type and FE_Nothing. For all other elements, an exception of
+ * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is
+ * thrown.
*/
virtual void
get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
#include <deal.II/fe/fe_face.h>
#include <deal.II/fe/fe_poly_face.templates.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/base/quadrature_lib.h>
#include <sstream>
DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+ std::vector<Point<1> >
+ get_QGaussLobatto_points (const unsigned int degree)
+ {
+ if (degree > 0)
+ {
+ QGaussLobatto<1> quad(degree+1);
+ return quad.get_points();
+ }
+ else
+ return std::vector<Point<1> >(1, Point<1>(0.5));
+ }
+}
+
template <int dim, int spacedim>
FE_FaceQ<dim,spacedim>::FE_FaceQ (const unsigned int degree)
:
FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim> (
- TensorProductPolynomials<dim-1>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
+ TensorProductPolynomials<dim-1>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
std::vector<bool>(1,true))
{
{
const double step = 1./this->degree;
Point<codim> p;
+ std::vector<Point<1> > points = get_QGaussLobatto_points(degree);
unsigned int k=0;
for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
for (unsigned int ix=0; ix<=this->degree; ++ix)
{
- p(0) = ix * step;
+ p(0) = points[ix][0];
if (codim>1)
- p(1) = iy * step;
+ p(1) = points[iy][0];
if (codim>2)
- p(2) = iz * step;
+ p(2) = points[iz][0];
this->unit_face_support_points[k++] = p;
}
AssertDimension (k, this->unit_face_support_points.size());
}
+
+ // initialize unit support points
+ this->unit_support_points.resize(GeometryInfo<dim>::faces_per_cell*
+ this->unit_face_support_points.size());
+ const unsigned int n_face_dofs = this->unit_face_support_points.size();
+ for (unsigned int i=0; i<n_face_dofs; ++i)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ for (unsigned int e=0, c=0; e<dim; ++e)
+ if (d!=e)
+ {
+ this->unit_support_points[n_face_dofs*2*d+i][e] =
+ this->unit_face_support_points[i][c];
+ this->unit_support_points[n_face_dofs*(2*d+1)+i][e] =
+ this->unit_face_support_points[i][c];
+ this->unit_support_points[n_face_dofs*(2*d+1)+i][d] = 1;
+ ++c;
+ }
+ }
}
+template <int dim, int spacedim>
+void
+FE_FaceQ<dim,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // this function is similar to the respective method in FE_Q
+
+ // this is only implemented, if the source FE is also a FE_FaceQ element
+ AssertThrow ((dynamic_cast<const FE_FaceQ<dim,spacedim> *>(&x_source_fe) != 0),
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented()));
+
+ Assert (interpolation_matrix.n() == this->dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.n(),
+ this->dofs_per_face));
+ Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_face));
+
+ // ok, source is a FaceQ element, so we will be able to do the work
+ const FE_FaceQ<dim,spacedim> &source_fe
+ = dynamic_cast<const FE_FaceQ<dim,spacedim>&>(x_source_fe);
+
+ // Make sure that the element for which the DoFs should be constrained is
+ // the one with the higher polynomial degree. Actually the procedure will
+ // work also if this assertion is not satisfied. But the matrices produced
+ // in that case might lead to problems in the hp procedures, which use this
+ // method.
+ Assert (this->dofs_per_face <= source_fe.dofs_per_face,
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented ()));
+
+ // generate a quadrature with the unit face support points.
+ Quadrature<dim-1> face_quadrature (source_fe.get_unit_face_support_points ());
+
+ // Rule of thumb for FP accuracy, that can be expected for a given
+ // polynomial degree. This value is used to cut off values close to zero.
+ const double eps = 2e-13*this->degree*(dim-1);
+
+ // compute the interpolation matrix by simply taking the value at the
+ // support points.
+ for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
+ {
+ const Point<dim-1> &p = face_quadrature.point (i);
+
+ for (unsigned int j=0; j<this->dofs_per_face; ++j)
+ {
+ double matrix_entry = this->poly_space.compute_value (j, p);
+
+ // Correct the interpolated value. I.e. if it is close to 1 or 0,
+ // make it exactly 1 or 0. Unfortunately, this is required to avoid
+ // problems with higher order elements.
+ if (std::fabs (matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs (matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i,j) = matrix_entry;
+ }
+ }
+
+ // make sure that the row sum of each of the matrices is 1 at this
+ // point. this must be so since the shape functions sum up to 1
+ for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ sum += interpolation_matrix(j,i);
+
+ Assert (std::fabs(sum-1) < 2e-13*this->degree*(dim-1),
+ ExcInternalError());
+ }
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_FaceQ<dim,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // this function is similar to the respective method in FE_Q
+
+ Assert (interpolation_matrix.n() == this->dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.n(),
+ this->dofs_per_face));
+ Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_face));
+
+ // see if source is a FaceQ element
+ if (const FE_FaceQ<dim,spacedim> *source_fe
+ = dynamic_cast<const FE_FaceQ<dim,spacedim> *>(&x_source_fe))
+ {
+
+ // Make sure that the element for which the DoFs should be constrained
+ // is the one with the higher polynomial degree. Actually the procedure
+ // will work also if this assertion is not satisfied. But the matrices
+ // produced in that case might lead to problems in the hp procedures,
+ // which use this method.
+ Assert (this->dofs_per_face <= source_fe->dofs_per_face,
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented ()));
+
+ // generate a quadrature with the unit face support points.
+ Quadrature<dim-1> face_quadrature (source_fe->get_unit_face_support_points ());
+
+ // Rule of thumb for FP accuracy, that can be expected for a given
+ // polynomial degree. This value is used to cut off values close to
+ // zero.
+ const double eps = 2e-13*this->degree*(dim-1);
+
+ // compute the interpolation matrix by simply taking the value at the
+ // support points.
+ for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
+ {
+ const Point<dim-1> p =
+ GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(i),
+ subface);
+
+ for (unsigned int j=0; j<this->dofs_per_face; ++j)
+ {
+ double matrix_entry = this->poly_space.compute_value (j, p);
+
+ // Correct the interpolated value. I.e. if it is close to 1 or 0,
+ // make it exactly 1 or 0. Unfortunately, this is required to avoid
+ // problems with higher order elements.
+ if (std::fabs (matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs (matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i,j) = matrix_entry;
+ }
+ }
+
+ // make sure that the row sum of each of the matrices is 1 at this
+ // point. this must be so since the shape functions sum up to 1
+ for (unsigned int j=0; j<source_fe->dofs_per_face; ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ sum += interpolation_matrix(j,i);
+
+ Assert (std::fabs(sum-1) < 2e-13*this->degree*(dim-1),
+ ExcInternalError());
+ }
+ }
+ else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != 0)
+ {
+ // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+ }
+ else
+ AssertThrow (false,(typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented()));
+}
+
+
+
template <int dim, int spacedim>
bool
FE_FaceQ<dim,spacedim>::has_support_on_face (
}
+
template <int dim, int spacedim>
std::vector<unsigned int>
FE_FaceQ<dim,spacedim>::get_dpo_vector (const unsigned int deg)
+template <int dim, int spacedim>
+bool
+FE_FaceQ<dim,spacedim>::hp_constraints_are_implemented () const
+{
+ return true;
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElementDomination::Domination
+FE_FaceQ<dim,spacedim>::
+compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ if (const FE_FaceQ<dim,spacedim> *fe_q_other
+ = dynamic_cast<const FE_FaceQ<dim,spacedim>*>(&fe_other))
+ {
+ if (this->degree < fe_q_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_q_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no degrees of freedom and it is typically used in
+ // a context where we don't require any continuity along the interface
+ return FiniteElementDomination::no_requirements;
+ }
+
+ Assert (false, ExcNotImplemented());
+ return FiniteElementDomination::neither_element_dominates;
+}
+
// explicit instantiations
#include "fe_face.inst"
{
// dofs are located along lines, so two dofs are identical if they are
// located at identical positions. if we had only equidistant points, we
- // could simple check for similarity like (i+1)*q == (j+1)*p, but we
+ // could simply check for similarity like (i+1)*q == (j+1)*p, but we
// might have other support points (e.g. Gauss-Lobatto
// points). Therefore, read the points in unit_support_points for the
// first coordinate direction. We take the lexicographic ordering of the
DEAL:Face2d-FaceQ2::
DEAL:Face2d-FaceQ2::
DEAL:Face2d-FaceQ3::0.00 0.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.00 0.12 1.12 2.05 0.79 1.04 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.00 0.25 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.00 0.38 1.04 0.79 2.05 1.12 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.00 0.12 1.05 1.99 0.94 1.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.00 0.25 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.00 0.38 1.02 0.94 1.99 1.05 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::0.00 0.50 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::0.50 0.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.06 1.00 1.00 1.00 1.00 1.44 1.80 0.69 1.06 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.12 1.00 1.00 1.00 1.00 1.12 2.05 0.79 1.04 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.19 1.00 1.00 1.00 1.00 0.97 1.92 1.13 0.98 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.06 1.00 1.00 1.00 1.00 1.40 1.73 0.81 1.06 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.12 1.00 1.00 1.00 1.00 1.05 1.99 0.94 1.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.19 1.00 1.00 1.00 1.00 0.89 1.91 1.26 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::
-DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.31 1.00 1.00 1.00 1.00 0.98 1.13 1.92 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.38 1.00 1.00 1.00 1.00 1.04 0.79 2.05 1.12 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.50 0.44 1.00 1.00 1.00 1.00 1.06 0.69 1.80 1.44 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.25 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.31 1.00 1.00 1.00 1.00 0.94 1.26 1.91 0.89 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.38 1.00 1.00 1.00 1.00 1.02 0.94 1.99 1.05 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.50 0.44 1.00 1.00 1.00 1.00 1.06 0.81 1.73 1.40 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::0.50 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::0.00 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.12 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.12 2.05 0.79 1.04 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.25 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94 1.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.38 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.04 0.79 2.05 1.12 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.12 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.05 1.99 0.94 1.02 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.25 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88 1.00 1.00 1.00 1.00
+DEAL:Face2d-FaceQ3::0.38 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.02 0.94 1.99 1.05 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::0.50 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00 1.00
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::0.00 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00 1.00 1.00 1.00
-DEAL:Face2d-FaceQ3::0.12 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.12 2.05 0.79 1.04
-DEAL:Face2d-FaceQ3::0.25 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.94 1.56 1.56 0.94
-DEAL:Face2d-FaceQ3::0.38 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.04 0.79 2.05 1.12
+DEAL:Face2d-FaceQ3::0.12 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.05 1.99 0.94 1.02
+DEAL:Face2d-FaceQ3::0.25 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.88 1.62 1.62 0.88
+DEAL:Face2d-FaceQ3::0.38 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.02 0.94 1.99 1.05
DEAL:Face2d-FaceQ3::0.50 0.50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 2.00
DEAL:Face2d-FaceQ3::
DEAL:Face2d-FaceQ3::