+/**
+ * This is a specialization of the general template for 1d. The implementation
+ * is sufficiently different for 1d to justify this specialization. The basic
+ * difference between 1d and all other space dimensions is that in 1d, there
+ * are no faces of cells, just the vertices between line segments, so we have
+ * to compute the jump terms differently. However, this class offers exactly
+ * the same public functions as the general template, so that a user will not
+ * see any difference.
+ *
+ * @author Wolfgang Bangerth, 1998, 2004.
+ */
+template <>
+class KellyErrorEstimator<1>
+{
+ public:
+ /**
+ * Implementation of the error
+ * estimator described above. You
+ * may give a coefficient, but
+ * there is a default value which
+ * denotes the constant
+ * coefficient with value
+ * one. The coefficient function
+ * may either be a scalar one, in
+ * which case it is used for all
+ * components of the finite
+ * element, or a vector-valued
+ * one with as many components as
+ * there are in the finite
+ * element; in the latter case,
+ * each component is weighted by
+ * the respective component in
+ * the coefficient.
+ *
+ * You might give a list of
+ * components you want to
+ * evaluate, in case the finite
+ * element used by the
+ * @ref{DoFHandler} object is
+ * vector-valued. You then have
+ * to set those entries to true
+ * in the bit-vector
+ * @p{component_mask} for which the
+ * respective component is to be
+ * used in the error
+ * estimator. The default is to
+ * use all components, which is
+ * done by either providing a
+ * bit-vector with all-set
+ * entries, or an empty
+ * bit-vector.
+ *
+ * The estimator supports
+ * multithreading and splits the
+ * cells to
+ * @p{multithread_info.n_default_threads}
+ * (default) threads. The number
+ * of threads to be used in
+ * multithreaded mode can be set
+ * with the last parameter of the
+ * error estimator.
+ * Multithreading is only
+ * implemented in two and three
+ * dimensions.
+ */
+ template <typename InputVector>
+ static void estimate (const Mapping<1> &mapping,
+ const DoFHandler<1> &dof,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask = std::vector<bool>(),
+ const Function<1> *coefficients = 0,
+ const unsigned int n_threads = multithread_info.n_default_threads);
+
+ /**
+ * Calls the @p{estimate}
+ * function, see above, with
+ * @p{mapping=MappingQ1<1>()}.
+ */
+ template <typename InputVector>
+ static void estimate (const DoFHandler<1> &dof,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask = std::vector<bool>(),
+ const Function<1> *coefficients = 0,
+ const unsigned int n_threads = multithread_info.n_default_threads);
+
+ /**
+ * Same function as above, but
+ * accepts more than one solution
+ * vectors and returns one error
+ * vector for each solution
+ * vector. For the reason of
+ * existence of this function,
+ * see the general documentation
+ * of this class.
+ *
+ * Since we do not want to force
+ * the user of this function to
+ * copy around their solution
+ * vectors, the vector of
+ * solution vectors takes
+ * pointers to the solutions,
+ * rather than being a vector of
+ * vectors. This makes it simpler
+ * to have the solution vectors
+ * somewhere in memory, rather
+ * than to have them collected
+ * somewhere special. (Note that
+ * it is not possible to
+ * construct of vector of
+ * references, so we had to use a
+ * vector of pointers.)
+ */
+ template <typename InputVector>
+ static void estimate (const Mapping<1> &mapping,
+ const DoFHandler<1> &dof,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const std::vector<const InputVector *> &solutions,
+ std::vector<Vector<float>*> &errors,
+ const std::vector<bool> &component_mask = std::vector<bool>(),
+ const Function<1> *coefficients = 0,
+ const unsigned int n_threads = multithread_info.n_default_threads);
+
+ /**
+ * Calls the @p{estimate}
+ * function, see above, with
+ * @p{mapping=MappingQ1<1>()}.
+ */
+ template <typename InputVector>
+ static void estimate (const DoFHandler<1> &dof,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const std::vector<const InputVector *> &solutions,
+ std::vector<Vector<float>*> &errors,
+ const std::vector<bool> &component_mask = std::vector<bool>(),
+ const Function<1> *coefficients = 0,
+ const unsigned int n_threads = multithread_info.n_default_threads);
+
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidBoundaryIndicator);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidComponentMask);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidCoefficient);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidBoundaryFunction);
+ /**
+ * Exception
+ */
+ DeclException2 (ExcIncompatibleNumberOfElements,
+ int, int,
+ << "The number of elements " << arg1 << " and " << arg2
+ << " of the vectors do not match!");
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidSolutionVector);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNoSolutions);
+};
+
+
+
#endif
}
-
-template <int dim>
-KellyErrorEstimator<dim>::PerThreadData::
-PerThreadData (const unsigned int n_solution_vectors,
- const unsigned int n_components,
- const unsigned int n_q_points)
-{
- // Init the size of a lot of vectors
- // needed in the calculations once
- // per thread.
- normal_vectors.resize(n_q_points);
- coefficient_values1.resize(n_q_points);
- coefficient_values.resize(n_q_points);
- JxW_values.resize(n_q_points);
-
- phi.resize(n_solution_vectors);
- psi.resize(n_solution_vectors);
- neighbor_psi.resize(n_solution_vectors);
-
- for (unsigned int i=0; i<n_solution_vectors; ++i)
- {
- phi[i].resize(n_q_points);
- psi[i].resize(n_q_points);
- neighbor_psi[i].resize(n_q_points);
-
- for (unsigned int qp=0;qp<n_q_points;++qp)
- {
- phi[i][qp].resize(n_components);
- psi[i][qp].resize(n_components);
- neighbor_psi[i][qp].resize(n_components);
- };
- };
-
- for (unsigned int qp=0;qp<n_q_points;++qp)
- coefficient_values[qp].reinit(n_components);
-}
-
-
+#if deal_II_dimension == 1
-// the following function is still independent of dimension, but it
-// calls dimension dependent functions
-template <int dim>
template <typename InputVector>
-void KellyErrorEstimator<dim>::estimate (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof_handler,
- const Quadrature<dim-1> &quadrature,
- const typename FunctionMap<dim>::type &neumann_bc,
- const InputVector &solution,
- Vector<float> &error,
- const std::vector<bool> &component_mask,
- const Function<dim> *coefficients,
- const unsigned int n_threads)
+void KellyErrorEstimator<1>::estimate (const Mapping<1> &mapping,
+ const DoFHandler<1> &dof_handler,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<1> *coefficients,
+ const unsigned int n_threads)
{
// just pass on to the other function
const std::vector<const InputVector *> solutions (1, &solution);
}
-template <int dim>
template <typename InputVector>
-void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim> &dof_handler,
- const Quadrature<dim-1> &quadrature,
- const typename FunctionMap<dim>::type &neumann_bc,
- const InputVector &solution,
- Vector<float> &error,
- const std::vector<bool> &component_mask,
- const Function<dim> *coefficients,
- const unsigned int n_threads)
+void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<1> *coefficients,
+ const unsigned int n_threads)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- static const MappingQ1<dim> mapping;
+ static const MappingQ1<1> mapping;
estimate(mapping, dof_handler, quadrature, neumann_bc, solution,
error, component_mask, coefficients, n_threads);
}
-#if deal_II_dimension == 1
-
-template <>
template <typename InputVector>
-void KellyErrorEstimator<1>::
-estimate_some (const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- const FunctionMap<1>::type &,
- const std::vector<const InputVector *> &,
- const std::vector<bool> &,
- const Function<1> *,
- const std::pair<unsigned int, unsigned int> ,
- FaceIntegrals &,
- PerThreadData &)
+void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const std::vector<const InputVector*> &solutions,
+ std::vector<Vector<float>*> &errors,
+ const std::vector<bool> &component_mask,
+ const Function<1> *coefficients,
+ const unsigned int n_threads)
{
- // in 1d, the @p{estimate} function
- // does all the work
- Assert (false, ExcInternalError() );
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ static const MappingQ1<1> mapping;
+ estimate(mapping, dof_handler, quadrature, neumann_bc, solutions,
+ errors, component_mask, coefficients, n_threads);
}
-template <>
template <typename InputVector>
void KellyErrorEstimator<1>::
estimate (const Mapping<1> &mapping,
// vertices of each cell.
QTrapez<1> quadrature;
FEValues<1> fe_values (mapping, dof_handler.get_fe(), quadrature, update_gradients);
- active_cell_iterator cell = dof_handler.begin_active();
+ DoFHandler<1>::active_cell_iterator cell = dof_handler.begin_active();
for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index)
{
for (unsigned int n=0; n<n_solution_vectors; ++n)
#else // #if deal_II_dimension !=1
+template <int dim>
+KellyErrorEstimator<dim>::PerThreadData::
+PerThreadData (const unsigned int n_solution_vectors,
+ const unsigned int n_components,
+ const unsigned int n_q_points)
+{
+ // Init the size of a lot of vectors
+ // needed in the calculations once
+ // per thread.
+ normal_vectors.resize(n_q_points);
+ coefficient_values1.resize(n_q_points);
+ coefficient_values.resize(n_q_points);
+ JxW_values.resize(n_q_points);
+
+ phi.resize(n_solution_vectors);
+ psi.resize(n_solution_vectors);
+ neighbor_psi.resize(n_solution_vectors);
+
+ for (unsigned int i=0; i<n_solution_vectors; ++i)
+ {
+ phi[i].resize(n_q_points);
+ psi[i].resize(n_q_points);
+ neighbor_psi[i].resize(n_q_points);
+
+ for (unsigned int qp=0;qp<n_q_points;++qp)
+ {
+ phi[i][qp].resize(n_components);
+ psi[i][qp].resize(n_components);
+ neighbor_psi[i][qp].resize(n_components);
+ };
+ };
+
+ for (unsigned int qp=0;qp<n_q_points;++qp)
+ coefficient_values[qp].reinit(n_components);
+}
+
+
+
+
+// the following function is still independent of dimension, but it
+// calls dimension dependent functions
+template <int dim>
+template <typename InputVector>
+void KellyErrorEstimator<dim>::estimate (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof_handler,
+ const Quadrature<dim-1> &quadrature,
+ const typename FunctionMap<dim>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<dim> *coefficients,
+ const unsigned int n_threads)
+{
+ // 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);
+}
+
+
+template <int dim>
+template <typename InputVector>
+void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim> &dof_handler,
+ const Quadrature<dim-1> &quadrature,
+ const typename FunctionMap<dim>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<dim> *coefficients,
+ const unsigned int n_threads)
+{
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ static const MappingQ1<dim> mapping;
+ estimate(mapping, dof_handler, quadrature, neumann_bc, solution,
+ error, component_mask, coefficients, n_threads);
+}
+
+
+
+
template <int dim>
template <typename InputVector>
void KellyErrorEstimator<dim>::
};
}
-#endif
template <int dim>
template <typename InputVector>
-#if deal_II_dimension == 1
-
-template <>
-template <typename InputVector>
-void KellyErrorEstimator<1>::
-integrate_over_regular_face (const DoFHandler<1> &,
- const Quadrature<0> &,
- const FunctionMap<1>::type &,
- const std::vector<const InputVector *> &,
- const std::vector<bool> &,
- const Function<1> *,
- FaceIntegrals &,
- PerThreadData &,
- const active_cell_iterator &,
- const unsigned int ,
- FEFaceValues<1> &,
- FEFaceValues<1> &)
-{
- Assert (false, ExcInternalError());
-}
-
-
-template <>
-template <typename InputVector>
-void KellyErrorEstimator<1>::
-integrate_over_irregular_face (const DoFHandler<1> &,
- const Quadrature<0> &,
- const std::vector<const InputVector *> &,
- const std::vector<bool> &,
- const Function<1> *,
- FaceIntegrals &,
- PerThreadData &,
- const active_cell_iterator &,
- const unsigned int ,
- FEFaceValues<1> &,
- FESubfaceValues<1> &)
-{
- Assert (false, ExcInternalError());
-}
-
-#endif
-
-
template <int dim>
template <typename InputVector>
void KellyErrorEstimator<dim>::
face_integrals[face] = sum;
}
+#endif
+
+
// explicit instantiations
+#if deal_II_dimension != 1
template class KellyErrorEstimator<deal_II_dimension>;
+#endif
// instantiate the externally visible functions. define a list of functions
// for vectors, where the vector/matrix can be replaced using a preprocessor