passed through the constructor. In debug mode, the accessor functions, which
return values from the different fields, check whether the required field
was initialized, thus avoiding use of unitialized data.
+
+
+ {\bf Member functions}
+
+ The functions of this class fall into different cathegories:
+ \begin{itemize}
+ \item #shape_value#, #shape_grad#, etc: return one of the values
+ of this object at a time. In many cases you will want to get
+ a whole bunch at a time for performance or convenience reasons,
+ then use the #get_*# functions.
+
+ \item #get_shape_values#, #get_shape_grads#, etc: these return
+ a reference to a whole field. Usually these fields contain
+ the values of all ansatz functions at all quadrature points.
+
+ \item #get_function_values#, #get_function_gradients#: these
+ two functions offer a simple way to avoid the detour of the
+ ansatz functions, if you have a finite solution (resp. the
+ vector of values associated with the different ansatz functions.)
+ Then you may want to get information from the restriction of
+ the finite element function to a certain cell, e.g. the values
+ of the function at the quadrature points or the values of its
+ gradient. These two functions provide the information needed:
+ you pass it a vector holding the finite element solution and the
+ functions return the values or gradients of the finite element
+ function restricted to the cell which was given last time the
+ #reinit# function was given.
+
+ Though possible in principle, these functions do not call the
+ #reinit# function, you have to do so yourself beforehand. On the
+ other hand, a copy of the cell iterator is stored which was used
+ last time the #reinit# function was called. This frees us from
+ the need to pass the cell iterator again to these two functions,
+ which guarantees that the cell used here is in sync with that used
+ for the #reinit# function. You should, however, make sure that
+ nothing substantial happens to the #DoFHandler# object or any
+ other involved instance between the #reinit# and the #get_function_*#
+ functions are called.
+
+ \item #reinit#: initialize the #FEValues# object for a certain cell.
+ See above for more information.
+ \end{itemize}
*/
template <int dim>
class FEValues {
*/
const dFMatrix & get_shape_values () const;
+ /**
+ * Return the values of the finite
+ * element function characterized
+ * by #fe_function# restricted to
+ * #cell# at the quadrature points.
+ *
+ * The function assumes that the
+ * #values# object already has the
+ * right size.
+ */
+ void get_function_values (const dVector &fe_function,
+ vector<double> &values) const;
+
/**
* Return the gradient of the #i#th shape
* function at the #j# quadrature point.
*/
const vector<vector<Point<dim> > > & get_shape_grads () const;
+ /**
+ * Return the gradients of the finite
+ * element function characterized
+ * by #fe_function# restricted to
+ * #cell# at the quadrature points.
+ *
+ * The function assumes that the
+ * #gradients# object already has the
+ * right size.
+ */
+ void get_function_grads (const dVector &fe_function,
+ vector<Point<dim> > &gradients) const;
+
/**
* Return the position of the #i#th
* quadrature point in real space.
* Exception
*/
DeclException0 (ExcCannotInitializeField);
+ /**
+ * Exception
+ */
+ DeclException2 (ExcWrongVectorSize,
+ int, int,
+ << "Vector has wrong size " << arg1
+ << ", expected size " << arg2);
private:
/**
* the reinit function.
*/
UpdateFlags update_flags;
+
+ /**
+ * Store the cell selected last time
+ * the #reinit# function was called
+ * to make access
+ * to the #get_function_*# functions
+ * safer.
+ */
+ DoFHandler<dim>::cell_iterator present_cell;
};
* documentation for the matrix itself.
*/
const dFMatrix & get_shape_values () const;
-
+
+ /**
+ * Return the values of the finite
+ * element function characterized
+ * by #fe_function# restricted to
+ * #cell# at the quadrature points.
+ *
+ * The function assumes that the
+ * #values# object already has the
+ * right size.
+ */
+ void get_function_values (const dVector &fe_function,
+ vector<double> &values) const;
+
/**
* Return the gradient of the #i#th shape
* function at the #j# quadrature point.
*/
const vector<vector<Point<dim> > > & get_shape_grads () const;
+ /**
+ * Return the gradients of the finite
+ * element function characterized
+ * by #fe_function# restricted to
+ * #cell# at the quadrature points.
+ *
+ * The function assumes that the
+ * #gradients# object already has the
+ * right size.
+ */
+ void get_function_grads (const dVector &fe_function,
+ vector<Point<dim> > &gradients) const;
+
/**
* Return the position of the #i#th
* quadrature point in real space.
* Exception
*/
DeclException0 (ExcNotImplemented);
-
+ /**
+ * Exception
+ */
+ DeclException2 (ExcWrongVectorSize,
+ int, int,
+ << "Vector has wrong size " << arg1
+ << ", expected size " << arg2);
private:
/**
* Store the values of the shape functions
*/
UpdateFlags update_flags;
+ /**
+ * Store the cell selected last time
+ * the #reinit# function was called
+ * to make access
+ * to the #get_function_*# functions
+ * safer.
+ */
+ DoFHandler<dim>::cell_iterator present_cell;
+
/**
* Store the number of the face selected
* last time the #reinit# function was
+template <int dim>
+void FEValues<dim>::get_function_values (const dVector &fe_function,
+ vector<double> &values) const {
+ Assert (values.size() == n_quadrature_points,
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ // get function values of dofs
+ // on this cell
+ vector<double> dof_values (total_dofs, 0);
+ present_cell->get_dof_values (fe_function, dof_values);
+
+ // initialize with zero
+ fill_n (values.begin(), n_quadrature_points, 0);
+
+ // add up contributions of ansatz
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+ values[point] += (dof_values[shape_func] *
+ shape_values(shape_func, point));
+};
+
+
+
template <int dim>
const Point<dim> &
FEValues<dim>::shape_grad (const unsigned int i,
+template <int dim>
+void FEValues<dim>::get_function_grads (const dVector &fe_function,
+ vector<Point<dim> > &gradients) const {
+ Assert (gradients.size() == n_quadrature_points,
+ ExcWrongVectorSize(gradients.size(), n_quadrature_points));
+
+ // get function values of dofs
+ // on this cell
+ vector<double> dof_values (total_dofs, 0);
+ present_cell->get_dof_values (fe_function, dof_values);
+
+ // initialize with zero
+ fill_n (gradients.begin(), n_quadrature_points, Point<dim>());
+
+ // add up contributions of ansatz
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+ gradients[point] += (dof_values[shape_func] *
+ shape_gradients[shape_func][point]);
+};
+
+
+
template <int dim>
const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
const FiniteElement<dim> &fe,
const Boundary<dim> &boundary) {
+ present_cell = cell;
// fill jacobi matrices and real
// quadrature points
if ((update_flags & update_jacobians) ||
case 0:
global_unit_quadrature_points[face][p]
= Point<dim>(unit_quadrature_points[p](0),0);
- cout << " face=" << face << ", p=" << p
- << ", uqp[p]=" << unit_quadrature_points[p]
- << ", guqp[face][p]="
- << global_unit_quadrature_points[face][p]
- << ", ??="
- << Point<dim>(unit_quadrature_points[p](0),0)
- << endl;
break;
case 1:
global_unit_quadrature_points[face][p]
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
{
- shape_values[face](i,j) = fe.shape_value(i, global_unit_quadrature_points[face][j]);
+ shape_values[face](i,j)
+ = fe.shape_value(i, global_unit_quadrature_points[face][j]);
unit_shape_gradients[face][i][j]
= fe.shape_grad(i, global_unit_quadrature_points[face][j]);
};
+template <int dim>
+void FEFaceValues<dim>::get_function_values (const dVector &fe_function,
+ vector<double> &values) const {
+ Assert (values.size() == n_quadrature_points,
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ // get function values of dofs
+ // on this cell
+ vector<double> dof_values (total_dofs, 0);
+ present_cell->get_dof_values (fe_function, dof_values);
+
+ // initialize with zero
+ fill_n (values.begin(), n_quadrature_points, 0);
+
+ // add up contributions of ansatz
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+ values[point] += (dof_values[shape_func] *
+ shape_values[selected_face](shape_func, point));
+};
+
+
+
template <int dim>
const Point<dim> &
FEFaceValues<dim>::shape_grad (const unsigned int i,
ExcInvalidIndex (i, shape_values[selected_face].m()));
Assert (j<shape_values[selected_face].n(),
ExcInvalidIndex (j, shape_values[selected_face].n()));
- Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
+ Assert (update_flags & update_gradients,
+ ExcAccessToUninitializedField());
return shape_gradients[i][j];
};
+template <int dim>
+void FEFaceValues<dim>::get_function_grads (const dVector &fe_function,
+ vector<Point<dim> > &gradients) const {
+ Assert (gradients.size() == n_quadrature_points,
+ ExcWrongVectorSize(gradients.size(), n_quadrature_points));
+
+ // get function values of dofs
+ // on this cell
+ vector<double> dof_values (total_dofs, 0);
+ present_cell->get_dof_values (fe_function, dof_values);
+
+ // initialize with zero
+ fill_n (gradients.begin(), n_quadrature_points, Point<dim>());
+
+ // add up contributions of ansatz
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+ gradients[point] += (dof_values[shape_func] *
+ shape_gradients[shape_func][point]);
+};
+
+
+
template <int dim>
const Point<dim> & FEFaceValues<dim>::quadrature_point (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
- Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
+ Assert (update_flags & update_q_points,
+ ExcAccessToUninitializedField());
return quadrature_points[i];
};
template <int dim>
const Point<dim> & FEFaceValues<dim>::ansatz_point (const unsigned int i) const {
Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
- Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
+ Assert (update_flags & update_ansatz_points,
+ ExcAccessToUninitializedField());
return ansatz_points[i];
};
template <int dim>
const Point<dim> & FEFaceValues<dim>::normal_vector (const unsigned int i) const {
Assert (i<normal_vectors.size(), ExcInvalidIndex(i, normal_vectors.size()));
- Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField());
+ Assert (update_flags & update_normal_vectors,
+ ExcAccessToUninitializedField());
return normal_vectors[i];
};
template <int dim>
double FEFaceValues<dim>::JxW (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
- Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
+ Assert (update_flags & update_JxW_values,
+ ExcAccessToUninitializedField());
return JxW_values[i];
};
const unsigned int face_no,
const FiniteElement<dim> &fe,
const Boundary<dim> &boundary) {
- cout << "ATTENTION: update_flags all set to " << (update_flags=255) << endl;
-
+ present_cell = cell;
selected_face = face_no;
// fill jacobi matrices and real
// quadrature points
Assert (update_flags & update_jacobians, ExcCannotInitializeField());
for (unsigned int i=0; i<fe.total_dofs; ++i)
- for (unsigned int j=0; j<n_quadrature_points; ++j)
- for (unsigned int s=0; s<dim; ++s)
- {
- shape_gradients[i][j](s) = 0;
-
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ shape_gradients[i][j] = Point<dim>();
+
+ for (unsigned int s=0; s<dim; ++s)
// (grad psi)_s =
// (grad_{\xi\eta})_b J_{bs}
// with J_{bs}=(d\xi_b)/(dx_s)
for (unsigned int b=0; b<dim; ++b)
shape_gradients[i][j](s)
- +=
- unit_shape_gradients[face_no][i][j](b) * jacobi_matrices[j](b,s);
- };
+ += (unit_shape_gradients[face_no][i][j](b) *
+ jacobi_matrices[j](b,s));
+ };
};
// determinant
if (update_flags & update_JxW_values)
{
- Assert (update_flags & update_jacobians, ExcCannotInitializeField());
+ Assert (update_flags & update_jacobians,
+ ExcCannotInitializeField());
for (unsigned int i=0; i<n_quadrature_points; ++i)
JxW_values[i] = weights[i] * face_jacobi_determinants[i];
};
-
- cout << "----------------------------------------" << endl;
- cout << "cell vertices: \n";
- for (unsigned int i=0; i<2*dim; ++i)
- cout << " " << cell->vertex(i) << endl;
- cout << "face no= " << face_no << endl;
-
- /*
- cout << "----------------------------------------" << endl;
- cout << "Unit face quadrature points:\n";
- for (unsigned i=0; i<unit_quadrature_points.size(); ++i)
- cout << " " << unit_quadrature_points[i];
- cout << endl;
- */
- cout << "----------------------------------------" << endl;
- cout << "Unit global quadrature points:\n";
- for (unsigned int face=0; face<2*dim; ++face)
- {
- cout << " face=" << face << endl;
- for (unsigned int i=0; i<unit_quadrature_points.size(); ++i)
- cout << global_unit_quadrature_points[face][i] << endl;
- };
- /*
- cout << "----------------------------------------" << endl;
- cout << "shape values: \n";
- for (unsigned int i=0; i<2*dim; ++i)
- {
- cout << "face=" << i << endl;
- shape_values[i].print_formatted(cout);
- };
- cout << endl;
- */
-
- cout << "----------------------------------------" << endl;
- Assert (false,ExcAccessToUninitializedField());
};