but the scalar type of the vector was always <code>double</code>. This
has been changed so that the two always need to match.
<br>
- (David Davydov, Wolfgang Bangerth, 2015/06/17)
+ (Denis Davydov, Wolfgang Bangerth, 2015/06/17)
</li>
<li>Changed: Functions such as FEValuesBase::get_function_values() that
complex-valued data types. This includes compiling PETSc with
complex scalars.
<br>
- (David Davydov, 2015/05/08)
+ (Denis Davydov, 2015/05/08)
</li>
<li>Removed: The generic, templated vmult, Tvmult, etc. -interfaces of
<ol>
+ <li> New: Implement a modified version of the Kelly error estimator, which
+ effectively provides the boundary residual term for the hp-FEM error estimators.
+ <br>
+ (Denis Davydov, 2015/06/17)
+ </li>
+
<li> New: DerivativeForm() now takes an additional optional
template argument specifying the type, similarly to Tensor() classes.
<br>
/**
* Implementation of the error indicator by Kelly, De S. R. Gago, Zienkiewicz and
- * Babuska. This error indicator tries to approximate the error per cell by
+ * Babuska and its modification for the hp-FEM.
+ * This error indicator tries to approximate the error per cell by
* integration of the jump of the gradient of the solution along the faces of
* each cell. It can be understood as a gradient recovery estimator; see the
* survey of Ainsworth and Oden, "A Posteriori Error Estimation in Finite Element
* Analysis" (Wiley, 2000) for a complete discussion.
*
- * @note In spite of the name, this is not truly an a posteriori error
+ * In the original Kelly error estimator, the contribution of each face to the
+ * cell error is scaled with the cell diagonal. In the modified version,
+ * however, we employ a scaling factor which depends on the face diagonal and
+ * polynomial degrees of the adjacent elements. The choice between the two is
+ * done by means of the enumerator, defined within the class.
+ *
+ * @note In spite of the name, Kelly estimator is not truly an a posteriori error
* estimator, even if applied to the Poisson problem only. It gives good hints
* for mesh refinement, but the estimate is not to be trusted. For higher
* order trial spaces the integrals computed here tend to zero faster than the
- * error itself, thus ruling out the values as error estimators.
+ * error itself, thus ruling out the values as error estimators. However, the
+ * modified version discussed below can be utilised to obtain the reliable
+ * error estimator by adding the residual (volume) part.
*
* The error estimator really only estimates the error for the generalized
* Poisson equation $-\nabla\cdot a(x) \nabla u = f$ with either Dirichlet
* <h3>Implementation</h3>
*
* In principle, the implementation of the error estimation is simple: let \f[
- * \eta_K^2 = \frac h{24} \int_{\partial K} \left[a \frac{\partial
+ * \eta_K^2 = \sum_{F\in\partial K} c_F \int_{\partial K_F} \left[a \frac{\partial
* u_h}{\partial n}\right]^2 do \f] be the error estimator for cell $K$.
* $[\cdot]$ denotes the jump of the argument at the face. In the paper of
- * Ainsworth, $h$ is divided by $24$, but this factor is a bit esoteric,
+ * Ainsworth $ c_F=\frac h{24} $, but this factor is a bit esoteric,
* stemming from interpolation estimates and stability constants which may
* hold for the Poisson problem, but may not hold for more general situations.
- * In the implementation, this factor is considered, but may lead to wrong
- * results. You may scale the vector appropriately afterwards.
+ * Alternatively, we consider the case when $ c_F=\frac {h_F}_{2p_F} $,
+ * where $ h_F $ is face diagonal and $ p_F=max(p^+,p^-) $ is the
+ * maximum polynomial degree of adjacent elements. The choice between the two is
+ * done by means of the enumerator, provided as the last argument in all functions.
*
* To perform the integration, use is made of the FEFaceValues and
* FESubfaceValues classes. The integration is performed by looping over all
* cells and integrating over faces that are not yet treated. This way we
* avoid integration on faces twice, once for each time we visit one of the
* adjacent cells. In a second loop over all cells, we sum up the
- * contributions of the faces (which are the integrated square of the jumps)
- * of each cell and take the square root.
+ * contributions of the faces (which are the integrated square of the jumps times
+ * some factor) of each cell and take the square root.
*
* The integration is done using a quadrature formula on the face. For linear
* trial functions (FEQ1), the QGauss2 or even the QMidpoint rule will
*
* We store the contribution of each face in a @p map, as provided by the C++
* standard library, with the iterator pointing to that face being the key
- * into the map. In fact, we do not store the indicator per face, but only the
- * integral listed above. When looping the second time over all cells, we have
- * to sum up the contributions of the faces, multiply them with $\frac h{24}$
- * and take the square root. By doing the multiplication with $h$ in the
- * second loop, we avoid problems to decide with which $h$ to multiply, that
- * of the cell on the one or that of the cell on the other side of the face.
+ * into the map. When looping the second time over all cells,
+ * we have to sum up the contributions of the faces and take the square root.
+ * For the Kelly estimator, the multiplication with $\frac h{24}$ is done
+ * in the second loop. By doing so we avoid problems to decide with which $h$
+ * to multiply, that of the cell on the one or that of the cell on the other
+ * side of the face. Whereas for the hp-estimator the @p map stores integrals
+ * multiplied by $\frac {h_f}_{2p_f}$, which are then summed in the second loop.
*
- * $h$ is taken to be the greatest length of the diagonals of the cell. For
- * more or less uniform cells without deformed angles, this coincides with the
- * diameter of the cell.
+ * $h$ ($h_f$) is taken to be the greatest length of the diagonals of the cell (face).
+ * For more or less uniform cells (faces) without deformed angles, this coincides
+ * with the diameter of the cell (face).
*
*
* <h3>Vector-valued functions</h3>
* different faces to the cells easier.
*
* <li> The face belongs to a Neumann boundary. In this case, the
- * contribution of the face $F\in\partial K$ looks like \f[ \int_F
+ * contribution of the face $F\in\partial K$ looks like \f[ n_F\int_F
* \left|g-a\frac{\partial u_h}{\partial n}\right|^2 ds \f] where $g$ is the
- * Neumann boundary function. If the finite element is vector-valued, then
+ * Neumann boundary function, $n_F=\frac h_{24}$ and $n_F=\frac {h_F}{p}$ for
+ * the Kelly and hp-estimator, respectively.
+ * If the finite element is vector-valued, then
* obviously the function denoting the Neumann boundary conditions needs to be
* vector-valued as well.
*
* that accepts several in- and output vectors at the same time.
*
* @ingroup numerics
- * @author Wolfgang Bangerth, 1998, 1999, 2000, 2004, 2006; parallelization by
- * Thomas Richter, 2000
+ * @author Wolfgang Bangerth, 1998, 1999, 2000, 2004, 2006, Denis Davydov, 2015;
+ * parallelization by Thomas Richter, 2000
*/
template <int dim, int spacedim=dim>
class KellyErrorEstimator
{
public:
+ /**
+ * The enum type given to the class functions to decide on the scaling factors
+ * of the facial integrals.
+ */
+ enum Strategy
+ {
+ //! Kelly error estimator with the factor $\frac h_{24}$.
+ cell_diameter_over_24 = 0,
+ //! the boundary residual estimator with the factor $\frac {h_F}_{2 max(p^+,p^-)$.
+ face_diameter_over_twice_max_degree
+ };
+
/**
* Implementation of the error estimator described above. You may give a
* coefficient, but there is a default value which denotes the constant
* the number of threads determined automatically. The parameter is retained
* for compatibility with old versions of the library.
*
+ * The @p strategy parameter is used to choose the scaling factor for
+ * the integral over cell's faces.
+ *
* @note If the DoFHandler object given as an argument to this function
* builds on a parallel::distributed::Triangulation, this function skips
* computations on all cells that are not locally owned. In that case, the
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
* Calls the @p estimate function, see above, with
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
* Same function as above, but accepts more than one solution vector and
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
* Calls the @p estimate function, see above, with
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
const Function<spacedim> *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
- const types::material_id material_id = numbers::invalid_material_id);
-
+ const types::material_id material_id = numbers::invalid_material_id,
+ const KellyErrorEstimator<dim,spacedim>::Strategy strategy = cell_diameter_over_24);
/**
* Exception
return face_integral;
}
+ /**
+ * A factor to scale the integral for the face at the boundary.
+ * Used for Neumann BC.
+ */
+ template <class DH>
+ double boundary_face_factor(const typename DH::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_cell,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
+ {
+ switch (strategy)
+ {
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::cell_diameter_over_24 :
+ return 1.0;
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::face_diameter_over_twice_max_degree :
+ const double cell_degree = fe_face_values_cell.get_fe_collection()[cell->active_fe_index()].degree;
+ return cell->face(face_no)->diameter() / cell_degree;
+ }
+ }
+
+
+ /**
+ * A factor to scale the integral for the regular face.
+ */
+ template <class DH>
+ double regular_face_factor(const typename DH::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_cell,
+ const dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_neighbor,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
+ {
+ switch (strategy)
+ {
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::cell_diameter_over_24 :
+ return 1.0;
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::face_diameter_over_twice_max_degree :
+ const double cell_degree = fe_face_values_cell.get_fe_collection()[cell->active_fe_index()].degree;
+ const double neighbor_degree = fe_face_values_neighbor.get_fe_collection()[cell->neighbor(face_no)->active_fe_index()].degree;
+ return cell->face(face_no)->diameter() / std::max(cell_degree,neighbor_degree) / 2.0;
+ }
+ }
+
+ /**
+ * A factor to scale the integral for the irregular face.
+ */
+ template <class DH>
+ double irregular_face_factor(const typename DH::active_cell_iterator &cell,
+ const typename DH::active_cell_iterator &neighbor_child,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values,
+ dealii::hp::FESubfaceValues<DH::dimension, DH::space_dimension> &fe_subface_values,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
+ {
+ switch (strategy)
+ {
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::cell_diameter_over_24 :
+ return 1.0;
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::face_diameter_over_twice_max_degree :
+ const double cell_degree = fe_face_values.get_fe_collection()[cell->active_fe_index()].degree;
+ const double neighbor_child_degree = fe_subface_values.get_fe_collection()[neighbor_child->active_fe_index()].degree;
+ return cell->face(face_no)->child(subface_no)->diameter()/std::max(neighbor_child_degree,cell_degree)/2.0;
+ }
+ }
+
+ /**
+ * A factor used when summing up all the contribution
+ * from different faces of each cell.
+ */
+ template <class DH>
+ double cell_factor(const typename DH::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const DH &dof_handler,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
+ {
+ switch (strategy)
+ {
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::cell_diameter_over_24 :
+ return cell->diameter()/24;
+ case KellyErrorEstimator<DH::dimension,DH::space_dimension>::face_diameter_over_twice_max_degree :
+ return 1.0;
+ }
+ }
+
/**
const typename DH::active_cell_iterator &cell,
const unsigned int face_no,
dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_cell,
- dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_neighbor)
+ dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_neighbor,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
{
const unsigned int dim = DH::dimension;
(void)dim;
fe_face_values_cell.get_present_fe_values()
.get_function_gradients (*solutions[n], parallel_data.psi[n]);
+ double factor;
// now compute over the other side of the face
if (face->at_boundary() == false)
// internal face; integrate jump of gradient across this face
fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor,
cell->active_fe_index());
+ factor = regular_face_factor<DH>(cell,face_no,
+ fe_face_values_cell,fe_face_values_neighbor,
+ strategy);
+
// get gradients on neighbor cell
for (unsigned int n=0; n<n_solution_vectors; ++n)
{
parallel_data.neighbor_psi[n]);
}
}
+ else
+ {
+ factor = boundary_face_factor<DH>(cell,face_no,
+ fe_face_values_cell,
+ strategy);
+ }
// now go to the generic function that does all the other things
local_face_integrals[face] =
integrate_over_face (parallel_data, face, fe_face_values_cell);
+ for (unsigned int i = 0; i < local_face_integrals[face].size(); i++)
+ local_face_integrals[face][i] *= factor;
}
const typename DH::active_cell_iterator &cell,
const unsigned int face_no,
dealii::hp::FEFaceValues<DH::dimension,DH::space_dimension> &fe_face_values,
- dealii::hp::FESubfaceValues<DH::dimension, DH::space_dimension> &fe_subface_values)
+ dealii::hp::FESubfaceValues<DH::dimension, DH::space_dimension> &fe_subface_values,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
{
const unsigned int dim = DH::dimension;
(void)dim;
fe_face_values.reinit (neighbor_child, neighbor_neighbor,
cell->active_fe_index());
+ const double factor = irregular_face_factor<DH>(cell,
+ neighbor_child,
+ face_no,
+ subface_no,
+ fe_face_values,
+ fe_subface_values,
+ strategy);
+
// store the gradient of the solution in psi
for (unsigned int n=0; n<n_solution_vectors; ++n)
fe_subface_values.get_present_fe_values()
// call generic evaluate function
local_face_integrals[neighbor_child->face(neighbor_neighbor)] =
integrate_over_face (parallel_data, face, fe_face_values);
+ for (unsigned int i = 0; i < local_face_integrals[neighbor_child->face(neighbor_neighbor)].size(); i++)
+ local_face_integrals[neighbor_child->face(neighbor_neighbor)][i] *= factor;
}
// finally loop over all subfaces to collect the contributions of the
estimate_one_cell (const typename DH::active_cell_iterator &cell,
ParallelData<DH, typename InputVector::value_type> ¶llel_data,
std::map<typename DH::face_iterator,std::vector<double> > &local_face_integrals,
- const std::vector<const InputVector *> &solutions)
+ const std::vector<const InputVector *> &solutions,
+ const typename KellyErrorEstimator<DH::dimension,DH::space_dimension>::Strategy strategy)
{
const unsigned int dim = DH::dimension;
const unsigned int n_solution_vectors = solutions.size();
local_face_integrals,
cell, face_no,
parallel_data.fe_face_values_cell,
- parallel_data.fe_face_values_neighbor);
+ parallel_data.fe_face_values_neighbor,
+ strategy);
else
// otherwise we need to do some special computations which do not
local_face_integrals,
cell, face_no,
parallel_data.fe_face_values_cell,
- parallel_data.fe_subface_values);
+ parallel_data.fe_subface_values,
+ strategy);
}
}
}
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
// just pass on to the other function
const std::vector<const InputVector *> solutions (1, &solution);
std::vector<Vector<float>*> errors (1, &error);
estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors,
- component_mask, coefficients, n_threads, subdomain_id, material_id);
+ component_mask, coefficients, n_threads, subdomain_id, material_id, strategy);
}
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
estimate(StaticMappingQ1<dim,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution,
error, component_mask, coefficients, n_threads,
- subdomain_id, material_id);
+ subdomain_id, material_id, strategy);
}
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
// just pass on to the other function
const std::vector<const InputVector *> solutions (1, &solution);
std::vector<Vector<float>*> errors (1, &error);
estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors,
- component_mask, coefficients, n_threads, subdomain_id, material_id);
+ component_mask, coefficients, n_threads, subdomain_id, material_id, strategy);
}
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
estimate(StaticMappingQ1<dim, spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution,
error, component_mask, coefficients, n_threads,
- subdomain_id, material_id);
+ subdomain_id, material_id, strategy);
}
const Function<spacedim> *coefficients,
const unsigned int ,
const types::subdomain_id subdomain_id_,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
#ifdef DEAL_II_WITH_P4EST
if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
WorkStream::run (dof_handler.begin_active(),
static_cast<typename DH::active_cell_iterator>(dof_handler.end()),
std_cxx11::bind (&internal::estimate_one_cell<InputVector,DH>,
- std_cxx11::_1, std_cxx11::_2, std_cxx11::_3, std_cxx11::ref(solutions)),
+ std_cxx11::_1, std_cxx11::_2, std_cxx11::_3, std_cxx11::ref(solutions),strategy),
std_cxx11::bind (&internal::copy_local_to_global<DH>,
std_cxx11::_1, std_cxx11::ref(face_integrals)),
parallel_data,
Assert(face_integrals.find(cell->face(face_no))
!= face_integrals.end(),
ExcInternalError());
- const double factor = cell->diameter() / 24;
+ const double factor = internal::cell_factor<DH>(cell,
+ face_no,
+ dof_handler,
+ strategy);
for (unsigned int n=0; n<n_solution_vectors; ++n)
{
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
// forward to the function with the QCollection
estimate (mapping, dof_handler,
hp::QCollection<dim-1>(quadrature),
neumann_bc, solutions,
errors, component_mask, coefficients,
- n_threads, subdomain_id, material_id);
+ n_threads, subdomain_id, material_id, strategy);
}
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
estimate(StaticMappingQ1<dim, spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions,
errors, component_mask, coefficients, n_threads,
- subdomain_id, material_id);
+ subdomain_id, material_id, strategy);
}
const Function<spacedim> *coefficients,
const unsigned int n_threads,
const types::subdomain_id subdomain_id,
- const types::material_id material_id)
+ const types::material_id material_id,
+ const typename KellyErrorEstimator<dim,spacedim>::Strategy strategy)
{
estimate(StaticMappingQ1<dim, spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions,
errors, component_mask, coefficients, n_threads,
- subdomain_id, material_id);
+ subdomain_id, material_id, strategy);
}
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
template
void
const Function<deal_II_space_dimension> *,
const unsigned int ,
const unsigned int ,
- const types::material_id);
+ const types::material_id,
+ const KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::Strategy);
#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2000 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+/* Author: Denis Davydov, University of Erlangen-Nuremberg, 2015 */
+// make sure that the modified Kelly error estimator with
+// strategy = face_diameter_over_twice_max_degree
+// returns correct results.
+// To that end consider 3 cases at domain [0,2]^dim and one at domain [0,1]^dim:
+//
+// 1)
+// Neuman BC g(x) = c; // constant
+// single element with p =3; u=0;
+// \eta = h * A * c^2 / p;
+//
+//
+// 2)
+//
+// ---------------
+// | | |
+// | | |
+// | p1 | p2 |
+// | | |
+// | | |
+// ---------------
+//
+// u (left) = 0;
+// u (right) = kx;
+//
+// \eta = h * A * k^2 / 2 max(p1,p2)
+//
+//
+// 3)
+// ----------------------
+// | | | |
+// | p3 | p3 f_1 |
+// |------------| p1 |
+// | p2 | p2 f_2 |
+// | | | |
+// ----------------------
+// solution is the same as above;
+//
+// \eta_1 = h * A * k^2 / 2 max(p3,p1)
+// \eta_2 = h * A * k^2 / 2 max(p2,p1)
+// \eta_3 = \eta_1 + \eta_2
+//
+// 4) just arbitrary mesh
+// -------------------
+// | | |
+// | 1 | 1 |
+// | | |
+// |------------------
+// | 3 | 2 | |
+// |--------| 1 |
+// | 3 | 2 | |
+// ------------------
+//
+// and evaluate error of the interpolated to it function.
+
+
+
+#include "../tests.h"
+
+// dealii
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+
+// c++
+#include <fstream>
+#include <iostream>
+
+using namespace dealii;
+
+template<int dim>
+class MyFunction : public dealii::Function<dim>
+{
+public:
+ MyFunction(const double k);
+
+ virtual double value(const dealii::Point<dim> & point,
+ const unsigned int component = 0 ) const;
+
+ double get_k() const;
+
+private:
+ const double k;
+};
+
+template<int dim>
+MyFunction<dim>::MyFunction(const double k)
+:
+dealii::Function<dim>(1),
+k(k)
+{
+
+}
+
+template<int dim>
+double MyFunction<dim>::value(const dealii::Point<dim> & point,
+ const unsigned int ) const
+{
+ const double x = point[0]-1.0;
+
+ if (x < 0)
+ return 0.0;
+ else
+ return k * x;
+}
+
+template<int dim>
+double MyFunction<dim>::get_k() const
+{
+ return k;
+}
+
+// neuman bc
+template<int dim>
+class NeumanBC : public dealii::Function<dim>
+{
+public:
+ NeumanBC(const double c);
+
+ virtual double value(const dealii::Point<dim> & point,
+ const unsigned int component = 0 ) const;
+
+ double get_c() const;
+
+private:
+ const double c;
+};
+
+template<int dim>
+NeumanBC<dim>::NeumanBC(const double c)
+:
+dealii::Function<dim>(1),
+c(c)
+{
+}
+
+template<int dim>
+double NeumanBC<dim>::value(const dealii::Point<dim> & point,
+ const unsigned int ) const
+{
+ return c;
+}
+
+template<int dim>
+double NeumanBC<dim>::get_c() const
+{
+ return c;
+}
+
+// helper function to get diagonal and
+// area of the squared element with lenght h
+template<int dim>
+void get_h_area(double &h, double &a, const double L);
+
+template<>
+void get_h_area<2>(double &h, double &a, const double L)
+{
+ h = L;
+ a = L;
+}
+
+template<>
+void get_h_area<3>(double &h, double &a, const double L)
+{
+ h = std::sqrt(2.0)*L;
+ a = L*L;
+}
+
+// helper function to get diagonal and area of the
+// h-refined face.
+template<int dim>
+void get_h_area_sub(double &h, double &a, const double L);
+
+template<>
+void get_h_area_sub<2>(double &h, double &a, const double L)
+{
+ h = L/2;
+ a = L/2;
+}
+
+template<>
+void get_h_area_sub<3>(double &h, double &a, const double L)
+{
+ h = std::sqrt(2.0)*L/2;
+ a = L*L/4.0;
+}
+
+// output for inspection
+template<int dim>
+void output(const std::string name,
+ const Triangulation<dim> &triangulation,
+ const hp::DoFHandler<dim> &dof_handler,
+ const Vector<double> &values,
+ const Vector<float> &error)
+{
+ dealii::Vector<double> fe_degrees(triangulation.n_active_cells());
+ {
+ typename dealii::hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (unsigned int index=0; cell!=endc; ++cell, ++index)
+ fe_degrees(index) = dof_handler.get_fe()[cell->active_fe_index()].degree;
+ }
+
+ dealii::DataOut<dim,dealii::hp::DoFHandler<dim> > data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector(values,
+ std::string("function_interpolation"));
+ data_out.add_data_vector(fe_degrees,
+ std::string("fe_degree"));
+ data_out.add_data_vector(error,
+ std::string("error"));
+ data_out.build_patches ();
+
+ std::ofstream output (name);
+ data_out.write_vtu (output);
+}
+
+// case 1)
+template<int dim>
+void test_neumann(const NeumanBC<dim> &func)
+{
+ deallog << "NeumanBC case:"<<std::endl;
+ deallog << "--------------"<<std::endl;
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_formula;
+ hp::QCollection<dim-1> face_quadrature_formula;
+ ConstraintMatrix constraints;
+
+ const unsigned int p = 3;
+
+ fe_collection.push_back(dealii::FE_Q<dim>(p));
+ quadrature_formula.push_back(dealii::QGauss<dim>(p+5));
+ face_quadrature_formula.push_back(dealii::QGauss<dim-1>(p+5));
+
+ const double L = 2.0;
+
+ // set-up domain
+ {
+ GridGenerator::hyper_cube (triangulation,.0,L,/*colorize*/true);
+ }
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ // constraints
+ constraints.clear();
+ dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ constraints.close ();
+
+ // interpolate some function
+ dealii::Vector<double> values(dof_handler.n_dofs());
+ values = 0.0;
+
+ dealii::deallog << "dof values:"<<std::endl;
+ for (unsigned int i = 0; i < values.size(); i++)
+ dealii::deallog << " " << values[i];
+ dealii::deallog << std::endl;
+
+ // call Kelly
+ typename dealii::FunctionMap<dim>::type function_map;
+ function_map[0] = &func;
+
+ dealii::Vector<float> error(dof_handler.n_dofs());
+ dealii::KellyErrorEstimator<dim>::estimate(dof_handler,
+ face_quadrature_formula,
+ function_map,
+ values,
+ error,
+ dealii::ComponentMask(),
+ 0,
+ dealii::numbers::invalid_unsigned_int,
+ dealii::numbers::invalid_subdomain_id,
+ dealii::numbers::invalid_material_id,
+ dealii::KellyErrorEstimator<dim>::face_diameter_over_twice_max_degree);
+
+ dealii::deallog <<"error:"<<std::endl;
+ for (unsigned int i = 0; i < error.size(); i++)
+ dealii::deallog << " " << error[i];
+ dealii::deallog << std::endl;
+
+ //output("neuman.vtu",
+ // triangulation,
+ // dof_handler,
+ // values,
+ // error);
+
+ double h,A;
+ get_h_area<dim>(h,A,L);
+ const double expected_value_squared = h*A*std::pow(func.get_c(),2)/p;
+ dealii::deallog << "expected:"<< std::endl <<" "<< std::sqrt(expected_value_squared) << std::endl;
+
+ AssertThrow (std::fabs(std::sqrt(expected_value_squared) - error[0] ) < 1e-5, dealii::ExcInternalError());
+
+ dof_handler.clear();
+}
+
+// case 2)
+template<int dim>
+void test_regular(const MyFunction<dim> &func)
+{
+ deallog << std::endl;
+ deallog << "Regular face:"<<std::endl;
+ deallog << "-------------"<<std::endl;
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_formula;
+ hp::QCollection<dim-1> face_quadrature_formula;
+ ConstraintMatrix constraints;
+
+ const unsigned int p1 = 1;
+ const unsigned int p2 = 2;
+ std::vector<unsigned int> p_degree;
+ p_degree.push_back(p1);
+ p_degree.push_back(p2);
+
+ for (unsigned int i=0;i<p_degree.size();i++)
+ {
+ const unsigned int &p = p_degree[i];
+ fe_collection.push_back(dealii::FE_Q<dim>(p));
+ quadrature_formula.push_back(dealii::QGauss<dim>(p+5));
+ face_quadrature_formula.push_back(dealii::QGauss<dim-1>(p+5));
+ }
+
+ const double L = 2.0;
+
+ // set-up domain
+ {
+ std::vector<unsigned int> repetitions(dim,1);
+ repetitions[0] = 2;
+ dealii::Point<dim> p1;
+ dealii::Point<dim> p2;
+ for (unsigned int d = 0; d < dim; d++)
+ {
+ p1[d] = 0.0;
+ p2[d] = L;
+ }
+ GridGenerator::subdivided_hyper_rectangle (triangulation,
+ repetitions,
+ p1,
+ p2,
+ /*colorize*/ false);
+
+ typename dealii::hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (;cell != endc; cell++)
+ if (cell->center()[0] > 1.0)
+ {
+ cell->set_active_fe_index(1);
+ break;
+ }
+
+ }
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ // constraints
+ constraints.clear();
+ dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ constraints.close ();
+
+ // interpolate some function
+ dealii::Vector<double> values(dof_handler.n_dofs());
+ dealii::VectorTools::interpolate (dof_handler,
+ func,
+ values);
+
+ dealii::deallog << "dof values:"<<std::endl;
+ for (unsigned int i = 0; i < values.size(); i++)
+ dealii::deallog << " " << values[i];
+ dealii::deallog << std::endl;
+
+ // call Kelly
+ dealii::Vector<float> error(dof_handler.n_dofs());
+ dealii::KellyErrorEstimator<dim>::estimate(dof_handler,
+ face_quadrature_formula,
+ typename dealii::FunctionMap<dim>::type(),
+ values,
+ error,
+ dealii::ComponentMask(),
+ 0,
+ dealii::numbers::invalid_unsigned_int,
+ dealii::numbers::invalid_subdomain_id,
+ dealii::numbers::invalid_material_id,
+ dealii::KellyErrorEstimator<dim>::face_diameter_over_twice_max_degree);
+
+ dealii::deallog <<"error:"<<std::endl;
+ for (unsigned int i = 0; i < error.size(); i++)
+ dealii::deallog << " " << error[i];
+ dealii::deallog << std::endl;
+
+ //output("regular.vtu",
+ // triangulation,
+ // dof_handler,
+ // values,
+ // error);
+
+ double h,A;
+ get_h_area<dim>(h,A,L);
+ const double expected_value_squared = h*A*std::pow(func.get_k(),2)/2.0/std::max(p1,p2);
+ dealii::deallog << "expected:"<< std::endl <<" "<< std::sqrt(expected_value_squared) << std::endl;
+ for (unsigned int i = 0; i < error.size(); i++)
+ AssertThrow (std::fabs(std::sqrt(expected_value_squared) - error[i] ) < 1e-6, dealii::ExcInternalError());
+
+ dof_handler.clear();
+}
+
+// case 3)
+template<int dim>
+void test_irregular(const MyFunction<dim> &func)
+{
+ deallog << std::endl;
+ deallog << "Irregular face:"<<std::endl;
+ deallog << "---------------"<<std::endl;
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_formula;
+ hp::QCollection<dim-1> face_quadrature_formula;
+ ConstraintMatrix constraints;
+
+ const unsigned int p1 = 1;
+ const unsigned int p2 = 2;
+ const unsigned int p3 = 3;
+ std::vector<unsigned int> p_degree;
+ p_degree.push_back(p1);
+ p_degree.push_back(p2);
+ p_degree.push_back(p3);
+
+ for (unsigned int i=0;i<p_degree.size();i++)
+ {
+ const unsigned int &p = p_degree[i];
+ fe_collection.push_back(dealii::FE_Q<dim>(p));
+ quadrature_formula.push_back(dealii::QGauss<dim>(p+5));
+ face_quadrature_formula.push_back(dealii::QGauss<dim-1>(p+5));
+ }
+
+ const double L = 2.0;
+
+ // set-up domain
+ {
+ std::vector<unsigned int> repetitions(dim,1);
+ repetitions[0] = 2;
+ dealii::Point<dim> p1;
+ dealii::Point<dim> p2;
+ for (unsigned int d = 0; d < dim; d++)
+ {
+ p1[d] = 0.0;
+ p2[d] = L;
+ }
+ GridGenerator::subdivided_hyper_rectangle (triangulation,
+ repetitions,
+ p1,
+ p2,
+ /*colorize*/ false);
+ // refine left side
+ {
+ typename dealii::hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+ typename dealii::hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (;cell != endc; cell++)
+ if (cell->center()[0] > 1.0) // right
+ {
+ cell->set_active_fe_index(0);
+ }
+ else if (cell->center()[1] > 1.0) // top
+ {
+ cell->set_active_fe_index(2);
+ }
+ else
+ {
+ cell->set_active_fe_index(1);
+ }
+ }
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ // constraints
+ constraints.clear();
+ dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ constraints.close ();
+
+ // interpolate some function
+ dealii::Vector<double> values(dof_handler.n_dofs());
+ dealii::VectorTools::interpolate (dof_handler,
+ func,
+ values);
+
+ dealii::deallog << "dof values:"<<std::endl;
+ for (unsigned int i = 0; i < values.size(); i++)
+ dealii::deallog << " " << values[i];
+ dealii::deallog << std::endl;
+
+ // call Kelly
+ dealii::Vector<float> error(dof_handler.n_dofs());
+ dealii::KellyErrorEstimator<dim>::estimate(dof_handler,
+ face_quadrature_formula,
+ typename dealii::FunctionMap<dim>::type(),
+ values,
+ error,
+ dealii::ComponentMask(),
+ 0,
+ dealii::numbers::invalid_unsigned_int,
+ dealii::numbers::invalid_subdomain_id,
+ dealii::numbers::invalid_material_id,
+ dealii::KellyErrorEstimator<dim>::face_diameter_over_twice_max_degree);
+
+ dealii::deallog <<"error:"<<std::endl;
+ for (unsigned int i = 0; i < error.size(); i++)
+ dealii::deallog << " " << error[i];
+ dealii::deallog << std::endl;
+
+ //output("irregular.vtu",
+ // triangulation,
+ // dof_handler,
+ // values,
+ // error);
+
+ double h,A;
+ get_h_area_sub<dim>(h,A,L);
+ //
+ const double expected_squared_1 = h*A*std::pow(func.get_k(),2)/2.0/std::max(p3,p1);
+ const double expected_squared_2 = h*A*std::pow(func.get_k(),2)/2.0/std::max(p2,p1);
+ const double expected_squared_3 = (dim==2) ?
+ expected_squared_1 + expected_squared_2:
+ 2*expected_squared_1 + 2*expected_squared_2;
+
+ std::vector<double> expected_error(error.size(),0.0);
+
+ expected_error[0] = std::sqrt(expected_squared_3);
+ expected_error[2] = std::sqrt(expected_squared_2);
+ expected_error[4] = std::sqrt(expected_squared_1);
+
+ if (dim ==3)
+ {
+ expected_error[6] = expected_error[2];
+ expected_error[8] = expected_error[4];
+ }
+
+ dealii::deallog << "expected:"<< std::endl;
+ for (unsigned int i = 0; i < expected_error.size(); i++)
+ deallog<<" " << expected_error[i];
+ deallog <<std::endl;
+
+ for (unsigned int i = 0; i < expected_error.size(); i++)
+ AssertThrow (std::fabs(expected_error[i] - error[i] ) < 1e-6, dealii::ExcInternalError());
+
+ dof_handler.clear();
+}
+
+template<int dim>
+class MySecondFunction : public dealii::Function<dim>
+{
+public:
+ MySecondFunction();
+
+ virtual double value(const dealii::Point<dim> & point,
+ const unsigned int component = 0 ) const;
+};
+
+template<int dim>
+MySecondFunction<dim>::MySecondFunction()
+:
+dealii::Function<dim>(1)
+{
+
+}
+
+template<int dim>
+double MySecondFunction<dim>::value(const dealii::Point<dim> & point,
+ const unsigned int ) const
+{
+ double f = 0.0;
+ const double &x = point[0];
+ Assert (dim>1, dealii::ExcNotImplemented());
+ const double &y = point[1];
+
+ return (1.-x)*(1.-y)*(1.-y)+std::pow(1.0-y,4)*std::exp(-x);
+}
+
+template<int dim>
+void test(const MySecondFunction<dim> &func)
+{
+ deallog << std::endl;
+ deallog << "More complicated mesh:"<<std::endl;
+ deallog << "----------------------"<<std::endl;
+
+ dealii::Triangulation<dim> triangulation;
+ dealii::hp::DoFHandler<dim> dof_handler(triangulation);
+ dealii::hp::FECollection<dim> fe_collection;
+ dealii::hp::QCollection<dim> quadrature_formula;
+ dealii::hp::QCollection<dim-1> face_quadrature_formula;
+ dealii::ConstraintMatrix constraints;
+ for (unsigned int p = 1; p <=3; p++)
+ {
+ fe_collection.push_back(dealii::FE_Q<dim>(p));
+ quadrature_formula.push_back(dealii::QGauss<dim>(p+1));
+ face_quadrature_formula.push_back(dealii::QGauss<dim-1>(p+1));
+ }
+ dealii::GridGenerator::hyper_cube (triangulation,0.0,1.0); // reference cell
+
+ // refine
+ {
+ triangulation.refine_global(1);
+
+ // otherwise the next set_active_fe_index
+ // will not carry to the child cells.
+ dof_handler.distribute_dofs (fe_collection);
+
+ typename dealii::hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (;cell != endc; cell++)
+ {
+ bool in_top_left = true;
+ for (unsigned int d=0; d< dim; d++)
+ in_top_left = in_top_left && (cell->center()[d] < 0.5);
+
+ if (in_top_left)
+ {
+ cell->set_active_fe_index(1);
+ cell->set_refine_flag();
+ break;
+ }
+ }
+
+ triangulation.prepare_coarsening_and_refinement();
+
+ triangulation.execute_coarsening_and_refinement();
+
+ cell = dof_handler.begin_active();
+
+ for (;cell != endc; cell++)
+ {
+ if (cell->center()[0] < 0.25)
+ {
+ cell->set_active_fe_index(2);
+ }
+ }
+ }
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ // constraints
+ constraints.clear();
+ dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ constraints.close ();
+
+ // interpolate some function
+ dealii::Vector<double> values(dof_handler.n_dofs());
+
+ dealii::VectorTools::interpolate (dof_handler,
+ func,
+ values);
+
+ dealii::deallog << "dof values:"<<std::endl;
+ for (unsigned int i = 0; i < values.size(); i++)
+ dealii::deallog << " " << values[i];
+ dealii::deallog << std::endl;
+
+ // call Kelly
+ dealii::Vector<float> error(dof_handler.n_dofs());
+ dealii::KellyErrorEstimator<dim>::estimate(dof_handler,
+ face_quadrature_formula,
+ typename dealii::FunctionMap<dim>::type (),
+ values,
+ error,
+ dealii::ComponentMask(),
+ 0,
+ dealii::numbers::invalid_unsigned_int,
+ dealii::numbers::invalid_subdomain_id,
+ dealii::numbers::invalid_material_id,
+ dealii::KellyErrorEstimator<dim>::face_diameter_over_twice_max_degree);
+
+ dealii::deallog <<"error:"<<std::endl;
+ for (unsigned int i = 0; i < error.size(); i++)
+ dealii::deallog << " " << error[i];
+ dealii::deallog << std::endl;
+
+ dof_handler.clear();
+}
+
+
+int main ()
+{
+ std::ofstream logfile("output");
+ dealii::deallog.attach(logfile);
+ dealii::deallog.depth_console(0);
+ dealii::deallog.threshold_double(1e-8);
+
+ {
+ NeumanBC<2> func(1.25);
+ test_neumann(func);
+ }
+
+ {
+ MyFunction<2> func(0.25);
+ test_regular(func);
+ }
+
+ {
+ MyFunction<2> func(0.75);
+ test_irregular(func);
+ }
+
+ deallog << "===3d==="<<std::endl;
+ {
+ NeumanBC<3> func(1.25);
+ test_neumann(func);
+ }
+
+ {
+ MyFunction<3> func(0.25);
+ test_regular(func);
+ }
+
+ {
+ MyFunction<3> func(0.75);
+ test_irregular(func);
+ }
+
+ {
+ MySecondFunction<2> function;
+ test(function);
+ }
+
+
+
+ dealii::deallog << "Ok"<<std::endl;
+
+}
--- /dev/null
+
+DEAL::NeumanBC case:
+DEAL::--------------
+DEAL::dof values:
+DEAL:: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::error:
+DEAL:: 1.44338
+DEAL::expected:
+DEAL:: 1.44338
+DEAL::
+DEAL::Regular face:
+DEAL::-------------
+DEAL::dof values:
+DEAL:: 0 0 0 0 0.250000 0.250000 0 0.250000 0.125000 0.125000 0.125000
+DEAL::error:
+DEAL:: 0.250000 0.250000
+DEAL::expected:
+DEAL:: 0.250000
+DEAL::
+DEAL::Irregular face:
+DEAL::---------------
+DEAL::dof values:
+DEAL:: 0 0.750000 0 0.750000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::error:
+DEAL:: 0.484123 0.00000 0.375000 0.00000 0.306186
+DEAL::expected:
+DEAL:: 0.484123 0 0.375000 0 0.306186
+DEAL::===3d===
+DEAL::NeumanBC case:
+DEAL::--------------
+DEAL::dof values:
+DEAL:: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::error:
+DEAL:: 2.42746
+DEAL::expected:
+DEAL:: 2.42746
+DEAL::
+DEAL::Regular face:
+DEAL::-------------
+DEAL::dof values:
+DEAL:: 0 0 0 0 0 0 0 0 0.250000 0.250000 0.250000 0.250000 0 0.250000 0.125000 0.125000 0 0.250000 0.125000 0.125000 0 0.250000 0 0.250000 0 0.250000 0.125000 0.125000 0.125000 0.125000 0.125000
+DEAL::error:
+DEAL:: 0.420448 0.420448
+DEAL::expected:
+DEAL:: 0.420448
+DEAL::
+DEAL::Irregular face:
+DEAL::---------------
+DEAL::dof values:
+DEAL:: 0 0.750000 0 0.750000 0 0.750000 0 0.750000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::error:
+DEAL:: 0.814194 0.00000 0.445953 0.00000 0.364119 0.00000 0.445953 0.00000 0.364119
+DEAL::expected:
+DEAL:: 0.814194 0 0.445953 0 0.364119 0 0.445953 0 0.364119
+DEAL::
+DEAL::More complicated mesh:
+DEAL::----------------------
+DEAL::dof values:
+DEAL:: 1.10653 0.367879 0.162908 0.0229925 0.312500 0 0 0 2.00000 0.878906 1.54634 1.17670 1.18009 0.896412 1.83671 1.67982 0.806733 0.736582 1.41987 1.29790 1.08027 0.986922 1.52880 0.668292 0.473160 1.03074 0.738350 1.31229 0.569025 0.881392 0.641975 0.456067 0.487171 0.345385 0.286669 0.261238 0.589145 0.537577 0.418452 0.381578 0.236175 0.411804 0.287862 0.199206 0.349013
+DEAL::error:
+DEAL:: 0.401222 0.102290 0.400943 0.000678580 0.0104936 0.0788644 0.0660732
+DEAL::Ok
projects / dealii.git / commitdiff