* integrals computed here show superconvergence properties, i.e. they tend
* to zero faster than the error itself, thus ruling out the values as error
* indicators.
- *
+ *
+ * The error estimator really only estimates the error for the generalized
+ * Poisson equation $-\nabla\cdot a(x) \nabla u = f$ with either Dirichlet
+ * boundary conditions or generalized Neumann boundary conditions involving
+ * the conormal derivative $a\frac{du}{dn} = g$.
+ *
* The error estimator returns a vector of estimated errors per cell which
* can be used to feed the #Triangulation<dim>::refine_*# functions.
*
*
* In principle, the implementation of the error estimation is simple: let
* $$ \eta_K^2 =
- * \frac h{24} \int_{\partial K} \left[\frac{\partial u_h}{\partial n}\right]^2 do
+ * \frac h{24} \int_{\partial K} \left[a \frac{\partial u_h}{\partial n}\right]^2 do
* $$
* 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$,
*
* \item The face belongs to a Neumann boundary. In this case, the
* contribution of the face $F\in\partial K$ looks like
- * $$ \int_F \left|g-\frac{\partial u_h}{\partial n}\right| ds $$
+ * $$ \int_F \left|g-a\frac{\partial u_h}{\partial n}\right|^2 ds $$
* where $g$ is the Neumann boundary function.
*
* \item No other boundary conditions are considered.
* guaranteed by the #map# data type.
*/
typedef map<unsigned char,const Function<dim>*> FunctionMap;
-
+
+ /**
+ * 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.
+ */
static void estimate_error (const DoFHandler<dim> &dof,
const Quadrature<dim-1> &quadrature,
const FiniteElement<dim> &fe,
const Boundary<dim> &boundary,
const FunctionMap &neumann_bc,
const dVector &solution,
- dVector &error);
+ dVector &error,
+ const Function<dim> *coefficient = 0);
/**
* Exception
FEFaceValues<dim> &fe_face_values_cell,
FEFaceValues<dim> &fe_face_values_neighbor,
FaceIntegrals &face_integrals,
- const dVector &solution);
+ const dVector &solution,
+ const Function<dim> *coefficient);
/**
* The same applies as for the function
FEFaceValues<dim> &fe_face_values,
FESubfaceValues<dim> &fe_subface_values,
FaceIntegrals &face_integrals,
- const dVector &solution);
+ const dVector &solution,
+ const Function<dim> *coefficient);
};
const Boundary<1> &,
const FunctionMap &,
const dVector &,
- dVector &) {
+ dVector &,
+ const Function<1> *) {
Assert(false, ExcNotImplemented());
};
const Boundary<dim> &boundary,
const FunctionMap &neumann_bc,
const dVector &solution,
- dVector &error) {
+ dVector &error,
+ const Function<dim> *coefficient) {
Assert (neumann_bc.find(255) == neumann_bc.end(),
ExcInvalidBoundaryIndicator());
fe_face_values_cell,
fe_face_values_neighbor,
face_integrals,
- solution);
+ solution,
+ coefficient);
else
// otherwise we need to do some
// special computations which do
n_q_points,
fe_face_values_cell,
fe_subface_values,
- face_integrals, solution);
+ face_integrals, solution,
+ coefficient);
};
FEFaceValues<1> &,
FEFaceValues<1> &,
FaceIntegrals &,
- const dVector &) {
+ const dVector &,
+ const Function<1> *) {
Assert (false, ExcInternalError());
};
FEFaceValues<1> &,
FESubfaceValues<1> &,
FaceIntegrals &,
- const dVector &) {
+ const dVector &,
+ const Function<1> *) {
Assert (false, ExcInternalError());
};
FEFaceValues<dim> &fe_face_values_cell,
FEFaceValues<dim> &fe_face_values_neighbor,
FaceIntegrals &face_integrals,
- const dVector &solution) {
+ const dVector &solution,
+ const Function<dim> *coefficient) {
const DoFHandler<dim>::face_iterator face = cell->face(face_no);
// initialize data of the restriction
// points
//
// let psi be a short name for
- // [grad u_h]
+ // [a grad u_h]
vector<Point<dim> > psi(n_q_points);
fe_face_values_cell.get_function_grads (solution, psi);
for (unsigned int point=0; point<n_q_points; ++point)
phi[point] = psi[point]*normal_vectors[point];
-
+
if (face->at_boundary() == true)
// neumann boundary face. compute
phi[point] -= g[point];
};
+ // if a coefficient was given: use that
+ // to scale the jump in the gradient
+ if (coefficient != 0)
+ {
+ vector<double> coefficient_values (n_q_points);
+ coefficient->value_list (fe_face_values_cell.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ phi[point] *= coefficient_values[point];
+ };
+
// now phi contains the following:
- // - for an internal face, phi=[du/dn]
+ // - for an internal face, phi=[a du/dn]
// - for a neumann boundary face,
- // phi=du/dn-g
+ // phi=a du/dn-g
// each component being the
// mentioned value at one of the
// quadrature points
FEFaceValues<dim> &fe_face_values,
FESubfaceValues<dim> &fe_subface_values,
FaceIntegrals &face_integrals,
- const dVector &solution) {
+ const dVector &solution,
+ const Function<dim> *coefficient) {
const DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
Assert (neighbor.state() == valid, ExcInternalError());
Assert (neighbor->has_children(), ExcInternalError());
// points
//
// let psi be a short name for
- // [grad u_h]
+ // [a grad u_h]
vector<Point<dim> > psi(n_q_points);
// store which number #cell# has in the
psi.begin(),
minus<Point<dim> >());
- // next we have to multiply this with
- // the normal vector. Since we have
- // taken the difference of gradients
- // for internal faces, we may chose
- // the normal vector of one cell,
- // taking that of the neighbor
- // would only change the sign. We take
- // the outward normal.
- //
- // let phi be the name of the integrand
+ // next we have to multiply this with
+ // the normal vector. Since we have
+ // taken the difference of gradients
+ // for internal faces, we may chose
+ // the normal vector of one cell,
+ // taking that of the neighbor
+ // would only change the sign. We take
+ // the outward normal.
+ //
+ // let phi be the name of the integrand
vector<double> phi(n_q_points,0);
const vector<Point<dim> > &normal_vectors(fe_face_values.
get_normal_vectors());
for (unsigned int point=0; point<n_q_points; ++point)
phi[point] = psi[point]*normal_vectors[point];
- // take the square of the phi[i]
+
+ // if a coefficient was given: use that
+ // to scale the jump in the gradient
+ if (coefficient != 0)
+ {
+ vector<double> coefficient_values (n_q_points);
+ coefficient->value_list (fe_face_values.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ phi[point] *= coefficient_values[point];
+ };
+
+ // take the square of the phi[i]
// for integration
transform (phi.begin(), phi.end(),
phi.begin(), ptr_fun(sqr));