// <code>setup_system</code>. Apart from this,
// everything is as before.
template <int dim>
-class LaplaceProblem
+class LaplaceProblem
{
public:
LaplaceProblem ();
void run ();
-
+
private:
void setup_system ();
void assemble_system ();
// here just as in the previous
// example.
template <int dim>
-class Coefficient : public Function<dim>
+class Coefficient : public Function<dim>
{
public:
Coefficient () : Function<dim>() {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component = 0) const;
// after all):
template <int dim>
double Coefficient<dim>::value (const Point<dim> &p,
- const unsigned int /*component*/) const
+ const unsigned int /*component*/) const
{
if (p.square() < 0.5*0.5)
return 20;
template <int dim>
void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
- const unsigned int component) const
+ const unsigned int component) const
{
- Assert (values.size() == points.size(),
+ Assert (values.size() == points.size(),
ExcDimensionMismatch (values.size(), points.size()));
// Since examples are not very good
// if they do not demonstrate their
// call stack to immediately find
// the place where the the array
// with the wrong size was set up.
-
+
// While we're at it, we can do
// another check: the coefficient
// is a scalar, but the
// of writing <code>v.size()-1</code> in
// many places, the range is
// defined as half-open.)
- Assert (component == 0,
+ Assert (component == 0,
ExcIndexRange (component, 0, 1));
// The rest of the function is
// are completely unchanged from
// before:
template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
+void LaplaceProblem<dim>::assemble_system ()
+{
QGauss<dim> quadrature_formula(2);
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
}
- cell->distribute_local_to_global (cell_rhs, system_rhs);
- cell->distribute_local_to_global (cell_matrix, system_matrix);
- /*
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
-
+
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
- */
}
// With the matrix so built, we use
// declared, and the CG solver will
// do the rest for us:
template <int dim>
-void LaplaceProblem<dim>::solve ()
+void LaplaceProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
// this is the first cycle, however,
// we first have to generate a mesh:
template <int dim>
-void LaplaceProblem<dim>::run ()
+void LaplaceProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<6; ++cycle)
{
// might be worth to show
// what not to do, after
// all.
-
+
// So if we got past the
// assertion, we know that
// dim==2, and we can now
// like the one in the previous
// example, so we won't comment on it
// further:
-int main ()
+int main ()
{
deallog.depth_console (0);
// Results section of the program
// to see what happens when the
// code is actually run:
-/*
+/*
Coefficient<2> coefficient;
std::vector<Point<2> > points (2);
std::vector<double> coefficient_values (1);
coefficient.value_list (points, coefficient_values);
*/
-
+
return 0;
}