// down to the base class's
// constructor, or are stored and
// used to generate a
- // ``DoFHandler'' object later.
+ // ``DoFHandler'' object
+ // later. Since finite elements and
+ // quadrature formula should match,
+ // it is also passed a quadrature
+ // object.
//
// The ``solve_problem'' sets up
// the data structures for the
public:
Solver (Triangulation<dim> &triangulation,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &boundary_values);
virtual
~Solver ();
// examples:
protected:
const SmartPointer<const FiniteElement<dim> > fe;
+ const SmartPointer<const Quadrature<dim> > quadrature;
DoFHandler<dim> dof_handler;
Vector<double> solution;
const SmartPointer<const Function<dim> > boundary_values;
template <int dim>
Solver<dim>::Solver (Triangulation<dim> &triangulation,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &boundary_values)
:
Base<dim> (triangulation),
fe (&fe),
+ quadrature (&quadrature),
dof_handler (triangulation),
boundary_values (&boundary_values)
{};
const typename DoFHandler<dim>::active_cell_iterator &end_cell,
Threads::ThreadMutex &mutex) const
{
- //TODO: adaptive
- QGauss4<dim> quadrature_formula;
-
- FEValues<dim> fe_values (*fe, quadrature_formula,
+ FEValues<dim> fe_values (*fe, *quadrature,
UpdateFlags(update_gradients |
update_JxW_values));
const unsigned int dofs_per_cell = fe->dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ const unsigned int n_q_points = quadrature->n_quadrature_points;
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
public:
PrimalSolver (Triangulation<dim> &triangulation,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
protected:
PrimalSolver<dim>::
PrimalSolver (Triangulation<dim> &triangulation,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
Base<dim> (triangulation),
- Solver<dim> (triangulation, fe, boundary_values),
+ Solver<dim> (triangulation, fe,
+ quadrature, boundary_values),
rhs_function (&rhs_function)
{};
PrimalSolver<dim>::
assemble_rhs (Vector<double> &rhs) const
{
- //TODO: adaptive
- QGauss4<dim> quadrature_formula;
-
- FEValues<dim> fe_values (*fe, quadrature_formula,
+ FEValues<dim> fe_values (*fe, *quadrature,
UpdateFlags(update_values |
update_q_points |
update_JxW_values));
const unsigned int dofs_per_cell = fe->dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ const unsigned int n_q_points = quadrature->n_quadrature_points;
Vector<double> cell_rhs (dofs_per_cell);
std::vector<double> rhs_values (n_q_points);
public:
RefinementGlobal (Triangulation<dim> &coarse_grid,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
RefinementGlobal<dim>::
RefinementGlobal (Triangulation<dim> &coarse_grid,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
Base<dim> (coarse_grid),
- PrimalSolver<dim> (coarse_grid, fe,
+ PrimalSolver<dim> (coarse_grid, fe, quadrature,
rhs_function, boundary_values)
{};
public:
RefinementKelly (Triangulation<dim> &coarse_grid,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
RefinementKelly<dim>::
RefinementKelly (Triangulation<dim> &coarse_grid,
const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
Base<dim> (coarse_grid),
- PrimalSolver<dim> (coarse_grid, fe, rhs_function, boundary_values)
+ PrimalSolver<dim> (coarse_grid, fe, quadrature,
+ rhs_function, boundary_values)
{};
Triangulation<dim> triangulation;
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (2);
- const FE_Q<dim> fe(1);
+ const FE_Q<dim> fe(1);
+ const QGauss4<dim> quadrature;
const RightHandSide<dim> rhs_function;
const Solution<dim> boundary_values;
LaplaceSolver::Base<dim> * solver = 0;
if (solver_name == "global")
solver = new LaplaceSolver::RefinementGlobal<dim> (triangulation, fe,
+ quadrature,
rhs_function,
boundary_values);
else if (solver_name == "kelly")
solver = new LaplaceSolver::RefinementKelly<dim> (triangulation, fe,
+ quadrature,
rhs_function,
boundary_values);
else