<ol>
+<li>New: VectorTools::project and a whole host of similar functions
+are now also available for objects of type hp::DoFHandler.
+<br>
+(Wolfgang Bangerth, 2013/07/21)
+</li>
+
<li>Fixed: hp::DoFHandler::n_boundary_dofs() had a bug that always led
to a failed assertion. This is now fixed.
<br>
/**
* Calls the project()
- * function, see above, with
+ * function above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
template <int dim, class VECTOR, int spacedim>
Quadrature<dim-1>(0)),
const bool project_to_boundary_first = false);
+ /**
+ * Same as above, but for arguments of type hp::DoFHandler,
+ * hp::QuadratureCollection, hp::MappingCollection
+ */
+ template <int dim, class VECTOR, int spacedim>
+ void project (const hp::MappingCollection<dim, spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const hp::QCollection<dim> &quadrature,
+ const Function<spacedim> &function,
+ VECTOR &vec,
+ const bool enforce_zero_boundary = false,
+ const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
+ QGauss<dim-1>(2) :
+ Quadrature<dim-1>(0)),
+ const bool project_to_boundary_first = false);
+
+ /**
+ * Calls the project()
+ * function above, with a collection of
+ * MappingQ1@<dim@>() objects.
+ */
+ template <int dim, class VECTOR, int spacedim>
+ void project (const hp::DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const hp::QCollection<dim> &quadrature,
+ const Function<spacedim> &function,
+ VECTOR &vec,
+ const bool enforce_zero_boundary = false,
+ const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
+ QGauss<dim-1>(2) :
+ Quadrature<dim-1>(0)),
+ const bool project_to_boundary_first = false);
+
/**
* Compute Dirichlet boundary
* conditions. This function makes up a map
std::map<types::global_dof_index,double> &boundary_values,
const ComponentMask &component_mask = ComponentMask());
+ /**
+ * Like the previous function, but take a mapping collection to go with
+ * the hp::DoFHandler object.
+ */
+ template <int dim, int spacedim>
+ void
+ interpolate_boundary_values (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &function_map,
+ std::map<types::global_dof_index,double> &boundary_values,
+ const ComponentMask &component_mask = ComponentMask());
+
/**
* Same function as above, but
* taking only one pair of
std::map<types::global_dof_index,double> &boundary_values,
std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
+ /**
+ * Same as above, but for objects of type hp::DoFHandler
+ */
+ template <int dim, int spacedim>
+ void project_boundary_values (const hp::MappingCollection<dim, spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_functions,
+ const hp::QCollection<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
+
+ /**
+ * Calls the project_boundary_values()
+ * function, see above, with
+ * <tt>mapping=MappingQ1@<dim@>()</tt>.
+ */
+ template <int dim, int spacedim>
+ void project_boundary_values (const hp::DoFHandler<dim,spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_function,
+ const hp::QCollection<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
+
/**
* Project a function to the boundary of
* the domain, using the given quadrature
namespace internal
{
+ /**
+ * Interpolate zero boundary values. We don't need to worry about a mapping
+ * here because the function we evaluate for the DoFs is zero in the mapped
+ * locations as well as in the original, unmapped locations
+ */
+ template <class DH>
void
- interpolate_zero_boundary_values (const dealii::DoFHandler<1> &dof_handler,
- std::map<types::global_dof_index,double> &boundary_values)
- {
- // we only need to find the
- // left-most and right-most
- // vertex and query its vertex
- // dof indices. that's easy :-)
- for (unsigned int direction=0; direction<2; ++direction)
- {
- dealii::DoFHandler<1>::cell_iterator
- cell = dof_handler.begin(0);
- while (!cell->at_boundary(direction))
- cell = cell->neighbor(direction);
-
- for (unsigned int i=0; i<dof_handler.get_fe().dofs_per_vertex; ++i)
- boundary_values[cell->vertex_dof_index (direction, i)] = 0.;
- }
- }
-
-
-
- // codimension 1
- void
- interpolate_zero_boundary_values (const dealii::DoFHandler<1,2> &dof_handler,
- std::map<types::global_dof_index,double> &boundary_values)
- {
- // we only need to find the
- // left-most and right-most
- // vertex and query its vertex
- // dof indices. that's easy :-)
- for (unsigned int direction=0; direction<2; ++direction)
- {
- dealii::DoFHandler<1,2>::cell_iterator
- cell = dof_handler.begin(0);
- while (!cell->at_boundary(direction))
- cell = cell->neighbor(direction);
-
- for (unsigned int i=0; i<dof_handler.get_fe().dofs_per_vertex; ++i)
- boundary_values[cell->vertex_dof_index (direction, i)] = 0.;
- }
- }
-
-
-
- template <int dim, int spacedim>
- void
- interpolate_zero_boundary_values (const dealii::DoFHandler<dim,spacedim> &dof_handler,
+ interpolate_zero_boundary_values (const DH &dof_handler,
std::map<types::global_dof_index,double> &boundary_values)
{
- const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
+ const unsigned int dim = DH::dimension;
+ const unsigned int spacedim = DH::space_dimension;
// loop over all boundary faces
// to get all dof indices of
// that is actually wholly on
// the boundary, not only by
// one line or one vertex
- typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+ typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- std::vector<types::global_dof_index> face_dof_indices (fe.dofs_per_face);
+ std::vector<types::global_dof_index> face_dof_indices;
for (; cell!=endc; ++cell)
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->at_boundary(f))
{
- cell->face(f)->get_dof_indices (face_dof_indices);
- for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+ face_dof_indices.resize (cell->get_fe().dofs_per_face);
+ cell->face(f)->get_dof_indices (face_dof_indices,
+ cell->active_fe_index());
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
// enter zero boundary values
// for all boundary nodes
//
/**
* Compute the boundary values to be used in the project() functions.
*/
- template <int dim, int spacedim>
- void project_compute_b_v (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim,spacedim> &dof,
- const Function<spacedim> &function,
- const bool enforce_zero_boundary,
- const Quadrature<dim-1> &q_boundary,
- const bool project_to_boundary_first,
- std::map<types::global_dof_index,double> &boundary_values)
+ template <int dim, int spacedim,
+ template <int,int> class DH,
+ template <int,int> class M_or_MC,
+ template <int> class Q_or_QC>
+ void project_compute_b_v (const M_or_MC<dim, spacedim> &mapping,
+ const DH<dim,spacedim> &dof,
+ const Function<spacedim> &function,
+ const bool enforce_zero_boundary,
+ const Q_or_QC<dim-1> &q_boundary,
+ const bool project_to_boundary_first,
+ std::map<types::global_dof_index,double> &boundary_values)
{
if (enforce_zero_boundary == true)
// no need to project boundary
return true;
}
+
+
+ /**
+ * Generic implementation of the project() function
+ */
+ template <int dim, int spacedim,
+ class Vector,
+ template <int,int> class DH,
+ template <int,int> class M_or_MC,
+ template <int> class Q_or_QC>
+ void do_project (const M_or_MC<dim, spacedim> &mapping,
+ const DH<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Q_or_QC<dim> &quadrature,
+ const Function<spacedim> &function,
+ Vector &vec_result,
+ const bool enforce_zero_boundary,
+ const Q_or_QC<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ Assert (dof.get_fe().n_components() == function.n_components,
+ ExcDimensionMismatch(dof.get_fe().n_components(),
+ function.n_components));
+ Assert (vec_result.size() == dof.n_dofs(),
+ ExcDimensionMismatch (vec_result.size(), dof.n_dofs()));
+
+ // make up boundary values
+ std::map<types::global_dof_index,double> boundary_values;
+ project_compute_b_v(mapping, dof, function, enforce_zero_boundary,
+ q_boundary, project_to_boundary_first, boundary_values);
+
+ // check if constraints are compatible (see below)
+ const bool constraints_are_compatible =
+ constraints_and_b_v_are_compatible(constraints, boundary_values);
+
+ // set up mass matrix and right hand side
+ dealii::Vector<double> vec (dof.n_dofs());
+ SparsityPattern sparsity;
+ {
+ CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs());
+ DoFTools::make_sparsity_pattern (dof, csp, constraints,
+ !constraints_are_compatible);
+
+ sparsity.copy_from (csp);
+ }
+ SparseMatrix<double> mass_matrix (sparsity);
+ dealii::Vector<double> tmp (mass_matrix.n());
+
+ // If the constraint matrix does not conflict with the given boundary
+ // values (i.e., it either does not contain boundary values or it contains
+ // the same as boundary_values), we can let it call
+ // distribute_local_to_global straight away, otherwise we need to first
+ // interpolate the boundary values and then condense the matrix and vector
+ if (constraints_are_compatible)
+ {
+ const Function<spacedim> *dummy = 0;
+ MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
+ mass_matrix, function, tmp,
+ dummy, constraints);
+ if (boundary_values.size() > 0)
+ MatrixTools::apply_boundary_values (boundary_values,
+ mass_matrix, vec, tmp,
+ true);
+ }
+ else
+ {
+ // create mass matrix and rhs at once, which is faster.
+ MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
+ mass_matrix, function, tmp);
+ MatrixTools::apply_boundary_values (boundary_values,
+ mass_matrix, vec, tmp,
+ true);
+ constraints.condense(mass_matrix, tmp);
+ }
+
+ // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
+ // steps may not be sufficient, since roundoff errors may accumulate for
+ // badly conditioned matrices
+ ReductionControl control(5*tmp.size(), 0., 1e-12, false, false);
+ GrowingVectorMemory<> memory;
+ SolverCG<> cg(control,memory);
+
+ PreconditionSSOR<> prec;
+ prec.initialize(mass_matrix, 1.2);
+
+ cg.solve (mass_matrix, vec, tmp, prec);
+
+ constraints.distribute (vec);
+
+ // copy vec into vec_result. we can't use vec_result itself above, since
+ // it may be of another type than Vector<double> and that wouldn't
+ // necessarily go together with the matrix and other functions
+ for (unsigned int i=0; i<vec.size(); ++i)
+ vec_result(i) = vec(i);
+ }
}
const Quadrature<dim-1> &q_boundary,
const bool project_to_boundary_first)
{
- Assert (dof.get_fe().n_components() == function.n_components,
- ExcDimensionMismatch(dof.get_fe().n_components(),
- function.n_components));
-
- Assert (vec_result.size() == dof.n_dofs(),
- ExcDimensionMismatch (vec_result.size(), dof.n_dofs()));
-
- // make up boundary values
- std::map<types::global_dof_index,double> boundary_values;
- project_compute_b_v(mapping, dof, function, enforce_zero_boundary,
- q_boundary, project_to_boundary_first, boundary_values);
-
- // check if constraints are compatible (see below)
- const bool constraints_are_compatible =
- constraints_and_b_v_are_compatible(constraints, boundary_values);
-
- // set up mass matrix and right hand side
- Vector<double> vec (dof.n_dofs());
- SparsityPattern sparsity;
- {
- CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs());
- DoFTools::make_sparsity_pattern (dof, csp, constraints,
- !constraints_are_compatible);
-
- sparsity.copy_from (csp);
- }
- SparseMatrix<double> mass_matrix (sparsity);
- Vector<double> tmp (mass_matrix.n());
-
- // If the constraint matrix does not conflict with the given boundary
- // values (i.e., it either does not contain boundary values or it contains
- // the same as boundary_values), we can let it call
- // distribute_local_to_global straight away, otherwise we need to first
- // interpolate the boundary values and then condense the matrix and vector
- if (constraints_are_compatible)
- {
- const Function<spacedim> *dummy = 0;
- MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
- mass_matrix, function, tmp,
- dummy, constraints);
- if (boundary_values.size() > 0)
- MatrixTools::apply_boundary_values (boundary_values,
- mass_matrix, vec, tmp,
- true);
- }
- else
- {
- // create mass matrix and rhs at once, which is faster.
- MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
- mass_matrix, function, tmp);
- MatrixTools::apply_boundary_values (boundary_values,
- mass_matrix, vec, tmp,
- true);
- constraints.condense(mass_matrix, tmp);
- }
-
- // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
- // steps may not be sufficient, since roundoff errors may accumulate for
- // badly conditioned matrices
- ReductionControl control(5*tmp.size(), 0., 1e-12, false, false);
- GrowingVectorMemory<> memory;
- SolverCG<> cg(control,memory);
-
- PreconditionSSOR<> prec;
- prec.initialize(mass_matrix, 1.2);
-
- cg.solve (mass_matrix, vec, tmp, prec);
-
- constraints.distribute (vec);
-
- // copy vec into vec_result. we can't use vec_result itself above, since
- // it may be of another type than Vector<double> and that wouldn't
- // necessarily go together with the matrix and other functions
- for (unsigned int i=0; i<vec.size(); ++i)
- vec_result(i) = vec(i);
+ do_project (mapping, dof, constraints, quadrature,
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
}
+ template <int dim, class VECTOR, int spacedim>
+ void project (const hp::MappingCollection<dim, spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const hp::QCollection<dim> &quadrature,
+ const Function<spacedim> &function,
+ VECTOR &vec_result,
+ const bool enforce_zero_boundary,
+ const hp::QCollection<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ do_project (mapping, dof, constraints, quadrature,
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
+ }
+
+
+ template <int dim, class VECTOR, int spacedim>
+ void project (const hp::DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const hp::QCollection<dim> &quadrature,
+ const Function<spacedim> &function,
+ VECTOR &vec,
+ const bool enforce_zero_boundary,
+ const hp::QCollection<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
+ {
+ project(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
+ dof, constraints, quadrature, function, vec,
+ enforce_zero_boundary, q_boundary, project_to_boundary_first);
+ }
+
template <int dim, int spacedim>
void create_right_hand_side (const Mapping<dim, spacedim> &mapping,
// ----------- interpolate_boundary_values for std::map --------------------
- namespace internal
+ namespace
{
// interpolate boundary values in
// 1d. in higher dimensions, we
// faces are points and it is far
// easier to simply work on
// individual vertices
- template <class DH>
+ template <class DH,
+ template <int,int> class M_or_MC>
static inline
- void interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension> &,
+ void do_interpolate_boundary_values (const M_or_MC<DH::dimension, DH::space_dimension> &,
const DH &dof,
const typename FunctionMap<DH::space_dimension>::type &function_map,
std::map<types::global_dof_index,double> &boundary_values,
- template <class DH>
+ template <class DH,
+ template <int,int> class M_or_MC>
static inline
void
- interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension> &mapping,
+ do_interpolate_boundary_values (const M_or_MC<DH::dimension, DH::space_dimension> &mapping,
const DH &dof,
const typename FunctionMap<DH::space_dimension>::type &function_map,
std::map<types::global_dof_index,double> &boundary_values,
std::map<types::global_dof_index,double> &boundary_values,
const ComponentMask &component_mask_)
{
- internal::
- interpolate_boundary_values (mapping, dof, function_map, boundary_values,
+ do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
component_mask_,
dealii::internal::int2type<DH::dimension>());
}
}
+ template <int dim, int spacedim>
+ void
+ interpolate_boundary_values (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &function_map,
+ std::map<types::global_dof_index,double> &boundary_values,
+ const ComponentMask &component_mask_)
+ {
+ do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
+ component_mask_,
+ dealii::internal::int2type<dim>());
+ }
+
+
template <class DH>
void
// -------- implementation for project_boundary_values with std::map --------
-// separate implementations for 1D because otherwise we get linking errors
-// because create_boundary_mass_matrix is not compiled for 1D
- template <>
- void
- project_boundary_values (const Mapping<1,1> &mapping,
- const DoFHandler<1,1> &dof,
- const FunctionMap<1>::type &boundary_functions,
- const Quadrature<0> &,
- std::map<types::global_dof_index,double> &boundary_values,
- std::vector<unsigned int> component_mapping)
- {
- Assert (component_mapping.size() == 0, ExcNotImplemented());
- // projection in 1d is equivalent
- // to interpolation
- interpolate_boundary_values (mapping, dof, boundary_functions,
- boundary_values, ComponentMask());
- }
-
-
- template <>
- void
- project_boundary_values (const Mapping<1,2> &mapping,
- const DoFHandler<1,2> &dof,
- const FunctionMap<2>::type &boundary_functions,
- const Quadrature<0> &,
- std::map<types::global_dof_index,double> &boundary_values,
- std::vector<unsigned int> component_mapping)
+ namespace
{
- Assert (component_mapping.size() == 0, ExcNotImplemented());
- // projection in 1d is equivalent
- // to interpolation
- interpolate_boundary_values (mapping, dof, boundary_functions,
- boundary_values, ComponentMask());
- }
+ template <int dim, int spacedim,
+ template <int,int> class DH,
+ template <int,int> class M_or_MC,
+ template <int> class Q_or_QC>
+ void
+ do_project_boundary_values (const M_or_MC<dim, spacedim> &mapping,
+ const DH<dim, spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_functions,
+ const Q_or_QC<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
+ {
+ // in 1d, projection onto the 0d end points == interpolation
+ if (dim == 1)
+ {
+ Assert (component_mapping.size() == 0, ExcNotImplemented());
+ interpolate_boundary_values (mapping, dof, boundary_functions,
+ boundary_values, ComponentMask());
+ return;
+ }
+ //TODO:[?] In project_boundary_values, no condensation of sparsity
+ // structures, matrices and right hand sides or distribution of
+ // solution vectors is performed. This is ok for dim<3 because then
+ // there are no constrained nodes on the boundary, but is not
+ // acceptable for higher dimensions. Fix this.
+ if (component_mapping.size() == 0)
+ {
+ AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components);
+ // I still do not see why i
+ // should create another copy
+ // here
+ component_mapping.resize(dof.get_fe().n_components());
+ for (unsigned int i= 0 ; i < component_mapping.size() ; ++i)
+ component_mapping[i] = i;
+ }
+ else
+ AssertDimension (dof.get_fe().n_components(), component_mapping.size());
+ std::vector<types::global_dof_index> dof_to_boundary_mapping;
+ std::set<types::boundary_id> selected_boundary_components;
+ for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
+ i!=boundary_functions.end(); ++i)
+ selected_boundary_components.insert (i->first);
- template <int dim, int spacedim>
- void
- project_boundary_values (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim, spacedim> &dof,
- const typename FunctionMap<spacedim>::type &boundary_functions,
- const Quadrature<dim-1> &q,
- std::map<types::global_dof_index,double> &boundary_values,
- std::vector<unsigned int> component_mapping)
- {
-//TODO:[?] In project_boundary_values, no condensation of sparsity
-// structures, matrices and right hand sides or distribution of
-// solution vectors is performed. This is ok for dim<3 because then
-// there are no constrained nodes on the boundary, but is not
-// acceptable for higher dimensions. Fix this.
+ DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
+ dof_to_boundary_mapping);
- if (component_mapping.size() == 0)
- {
- AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components);
- // I still do not see why i
- // should create another copy
- // here
- component_mapping.resize(dof.get_fe().n_components());
- for (unsigned int i= 0 ; i < component_mapping.size() ; ++i)
- component_mapping[i] = i;
- }
- else
- AssertDimension (dof.get_fe().n_components(), component_mapping.size());
-
- std::vector<types::global_dof_index> dof_to_boundary_mapping;
- std::set<types::boundary_id> selected_boundary_components;
- for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
- i!=boundary_functions.end(); ++i)
- selected_boundary_components.insert (i->first);
-
- DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
- dof_to_boundary_mapping);
-
- // Done if no degrees of freedom on the boundary
- if (dof.n_boundary_dofs (boundary_functions) == 0)
- return;
-
- // set up sparsity structure
- CompressedSparsityPattern c_sparsity(dof.n_boundary_dofs (boundary_functions),
- dof.n_boundary_dofs (boundary_functions));
- DoFTools::make_boundary_sparsity_pattern (dof,
- boundary_functions,
- dof_to_boundary_mapping,
- c_sparsity);
- SparsityPattern sparsity;
- sparsity.copy_from(c_sparsity);
-
-
-
- // note: for three or more dimensions, there
- // may be constrained nodes on the boundary
- // in this case the boundary mass matrix has
- // to be condensed and the solution is to
- // be distributed afterwards, which is not
- // yet implemented. The reason for this is
- // that we cannot simply use the condense
- // family of functions, since the matrices
- // and vectors do not use the global
- // numbering but rather the boundary
- // numbering, i.e. the condense function
- // needs to use another indirection. There
- // should be not many technical problems,
- // but it needs to be implemented
- if (dim>=3)
- {
-#ifdef DEBUG
-// Assert that there are no hanging nodes at the boundary
- int level = -1;
- for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
- cell != dof.end(); ++cell)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- {
- if (cell->at_boundary(f))
- {
- if (level == -1)
- level = cell->level();
- else
- {
- Assert (level == cell->level(), ExcNotImplemented());
- }
- }
- }
-#endif
- }
- sparsity.compress();
+ // Done if no degrees of freedom on the boundary
+ if (dof.n_boundary_dofs (boundary_functions) == 0)
+ return;
+
+ // set up sparsity structure
+ CompressedSparsityPattern c_sparsity(dof.n_boundary_dofs (boundary_functions),
+ dof.n_boundary_dofs (boundary_functions));
+ DoFTools::make_boundary_sparsity_pattern (dof,
+ boundary_functions,
+ dof_to_boundary_mapping,
+ c_sparsity);
+ SparsityPattern sparsity;
+ sparsity.copy_from(c_sparsity);
+
+
+
+ // note: for three or more dimensions, there
+ // may be constrained nodes on the boundary
+ // in this case the boundary mass matrix has
+ // to be condensed and the solution is to
+ // be distributed afterwards, which is not
+ // yet implemented. The reason for this is
+ // that we cannot simply use the condense
+ // family of functions, since the matrices
+ // and vectors do not use the global
+ // numbering but rather the boundary
+ // numbering, i.e. the condense function
+ // needs to use another indirection. There
+ // should be not many technical problems,
+ // but it needs to be implemented
+ if (dim>=3)
+ {
+ #ifdef DEBUG
+ // Assert that there are no hanging nodes at the boundary
+ int level = -1;
+ for (typename DH<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
+ cell != dof.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ {
+ if (cell->at_boundary(f))
+ {
+ if (level == -1)
+ level = cell->level();
+ else
+ {
+ Assert (level == cell->level(), ExcNotImplemented());
+ }
+ }
+ }
+ #endif
+ }
+ sparsity.compress();
- // make mass matrix and right hand side
- SparseMatrix<double> mass_matrix(sparsity);
- Vector<double> rhs(sparsity.n_rows());
+ // make mass matrix and right hand side
+ SparseMatrix<double> mass_matrix(sparsity);
+ Vector<double> rhs(sparsity.n_rows());
- MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
- mass_matrix, boundary_functions,
- rhs, dof_to_boundary_mapping, (const Function<spacedim> *) 0,
- component_mapping);
+ MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
+ mass_matrix, boundary_functions,
+ rhs, dof_to_boundary_mapping, (const Function<spacedim> *) 0,
+ component_mapping);
- // For certain weird elements,
- // there might be degrees of
- // freedom on the boundary, but
- // their shape functions do not
- // have support there. Let's
- // eliminate them here.
+ // For certain weird elements,
+ // there might be degrees of
+ // freedom on the boundary, but
+ // their shape functions do not
+ // have support there. Let's
+ // eliminate them here.
- // The Bogner-Fox-Schmidt element
- // is an example for those.
+ // The Bogner-Fox-Schmidt element
+ // is an example for those.
//TODO: Maybe we should figure out if the element really needs this
- FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
- FilteredMatrix<Vector<double> > filtered_precondition;
- std::vector<bool> excluded_dofs(mass_matrix.m(), false);
+ FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
+ FilteredMatrix<Vector<double> > filtered_precondition;
+ std::vector<bool> excluded_dofs(mass_matrix.m(), false);
- double max_element = 0.;
- for (unsigned int i=0; i<mass_matrix.m(); ++i)
- if (mass_matrix.diag_element(i) > max_element)
- max_element = mass_matrix.diag_element(i);
+ double max_element = 0.;
+ for (unsigned int i=0; i<mass_matrix.m(); ++i)
+ if (mass_matrix.diag_element(i) > max_element)
+ max_element = mass_matrix.diag_element(i);
- for (unsigned int i=0; i<mass_matrix.m(); ++i)
- if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
- {
- filtered_mass_matrix.add_constraint(i, 0.);
- filtered_precondition.add_constraint(i, 0.);
- mass_matrix.diag_element(i) = 1.;
- excluded_dofs[i] = true;
- }
+ for (unsigned int i=0; i<mass_matrix.m(); ++i)
+ if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
+ {
+ filtered_mass_matrix.add_constraint(i, 0.);
+ filtered_precondition.add_constraint(i, 0.);
+ mass_matrix.diag_element(i) = 1.;
+ excluded_dofs[i] = true;
+ }
- Vector<double> boundary_projection (rhs.size());
-
- // Allow for a maximum of 5*n
- // steps to reduce the residual by
- // 10^-12. n steps may not be
- // sufficient, since roundoff
- // errors may accumulate for badly
- // conditioned matrices
- ReductionControl control(5*rhs.size(), 0., 1.e-12, false, false);
- GrowingVectorMemory<> memory;
- SolverCG<> cg(control,memory);
-
- PreconditionSSOR<> prec;
- prec.initialize(mass_matrix, 1.2);
- filtered_precondition.initialize(prec, true);
- // solve
- cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
- filtered_precondition.apply_constraints(boundary_projection, true);
- filtered_precondition.clear();
- // fill in boundary values
- for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
- if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
- && ! excluded_dofs[dof_to_boundary_mapping[i]])
- {
- Assert(numbers::is_finite(boundary_projection(dof_to_boundary_mapping[i])), ExcNumberNotFinite());
+ Vector<double> boundary_projection (rhs.size());
+
+ // Allow for a maximum of 5*n
+ // steps to reduce the residual by
+ // 10^-12. n steps may not be
+ // sufficient, since roundoff
+ // errors may accumulate for badly
+ // conditioned matrices
+ ReductionControl control(5*rhs.size(), 0., 1.e-12, false, false);
+ GrowingVectorMemory<> memory;
+ SolverCG<> cg(control,memory);
+
+ PreconditionSSOR<> prec;
+ prec.initialize(mass_matrix, 1.2);
+ filtered_precondition.initialize(prec, true);
+ // solve
+ cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
+ filtered_precondition.apply_constraints(boundary_projection, true);
+ filtered_precondition.clear();
+ // fill in boundary values
+ for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
+ if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
+ && ! excluded_dofs[dof_to_boundary_mapping[i]])
+ {
+ Assert(numbers::is_finite(boundary_projection(dof_to_boundary_mapping[i])), ExcNumberNotFinite());
- // this dof is on one of the
- // interesting boundary parts
- //
- // remember: i is the global dof
- // number, dof_to_boundary_mapping[i]
- // is the number on the boundary and
- // thus in the solution vector
- boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
- }
+ // this dof is on one of the
+ // interesting boundary parts
+ //
+ // remember: i is the global dof
+ // number, dof_to_boundary_mapping[i]
+ // is the number on the boundary and
+ // thus in the solution vector
+ boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
+ }
+ }
+ }
+
+ template <int dim, int spacedim>
+ void
+ project_boundary_values (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim, spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_functions,
+ const Quadrature<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
+ {
+ do_project_boundary_values(mapping, dof, boundary_functions, q,
+ boundary_values, component_mapping);
}
-// ----- implementation for project_boundary_values with ConstraintMatrix -----
+ template <int dim, int spacedim>
+ void project_boundary_values (const hp::MappingCollection<dim, spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_functions,
+ const hp::QCollection<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
+ {
+ do_project_boundary_values (hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof,
+ boundary_functions,
+ q, boundary_values,
+ component_mapping);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void project_boundary_values (const hp::DoFHandler<dim,spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_function,
+ const hp::QCollection<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
+ {
+ project_boundary_values (hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof,
+ boundary_function,
+ q, boundary_values,
+ component_mapping);
+ }
+
+
+ // ----- implementation for project_boundary_values with ConstraintMatrix -----
const hp::DoFHandler<deal_II_dimension,deal_II_dimension> &dof_col,
SP &sparsity);
- // #if deal_II_dimension > 1
template void
DoFTools::make_boundary_sparsity_pattern<DoFHandler<deal_II_dimension>,SP>
(const DoFHandler<deal_II_dimension>& dof,
const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
SP &sparsity);
+#if deal_II_dimension < 3
+ template void
+ DoFTools::make_boundary_sparsity_pattern<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1>,SP>
+ (const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1>& dof,
+ const FunctionMap<deal_II_dimension+1>::type &boundary_indicators,
+ const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+ SP &sparsity);
+ #endif
+
template void
DoFTools::make_flux_sparsity_pattern<DoFHandler<deal_II_dimension>,SP>
(const DoFHandler<deal_II_dimension> &dof,
UpdateFlags update_flags = UpdateFlags (update_values |
update_JxW_values |
update_quadrature_points);
- hp::FEFaceValues<dim> x_fe_values (mapping, fe_collection, q, update_flags);
+ hp::FEFaceValues<dim,spacedim> x_fe_values (mapping, fe_collection, q, update_flags);
// two variables for the coefficient,
// one for the two cases indicated in
const Function<deal_II_space_dimension> * const a,
std::vector<unsigned int>);
+ template
+ void
+ MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
+ const hp::QCollection<deal_II_dimension-1>&,
+ SparseMatrix<double>&,
+ const FunctionMap<deal_II_space_dimension>::type&,
+ Vector<double>&,
+ std::vector<types::global_dof_index>&,
+ const Function<deal_II_space_dimension> * const,
+ std::vector<unsigned int>);
-
+ template
+ void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
+ const hp::QCollection<deal_II_dimension-1>&,
+ SparseMatrix<double>&,
+ const FunctionMap<deal_II_space_dimension>::type&,
+ Vector<double>&,
+ std::vector<types::global_dof_index>&,
+ const Function<deal_II_space_dimension> * const,
+ std::vector<unsigned int>);
#endif
}
// where applicable and move to codimension cases above also when applicable
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
{
-#if deal_II_dimension == deal_II_space_dimension
-
+#if deal_II_dimension <= deal_II_space_dimension
// hp versions of functions
template
- void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::MappingCollection<deal_II_dimension> &mapping,
- const hp::DoFHandler<deal_II_dimension> &dof,
+ void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+
template
- void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::DoFHandler<deal_II_dimension> &dof,
+ void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient,
+ const Function<deal_II_space_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::MappingCollection<deal_II_dimension> &mapping,
- const hp::DoFHandler<deal_II_dimension> &dof,
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
- SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> &rhs,
+ SparseMatrix<float> &matrix,
+ const Function<deal_II_space_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::DoFHandler<deal_II_dimension> &dof,
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
- SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> &rhs,
+ SparseMatrix<float> &matrix,
+ const Function<deal_II_space_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+#endif
+#if deal_II_dimension == deal_II_space_dimension
+
+
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::MappingCollection<deal_II_dimension> &mapping,
- const hp::DoFHandler<deal_II_dimension> &dof,
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
- SparseMatrix<float> &matrix,
- const Function<deal_II_dimension> * const coefficient,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::DoFHandler<deal_II_dimension> &dof,
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
- SparseMatrix<float> &matrix,
- const Function<deal_II_dimension> * const coefficient,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_space_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::MappingCollection<deal_II_dimension> &mapping,
- const hp::DoFHandler<deal_II_dimension> &dof,
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<float> &matrix,
- const Function<deal_II_dimension> &rhs,
- Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient,
+ const Function<deal_II_space_dimension> * const coefficient,
const ConstraintMatrix &constraints);
+
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
- (const hp::DoFHandler<deal_II_dimension> &dof,
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<float> &matrix,
- const Function<deal_II_dimension> &rhs,
- Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient,
- const ConstraintMatrix &constraints);
-
-
-#if deal_II_dimension != 1
-
-
-
- template
- void
- MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
- (const hp::MappingCollection<deal_II_dimension>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const hp::QCollection<deal_II_dimension-1>&,
- SparseMatrix<double>&,
- const FunctionMap<deal_II_dimension>::type&,
- Vector<double>&,
- std::vector<types::global_dof_index>&,
- const Function<deal_II_dimension> * const,
- std::vector<unsigned int>);
-
- template
- void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
- (const hp::DoFHandler<deal_II_dimension>&,
- const hp::QCollection<deal_II_dimension-1>&,
- SparseMatrix<double>&,
- const FunctionMap<deal_II_dimension>::type&,
- Vector<double>&,
- std::vector<types::global_dof_index>&,
- const Function<deal_II_dimension> * const,
- std::vector<unsigned int>);
-
-#endif
-
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
// non-hp versions of create_laplace_matrix
template
const bool,
const Quadrature<deal_II_dimension-1> &,
const bool);
+
+ template
+ void project
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+ const ConstraintMatrix &,
+ const hp::QCollection<deal_II_dimension> &,
+ const Function<deal_II_space_dimension> &,
+ VEC &,
+ const bool,
+ const hp::QCollection<deal_II_dimension-1> &,
+ const bool);
+
+ template
+ void project
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+ const ConstraintMatrix &,
+ const hp::QCollection<deal_II_dimension> &,
+ const Function<deal_II_space_dimension> &,
+ VEC &,
+ const bool,
+ const hp::QCollection<deal_II_dimension-1> &,
+ const bool);
\}
#endif
}
// where applicable and move to codimension cases above also when applicable
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
{
-#if deal_II_dimension == deal_II_space_dimension
namespace VectorTools \{
+
+#if deal_II_dimension <= deal_II_space_dimension
+ template
+ void project_boundary_values<deal_II_dimension,deal_II_space_dimension>
+ (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &,
+ const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+ const FunctionMap<deal_II_space_dimension>::type &,
+ const hp::QCollection<deal_II_dimension-1>&,
+ std::map<types::global_dof_index,double>&,
+ std::vector<unsigned int>);
+
+ template
+ void project_boundary_values<deal_II_dimension,deal_II_space_dimension>
+ (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+ const FunctionMap<deal_II_space_dimension>::type &,
+ const hp::QCollection<deal_II_dimension-1>&,
+ std::map<types::global_dof_index,double>&,
+ std::vector<unsigned int>);
+#endif
+
+#if deal_II_dimension == deal_II_space_dimension
template
void create_right_hand_side<deal_II_dimension>
(const hp::MappingCollection<deal_II_dimension> &,
-#if deal_II_dimension > 1
template
void project_boundary_values<deal_II_dimension>
(const Mapping<deal_II_dimension> &,
const FunctionMap<deal_II_dimension>::type &,
const Quadrature<deal_II_dimension-1>&,
std::map<types::global_dof_index,double>&, std::vector<unsigned int>);
-#endif
template
void project_boundary_values<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &,
const FunctionMap<deal_II_dimension>::type &,
const Quadrature<deal_II_dimension-1>&,
- std::map<types::global_dof_index,double>&, std::vector<unsigned int>);
+ std::map<types::global_dof_index,double>&,
+ std::vector<unsigned int>);
- template
- void project_boundary_values<deal_II_dimension>
- (const Mapping<deal_II_dimension> &,
- const DoFHandler<deal_II_dimension> &,
- const FunctionMap<deal_II_dimension>::type &,
- const Quadrature<deal_II_dimension-1>&,
- ConstraintMatrix&, std::vector<unsigned int>);
-
- template
- void project_boundary_values<deal_II_dimension>
- (const DoFHandler<deal_II_dimension> &,
- const FunctionMap<deal_II_dimension>::type &,
- const Quadrature<deal_II_dimension-1>&,
- ConstraintMatrix&, std::vector<unsigned int>);
-
#if deal_II_dimension != 1
template
ConstraintMatrix &constraints,
const Mapping<deal_II_dimension> &mapping);
#endif
- \}
#endif
+ \}
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2006 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::project for Vector<double> arguments and hp
+
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <fstream>
+
+
+// define the multi-linear function x or x*y or x*y*z that we will
+// subsequently project onto the ansatz space
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int = 0) const
+ {
+ double s = 1;
+ for (unsigned int i=0; i<dim; ++i)
+ s *= p[i];
+ return s;
+ }
+};
+
+
+template<int dim>
+void test()
+{
+ // create 2 triangulations with the
+ // same coarse grid, and refine
+ // them differently
+ Triangulation<dim> tria;
+
+ GridGenerator::hyper_cube (tria);
+ tria.refine_global (2);
+
+ hp::FECollection<dim> fe;
+ fe.push_back(FE_Q<dim>(1));
+ hp::DoFHandler<dim> dh (tria);
+ dh.distribute_dofs (fe);
+
+ Vector<double> v(dh.n_dofs());
+
+ ConstraintMatrix cm;
+ cm.close ();
+ VectorTools::project (dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v);
+
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=dh.begin_active();
+ cell != dh.end(); ++cell)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ // check that the error is
+ // somewhat small. it won't
+ // be zero since we project
+ // and do not interpolate
+ Assert (std::fabs (v(cell->vertex_dof_index(i,0,cell->active_fe_index())) -
+ F<dim>().value (cell->vertex(i)))
+ < 1e-4,
+ ExcInternalError());
+ deallog << cell->vertex(i) << ' ' << v(cell->vertex_dof_index(i,0,cell->active_fe_index()))
+ << std::endl;
+ }
+}
+
+
+int main()
+{
+ std::ofstream logfile ("project_01/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1>();
+ test<2>();
+ test<3>();
+}
+
--- /dev/null
+
+DEAL::0.00000 0
+DEAL::0.250000 0.250000
+DEAL::0.250000 0.250000
+DEAL::0.500000 0.500000
+DEAL::0.500000 0.500000
+DEAL::0.750000 0.750000
+DEAL::0.750000 0.750000
+DEAL::1.00000 1.00000
+DEAL::0.00000 0.00000 0
+DEAL::0.250000 0.00000 0
+DEAL::0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.250000 0.00000 0
+DEAL::0.500000 0.00000 0
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.00000 0.500000 0
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.500000 0.00000 0
+DEAL::0.750000 0.00000 0
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.750000 0.250000 0.187500
+DEAL::0.750000 0.00000 0
+DEAL::1.00000 0.00000 0
+DEAL::0.750000 0.250000 0.187500
+DEAL::1.00000 0.250000 0.250000
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.750000 0.250000 0.187500
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.750000 0.500000 0.375000
+DEAL::0.750000 0.250000 0.187500
+DEAL::1.00000 0.250000 0.250000
+DEAL::0.750000 0.500000 0.375000
+DEAL::1.00000 0.500000 0.500000
+DEAL::0.00000 0.500000 0
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.00000 0.750000 0
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.00000 0.750000 0
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.00000 1.00000 0
+DEAL::0.250000 1.00000 0.250000
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.250000 1.00000 0.250000
+DEAL::0.500000 1.00000 0.500000
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.750000 0.500000 0.375000
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.562500
+DEAL::0.750000 0.500000 0.375000
+DEAL::1.00000 0.500000 0.500000
+DEAL::0.750000 0.750000 0.562500
+DEAL::1.00000 0.750000 0.750000
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.562500
+DEAL::0.500000 1.00000 0.500000
+DEAL::0.750000 1.00000 0.750000
+DEAL::0.750000 0.750000 0.562500
+DEAL::1.00000 0.750000 0.750000
+DEAL::0.750000 1.00000 0.750000
+DEAL::1.00000 1.00000 1.00000
+DEAL::0.00000 0.00000 0.00000 0
+DEAL::0.250000 0.00000 0.00000 0
+DEAL::0.00000 0.250000 0.00000 0
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.00000 0.00000 0.250000 0
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.250000 0.00000 0.00000 0
+DEAL::0.500000 0.00000 0.00000 0
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.00000 0.250000 0.00000 0
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.00000 0.500000 0.00000 0
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.00000 0.00000 0.250000 0
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.00000 0.00000 0.500000 0
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.500000 0.00000 0.00000 0
+DEAL::0.750000 0.00000 0.00000 0
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.750000 0.00000 0.00000 0
+DEAL::1.00000 0.00000 0.00000 0
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::1.00000 0.250000 0.00000 0
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::1.00000 0.00000 0.250000 0
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::1.00000 0.250000 0.00000 0
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::1.00000 0.500000 0.00000 0
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::1.00000 0.00000 0.250000 0
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::1.00000 0.00000 0.500000 0
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.00000 0.500000 0.00000 0
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.00000 0.750000 0.00000 0
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.00000 0.750000 0.00000 0
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.00000 1.00000 0.00000 0
+DEAL::0.250000 1.00000 0.00000 0
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.00000 1.00000 0.250000 0
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.250000 1.00000 0.00000 0
+DEAL::0.500000 1.00000 0.00000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.00000 1.00000 0.250000 0
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.00000 1.00000 0.500000 0
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::1.00000 0.500000 0.00000 0
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::1.00000 0.750000 0.00000 0
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::0.500000 1.00000 0.00000 0
+DEAL::0.750000 1.00000 0.00000 0
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::1.00000 0.750000 0.00000 0
+DEAL::0.750000 1.00000 0.00000 0
+DEAL::1.00000 1.00000 0.00000 0
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::1.00000 1.00000 0.250000 0.250000
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::1.00000 1.00000 0.250000 0.250000
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::1.00000 1.00000 0.500000 0.500000
+DEAL::0.00000 0.00000 0.500000 0
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.00000 0.00000 0.750000 0
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.00000 0.00000 0.750000 0
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.00000 0.00000 1.00000 0
+DEAL::0.250000 0.00000 1.00000 0
+DEAL::0.00000 0.250000 1.00000 0
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.250000 0.00000 1.00000 0
+DEAL::0.500000 0.00000 1.00000 0
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.00000 0.250000 1.00000 0
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.00000 0.500000 1.00000 0
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::1.00000 0.00000 0.500000 0
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::1.00000 0.00000 0.750000 0
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.500000 0.00000 1.00000 0
+DEAL::0.750000 0.00000 1.00000 0
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::1.00000 0.00000 0.750000 0
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.750000 0.00000 1.00000 0
+DEAL::1.00000 0.00000 1.00000 0
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::1.00000 0.250000 1.00000 0.250000
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::1.00000 0.250000 1.00000 0.250000
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::1.00000 0.500000 1.00000 0.500000
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.00000 1.00000 0.500000 0
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.00000 1.00000 0.750000 0
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.00000 0.500000 1.00000 0
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.00000 0.750000 1.00000 0
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.00000 1.00000 0.750000 0
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.00000 0.750000 1.00000 0
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.00000 1.00000 1.00000 0
+DEAL::0.250000 1.00000 1.00000 0.250000
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.250000 1.00000 1.00000 0.250000
+DEAL::0.500000 1.00000 1.00000 0.500000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::1.00000 1.00000 0.500000 0.500000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::1.00000 1.00000 0.750000 0.750000
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::1.00000 0.500000 1.00000 0.500000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::1.00000 0.750000 1.00000 0.750000
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::0.500000 1.00000 1.00000 0.500000
+DEAL::0.750000 1.00000 1.00000 0.750000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::1.00000 1.00000 0.750000 0.750000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::1.00000 0.750000 1.00000 0.750000
+DEAL::0.750000 1.00000 1.00000 0.750000
+DEAL::1.00000 1.00000 1.00000 1.00000
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2006 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::project for Vector<double> arguments and hp
+
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <fstream>
+
+
+// define a function that is zero within the square
+// [-sqrt(2)/2,sqrt(2)/2]^d
+// which is the domain produced by GridGenerator::hyper_sphere
+// when using a linear mapping. we will want to see a solution
+// that is not zero when using a higher order mapping and a Q2
+// element
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int = 0) const
+ {
+ // compute the linfty norm of p
+ double m = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ m = std::max(m, std::fabs(p[d]));
+
+ // let the value increase linearly away from the square
+ return std::max (0., m-std::sqrt(2.)/2);
+ }
+};
+
+
+template<int dim>
+void test()
+{
+ HyperBallBoundary<dim> boundary;
+
+ // create 2 triangulations with the
+ // same coarse grid, and refine
+ // them differently
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball (tria);
+ tria.set_boundary (0, boundary);
+
+
+ hp::FECollection<dim> fe;
+ fe.push_back(FE_Q<dim>(1));
+ hp::DoFHandler<dim> dh (tria);
+ dh.distribute_dofs (fe);
+
+ Vector<double> v(dh.n_dofs());
+
+ ConstraintMatrix cm;
+ cm.close ();
+
+ // use the implicit Q1 mapping. this will yield a zero solution
+ {
+ VectorTools::project (dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v);
+ deallog << v.l2_norm() << std::endl;
+ Assert (v.l2_norm() == 0, ExcInternalError());
+ }
+
+ // use an explicit Q1 mapping. this will yield a zero solution
+ {
+ VectorTools::project (hp::MappingCollection<dim>(MappingQ1<dim>()),
+ dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v);
+ deallog << v.l2_norm() << std::endl;
+ Assert (v.l2_norm() == 0, ExcInternalError());
+ }
+
+ // use an explicit Q2 mapping. this will yield a nonzero solution if it's a
+ // straight projection since some of the quadrature points are lying outside
+ // the area where the function is zero
+ {
+ VectorTools::project (hp::MappingCollection<dim>(MappingQ<dim>(2)),
+ dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v);
+ deallog << v.l2_norm() << std::endl;
+ Assert (v.l2_norm() != 0, ExcInternalError());
+ }
+
+ // use an explicit Q2 mapping but enforce zero boundary values. this will
+ // yield a nonzero solution since some of the quadrature points are lying
+ // outside the area where the function is zero, even though the values of
+ // the DoFs at the boundary are in fact zero (they are interpolated only at
+ // points where the function is zero)
+ {
+ VectorTools::project (hp::MappingCollection<dim>(MappingQ<dim>(2)),
+ dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v, true);
+ deallog << v.l2_norm() << std::endl;
+ Assert (v.l2_norm() != 0, ExcInternalError());
+ }
+
+ // use an explicit Q2 mapping and project to the boundary first. this will
+ // yield a nonzero solution since some of the quadrature points are lying
+ // outside the area where the function is zero. furthermore, the values on
+ // the boundary will be nonzero since we project, even though the
+ // *interpolation* of boundary values onto the trace of the Q1 space is zero
+ {
+ VectorTools::project (hp::MappingCollection<dim>(MappingQ<dim>(2)),
+ dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v, false,
+ hp::QCollection<dim-1>(QGauss<dim-1>(2)), true);
+ deallog << v.l2_norm() << std::endl;
+ Assert (v.l2_norm() != 0, ExcInternalError());
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=dh.begin_active();
+ cell != dh.end(); ++cell)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ deallog << cell->vertex(i) << ' ' << v(cell->vertex_dof_index(i,0,cell->active_fe_index()))
+ << std::endl;
+ }
+
+
+ // same as above, but use a projection with a QTrapez formula. this happens
+ // to evaluate the function only at points where it is zero, and
+ // consequently the values at the boundary should be zero
+ {
+ VectorTools::project (hp::MappingCollection<dim>(MappingQ<dim>(2)),
+ dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v, false,
+ hp::QCollection<dim-1>(QTrapez<dim-1>()), true);
+ deallog << v.l2_norm() << std::endl;
+ Assert (v.l2_norm() != 0, ExcInternalError());
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=dh.begin_active();
+ cell != dh.end(); ++cell)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ deallog << cell->vertex(i) << ' ' << v(cell->vertex_dof_index(i,0,cell->active_fe_index()))
+ << std::endl;
+ }
+}
+
+
+int main()
+{
+ std::ofstream logfile ("project_01_curved_boundary/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<2>();
+ test<3>();
+}
+
--- /dev/null
+
+DEAL::0
+DEAL::0
+DEAL::0.227756
+DEAL::0.0359680
+DEAL::0.421875
+DEAL::-0.707107 -0.707107 0.195262
+DEAL::0.707107 -0.707107 0.195262
+DEAL::-0.292893 -0.292893 -0.0797954
+DEAL::0.292893 -0.292893 -0.0797954
+DEAL::-0.707107 -0.707107 0.195262
+DEAL::-0.292893 -0.292893 -0.0797954
+DEAL::-0.707107 0.707107 0.195262
+DEAL::-0.292893 0.292893 -0.0797954
+DEAL::-0.292893 -0.292893 -0.0797954
+DEAL::0.292893 -0.292893 -0.0797954
+DEAL::-0.292893 0.292893 -0.0797954
+DEAL::0.292893 0.292893 -0.0797954
+DEAL::0.707107 -0.707107 0.195262
+DEAL::0.707107 0.707107 0.195262
+DEAL::0.292893 -0.292893 -0.0797954
+DEAL::0.292893 0.292893 -0.0797954
+DEAL::-0.707107 0.707107 0.195262
+DEAL::-0.292893 0.292893 -0.0797954
+DEAL::0.707107 0.707107 0.195262
+DEAL::0.292893 0.292893 -0.0797954
+DEAL::0.0359680
+DEAL::-0.707107 -0.707107 0
+DEAL::0.707107 -0.707107 0
+DEAL::-0.292893 -0.292893 0.0179840
+DEAL::0.292893 -0.292893 0.0179840
+DEAL::-0.707107 -0.707107 0
+DEAL::-0.292893 -0.292893 0.0179840
+DEAL::-0.707107 0.707107 0
+DEAL::-0.292893 0.292893 0.0179840
+DEAL::-0.292893 -0.292893 0.0179840
+DEAL::0.292893 -0.292893 0.0179840
+DEAL::-0.292893 0.292893 0.0179840
+DEAL::0.292893 0.292893 0.0179840
+DEAL::0.707107 -0.707107 0
+DEAL::0.707107 0.707107 0
+DEAL::0.292893 -0.292893 0.0179840
+DEAL::0.292893 0.292893 0.0179840
+DEAL::-0.707107 0.707107 0
+DEAL::-0.292893 0.292893 0.0179840
+DEAL::0.707107 0.707107 0
+DEAL::0.292893 0.292893 0.0179840
+DEAL::0
+DEAL::0
+DEAL::0.166690
+DEAL::0.0709393
+DEAL::0.440214
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 -0.577350 0.115757
+DEAL::0.577350 -0.577350 -0.577350 0.115757
+DEAL::-0.577350 0.577350 -0.577350 0.115757
+DEAL::0.577350 0.577350 -0.577350 0.115757
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.577350 -0.577350 -0.577350 0.115757
+DEAL::0.577350 0.577350 -0.577350 0.115757
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.577350 -0.577350 0.577350 0.115757
+DEAL::0.577350 0.577350 0.577350 0.115757
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 0.577350 0.115757
+DEAL::0.577350 -0.577350 0.577350 0.115757
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.577350 0.577350 0.577350 0.115757
+DEAL::0.577350 0.577350 0.577350 0.115757
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 -0.577350 0.115757
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.577350 0.577350 -0.577350 0.115757
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::-0.577350 -0.577350 0.577350 0.115757
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.577350 0.577350 0.577350 0.115757
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 -0.577350 0.115757
+DEAL::0.577350 -0.577350 -0.577350 0.115757
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.577350 -0.577350 0.577350 0.115757
+DEAL::0.577350 -0.577350 0.577350 0.115757
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.577350 0.577350 -0.577350 0.115757
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.577350 0.577350 -0.577350 0.115757
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::-0.577350 0.577350 0.577350 0.115757
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::0.577350 0.577350 0.577350 0.115757
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::0.0709393
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 -0.577350 0
+DEAL::0.577350 -0.577350 -0.577350 0
+DEAL::-0.577350 0.577350 -0.577350 0
+DEAL::0.577350 0.577350 -0.577350 0
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.577350 -0.577350 -0.577350 0
+DEAL::0.577350 0.577350 -0.577350 0
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.577350 -0.577350 0.577350 0
+DEAL::0.577350 0.577350 0.577350 0
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 0.577350 0
+DEAL::0.577350 -0.577350 0.577350 0
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.577350 0.577350 0.577350 0
+DEAL::0.577350 0.577350 0.577350 0
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 -0.577350 0
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.577350 0.577350 -0.577350 0
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::-0.577350 -0.577350 0.577350 0
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.577350 0.577350 0.577350 0
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 -0.577350 0
+DEAL::0.577350 -0.577350 -0.577350 0
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.577350 -0.577350 0.577350 0
+DEAL::0.577350 -0.577350 0.577350 0
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.577350 0.577350 -0.577350 0
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.577350 0.577350 -0.577350 0
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::-0.577350 0.577350 0.577350 0
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::0.577350 0.577350 0.577350 0
+DEAL::0.211325 0.211325 0.211325 0.0250808
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2006 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::project for BlockVector<double> arguments and hp::
+
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <fstream>
+
+
+// define the multi-linear function x or x*y or x*y*z that we will
+// subsequently project onto the ansatz space
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int = 0) const
+ {
+ double s = 1;
+ for (unsigned int i=0; i<dim; ++i)
+ s *= p[i];
+ return s;
+ }
+};
+
+
+template<int dim>
+void test()
+{
+ // create 2 triangulations with the
+ // same coarse grid, and refine
+ // them differently
+ Triangulation<dim> tria;
+
+ GridGenerator::hyper_cube (tria);
+ tria.refine_global (2);
+
+ hp::FECollection<dim> fe;
+ fe.push_back(FE_Q<dim>(1));
+ hp::DoFHandler<dim> dh (tria);
+ dh.distribute_dofs (fe);
+
+ BlockVector<double> v(2);
+ v.block(0).reinit (dh.n_dofs()/2);
+ v.block(1).reinit (dh.n_dofs()-dh.n_dofs()/2);
+ v.collect_sizes();
+
+ ConstraintMatrix cm;
+ cm.close ();
+ VectorTools::project (dh, cm, hp::QCollection<dim>(QGauss<dim>(3)), F<dim>(),
+ v);
+
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=dh.begin_active();
+ cell != dh.end(); ++cell)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ // check that the error is
+ // somewhat small. it won't
+ // be zero since we project
+ // and do not interpolate
+ Assert (std::fabs (v(cell->vertex_dof_index(i,0,cell->active_fe_index())) -
+ F<dim>().value (cell->vertex(i)))
+ < 1e-4,
+ ExcInternalError());
+ deallog << cell->vertex(i) << ' ' << v(cell->vertex_dof_index(i,0,cell->active_fe_index()))
+ << std::endl;
+ }
+}
+
+
+int main()
+{
+ std::ofstream logfile ("project_02/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1>();
+ test<2>();
+ test<3>();
+}
+
--- /dev/null
+
+DEAL::0.00000 0
+DEAL::0.250000 0.250000
+DEAL::0.250000 0.250000
+DEAL::0.500000 0.500000
+DEAL::0.500000 0.500000
+DEAL::0.750000 0.750000
+DEAL::0.750000 0.750000
+DEAL::1.00000 1.00000
+DEAL::0.00000 0.00000 0
+DEAL::0.250000 0.00000 0
+DEAL::0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.250000 0.00000 0
+DEAL::0.500000 0.00000 0
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.00000 0.500000 0
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.250000 0.250000 0.0625000
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.500000 0.00000 0
+DEAL::0.750000 0.00000 0
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.750000 0.250000 0.187500
+DEAL::0.750000 0.00000 0
+DEAL::1.00000 0.00000 0
+DEAL::0.750000 0.250000 0.187500
+DEAL::1.00000 0.250000 0.250000
+DEAL::0.500000 0.250000 0.125000
+DEAL::0.750000 0.250000 0.187500
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.750000 0.500000 0.375000
+DEAL::0.750000 0.250000 0.187500
+DEAL::1.00000 0.250000 0.250000
+DEAL::0.750000 0.500000 0.375000
+DEAL::1.00000 0.500000 0.500000
+DEAL::0.00000 0.500000 0
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.00000 0.750000 0
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.250000 0.500000 0.125000
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.00000 0.750000 0
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.00000 1.00000 0
+DEAL::0.250000 1.00000 0.250000
+DEAL::0.250000 0.750000 0.187500
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.250000 1.00000 0.250000
+DEAL::0.500000 1.00000 0.500000
+DEAL::0.500000 0.500000 0.250000
+DEAL::0.750000 0.500000 0.375000
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.562500
+DEAL::0.750000 0.500000 0.375000
+DEAL::1.00000 0.500000 0.500000
+DEAL::0.750000 0.750000 0.562500
+DEAL::1.00000 0.750000 0.750000
+DEAL::0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.562500
+DEAL::0.500000 1.00000 0.500000
+DEAL::0.750000 1.00000 0.750000
+DEAL::0.750000 0.750000 0.562500
+DEAL::1.00000 0.750000 0.750000
+DEAL::0.750000 1.00000 0.750000
+DEAL::1.00000 1.00000 1.00000
+DEAL::0.00000 0.00000 0.00000 0
+DEAL::0.250000 0.00000 0.00000 0
+DEAL::0.00000 0.250000 0.00000 0
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.00000 0.00000 0.250000 0
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.250000 0.00000 0.00000 0
+DEAL::0.500000 0.00000 0.00000 0
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.00000 0.250000 0.00000 0
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.00000 0.500000 0.00000 0
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.250000 0.250000 0.00000 0
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.00000 0.00000 0.250000 0
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.00000 0.00000 0.500000 0
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.250000 0.00000 0.250000 0
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.00000 0.250000 0.250000 0
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.250000 0.250000 0.250000 0.0156250
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.500000 0.00000 0.00000 0
+DEAL::0.750000 0.00000 0.00000 0
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.750000 0.00000 0.00000 0
+DEAL::1.00000 0.00000 0.00000 0
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::1.00000 0.250000 0.00000 0
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::1.00000 0.00000 0.250000 0
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.500000 0.250000 0.00000 0
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.750000 0.250000 0.00000 0
+DEAL::1.00000 0.250000 0.00000 0
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::1.00000 0.500000 0.00000 0
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.500000 0.00000 0.250000 0
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.750000 0.00000 0.250000 0
+DEAL::1.00000 0.00000 0.250000 0
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::1.00000 0.00000 0.500000 0
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.500000 0.250000 0.250000 0.0312500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.750000 0.250000 0.250000 0.0468750
+DEAL::1.00000 0.250000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.00000 0.500000 0.00000 0
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.00000 0.750000 0.00000 0
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.250000 0.500000 0.00000 0
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.00000 0.750000 0.00000 0
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.00000 1.00000 0.00000 0
+DEAL::0.250000 1.00000 0.00000 0
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.00000 1.00000 0.250000 0
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.250000 0.750000 0.00000 0
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.250000 1.00000 0.00000 0
+DEAL::0.500000 1.00000 0.00000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.00000 0.500000 0.250000 0
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.250000 0.500000 0.250000 0.0312500
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.00000 0.750000 0.250000 0
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.00000 1.00000 0.250000 0
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.00000 1.00000 0.500000 0
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.250000 0.750000 0.250000 0.0468750
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.250000 1.00000 0.250000 0.0625000
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.500000 0.500000 0.00000 0
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.750000 0.500000 0.00000 0
+DEAL::1.00000 0.500000 0.00000 0
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::1.00000 0.750000 0.00000 0
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.500000 0.750000 0.00000 0
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::0.500000 1.00000 0.00000 0
+DEAL::0.750000 1.00000 0.00000 0
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::0.750000 0.750000 0.00000 0
+DEAL::1.00000 0.750000 0.00000 0
+DEAL::0.750000 1.00000 0.00000 0
+DEAL::1.00000 1.00000 0.00000 0
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::1.00000 1.00000 0.250000 0.250000
+DEAL::0.500000 0.500000 0.250000 0.0625000
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.750000 0.500000 0.250000 0.0937500
+DEAL::1.00000 0.500000 0.250000 0.125000
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.500000 0.750000 0.250000 0.0937500
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::0.500000 1.00000 0.250000 0.125000
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::0.750000 0.750000 0.250000 0.140625
+DEAL::1.00000 0.750000 0.250000 0.187500
+DEAL::0.750000 1.00000 0.250000 0.187500
+DEAL::1.00000 1.00000 0.250000 0.250000
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::1.00000 1.00000 0.500000 0.500000
+DEAL::0.00000 0.00000 0.500000 0
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.00000 0.00000 0.750000 0
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.250000 0.00000 0.500000 0
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.00000 0.250000 0.500000 0
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.250000 0.250000 0.500000 0.0312500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.00000 0.00000 0.750000 0
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.00000 0.00000 1.00000 0
+DEAL::0.250000 0.00000 1.00000 0
+DEAL::0.00000 0.250000 1.00000 0
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.250000 0.00000 0.750000 0
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.250000 0.00000 1.00000 0
+DEAL::0.500000 0.00000 1.00000 0
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.00000 0.250000 0.750000 0
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.00000 0.250000 1.00000 0
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.00000 0.500000 1.00000 0
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.250000 0.250000 0.750000 0.0468750
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.250000 0.250000 1.00000 0.0625000
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.500000 0.00000 0.500000 0
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.750000 0.00000 0.500000 0
+DEAL::1.00000 0.00000 0.500000 0
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::1.00000 0.00000 0.750000 0
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.500000 0.250000 0.500000 0.0625000
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.750000 0.250000 0.500000 0.0937500
+DEAL::1.00000 0.250000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.500000 0.00000 0.750000 0
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.500000 0.00000 1.00000 0
+DEAL::0.750000 0.00000 1.00000 0
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::0.750000 0.00000 0.750000 0
+DEAL::1.00000 0.00000 0.750000 0
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.750000 0.00000 1.00000 0
+DEAL::1.00000 0.00000 1.00000 0
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::1.00000 0.250000 1.00000 0.250000
+DEAL::0.500000 0.250000 0.750000 0.0937500
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.500000 0.250000 1.00000 0.125000
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::0.750000 0.250000 0.750000 0.140625
+DEAL::1.00000 0.250000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.750000 0.250000 1.00000 0.187500
+DEAL::1.00000 0.250000 1.00000 0.250000
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::1.00000 0.500000 1.00000 0.500000
+DEAL::0.00000 0.500000 0.500000 0
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.250000 0.500000 0.500000 0.0625000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.00000 0.750000 0.500000 0
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.00000 1.00000 0.500000 0
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.00000 1.00000 0.750000 0
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.250000 0.750000 0.500000 0.0937500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.250000 1.00000 0.500000 0.125000
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.00000 0.500000 0.750000 0
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.00000 0.500000 1.00000 0
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.00000 0.750000 1.00000 0
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.250000 0.500000 0.750000 0.0937500
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.250000 0.500000 1.00000 0.125000
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.00000 0.750000 0.750000 0
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.00000 1.00000 0.750000 0
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.00000 0.750000 1.00000 0
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.00000 1.00000 1.00000 0
+DEAL::0.250000 1.00000 1.00000 0.250000
+DEAL::0.250000 0.750000 0.750000 0.140625
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.250000 1.00000 0.750000 0.187500
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.250000 0.750000 1.00000 0.187500
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.250000 1.00000 1.00000 0.250000
+DEAL::0.500000 1.00000 1.00000 0.500000
+DEAL::0.500000 0.500000 0.500000 0.125000
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.750000 0.500000 0.500000 0.187500
+DEAL::1.00000 0.500000 0.500000 0.250000
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.500000 0.750000 0.500000 0.187500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::0.500000 1.00000 0.500000 0.250000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::0.750000 0.750000 0.500000 0.281250
+DEAL::1.00000 0.750000 0.500000 0.375000
+DEAL::0.750000 1.00000 0.500000 0.375000
+DEAL::1.00000 1.00000 0.500000 0.500000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::1.00000 1.00000 0.750000 0.750000
+DEAL::0.500000 0.500000 0.750000 0.187500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.500000 0.500000 1.00000 0.250000
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::0.750000 0.500000 0.750000 0.281250
+DEAL::1.00000 0.500000 0.750000 0.375000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.750000 0.500000 1.00000 0.375000
+DEAL::1.00000 0.500000 1.00000 0.500000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::1.00000 0.750000 1.00000 0.750000
+DEAL::0.500000 0.750000 0.750000 0.281250
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::0.500000 1.00000 0.750000 0.375000
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::0.500000 0.750000 1.00000 0.375000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::0.500000 1.00000 1.00000 0.500000
+DEAL::0.750000 1.00000 1.00000 0.750000
+DEAL::0.750000 0.750000 0.750000 0.421875
+DEAL::1.00000 0.750000 0.750000 0.562500
+DEAL::0.750000 1.00000 0.750000 0.562500
+DEAL::1.00000 1.00000 0.750000 0.750000
+DEAL::0.750000 0.750000 1.00000 0.562500
+DEAL::1.00000 0.750000 1.00000 0.750000
+DEAL::0.750000 1.00000 1.00000 0.750000
+DEAL::1.00000 1.00000 1.00000 1.00000