// tria.create_hyper_ball(Point<2>(2,3),4);
// tria.set_boundary (&boundary);
-/* tria.refine_global (1);
- (--tria.last_active())->set_refine_flag();
- tria.execute_refinement ();
- tria.begin_active(2)->set_refine_flag();
- tria.execute_refinement ();
+ tria.refine_global (1);
+// (--tria.last_active())->set_refine_flag();
+// tria.execute_refinement ();
+// tria.begin_active(2)->set_refine_flag();
+// tria.execute_refinement ();
tria.refine_global (2);
-*/
+/*
const unsigned int dim=2;
tria.refine_global (1);
tria.execute_refinement ();
};
tria.refine_global (1);
-
+*/
cout << "Distributing dofs... ";
dof.distribute_dofs (fe);
$dx dy$ from the real to the unit cell, we have to take the determinant of
the inverse matrix, which is the reciprocal value of the determinant of the
matrix defined above.
+
+ The #FEValues# object keeps track of those fields which really need to
+ be computed, since the computation of the gradients of the ansatz functions
+ on each real cell can be quite an expensive thing if it is not needed. The
+ object knows about which fields are needed by the #UpdateStruct# object
+ 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.
*/
template <int dim>
class FEValues {
public:
+ /**
+ * Provide a structure which tells the
+ * #reinit# function, which fields are
+ * to be updated for each cell. E.g. if
+ * you do not need the gradients since
+ * you want to assemble the mass matrix,
+ * you can switch that off. By default,
+ * all flags are off, i.e. no
+ * reinitialization will be done.
+ *
+ * A structure of this type has to be
+ * passed to the constructor of the
+ * #FEValues# object.
+ */
+ struct UpdateStruct {
+ /**
+ * Constructor. Sets all fields to
+ * false.
+ */
+ UpdateStruct ();
+ /**
+ * Compute quadrature points in real
+ * space (not on unit cell).
+ */
+ bool update_q_points;
+ /**
+ * Transform gradients on unit cell to
+ * gradients on real cell.
+ */
+ bool update_gradients;
+ /**
+ * Compute jacobian matrices of the
+ * transform between unit and real cell
+ * in the evaluation points.
+ */
+ bool update_jacobians;
+ /**
+ * Compute the JxW values (Jacobian
+ * determinant at the quadrature point
+ * times the weight of this point).
+ */
+ bool update_JxW_values;
+ /**
+ * Compute the points on the real cell
+ * on which the ansatz functions are
+ * located.
+ */
+ bool update_ansatz_points;
+ };
+
+
+
/**
* Number of quadrature points.
*/
* quadrature rule.
*/
FEValues (const FiniteElement<dim> &,
- const Quadrature<dim> &);
+ const Quadrature<dim> &,
+ const UpdateStruct &);
/**
* Return the value of the #i#th shape
int, int,
<< "The index " << arg1
<< " is out of range, it should be less than " << arg2);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcAccessToUninitializedField);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcCannotInitializeField);
private:
/**
* is set each time #reinit# is called.
*/
vector<dFMatrix> jacobi_matrices;
+
+ /**
+ * Store which fields are to be updated by
+ * the reinit function.
+ */
+ UpdateStruct update_flags;
};
* elements need different transformations
* of the unit cell to a real cell.
*
+ * The computation of the two fields may
+ * share some common code, which is why we
+ * put it in one function. However, it may
+ * not always be necessary to really
+ * compute both fields, so there are two
+ * bool flags which tell the function which
+ * of the fields to actually compute.
+ *
* Refer to the documentation of the
* \Ref{FEValues} class for a definition
* of the Jacobi matrix.
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<dim> > &points) const;
+ const bool compute_jacobians,
+ vector<Point<dim> > &points,
+ const bool compute_points) const;
/**
* Comparison operator. We also check for
virtual void fill_fe_values (const Triangulation<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<1> > &points) const;
+ const bool compute_jacobians,
+ vector<Point<1> > &points,
+ const bool compute_points) const;
};
virtual void fill_fe_values (const Triangulation<2>::cell_iterator &cell,
const vector<Point<2> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<2> > &points) const;
+ const bool compute_jacobians,
+ vector<Point<2> > &points,
+ const bool compute_points) const;
};
inline
const vector<vector<Point<dim> > > &
FEValues<dim>::get_shape_grads () const {
+ Assert (update_flags.update_gradients, ExcAccessToUninitializedField());
return shape_gradients;
};
inline
const vector<Point<dim> > &
FEValues<dim>::get_quadrature_points () const {
+ Assert (update_flags.update_q_points, ExcAccessToUninitializedField());
return quadrature_points;
};
inline
const vector<double> &
FEValues<dim>::get_JxW_values () const {
+ Assert (update_flags.update_JxW_values, ExcAccessToUninitializedField());
return JxW_values;
};
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<dim> > &points) const;
+ const bool compute_jacobians,
+ vector<Point<dim> > &points,
+ const bool compute_points) const;
};
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<dim> > &points) const;
+ const bool compute_jacobians,
+ vector<Point<dim> > &points,
+ const bool compute_points) const;
};
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<dim> > &points) const;
+ const bool compute_jacobians,
+ vector<Point<dim> > &points,
+ const bool compute_points) const;
};
const bool assemble_rhs,
ProblemBase<dim> &problem,
const Quadrature<dim> &quadrature,
- const FiniteElement<dim> &fe);
+ const FiniteElement<dim> &fe,
+ const FEValues<dim>::UpdateStruct &update_flags);
/**
* Pointer to the dof handler object
* object.
*/
const FiniteElement<dim> &fe;
+ /**
+ * Store which of the fields of the
+ * FEValues object need to be reinitialized
+ * on each cell.
+ */
+ const FEValues<dim>::UpdateStruct update_flags;
};
* For what exactly happens here, refer to
* the general doc of this class.
*/
- virtual void assemble (const Equation<dim> &equation,
- const Quadrature<dim> &q,
- const FiniteElement<dim> &fe,
- const DirichletBC &dirichlet_bc = DirichletBC());
+ virtual void assemble (const Equation<dim> &equation,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const FEValues<dim>::UpdateStruct &update_flags,
+ const DirichletBC &dirichlet_bc = DirichletBC());
/**
* Solve the system of equations.
+template <int dim>
+FEValues<dim>::UpdateStruct::UpdateStruct () :
+ update_q_points(false),
+ update_gradients(false),
+ update_jacobians(false),
+ update_JxW_values(false),
+ update_ansatz_points(false) {};
+
+
+
+
+
/*------------------------------- FEValues -------------------------------*/
template <int dim>
FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature) :
+ const Quadrature<dim> &quadrature,
+ const UpdateStruct &update_flags) :
n_quadrature_points(quadrature.n_quadrature_points),
total_dofs(fe.total_dofs),
shape_values(fe.total_dofs, quadrature.n_quadrature_points),
quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
jacobi_matrices (quadrature.n_quadrature_points,
- dFMatrix(dim,dim))
+ dFMatrix(dim,dim)),
+ update_flags (update_flags)
{
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
const unsigned int j) const {
Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+ Assert (update_flags.update_gradients, ExcAccessToUninitializedField());
return shape_gradients[i][j];
};
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));
+ Assert (update_flags.update_q_points, ExcAccessToUninitializedField());
return quadrature_points[i];
};
template <int dim>
double FEValues<dim>::JxW (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+ Assert (update_flags.update_JxW_values, ExcAccessToUninitializedField());
return JxW_values[i];
};
void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
const FiniteElement<dim> &fe) {
// fill jacobi matrices and real
- //quadrature points
- fe.fill_fe_values (cell,
- unit_quadrature_points,
- jacobi_matrices,
- quadrature_points);
-
- // compute gradients on real element
- 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;
-
- // (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[i][j](b) * jacobi_matrices[j](b,s);
- };
+ // quadrature points
+ if (update_flags.update_jacobians || update_flags.update_q_points)
+ fe.fill_fe_values (cell,
+ unit_quadrature_points,
+ jacobi_matrices,
+ update_flags.update_jacobians,
+ quadrature_points,
+ update_flags.update_q_points);
+
+ // compute gradients on real element if
+ // requested
+ if (update_flags.update_gradients)
+ {
+ 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;
+
+ // (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[i][j](b) * jacobi_matrices[j](b,s);
+ };
+ };
+
// compute Jacobi determinants in
// quadrature points.
// refer to the general doc for
// why we take the inverse of the
// determinant
- for (unsigned int i=0; i<n_quadrature_points; ++i)
- JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+ if (update_flags.update_JxW_values)
+ {
+ Assert (update_flags.update_jacobians,
+ ExcCannotInitializeField());
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+ };
};
void FiniteElementBase<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator &,
const vector<Point<dim> > &,
vector<dFMatrix> &,
- vector<Point<dim> > &) const {
+ const bool,
+ vector<Point<dim> > &,
+ const bool) const {
Assert (false, ExcPureFunctionCalled());
};
void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<1> > &points) const {
+ const bool compute_jacobians,
+ vector<Point<1> > &points,
+ const bool compute_points) const {
// local mesh width
const double h=(cell->vertex(1)(0) - cell->vertex(0)(0));
for (unsigned int i=0; i<points.size(); ++i)
{
- jacobians[i](0,0) = 1./h;
- points[i] = cell->vertex(0) + h*unit_points[i];
+ if (compute_jacobians)
+ jacobians[i](0,0) = 1./h;
+ if (compute_points)
+ points[i] = cell->vertex(0) + h*unit_points[i];
};
};
void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
const vector<Point<2> > &,
vector<dFMatrix> &,
- vector<Point<2> > &) const {
+ const bool,
+ vector<Point<2> > &,
+ const bool) const {
Assert (false, ExcPureFunctionCalled());
};
/*------------------------------- Explicit Instantiations -------------*/
+template struct FEValues<1>::UpdateStruct;
+template struct FEValues<2>::UpdateStruct;
+
template class FEValues<1>;
template class FEValues<2>;
void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<1> > &points) const {
+ const bool compute_jacobians,
+ vector<Point<1> > &points,
+ const bool compute_points) const {
// simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+ FiniteElement<1>::fill_fe_values (cell, unit_points,
+ jacobians, compute_jacobians,
+ points, compute_points);
};
void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
const vector<Point<2> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<2> > &points) const {
+ const bool compute_jacobians,
+ vector<Point<2> > &points,
+ const bool compute_points) const {
const unsigned int dim=2;
const unsigned int n_vertices=4;
unsigned int n_points=unit_points.size();
- // recomputation of values
Point<dim> vertices[n_vertices];
for (unsigned int l=0; l<n_vertices; ++l)
vertices[l] = cell->vertex(l);
+
+ if (compute_points)
+ {
+ // initialize points to zero
+ for (unsigned int i=0; i<n_points; ++i)
+ points[i] = Point<dim> ();
+
+ // note: let x_l be the vector of the
+ // lth quadrature point in real space and
+ // xi_l that on the unit cell, let further
+ // p_j be the vector of the jth vertex
+ // of the cell in real space and
+ // N_j(xi_l) be the value of the associated
+ // basis function at xi_l, then
+ // x_l(xi_l) = sum_j p_j N_j(xi_l)
+ for (unsigned int j=0; j<n_vertices; ++j)
+ for (unsigned int l=0; l<n_points; ++l)
+ points[l] += vertices[j] * shape_value(j, unit_points[l]);
+ };
- // initialize points to zero
- for (unsigned int i=0; i<n_points; ++i)
- points[i] = Point<dim> ();
-
- // note: let x_l be the vector of the
- // lth quadrature point in real space and
- // xi_l that on the unit cell, let further
- // p_j be the vector of the jth vertex
- // of the cell in real space and
- // N_j(xi_l) be the value of the associated
- // basis function at xi_l, then
- // x_l(xi_l) = sum_j p_j N_j(xi_l)
- for (unsigned int j=0; j<n_vertices; ++j)
- for (unsigned int l=0; l<n_points; ++l)
- points[l] += vertices[j] * shape_value(j, unit_points[l]);
/* jacobi matrices: compute d(x)/d(xi) and invert this
Let M(l) be the inverse of J at the quadrature point l, then
However, we rewrite the loops to only compute the gradient once for
each integration point and basis function.
*/
- dFMatrix M(dim,dim);
- for (unsigned int l=0; l<n_points; ++l)
+ if (compute_jacobians)
{
- M.clear ();
- for (unsigned int s=0; s<n_vertices; ++s)
+ dFMatrix M(dim,dim);
+ for (unsigned int l=0; l<n_points; ++l)
{
- // we want a linear transform and
- // if we prepend the class name in
- // front of the #shape_grad#, we
- // need not use virtual function
- // calls.
- const Point<dim> gradient
- = FELinear<dim>::shape_grad (s, unit_points[l]);
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- M(i,j) += vertices[s](i) * gradient(j);
+ M.clear ();
+ for (unsigned int s=0; s<n_vertices; ++s)
+ {
+ // we want a linear transform and
+ // if we prepend the class name in
+ // front of the #shape_grad#, we
+ // need not use virtual function
+ // calls.
+ const Point<dim> gradient
+ = FELinear<dim>::shape_grad (s, unit_points[l]);
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ M(i,j) += vertices[s](i) * gradient(j);
+ };
+ jacobians[l].invert(M);
};
- jacobians[l].invert(M);
};
};
void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<1> > &points) const {
+ const bool compute_jacobians,
+ vector<Point<1> > &points,
+ const bool compute_points) const {
// simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+ FiniteElement<1>::fill_fe_values (cell, unit_points,
+ jacobians, compute_jacobians,
+ points, compute_points);
};
void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
const vector<Point<2> > &,
vector<dFMatrix> &,
- vector<Point<2> > &) const {
+ const bool,
+ vector<Point<2> > &,
+ const bool) const {
Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
vector<dFMatrix> &jacobians,
- vector<Point<1> > &points) const {
+ const bool compute_jacobians,
+ vector<Point<1> > &points,
+ const bool compute_points) const {
// simply pass down
- FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+ FiniteElement<1>::fill_fe_values (cell, unit_points,
+ jacobians, compute_jacobians,
+ points, compute_points);
};
void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
const vector<Point<2> > &,
vector<dFMatrix> &,
- vector<Point<2> > &) const {
+ const bool,
+ vector<Point<2> > &,
+ const bool) const {
Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
const bool assemble_rhs,
ProblemBase<dim> &problem,
const Quadrature<dim> &quadrature,
- const FiniteElement<dim> &fe) :
+ const FiniteElement<dim> &fe,
+ const FEValues<dim>::UpdateStruct &update_flags) :
dof(dof),
assemble_matrix(assemble_matrix),
assemble_rhs(assemble_rhs),
problem(problem),
quadrature(quadrature),
- fe(fe) {};
+ fe(fe),
+ update_flags(update_flags) {};
problem (((AssemblerData<dim>*)local_data)->problem),
fe(((AssemblerData<dim>*)local_data)->fe),
fe_values (((AssemblerData<dim>*)local_data)->fe,
- ((AssemblerData<dim>*)local_data)->quadrature)
+ ((AssemblerData<dim>*)local_data)->quadrature,
+ ((AssemblerData<dim>*)local_data)->update_flags)
{
Assert ((unsigned int)problem.system_matrix.m() == dof_handler->n_dofs(),
ExcInvalidData());
template <int dim>
-void ProblemBase<dim>::assemble (const Equation<dim> &equation,
- const Quadrature<dim> &quadrature,
- const FiniteElement<dim> &fe,
- const DirichletBC &dirichlet_bc) {
+void ProblemBase<dim>::assemble (const Equation<dim> &equation,
+ const Quadrature<dim> &quadrature,
+ const FiniteElement<dim> &fe,
+ const FEValues<dim>::UpdateStruct &update_flags,
+ const DirichletBC &dirichlet_bc) {
system_sparsity.reinit (dof_handler->n_dofs(),
dof_handler->n_dofs(),
dof_handler->max_couplings_between_dofs());
// reinit solution vector, preset
// with zeroes.
solution.reinit (dof_handler->n_dofs());
-
+
// create assembler
AssemblerData<dim> data (*dof_handler,
true, true, //assemble matrix and rhs
*this,
quadrature,
- fe);
+ fe,
+ update_flags);
TriaActiveIterator<dim, Assembler<dim> > assembler (tria,
tria->begin_active()->level(),
tria->begin_active()->index(),
&data);
-
// loop over all cells, fill matrix and rhs
do
{
difference.erase (difference.begin(), difference.end());
difference.reserve (tria->n_cells());
- FEValues<dim> fe_values(fe, q);
+ FEValues<dim>::UpdateStruct update_flags;
+ update_flags.update_q_points = true;
+ update_flags.update_JxW_values = true;
+ FEValues<dim> fe_values(fe, q, update_flags);
// loop over all cells
// (we need an iterator on the triangulation for
psi);
// then subtract finite element
// solution
- for (unsigned int j=0; j<n_dofs; ++j)
+ const vector<double> &JxW_values = fe_values.get_JxW_values();
+ for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int i=0; i<n_dofs; ++i)
- psi[j] -= dof_values[i]*shape_values(i,j);
+ psi[j] -= dof_values[i]*shape_values(i,j)*JxW_values[j];
// for L1_norm and Linfty_norm:
// take absolute
// tria.create_hyper_ball(Point<2>(2,3),4);
// tria.set_boundary (&boundary);
-/* tria.refine_global (1);
- (--tria.last_active())->set_refine_flag();
- tria.execute_refinement ();
- tria.begin_active(2)->set_refine_flag();
- tria.execute_refinement ();
+ tria.refine_global (1);
+// (--tria.last_active())->set_refine_flag();
+// tria.execute_refinement ();
+// tria.begin_active(2)->set_refine_flag();
+// tria.execute_refinement ();
tria.refine_global (2);
-*/
+/*
const unsigned int dim=2;
tria.refine_global (1);
tria.execute_refinement ();
};
tria.refine_global (1);
-
+*/
cout << "Distributing dofs... ";
dof.distribute_dofs (fe);