* coefficient, but there is a default
* value which denotes the constant
* coefficient with value one.
+ *
+ * You must give the component if the
+ * finite element in use by the #dof#
+ * object has more than one component.
+ * This number shall be between zero
+ * and the number of components within
+ * the finite element. If the finite
+ * element has only one component,
+ * then the parameter selecting the
+ * component shall be zero, which is
+ * also the default value.
*/
static void estimate (const DoFHandler<dim> &dof,
const Quadrature<dim-1> &quadrature,
const FunctionMap &neumann_bc,
const Vector<double> &solution,
Vector<float> &error,
- const Function<dim> *coefficient = 0);
+ const Function<dim> *coefficient = 0,
+ const unsigned int selected_component = 0);
/**
* Exception
* Exception
*/
DeclException0 (ExcInvalidBoundaryIndicator);
-
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidComponent,
+ int, int,
+ << "The component you gave (" << arg1 << ") "
+ << "is invalid with respect to the number "
+ << "of components in the finite element "
+ << "(" << arg2 << ")");
private:
/**
* Declare a data type to represent the
FEFaceValues<dim> &fe_face_values_neighbor,
FaceIntegrals &face_integrals,
const Vector<double>&solution,
+ const unsigned int n_components,
+ const unsigned int selected_component,
const Function<dim> *coefficient);
/**
FESubfaceValues<dim> &fe_subface_values,
FaceIntegrals &face_integrals,
const Vector<double> &solution,
- const Function<dim> *coefficient);
+ const unsigned int n_components,
+ const unsigned int selected_component,
+ const Function<dim> *coefficient);
};
const FunctionMap &,
const Vector<double> &,
Vector<float> &,
- const Function<1> *) {
+ const Function<1> *,
+ const unsigned int) {
Assert(false, ExcNotImplemented());
};
const FunctionMap &neumann_bc,
const Vector<double> &solution,
Vector<float> &error,
- const Function<dim> *coefficient)
+ const Function<dim> *coefficient,
+ const unsigned int selected_component)
{
Assert (neumann_bc.find(255) == neumann_bc.end(),
ExcInvalidBoundaryIndicator());
+ Assert (selected_component < dof.get_fe().n_components,
+ ExcInvalidComponent (selected_component, dof.get_fe().n_components));
// create a map of integrals indexed by
// the corresponding face. In this map
fe_face_values_neighbor,
face_integrals,
solution,
+ dof.get_fe().n_components,
+ selected_component,
coefficient);
else
// otherwise we need to do some
fe_face_values_cell,
fe_subface_values,
face_integrals, solution,
+ dof.get_fe().n_components,
+ selected_component,
coefficient);
};
FEFaceValues<1> &,
FaceIntegrals &,
const Vector<double> &,
+ const unsigned int ,
+ const unsigned int ,
const Function<1> *) {
Assert (false, ExcInternalError());
};
FESubfaceValues<1> &,
FaceIntegrals &,
const Vector<double> &,
+ const unsigned int ,
+ const unsigned int ,
const Function<1> *) {
Assert (false, ExcInternalError());
};
FEFaceValues<dim> &fe_face_values_neighbor,
FaceIntegrals &face_integrals,
const Vector<double> &solution,
+ const unsigned int n_components,
+ const unsigned int selected_component,
const Function<dim> *coefficient) {
const DoFHandler<dim>::face_iterator face = cell->face(face_no);
// let psi be a short name for
// [a grad u_h]
vector<Tensor<1,dim> > psi(n_q_points);
- fe_face_values_cell.get_function_grads (solution, psi);
+ if (n_components == 1)
+ fe_face_values_cell.get_function_grads (solution, psi);
+ else
+ {
+ vector<vector<Tensor<1,dim> > > tmp (n_components,
+ psi);
+ fe_face_values_cell.get_function_grads (solution, tmp);
+
+ psi.swap (tmp[selected_component]);
+ };
+
// now compute over the other side of
// the finite element solution
// restricted to the neighbor cell
vector<Tensor<1,dim> > neighbor_psi (n_q_points);
- fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
+ if (n_components == 1)
+ fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
+ else
+ {
+ vector<vector<Tensor<1,dim> > > tmp (n_components,
+ neighbor_psi);
+ fe_face_values_neighbor.get_function_grads (solution, tmp);
+ neighbor_psi.swap (tmp[selected_component]);
+ };
+
// compute the jump in the gradients
transform (psi.begin(), psi.end(),
FESubfaceValues<dim> &fe_subface_values,
FaceIntegrals &face_integrals,
const Vector<double> &solution,
+ const unsigned int n_components,
+ const unsigned int selected_component,
const Function<dim> *coefficient) {
const DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
Assert (neighbor.state() == valid, ExcInternalError());
// store the gradient of the solution
// in psi
fe_subface_values.reinit (cell, face_no, subface_no);
- fe_subface_values.get_function_grads (solution, psi);
+ if (n_components == 1)
+ fe_subface_values.get_function_grads (solution, psi);
+ else
+ {
+ vector<vector<Tensor<1,dim> > > tmp (n_components,
+ psi);
+ fe_subface_values.get_function_grads (solution, tmp);
+ psi.swap (tmp[selected_component]);
+ };
// restrict the finite element on the
// neighbor cell to the common #subface#.
// store the gradient in #neighbor_psi#
vector<Tensor<1,dim> > neighbor_psi (n_q_points);
fe_face_values.reinit (neighbor_child, neighbor_neighbor);
- fe_face_values.get_function_grads (solution, neighbor_psi);
+ if (n_components == 1)
+ fe_face_values.get_function_grads (solution, neighbor_psi);
+ else
+ {
+ vector<vector<Tensor<1,dim> > > tmp (n_components,
+ neighbor_psi);
+ fe_face_values.get_function_grads (solution, tmp);
+ neighbor_psi.swap (tmp[selected_component]);
+ };
// compute the jump in the gradients
transform (psi.begin(), psi.end(),