/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-#include <grid/tria.h>
-#include <dofs/mg_dof_handler.h>
-#include <grid/tria_accessor.h>
-#include <dofs/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-#include <grid/grid_generator.h>
-#include <base/function.h>
-#include <numerics/data_out.h>
-#include <fe/fe_lib.lagrange.h>
-#include <fe/fe_lib.criss_cross.h>
-#include <fe/fe_update_flags.h>
-#include <base/quadrature_lib.h>
-#include <numerics/assembler.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
#include <lac/vector.h>
+#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/vector_memory.h>
#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/fe_lib.lagrange.h>
+#include <grid/grid_out.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
-#include <map>
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <cstdlib>
-
-
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation (const Function<dim> &rhs) :
- Equation<dim>(1),
- right_hand_side (rhs) {};
-
- virtual void assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (FullMatrix<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- protected:
- const Function<dim> &right_hand_side;
-};
-
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_dof_tools.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_smoother.h>
+#include <multigrid/multigrid.h>
+#include <lac/solver_richardson.h>
+#include <fstream>
template <int dim>
-class PoissonProblem {
+class LaplaceProblem
+{
public:
- /**
- * Declare a data type which denotes a
- * mapping between a boundary indicator
- * and the function denoting the boundary
- * values on this part of the boundary.
- * Only one boundary function may be given
- * for each boundary indicator, which is
- * guaranteed by the #map# data type.
- *
- * See the general documentation of this
- * class for more detail.
- */
- typedef map<unsigned char,const Function<dim>*> FunctionMap;
- /**
- * Typdedef an iterator which assembles
- * matrices and vectors.
- */
- typedef TriaActiveIterator<dim, Assembler<dim> > active_assemble_iterator;
-
- PoissonProblem (unsigned int order);
- virtual ~PoissonProblem ();
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+ void run ();
- void clear ();
- void create_new ();
+ private:
+ void setup_system ();
+ void assemble_system ();
void solve ();
-
- /**
- * Initiate the process of assemblage of
- * vectors and system matrix. Use the
- * given equation object and the given
- * quadrature formula. Also use the list
- * of dirichlet boundary value functions
- * (by default, no dirichlet bc are assumed
- * which means that all bc are included
- * into the weak formulation).
- *
- * 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 UpdateFlags update_flags,
- const FunctionMap &dirichlet_bc = FunctionMap());
-
- int run (unsigned int level);
- void print_history (string filename) const;
-
- /**
- * Exception
- */
- DeclException0 (ExcNoTriaSelected);
-
- protected:
- Triangulation<dim> *tria;
- MGDoFHandler<dim> *dof;
-
-
- /**
- * Sparsity pattern of the system matrix.
- */
- SparsityPattern system_sparsity;
-
- /**
- * System matrix.
- */
- SparseMatrix<double> system_matrix;
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
- /**
- * Vector storing the right hand side.
- */
- Vector<double> right_hand_side;
+ Triangulation<dim> triangulation;
+ MGDoFHandler<dim> mg_dof_handler;
- /**
- * Solution vector.
- */
- Vector<double> solution;
+ FEQ1<dim> fe;
- /**
- * List of constraints introduced by
- * hanging nodes.
- */
- ConstraintMatrix constraints;
+ ConstraintMatrix hanging_node_constraints;
- Function<dim> *rhs;
- Function<dim> *boundary_values;
+ SparsityPattern global_sparsity_pattern;
+ SparseMatrix<double> global_system_matrix;
- vector<double> l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error;
- vector<int> n_dofs;
-
- unsigned int order;
+ MGLevelObject<SparsityPattern> level_sparsity_patterns;
+ MGLevelObject<SparseMatrix<double> > level_system_matrices;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
};
-
-
-/**
- Right hand side constructed such that the exact solution is
- $sin(2 pi x) + sin(2 pi y)$
- */
template <int dim>
-class RHSPoly : public Function<dim> {
+class Coefficient : public Function<dim>
+{
public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int component = 0) const;
};
template <int dim>
-class Solution : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- /**
- * Return the gradient of the function
- * at the given point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component) const;
+double Coefficient<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+{
+ if (p.square() < 0.5*0.5)
+ return 20;
+ else
+ return 1;
};
+template <int dim>
+void Coefficient<dim>::value_list (const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int component) const
+{
+ const unsigned int n_points = points.size();
-template <>
-double RHSPoly<2>::value (const Point<2> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
+ Assert (values.size() == n_points,
+ ExcVectorHasWrongSize (values.size(), n_points));
+
+ Assert (component == 0,
+ ExcWrongComponent (component, 1));
- const double x = p(0),
- y = p(1);
- const double pi= 3.1415926536;
- return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y));
+ for (unsigned int i=0; i<n_points; ++i)
+ {
+ if (points[i].square() < 0.5*0.5)
+ values[i] = 20;
+ else
+ values[i] = 1;
+ };
};
-template <>
-double Solution<2>::value (const Point<2> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
- const double x = p(0),
- y = p(1);
- const double pi= 3.1415926536;
- return sin(2*pi*x)+sin(2*pi*y);
+class MGSmootherLAC : public MGSmootherBase
+{
+ private:
+ SmartPointer<MGLevelObject<SparseMatrix<double> > >matrices;
+ public:
+ MGSmootherLAC(MGLevelObject<SparseMatrix<double> >&);
+
+ virtual void smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
};
-template <>
-Tensor<1,2> Solution<2>::gradient (const Point<2> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
+MGSmootherLAC::MGSmootherLAC(MGLevelObject<SparseMatrix<double> >& matrix)
+ :
+ matrices(&matrix)
+{}
- const double x = p(0),
- y = p(1);
- const double pi= 3.1415926536;
- return Point<2> (2*pi*cos(2*pi*x),
- 2*pi*cos(2*pi*y));
-};
-
+void
+MGSmootherLAC::smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const
+{
+ SolverControl control(2,1.e-300,false,false);
+ PrimitiveVectorMemory<> mem;
+ SolverRichardson<> rich(control, mem);
+ PreconditionRelaxation<>
+ prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
+ rich.solve((*matrices)[level], u, rhs, prec);
+}
-template <>
-void PoissonEquation<2>::assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<2> &fe_values,
- const DoFHandler<2>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side.value(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+ mg_dof_handler (triangulation)
+{};
template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, ExcPureVirtualFunctionCalled());
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+ mg_dof_handler.clear ();
};
template <int dim>
-void PoissonEquation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, ExcPureVirtualFunctionCalled());
-};
-
+void LaplaceProblem<dim>::setup_system ()
+{
+ mg_dof_handler.distribute_dofs (fe);
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (mg_dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+ global_sparsity_pattern.reinit (mg_dof_handler.DoFHandler<dim>::n_dofs(),
+ mg_dof_handler.DoFHandler<dim>::n_dofs(),
+ mg_dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (mg_dof_handler, global_sparsity_pattern);
+ hanging_node_constraints.condense (global_sparsity_pattern);
+ global_sparsity_pattern.compress();
+ global_system_matrix.reinit (global_sparsity_pattern);
+ solution.reinit (mg_dof_handler.DoFHandler<dim>::n_dofs());
+ system_rhs.reinit (mg_dof_handler.DoFHandler<dim>::n_dofs());
+ const unsigned int n_levels = triangulation.n_levels();
+ level_system_matrices.resize (0, n_levels);
+ level_sparsity_patterns.resize (0, n_levels);
+
+ for (unsigned int level=0; level<n_levels; ++level)
+ {
+ level_sparsity_patterns[level].reinit (mg_dof_handler.n_dofs(level),
+ mg_dof_handler.n_dofs(level),
+ mg_dof_handler.max_couplings_between_dofs()); //xxx
+ MGDoFTools::make_sparsity_pattern (mg_dof_handler,
+ level,
+ level_sparsity_patterns[level]);
+ level_sparsity_patterns[level].compress();
+
+ level_system_matrices[level].reinit (level_sparsity_patterns[level]);
+ };
+};
template <int dim>
-PoissonProblem<dim>::PoissonProblem (unsigned int order) :
- tria(0), dof(0), rhs(0),
- boundary_values(0), order(order) {};
+void LaplaceProblem<dim>::assemble_system ()
+{
+ const Coefficient<dim> coefficient;
+ QGauss2<dim> quadrature_formula;
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
-template <int dim>
-PoissonProblem<dim>::~PoissonProblem ()
-{
- clear ();
-};
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+ vector<int> local_dof_indices (dofs_per_cell);
-template <int dim>
-void PoissonProblem<dim>::clear () {
- if (dof != 0) {
- delete dof;
- dof = 0;
- };
-
- if (tria != 0) {
- delete tria;
- tria = 0;
- };
-
- if (rhs != 0)
- {
- delete rhs;
- rhs = 0;
- };
+ // FIX
+ vector<double> coefficient_values (n_q_points, 1.0);
- if (boundary_values != 0)
+ // not only active cells
+ MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+ endc = mg_dof_handler.end();
+ for (; cell!=endc; ++cell)
{
- delete boundary_values;
- boundary_values = 0;
+ cell_matrix.clear ();
+ cell_rhs.clear ();
+
+ fe_values.reinit (cell);
+ const FullMatrix<double>
+ & shape_values = fe_values.get_shape_values();
+ const vector<vector<Tensor<1,dim> > >
+ & shape_grads = fe_values.get_shape_grads();
+ const vector<double>
+ & JxW_values = fe_values.get_JxW_values();
+ const vector<Point<dim> >
+ & q_points = fe_values.get_quadrature_points();
+
+ // FIX
+// coefficient.value_list (q_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)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (coefficient_values[q_point] *
+ (shape_grads[i][q_point] *
+ shape_grads[j][q_point] +
+ shape_values(i,q_point) *
+ shape_values(j,q_point) ) *
+ JxW_values[q_point]);
+
+ cell_rhs(i) += (shape_values (i,q_point) *
+ sin(4*sqrt(q_points[q_point].square())) *
+ fe_values.JxW (q_point));
+ };
+
+
+ cell->get_mg_dof_indices (local_dof_indices);
+ const unsigned int level = cell->level();
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ level_system_matrices[level].add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ // if active, then also into
+ // global matrix
+ if (cell->active())
+ {
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ global_system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ };
+ };
};
- system_sparsity.reinit (0,0,1);
- system_matrix.clear ();
- right_hand_side.reinit (1);
- solution.reinit (1);
- constraints.clear ();
+ hanging_node_constraints.condense (global_system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+
+// map<int,double> boundary_values;
+// VectorTools::interpolate_boundary_values (mg_dof_handler,
+// 0,
+// ZeroFunction<dim>(),
+// boundary_values);
+// MatrixTools<dim>::apply_boundary_values (boundary_values,
+// global_system_matrix,
+// solution,
+// system_rhs);
};
-
template <int dim>
-void PoissonProblem<dim>::create_new () {
- clear ();
-
- tria = new Triangulation<dim>();
- dof = new MGDoFHandler<dim> (*tria);
-};
+void LaplaceProblem<dim>::solve ()
+{
+ {
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+ SolverControl coarse_grid_solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> coarse_grid_vector_memory;
+
+ SolverCG<> coarse_grid_cg (coarse_grid_solver_control,
+ coarse_grid_vector_memory);
+
+// PreconditionRelaxation<>
+// coarse_grid_solver_preconditioner(level_system_matrices[level_system_matrices.get_minlevel()],
+// &SparseMatrix<double>::template precondition_SSOR<double>,
+// 1.2);
+ PreconditionIdentity coarse_grid_solver_preconditioner;
+
+ MGCoarseGridLACIteration<SolverCG<>, SparseMatrix<double>, PreconditionIdentity>
+ coarse_grid_solver (coarse_grid_cg,
+ level_system_matrices[level_system_matrices.get_minlevel()],
+ coarse_grid_solver_preconditioner);
+
+ MGSmootherLAC smoother (level_system_matrices);
+ MGTransferPrebuilt grid_transfer;
+ grid_transfer.build_matrices (mg_dof_handler);
+
+ Multigrid<2> multigrid (mg_dof_handler,
+ hanging_node_constraints,
+ level_sparsity_patterns,
+ level_system_matrices,
+ grid_transfer);
+
+ PreconditionMG<Multigrid<2> >
+ mg_precondition (multigrid, smoother, smoother, coarse_grid_solver);
-template <int dim>
-void PoissonProblem<dim>::assemble (const Equation<dim> &equation,
- const Quadrature<dim> &quadrature,
- const UpdateFlags update_flags,
- const FunctionMap &dirichlet_bc) {
- Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected());
-
- system_sparsity.reinit (dof->DoFHandler<dim>::n_dofs(),
- dof->DoFHandler<dim>::n_dofs(),
- dof->max_couplings_between_dofs());
- right_hand_side.reinit (dof->DoFHandler<dim>::n_dofs());
-
- // make up sparsity pattern and
- // compress with constraints
- constraints.clear ();
- DoFTools::make_hanging_node_constraints (*dof, constraints);
- constraints.close ();
- DoFTools::make_sparsity_pattern (*dof, system_sparsity);
- constraints.condense (system_sparsity);
-
- MGTransferPrebuilt p;
- p.build_matrices (*dof);
-// MGSmoother smoother(*dof);
-
+ solution.clear ();
+ cg.solve (global_system_matrix, solution, system_rhs,
+ mg_precondition);
- // reinite system matrix
- system_matrix.reinit (system_sparsity);
- // reinit solution vector, preset
- // with zeroes.
- solution.reinit (dof->DoFHandler<dim>::n_dofs());
+ cout << " MG Outer iterations: " << solver_control.last_step()
+ << endl;
+
+ cout << " MG Total inner iterations: " << coarse_grid_solver_control.last_step()
+ << endl;
+ };
- // create assembler
- Assembler<dim>::AssemblerData data (*dof,
- true, true, //assemble matrix and rhs
- system_matrix,
- right_hand_side,
- quadrature,
- update_flags);
- active_assemble_iterator assembler (tria,
- tria->begin_active()->level(),
- tria->begin_active()->index(),
- &data);
- // loop over all cells, fill matrix and rhs
- do
{
- assembler->assemble (equation);
- }
- while ((++assembler).state() == valid);
-
- // condense system matrix in-place
- constraints.condense (system_matrix);
-
- // condense right hand side in-place
- constraints.condense (right_hand_side);
-
- // apply Dirichlet bc as described
- // in the docs
- map<int, double> boundary_value_list;
-
- for (FunctionMap::const_iterator bc=dirichlet_bc.begin();
- bc != dirichlet_bc.end(); ++bc)
- VectorTools::interpolate_boundary_values (*dof,
- bc->first,
- *bc->second,
- boundary_value_list);
- MatrixTools<dim>::apply_boundary_values (boundary_value_list,
- system_matrix, solution,
- right_hand_side);
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+
+ PreconditionRelaxation<>
+ preconditioner(global_system_matrix,
+ &SparseMatrix<double>::template precondition_SSOR<double>,
+ 1.2);
+
+ solution.clear ();
+ cg.solve (global_system_matrix, solution, system_rhs,
+ preconditioner);
+
+ cout << " CG Outer iterations: " << solver_control.last_step()
+ << endl;
+ };
+
+ hanging_node_constraints.distribute (solution);
+};
+
+
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (mg_dof_handler,
+ QGauss3<dim-1>(),
+ neumann_boundary,
+ solution,
+ estimated_error_per_cell);
+
+ triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+ 0.3, 0.03);
+ triangulation.execute_coarsening_and_refinement ();
};
template <int dim>
-void PoissonProblem<dim>::solve () {
- Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected());
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+ string filename = "grid-";
+ filename += ('0' + cycle);
+ Assert (cycle < 10, ExcInternalError());
- SolverControl control(4000, 1e-16);
- PrimitiveVectorMemory<Vector<double> > memory;
- SolverCG<SparseMatrix<double>,Vector<double> > cg(control,memory);
-
- // solve
- cg.solve (system_matrix, solution, right_hand_side,
- PreconditionIdentity());
- // distribute solution
- constraints.distribute (solution);
+ filename += ".eps";
+ ofstream output (filename.c_str());
+
+ GridOut grid_out;
+ grid_out.write_eps (triangulation, output);
};
template <int dim>
-int PoissonProblem<dim>::run (const unsigned int level) {
- create_new ();
-
- cout << "Refinement level = " << level
- << ", using elements of type <";
- switch (order)
+void LaplaceProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<8; ++cycle)
{
- case 0:
- cout << "criss-cross";
- break;
- default:
- cout << "Lagrange-" << order;
- break;
- };
- cout << ">" << endl;
-
- cout << " Making grid... ";
- GridGenerator::hyper_cube (*tria);
- tria->refine_global (level+1);
- tria->begin_active()->set_refine_flag();
- (++(++(tria->begin_active())))->set_refine_flag();
- tria->execute_coarsening_and_refinement ();
- cout << tria->n_active_cells() << " active cells." << endl;
-
- rhs = new RHSPoly<dim>();
- boundary_values = new Solution<dim> ();
-
+ cout << "Cycle " << cycle << ':' << endl;
- FiniteElement<dim> *fe;
- PoissonEquation<dim> equation (*rhs);
- Quadrature<dim> *quadrature;
- Quadrature<dim-1> *boundary_quadrature;
- switch (order) {
- case 0:
- fe = new FECrissCross<dim>();
- quadrature = new QCrissCross1<dim>();
- boundary_quadrature = new QGauss2<dim-1>();
- break;
- case 1:
- fe = new FEQ1<dim>();
- quadrature = new QGauss3<dim>();
- boundary_quadrature = new QGauss2<dim-1>();
- break;
- case 2:
- fe = new FEQ2<dim>();
- quadrature = new QGauss4<dim>();
- boundary_quadrature = new QGauss3<dim-1>();
- break;
- case 3:
- fe = new FEQ3<dim>();
- quadrature = new QGauss5<dim>();
- boundary_quadrature = new QGauss4<dim-1>();
- break;
- case 4:
- fe = new FEQ4<dim>();
- quadrature = new QGauss6<dim>();
- boundary_quadrature = new QGauss5<dim-1>();
- break;
- default:
- return 100000;
- };
-
- cout << " Distributing dofs... ";
- dof->distribute_dofs (*fe);
- cout << dof->DoFHandler<dim>::n_dofs() << " degrees of freedom." << endl;
- n_dofs.push_back (dof->DoFHandler<dim>::n_dofs());
-
- cout << " Assembling matrices..." << endl;
- UpdateFlags update_flags = UpdateFlags(update_values | update_q_points |
- update_gradients | update_JxW_values);
-
- map<unsigned char,const Function<dim>*> dirichlet_bc;
- dirichlet_bc[0] = boundary_values;
- assemble (equation, *quadrature, update_flags, dirichlet_bc);
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (1);
+ }
+ else
+ {
+ refine_grid ();
+ };
+
- cout << " Solving..." << endl;
- solve ();
+ cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << endl;
- Solution<dim> sol;
- Vector<float> l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell;
- Vector<float> h1_seminorm_error_per_cell, h1_error_per_cell;
-
- cout << " Calculating L1 error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- l1_error_per_cell,
- *quadrature, L1_norm);
- cout << l1_error_per_cell.l1_norm() << endl;
- l1_error.push_back (l1_error_per_cell.l1_norm());
-
- cout << " Calculating L2 error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- l2_error_per_cell,
- *quadrature, L2_norm);
- cout << l2_error_per_cell.l2_norm() << endl;
- l2_error.push_back (l2_error_per_cell.l2_norm());
-
- cout << " Calculating L-infinity error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- linfty_error_per_cell,
- *quadrature, Linfty_norm);
- cout << linfty_error_per_cell.linfty_norm() << endl;
- linfty_error.push_back (linfty_error_per_cell.linfty_norm());
-
- cout << " Calculating H1-seminorm error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- h1_seminorm_error_per_cell,
- *quadrature, H1_seminorm);
- cout << h1_seminorm_error_per_cell.l2_norm() << endl;
- h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm());
-
- cout << " Calculating H1 error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- h1_error_per_cell,
- *quadrature, H1_norm);
- cout << h1_error_per_cell.l2_norm() << endl;
- h1_error.push_back (h1_error_per_cell.l2_norm());
-
- if (dof->DoFHandler<dim>::n_dofs()<=5000)
- {
- Vector<double> l1_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> l2_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> linfty_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> h1_seminorm_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> h1_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- DoFTools::distribute_cell_to_dof_vector (*dof, l1_error_per_cell, l1_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, l2_error_per_cell, l2_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, linfty_error_per_cell,
- linfty_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, h1_seminorm_error_per_cell,
- h1_seminorm_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, h1_error_per_cell, h1_error_per_dof);
-
-// Vector<double> projected_solution;
-// ConstraintMatrix constraints;
-// constraints.close ();
-// VectorTools::project (*dof, constraints, *fe,
-// StraightBoundary<dim>(), *quadrature,
-// sol, projected_solution, false,
-// *boundary_quadrature);
-// cout << " Calculating L2 error of projected solution... ";
-// VectorTools::integrate_difference (*dof,
-// projected_solution, sol,
-// l2_error_per_cell,
-// *quadrature, *fe, L2_norm);
-// cout << l2_error_per_cell.l2_norm() << endl;
-
-
- string filename;
- filename = ('0'+order);
- filename += ".";
- filename += ('0'+level);
- filename += ".ucd";
- cout << " Writing error plots to <" << filename << ">..." << endl;
+ setup_system ();
+
+ cout << " Number of degrees of freedom: "
+ << mg_dof_handler.DoFHandler<dim>::n_dofs()
+ << endl;
- DataOut<dim> out;
- ofstream o(filename.c_str());
- out.attach_dof_handler (*dof);
-
- out.add_data_vector (solution, "u");
- out.add_data_vector (l1_error_per_dof, "L1_Error");
- out.add_data_vector (l2_error_per_dof, "L2_Error");
- out.add_data_vector (linfty_error_per_dof, "Linfty_Error");
- out.add_data_vector (h1_seminorm_error_per_dof, "H1_seminorm_Error");
- out.add_data_vector (h1_error_per_dof, "H1_Error");
- out.write_ucd (o);
- o.close ();
- }
- else
- cout << " Not writing error as grid." << endl;
-
- cout << endl;
-
- const unsigned int n_dofs = dof->DoFHandler<dim>::n_dofs();
- dof->clear ();
- tria->set_boundary (0);
- delete fe;
- delete quadrature;
- delete boundary_quadrature;
-
- return n_dofs;
-};
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ DataOut<dim>::EpsFlags eps_flags;
+ eps_flags.z_scaling = 4;
+
+ DataOut<dim> data_out;
+ data_out.set_flags (eps_flags);
-template <int dim>
-void PoissonProblem<dim>::print_history (string filename) const {
- ofstream out(filename.c_str());
- out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error"
- << endl;
- for (unsigned int i=0; i<n_dofs.size(); ++i)
- out << n_dofs[i]
- << " "
- << l1_error[i] << " "
- << l2_error[i] << " "
- << linfty_error[i] << " "
- << h1_seminorm_error[i] << " "
- << h1_error[i] << endl;
-
- double average_l1=0,
- average_l2=0,
- average_linfty=0,
- average_h1_semi=0,
- average_h1=0;
- for (unsigned int i=1; i<n_dofs.size(); ++i)
- {
- average_l1 += l1_error[i]/l1_error[i-1];
- average_l2 += l2_error[i]/l2_error[i-1];
- average_linfty += linfty_error[i]/linfty_error[i-1];
- average_h1_semi += h1_seminorm_error[i]/h1_seminorm_error[i-1];
- average_h1 += h1_error[i]/h1_error[i-1];
+ data_out.attach_dof_handler (mg_dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ ofstream output ("final-solution.eps");
+ data_out.write_eps (output);
};
-
- average_l1 /= (l1_error.size()-1);
- average_l2 /= (l1_error.size()-1);
- average_linfty /= (l1_error.size()-1);
- average_h1_semi /= (l1_error.size()-1);
- average_h1 /= (l1_error.size()-1);
-
- cout << "==========================================================\n";
- cout << "Average error reduction rates for h->h/2:" << endl;
- cout << " L1 error : " << 1./average_l1 << endl
- << " L2 error : " << 1./average_l2 << endl
- << " Linfty error : " << 1./average_linfty << endl
- << " H1 seminorm error: " << 1./average_h1_semi << endl
- << " H1 error : " << 1./average_h1 << endl;
- cout << "==========================================================\n";
- cout << "==========================================================\n";
};
+
+int main ()
+{
+ try
+ {
+ deallog.depth_console (0);
-
-int main () {
- for (unsigned int order=0; order<5; ++order)
+ LaplaceProblem<2> laplace_problem_2d;
+ laplace_problem_2d.run ();
+ }
+ catch (exception &exc)
{
- PoissonProblem<2> problem (order);
-
- unsigned int level=0;
- unsigned int n_dofs;
- do
- n_dofs = problem.run (level++);
- while (n_dofs<25000);
-
- string filename;
- switch (order)
- {
- case 0:
- filename = "criss_cross";
- break;
- case 1:
- filename = "linear";
- break;
- case 2:
- filename = "quadratic";
- break;
- case 3:
- filename = "cubic";
- break;
- case 4:
- filename = "quartic";
- break;
- };
- filename += ".history";
-
- cout << endl << "Printing convergence history to <"
- << filename << ">..." << endl;
- problem.print_history (filename);
- cout << endl << endl << endl;
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Exception on processing: " << endl
+ << exc.what() << endl
+ << "Aborting!" << endl
+ << "----------------------------------------------------"
+ << endl;
+ return 1;
+ }
+ catch (...)
+ {
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Unknown exception!" << endl
+ << "Aborting!" << endl
+ << "----------------------------------------------------"
+ << endl;
+ return 1;
};
-
+
return 0;
};
+++ /dev/null
-/*---------------------------- mg_dof_accessor.h ---------------------------*/
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-#ifndef __mg_dof_accessor_H
-#define __mg_dof_accessor_H
-/*---------------------------- mg_dof_accessor.h ---------------------------*/
-
-
-#include <dofs/dof_accessor.h>
-
-
-
-
-/**
- * Define the basis for accessors to the degrees of freedom for
- * a multigrid DoF handler object.
- *
- * Note that it is allowed to construct an object of which the
- * #mg_dof_handler# pointer is a Null pointer. Such an object would
- * result in a strange kind of behaviour, though every reasonable
- * operating system should disallow access through that pointer.
- * The reason we do not check for the null pointer in the
- * constructor which gets passed the #DoFHandler# pointer is that
- * if we did we could not make dof iterators member of other classes
- * (like in the #FEValues# class) if we did not know about the
- * #DoFHandler# object to be used upon construction of that object.
- * Through the way this class is implemented here, we allow the
- * creation of a kind of virgin object which only gets useful if
- * assigned to from another object before first usage.
- *
- * Opposite to construction, it is not possible to copy an object
- * which has an invalid dof handler pointer. This is to guarantee
- * that every iterator which is once assigned to is a valid
- * object. However, this assertion only holds in debug mode, when
- * the #Assert# macro is switched on.
- *
- * @author Wolfgang Bangerth, 1998
- */
-template <int dim>
-class MGDoFAccessor {
- public:
- /**
- * Constructor
- */
- MGDoFAccessor () : mg_dof_handler(0) {
- Assert (false, ExcInvalidObject());
- };
-
- /**
- * This should be the default constructor.
- * We cast away the #const#ness of the
- * pointer which clearly is EVIL but
- * we can't help without making all
- * functions which could somehow use
- * iterators (directly or indirectly) make
- * non-const, even if they preserve
- * constness.
- */
- MGDoFAccessor (const MGDoFHandler<dim> *mg_dof_handler) :
- mg_dof_handler(const_cast<MGDoFHandler<dim>*>(mg_dof_handler)) {};
-
- /**
- * Reset the DoF handler pointer.
- */
- void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
- Assert (dh != 0, DoFAccessor<dim>::ExcInvalidObject());
- mg_dof_handler = dh;
- };
-
- /**
- * Exception for child classes
- */
- DeclException0 (ExcInvalidObject);
-
- /**
- * Copy operator.
- */
- MGDoFAccessor<dim> & operator = (const MGDoFAccessor<dim> &da) {
- set_dof_handler (da.mg_dof_handler);
- return *this;
- };
-
- protected:
- /**
- * Store the address of the #MGDoFHandler# object
- * to be accessed.
- */
- MGDoFHandler<dim> *mg_dof_handler;
-};
-
-
-/* -------------------------------------------------------------------------- */
-
-/**
- * This is a switch class which only declares a #typdef#. It is meant to
- * determine which class a #MGDoFObjectAccessor# class is to be derived
- * from. By default, #MGDoFObjectAccessor<celldim,dim># derives from
- * the #typedef# in the general #MGDoFObjectAccessor_Inheritance<celldim,dim>#
- * class, which is #DoFObjectAccessor<celldim,dim>#,
- * but if #celldim==dim#, then the specialization #MGDoFObjectAccessor_Inheritance<dim,dim>#
- * is used which declares its local type to be #DoFCellAccessor<dim>#. Therefore,
- * the inheritance is automatically chosen to be from #DoFCellAccessor# if the
- * object under consideration has full dimension, i.e. constitutes a cell.
- *
- * @author Wolfgang Bangerth, 1999
- */
-template <int celldim, int dim>
-class MGDoFObjectAccessor_Inheritance
-{
- /**
- * Declaration of the #typedef#.
- * See the full documentation for
- * more information.
- */
- typedef DoFObjectAccessor<celldim,dim> BaseClass;
-};
-
-
-
-/**
- * This is a switch class which only declares a #typdef#. It is meant to
- * determine which class a #DoFObjectAccessor# class is to be derived
- * from. By default, #DoFObjectAccessor<celldim,dim># derives from
- * the #typedef# in the general #DoFObjectAccessor_Inheritance<celldim,dim>#
- * class, which is #TriaObjectAccessor<celldim,dim>#,
- * but if #celldim==dim#, then the specialization #TriaObjectAccessor<dim,dim>#
- * is used which declares its local type to be #CellAccessor<dim>#. Therefore,
- * the inheritance is automatically chosen to be from #CellAccessor# if the
- * object under consideration has full dimension, i.e. constitutes a cell.
- *
- * @author Wolfgang Bangerth, 1999
- */
-template <int dim>
-class MGDoFObjectAccessor_Inheritance<dim,dim>
-{
- /**
- * Declaration of the #typedef#.
- * See the full documentation for
- * more information.
- */
- typedef DoFCellAccessor<dim> BaseClass;
-};
-
-
-
-
-/* -------------------------------------------------------------------------- */
-
-
-
-/**
- * Common template for line, quad, hex.
- *
- * Internal: inheritance is necessary for the general template due to
- * a compiler error.
- * @author Guido Kanschat, 1999
- */
-template <int celldim, int dim>
-class MGDoFObjectAccessor : public MGDoFAccessor<dim>,
- public MGDoFObjectAccessor_Inheritance<celldim, dim>::BaseClass
-{};
-
-
-
-/**
- * Closure class.
- */
-template<int dim>
-class MGDoFObjectAccessor<0, dim>
-{
- public:
- typedef void AccessorData;
-
- /**
- * Constructor. Should never be called
- * and thus throws an error.
- */
- MGDoFObjectAccessor (Triangulation<dim> *,
- const int,
- const int,
- const AccessorData *)
- {
- Assert (false, ExcInternalError());
- }
-};
-
-
-
-/**
- * Grant access to the degrees of freedom located on lines.
- * This class follows mainly the route laid out by the accessor library
- * declared in the triangulation library (\Ref{TriaAccessor}). It enables
- * the user to access the degrees of freedom on the lines (there are similar
- * versions for the DoFs on quads, etc), where the dimension of the underlying
- * triangulation does not really matter (i.e. this accessor works with the
- * lines in 1D-, 2D-, etc dimensions).
- *
- *
- * \subsection{Usage}
- *
- * The \Ref{DoFDimensionInfo} classes inherited by the \Ref{DoFHandler} classes
- * declare typedefs to iterators using the accessors declared in this class
- * hierarchy tree. Usage is best to happens through these typedefs, since they
- * are more secure to changes in the class naming and template interface as well
- * as they provide easier typing (much less complicated names!).
- *
- *
- * \subsection{Notes about the class hierarchy structure}
- *
- * Inheritance from #MGDoFObjectAccessor_Inheritance<1,dim>::BaseClass# yields
- * inheritance from #DoFCellAccessor<1># if #dim==1# and from
- * #TriaObjectAccessor<1,dim># for all other #dim# values. Thus, an object
- * of this class shares all features of cells in one dimension, but behaves
- * like an ordinary line in all other cases.
- *
- * @author Wolfgang Bangerth, 1998
- */
-template <int dim>
-class MGDoFObjectAccessor<1, dim> : public MGDoFAccessor<dim>,
- public MGDoFObjectAccessor_Inheritance<1, dim>::BaseClass
-{
- public:
- /**
- * Declare the data type that this accessor
- * class expects to get passed from the
- * iterator classes.
- */
- typedef MGDoFHandler<dim> AccessorData;
-
- /**
- * Default constructor, unused thus
- * not implemented.
- */
- MGDoFObjectAccessor ();
-
- /**
- * Constructor. The #local_data#
- * argument is assumed to be a pointer
- * to a #MGDoFHandler<dim># object.
- */
- MGDoFObjectAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data);
-
- /**
- * Return the index of the #i#th degree
- * of freedom of this line on the level
- * this line lives on.
- */
- int mg_dof_index (const unsigned int i) const;
-
- /**
- * Set the index of the #i#th degree
- * of freedom of this line on the
- * level this line lives on to #index#.
- */
- void set_mg_dof_index (const unsigned int i, const int index) const;
-
- /**
- * Return the index of the #i#th degree
- * on the #vertex#th vertex for the
- * level this line lives on.
- */
- int mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i) const;
-
- /**
- * Set the index of the #i#th degree
- * on the #vertex#th vertex to #index#
- * for the level this line lives on.
- */
- void set_mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i,
- const int index) const;
-
- /**
- * Return the indices of the dofs of this
- * line in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1,
- * dofs on line.
- *
- * It is assumed that the vector already
- * has the right size beforehand. The
- * indices refer to the local numbering
- * for the level this line lives on.
- */
- void get_mg_dof_indices (vector<int> &dof_indices) const;
-
- /**
- * Return the value of the given vector
- * restricted to the dofs of this
- * cell in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1, etc,
- * dofs on line 0, dofs on line 1, etc,
- * dofs on quad 0, etc.
- *
- * It is assumed that the vector already
- * has the right size beforehand. The
- * indices refer to the multilevel
- * numbering local to the present
- * level of this cell. The vector shall
- * therefore have the same number of
- * entries as there are degrees of
- * freedom on this level.
- */
- template <typename number>
- void get_mg_dof_values (const Vector<number> &values,
- Vector<number> &dof_values) const;
-
- /**
- * Return the #i#th child as a MGDoF line
- * iterator. This function is needed since
- * the child function of the base
- * class returns a line accessor without
- * access to the MG data.
- */
- TriaIterator<dim,MGDoFObjectAccessor<1, dim> > child (const unsigned int) const;
-
- /**
- * Implement the copy operator needed
- * for the iterator classes.
- */
- void copy_from (const MGDoFObjectAccessor<1, dim> &a);
-};
-
-
-
-
-
-
-/**
- * Grant access to the multilevel degrees of freedom located on quads.
- *
- * @see DoFLineAccessor
- */
-template <int dim>
-class MGDoFObjectAccessor<2, dim> : public MGDoFAccessor<dim>,
- public MGDoFObjectAccessor_Inheritance<2, dim>::BaseClass
-{
- public:
- /**
- * Declare the data type that this accessor
- * class expects to get passed from the
- * iterator classes.
- */
- typedef MGDoFHandler<dim> AccessorData;
-
- /**
- * Default constructor, unused thus
- * not implemented.
- */
- MGDoFObjectAccessor ();
-
- /**
- * Constructor. The #local_data#
- * argument is assumed to be a pointer
- * to a #DoFHandler<dim># object.
- */
- MGDoFObjectAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data);
-
- /**
- * Return the index of the #i#th degree
- * of freedom of this quad on the level
- * this quad lives on.
- */
- int mg_dof_index (const unsigned int i) const;
-
- /**
- * Set the index of the #i#th degree
- * of freedom of this quad on the
- * level this quad lives on to #index#.
- */
- void set_mg_dof_index (const unsigned int i, const int index) const;
-
- /**
- * Return the index of the #i#th degree
- * on the #vertex#th vertex for the level
- * this quad lives on.
- */
- int mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i) const;
-
- /**
- * Set the index of the #i#th degree
- * on the #vertex#th vertex for the
- * level this quad lives on to #index#.
- */
- void set_mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i,
- const int index) const;
-
- /**
- * Return the indices of the dofs of this
- * quad in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1, etc,
- * dofs on line 0, dofs on line 1, etc,
- * dofs on quad 0, etc.
- *
- * It is assumed that the vector already
- * has the right size beforehand. The
- * indices refer to the local numbering
- * for the level this quad lives on.
- */
- void get_mg_dof_indices (vector<int> &dof_indices) const;
-
- /**
- * Return the value of the given vector
- * restricted to the dofs of this
- * cell in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1, etc,
- * dofs on line 0, dofs on line 1, etc,
- * dofs on quad 0, etc.
- *
- * It is assumed that the vector already
- * has the right size beforehand. The
- * indices refer to the multilevel
- * numbering local to the present
- * level of this cell. The vector shall
- * therefore have the same number of
- * entries as there are degrees of
- * freedom on this level.
- */
- template <typename number>
- void get_mg_dof_values (const Vector<number> &values,
- Vector<number> &dof_values) const;
-
- /**
- * Return a pointer to the #i#th line
- * bounding this #Quad#.
- */
- TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
- line (const unsigned int i) const;
-
- /**
- * Return the #i#th child as a DoF quad
- * iterator. This function is needed since
- * the child function of the base
- * class returns a quad accessor without
- * access to the DoF data.
- */
- TriaIterator<dim,MGDoFObjectAccessor<2, dim> >
- child (const unsigned int) const;
-
- /**
- * Implement the copy operator needed
- * for the iterator classes.
- */
- void copy_from (const MGDoFObjectAccessor<2, dim> &a);
-};
-
-
-
-
-/**
- * Grant access to the multilevel degrees of freedom located on hexhedra.
- *
- * @see DoFLineAccessor
- */
-template <int dim>
-class MGDoFObjectAccessor<3, dim> : public MGDoFAccessor<dim>,
- public MGDoFObjectAccessor_Inheritance<3, dim>::BaseClass
-{
- public:
- /**
- * Declare the data type that this accessor
- * class expects to get passed from the
- * iterator classes.
- */
- typedef MGDoFHandler<dim> AccessorData;
-
- /**
- * Default constructor, unused thus
- * not implemented.
- */
- MGDoFObjectAccessor ();
-
- /**
- * Constructor. The #local_data#
- * argument is assumed to be a pointer
- * to a #DoFHandler<dim># object.
- */
- MGDoFObjectAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data);
-
- /**
- * Return the index of the #i#th degree
- * of freedom of this hex on the level
- * this quad lives on.
- */
- int mg_dof_index (const unsigned int i) const;
-
- /**
- * Set the index of the #i#th degree
- * of freedom of this hex on the
- * level this hex lives on to #index#.
- */
- void set_mg_dof_index (const unsigned int i, const int index) const;
-
- /**
- * Return the index of the #i#th degree
- * on the #vertex#th vertex for the level
- * this hex lives on.
- */
- int mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i) const;
-
- /**
- * Set the index of the #i#th degree
- * on the #vertex#th vertex for the
- * level this hex lives on to #index#.
- */
- void set_mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i,
- const int index) const;
-
- /**
- * Return the indices of the dofs of this
- * hex in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1, etc,
- * dofs on line 0, dofs on line 1, etc,
- * dofs on quad 0, etc.
- *
- * It is assumed that the vector already
- * has the right size beforehand. The
- * indices refer to the local numbering
- * for the level this hex lives on.
- */
- void get_mg_dof_indices (vector<int> &dof_indices) const;
-
- /**
- * Return the value of the given vector
- * restricted to the dofs of this
- * cell in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1, etc,
- * dofs on line 0, dofs on line 1, etc,
- * dofs on quad 0, etc.
- *
- * It is assumed that the vector already
- * has the right size beforehand. The
- * indices refer to the multilevel
- * numbering local to the present
- * level of this cell. The vector shall
- * therefore have the same number of
- * entries as there are degrees of
- * freedom on this level.
- */
- template <typename number>
- void get_mg_dof_values (const Vector<number> &values,
- Vector<number> &dof_values) const;
-
- /**
- * Return a pointer to the #i#th line
- * bounding this #Hex#.
- */
- TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
- line (const unsigned int i) const;
-
- /**
- * Return a pointer to the #i#th quad
- * bounding this #Hex#.
- */
- TriaIterator<dim,MGDoFObjectAccessor<2, dim> >
- quad (const unsigned int i) const;
-
- /**
- * Return the #i#th child as a DoF quad
- * iterator. This function is needed since
- * the child function of the base
- * class returns a hex accessor without
- * access to the DoF data.
- */
- TriaIterator<dim,MGDoFObjectAccessor<3, dim> > child (const unsigned int) const;
-
- /**
- * Implement the copy operator needed
- * for the iterator classes.
- */
- void copy_from (const MGDoFObjectAccessor<3, dim> &a);
-};
-
-
-
-/**
- * Grant access to the degrees of freedom on a cell. In fact, since all
- * access to the degrees of freedom has been enabled by the classes
- * #DoFObjectAccessor<1, 1># and #DoFObjectAccessor<2, 2># for the space dimension
- * one and two, respectively, this class only collects the pieces
- * together by deriving from the appropriate #DoF*Accessor# and the
- * right #CellAccessor<dim># and finally adding two functions which give
- * access to the neighbors and children as #DoFCellAccessor# objects
- * rather than #CellAccessor# objects (the latter function was inherited
- * from the #CellAccessor<dim># class).
- *
- * @author Wolfgang Bangerth, 1998
- */
-template <int dim>
-class MGDoFCellAccessor : public MGDoFObjectAccessor<dim, dim> {
- public:
- /**
- * Type of faces.
- */
- typedef
- TriaIterator<dim, MGDoFObjectAccessor<dim-1, dim> >
- face_iterator;
- /**
- * Declare the data type that this accessor
- * class expects to get passed from the
- * iterator classes.
- */
- typedef typename MGDoFObjectAccessor<dim, dim>::AccessorData AccessorData;
-
- /**
- * Constructor
- */
- MGDoFCellAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data) :
- MGDoFObjectAccessor<dim, dim> (tria,level,index,local_data) {};
-
- /**
- * Return the #i#th neighbor as a MGDoF cell
- * iterator. This function is needed since
- * the neighbor function of the base
- * class returns a cell accessor without
- * access to the MGDoF data.
- */
- TriaIterator<dim,MGDoFCellAccessor<dim> > neighbor (const unsigned int) const;
-
- /**
- * Return the #i#th child as a MGDoF cell
- * iterator. This function is needed since
- * the child function of the base
- * class returns a cell accessor without
- * access to the DoF data.
- */
- TriaIterator<dim,MGDoFCellAccessor<dim> > child (const unsigned int) const;
-
- /**
- * Return an iterator to the #i#th face
- * of this cell.
- *
- * This function is not implemented in 1D,
- * and maps to MGDoFObjectAccessor<2, dim>::line in 2D.
- */
- face_iterator
- face (const unsigned int i) const;
-
- /**
- * Exception
- */
- DeclException0 (ExcNotUsefulForThisDimension);
-};
-
-
-
-
-
-
-
-/*---------------------------- mg_dof_accessor.h ---------------------------*/
-/* end of #ifndef __mg_dof_accessor_H */
-#endif
-/*---------------------------- mg_dof_accessor.h ---------------------------*/
+++ /dev/null
-/*---------------------------- mg_dof.h ---------------------------*/
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-#ifndef __mg_dof_H
-#define __mg_dof_H
-/*---------------------------- mg_dof.h ---------------------------*/
-
-
-#include <dofs/dof_handler.h>
-
-
-
-/**
- * Define some types which differ between the dimensions. This class
- * is analogous to the \Ref{TriaDimensionInfo} class hierarchy.
- *
- * @see MGDoFDimensionInfo<1>
- * @see MGDoFDimensionInfo<2>
- */
-template <int dim>
-class MGDoFDimensionInfo;
-
-
-
-
-
-/**
- * Define some types for the DoF handling in one dimension.
- *
- * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
- */
-class MGDoFDimensionInfo<1> {
- public:
- typedef TriaRawIterator<1,MGDoFCellAccessor<1> > raw_line_iterator;
- typedef TriaIterator<1,MGDoFCellAccessor<1> > line_iterator;
- typedef TriaActiveIterator<1,MGDoFCellAccessor<1> > active_line_iterator;
-
- typedef void * raw_quad_iterator;
- typedef void * quad_iterator;
- typedef void * active_quad_iterator;
-
- typedef void * raw_hex_iterator;
- typedef void * hex_iterator;
- typedef void * active_hex_iterator;
-
- typedef raw_line_iterator raw_cell_iterator;
- typedef line_iterator cell_iterator;
- typedef active_line_iterator active_cell_iterator;
-
- typedef void * raw_face_iterator;
- typedef void * face_iterator;
- typedef void * active_face_iterator;
-};
-
-
-
-
-
-
-/**
- * Define some types for the DoF handling in two dimensions.
- *
- * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
- */
-class MGDoFDimensionInfo<2> {
- public:
- typedef TriaRawIterator<2,MGDoFObjectAccessor<1, 2> > raw_line_iterator;
- typedef TriaIterator<2,MGDoFObjectAccessor<1, 2> > line_iterator;
- typedef TriaActiveIterator<2,MGDoFObjectAccessor<1, 2> > active_line_iterator;
-
- typedef TriaRawIterator<2,MGDoFCellAccessor<2> > raw_quad_iterator;
- typedef TriaIterator<2,MGDoFCellAccessor<2> > quad_iterator;
- typedef TriaActiveIterator<2,MGDoFCellAccessor<2> > active_quad_iterator;
-
- typedef void * raw_hex_iterator;
- typedef void * hex_iterator;
- typedef void * active_hex_iterator;
-
- typedef raw_quad_iterator raw_cell_iterator;
- typedef quad_iterator cell_iterator;
- typedef active_quad_iterator active_cell_iterator;
-
- typedef raw_line_iterator raw_face_iterator;
- typedef line_iterator face_iterator;
- typedef active_line_iterator active_face_iterator;
-};
-
-
-
-
-/**
- * Define some types for the DoF handling in two dimensions.
- *
- * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
- */
-class MGDoFDimensionInfo<3> {
- public:
- typedef TriaRawIterator<3,MGDoFObjectAccessor<1, 3> > raw_line_iterator;
- typedef TriaIterator<3,MGDoFObjectAccessor<1, 3> > line_iterator;
- typedef TriaActiveIterator<3,MGDoFObjectAccessor<1, 3> > active_line_iterator;
-
- typedef TriaRawIterator<3,MGDoFObjectAccessor<2, 3> > raw_quad_iterator;
- typedef TriaIterator<3,MGDoFObjectAccessor<2, 3> > quad_iterator;
- typedef TriaActiveIterator<3,MGDoFObjectAccessor<2, 3> > active_quad_iterator;
-
- typedef TriaRawIterator<3,MGDoFCellAccessor<3> > raw_hex_iterator;
- typedef TriaIterator<3,MGDoFCellAccessor<3> > hex_iterator;
- typedef TriaActiveIterator<3,MGDoFCellAccessor<3> > active_hex_iterator;
-
- typedef raw_hex_iterator raw_cell_iterator;
- typedef hex_iterator cell_iterator;
- typedef active_hex_iterator active_cell_iterator;
-
- typedef raw_quad_iterator raw_face_iterator;
- typedef quad_iterator face_iterator;
- typedef active_quad_iterator active_face_iterator;
-};
-
-
-
-
-/**
- * This class manages degrees of freedom for a multilevel hierarchy of
- * grids. It does mostly the same as does the #DoDHandler# class, but
- * it uses a separate enumeration of the degrees of freedom on each
- * level. For example, a vertex has several DoF numbers, one for each
- * level of the triangulation on which it exists.
- *
- * At present, multilevel algorithms are not fully functional, so this
- * documentation is still very brief.
-//TODO: Extend MGDoFHandler doc
- *
- * @author Wolfgang Bangerth, 1998, 1999
- */
-template <int dim>
-class MGDoFHandler : public DoFHandler<dim>
-{
- public:
- typedef typename MGDoFDimensionInfo<dim>::raw_line_iterator raw_line_iterator;
- typedef typename MGDoFDimensionInfo<dim>::line_iterator line_iterator;
- typedef typename MGDoFDimensionInfo<dim>::active_line_iterator active_line_iterator;
-
- typedef typename MGDoFDimensionInfo<dim>::raw_quad_iterator raw_quad_iterator;
- typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
- typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
-
- typedef typename MGDoFDimensionInfo<dim>::raw_hex_iterator raw_hex_iterator;
- typedef typename MGDoFDimensionInfo<dim>::hex_iterator hex_iterator;
- typedef typename MGDoFDimensionInfo<dim>::active_hex_iterator active_hex_iterator;
-
- typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
- typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
- typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
-
- typedef typename MGDoFDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
- typedef typename MGDoFDimensionInfo<dim>::face_iterator face_iterator;
- typedef typename MGDoFDimensionInfo<dim>::active_face_iterator active_face_iterator;
-
- /**
- * Constructor. Take #tria# as the
- * triangulation to work on.
- */
- MGDoFHandler (Triangulation<dim> &tria);
-
- /**
- * Destructor
- */
- virtual ~MGDoFHandler ();
-
- /**
- * Go through the triangulation and
- * distribute the degrees of freedoms
- * needed for the given finite element
- * according to the given distribution
- * method. We first call the #DoFHandler#'s
- * function and then distribute the
- * levelwise numbers.
- *
- * A copy of the transferred finite
- * element is stored.
- *
- * This function uses the user flags of the
- * triangulation object, so make sure you
- * don't need them after calling this
- * function, or if so store them.
- */
- virtual void distribute_dofs (const FiniteElement<dim> &, unsigned int offset=0);
-
- /**
- * Actually do the renumbering based on
- * a list of new dof numbers for all the
- * dofs.
- *
- * #new_numbers# is an array of integers
- * with size equal to the number of dofs
- * on the present level. It stores the new
- * indices after renumbering in the
- * order of the old indices.
- */
- void renumber_dofs (const unsigned int level,
- const vector<int> &new_numbers);
-
- /*--------------------------------------*/
-
- /**
- * @name Cell iterator functions
- */
- /*@{*/
- /**
- * Iterator to the first cell, used
- * or not, on level #level#. If a level
- * has no cells, a past-the-end iterator
- * is returned.
- *
- * This function calls #begin_raw_line#
- * in 1D and #begin_raw_quad# in 2D.
- */
- raw_cell_iterator begin_raw (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first used cell
- * on level #level#.
- *
- * This function calls #begin_line#
- * in 1D and #begin_quad# in 2D.
- */
- cell_iterator begin (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first active
- * cell on level #level#.
- *
- * This function calls #begin_active_line#
- * in 1D and #begin_active_quad# in 2D.
- */
- active_cell_iterator begin_active(const unsigned int level = 0) const;
-
- /**
- * Iterator past the end; this
- * iterator serves for comparisons of
- * iterators with past-the-end or
- * before-the-beginning states.
- *
- * This function calls #end_line#
- * in 1D and #end_quad# in 2D.
- */
- raw_cell_iterator end () const;
-
- /**
- * Return an iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- cell_iterator end (const unsigned int level) const;
-
- /**
- * Return a raw iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- raw_cell_iterator end_raw (const unsigned int level) const;
-
- /**
- * Return an active iterator which is the
- * first iterator not on level. If #level#
- * is the last level, then this returns
- * #end()#.
- */
- active_cell_iterator end_active (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the
- * last cell, used or not.
- *
- * This function calls #last_raw_line#
- * in 1D and #last_raw_quad# in 2D.
- */
- raw_cell_iterator last_raw () const;
-
- /**
- * Return an iterator pointing to the last
- * cell of the level #level#, used or not.
- *
- * This function calls #last_raw_line#
- * in 1D and #last_raw_quad# in 2D.
- */
- raw_cell_iterator last_raw (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * used cell.
- *
- * This function calls #last_line#
- * in 1D and #last_quad# in 2D.
- */
- cell_iterator last () const;
-
- /**
- * Return an iterator pointing to the last
- * used cell on level #level#.
- *
- * This function calls #last_line#
- * in 1D and #last_quad# in 2D.
- */
- cell_iterator last (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * active cell.
- *
- * This function calls #last_active_line#
- * in 1D and #last_active_quad# in 2D.
- */
- active_cell_iterator last_active () const;
-
- /**
- * Return an iterator pointing to the last
- * active cell on level #level#.
- *
- * This function calls #last_active_line#
- * in 1D and #last_active_quad# in 2D.
- */
- active_cell_iterator last_active (const unsigned int level) const;
- //@}
-
- /*---------------------------------------*/
-
- /**
- * @name Face iterator functions
- */
- /*@{*/
- /**
- * Iterator to the first face, used
- * or not, on level #level#. If a level
- * has no faces, a past-the-end iterator
- * is returned.
- *
- * This function calls #begin_raw_line#
- * in 2D and #begin_raw_quad# in 3D.
- */
- raw_face_iterator begin_raw_face (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first used face
- * on level #level#.
- *
- * This function calls #begin_line#
- * in 2D and #begin_quad# in 3D.
- */
- face_iterator begin_face (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first active
- * face on level #level#.
- *
- * This function calls #begin_active_line#
- * in 2D and #begin_active_quad# in 3D.
- */
- active_face_iterator begin_active_face(const unsigned int level = 0) const;
-
- /**
- * Iterator past the end; this
- * iterator serves for comparisons of
- * iterators with past-the-end or
- * before-the-beginning states.
- *
- * This function calls #end_line#
- * in 2D and #end_quad# in 3D.
- */
- raw_face_iterator end_face () const;
-
- /**
- * Return an iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- face_iterator end_face (const unsigned int level) const;
-
- /**
- * Return a raw iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- raw_face_iterator end_raw_face (const unsigned int level) const;
-
- /**
- * Return an active iterator which is the
- * first iterator not on level. If #level#
- * is the last level, then this returns
- * #end()#.
- */
- active_face_iterator end_active_face (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the
- * last face, used or not.
- *
- * This function calls #last_raw_line#
- * in 2D and #last_raw_quad# in 3D.
- */
- raw_face_iterator last_raw_face () const;
-
- /**
- * Return an iterator pointing to the last
- * face of the level #level#, used or not.
- *
- * This function calls #last_raw_line#
- * in 2D and #last_raw_quad# in 3D.
- */
- raw_face_iterator last_raw_face (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * used face.
- *
- * This function calls #last_line#
- * in 2D and #last_quad# in 3D.
- */
- face_iterator last_face () const;
-
- /**
- * Return an iterator pointing to the last
- * used face on level #level#.
- *
- * This function calls #last_line#
- * in 2D and #last_quad# in 3D.
- */
- face_iterator last_face (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * active face.
- *
- * This function calls #last_active_line#
- * in 2D and #last_active_quad# in 3D.
- */
- active_face_iterator last_active_face () const;
-
- /**
- * Return an iterator pointing to the last
- * active face on level #level#.
- *
- * This function calls #last_active_line#
- * in 2D and #last_active_quad# in 3D.
- */
- active_face_iterator last_active_face (const unsigned int level) const;
- //@}
-
-
- /*---------------------------------------*/
-
- /**
- * @name Line iterator functions
- */
- /*@{*/
- /**
- * Iterator to the first line, used
- * or not, on level #level#. If a level
- * has no lines, a past-the-end iterator
- * is returned.
- */
- raw_line_iterator begin_raw_line (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first used line
- * on level #level#.
- */
- line_iterator begin_line (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first active
- * line on level #level#.
- */
- active_line_iterator begin_active_line(const unsigned int level = 0) const;
-
- /**
- * Iterator past the end; this
- * iterator serves for comparisons of
- * iterators with past-the-end or
- * before-the-beginning states.
- */
- raw_line_iterator end_line () const;
-
- /**
- * Return an iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- line_iterator end_line (const unsigned int level) const;
-
- /**
- * Return a raw iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- raw_line_iterator end_raw_line (const unsigned int level) const;
-
- /**
- * Return an active iterator which is the
- * first iterator not on level. If #level#
- * is the last level, then this returns
- * #end()#.
- */
- active_line_iterator end_active_line (const unsigned int level) const;
-
-
- /**
- * Return an iterator pointing to the
- * last line, used or not.
- */
- raw_line_iterator last_raw_line () const;
-
- /**
- * Return an iterator pointing to the last
- * line of the level #level#, used or not.
-
- */
- raw_line_iterator last_raw_line (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * used line.
- */
- line_iterator last_line () const;
-
- /**
- * Return an iterator pointing to the last
- * used line on level #level#.
- */
- line_iterator last_line (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * active line.
- */
- active_line_iterator last_active_line () const;
-
- /**
- * Return an iterator pointing to the last
- * active line on level #level#.
- */
- active_line_iterator last_active_line (const unsigned int level) const;
- /*@}*/
-
- /*---------------------------------------*/
-
- /**
- * @name Quad iterator functions*/
- /*@{
- */
- /**
- * Iterator to the first quad, used
- * or not, on level #level#. If a level
- * has no quads, a past-the-end iterator
- * is returned.
- */
- raw_quad_iterator begin_raw_quad (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first used quad
- * on level #level#.
- */
- quad_iterator begin_quad (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first active
- * quad on level #level#.
- */
- active_quad_iterator begin_active_quad(const unsigned int level = 0) const;
-
- /**
- * Iterator past the end; this
- * iterator serves for comparisons of
- * iterators with past-the-end or
- * before-the-beginning states.
- */
- raw_quad_iterator end_quad () const;
-
- /**
- * Return an iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- quad_iterator end_quad (const unsigned int level) const;
-
- /**
- * Return a raw iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- raw_quad_iterator end_raw_quad (const unsigned int level) const;
-
- /**
- * Return an active iterator which is the
- * first iterator not on level. If #level#
- * is the last level, then this returns
- * #end()#.
- */
- active_quad_iterator end_active_quad (const unsigned int level) const;
-
-
- /**
- * Return an iterator pointing to the
- * last quad, used or not.
- */
- raw_quad_iterator last_raw_quad () const;
-
- /**
- * Return an iterator pointing to the last
- * quad of the level #level#, used or not.
-
- */
- raw_quad_iterator last_raw_quad (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * used quad.
- */
- quad_iterator last_quad () const;
-
- /**
- * Return an iterator pointing to the last
- * used quad on level #level#.
- */
- quad_iterator last_quad (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * active quad.
- */
- active_quad_iterator last_active_quad () const;
-
- /**
- * Return an iterator pointing to the last
- * active quad on level #level#.
- */
- active_quad_iterator last_active_quad (const unsigned int level) const;
- /*@}*/
-
- /*---------------------------------------*/
-
- /**
- * @name Hex iterator functions*/
- /*@{
- */
- /**
- * Iterator to the first hex, used
- * or not, on level #level#. If a level
- * has no hexs, a past-the-end iterator
- * is returned.
- */
- raw_hex_iterator begin_raw_hex (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first used hex
- * on level #level#.
- */
- hex_iterator begin_hex (const unsigned int level = 0) const;
-
- /**
- * Iterator to the first active
- * hex on level #level#.
- */
- active_hex_iterator begin_active_hex(const unsigned int level = 0) const;
-
- /**
- * Iterator past the end; this
- * iterator serves for comparisons of
- * iterators with past-the-end or
- * before-the-beginning states.
- */
- raw_hex_iterator end_hex () const;
-
- /**
- * Return an iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- hex_iterator end_hex (const unsigned int level) const;
-
- /**
- * Return a raw iterator which is the first
- * iterator not on level. If #level# is
- * the last level, then this returns
- * #end()#.
- */
- raw_hex_iterator end_raw_hex (const unsigned int level) const;
-
- /**
- * Return an active iterator which is the
- * first iterator not on level. If #level#
- * is the last level, then this returns
- * #end()#.
- */
- active_hex_iterator end_active_hex (const unsigned int level) const;
-
-
- /**
- * Return an iterator pointing to the
- * last hex, used or not.
- */
- raw_hex_iterator last_raw_hex () const;
-
- /**
- * Return an iterator pointing to the last
- * hex of the level #level#, used or not.
-
- */
- raw_hex_iterator last_raw_hex (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * used hex.
- */
- hex_iterator last_hex () const;
-
- /**
- * Return an iterator pointing to the last
- * used hex on level #level#.
- */
- hex_iterator last_hex (const unsigned int level) const;
-
- /**
- * Return an iterator pointing to the last
- * active hex.
- */
- active_hex_iterator last_active_hex () const;
-
- /**
- * Return an iterator pointing to the last
- * active hex on level #level#.
- */
- active_hex_iterator last_active_hex (const unsigned int level) const;
- /*@}*/
-
- /*---------------------------------------*/
-
- /**
- * Return the number of degrees of freedom
- * on the specified level.
- * Included in this number are those
- * DoFs which are constrained by
- * hanging nodes.
- */
- unsigned int n_dofs (const unsigned int level) const;
-
- /**
- * Exception.
- */
- DeclException1 (ExcInvalidLevel,
- int,
- << "The level index " << arg1 << "was not in the valid range.");
-
- private:
-
- /** We need for each vertex a list of the
- * degree of freedom indices on each of the
- * levels this vertex lives on. Since most
- * vertices live only on a few levels, it
- * is not economical to reserve indices for
- * all the levels there are; rather, we
- * create an object which holds the indices
- * on those levels only where the vertex
- * lives. To construct such an array, it
- * is necessary to know beforehand which
- * is the coarsest level the vertex lives
- * on, how many levels it lives on and
- * how many dofs there are on each vertex.
- * If we have this information, we can
- * allocate exactly the amount of memory
- * which is needed and need not handle
- * growing arrays and the like.
- */
- class MGVertexDoFs {
- public:
-
- /**
- * Constructor. This one is empty
- * because it is difficult to make it
- * efficient to use vector<>'s and
- * still construct the object using
- * the constructor. Use the #init#
- * function to really allocate
- * memory.
- */
- MGVertexDoFs ();
-
- /**
- * Allocate memory and
- * set all indices to #-1#.
- */
- void init (const unsigned int coarsest_level,
- const unsigned int finest_level,
- const unsigned int dofs_per_vertex);
-
- /**
- * Destructor
- */
- ~MGVertexDoFs ();
-
- /**
- * Set the index with number
- * #dof_number# of this vertex on
- * #level# to the given index. To
- * compute the position in the array,
- * one has to specify how many
- * dofs per vertex there are. It is
- * not checked that the level number
- * is below the number of the finest
- * level this vertex lives on.
- *
- * The function is inline, so should
- * be reasonably fast.
- */
- void set_index (const unsigned int level,
- const unsigned int dof_number,
- const unsigned int dofs_per_vertex,
- const unsigned int index);
-
- /**
- * Return the index with number
- * #dof_number# of this vertex on
- * #level#. To
- * compute the position in the array,
- * one has to specify how many
- * dofs per vertex there are. It is
- * not checked that the level number
- * is below the number of the finest
- * level this vertex lives on.
- *
- * The function is inline, so should
- * be reasonably fast.
- */
- int get_index (const unsigned int level,
- const unsigned int dof_number,
- const unsigned int dofs_per_vertex) const;
-
- /**
- * Return the index of the coarsest
- * level this vertex lives on.
- */
- unsigned int get_coarsest_level () const;
-
- /**
- * Return the index of the finest
- * level this vertex lives on.
- */
- unsigned int get_finest_level () const;
-
- /**
- * Exception.
- */
- DeclException0 (ExcNoMemory);
- /**
- * Exception.
- */
- DeclException1 (ExcInvalidLevel,
- int,
- << "The given level number " << arg1 << " is outside "
- << "the range of levels this vertex lives on.");
- /**
- * Exception.
- */
- DeclException0 (ExcInternalError);
-
- private:
- /**
- * Store the coarsest level this
- * vertex lives on. This is used
- * as an offset when accessing the
- * dofs of a specific level.
- */
- unsigned int coarsest_level;
-
- /**
- * Finest level this level lives on.
- * This is mostly used for error
- * checking but can also be accessed
- * by the function #get_finest_level#.
- */
- unsigned int finest_level;
-
- /**
- * Array holding the indices.
- */
- int *indices;
- };
-
-
- /**
- * Distribute dofs on the given cell,
- * with new dofs starting with index
- * #next_free_dof#. Return the next
- * unused index number. The finite
- * element used is the one given to
- * #distribute_dofs#, which is copied
- * to #selected_fe#.
- *
- * This function is excluded from the
- * #distribute_dofs# function since
- * it can not be implemented dimension
- * independent.
- *
- * Note that unlike for the usual dofs,
- * here all cells and not only active
- * ones are allowed.
- */
- unsigned int distribute_dofs_on_cell (cell_iterator &cell,
- unsigned int next_free_dof);
-
- /**
- * Reserve enough space for the MG dof
- * indices for a given triangulation.
- */
- void reserve_space ();
-
- /**
- * Space to store the DoF numbers for the
- * different levels. Unlike the #levels#
- * object in the #DoFHandler#, these are
- * not global numbers but rather are
- * numbers which start from zero on each
- * level.
- */
- vector<DoFLevel<dim>*> mg_levels;
-
- /**
- * For each vertex there is a list of
- * indices of the degrees of freedom indices
- * on the different levels it lives on and
- * which are these levels.
- */
- vector<MGVertexDoFs> mg_vertex_dofs;
-
- /**
- * Vectors storing the number of degrees of
- * freedom on each level.
- */
- vector<unsigned int> mg_used_dofs;
-
-#if (__GNUC__==2) && (__GNUC_MINOR__ < 95)
- // this seems to be disallowed
- // by the standard, so gcc2.95
- // does not accept it
- friend class MGDoFObjectAccessor<1, dim>;
- friend class MGDoFObjectAccessor<2, dim>;
- friend class MGDoFObjectAccessor<3, dim>;
-#else
- // this, however, may grant
- // access to too many classes,
- // but ...
- template <int dim1, int dim2> friend class MGDoFObjectAccessor;
-#endif
-};
-
-
-
-
-
-
-
-/* ----------------------- Inline functions of MGVertexDoFs -------------------*/
-
-
-template <int dim>
-inline
-void MGDoFHandler<dim>::MGVertexDoFs::set_index (const unsigned int level,
- const unsigned int dof_number,
- const unsigned int dofs_per_vertex,
- const unsigned int index) {
- Assert ((level >= coarsest_level) && (level <= finest_level),
- ExcInvalidLevel(level));
- Assert (dof_number < dofs_per_vertex,
- ExcIndexRange(dof_number, 0, dofs_per_vertex));
-
- indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index;
-};
-
-
-
-
-template <int dim>
-inline
-int MGDoFHandler<dim>::MGVertexDoFs::get_index (const unsigned int level,
- const unsigned int dof_number,
- const unsigned int dofs_per_vertex) const {
- Assert ((level >= coarsest_level) && (level <= finest_level),
- ExcInvalidLevel(level));
- Assert (dof_number < dofs_per_vertex,
- ExcIndexRange (dof_number, 0, dofs_per_vertex));
-
- return indices[(level-coarsest_level)*dofs_per_vertex + dof_number];
-};
-
-
-
-/*---------------------------- mg_dof.h ---------------------------*/
-/* end of #ifndef __mg_dof_H */
-#endif
-/*---------------------------- mg_dof.h ---------------------------*/
+++ /dev/null
-/*---------------------------- mg_dof_tools.h ---------------------------*/
-/* $Id$ */
-#ifndef __mg_dof_tools_H
-#define __mg_dof_tools_H
-/*---------------------------- mg_dof_tools.h ---------------------------*/
-
-
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
-
-
-/**
- * This is a collection of functions operating on, and manipulating
- * the numbers of degrees of freedom in a multilevel triangulation. It
- * is similar in purpose and function to the #DoFTools# class, but
- * operates on #MGDoFHandler# objects instead of #DoFHandler#
- * objects. See there and the documentation of the member functions
- * for more information.
- *
- * All member functions are static, so there is no need to create an
- * object of class #DoFTools#.
- *
- * @author Wolfgang Bangerth and others, 1999
- */
-class MGDoFTools
-{
- public:
- /**
- * Write the sparsity structure
- * of the matrix belonging to the
- * specified #level# including
- * constrained degrees of freedom
- * into the matrix structure. The
- * sparsity pattern does not
- * include entries introduced by
- * the elimination of constrained
- * nodes. The sparsity pattern
- * is not compressed, since if
- * you want to call
- * #ConstraintMatrix::condense(1)#
- * afterwards, new entries have
- * to be added. However, if you
- * don't want to call
- * #ConstraintMatrix::condense(1)#,
- * you have to compress the
- * matrix yourself, using
- * #SparsityPattern::compress()#.
- */
- template <int dim>
- static void
- make_sparsity_pattern (const MGDoFHandler<dim> &dof_handler,
- const unsigned int level,
- SparsityPattern &sparsity);
-};
-
-
-
-/*---------------------------- mg_dof_tools.h ---------------------------*/
-/* end of #ifndef __mg_dof_tools_H */
-#endif
-/*---------------------------- mg_dof_tools.h ---------------------------*/
-/*---------------------------- mg_base.h ---------------------------*/
+/*---------------------------- mgbase.h ---------------------------*/
/* $Id$ */
-#ifndef __mg_base_H
-#define __mg_base_H
-/*---------------------------- mg_base.h ---------------------------*/
+#ifndef __mgbase_H
+#define __mgbase_H
+/*---------------------------- mgbase.h ---------------------------*/
#include <base/subscriptor.h>
#include <base/smartpointer.h>
-#include <lac/forward_declarations.h>
#include <lac/vector.h>
+#include <lac/forward_declarations.h>
+#include <grid/forward_declarations.h>
+
#include <vector>
+
+/**
+ * An array with an object for each level. The purpose of this class
+ * is mostly to allow access by level number, even if the lower levels
+ * are not used and therefore have no object at all; this is done by
+ * simply shifting the given index by the minimum level we have
+ * stored.
+ *
+ * In most cases, the objects which are stored on each levels, are
+ * either matrices or vectors.
+ *
+ * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ */
+template<class Object>
+class MGLevelObject : public Subscriptor
+{
+ public:
+ /**
+ * Constructor allowing to
+ * initialize the number of
+ * levels. By default, the object
+ * is created empty.
+ */
+ MGLevelObject (const unsigned int minlevel = 0,
+ const unsigned int maxlevel = 0);
+
+ /**
+ * Access object on level #level#.
+ */
+ Object & operator[] (const unsigned int level);
+
+ /**
+ * Access object on level
+ * #level#. Constant version.
+ */
+ const Object & operator[] (const unsigned int level) const;
+
+ /**
+ * Delete all previous contents
+ * of this object and reset its
+ * size according to the values
+ * of #new_minlevel# and
+ * #new_maxlevel#.
+ */
+ void resize (const unsigned int new_minlevel,
+ const unsigned int new_maxlevel);
+
+ /**
+ * Call #clear# on all objects
+ * stored by this object. This
+ * function is only implemented
+ * for some #Object# classes,
+ * most notably for vectors and
+ * matrices. Note that if
+ * #Object==Vector<T>#, #clear#
+ * will set all entries to zero,
+ * while if
+ * #Object==std::vector<T>#,
+ * #clear# deletes the elements
+ * of the vectors. This class
+ * might therefore not be useful
+ * for STL vectors.
+ */
+ void clear();
+
+ //////
+ unsigned int get_minlevel () const;
+ unsigned int get_maxlevel () const;
+
+ private:
+ /**
+ * Level of first component.
+ */
+ unsigned int minlevel;
+
+ /**
+ * Array of the objects to be held.
+ */
+ vector<Object> objects;
+};
+
+
+
+
+
/**
* Abstract base class for coarse grid solvers. The interface of a
* function call operator is defined to execute coarse grid solution
-
/**
- * Abstract base class for multigrid smoothers.
- * In fact, this class only provides the interface of the smoothing function.
- * Using #deal.II# grid handling, #MGSmoother# is a good point to start.
+ * Coarse grid solver using LAC iterative methods.
+ * This is a little wrapper, transforming a triplet of iterative
+ * solver, matrix and preconditioner into a coarse grid solver.
*
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ * The type of the matrix (i.e. the template parameter #MATRIX#)
+ * should be derived from #Subscriptor# to allow for the use of a
+ * smart pointer to it.
+ *
+ * @author Guido Kanschat, 1999
*/
-class MGSmootherBase : public Subscriptor
-{
+template<class SOLVER, class MATRIX, class PRECOND>
+class MGCoarseGridLACIteration : public MGCoarseGridSolver
+{
public:
/**
- * Virtual destructor. Does
- * nothing in particular, but
- * needs to be declared anyhow.
+ * Constructor.
+ * Store solver, matrix and
+ * preconditioning method for later
+ * use.
*/
- virtual ~MGSmootherBase();
+ MGCoarseGridLACIteration (SOLVER &,
+ const MATRIX &,
+ const PRECOND &);
/**
- * Smooth the residual of #u# on
- * the given level. If $S$ is the
- * smoothing operator, then this
- * function should do the
- * following operation:
- * $u += S (rhs - Au)$, where #u# and
- * #rhs# are the input parameters.
- *
- * This function should keep the
- * interior level boundary
- * values, so you may want to
- * call #set_zero_interior_boundary#
- * in the derived class
- * #MGSmoother# somewhere in your
- * derived function, or another
- * function doing similar things.
- */
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const = 0;
+ * Implementation of the abstract
+ * function.
+ * Calls the solver method with
+ * matrix, vectors and
+ * preconditioner.
+ */
+ virtual void operator() (const unsigned int level,
+ Vector<double> &dst,
+ const Vector<double> &src) const;
+ private:
+ /**
+ * Reference to the solver.
+ */
+ SOLVER& solver;
+
+ /**
+ * Reference to the matrix.
+ */
+ const SmartPointer<const MATRIX> matrix;
+
+ /**
+ * Reference to the preconditioner.
+ */
+ const PRECOND& precondition;
};
-/**
- * Smoother doing nothing. This class is not useful for many applications other
- * than for testing some multigrid procedures. Also some applications might
- * get convergence without smoothing and then this class brings you the
- * cheapest possible multigrid.
- *
- * @author Guido Kanschat, 1999
- */
-class MGSmootherIdentity : public MGSmootherBase
-{
- public:
- /**
- * Implementation of the
- * interface in #MGSmootherBase#.
- * This function does nothing,
- * which by comparison with the
- * definition of this function
- * means that the the smoothing
- * operator equals the null
- * operator.
- */
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const;
-};
-
-
/**
-/**
- * An array with a vector for each level. The purpose of this class
- * is mostly to allow access by level number, even if the lower levels
- * are not used and therefore have no vector at all; this is done by
- * simply shifting the given index by the minimum level we have
- * stored.
- *
- * @author Guido Kanschat, 1999
- */
-template<class VECTOR = Vector<double> >
-class MGVector : public Subscriptor
-{
- public:
- /**
- * Constructor allowing to
- * initialize the number of
- * levels.
- */
- MGVector (const unsigned int minlevel,
- const unsigned int maxlevel);
-
- /**
- * Access vector on level #level#.
- */
- VECTOR & operator[] (const unsigned int level);
-
- /**
- * Access vector on level
- * #level#. Constant version.
- */
- const VECTOR & operator[] (const unsigned int level) const;
-
- /**
- * Call #clear# on all vectors
- * stored by this object. Note
- * that if #VECTOR==Vector<T>#,
- * #clear# will set all entries
- * to zero, while if
- * #VECTOR==std::vector<T>#,
- * #clear# deletes the elements
- * of the vectors. This class
- * might therefore not be useful
- * for STL vectors.
- */
- void clear();
-
- private:
- /**
- * Level of first component.
- */
- const unsigned int minlevel;
-
- /**
- * Array of the vectors to be held.
- */
- vector<VECTOR> vectors;
-};
-
-
-
-/**
- * An array of matrices for each level.
- * This class provides the functionality of #vector<MATRIX># combined
- * with a subscriptor for smart pointers.
- * @author Guido Kanschat, 1999
- */
-template <class MATRIX = SparseMatrix<double> >
-class MGMatrix : public Subscriptor,
- public vector<MATRIX>
-{
- public:
- /**
- * Constructor allowing to initialize the number of levels.
- */
- MGMatrix(unsigned int minlevel, unsigned int maxlevel);
-
- /**
- * Safe access operator.
- */
- MATRIX& operator[](unsigned int);
-
- /**
- * Safe access operator.
- */
- const MATRIX& operator[](unsigned int) const;
-
- private:
- /**
- * Level of first component.
- */
- unsigned int minlevel;
-};
-
-
-
-
-/**
- * Coarse grid solver using LAC iterative methods.
- * This is a little wrapper, transforming a triplet of iterative
- * solver, matrix and preconditioner into a coarse grid solver.
- *
- * The type of the matrix (i.e. the template parameter #MATRIX#)
- * should be derived from #Subscriptor# to allow for the use of a
- * smart pointer to it.
- *
- * @author Guido Kanschat, 1999
- */
-template<class SOLVER, class MATRIX, class PRECOND>
-class MGCoarseGridLACIteration : public MGCoarseGridSolver
-{
- public:
- /**
- * Constructor.
- * Store solver, matrix and
- * preconditioning method for later
- * use.
- */
- MGCoarseGridLACIteration (SOLVER &,
- const MATRIX &,
- const PRECOND &);
-
- /**
- * Implementation of the abstract
- * function.
- * Calls the solver method with
- * matrix, vectors and
- * preconditioner.
- */
- virtual void operator() (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src) const;
- private:
- /**
- * Reference to the solver.
- */
- SOLVER& solver;
-
- /**
- * Reference to the matrix.
- */
- const SmartPointer<const MATRIX> matrix;
-
- /**
- * Reference to the preconditioner.
- */
- const PRECOND& precondition;
-};
/**
* Auxiliary vector.
*/
- MGVector<Vector<double> > defect;
+ MGLevelObject<Vector<double> > defect;
/**
* Auxiliary vector.
*/
- MGVector<Vector<double> > solution;
+ MGLevelObject<Vector<double> > solution;
/**
* Prolongation and restriction object.
/* ------------------------------------------------------------------- */
-template<class VECTOR>
-MGVector<VECTOR>::MGVector(const unsigned int min,
- const unsigned int max)
+
+template<class Object>
+MGLevelObject<Object>::MGLevelObject(const unsigned int min,
+ const unsigned int max)
:
- minlevel(min),
- vectors(max-min+1)
-{};
+ minlevel(0)
+{
+ resize (min, max);
+};
-template<class VECTOR>
-VECTOR &
-MGVector<VECTOR>::operator[] (const unsigned int i)
+template<class Object>
+Object &
+MGLevelObject<Object>::operator[] (const unsigned int i)
{
- Assert((i>=minlevel) && (i<minlevel+vectors.size()),
- ExcIndexRange (i, minlevel, minlevel+vectors.size()));
- return vectors[i-minlevel];
+ Assert((i>=minlevel) && (i<minlevel+objects.size()),
+ ExcIndexRange (i, minlevel, minlevel+objects.size()));
+ return objects[i-minlevel];
}
-template<class VECTOR>
-const VECTOR &
-MGVector<VECTOR>::operator[] (const unsigned int i) const
+template<class Object>
+const Object &
+MGLevelObject<Object>::operator[] (const unsigned int i) const
{
- Assert((i>=minlevel) && (i<minlevel+vectors.size()),
- ExcIndexRange (i, minlevel, minlevel+vectors.size()));
- return vectors[i-minlevel];
+ Assert((i>=minlevel) && (i<minlevel+objects.size()),
+ ExcIndexRange (i, minlevel, minlevel+objects.size()));
+ return objects[i-minlevel];
}
-template<class VECTOR>
+
+template<class Object>
void
-MGVector<VECTOR>::clear()
+MGLevelObject<Object>::resize (const unsigned int new_minlevel,
+ const unsigned int new_maxlevel)
{
- typename vector<VECTOR>::iterator v;
- for (v = vectors.begin(); v != vectors.end(); ++v)
- v->clear();
-
-}
+ Assert (new_minlevel <= new_maxlevel, ExcInternalError());
+ objects.clear ();
+ minlevel = new_minlevel;
+ objects.resize (new_maxlevel - new_minlevel + 1);
+};
-/* ------------------------------------------------------------------- */
-template<class MATRIX>
-MGMatrix<MATRIX>::MGMatrix(unsigned int min, unsigned int max)
- :
- vector<MATRIX>(max-min+1),
- minlevel(min)
-{}
+template<class Object>
+void
+MGLevelObject<Object>::clear ()
+{
+ typename vector<Object>::iterator v;
+ for (v = objects.begin(); v != objects.end(); ++v)
+ v->clear();
+};
-template<class MATRIX>
-MATRIX&
-MGMatrix<MATRIX>::operator[](unsigned int i)
+template<class Object>
+unsigned int
+MGLevelObject<Object>::get_minlevel () const
{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<MATRIX>::operator[](i-minlevel);
+ return minlevel;
};
-template<class MATRIX>
-const MATRIX&
-MGMatrix<MATRIX>::operator[](unsigned int i) const
+template<class Object>
+unsigned int
+MGLevelObject<Object>::get_maxlevel () const
{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<MATRIX>::operator[](i-minlevel);
+ return minlevel + objects.size() - 1;
};
-/*---------------------------- mg_base.h ---------------------------*/
-/* end of #ifndef __mg_base_H */
+/*---------------------------- mgbase.h ---------------------------*/
+/* end of #ifndef __mgbase_H */
#endif
-/*---------------------------- mg_base.h ---------------------------*/
+/*---------------------------- mgbase.h ---------------------------*/
* #ConstraintMatrix::condense(1)#,
* you have to compress the
* matrix yourself, using
- * #SparsityPattern::compress()#.
+ * #SparseMatrixStruct::compress()#.
*/
template <int dim>
static void
#include <lac/forward_declarations.h>
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
#include <grid/forward_declarations.h>
#include <base/smartpointer.h>
#include <vector>
+
+/**
+ * Abstract base class for multigrid smoothers.
+ * In fact, this class only provides the interface of the smoothing function.
+ * Using #deal.II# grid handling, #MGSmoother# is a good point to start.
+ *
+ * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ */
+class MGSmootherBase : public Subscriptor
+{
+ public:
+ /**
+ * Virtual destructor. Does
+ * nothing in particular, but
+ * needs to be declared anyhow.
+ */
+ virtual ~MGSmootherBase();
+
+ /**
+ * Smooth the residual of #u# on
+ * the given level. If $S$ is the
+ * smoothing operator, then this
+ * function should do the
+ * following operation:
+ * $u += S (rhs - Au)$, where #u# and
+ * #rhs# are the input parameters.
+ *
+ * This function should keep the
+ * interior level boundary
+ * values, so you may want to
+ * call #set_zero_interior_boundary#
+ * in the derived class
+ * #MGSmoother# somewhere in your
+ * derived function, or another
+ * function doing similar things.
+ */
+ virtual void smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const = 0;
+};
+
+
+
+
+
+/**
+ * Smoother doing nothing. This class is not useful for many applications other
+ * than for testing some multigrid procedures. Also some applications might
+ * get convergence without smoothing and then this class brings you the
+ * cheapest possible multigrid.
+ *
+ * @author Guido Kanschat, 1999
+ */
+class MGSmootherIdentity : public MGSmootherBase
+{
+ public:
+ /**
+ * Implementation of the
+ * interface in #MGSmootherBase#.
+ * This function does nothing,
+ * which by comparison with the
+ * definition of this function
+ * means that the the smoothing
+ * operator equals the null
+ * operator.
+ */
+ virtual void smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
+};
+
+
+
/**
* Base class for multigrid smoothers. It declares the interface
* to smoothers by inheriting #MGSmootherBase# and implements some functionality for setting
vector<vector<int> > interior_boundary_dofs;
};
+
+
+
/**
* Implementation of a smoother using matrix builtin relaxation methods.
*
*/
template<int dim>
- MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<number> >& matrix,
- function_ptr function,
- unsigned int steps,
- double omega = 1.);
+ MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+ const MGLevelObject<SparseMatrix<number> > &matrix,
+ const function_ptr function,
+ const unsigned int steps,
+ const double omega = 1.);
/**
* Implementation of the interface in #MGSmootherBase#.
/**
* Pointer to the matrices.
*/
- SmartPointer< const MGMatrix< SparseMatrix<number> > > matrix;
+ SmartPointer<const MGLevelObject<SparseMatrix<number> > > matrix;
+
/**
* Pointer to the relaxation function.
*/
function_ptr relaxation;
+
/**
* Relaxation parameter.
*/
double omega;
+
/**
* Auxiliary vector.
*/
mutable Vector<double> h1;
+
/**
* Auxiliary vector.
*/
mutable Vector<double> h2;
};
+
+
+/* ------------------------------- Inline functions -------------------------- */
+
inline
void
MGSmoother::set_steps(unsigned int i)
// $Id$
+
template<typename number>
template<int dim>
MGSmootherRelaxation<number>::
-MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<number> >& matrix,
- function_ptr relaxation,
- unsigned int steps,
- double omega)
+MGSmootherRelaxation (const MGDoFHandler<dim> &mg_dof,
+ const MGLevelObject<SparseMatrix<number> > &matrix,
+ const function_ptr relaxation,
+ const unsigned int steps,
+ const double omega)
:
MGSmoother(mg_dof, steps),
matrix(&matrix),
relaxation(relaxation),
omega(omega)
-{}
+{};
+
+
template<typename number>
void
MGSmootherRelaxation<number>::smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const
+ Vector<double> &u,
+ const Vector<double> &rhs) const
{
h1.reinit(u);
h2.reinit(u);
#include <base/smartpointer.h>
#include <lac/forward_declarations.h>
#include <grid/forward_declarations.h>
+#include <lac/sparse_matrix.h>
#include <lac/vector.h>
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
#include <vector>
/**
* Implementation of multigrid with matrices.
- * While #MGBase# was only an abstract framework for the v-cycle,
+ * While #MGBase# was only an abstract framework for the V-cycle,
* here we have the implementation of the pure virtual functions defined there.
* Furthermore, level information is obtained from a triangulation object.
*
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
-template<int dim>
-class MG : public MGBase
+template <int dim>
+class Multigrid : public MGBase
{
public:
/**
* (for example for hanging
* nodes), a set of matrices
* which describe the equation
- * discretized on each level, and
- * an object mediating the
- * transfer of solutions between
- * different levels.
+ * discretized on each level,
+ * their respective sparsity
+ * patterns, and an object
+ * mediating the transfer of
+ * solutions between different
+ * levels.
*
* This function already
* initializes the vectors which
* this type as late as possible.
*/
//TODO: minlevel, maxlevel?
- MG(const MGDoFHandler<dim> &mg_dof_handler,
- ConstraintMatrix &constraints,
- const MGMatrix<SparseMatrix<double> > &level_matrices,
- const MGTransferBase &transfer,
- const unsigned int minlevel = 0,
- const unsigned int maxlevel = 1000000);
+ Multigrid(const MGDoFHandler<dim> &mg_dof_handler,
+ const ConstraintMatrix &constraints,
+ const MGLevelObject<SparsityPattern> &level_sparsities,
+ const MGLevelObject<SparseMatrix<double> > &level_matrices,
+ const MGTransferBase &transfer,
+ const unsigned int minlevel = 0,
+ const unsigned int maxlevel = 1000000);
/**
* Transfer from a vector on the
* the respective positions of a
* #Vector<double>#. In order to
* keep the result consistent,
- * constrained degrees of freedom are set to zero.
+ * constrained degrees of freedom
+ * are set to zero.
*/
template<typename number>
void copy_from_mg (Vector<number> &dst) const;
*/
SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
+
+ /**
+ * Sparsity patterns for each level.
+ */
+ SmartPointer<const MGLevelObject<SparsityPattern> > level_sparsities;
+
/**
- * Matrices for each level.
- * The matrices are prepared by
- * the constructor of #MG# and can
- * be accessed for assembling.
+ * Matrices for each level. The
+ * matrices are prepared by the
+ * constructor of #Multigrid# and
+ * can be accessed for
+ * assembling.
*/
- SmartPointer<const MGMatrix<SparseMatrix<double> > > level_matrices;
+ SmartPointer<const MGLevelObject<SparseMatrix<double> > > level_matrices;
/**
* Pointer to the object holding
private:
- /**
- * Sparsity patterns for the
- * prolongation matrices on each
- * level.
- */
- vector<SparsityPattern> prolongation_sparsities;
+ vector<SparsityPattern> prolongation_sparsities;
- /**
- * The actual prolongation matrix.
- * column indices belong to the
- * dof indices of the mother cell,
- * i.e. the coarse level.
- * while row indices belong to the
- * child cell, i.e. the fine level.
- */
+ /**
+ * The actual prolongation matrix.
+ * column indices belong to the
+ * dof indices of the mother cell,
+ * i.e. the coarse level.
+ * while row indices belong to the
+ * child cell, i.e. the fine level.
+ */
vector<SparseMatrix<float> > prolongation_matrices;
};
template<int dim>
inline
const SparseMatrix<double>&
-MG<dim>::get_matrix(unsigned int level) const
+Multigrid<dim>::get_matrix (const unsigned int level) const
{
Assert((level>=minlevel) && (level<maxlevel),
ExcIndexRange(level, minlevel, maxlevel));
#include <dofs/dof_constraints.h>
-#include <numerics/multigrid.h>
+#include <multigrid/multigrid.h>
#include <algorithm>
#include <fstream>
+#include <lac/sparse_matrix.h>
+//TODO: This function needs to be specially implemented, since in 2d mode we use faces
+#if deal_II_dimension != 1
template <int dim>
template <typename number>
void
-MG<dim>::copy_to_mg(const Vector<number>& src)
+Multigrid<dim>::copy_to_mg(const Vector<number>& src)
{
const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell;
const unsigned int dofs_per_face = mg_dof_handler->get_fe().dofs_per_face;
// constraints->condense(src);
vector<int> global_dof_indices (dofs_per_cell);
- vector<int> level_dof_indices (dofs_per_cell);
+ vector<int> level_dof_indices (dofs_per_cell);
vector<int> level_face_indices (dofs_per_face);
+ // initialize the objects with
+ // their correct size
+ for (unsigned int level=minlevel; level<=maxlevel ; ++level)
+ {
+ const unsigned int system_size = (*level_matrices)[level].m();
+
+ solution[level].reinit(system_size);
+ defect[level].reinit(system_size);
+ };
+
// traverse the grid top-down
// (i.e. starting with the most
// refined grid). this way, we can
};
};
+#endif
template <int dim>
template <typename number>
void
-MG<dim>::copy_from_mg(Vector<number> &dst) const
+Multigrid<dim>::copy_from_mg(Vector<number> &dst) const
{
const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell;
+++ /dev/null
-/*---------------------------- mg_smoother.h ---------------------------*/
-/* $Id$ */
-#ifndef __mg_smoother_H
-#define __mg_smoother_H
-/*---------------------------- mg_smoother.h ---------------------------*/
-
-
-#include <lac/forward_declarations.h>
-#include <lac/mgbase.h>
-#include <grid/forward_declarations.h>
-#include <base/smartpointer.h>
-#include <vector>
-
-
-/**
- * Base class for multigrid smoothers. It declares the interface
- * to smoothers by inheriting #MGSmootherBase# and implements some functionality for setting
- * the values
- * of vectors at interior boundaries (i.e. boundaries between differing
- * levels of the triangulation) to zero, by building a list of these degrees
- * of freedom's indices at construction time.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-class MGSmoother : public MGSmootherBase
-{
- private:
- /**
- * Default constructor. Made private
- * to prevent it being called, which
- * is necessary since this could
- * circumpass the set-up of the list
- * if interior boundary dofs.
- */
- MGSmoother ();
-
- public:
-
- /**
- * Constructor. This one collects
- * the indices of the degrees of freedom
- * on the interior boundaries between
- * the different levels, which are
- * needed by the function
- * #set_zero_interior_boundaries#.
- *
- * Since this function is implemented
- * a bit different in 1d (there are no
- * faces of cells, just vertices),
- * there are actually two sets of
- * constructors, namely this one for 1d
- * and the following one for all other
- * dimensions.
- */
- MGSmoother (const MGDoFHandler<1> &mg_dof,
- unsigned int steps);
-
- /**
- * Constructor. This one collects
- * the indices of the degrees of freedom
- * on the interior boundaries between
- * the different levels, which are
- * needed by the function
- * #set_zero_interior_boundaries#.
- *
- * The parameter steps indicates the number of smoothing
- * steps to be executed by #smooth#.
- */
- template <int dim>
- MGSmoother (const MGDoFHandler<dim> &mg_dof,
- unsigned int steps);
-
- /**
- * Reset the values of the degrees of
- * freedom on interior boundaries between
- * different levels to zero in the given
- * data vector #u#.
- *
- * Since the coarsest level (#level==0#)
- * has no interior boundaries, this
- * function does nothing in this case.
- */
- void set_zero_interior_boundary (const unsigned int level,
- Vector<double> &u) const;
- /**
- * Modify the number of smoothing steps.
- */
- void set_steps(unsigned int steps);
- /**
- * How many steps should be used?
- */
- unsigned int get_steps() const;
-
- private:
- /**
- * Number of smoothing steps.
- */
- unsigned int steps;
- /**
- * For each level, we store a list of
- * degree of freedom indices which are
- * located on interior boundaries between
- * differing levels of the triangulation.
- * Since the coarsest level has no
- * interior boundary dofs, the first
- * entry refers to the second level.
- *
- * These arrays are set by the constructor.
- * The entries for each level are sorted ascendingly.
- */
- vector<vector<int> > interior_boundary_dofs;
-};
-
-/**
- * Implementation of a smoother using matrix builtin relaxation methods.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-template<typename number>
-class MGSmootherRelaxation : public MGSmoother
-{
- public:
- /**
- * Type of the smoothing
- * function of the matrix.
- */
- typedef void
- (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
- typename SparseMatrix<number>::value_type) const;
-
- /**
- * Constructor.
- * The constructor uses an #MGDoFHandler# to initialize
- * data structures for interior level boundary handling
- * in #MGSmoother#.
- *
- * Furthermore, it takes a pointer to the
- * level matrices and their smoothing function.
- * This function must perform one relaxation step
- * like #SparseMatrix<number>::SOR_step# does. Do not
- * use the preconditioning methods because they apply
- * a preconditioning matrix to the residual.
- *
- * The final two parameters are the number of relaxation
- * steps and the relaxation parameter.
- */
-
- template<int dim>
- MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<number> >& matrix,
- function_ptr function,
- unsigned int steps,
- double omega = 1.);
-
- /**
- * Implementation of the interface in #MGSmootherBase#.
- * We use the SOR method in #SparseMatrix# for the real work
- * and find a way to keep the boundary values.
- */
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const;
- private:
- /**
- * Pointer to the matrices.
- */
- SmartPointer< const MGMatrix< SparseMatrix<number> > > matrix;
- /**
- * Pointer to the relaxation function.
- */
- function_ptr relaxation;
- /**
- * Relaxation parameter.
- */
- double omega;
- /**
- * Auxiliary vector.
- */
- mutable Vector<double> h1;
- /**
- * Auxiliary vector.
- */
- mutable Vector<double> h2;
-};
-
-inline
-void
-MGSmoother::set_steps(unsigned int i)
-{
- steps = i;
-}
-
-inline
-unsigned int
-MGSmoother::get_steps() const
-{
- return steps;
-}
-
-/*---------------------------- mg_smoother.h ---------------------------*/
-/* end of #ifndef __mg_smoother_H */
-#endif
-/*---------------------------- mg_smoother.h ---------------------------*/
-
+++ /dev/null
-// $Id$
-
-template<typename number>
-template<int dim>
-MGSmootherRelaxation<number>::
-MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<number> >& matrix,
- function_ptr relaxation,
- unsigned int steps,
- double omega)
- :
- MGSmoother(mg_dof, steps),
- matrix(&matrix),
- relaxation(relaxation),
- omega(omega)
-{}
-
-template<typename number>
-void
-MGSmootherRelaxation<number>::smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const
-{
- h1.reinit(u);
- h2.reinit(u);
- for(unsigned i=0;i<get_steps();++i)
- {
- (*matrix)[level].residual(h1, u, rhs);
- ((*matrix)[level].*relaxation)(h2,h1,omega);
- set_zero_interior_boundary(level,h2);
- u.add(h2);
- }
-}
-
+++ /dev/null
-/*---------------------------- multigrid.h ---------------------------*/
-/* $Id$ */
-#ifndef __multigrid_H
-#define __multigrid_H
-/*---------------------------- multigrid.h ---------------------------*/
-
-#include <base/subscriptor.h>
-#include <base/smartpointer.h>
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
-#include <lac/vector.h>
-#include <lac/mgbase.h>
-
-#include <vector>
-
-
-
-/**
- * Implementation of multigrid with matrices.
- * While #MGBase# was only an abstract framework for the v-cycle,
- * here we have the implementation of the pure virtual functions defined there.
- * Furthermore, level information is obtained from a triangulation object.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-template<int dim>
-class MG : public MGBase
-{
- public:
- /**
- * Constructor. Take the object
- * describing the multilevel
- * hierarchy of degrees of
- * freedom, an object describing
- * constraints to degrees of
- * freedom on the global grid
- * (for example for hanging
- * nodes), a set of matrices
- * which describe the equation
- * discretized on each level, and
- * an object mediating the
- * transfer of solutions between
- * different levels.
- *
- * This function already
- * initializes the vectors which
- * will be used later on in the
- * course of the
- * computations. You should
- * therefore create objects of
- * this type as late as possible.
- */
-//TODO: minlevel, maxlevel?
- MG(const MGDoFHandler<dim> &mg_dof_handler,
- ConstraintMatrix &constraints,
- const MGMatrix<SparseMatrix<double> > &level_matrices,
- const MGTransferBase &transfer,
- const unsigned int minlevel = 0,
- const unsigned int maxlevel = 1000000);
-
- /**
- * Transfer from a vector on the
- * global grid to vectors defined
- * on each of the levels
- * separately, i.a. an #MGVector#.
- */
- template<typename number>
- void copy_to_mg (const Vector<number> &src);
-
- /**
- * Transfer from multi-level vector to
- * normal vector.
- *
- * Copies data from active
- * portions of an MGVector into
- * the respective positions of a
- * #Vector<double>#. In order to
- * keep the result consistent,
- * constrained degrees of freedom are set to zero.
- */
- template<typename number>
- void copy_from_mg (Vector<number> &dst) const;
-
- /**
- * Negative #vmult# on a level.
-//TODO: what does this function do?
- * @see MGBase.
- */
- virtual void level_vmult (const unsigned int level,
- Vector<double> &dest,
- const Vector<double> &src,
- const Vector<double> &rhs);
-
- /**
- * Read-only access to level matrices.
- */
- const SparseMatrix<double>& get_matrix (const unsigned int level) const;
-
- private:
- /**
- * Associated #MGDoFHandler#.
- */
- SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
-
- /**
- * Matrices for each level.
- * The matrices are prepared by
- * the constructor of #MG# and can
- * be accessed for assembling.
- */
- SmartPointer<const MGMatrix<SparseMatrix<double> > > level_matrices;
-
- /**
- * Pointer to the object holding
- * constraints.
- */
- SmartPointer<const ConstraintMatrix> constraints;
-};
-
-
-
-
-
-/**
- * Implementation of the #MGTransferBase# interface for which the transfer
- * operations are prebuilt upon construction of the object of this class as
- * matrices. This is the fast way, since it only needs to build the operation
- * once by looping over all cells and storing the result in a matrix for
- * each level, but requires additional memory.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-class MGTransferPrebuilt : public MGTransferBase
-{
- public:
- /**
- * Destructor.
- */
- virtual ~MGTransferPrebuilt ();
-
- /**
- * Actually build the prolongation
- * matrices for each level.
- */
- template <int dim>
- void build_matrices (const MGDoFHandler<dim> &mg_dof);
-
- /**
- * Prolongate a vector from level
- * #to_level-1# to level
- * #to_level#. The previous
- * content of #dst# is
- * overwritten.
- *
- * #src# is assumed to be a vector with
- * as many elements as there are degrees
- * of freedom on the coarser level of
- * the two involved levels, while #src#
- * shall have as many elements as there
- * are degrees of freedom on the finer
- * level.
- */
- virtual void prolongate (const unsigned int to_level,
- Vector<double> &dst,
- const Vector<double> &src) const;
-
- /**
- * Restrict a vector from level
- * #from_level# to level
- * #from_level-1# and add this
- * restriction to
- * #dst#. Obviously, if the
- * refined region on level
- * #from_level# is smaller than
- * that on level #from_level-1#,
- * some degrees of freedom in
- * #dst# are not covered and will
- * not be altered. For the other
- * degress of freedom, the result
- * of the restriction is added.
- *
- * #src# is assumed to be a vector with
- * as many elements as there are degrees
- * of freedom on the finer level of
- * the two involved levels, while #src#
- * shall have as many elements as there
- * are degrees of freedom on the coarser
- * level.
- */
- virtual void restrict_and_add (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const;
-
- private:
-
- /**
- * Sparsity patterns for the
- * prolongation matrices on each
- * level.
- */
- vector<SparsityPattern> prolongation_sparsities;
-
- /**
- * The actual prolongation matrix.
- * column indices belong to the
- * dof indices of the mother cell,
- * i.e. the coarse level.
- * while row indices belong to the
- * child cell, i.e. the fine level.
- */
- vector<SparseMatrix<float> > prolongation_matrices;
-};
-
-
-
-
-/* --------------------------- inline functions --------------------- */
-
-
-template<int dim>
-inline
-const SparseMatrix<double>&
-MG<dim>::get_matrix(unsigned int level) const
-{
- Assert((level>=minlevel) && (level<maxlevel),
- ExcIndexRange(level, minlevel, maxlevel));
-
- return (*level_matrices)[level];
-}
-
-
-/*---------------------------- multigrid.h ---------------------------*/
-/* end of #ifndef __multigrid_H */
-#endif
-/*---------------------------- multigrid.h ---------------------------*/
+++ /dev/null
-/*---------------------------- multigrid.templates.h ---------------------------*/
-/* $Id$ */
-#ifndef __multigrid_templates_H
-#define __multigrid_templates_H
-/*---------------------------- multigrid.templates.h ---------------------------*/
-
-
-#include <dofs/dof_constraints.h>
-#include <numerics/multigrid.h>
-#include <algorithm>
-#include <fstream>
-
-
-
-
-template <int dim>
-template <typename number>
-void
-MG<dim>::copy_to_mg(const Vector<number>& src)
-{
- const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell;
- const unsigned int dofs_per_face = mg_dof_handler->get_fe().dofs_per_face;
-
- // set the elements of the vectors
- // on all levels to zero
- defect.clear();
-// constraints->condense(src);
-
- vector<int> global_dof_indices (dofs_per_cell);
- vector<int> level_dof_indices (dofs_per_cell);
- vector<int> level_face_indices (dofs_per_face);
-
- // traverse the grid top-down
- // (i.e. starting with the most
- // refined grid). this way, we can
- // always get that part of one
- // level of the output vector which
- // corresponds to a region which is
- // more refined, by restriction of
- // the respective vector on the
- // next finer level, which we then
- // already have built.
- for (int level=maxlevel; level>=static_cast<int>(minlevel); --level)
- {
- DoFHandler<dim>::active_cell_iterator
- global_cell = mg_dof_handler->DoFHandler<dim>::begin_active(level);
- MGDoFHandler<dim>::active_cell_iterator
- level_cell = mg_dof_handler->begin_active(level);
- const MGDoFHandler<dim>::active_cell_iterator
- level_end = mg_dof_handler->end_active(level);
-
- for (; level_cell!=level_end; ++level_cell, ++global_cell)
- {
- // get the dof numbers of
- // this cell for the global
- // and the level-wise
- // numbering
- global_cell->get_dof_indices(global_dof_indices);
- level_cell->get_mg_dof_indices (level_dof_indices);
-
- // transfer the global
- // defect in the vector
- // into the level-wise one
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- defect[level](level_dof_indices[i]) = src(global_dof_indices[i]);
-
-//TODO: what happens here?
- // Delete values on refinement edge
- for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
- {
- const MGDoFHandler<dim>::face_iterator face = level_cell->face(face_n);
- if (face->has_children())
- {
- face->get_mg_dof_indices(level_face_indices);
- for (unsigned int i=0; i<dofs_per_face; ++i)
- defect[level](level_face_indices[i]) = 0.;
- };
- };
- };
-
- // for that part of the level
- // which is further refined:
- // get the defect by
- // restriction of the defect on
- // one level higher
- if (static_cast<unsigned int>(level) < maxlevel)
- transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
- };
-};
-
-
-
-template <int dim>
-template <typename number>
-void
-MG<dim>::copy_from_mg(Vector<number> &dst) const
-{
- const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell;
-
- vector<int> global_dof_indices (dofs_per_cell);
- vector<int> level_dof_indices (dofs_per_cell);
-
- DoFHandler<dim>::active_cell_iterator
- global_cell = mg_dof_handler->DoFHandler<dim>::begin_active();
- MGDoFHandler<dim>::active_cell_iterator
- level_cell = mg_dof_handler->begin_active();
- const MGDoFHandler<dim>::active_cell_iterator
- endc = mg_dof_handler->end();
-
- // traverse all cells and copy the
- // data appropriately to the output
- // vector
- for (; level_cell != endc; ++level_cell, ++global_cell)
- {
- const unsigned int level = level_cell->level();
-
- // get the dof numbers of
- // this cell for the global
- // and the level-wise
- // numbering
- global_cell->get_dof_indices (global_dof_indices);
- level_cell->get_mg_dof_indices(level_dof_indices);
-
- // copy level-wise data to
- // global vector
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- dst(global_dof_indices[i]) = solution[level](level_dof_indices[i]);
- };
-
- // clear constrained nodes
- constraints->set_zero(dst);
-}
-
-
-
-
-/*---------------------------- multigrid.templates.h ---------------------------*/
-/* end of #ifndef __multigrid.templates_H */
-#endif
-/*---------------------------- multigrid.templates.h ---------------------------*/
#include <dofs/dof_accessor.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_handler.h>
-#include <dofs/mg_dof_handler.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_tools.h>
-#include <dofs/mg_dof_tools.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_tools.h>
#include <fe/fe.h>
#include <numerics/dof_renumbering.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_handler.h>
-#include <dofs/mg_dof_handler.h>
#include <dofs/dof_accessor.h>
-#include <dofs/mg_dof_accessor.h>
#include <dofs/dof_constraints.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <fe/fe.h>
#include <fe/fe_system.h>
#include <dofs/dof_tools.h>
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-
-#include <dofs/mg_dof_accessor.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/dof_levels.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_iterator.templates.h>
-#include <fe/fe.h>
-
-#include <lac/vector.h>
-
-
-
-
-
-
-
-
-/* ------------------------ MGDoFLineAccessor --------------------------- */
-
-template <int dim>
-MGDoFObjectAccessor<1, dim>::MGDoFObjectAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data) :
- MGDoFAccessor<dim> (local_data),
- MGDoFObjectAccessor_Inheritance<1,dim>::BaseClass(tria,level,index,local_data)
-{};
-
-
-
-template <int dim>
-int MGDoFObjectAccessor<1, dim>::mg_dof_index (const unsigned int i) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- // make sure a FE has been selected
- // and enough room was reserved
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (i<dof_handler->get_fe().dofs_per_line,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_line));
-
- return mg_dof_handler->mg_levels[present_level]
- ->line_dofs[present_index*dof_handler->get_fe().dofs_per_line+i];
-};
-
-
-
-
-
-template <int dim>
-void MGDoFObjectAccessor<1, dim>::set_mg_dof_index (const unsigned int i,
- const int index) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- // make sure a FE has been selected
- // and enough room was reserved
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (i<dof_handler->get_fe().dofs_per_line,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_line));
-
- mg_dof_handler->mg_levels[present_level]
- ->line_dofs[present_index*dof_handler->get_fe().dofs_per_line+i] = index;
-};
-
-
-
-
-template <int dim>
-int MGDoFObjectAccessor<1, dim>::mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (vertex<2, ExcIndexRange (i,0,2));
- Assert (i<dof_handler->get_fe().dofs_per_vertex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex));
-
- return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
- .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
-};
-
-
-
-template <int dim>
-void MGDoFObjectAccessor<1, dim>::set_mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i,
- const int index) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (vertex<2, ExcIndexRange (i,0,2));
- Assert (i<dof_handler->get_fe().dofs_per_vertex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex));
-
- mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
- .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
-};
-
-
-
-template <int dim>
-void
-MGDoFObjectAccessor<1, dim>::get_mg_dof_indices (vector<int> &dof_indices) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (dof_indices.size() == (2*dof_handler->get_fe().dofs_per_vertex +
- dof_handler->get_fe().dofs_per_line),
- DoFAccessor<dim>::ExcVectorDoesNotMatch());
-
- const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
- dofs_per_line = dof_handler->get_fe().dofs_per_line;
- vector<int>::iterator next = dof_indices.begin();
- for (unsigned int vertex=0; vertex<2; ++vertex)
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- *next++ = mg_vertex_dof_index(vertex,d);
- for (unsigned int d=0; d<dofs_per_line; ++d)
- *next++ = mg_dof_index(d);
-
- Assert (next == dof_indices.end(),
- ExcInternalError());
-};
-
-
-
-template <int dim>
-template <typename number>
-void
-MGDoFObjectAccessor<1,dim>::get_mg_dof_values (const Vector<number> &values,
- Vector<number> &dof_values) const
-{
- Assert (dof_handler != 0, DoFAccessor<1>::DoFAccessor<1>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<1>::DoFAccessor<1>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<1>::DoFAccessor<1>::ExcInvalidObject());
- Assert (dof_values.size() == dof_handler->get_fe().dofs_per_cell,
- DoFAccessor<1>::ExcVectorDoesNotMatch());
- Assert (values.size() == dof_handler->n_dofs(),
- DoFAccessor<1>::ExcVectorDoesNotMatch());
-
- const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
- dofs_per_line = dof_handler->get_fe().dofs_per_line;
- vector<number>::iterator next_dof_value=dof_values.begin();
- for (unsigned int vertex=0; vertex<2; ++vertex)
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
- for (unsigned int d=0; d<dofs_per_line; ++d)
- *next_dof_value++ = values(mg_dof_index(d));
-
- Assert (next_dof_value == dof_values.end(),
- ExcInternalError());
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
-MGDoFObjectAccessor<1, dim>::child (const unsigned int i) const
-{
- TriaIterator<dim,MGDoFObjectAccessor<1, dim> > q (tria,
- present_level+1,
- child_index (i),
- mg_dof_handler);
-
-#ifdef DEBUG
- if (q.state() != past_the_end)
- Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
-#endif
- return q;
-};
-
-
-
-template <int dim>
-void
-MGDoFObjectAccessor<1, dim>::copy_from (const MGDoFObjectAccessor<1, dim> &a)
-{
- DoFObjectAccessor<1, dim>::copy_from (a);
- set_mg_dof_handler (a.mg_dof_handler);
-};
-
-
-
-
-/* ------------------------ MGDoFQuadAccessor --------------------------- */
-
-template <int dim>
-MGDoFObjectAccessor<2, dim>::MGDoFObjectAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data) :
- MGDoFAccessor<dim> (local_data),
- MGDoFObjectAccessor_Inheritance<2,dim>::BaseClass(tria,level,index,local_data)
-{};
-
-
-
-template <int dim>
-inline
-int MGDoFObjectAccessor<2, dim>::mg_dof_index (const unsigned int i) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- // make sure a FE has been selected
- // and enough room was reserved
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (i<dof_handler->get_fe().dofs_per_quad,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_quad));
-
- return mg_dof_handler->mg_levels[present_level]
- ->quad_dofs[present_index*dof_handler->get_fe().dofs_per_quad+i];
-};
-
-
-
-template <int dim>
-void MGDoFObjectAccessor<2, dim>::set_mg_dof_index (const unsigned int i,
- const int index) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- // make sure a FE has been selected
- // and enough room was reserved
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (i<dof_handler->get_fe().dofs_per_quad,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_quad));
-
- mg_dof_handler->mg_levels[present_level]
- ->quad_dofs[present_index*dof_handler->get_fe().dofs_per_quad+i] = index;
-};
-
-
-
-template <int dim>
-inline
-int MGDoFObjectAccessor<2, dim>::mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (vertex<4, ExcIndexRange (i,0,4));
- Assert (i<dof_handler->get_fe().dofs_per_vertex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex));
-
- return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
- .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
-};
-
-
-
-template <int dim>
-void MGDoFObjectAccessor<2, dim>::set_mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i,
- const int index) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (vertex<4, ExcIndexRange (i,0,4));
- Assert (i<dof_handler->get_fe().dofs_per_vertex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex));
-
- mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
- .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
-};
-
-
-
-template <int dim>
-void
-MGDoFObjectAccessor<2, dim>::get_mg_dof_indices (vector<int> &dof_indices) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (dof_indices.size() == (4*dof_handler->get_fe().dofs_per_vertex +
- 4*dof_handler->get_fe().dofs_per_line +
- dof_handler->get_fe().dofs_per_quad),
- DoFAccessor<2>::ExcVectorDoesNotMatch());
-
- const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
- dofs_per_line = dof_handler->get_fe().dofs_per_line,
- dofs_per_quad = dof_handler->get_fe().dofs_per_quad;
- vector<int>::iterator next = dof_indices.begin();
- for (unsigned int vertex=0; vertex<4; ++vertex)
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- *next++ = mg_vertex_dof_index(vertex,d);
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int d=0; d<dofs_per_line; ++d)
- *next++ = this->line(line)->mg_dof_index(d);
- for (unsigned int d=0; d<dofs_per_quad; ++d)
- *next++ = mg_dof_index(d);
-
- Assert (next == dof_indices.end(),
- ExcInternalError());
-};
-
-
-
-template <int dim>
-template <typename number>
-void
-MGDoFObjectAccessor<2,dim>::get_mg_dof_values (const Vector<number> &values,
- Vector<number> &dof_values) const
-{
- Assert (dof_handler != 0, DoFAccessor<2>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<2>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<2>::ExcInvalidObject());
- Assert (dof_values.size() == dof_handler->get_fe().dofs_per_cell,
- DoFAccessor<2>::ExcVectorDoesNotMatch());
- Assert (values.size() == mg_dof_handler->n_dofs(present_level),
- DoFAccessor<2>::ExcVectorDoesNotMatch());
-
- const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
- dofs_per_line = dof_handler->get_fe().dofs_per_line,
- dofs_per_quad = dof_handler->get_fe().dofs_per_quad;
- vector<number>::iterator next_dof_value=dof_values.begin();
- for (unsigned int vertex=0; vertex<4; ++vertex)
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int d=0; d<dofs_per_line; ++d)
- *next_dof_value++ = values(this->line(line)->mg_dof_index(d));
- for (unsigned int d=0; d<dofs_per_quad; ++d)
- *next_dof_value++ = values(mg_dof_index(d));
-
- Assert (next_dof_value == dof_values.end(),
- ExcInternalError());
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
-MGDoFObjectAccessor<2, dim>::line (const unsigned int i) const
-{
- Assert (i<4, ExcIndexRange (i, 0, 4));
-
- return TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
- (
- tria,
- present_level,
- line_index (i),
- mg_dof_handler
- );
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFObjectAccessor<2, dim> >
-MGDoFObjectAccessor<2, dim>::child (const unsigned int i) const
-{
- TriaIterator<dim,MGDoFObjectAccessor<2, dim> > q (tria,
- present_level+1,
- child_index (i),
- mg_dof_handler);
-
-#ifdef DEBUG
- if (q.state() != past_the_end)
- Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
-#endif
- return q;
-};
-
-
-
-template <int dim>
-void
-MGDoFObjectAccessor<2, dim>::copy_from (const MGDoFObjectAccessor<2, dim> &a)
-{
- DoFObjectAccessor<2, dim>::copy_from (a);
- set_mg_dof_handler (a.mg_dof_handler);
-};
-
-
-
-
-/* ------------------------ MGDoFHexAccessor --------------------------- */
-
-template <int dim>
-MGDoFObjectAccessor<3, dim>::MGDoFObjectAccessor (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AccessorData *local_data) :
- MGDoFAccessor<dim> (local_data),
- MGDoFObjectAccessor_Inheritance<3,dim>::BaseClass(tria,level,index,local_data)
-{};
-
-
-
-template <int dim>
-inline
-int MGDoFObjectAccessor<3, dim>::mg_dof_index (const unsigned int i) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- // make sure a FE has been selected
- // and enough room was reserved
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (i<dof_handler->get_fe().dofs_per_hex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_hex));
-
- return mg_dof_handler->mg_levels[present_level]
- ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i];
-};
-
-
-
-template <int dim>
-void MGDoFObjectAccessor<3, dim>::set_mg_dof_index (const unsigned int i,
- const int index) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- // make sure a FE has been selected
- // and enough room was reserved
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (i<dof_handler->get_fe().dofs_per_hex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_hex));
-
- mg_dof_handler->mg_levels[present_level]
- ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i] = index;
-};
-
-
-
-template <int dim>
-inline
-int MGDoFObjectAccessor<3, dim>::mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (vertex<8, ExcIndexRange (i,0,8));
- Assert (i<dof_handler->get_fe().dofs_per_vertex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex));
-
- return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
- .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
-};
-
-
-
-template <int dim>
-void MGDoFObjectAccessor<3, dim>::set_mg_vertex_dof_index (const unsigned int vertex,
- const unsigned int i,
- const int index) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (vertex<4, ExcIndexRange (i,0,4));
- Assert (i<dof_handler->get_fe().dofs_per_vertex,
- ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex));
-
- mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
- .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
-};
-
-
-
-template <int dim>
-void
-MGDoFObjectAccessor<3, dim>::get_mg_dof_indices (vector<int> &dof_indices) const
-{
- Assert (dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<dim>::ExcInvalidObject());
- Assert (dof_indices.size() == (8*dof_handler->get_fe().dofs_per_vertex +
- 12*dof_handler->get_fe().dofs_per_line +
- 6*dof_handler->get_fe().dofs_per_quad +
- dof_handler->get_fe().dofs_per_hex),
- DoFAccessor<3>::ExcVectorDoesNotMatch());
-
- const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
- dofs_per_line = dof_handler->get_fe().dofs_per_line,
- dofs_per_quad = dof_handler->get_fe().dofs_per_quad,
- dofs_per_hex = dof_handler->get_fe().dofs_per_hex;
- vector<int>::iterator next = dof_indices.begin();
- for (unsigned int vertex=0; vertex<8; ++vertex)
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- *next++ = mg_vertex_dof_index(vertex,d);
- for (unsigned int line=0; line<12; ++line)
- for (unsigned int d=0; d<dofs_per_line; ++d)
- *next++ = this->line(line)->mg_dof_index(d);
- for (unsigned int quad=0; quad<12; ++quad)
- for (unsigned int d=0; d<dofs_per_quad; ++d)
- *next++ = this->quad(quad)->mg_dof_index(d);
- for (unsigned int d=0; d<dofs_per_hex; ++d)
- *next++ = mg_dof_index(d);
-
- Assert (next == dof_indices.end(),
- ExcInternalError());
-};
-
-
-
-template <int dim>
-template <typename number>
-void
-MGDoFObjectAccessor<3,dim>::get_mg_dof_values (const Vector<number> &values,
- Vector<number> &dof_values) const
-{
- Assert (dof_handler != 0, DoFAccessor<2>::ExcInvalidObject());
- Assert (mg_dof_handler != 0, DoFAccessor<2>::ExcInvalidObject());
- Assert (&dof_handler->get_fe() != 0, DoFAccessor<2>::ExcInvalidObject());
- Assert (dof_values.size() == dof_handler->get_fe().dofs_per_cell,
- DoFAccessor<3>::ExcVectorDoesNotMatch());
- Assert (values.size() == mg_dof_handler->n_dofs(present_level),
- DoFAccessor<3>::ExcVectorDoesNotMatch());
-
- const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
- dofs_per_line = dof_handler->get_fe().dofs_per_line,
- dofs_per_quad = dof_handler->get_fe().dofs_per_quad,
- dofs_per_hex = dof_handler->get_fe().dofs_per_hex;
- vector<number>::iterator next_dof_value=dof_values.begin();
- for (unsigned int vertex=0; vertex<8; ++vertex)
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
- for (unsigned int line=0; line<12; ++line)
- for (unsigned int d=0; d<dofs_per_line; ++d)
- *next_dof_value++ = values(this->line(line)->mg_dof_index(d));
- for (unsigned int quad=0; quad<6; ++quad)
- for (unsigned int d=0; d<dofs_per_quad; ++d)
- *next_dof_value++ = values(this->quad(quad)->mg_dof_index(d));
- for (unsigned int d=0; d<dofs_per_hex; ++d)
- *next_dof_value++ = values(dof_index(d));
-
- Assert (next_dof_value == dof_values.end(),
- ExcInternalError());
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
-MGDoFObjectAccessor<3, dim>::line (const unsigned int i) const {
- Assert (i<12, ExcIndexRange (i, 0, 12));
-
- return TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
- (
- tria,
- present_level,
- line_index (i),
- mg_dof_handler
- );
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFObjectAccessor<2, dim> >
-MGDoFObjectAccessor<3, dim>::quad (const unsigned int i) const {
- Assert (i<12, ExcIndexRange (i, 0, 6));
-
- return TriaIterator<dim,MGDoFObjectAccessor<2, dim> >
- (
- tria,
- present_level,
- quad_index (i),
- mg_dof_handler
- );
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFObjectAccessor<3, dim> >
-MGDoFObjectAccessor<3, dim>::child (const unsigned int i) const {
- TriaIterator<dim,MGDoFObjectAccessor<3, dim> > q (tria,
- present_level+1,
- child_index (i),
- mg_dof_handler);
-
-#ifdef DEBUG
- if (q.state() != past_the_end)
- Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
-#endif
- return q;
-};
-
-
-
-template <int dim>
-void
-MGDoFObjectAccessor<3, dim>::copy_from (const MGDoFObjectAccessor<3, dim> &a) {
- DoFObjectAccessor<3, dim>::copy_from (a);
- set_mg_dof_handler (a.mg_dof_handler);
-};
-
-
-
-
-
-/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
-
-template <int dim>
-TriaIterator<dim,MGDoFCellAccessor<dim> >
-MGDoFCellAccessor<dim>::neighbor (const unsigned int i) const {
- TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
- neighbor_level (i),
- neighbor_index (i),
- mg_dof_handler);
-
-#ifdef DEBUG
- if (q.state() != past_the_end)
- Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsNeighbor());
-#endif
- return q;
-};
-
-
-
-template <int dim>
-TriaIterator<dim,MGDoFCellAccessor<dim> >
-MGDoFCellAccessor<dim>::child (const unsigned int i) const {
- TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
- present_level+1,
- child_index (i),
- mg_dof_handler);
-
-#ifdef DEBUG
- if (q.state() != past_the_end)
- Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
-#endif
- return q;
-};
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-MGDoFCellAccessor<1>::face_iterator
-MGDoFCellAccessor<1>::face (const unsigned int) const
-{
- Assert (false, ExcNotUsefulForThisDimension());
- return face_iterator();
-};
-
-#endif
-
-
-
-#if deal_II_dimension == 2
-
-template <>
-MGDoFCellAccessor<2>::face_iterator
-MGDoFCellAccessor<2>::face (const unsigned int i) const
-{
- return line(i);
-};
-
-#endif
-
-
-
-#if deal_II_dimension == 3
-
-template <>
-MGDoFCellAccessor<3>::face_iterator
-MGDoFCellAccessor<3>::face (const unsigned int i) const
-{
- return quad(i);
-};
-
-#endif
-
-
-
-
-
-// explicit instantiations
-
-template
-void
-MGDoFObjectAccessor<1,deal_II_dimension>::
-get_mg_dof_values (const Vector<double> &values,
- Vector<double> &dof_values) const;
-
-template
-void
-MGDoFObjectAccessor<1,deal_II_dimension>::
-get_mg_dof_values (const Vector<float> &values,
- Vector<float> &dof_values) const;
-
-
-
-#if deal_II_dimension >= 2
-
-template
-void
-MGDoFObjectAccessor<2,deal_II_dimension>::
-get_mg_dof_values (const Vector<double> &values,
- Vector<double> &dof_values) const;
-
-template
-void
-MGDoFObjectAccessor<2,deal_II_dimension>::
-get_mg_dof_values (const Vector<float> &values,
- Vector<float> &dof_values) const;
-
-#endif
-
-
-#if deal_II_dimension >= 3
-
-template
-void
-MGDoFObjectAccessor<3,deal_II_dimension>::
-get_mg_dof_values (const Vector<double> &values,
- Vector<double> &dof_values) const;
-
-template
-void
-MGDoFObjectAccessor<3,deal_II_dimension>::
-get_mg_dof_values (const Vector<float> &values,
- Vector<float> &dof_values) const;
-
-#endif
-
-
-#if deal_II_dimension == 1
-template class MGDoFObjectAccessor<1, 1>;
-template class MGDoFCellAccessor<1>;
-
-template class TriaRawIterator<1,MGDoFCellAccessor<1> >;
-template class TriaIterator<1,MGDoFCellAccessor<1> >;
-template class TriaActiveIterator<1,MGDoFCellAccessor<1> >;
-#endif
-
-#if deal_II_dimension == 2
-template class MGDoFObjectAccessor<1, 2>;
-template class MGDoFObjectAccessor<2, 2>;
-template class MGDoFCellAccessor<2>;
-
-template class TriaRawIterator <2,MGDoFObjectAccessor<1, 2> >;
-template class TriaIterator <2,MGDoFObjectAccessor<1, 2> >;
-template class TriaActiveIterator<2,MGDoFObjectAccessor<1, 2> >;
-
-template class TriaRawIterator <2,MGDoFCellAccessor<2> >;
-template class TriaIterator <2,MGDoFCellAccessor<2> >;
-template class TriaActiveIterator<2,MGDoFCellAccessor<2> >;
-#endif
-
-
-#if deal_II_dimension == 3
-template class MGDoFObjectAccessor<1, 3>;
-template class MGDoFObjectAccessor<2, 3>;
-template class MGDoFObjectAccessor<3, 3>;
-template class MGDoFCellAccessor<3>;
-
-template class TriaRawIterator <3,MGDoFObjectAccessor<1, 3> >;
-template class TriaIterator <3,MGDoFObjectAccessor<1, 3> >;
-template class TriaActiveIterator<3,MGDoFObjectAccessor<1, 3> >;
-
-template class TriaRawIterator <3,MGDoFObjectAccessor<2, 3> >;
-template class TriaIterator <3,MGDoFObjectAccessor<2, 3> >;
-template class TriaActiveIterator<3,MGDoFObjectAccessor<2, 3> >;
-
-template class TriaRawIterator <3,MGDoFCellAccessor<3> >;
-template class TriaIterator <3,MGDoFCellAccessor<3> >;
-template class TriaActiveIterator<3,MGDoFCellAccessor<3> >;
-#endif
-
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-#include <dofs/mg_dof_handler.h>
-#include <dofs/dof_levels.h>
-#include <dofs/mg_dof_accessor.h>
-#include <dofs/dof_constraints.h>
-#include <grid/tria_levels.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria.h>
-#include <grid/geometry_info.h>
-#include <fe/fe.h>
-
-#include <algorithm>
-
-
-
-/* ------------------------ MGVertexDoFs ----------------------------------- */
-
-template <int dim>
-MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs () :
- coarsest_level (1<<30),
- finest_level (0),
- indices (0)
-{};
-
-
-
-template <int dim>
-void MGDoFHandler<dim>::MGVertexDoFs::init (const unsigned int cl,
- const unsigned int fl,
- const unsigned int dofs_per_vertex) {
- Assert (indices == 0, ExcInternalError());
-
- coarsest_level = cl;
- finest_level = fl;
-
- const unsigned int n_levels = finest_level-coarsest_level + 1;
-
- indices = new int[n_levels * dofs_per_vertex];
- Assert (indices != 0, ExcNoMemory ());
-
- for (unsigned int i=0; i<n_levels*dofs_per_vertex; ++i)
- indices[i] = -1;
-};
-
-
-
-template <int dim>
-MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
- delete[] indices;
-};
-
-
-
-template <int dim>
-inline
-unsigned int MGDoFHandler<dim>::MGVertexDoFs::get_coarsest_level () const {
- return coarsest_level;
-};
-
-
-
-template <int dim>
-inline
-unsigned int MGDoFHandler<dim>::MGVertexDoFs::get_finest_level () const {
- return finest_level;
-};
-
-
-
-
-
-
-/* ------------------------ MGDoFHandler ------------------------------------- */
-
-template <int dim>
-MGDoFHandler<dim>::MGDoFHandler (Triangulation<dim> &tria) :
- DoFHandler<dim> (tria)
-{};
-
-
-
-template <int dim>
-MGDoFHandler<dim>::~MGDoFHandler () {};
-
-
-
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-MGDoFHandler<1>::raw_cell_iterator
-MGDoFHandler<1>::begin_raw (const unsigned int level) const {
- return begin_raw_line (level);
-};
-
-
-
-template <>
-MGDoFHandler<1>::cell_iterator
-MGDoFHandler<1>::begin (const unsigned int level) const {
- return begin_line (level);
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_cell_iterator
-MGDoFHandler<1>::begin_active (const unsigned int level) const {
- return begin_active_line (level);
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_cell_iterator
-MGDoFHandler<1>::end () const {
- return end_line ();
-};
-
-
-
-
-template <>
-MGDoFHandler<1>::raw_cell_iterator
-MGDoFHandler<1>::last_raw () const {
- return last_raw_line ();
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_cell_iterator
-MGDoFHandler<1>::last_raw (const unsigned int level) const {
- return last_raw_line (level);
-};
-
-
-
-template <>
-MGDoFHandler<1>::cell_iterator
-MGDoFHandler<1>::last () const {
- return last_line ();
-};
-
-
-
-template <>
-MGDoFHandler<1>::cell_iterator
-MGDoFHandler<1>::last (const unsigned int level) const {
- return last_line (level);
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_cell_iterator
-MGDoFHandler<1>::last_active () const {
- return last_active_line ();
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_cell_iterator
-MGDoFHandler<1>::last_active (const unsigned int level) const {
- return last_active_line (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::raw_face_iterator
-MGDoFHandler<1>::begin_raw_face (const unsigned int) const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::face_iterator
-MGDoFHandler<1>::begin_face (const unsigned int) const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::active_face_iterator
-MGDoFHandler<1>::begin_active_face (const unsigned int) const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::raw_face_iterator
-MGDoFHandler<1>::end_face () const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::raw_face_iterator
-MGDoFHandler<1>::last_raw_face () const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::raw_face_iterator
-MGDoFHandler<1>::last_raw_face (const unsigned int) const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::face_iterator
-MGDoFHandler<1>::last_face () const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::face_iterator
-MGDoFHandler<1>::last_face (const unsigned int) const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::active_face_iterator
-MGDoFHandler<1>::last_active_face () const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-template <>
-MGDoFDimensionInfo<1>::active_face_iterator
-MGDoFHandler<1>::last_active_face (const unsigned int) const {
- Assert (false, ExcFunctionNotUseful());
- return 0;
-};
-
-
-
-
-template <>
-MGDoFHandler<1>::raw_quad_iterator
-MGDoFHandler<1>::begin_raw_quad (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::quad_iterator
-MGDoFHandler<1>::begin_quad (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_quad_iterator
-MGDoFHandler<1>::begin_active_quad (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_quad_iterator
-MGDoFHandler<1>::end_quad () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_quad_iterator
-MGDoFHandler<1>::last_raw_quad (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-template <>
-MGDoFHandler<1>::quad_iterator
-MGDoFHandler<1>::last_quad (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_quad_iterator
-MGDoFHandler<1>::last_active_quad (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_quad_iterator
-MGDoFHandler<1>::last_raw_quad () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::quad_iterator
-MGDoFHandler<1>::last_quad () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_quad_iterator
-MGDoFHandler<1>::last_active_quad () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-template <>
-MGDoFHandler<1>::raw_hex_iterator
-MGDoFHandler<1>::begin_raw_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::hex_iterator
-MGDoFHandler<1>::begin_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_hex_iterator
-MGDoFHandler<1>::begin_active_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_hex_iterator
-MGDoFHandler<1>::end_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_hex_iterator
-MGDoFHandler<1>::last_raw_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-template <>
-MGDoFHandler<1>::hex_iterator
-MGDoFHandler<1>::last_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_hex_iterator
-MGDoFHandler<1>::last_active_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::raw_hex_iterator
-MGDoFHandler<1>::last_raw_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::hex_iterator
-MGDoFHandler<1>::last_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<1>::active_hex_iterator
-MGDoFHandler<1>::last_active_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-#endif
-
-
-
-#if deal_II_dimension == 2
-
-template <>
-MGDoFHandler<2>::raw_cell_iterator
-MGDoFHandler<2>::begin_raw (const unsigned int level) const {
- return begin_raw_quad (level);
-};
-
-
-
-template <>
-MGDoFHandler<2>::cell_iterator
-MGDoFHandler<2>::begin (const unsigned int level) const {
- return begin_quad (level);
-};
-
-
-
-template <>
-MGDoFHandler<2>::active_cell_iterator
-MGDoFHandler<2>::begin_active (const unsigned int level) const {
- return begin_active_quad (level);
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_cell_iterator
-MGDoFHandler<2>::end () const {
- return end_quad ();
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_cell_iterator
-MGDoFHandler<2>::last_raw () const {
- return last_raw_quad ();
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_cell_iterator
-MGDoFHandler<2>::last_raw (const unsigned int level) const {
- return last_raw_quad (level);
-};
-
-
-
-template <>
-MGDoFHandler<2>::cell_iterator
-MGDoFHandler<2>::last () const {
- return last_quad ();
-};
-
-
-
-template <>
-MGDoFHandler<2>::cell_iterator
-MGDoFHandler<2>::last (const unsigned int level) const {
- return last_quad (level);
-};
-
-
-
-template <>
-MGDoFHandler<2>::active_cell_iterator
-MGDoFHandler<2>::last_active () const {
- return last_active_quad ();
-};
-
-
-
-template <>
-MGDoFHandler<2>::active_cell_iterator
-MGDoFHandler<2>::last_active (const unsigned int level) const {
- return last_active_quad (level);
-};
-
-
-template <>
-MGDoFDimensionInfo<2>::raw_face_iterator
-MGDoFHandler<2>::begin_raw_face (const unsigned int level) const {
- return begin_raw_line (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::face_iterator
-MGDoFHandler<2>::begin_face (const unsigned int level) const {
- return begin_line (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::active_face_iterator
-MGDoFHandler<2>::begin_active_face (const unsigned int level) const {
- return begin_active_line (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::raw_face_iterator
-MGDoFHandler<2>::end_face () const {
- return end_line ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::raw_face_iterator
-MGDoFHandler<2>::last_raw_face () const {
- return last_raw_line ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::raw_face_iterator
-MGDoFHandler<2>::last_raw_face (const unsigned int level) const {
- return last_raw_line (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::face_iterator
-MGDoFHandler<2>::last_face () const {
- return last_line ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::face_iterator
-MGDoFHandler<2>::last_face (const unsigned int level) const {
- return last_line (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::active_face_iterator
-MGDoFHandler<2>::last_active_face () const {
- return last_active_line ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<2>::active_face_iterator
-MGDoFHandler<2>::last_active_face (const unsigned int level) const {
- return last_active_line (level);
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_hex_iterator
-MGDoFHandler<2>::begin_raw_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::hex_iterator
-MGDoFHandler<2>::begin_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::active_hex_iterator
-MGDoFHandler<2>::begin_active_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_hex_iterator
-MGDoFHandler<2>::end_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_hex_iterator
-MGDoFHandler<2>::last_raw_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-template <>
-MGDoFHandler<2>::hex_iterator
-MGDoFHandler<2>::last_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::active_hex_iterator
-MGDoFHandler<2>::last_active_hex (const unsigned int) const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::raw_hex_iterator
-MGDoFHandler<2>::last_raw_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::hex_iterator
-MGDoFHandler<2>::last_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <>
-MGDoFHandler<2>::active_hex_iterator
-MGDoFHandler<2>::last_active_hex () const {
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-#endif
-
-
-
-
-#if deal_II_dimension == 3
-
-template <>
-MGDoFHandler<3>::raw_cell_iterator
-MGDoFHandler<3>::begin_raw (const unsigned int level) const {
- return begin_raw_hex (level);
-};
-
-
-
-template <>
-MGDoFHandler<3>::cell_iterator
-MGDoFHandler<3>::begin (const unsigned int level) const {
- return begin_hex (level);
-};
-
-
-
-template <>
-MGDoFHandler<3>::active_cell_iterator
-MGDoFHandler<3>::begin_active (const unsigned int level) const {
- return begin_active_hex (level);
-};
-
-
-
-template <>
-MGDoFHandler<3>::raw_cell_iterator
-MGDoFHandler<3>::end () const {
- return end_hex ();
-};
-
-
-
-template <>
-MGDoFHandler<3>::raw_cell_iterator
-MGDoFHandler<3>::last_raw () const {
- return last_raw_hex ();
-};
-
-
-
-template <>
-MGDoFHandler<3>::raw_cell_iterator
-MGDoFHandler<3>::last_raw (const unsigned int level) const {
- return last_raw_hex (level);
-};
-
-
-
-template <>
-MGDoFHandler<3>::cell_iterator
-MGDoFHandler<3>::last () const {
- return last_hex ();
-};
-
-
-
-template <>
-MGDoFHandler<3>::cell_iterator
-MGDoFHandler<3>::last (const unsigned int level) const {
- return last_hex (level);
-};
-
-
-
-template <>
-MGDoFHandler<3>::active_cell_iterator
-MGDoFHandler<3>::last_active () const {
- return last_active_hex ();
-};
-
-
-
-template <>
-MGDoFHandler<3>::active_cell_iterator
-MGDoFHandler<3>::last_active (const unsigned int level) const {
- return last_active_hex (level);
-};
-
-
-template <>
-MGDoFDimensionInfo<3>::raw_face_iterator
-MGDoFHandler<3>::begin_raw_face (const unsigned int level) const {
- return begin_raw_quad (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::face_iterator
-MGDoFHandler<3>::begin_face (const unsigned int level) const {
- return begin_quad (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::active_face_iterator
-MGDoFHandler<3>::begin_active_face (const unsigned int level) const {
- return begin_active_quad (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::raw_face_iterator
-MGDoFHandler<3>::end_face () const {
- return end_quad ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::raw_face_iterator
-MGDoFHandler<3>::last_raw_face () const {
- return last_raw_quad ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::raw_face_iterator
-MGDoFHandler<3>::last_raw_face (const unsigned int level) const {
- return last_raw_quad (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::face_iterator
-MGDoFHandler<3>::last_face () const {
- return last_quad ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::face_iterator
-MGDoFHandler<3>::last_face (const unsigned int level) const {
- return last_quad (level);
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::active_face_iterator
-MGDoFHandler<3>::last_active_face () const {
- return last_active_quad ();
-};
-
-
-
-template <>
-MGDoFDimensionInfo<3>::active_face_iterator
-MGDoFHandler<3>::last_active_face (const unsigned int level) const {
- return last_active_quad (level);
-};
-
-#endif
-
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_line_iterator
-MGDoFHandler<dim>::begin_raw_line (const unsigned int level) const {
- return raw_line_iterator (tria,
- tria->begin_raw_line(level)->level(),
- tria->begin_raw_line(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::line_iterator
-MGDoFHandler<dim>::begin_line (const unsigned int level) const {
- return line_iterator (tria,
- tria->begin_line(level)->level(),
- tria->begin_line(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_line_iterator
-MGDoFHandler<dim>::begin_active_line (const unsigned int level) const {
- return active_line_iterator (tria,
- tria->begin_active_line(level)->level(),
- tria->begin_active_line(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_quad_iterator
-MGDoFHandler<dim>::begin_raw_quad (const unsigned int level) const {
- return raw_quad_iterator (tria,
- tria->begin_raw_quad(level)->level(),
- tria->begin_raw_quad(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::quad_iterator
-MGDoFHandler<dim>::begin_quad (const unsigned int level) const {
- return quad_iterator (tria,
- tria->begin_quad(level)->level(),
- tria->begin_quad(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_quad_iterator
-MGDoFHandler<dim>::begin_active_quad (const unsigned int level) const {
- return active_quad_iterator (tria,
- tria->begin_active_quad(level)->level(),
- tria->begin_active_quad(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_hex_iterator
-MGDoFHandler<dim>::begin_raw_hex (const unsigned int level) const {
- return raw_hex_iterator (tria,
- tria->begin_raw_hex(level)->level(),
- tria->begin_raw_hex(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::hex_iterator
-MGDoFHandler<dim>::begin_hex (const unsigned int level) const {
- return hex_iterator (tria,
- tria->begin_hex(level)->level(),
- tria->begin_hex(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_hex_iterator
-MGDoFHandler<dim>::begin_active_hex (const unsigned int level) const {
- return active_hex_iterator (tria,
- tria->begin_active_hex(level)->level(),
- tria->begin_active_hex(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_line_iterator
-MGDoFHandler<dim>::end_line () const {
- return raw_line_iterator (tria, -1, -1, this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_quad_iterator
-MGDoFHandler<dim>::end_quad () const {
- return raw_quad_iterator (tria, -1, -1, this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_hex_iterator
-MGDoFHandler<dim>::end_hex () const {
- return raw_hex_iterator (tria, -1, -1, this);
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_cell_iterator
-MGDoFHandler<dim>::end_raw (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- end() :
- begin_raw (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::cell_iterator
-MGDoFHandler<dim>::end (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- cell_iterator(end()) :
- begin (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::active_cell_iterator
-MGDoFHandler<dim>::end_active (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- active_cell_iterator(end()) :
- begin_active (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_face_iterator
-MGDoFHandler<dim>::end_raw_face (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- end_face() :
- begin_raw_face (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::face_iterator
-MGDoFHandler<dim>::end_face (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- face_iterator(end_face()) :
- begin_face (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::active_face_iterator
-MGDoFHandler<dim>::end_active_face (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- active_face_iterator(end_face()) :
- begin_active_face (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_line_iterator
-MGDoFHandler<dim>::end_raw_line (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- end_line() :
- begin_raw_line (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::line_iterator
-MGDoFHandler<dim>::end_line (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- line_iterator(end_line()) :
- begin_line (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::active_line_iterator
-MGDoFHandler<dim>::end_active_line (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- active_line_iterator(end_line()) :
- begin_active_line (level+1));
-};
-
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_quad_iterator
-MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- end_quad() :
- begin_raw_quad (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::quad_iterator
-MGDoFHandler<dim>::end_quad (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- quad_iterator(end_quad()) :
- begin_quad (level+1));
-};
-
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::active_quad_iterator
-MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- active_quad_iterator(end_quad()) :
- begin_active_quad (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_hex_iterator
-MGDoFHandler<dim>::end_raw_hex (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- end_hex() :
- begin_raw_hex (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::hex_iterator
-MGDoFHandler<dim>::end_hex (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- hex_iterator(end_hex()) :
- begin_hex (level+1));
-};
-
-
-
-
-template <int dim>
-typename MGDoFDimensionInfo<dim>::active_hex_iterator
-MGDoFHandler<dim>::end_active_hex (const unsigned int level) const {
- return (level == mg_levels.size()-1 ?
- active_hex_iterator(end_hex()) :
- begin_active_hex (level+1));
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_line_iterator
-MGDoFHandler<dim>::last_raw_line (const unsigned int level) const {
- return raw_line_iterator (tria,
- tria->last_raw_line(level)->level(),
- tria->last_raw_line(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::line_iterator
-MGDoFHandler<dim>::last_line (const unsigned int level) const {
- return line_iterator (tria,
- tria->last_line(level)->level(),
- tria->last_line(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_line_iterator
-MGDoFHandler<dim>::last_active_line (const unsigned int level) const {
- return active_line_iterator (tria,
- tria->last_active_line(level)->level(),
- tria->last_active_line(level)->index(),
- this);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_quad_iterator
-MGDoFHandler<dim>::last_raw_quad (const unsigned int level) const {
- return raw_quad_iterator (tria,
- tria->last_raw_quad(level)->level(),
- tria->last_raw_quad(level)->index(),
- this);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::quad_iterator
-MGDoFHandler<dim>::last_quad (const unsigned int level) const {
- return quad_iterator (tria,
- tria->last_quad(level)->level(),
- tria->last_quad(level)->index(),
- this);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_quad_iterator
-MGDoFHandler<dim>::last_active_quad (const unsigned int level) const {
- return active_quad_iterator (tria,
- tria->last_active_quad(level)->level(),
- tria->last_active_quad(level)->index(),
- this);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_hex_iterator
-MGDoFHandler<dim>::last_raw_hex (const unsigned int level) const {
- return raw_hex_iterator (tria,
- tria->last_raw_hex(level)->level(),
- tria->last_raw_hex(level)->index(),
- this);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::hex_iterator
-MGDoFHandler<dim>::last_hex (const unsigned int level) const {
- return hex_iterator (tria,
- tria->last_hex(level)->level(),
- tria->last_hex(level)->index(),
- this);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_hex_iterator
-MGDoFHandler<dim>::last_active_hex (const unsigned int level) const {
- return active_hex_iterator (tria,
- tria->last_active_hex(level)->level(),
- tria->last_active_hex(level)->index(),
- this);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_line_iterator
-MGDoFHandler<dim>::last_raw_line () const {
- return last_raw_line (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_quad_iterator
-MGDoFHandler<dim>::last_raw_quad () const {
- return last_raw_quad (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::raw_hex_iterator
-MGDoFHandler<dim>::last_raw_hex () const {
- return last_raw_hex (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::line_iterator
-MGDoFHandler<dim>::last_line () const {
- return last_line (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::quad_iterator
-MGDoFHandler<dim>::last_quad () const {
- return last_quad (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::hex_iterator
-MGDoFHandler<dim>::last_hex () const {
- return last_hex (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_line_iterator
-MGDoFHandler<dim>::last_active_line () const {
- return last_active_line (mg_levels.size()-1);
-};
-
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_quad_iterator
-MGDoFHandler<dim>::last_active_quad () const {
- return last_active_quad (mg_levels.size()-1);
-};
-
-
-
-template <int dim>
-typename MGDoFHandler<dim>::active_hex_iterator
-MGDoFHandler<dim>::last_active_hex () const {
- return last_active_hex (mg_levels.size()-1);
-};
-
-
-
-
-
-
-//------------------------------------------------------------------
-
-
-
-
-
-
-
-template <int dim>
-void MGDoFHandler<dim>::distribute_dofs (const FiniteElement<dim> &fe,
- unsigned int offset)
-{
- // first distribute global dofs
- DoFHandler<dim>::distribute_dofs (fe, offset);
-
-
- // reserve space for the MG dof numbers
- reserve_space ();
- mg_used_dofs.resize (tria->n_levels(), 0);
-
- // clear user flags because we will
- // need them
- tria->clear_user_flags ();
-
- // now distribute indices on each level
- // separately
- for (unsigned int level=0; level<tria->n_levels(); ++level)
- {
- unsigned int next_free_dof = offset;
- cell_iterator cell = begin(level),
- endc = end(level);
-
- for (; cell != endc; ++cell)
- next_free_dof = distribute_dofs_on_cell (cell, next_free_dof);
-
- mg_used_dofs[level] = next_free_dof;
- };
-};
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-unsigned int
-MGDoFHandler<1>::distribute_dofs_on_cell (cell_iterator &cell,
- unsigned int next_free_dof) {
-
- // distribute dofs of vertices
- if (selected_fe->dofs_per_vertex > 0)
- for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
- {
- cell_iterator neighbor = cell->neighbor(v);
-
- if (neighbor.state() == valid)
- {
- // has neighbor already been processed?
- if (neighbor->user_flag_set() &&
- (neighbor->level() == cell->level()))
- // copy dofs if the neighbor is on
- // the same level (only then are
- // mg dofs the same)
- {
- if (v==0)
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- cell->set_mg_vertex_dof_index (0, d,
- neighbor->mg_vertex_dof_index (1, d));
- else
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- cell->set_mg_vertex_dof_index (1, d,
- neighbor->mg_vertex_dof_index (0, d));
-
- // next neighbor
- continue;
- };
- };
-
- // otherwise: create dofs newly
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- cell->set_mg_vertex_dof_index (v, d, next_free_dof++);
- };
-
- // dofs of line
- if (selected_fe->dofs_per_line > 0)
- for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
- cell->set_mg_dof_index (d, next_free_dof++);
-
- // note that this cell has been processed
- cell->set_user_flag ();
-
- return next_free_dof;
-};
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-unsigned int
-MGDoFHandler<2>::distribute_dofs_on_cell (cell_iterator &cell,
- unsigned int next_free_dof) {
- if (selected_fe->dofs_per_vertex > 0)
- // number dofs on vertices
- for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
- // check whether dofs for this
- // vertex have been distributed
- // (only check the first dof)
- if (cell->mg_vertex_dof_index(vertex, 0) == -1)
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
-
- // for the four sides
- if (selected_fe->dofs_per_line > 0)
- for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
- {
- line_iterator line = cell->line(side);
-
- // distribute dofs if necessary:
- // check whether line dof is already
- // numbered (check only first dof)
- if (line->mg_dof_index(0) == -1)
- // if not: distribute dofs
- for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
- line->set_mg_dof_index (d, next_free_dof++);
- };
-
-
- // dofs of quad
- if (selected_fe->dofs_per_quad > 0)
- for (unsigned int d=0; d<selected_fe->dofs_per_quad; ++d)
- cell->set_mg_dof_index (d, next_free_dof++);
-
-
- // note that this cell has been processed
- cell->set_user_flag ();
-
- return next_free_dof;
-};
-
-#endif
-
-
-
-#if deal_II_dimension == 3
-
-template <>
-unsigned int
-MGDoFHandler<3>::distribute_dofs_on_cell (cell_iterator &cell,
- unsigned int next_free_dof) {
- if (selected_fe->dofs_per_vertex > 0)
- // number dofs on vertices
- for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
- // check whether dofs for this
- // vertex have been distributed
- // (only check the first dof)
- if (cell->mg_vertex_dof_index(vertex, 0) == -1)
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
-
- // for the lines
- if (selected_fe->dofs_per_line > 0)
- for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
- {
- line_iterator line = cell->line(l);
-
- // distribute dofs if necessary:
- // check whether line dof is already
- // numbered (check only first dof)
- if (line->mg_dof_index(0) == -1)
- // if not: distribute dofs
- for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
- line->set_mg_dof_index (d, next_free_dof++);
- };
-
- // for the quads
- if (selected_fe->dofs_per_quad > 0)
- for (unsigned int q=0; q<GeometryInfo<3>::lines_per_cell; ++q)
- {
- quad_iterator quad = cell->quad(q);
-
- // distribute dofs if necessary:
- // check whether line dof is already
- // numbered (check only first dof)
- if (quad->mg_dof_index(0) == -1)
- // if not: distribute dofs
- for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
- quad->set_mg_dof_index (d, next_free_dof++);
- };
-
-
- // dofs of cell
- if (selected_fe->dofs_per_hex > 0)
- for (unsigned int d=0; d<selected_fe->dofs_per_hex; ++d)
- cell->set_mg_dof_index (d, next_free_dof++);
-
-
- // note that this cell has been processed
- cell->set_user_flag ();
-
- return next_free_dof;
-};
-
-#endif
-
-
-
-
-template <int dim>
-unsigned int MGDoFHandler<dim>::n_dofs (const unsigned int level) const {
- Assert (level < mg_used_dofs.size(), ExcInvalidLevel(level));
-
- return mg_used_dofs[level];
-};
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-void MGDoFHandler<1>::renumber_dofs (const unsigned int level,
- const vector<int> &new_numbers) {
- Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
-
- // note that we can not use cell iterators
- // in this function since then we would
- // renumber the dofs on the interface of
- // two cells more than once. Anyway, this
- // ways it's not only more correct but also
- // faster
- for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
- i!=mg_vertex_dofs.end(); ++i)
- // if the present vertex lives on
- // the present level
- if ((i->get_coarsest_level() <= level) &&
- (i->get_finest_level() >= level))
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- i->set_index (level, d, selected_fe->dofs_per_vertex,
- new_numbers[i->get_index (level, d,
- selected_fe->dofs_per_vertex)]);
-
- for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
- i!=mg_levels[level]->line_dofs.end(); ++i)
- {
- Assert (*i != -1, ExcInternalError());
- *i = new_numbers[*i];
- };
-};
-
-#endif
-
-
-
-
-#if deal_II_dimension == 2
-
-template <>
-void MGDoFHandler<2>::renumber_dofs (const unsigned int level,
- const vector<int> &new_numbers) {
- Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
-
- for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
- i!=mg_vertex_dofs.end(); ++i)
- // if the present vertex lives on
- // the present level
- if ((i->get_coarsest_level() <= level) &&
- (i->get_finest_level() >= level))
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- i->set_index (level, d, selected_fe->dofs_per_vertex,
- new_numbers[i->get_index (level, d,
- selected_fe->dofs_per_vertex)]);
-
- for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
- i!=mg_levels[level]->line_dofs.end(); ++i)
- {
- Assert (*i != -1, ExcInternalError());
- *i = new_numbers[*i];
- };
-
- for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
- i!=mg_levels[level]->quad_dofs.end(); ++i)
- {
- Assert (*i != -1, ExcInternalError());
- *i = new_numbers[*i];
- };
-};
-
-#endif
-
-
-
-#if deal_II_dimension == 3
-
-template <>
-void MGDoFHandler<3>::renumber_dofs (const unsigned int level,
- const vector<int> &new_numbers) {
- Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
-
- for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
- i!=mg_vertex_dofs.end(); ++i)
- // if the present vertex lives on
- // the present level
- if ((i->get_coarsest_level() <= level) &&
- (i->get_finest_level() >= level))
- for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
- i->set_index (level, d, selected_fe->dofs_per_vertex,
- new_numbers[i->get_index (level, d,
- selected_fe->dofs_per_vertex)]);
-
- for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
- i!=mg_levels[level]->line_dofs.end(); ++i)
- {
- Assert (*i != -1, ExcInternalError());
- *i = new_numbers[*i];
- };
-
- for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
- i!=mg_levels[level]->quad_dofs.end(); ++i)
- {
- Assert (*i != -1, ExcInternalError());
- *i = new_numbers[*i];
- };
-
- for (vector<int>::iterator i=mg_levels[level]->hex_dofs.begin();
- i!=mg_levels[level]->hex_dofs.end(); ++i)
- {
- Assert (*i != -1, ExcInternalError());
- *i = new_numbers[*i];
- };
-};
-
-#endif
-
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-void MGDoFHandler<1>::reserve_space () {
- const unsigned int dim = 1;
-
- Assert (selected_fe != 0, ExcNoFESelected());
- Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
-
- //////////////////////////
- // DESTRUCTION
-
- // delete all levels and set them up
- // newly, since vectors are
- // troublesome if you want to change
- // their size
- for (unsigned int i=0; i<mg_levels.size(); ++i)
- delete mg_levels[i];
- mg_levels.resize (0);
-
- // also delete vector of vertex indices
- // this calls the destructor which
- // must free the space
- mg_vertex_dofs.resize (0);
-
- ////////////////////////////
- // CONSTRUCTION
-
- // first allocate space for the
- // lines on each level
- for (unsigned int i=0; i<tria->n_levels(); ++i)
- {
- mg_levels.push_back (new DoFLevel<1>);
-
- mg_levels.back()->line_dofs = vector<int>(tria->levels[i]->lines.lines.size() *
- selected_fe->dofs_per_line,
- -1);
- };
-
- // now allocate space for the
- // vertices. To this end, we need
- // to construct as many objects as
- // there are vertices and let them
- // allocate enough space for their
- // vertex indices on the levels they
- // live on. We need therefore to
- // count to how many levels a cell
- // belongs to, which we do by looping
- // over all cells and storing the
- // maximum and minimum level each
- // vertex we pass by belongs to
- mg_vertex_dofs.resize (tria->vertices.size());
-
- vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
- vector<unsigned int> max_level (tria->vertices.size(), 0);
-
- Triangulation<dim>::cell_iterator cell = tria->begin(),
- endc = tria->end();
- for (; cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- {
- const unsigned int vertex_index = cell->vertex_index(vertex);
- if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
- min_level[vertex_index] = cell->level();
- if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
- max_level[vertex_index] = cell->level();
- };
-
- // now allocate the needed space
- for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
- {
- Assert (min_level[vertex] < tria->n_levels(), ExcInternalError());
- Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
-
- mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex],
- selected_fe->dofs_per_vertex);
- };
-};
-
-#endif
-
-
-
-#if deal_II_dimension == 2
-
-template <>
-void MGDoFHandler<2>::reserve_space () {
- const unsigned int dim = 2;
-
- Assert (selected_fe != 0, ExcNoFESelected());
- Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
-
- ////////////////////////////
- // DESTRUCTION
-
- // delete all levels and set them up
- // newly, since vectors are
- // troublesome if you want to change
- // their size
- for (unsigned int i=0; i<mg_levels.size(); ++i)
- delete mg_levels[i];
- mg_levels.clear ();
-
- // also delete vector of vertex indices
- // this calls the destructor which
- // must free the space
- mg_vertex_dofs.clear ();
-
-
- ////////////////////////////
- // CONSTRUCTION
-
- // first allocate space for the
- // lines and quads on each level
- for (unsigned int i=0; i<tria->n_levels(); ++i)
- {
- mg_levels.push_back (new DoFLevel<2>);
-
- mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
- selected_fe->dofs_per_line,
- -1);
- mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
- selected_fe->dofs_per_quad,
- -1);
- };
-
-
- // now allocate space for the
- // vertices. To this end, we need
- // to construct as many objects as
- // there are vertices and let them
- // allocate enough space for their
- // vertex indices on the levels they
- // live on. We need therefore to
- // count to how many levels a cell
- // belongs to, which we do by looping
- // over all cells and storing the
- // maximum and minimum level each
- // vertex we pass by belongs to
- mg_vertex_dofs.resize (tria->vertices.size());
-
- vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
- vector<unsigned int> max_level (tria->vertices.size(), 0);
-
- Triangulation<dim>::cell_iterator cell = tria->begin(),
- endc = tria->end();
- for (; cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- {
- const unsigned int vertex_index = cell->vertex_index(vertex);
- if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
- min_level[vertex_index] = cell->level();
- if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
- max_level[vertex_index] = cell->level();
- };
-
-
- // now allocate the needed space
- for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
- {
- Assert (min_level[vertex] < tria->n_levels(), ExcInternalError());
- Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
-
- mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex],
- selected_fe->dofs_per_vertex);
- };
-};
-
-#endif
-
-
-
-
-#if deal_II_dimension == 3
-
-template <>
-void MGDoFHandler<3>::reserve_space () {
- const unsigned int dim = 3;
-
- Assert (selected_fe != 0, ExcNoFESelected());
- Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
-
- ////////////////////////////
- // DESTRUCTION
-
- // delete all levels and set them up
- // newly, since vectors are
- // troublesome if you want to change
- // their size
- for (unsigned int i=0; i<mg_levels.size(); ++i)
- delete mg_levels[i];
- mg_levels.clear ();
-
- // also delete vector of vertex indices
- // this calls the destructor which
- // must free the space
- mg_vertex_dofs.clear ();
-
-
- ////////////////////////////
- // CONSTRUCTION
-
- // first allocate space for the
- // lines and quads on each level
- for (unsigned int i=0; i<tria->n_levels(); ++i)
- {
- mg_levels.push_back (new DoFLevel<3>);
-
- mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
- selected_fe->dofs_per_line,
- -1);
- mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
- selected_fe->dofs_per_quad,
- -1);
- mg_levels.back()->hex_dofs = vector<int> (tria->levels[i]->hexes.hexes.size() *
- selected_fe->dofs_per_hex,
- -1);
- };
-
-
- // now allocate space for the
- // vertices. To this end, we need
- // to construct as many objects as
- // there are vertices and let them
- // allocate enough space for their
- // vertex indices on the levels they
- // live on. We need therefore to
- // count to how many levels a cell
- // belongs to, which we do by looping
- // over all cells and storing the
- // maximum and minimum level each
- // vertex we pass by belongs to
- mg_vertex_dofs.resize (tria->vertices.size());
-
- vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
- vector<unsigned int> max_level (tria->vertices.size(), 0);
-
- Triangulation<dim>::cell_iterator cell = tria->begin(),
- endc = tria->end();
- for (; cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- {
- const unsigned int vertex_index = cell->vertex_index(vertex);
- if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
- min_level[vertex_index] = cell->level();
- if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
- max_level[vertex_index] = cell->level();
- };
-
-
- // now allocate the needed space
- for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
- {
- Assert (min_level[vertex] < tria->n_levels(), ExcInternalError());
- Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
-
- mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex],
- selected_fe->dofs_per_vertex);
- };
-};
-
-#endif
-
-
-// explicite instantiations
-template class MGDoFHandler<deal_II_dimension>;
+++ /dev/null
-/* $Id$ */
-
-#include <lac/sparsity_pattern.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <dofs/mg_dof_tools.h>
-#include <fe/fe.h>
-
-
-template <int dim>
-void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<dim> &mg_dof_handler,
- const unsigned int level,
- SparsityPattern &sparsity)
-{
- Assert (sparsity.n_rows() == mg_dof_handler.n_dofs(level),
- ExcDimensionMismatch (sparsity.n_rows(), mg_dof_handler.n_dofs(level)));
- Assert (sparsity.n_cols() == mg_dof_handler.n_dofs(level),
- ExcDimensionMismatch (sparsity.n_cols(), mg_dof_handler.n_dofs(level)));
-
- const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
- vector<int> dofs_on_this_cell(dofs_per_cell);
- MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(level),
- endc = mg_dof_handler.end(level);
- for (; cell!=endc; ++cell)
- {
- cell->get_mg_dof_indices (dofs_on_this_cell);
- // make sparsity pattern for this cell
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- };
-};
-
-
-
-// explicit instantiations
-template void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
- const unsigned int,
- SparsityPattern &);
-
#include <grid/tria.h>
#include <dofs/dof_handler.h>
-#include <dofs/mg_dof_handler.h>
#include <grid/tria_accessor.h>
#include <dofs/dof_accessor.h>
-#include <dofs/mg_dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/intergrid_map.h>
// $Id$
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_smoother.h>
#include <iostream>
#include <cmath>
-MGBase::~MGBase ()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-MGSmootherBase::~MGSmootherBase()
-{};
-
-
//////////////////////////////////////////////////////////////////////
-void
-MGSmootherIdentity::smooth (const unsigned int,
- Vector<double> &,
- const Vector<double> &) const
+MGBase::~MGBase ()
{};
-//////////////////////////////////////////////////////////////////////
-
-
-
MGBase::MGBase(const MGTransferBase &transfer,
const unsigned minlevel,
const unsigned maxlevel)
-#include <dofs/mg_dof_accessor.h>
-#include <dofs/mg_dof_handler.h>
#include <dofs/dof_levels.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
#include <grid/tria_iterator.h>
#include <grid/tria_iterator.templates.h>
#include <fe/fe.h>
#include <lac/vector.h>
-
-
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+
#if deal_II_dimension == 3
template <>
/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-#include <dofs/mg_dof_handler.h>
#include <dofs/dof_levels.h>
-#include <dofs/mg_dof_accessor.h>
#include <dofs/dof_constraints.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_levels.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria.h>
#include <grid/geometry_info.h>
#include <fe/fe.h>
+#include <lac/sparse_matrix.h>
#include <algorithm>
/* $Id$ */
+
#include <lac/sparsity_pattern.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/mg_dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_iterator.h>
-#include <dofs/mg_dof_tools.h>
+#include <multigrid/mg_dof_tools.h>
#include <fe/fe.h>
/* $Id$ */
+
#include <grid/tria.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/mg_dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_iterator.h>
#include <fe/fe.h>
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
+#include <multigrid/multigrid.h>
+#include <multigrid/mg_smoother.h>
+#include <multigrid/mg_smoother.templates.h>
#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
#include <algorithm>
+
+//////////////////////////////////////////////////////////////////////
+
+
+MGSmootherBase::~MGSmootherBase()
+{};
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+ Vector<double> &,
+ const Vector<double> &) const
+{};
+
+
+
+
#if deal_II_dimension == 1
MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps)
//////////////////////////////////////////////////////////////////////
-#include <numerics/mg_smoother.templates.h>
// explicit instantiations
// don't do the following instantiation in 1d, since there is a specialized
template
MGSmootherRelaxation<float>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- function_ptr relaxation,
- unsigned int steps,
- double omega);
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const MGLevelObject<SparseMatrix<float> > &matrix,
+ const function_ptr relaxation,
+ const unsigned int steps,
+ const double omega);
+
template
MGSmootherRelaxation<double>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGMatrix<SparseMatrix<double> >& matrix,
- function_ptr relaxation,
- unsigned int steps,
- double omega);
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const MGLevelObject<SparseMatrix<double> > &matrix,
+ const function_ptr relaxation,
+ const unsigned int steps,
+ const double omega);
// $Id$
-// Copyright Guido Kanschat, Universitaet Heidelberg, 1999
+
#include <grid/tria.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/mg_dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_iterator.h>
#include <fe/fe.h>
-#include <numerics/multigrid.h>
-#include <numerics/multigrid.templates.h>
-#include <numerics/mg_smoother.h>
+#include <multigrid/multigrid.h>
+#include <multigrid/multigrid.templates.h>
+#include <multigrid/mg_smoother.h>
#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
+
/* --------------------------------- MG ------------------------------------ */
template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim> &mg_dof_handler,
- ConstraintMatrix &constraints,
- const MGMatrix<SparseMatrix<double> > &level_matrices,
- const MGTransferBase &transfer,
- unsigned int minl, unsigned int maxl)
+Multigrid<dim>::Multigrid (const MGDoFHandler<dim> &mg_dof_handler,
+ const ConstraintMatrix &constraints,
+ const MGLevelObject<SparsityPattern> &level_sparsities,
+ const MGLevelObject<SparseMatrix<double> > &level_matrices,
+ const MGTransferBase &transfer,
+ const unsigned int minlevel,
+ const unsigned int maxlevel)
:
- MGBase(transfer, minl,
- min(mg_dof_handler.get_tria().n_levels()-1, maxl)),
+ MGBase(transfer, minlevel,
+ min(mg_dof_handler.get_tria().n_levels()-1,
+ maxlevel)),
mg_dof_handler(&mg_dof_handler),
+ level_sparsities(&level_sparsities),
level_matrices(&level_matrices),
constraints(&constraints)
-{
- // initialize the objects with
- // their correct size
- for (unsigned int level=minlevel; level<=maxlevel ; ++level)
- {
- solution[level].reinit(level_matrices[level].m());
- defect[level].reinit(level_matrices[level].m());
- };
-};
+{};
template <int dim>
void
-MG<dim>::level_vmult (const unsigned int level,
- Vector<double> &result,
- const Vector<double> &u,
- const Vector<double> &/* rhs */)
+Multigrid<dim>::level_vmult (const unsigned int level,
+ Vector<double> &result,
+ const Vector<double> &u,
+ const Vector<double> &/* rhs */)
{
(*level_matrices)[level].vmult(result,u);
result.scale(-1.);
// level which have children
for (unsigned int level=0; level<n_levels-1; ++level)
{
- prolongation_sparsities.push_back (SparsityPattern ());
+ prolongation_sparsities.push_back (SparsityPattern());
// reset the dimension of the structure.
// note that for the number of entries
// per row, the number of mother dofs
// explicit instatations
-template class MG<deal_II_dimension>;
+template class Multigrid<deal_II_dimension>;
template
void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
+template
+void Multigrid<deal_II_dimension>::copy_to_mg (const Vector<float>& src);
+
+template
+void Multigrid<deal_II_dimension>::copy_to_mg (const Vector<double>& src);
+
+template
+void Multigrid<deal_II_dimension>::copy_from_mg (Vector<float>& src) const;
+
+template
+void Multigrid<deal_II_dimension>::copy_from_mg (Vector<double>& src) const;
+
+
+++ /dev/null
-/* $Id$ */
-
-#include <grid/tria.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <fe/fe.h>
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-
-#include <algorithm>
-
-
-#if deal_II_dimension == 1
-
-MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps)
- :
- steps(steps)
-{
- Assert (false, ExcNotImplemented());
-};
-
-#endif
-
-
-#if deal_II_dimension > 1
-
-template <int dim>
-MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
- :
- steps(steps)
-{
- const unsigned int n_levels = mg_dof.get_tria().n_levels();
-
- // allocate the right number of
- // elements
- interior_boundary_dofs.resize (n_levels-1);
-
- // use a temporary to store the
- // indices. this allows to not allocate
- // to much memory by later copying the
- // content of this vector to its final
- // destination
- vector<int> boundary_dofs;
-
- // temporary to hold the dof indices
- // on a face between to levels
- vector<int> dofs_on_face (mg_dof.get_fe().dofs_per_face);
-
- for (unsigned int level=1; level<n_levels; ++level)
- {
- boundary_dofs.clear ();
-
- // for each cell on this level:
- // find out whether a face is
- // at the boundary of this level's
- // cells and if so add the dofs
- // to the interior boundary dofs
- for (MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
- cell != mg_dof.end(level); ++cell)
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if ((cell->neighbor(face).state() == valid) &&
- (static_cast<unsigned int>(cell->neighbor(face)->level())
- == level-1))
- {
- // get indices of this face
- cell->face(face)->get_mg_dof_indices (dofs_on_face);
- // append them to the levelwise
- // list
- boundary_dofs.insert (boundary_dofs.end(),
- dofs_on_face.begin(),
- dofs_on_face.end());
- };
-
- // now sort the list of interior boundary
- // dofs and eliminate duplicates
- sort (boundary_dofs.begin(), boundary_dofs.end());
- boundary_dofs.erase (unique (boundary_dofs.begin(),
- boundary_dofs.end()),
- boundary_dofs.end());
-
- // now finally copy the result
- // for this level its destination
- interior_boundary_dofs[level-1] = boundary_dofs;
- };
-};
-
-#endif
-
-//////////////////////////////////////////////////////////////////////
-
-void
-MGSmoother::set_zero_interior_boundary (const unsigned int level,
- Vector<double> &u) const
-{
- if (level==0)
- return;
- else
- for (vector<int>::const_iterator p=interior_boundary_dofs[level-1].begin();
- p!=interior_boundary_dofs[level-1].end(); ++p)
- u(*p) = 0;
-};
-
-//////////////////////////////////////////////////////////////////////
-
-#include <numerics/mg_smoother.templates.h>
-
-// explicit instantiations
-// don't do the following instantiation in 1d, since there is a specialized
-// function there
-#if deal_II_dimension > 1
-template MGSmoother::MGSmoother (const MGDoFHandler<deal_II_dimension>&, unsigned int);
-#endif
-
-template
-MGSmootherRelaxation<float>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- function_ptr relaxation,
- unsigned int steps,
- double omega);
-template
-MGSmootherRelaxation<double>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGMatrix<SparseMatrix<double> >& matrix,
- function_ptr relaxation,
- unsigned int steps,
- double omega);
-
+++ /dev/null
-// $Id$
-// Copyright Guido Kanschat, Universitaet Heidelberg, 1999
-
-#include <grid/tria.h>
-#include <dofs/mg_dof_handler.h>
-#include <dofs/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <fe/fe.h>
-#include <numerics/multigrid.h>
-#include <numerics/multigrid.templates.h>
-#include <numerics/mg_smoother.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-
-
-/* --------------------------------- MG ------------------------------------ */
-
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim> &mg_dof_handler,
- ConstraintMatrix &constraints,
- const MGMatrix<SparseMatrix<double> > &level_matrices,
- const MGTransferBase &transfer,
- unsigned int minl, unsigned int maxl)
- :
- MGBase(transfer, minl,
- min(mg_dof_handler.get_tria().n_levels()-1, maxl)),
- mg_dof_handler(&mg_dof_handler),
- level_matrices(&level_matrices),
- constraints(&constraints)
-{
- // initialize the objects with
- // their correct size
- for (unsigned int level=minlevel; level<=maxlevel ; ++level)
- {
- solution[level].reinit(level_matrices[level].m());
- defect[level].reinit(level_matrices[level].m());
- };
-};
-
-
-
-template <int dim>
-void
-MG<dim>::level_vmult (const unsigned int level,
- Vector<double> &result,
- const Vector<double> &u,
- const Vector<double> &/* rhs */)
-{
- (*level_matrices)[level].vmult(result,u);
- result.scale(-1.);
-}
-
-
-
-
-
-
-/* ----------------------------- MGTransferPrebuilt ------------------------ */
-
-
-
-MGTransferPrebuilt::~MGTransferPrebuilt ()
-{};
-
-
-
-template <int dim>
-void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
-{
- const unsigned int n_levels = mg_dof.get_tria().n_levels();
- const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
-
- // reset the size of the array of
- // matrices
- prolongation_sparsities.clear ();
- prolongation_matrices.clear ();
- prolongation_sparsities.reserve (n_levels-1);
- prolongation_matrices.reserve (n_levels-1);
-
-
- // two fields which will store the
- // indices of the multigrid dofs
- // for a cell and one of its children
- vector<int> dof_indices_mother (dofs_per_cell);
- vector<int> dof_indices_child (dofs_per_cell);
-
- // for each level: first build the sparsity
- // pattern of the matrices and then build the
- // matrices themselves. note that we only
- // need to take care of cells on the coarser
- // level which have children
- for (unsigned int level=0; level<n_levels-1; ++level)
- {
- prolongation_sparsities.push_back (SparsityPattern ());
- // reset the dimension of the structure.
- // note that for the number of entries
- // per row, the number of mother dofs
- // coupling to a child dof is
- // necessary. this, of course, is the
- // number of degrees of freedom per
- // cell
- prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
- mg_dof.n_dofs(level),
-//TODO: evil hack, must be corrected!
- dofs_per_cell+1);
-
- for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
- cell != mg_dof.end(level); ++cell)
- if (cell->has_children())
- {
- cell->get_mg_dof_indices (dof_indices_mother);
-
- for (unsigned int child=0;
- child<GeometryInfo<dim>::children_per_cell; ++child)
- {
- // set an alias to the
- // prolongation matrix for
- // this child
- const FullMatrix<double> &prolongation
- = mg_dof.get_fe().prolongate(child);
-
- cell->child(child)->get_mg_dof_indices (dof_indices_child);
-
- // now tag the entries in the
- // matrix which will be used
- // for this pair of mother/child
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (prolongation(i,j) != 0)
- {
- prolongation_sparsities[level].add
- (dof_indices_child[i],
- dof_indices_mother[j]);
- };
- };
- };
- prolongation_sparsities[level].compress ();
-
- prolongation_matrices.push_back (SparseMatrix<float>());
- prolongation_matrices[level].reinit (prolongation_sparsities[level]);
-
-
- // now actually build the matrices
- for (MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
- cell != mg_dof.end(level); ++cell)
- if (cell->has_children())
- {
- cell->get_mg_dof_indices (dof_indices_mother);
-
- for (unsigned int child=0;
- child<GeometryInfo<dim>::children_per_cell; ++child)
- {
- // set an alias to the
- // prolongation matrix for
- // this child
- const FullMatrix<double> &prolongation
- = mg_dof.get_fe().prolongate(child);
-
- cell->child(child)->get_mg_dof_indices (dof_indices_child);
-
- // now set the entries in the
- // matrix
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (prolongation(i,j) != 0)
- {
- prolongation_matrices[level].set
- (dof_indices_child[i],
- dof_indices_mother[j],
- prolongation(i,j));
- };
- };
- };
- };
-};
-
-
-
-void MGTransferPrebuilt::prolongate (const unsigned int to_level,
- Vector<double> &dst,
- const Vector<double> &src) const
-{
- Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
- ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
-
- prolongation_matrices[to_level-1].vmult (dst, src);
-};
-
-
-
-void MGTransferPrebuilt::restrict_and_add (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const
-{
- Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
- ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
-
- prolongation_matrices[from_level-1].Tvmult_add (dst, src);
-};
-
-
-
-// explicit instatations
-template class MG<deal_II_dimension>;
-template
-void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
-
+++ /dev/null
-/*---------------------------- mgbase.h ---------------------------*/
-/* $Id$ */
-#ifndef __mgbase_H
-#define __mgbase_H
-/*---------------------------- mgbase.h ---------------------------*/
-
-
-
-#include <base/subscriptor.h>
-#include <base/smartpointer.h>
-#include <lac/forward_declarations.h>
-#include <lac/vector.h>
-
-#include <vector>
-
-
-
-/**
- * Abstract base class for coarse grid solvers. The interface of a
- * function call operator is defined to execute coarse grid solution
- * in a derived class.
- *
- * @author Guido Kanschat, 1999
- */
-class MGCoarseGridSolver : public Subscriptor
-{
- public:
- /**
- * Virtual destructor. Does
- * nothing in particular, but
- * needs to be declared anyhow.
- */
- virtual ~MGCoarseGridSolver();
-
- /**
- * Coarse grid solving method.
- * This is only the interface for
- * a function defined in a derived
- * class like
- * #MGCoarseGridLACIteration#.
- *
- * Note that the information
- * about the matrix is removed to
- * that class.
- */
- virtual void operator() (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src) const = 0;
-};
-
-
-
-
-/**
- * Abstract base class for multigrid smoothers.
- * In fact, this class only provides the interface of the smoothing function.
- * Using #deal.II# grid handling, #MGSmoother# is a good point to start.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-class MGSmootherBase : public Subscriptor
-{
- public:
- /**
- * Virtual destructor. Does
- * nothing in particular, but
- * needs to be declared anyhow.
- */
- virtual ~MGSmootherBase();
-
- /**
- * Smooth the residual of #u# on
- * the given level. If $S$ is the
- * smoothing operator, then this
- * function should do the
- * following operation:
- * $u += S (rhs - Au)$, where #u# and
- * #rhs# are the input parameters.
- *
- * This function should keep the
- * interior level boundary
- * values, so you may want to
- * call #set_zero_interior_boundary#
- * in the derived class
- * #MGSmoother# somewhere in your
- * derived function, or another
- * function doing similar things.
- */
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const = 0;
-};
-
-
-
-
-
-/**
- * Smoother doing nothing. This class is not useful for many applications other
- * than for testing some multigrid procedures. Also some applications might
- * get convergence without smoothing and then this class brings you the
- * cheapest possible multigrid.
- *
- * @author Guido Kanschat, 1999
- */
-class MGSmootherIdentity : public MGSmootherBase
-{
- public:
- /**
- * Implementation of the
- * interface in #MGSmootherBase#.
- * This function does nothing,
- * which by comparison with the
- * definition of this function
- * means that the the smoothing
- * operator equals the null
- * operator.
- */
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const;
-};
-
-
-
-
-/**
- * Base class used to declare the operations needed by a concrete class
- * implementing prolongation and restriction of vectors in the multigrid
- * context. This class is an abstract one and has no implementations of
- * possible algorithms for these operations.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-class MGTransferBase : public Subscriptor
-{
- public:
- /**
- * Destructor. Does nothing here, but
- * needs to be declared virtual anyway.
- */
- virtual ~MGTransferBase();
-
- /**
- * Prolongate a vector from level
- * #to_level-1# to level
- * #to_level#. The previous
- * content of #dst# is
- * overwritten.
- *
- * #src# is assumed to be a vector with
- * as many elements as there are degrees
- * of freedom on the coarser level of
- * the two involved levels, while #src#
- * shall have as many elements as there
- * are degrees of freedom on the finer
- * level.
- */
- virtual void prolongate (const unsigned int to_level,
- Vector<double> &dst,
- const Vector<double> &src) const = 0;
-
- /**
- * Restrict a vector from level
- * #from_level# to level
- * #from_level-1# and add this
- * restriction to
- * #dst#. Obviously, if the
- * refined region on level
- * #from_level# is smaller than
- * that on level #from_level-1#,
- * some degrees of freedom in
- * #dst# are not covered and will
- * not be altered. For the other
- * degress of freedom, the result
- * of the restriction is added.
- *
- * #src# is assumed to be a vector with
- * as many elements as there are degrees
- * of freedom on the finer level of
- * the two involved levels, while #src#
- * shall have as many elements as there
- * are degrees of freedom on the coarser
- * level.
- */
- virtual void restrict_and_add (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const = 0;
-};
-
-
-
-/**
- * An array with a vector for each level. The purpose of this class
- * is mostly to allow access by level number, even if the lower levels
- * are not used and therefore have no vector at all; this is done by
- * simply shifting the given index by the minimum level we have
- * stored.
- *
- * @author Guido Kanschat, 1999
- */
-template<class VECTOR = Vector<double> >
-class MGVector : public Subscriptor
-{
- public:
- /**
- * Constructor allowing to
- * initialize the number of
- * levels.
- */
- MGVector (const unsigned int minlevel,
- const unsigned int maxlevel);
-
- /**
- * Access vector on level #level#.
- */
- VECTOR & operator[] (const unsigned int level);
-
- /**
- * Access vector on level
- * #level#. Constant version.
- */
- const VECTOR & operator[] (const unsigned int level) const;
-
- /**
- * Call #clear# on all vectors
- * stored by this object. Note
- * that if #VECTOR==Vector<T>#,
- * #clear# will set all entries
- * to zero, while if
- * #VECTOR==std::vector<T>#,
- * #clear# deletes the elements
- * of the vectors. This class
- * might therefore not be useful
- * for STL vectors.
- */
- void clear();
-
- private:
- /**
- * Level of first component.
- */
- const unsigned int minlevel;
-
- /**
- * Array of the vectors to be held.
- */
- vector<VECTOR> vectors;
-};
-
-
-
-
-/**
- * An array of matrices for each level.
- * This class provides the functionality of #vector<MATRIX># combined
- * with a subscriptor for smart pointers.
- * @author Guido Kanschat, 1999
- */
-template <class MATRIX = SparseMatrix<double> >
-class MGMatrix : public Subscriptor,
- public vector<MATRIX>
-{
- public:
- /**
- * Constructor allowing to initialize the number of levels.
- */
- MGMatrix(unsigned int minlevel, unsigned int maxlevel);
-
- /**
- * Safe access operator.
- */
- MATRIX& operator[](unsigned int);
-
- /**
- * Safe access operator.
- */
- const MATRIX& operator[](unsigned int) const;
-
- private:
- /**
- * Level of first component.
- */
- unsigned int minlevel;
-};
-
-
-
-
-/**
- * Coarse grid solver using LAC iterative methods.
- * This is a little wrapper, transforming a triplet of iterative
- * solver, matrix and preconditioner into a coarse grid solver.
- *
- * The type of the matrix (i.e. the template parameter #MATRIX#)
- * should be derived from #Subscriptor# to allow for the use of a
- * smart pointer to it.
- *
- * @author Guido Kanschat, 1999
- */
-template<class SOLVER, class MATRIX, class PRECOND>
-class MGCoarseGridLACIteration : public MGCoarseGridSolver
-{
- public:
- /**
- * Constructor.
- * Store solver, matrix and
- * preconditioning method for later
- * use.
- */
- MGCoarseGridLACIteration (SOLVER &,
- const MATRIX &,
- const PRECOND &);
-
- /**
- * Implementation of the abstract
- * function.
- * Calls the solver method with
- * matrix, vectors and
- * preconditioner.
- */
- virtual void operator() (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src) const;
- private:
- /**
- * Reference to the solver.
- */
- SOLVER& solver;
-
- /**
- * Reference to the matrix.
- */
- const SmartPointer<const MATRIX> matrix;
-
- /**
- * Reference to the preconditioner.
- */
- const PRECOND& precondition;
-};
-
-
-
-
-/**
- * Basic class for preconditioning by multigrid.
- *
- * The functionality of the multigrid method is restricted to defect
- * correction. It is {\bf not} iterative and the start solution is
- * always zero. Since by this $u^E_l$ and $u^A_l$ (see report on
- * multigrid) are always zero, restriction is simplified a lot and
- * maybe even the seam condition on grids is oblivious. Still, I am
- * not sure that these restrictions on the class might cause numerical
- * inefficiencies.
- *
- * The function #precondition# is the actual multigrid method and
- * makes use of several operations to be implemented in derived
- * classes. It takes a defect in #src# and the result is the multigrid
- * preconditioned defect in #dst#.
- *
- * @author Guido Kanschat, 1999
- */
-class MGBase : public Subscriptor
-{
- private:
- /**
- * Copy constructor. Made private
- * to prevent use.
- */
- MGBase(const MGBase&);
-
- /**
- * Copy operator. Made private
- * to prevent use.
- */
- const MGBase& operator=(const MGBase&);
-
- public:
- /**
- * Constructor, subject to change.
- */
- MGBase (const MGTransferBase &transfer,
- const unsigned int minlevel,
- const unsigned int maxlevel);
-
- /**
- * Virtual destructor.
- */
- virtual ~MGBase();
-
- /**
- * Execute one step of the
- * v-cycle algorithm. This
- * function assumes, that the
- * vector #d# is properly filled
- * with the residual in the outer
- * defect correction
- * scheme. After execution of
- * #vcycle()#, the result is in
- * the vector #s#. We propose to
- * write functions #copy_from_mg#
- * and #copy_to_mg# in derived
- * classes of #MGBase# to access
- * #d# and #s#.
- *
- * The actual work for this
- * function is done in
- * #level_mgstep#.
- */
- void vcycle(const MGSmootherBase &pre_smooth,
- const MGSmootherBase &post_smooth,
- const MGCoarseGridSolver &cgs);
-
- /**
- * Exception.
- */
- DeclException2(ExcSwitchedLevels, int, int,
- << "minlevel and maxlevel switched, should be: "
- << arg1 << "<=" << arg2);
-
- protected:
- /**
- * Highest level of cells.
- */
- unsigned int maxlevel;
-
- /**
- * Level for coarse grid solution.
- */
- unsigned int minlevel;
-
- /**
- * Auxiliary vector.
- */
- MGVector<Vector<double> > defect;
-
- /**
- * Auxiliary vector.
- */
- MGVector<Vector<double> > solution;
-
- /**
- * Prolongation and restriction object.
- */
- SmartPointer<const MGTransferBase> transfer;
-
- /**
- * The actual v-cycle multigrid method.
- * #level# is the level the
- * function should work on. It
- * will usually be called for the
- * highest level from outside,
- * but will then call itself
- * recursively for #level-1#,
- * unless we are on #minlevel#
- * where instead of recursing
- * deeper, the coarse grid solver
- * is used to solve the matrix of
- * this level.
- */
- void level_mgstep (const unsigned int level,
- const MGSmootherBase &pre_smooth,
- const MGSmootherBase &post_smooth,
- const MGCoarseGridSolver &cgs);
-
- /**
- * Apply negative #vmult# on all
- * cells of a level.
- * This is implemented in a
- * derived class.
- */
- virtual void level_vmult (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src,
- const Vector<double> &rhs) = 0;
-
- private:
- /**
- * Auxiliary vector.
- */
- Vector<double> t;
-};
-
-
-
-/**
- * Multi-level preconditioner.
- * Here, we collect all information needed for multi-level preconditioning
- * and provide the standard interface for LAC iterative methods.
- *
- * The template parameter class #MG# is required to inherit #MGBase#.
- * Furthermore, it needs functions #void copy_to_mg(const VECTOR&)#
- * to store #src# in the right hand side of the multi-level method and
- * #void copy_from_mg(VECTOR&)# to store the result of the v-cycle in #dst#.
- *
- * @author Guido Kanschat, 1999
- */
-template<class MG, class VECTOR = Vector<double> >
-class PreconditionMG
-{
- public:
- /**
- * Constructor.
- * Arguments are the multigrid object,
- * pre-smoother, post-smoother and
- * coarse grid solver.
- */
- PreconditionMG(MG &mg,
- const MGSmootherBase &pre,
- const MGSmootherBase &post,
- const MGCoarseGridSolver &coarse);
-
- /**
- * Preconditioning operator,
- * calling the #vcycle# function
- * of the #MG# object passed to
- * the constructor.
- *
- * This is the operator used by
- * LAC iterative solvers.
- */
- void operator() (VECTOR &dst,
- const VECTOR &src) const;
-
- private:
- /**
- * The multigrid object.
- */
- SmartPointer<MG> multigrid;
-
- /**
- * The pre-smoothing object.
- */
- SmartPointer<const MGSmootherBase> pre;
-
- /**
- * The post-smoothing object.
- */
- SmartPointer<const MGSmootherBase> post;
-
- /**
- * The coarse grid solver.
- */
- SmartPointer<const MGCoarseGridSolver> coarse;
-};
-
-
-
-
-/* ------------------------------------------------------------------- */
-
-template<class VECTOR>
-MGVector<VECTOR>::MGVector(const unsigned int min,
- const unsigned int max)
- :
- minlevel(min),
- vectors(max-min+1)
-{};
-
-
-
-template<class VECTOR>
-VECTOR &
-MGVector<VECTOR>::operator[] (const unsigned int i)
-{
- Assert((i>=minlevel) && (i<minlevel+vectors.size()),
- ExcIndexRange (i, minlevel, minlevel+vectors.size()));
- return vectors[i-minlevel];
-}
-
-
-template<class VECTOR>
-const VECTOR &
-MGVector<VECTOR>::operator[] (const unsigned int i) const
-{
- Assert((i>=minlevel) && (i<minlevel+vectors.size()),
- ExcIndexRange (i, minlevel, minlevel+vectors.size()));
- return vectors[i-minlevel];
-}
-
-
-template<class VECTOR>
-void
-MGVector<VECTOR>::clear()
-{
- typename vector<VECTOR>::iterator v;
- for (v = vectors.begin(); v != vectors.end(); ++v)
- v->clear();
-
-}
-
-
-/* ------------------------------------------------------------------- */
-
-
-template<class MATRIX>
-MGMatrix<MATRIX>::MGMatrix(unsigned int min, unsigned int max)
- :
- vector<MATRIX>(max-min+1),
- minlevel(min)
-{}
-
-
-template<class MATRIX>
-MATRIX&
-MGMatrix<MATRIX>::operator[](unsigned int i)
-{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<MATRIX>::operator[](i-minlevel);
-};
-
-
-
-template<class MATRIX>
-const MATRIX&
-MGMatrix<MATRIX>::operator[](unsigned int i) const
-{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<MATRIX>::operator[](i-minlevel);
-};
-
-
-
-/* ------------------ Functions for MGCoarseGridLACIteration ------------ */
-
-
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::
-MGCoarseGridLACIteration(SOLVER &s,
- const MATRIX &m,
- const PRECOND &p)
- :
- solver(s),
- matrix(&m),
- precondition(p)
-{};
-
-
-
-template<class SOLVER, class MATRIX, class PRECOND>
-void
-MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::
-operator() (const unsigned int /* level */,
- Vector<double> &dst,
- const Vector<double> &src) const
-{
- solver.solve(*matrix, dst, src, precondition);
-}
-
-
-
-
-
-
-/* ------------------------------------------------------------------- */
-
-
-template<class MG, class VECTOR>
-PreconditionMG<MG, VECTOR>::PreconditionMG(MG &mg,
- const MGSmootherBase &pre,
- const MGSmootherBase &post,
- const MGCoarseGridSolver &coarse)
- :
- multigrid(&mg),
- pre(&pre),
- post(&post),
- coarse(&coarse)
-{}
-
-
-
-template<class MG, class VECTOR>
-void
-PreconditionMG<MG,VECTOR>::operator() (VECTOR &dst,
- const VECTOR &src) const
-{
- multigrid->copy_to_mg(src);
- multigrid->vcycle(*pre, *post, *coarse);
- multigrid->copy_from_mg(dst);
-}
-
-
-
-/*---------------------------- mgbase.h ---------------------------*/
-/* end of #ifndef __mgbase_H */
-#endif
-/*---------------------------- mgbase.h ---------------------------*/
+++ /dev/null
-// $Id$
-
-#include <lac/mgbase.h>
-#include <iostream>
-#include <cmath>
-
-
-
-//TODO: this function is only for debugging purposes and should be removed sometimes
-template<class VECTOR>
-static
-void print_vector(ostream& s, const VECTOR& v)
-{
- const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
- unsigned int k=0;
-
- // write the vector entries in a
- // kind of square
- for (unsigned int i=0;i<n;++i)
- {
- for (unsigned int j=0;j<n;++j)
- s << '\t' << v(k++);
- s << endl;
- }
- s << endl;
-}
-
-
-
-MGBase::~MGBase ()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-MGSmootherBase::~MGSmootherBase()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
- Vector<double> &,
- const Vector<double> &) const
-{};
-
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-
-MGBase::MGBase(const MGTransferBase &transfer,
- const unsigned minlevel,
- const unsigned maxlevel)
- :
- maxlevel(maxlevel),
- minlevel(minlevel),
- defect(minlevel,maxlevel),
- solution(minlevel,maxlevel),
- transfer(&transfer)
-{
- Assert(minlevel <= maxlevel,
- ExcSwitchedLevels(minlevel, maxlevel));
-};
-
-
-
-void
-MGBase::vcycle(const MGSmootherBase &pre_smooth,
- const MGSmootherBase &post_smooth,
- const MGCoarseGridSolver &coarse_grid_solver)
-{
- level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
-};
-
-
-
-void
-MGBase::level_mgstep(const unsigned int level,
- const MGSmootherBase &pre_smooth,
- const MGSmootherBase &post_smooth,
- const MGCoarseGridSolver &coarse_grid_solver)
-{
- solution[level] = 0.;
-
- if (level == minlevel)
- {
- coarse_grid_solver(level, solution[level], defect[level]);
- return;
- }
-
- // smoothing of the residual by modifying s
- pre_smooth.smooth(level, solution[level], defect[level]);
- // t = d-As
- t.reinit(solution[level].size());
- level_vmult(level, t, solution[level], defect[level]);
-
- // make t rhs of lower level
-//TODO: this function adds the restricted t to defect[level-1].
-//TODO: why don't we have to clear it before?
- transfer->restrict_and_add (level, defect[level-1], t);
-
- // do recursion
- level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
-
- // reset size of the auxiliary
- // vector, since it has been
- // resized in the recursive call to
- // level_mgstep directly above
- t.reinit(solution[level].size());
-
- // do coarse grid correction
- transfer->prolongate(level, t, solution[level-1]);
- solution[level] += t;
-
- // smoothing (modify solution again)
-//TODO: what happens here? smooth overwrites the solution[level],
-//TODO: so the previous two statements should have no effect. No?
- post_smooth.smooth(level, solution[level], defect[level]);
-};
-
-
-//////////////////////////////////////////////////////////////////////
-
-MGCoarseGridSolver::~MGCoarseGridSolver()
-{};
-
-
-
-//////////////////////////////////////////////////////////////////////
-
-MGTransferBase::~MGTransferBase()
-{};
-
-
-
-
-
/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-#include <grid/tria.h>
-#include <dofs/mg_dof_handler.h>
-#include <grid/tria_accessor.h>
-#include <dofs/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-#include <grid/grid_generator.h>
-#include <base/function.h>
-#include <numerics/data_out.h>
-#include <fe/fe_lib.lagrange.h>
-#include <fe/fe_lib.criss_cross.h>
-#include <fe/fe_update_flags.h>
-#include <base/quadrature_lib.h>
-#include <numerics/assembler.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
#include <lac/vector.h>
+#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/vector_memory.h>
#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/fe_lib.lagrange.h>
+#include <grid/grid_out.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
-#include <map>
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <cstdlib>
-
-
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation (const Function<dim> &rhs) :
- Equation<dim>(1),
- right_hand_side (rhs) {};
-
- virtual void assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (FullMatrix<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- protected:
- const Function<dim> &right_hand_side;
-};
-
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_dof_tools.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_smoother.h>
+#include <multigrid/multigrid.h>
+#include <lac/solver_richardson.h>
+#include <fstream>
template <int dim>
-class PoissonProblem {
+class LaplaceProblem
+{
public:
- /**
- * Declare a data type which denotes a
- * mapping between a boundary indicator
- * and the function denoting the boundary
- * values on this part of the boundary.
- * Only one boundary function may be given
- * for each boundary indicator, which is
- * guaranteed by the #map# data type.
- *
- * See the general documentation of this
- * class for more detail.
- */
- typedef map<unsigned char,const Function<dim>*> FunctionMap;
- /**
- * Typdedef an iterator which assembles
- * matrices and vectors.
- */
- typedef TriaActiveIterator<dim, Assembler<dim> > active_assemble_iterator;
-
- PoissonProblem (unsigned int order);
- virtual ~PoissonProblem ();
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+ void run ();
- void clear ();
- void create_new ();
+ private:
+ void setup_system ();
+ void assemble_system ();
void solve ();
-
- /**
- * Initiate the process of assemblage of
- * vectors and system matrix. Use the
- * given equation object and the given
- * quadrature formula. Also use the list
- * of dirichlet boundary value functions
- * (by default, no dirichlet bc are assumed
- * which means that all bc are included
- * into the weak formulation).
- *
- * 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 UpdateFlags update_flags,
- const FunctionMap &dirichlet_bc = FunctionMap());
-
- int run (unsigned int level);
- void print_history (string filename) const;
-
- /**
- * Exception
- */
- DeclException0 (ExcNoTriaSelected);
-
- protected:
- Triangulation<dim> *tria;
- MGDoFHandler<dim> *dof;
-
-
- /**
- * Sparsity pattern of the system matrix.
- */
- SparsityPattern system_sparsity;
-
- /**
- * System matrix.
- */
- SparseMatrix<double> system_matrix;
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
- /**
- * Vector storing the right hand side.
- */
- Vector<double> right_hand_side;
+ Triangulation<dim> triangulation;
+ MGDoFHandler<dim> mg_dof_handler;
- /**
- * Solution vector.
- */
- Vector<double> solution;
+ FEQ1<dim> fe;
- /**
- * List of constraints introduced by
- * hanging nodes.
- */
- ConstraintMatrix constraints;
+ ConstraintMatrix hanging_node_constraints;
- Function<dim> *rhs;
- Function<dim> *boundary_values;
+ SparsityPattern global_sparsity_pattern;
+ SparseMatrix<double> global_system_matrix;
- vector<double> l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error;
- vector<int> n_dofs;
-
- unsigned int order;
+ MGLevelObject<SparsityPattern> level_sparsity_patterns;
+ MGLevelObject<SparseMatrix<double> > level_system_matrices;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
};
-
-
-/**
- Right hand side constructed such that the exact solution is
- $sin(2 pi x) + sin(2 pi y)$
- */
template <int dim>
-class RHSPoly : public Function<dim> {
+class Coefficient : public Function<dim>
+{
public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int component = 0) const;
};
template <int dim>
-class Solution : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- /**
- * Return the gradient of the function
- * at the given point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component) const;
+double Coefficient<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+{
+ if (p.square() < 0.5*0.5)
+ return 20;
+ else
+ return 1;
};
+template <int dim>
+void Coefficient<dim>::value_list (const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int component) const
+{
+ const unsigned int n_points = points.size();
-template <>
-double RHSPoly<2>::value (const Point<2> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
+ Assert (values.size() == n_points,
+ ExcVectorHasWrongSize (values.size(), n_points));
+
+ Assert (component == 0,
+ ExcWrongComponent (component, 1));
- const double x = p(0),
- y = p(1);
- const double pi= 3.1415926536;
- return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y));
+ for (unsigned int i=0; i<n_points; ++i)
+ {
+ if (points[i].square() < 0.5*0.5)
+ values[i] = 20;
+ else
+ values[i] = 1;
+ };
};
-template <>
-double Solution<2>::value (const Point<2> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
- const double x = p(0),
- y = p(1);
- const double pi= 3.1415926536;
- return sin(2*pi*x)+sin(2*pi*y);
+class MGSmootherLAC : public MGSmootherBase
+{
+ private:
+ SmartPointer<MGLevelObject<SparseMatrix<double> > >matrices;
+ public:
+ MGSmootherLAC(MGLevelObject<SparseMatrix<double> >&);
+
+ virtual void smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
};
-template <>
-Tensor<1,2> Solution<2>::gradient (const Point<2> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
+MGSmootherLAC::MGSmootherLAC(MGLevelObject<SparseMatrix<double> >& matrix)
+ :
+ matrices(&matrix)
+{}
- const double x = p(0),
- y = p(1);
- const double pi= 3.1415926536;
- return Point<2> (2*pi*cos(2*pi*x),
- 2*pi*cos(2*pi*y));
-};
-
+void
+MGSmootherLAC::smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const
+{
+ SolverControl control(2,1.e-300,false,false);
+ PrimitiveVectorMemory<> mem;
+ SolverRichardson<> rich(control, mem);
+ PreconditionRelaxation<>
+ prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
+ rich.solve((*matrices)[level], u, rhs, prec);
+}
-template <>
-void PoissonEquation<2>::assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<2> &fe_values,
- const DoFHandler<2>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side.value(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+ mg_dof_handler (triangulation)
+{};
template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, ExcPureVirtualFunctionCalled());
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+ mg_dof_handler.clear ();
};
template <int dim>
-void PoissonEquation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, ExcPureVirtualFunctionCalled());
-};
-
+void LaplaceProblem<dim>::setup_system ()
+{
+ mg_dof_handler.distribute_dofs (fe);
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (mg_dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+ global_sparsity_pattern.reinit (mg_dof_handler.DoFHandler<dim>::n_dofs(),
+ mg_dof_handler.DoFHandler<dim>::n_dofs(),
+ mg_dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (mg_dof_handler, global_sparsity_pattern);
+ hanging_node_constraints.condense (global_sparsity_pattern);
+ global_sparsity_pattern.compress();
+ global_system_matrix.reinit (global_sparsity_pattern);
+ solution.reinit (mg_dof_handler.DoFHandler<dim>::n_dofs());
+ system_rhs.reinit (mg_dof_handler.DoFHandler<dim>::n_dofs());
+ const unsigned int n_levels = triangulation.n_levels();
+ level_system_matrices.resize (0, n_levels);
+ level_sparsity_patterns.resize (0, n_levels);
+
+ for (unsigned int level=0; level<n_levels; ++level)
+ {
+ level_sparsity_patterns[level].reinit (mg_dof_handler.n_dofs(level),
+ mg_dof_handler.n_dofs(level),
+ mg_dof_handler.max_couplings_between_dofs()); //xxx
+ MGDoFTools::make_sparsity_pattern (mg_dof_handler,
+ level,
+ level_sparsity_patterns[level]);
+ level_sparsity_patterns[level].compress();
+
+ level_system_matrices[level].reinit (level_sparsity_patterns[level]);
+ };
+};
template <int dim>
-PoissonProblem<dim>::PoissonProblem (unsigned int order) :
- tria(0), dof(0), rhs(0),
- boundary_values(0), order(order) {};
+void LaplaceProblem<dim>::assemble_system ()
+{
+ const Coefficient<dim> coefficient;
+ QGauss2<dim> quadrature_formula;
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
-template <int dim>
-PoissonProblem<dim>::~PoissonProblem ()
-{
- clear ();
-};
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+ vector<int> local_dof_indices (dofs_per_cell);
-template <int dim>
-void PoissonProblem<dim>::clear () {
- if (dof != 0) {
- delete dof;
- dof = 0;
- };
-
- if (tria != 0) {
- delete tria;
- tria = 0;
- };
-
- if (rhs != 0)
- {
- delete rhs;
- rhs = 0;
- };
+ // FIX
+ vector<double> coefficient_values (n_q_points, 1.0);
- if (boundary_values != 0)
+ // not only active cells
+ MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+ endc = mg_dof_handler.end();
+ for (; cell!=endc; ++cell)
{
- delete boundary_values;
- boundary_values = 0;
+ cell_matrix.clear ();
+ cell_rhs.clear ();
+
+ fe_values.reinit (cell);
+ const FullMatrix<double>
+ & shape_values = fe_values.get_shape_values();
+ const vector<vector<Tensor<1,dim> > >
+ & shape_grads = fe_values.get_shape_grads();
+ const vector<double>
+ & JxW_values = fe_values.get_JxW_values();
+ const vector<Point<dim> >
+ & q_points = fe_values.get_quadrature_points();
+
+ // FIX
+// coefficient.value_list (q_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)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (coefficient_values[q_point] *
+ (shape_grads[i][q_point] *
+ shape_grads[j][q_point] +
+ shape_values(i,q_point) *
+ shape_values(j,q_point) ) *
+ JxW_values[q_point]);
+
+ cell_rhs(i) += (shape_values (i,q_point) *
+ sin(4*sqrt(q_points[q_point].square())) *
+ fe_values.JxW (q_point));
+ };
+
+
+ cell->get_mg_dof_indices (local_dof_indices);
+ const unsigned int level = cell->level();
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ level_system_matrices[level].add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ // if active, then also into
+ // global matrix
+ if (cell->active())
+ {
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ global_system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ };
+ };
};
- system_sparsity.reinit (0,0,1);
- system_matrix.clear ();
- right_hand_side.reinit (1);
- solution.reinit (1);
- constraints.clear ();
+ hanging_node_constraints.condense (global_system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+
+// map<int,double> boundary_values;
+// VectorTools::interpolate_boundary_values (mg_dof_handler,
+// 0,
+// ZeroFunction<dim>(),
+// boundary_values);
+// MatrixTools<dim>::apply_boundary_values (boundary_values,
+// global_system_matrix,
+// solution,
+// system_rhs);
};
-
template <int dim>
-void PoissonProblem<dim>::create_new () {
- clear ();
-
- tria = new Triangulation<dim>();
- dof = new MGDoFHandler<dim> (*tria);
-};
+void LaplaceProblem<dim>::solve ()
+{
+ {
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+ SolverControl coarse_grid_solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> coarse_grid_vector_memory;
+
+ SolverCG<> coarse_grid_cg (coarse_grid_solver_control,
+ coarse_grid_vector_memory);
+
+// PreconditionRelaxation<>
+// coarse_grid_solver_preconditioner(level_system_matrices[level_system_matrices.get_minlevel()],
+// &SparseMatrix<double>::template precondition_SSOR<double>,
+// 1.2);
+ PreconditionIdentity coarse_grid_solver_preconditioner;
+
+ MGCoarseGridLACIteration<SolverCG<>, SparseMatrix<double>, PreconditionIdentity>
+ coarse_grid_solver (coarse_grid_cg,
+ level_system_matrices[level_system_matrices.get_minlevel()],
+ coarse_grid_solver_preconditioner);
+
+ MGSmootherLAC smoother (level_system_matrices);
+ MGTransferPrebuilt grid_transfer;
+ grid_transfer.build_matrices (mg_dof_handler);
+
+ Multigrid<2> multigrid (mg_dof_handler,
+ hanging_node_constraints,
+ level_sparsity_patterns,
+ level_system_matrices,
+ grid_transfer);
+
+ PreconditionMG<Multigrid<2> >
+ mg_precondition (multigrid, smoother, smoother, coarse_grid_solver);
-template <int dim>
-void PoissonProblem<dim>::assemble (const Equation<dim> &equation,
- const Quadrature<dim> &quadrature,
- const UpdateFlags update_flags,
- const FunctionMap &dirichlet_bc) {
- Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected());
-
- system_sparsity.reinit (dof->DoFHandler<dim>::n_dofs(),
- dof->DoFHandler<dim>::n_dofs(),
- dof->max_couplings_between_dofs());
- right_hand_side.reinit (dof->DoFHandler<dim>::n_dofs());
-
- // make up sparsity pattern and
- // compress with constraints
- constraints.clear ();
- DoFTools::make_hanging_node_constraints (*dof, constraints);
- constraints.close ();
- DoFTools::make_sparsity_pattern (*dof, system_sparsity);
- constraints.condense (system_sparsity);
-
- MGTransferPrebuilt p;
- p.build_matrices (*dof);
-// MGSmoother smoother(*dof);
-
+ solution.clear ();
+ cg.solve (global_system_matrix, solution, system_rhs,
+ mg_precondition);
- // reinite system matrix
- system_matrix.reinit (system_sparsity);
- // reinit solution vector, preset
- // with zeroes.
- solution.reinit (dof->DoFHandler<dim>::n_dofs());
+ cout << " MG Outer iterations: " << solver_control.last_step()
+ << endl;
+
+ cout << " MG Total inner iterations: " << coarse_grid_solver_control.last_step()
+ << endl;
+ };
- // create assembler
- Assembler<dim>::AssemblerData data (*dof,
- true, true, //assemble matrix and rhs
- system_matrix,
- right_hand_side,
- quadrature,
- update_flags);
- active_assemble_iterator assembler (tria,
- tria->begin_active()->level(),
- tria->begin_active()->index(),
- &data);
- // loop over all cells, fill matrix and rhs
- do
{
- assembler->assemble (equation);
- }
- while ((++assembler).state() == valid);
-
- // condense system matrix in-place
- constraints.condense (system_matrix);
-
- // condense right hand side in-place
- constraints.condense (right_hand_side);
-
- // apply Dirichlet bc as described
- // in the docs
- map<int, double> boundary_value_list;
-
- for (FunctionMap::const_iterator bc=dirichlet_bc.begin();
- bc != dirichlet_bc.end(); ++bc)
- VectorTools::interpolate_boundary_values (*dof,
- bc->first,
- *bc->second,
- boundary_value_list);
- MatrixTools<dim>::apply_boundary_values (boundary_value_list,
- system_matrix, solution,
- right_hand_side);
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+
+ PreconditionRelaxation<>
+ preconditioner(global_system_matrix,
+ &SparseMatrix<double>::template precondition_SSOR<double>,
+ 1.2);
+
+ solution.clear ();
+ cg.solve (global_system_matrix, solution, system_rhs,
+ preconditioner);
+
+ cout << " CG Outer iterations: " << solver_control.last_step()
+ << endl;
+ };
+
+ hanging_node_constraints.distribute (solution);
+};
+
+
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (mg_dof_handler,
+ QGauss3<dim-1>(),
+ neumann_boundary,
+ solution,
+ estimated_error_per_cell);
+
+ triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+ 0.3, 0.03);
+ triangulation.execute_coarsening_and_refinement ();
};
template <int dim>
-void PoissonProblem<dim>::solve () {
- Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected());
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+ string filename = "grid-";
+ filename += ('0' + cycle);
+ Assert (cycle < 10, ExcInternalError());
- SolverControl control(4000, 1e-16);
- PrimitiveVectorMemory<Vector<double> > memory;
- SolverCG<SparseMatrix<double>,Vector<double> > cg(control,memory);
-
- // solve
- cg.solve (system_matrix, solution, right_hand_side,
- PreconditionIdentity());
- // distribute solution
- constraints.distribute (solution);
+ filename += ".eps";
+ ofstream output (filename.c_str());
+
+ GridOut grid_out;
+ grid_out.write_eps (triangulation, output);
};
template <int dim>
-int PoissonProblem<dim>::run (const unsigned int level) {
- create_new ();
-
- cout << "Refinement level = " << level
- << ", using elements of type <";
- switch (order)
+void LaplaceProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<8; ++cycle)
{
- case 0:
- cout << "criss-cross";
- break;
- default:
- cout << "Lagrange-" << order;
- break;
- };
- cout << ">" << endl;
-
- cout << " Making grid... ";
- GridGenerator::hyper_cube (*tria);
- tria->refine_global (level+1);
- tria->begin_active()->set_refine_flag();
- (++(++(tria->begin_active())))->set_refine_flag();
- tria->execute_coarsening_and_refinement ();
- cout << tria->n_active_cells() << " active cells." << endl;
-
- rhs = new RHSPoly<dim>();
- boundary_values = new Solution<dim> ();
-
+ cout << "Cycle " << cycle << ':' << endl;
- FiniteElement<dim> *fe;
- PoissonEquation<dim> equation (*rhs);
- Quadrature<dim> *quadrature;
- Quadrature<dim-1> *boundary_quadrature;
- switch (order) {
- case 0:
- fe = new FECrissCross<dim>();
- quadrature = new QCrissCross1<dim>();
- boundary_quadrature = new QGauss2<dim-1>();
- break;
- case 1:
- fe = new FEQ1<dim>();
- quadrature = new QGauss3<dim>();
- boundary_quadrature = new QGauss2<dim-1>();
- break;
- case 2:
- fe = new FEQ2<dim>();
- quadrature = new QGauss4<dim>();
- boundary_quadrature = new QGauss3<dim-1>();
- break;
- case 3:
- fe = new FEQ3<dim>();
- quadrature = new QGauss5<dim>();
- boundary_quadrature = new QGauss4<dim-1>();
- break;
- case 4:
- fe = new FEQ4<dim>();
- quadrature = new QGauss6<dim>();
- boundary_quadrature = new QGauss5<dim-1>();
- break;
- default:
- return 100000;
- };
-
- cout << " Distributing dofs... ";
- dof->distribute_dofs (*fe);
- cout << dof->DoFHandler<dim>::n_dofs() << " degrees of freedom." << endl;
- n_dofs.push_back (dof->DoFHandler<dim>::n_dofs());
-
- cout << " Assembling matrices..." << endl;
- UpdateFlags update_flags = UpdateFlags(update_values | update_q_points |
- update_gradients | update_JxW_values);
-
- map<unsigned char,const Function<dim>*> dirichlet_bc;
- dirichlet_bc[0] = boundary_values;
- assemble (equation, *quadrature, update_flags, dirichlet_bc);
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (1);
+ }
+ else
+ {
+ refine_grid ();
+ };
+
- cout << " Solving..." << endl;
- solve ();
+ cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << endl;
- Solution<dim> sol;
- Vector<float> l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell;
- Vector<float> h1_seminorm_error_per_cell, h1_error_per_cell;
-
- cout << " Calculating L1 error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- l1_error_per_cell,
- *quadrature, L1_norm);
- cout << l1_error_per_cell.l1_norm() << endl;
- l1_error.push_back (l1_error_per_cell.l1_norm());
-
- cout << " Calculating L2 error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- l2_error_per_cell,
- *quadrature, L2_norm);
- cout << l2_error_per_cell.l2_norm() << endl;
- l2_error.push_back (l2_error_per_cell.l2_norm());
-
- cout << " Calculating L-infinity error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- linfty_error_per_cell,
- *quadrature, Linfty_norm);
- cout << linfty_error_per_cell.linfty_norm() << endl;
- linfty_error.push_back (linfty_error_per_cell.linfty_norm());
-
- cout << " Calculating H1-seminorm error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- h1_seminorm_error_per_cell,
- *quadrature, H1_seminorm);
- cout << h1_seminorm_error_per_cell.l2_norm() << endl;
- h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm());
-
- cout << " Calculating H1 error... ";
- VectorTools::integrate_difference (*dof,
- solution, sol,
- h1_error_per_cell,
- *quadrature, H1_norm);
- cout << h1_error_per_cell.l2_norm() << endl;
- h1_error.push_back (h1_error_per_cell.l2_norm());
-
- if (dof->DoFHandler<dim>::n_dofs()<=5000)
- {
- Vector<double> l1_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> l2_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> linfty_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> h1_seminorm_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- Vector<double> h1_error_per_dof (dof->DoFHandler<dim>::n_dofs());
- DoFTools::distribute_cell_to_dof_vector (*dof, l1_error_per_cell, l1_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, l2_error_per_cell, l2_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, linfty_error_per_cell,
- linfty_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, h1_seminorm_error_per_cell,
- h1_seminorm_error_per_dof);
- DoFTools::distribute_cell_to_dof_vector (*dof, h1_error_per_cell, h1_error_per_dof);
-
-// Vector<double> projected_solution;
-// ConstraintMatrix constraints;
-// constraints.close ();
-// VectorTools::project (*dof, constraints, *fe,
-// StraightBoundary<dim>(), *quadrature,
-// sol, projected_solution, false,
-// *boundary_quadrature);
-// cout << " Calculating L2 error of projected solution... ";
-// VectorTools::integrate_difference (*dof,
-// projected_solution, sol,
-// l2_error_per_cell,
-// *quadrature, *fe, L2_norm);
-// cout << l2_error_per_cell.l2_norm() << endl;
-
-
- string filename;
- filename = ('0'+order);
- filename += ".";
- filename += ('0'+level);
- filename += ".ucd";
- cout << " Writing error plots to <" << filename << ">..." << endl;
+ setup_system ();
+
+ cout << " Number of degrees of freedom: "
+ << mg_dof_handler.DoFHandler<dim>::n_dofs()
+ << endl;
- DataOut<dim> out;
- ofstream o(filename.c_str());
- out.attach_dof_handler (*dof);
-
- out.add_data_vector (solution, "u");
- out.add_data_vector (l1_error_per_dof, "L1_Error");
- out.add_data_vector (l2_error_per_dof, "L2_Error");
- out.add_data_vector (linfty_error_per_dof, "Linfty_Error");
- out.add_data_vector (h1_seminorm_error_per_dof, "H1_seminorm_Error");
- out.add_data_vector (h1_error_per_dof, "H1_Error");
- out.write_ucd (o);
- o.close ();
- }
- else
- cout << " Not writing error as grid." << endl;
-
- cout << endl;
-
- const unsigned int n_dofs = dof->DoFHandler<dim>::n_dofs();
- dof->clear ();
- tria->set_boundary (0);
- delete fe;
- delete quadrature;
- delete boundary_quadrature;
-
- return n_dofs;
-};
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ DataOut<dim>::EpsFlags eps_flags;
+ eps_flags.z_scaling = 4;
+
+ DataOut<dim> data_out;
+ data_out.set_flags (eps_flags);
-template <int dim>
-void PoissonProblem<dim>::print_history (string filename) const {
- ofstream out(filename.c_str());
- out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error"
- << endl;
- for (unsigned int i=0; i<n_dofs.size(); ++i)
- out << n_dofs[i]
- << " "
- << l1_error[i] << " "
- << l2_error[i] << " "
- << linfty_error[i] << " "
- << h1_seminorm_error[i] << " "
- << h1_error[i] << endl;
-
- double average_l1=0,
- average_l2=0,
- average_linfty=0,
- average_h1_semi=0,
- average_h1=0;
- for (unsigned int i=1; i<n_dofs.size(); ++i)
- {
- average_l1 += l1_error[i]/l1_error[i-1];
- average_l2 += l2_error[i]/l2_error[i-1];
- average_linfty += linfty_error[i]/linfty_error[i-1];
- average_h1_semi += h1_seminorm_error[i]/h1_seminorm_error[i-1];
- average_h1 += h1_error[i]/h1_error[i-1];
+ data_out.attach_dof_handler (mg_dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ ofstream output ("final-solution.eps");
+ data_out.write_eps (output);
};
-
- average_l1 /= (l1_error.size()-1);
- average_l2 /= (l1_error.size()-1);
- average_linfty /= (l1_error.size()-1);
- average_h1_semi /= (l1_error.size()-1);
- average_h1 /= (l1_error.size()-1);
-
- cout << "==========================================================\n";
- cout << "Average error reduction rates for h->h/2:" << endl;
- cout << " L1 error : " << 1./average_l1 << endl
- << " L2 error : " << 1./average_l2 << endl
- << " Linfty error : " << 1./average_linfty << endl
- << " H1 seminorm error: " << 1./average_h1_semi << endl
- << " H1 error : " << 1./average_h1 << endl;
- cout << "==========================================================\n";
- cout << "==========================================================\n";
};
+
+int main ()
+{
+ try
+ {
+ deallog.depth_console (0);
-
-int main () {
- for (unsigned int order=0; order<5; ++order)
+ LaplaceProblem<2> laplace_problem_2d;
+ laplace_problem_2d.run ();
+ }
+ catch (exception &exc)
{
- PoissonProblem<2> problem (order);
-
- unsigned int level=0;
- unsigned int n_dofs;
- do
- n_dofs = problem.run (level++);
- while (n_dofs<25000);
-
- string filename;
- switch (order)
- {
- case 0:
- filename = "criss_cross";
- break;
- case 1:
- filename = "linear";
- break;
- case 2:
- filename = "quadratic";
- break;
- case 3:
- filename = "cubic";
- break;
- case 4:
- filename = "quartic";
- break;
- };
- filename += ".history";
-
- cout << endl << "Printing convergence history to <"
- << filename << ">..." << endl;
- problem.print_history (filename);
- cout << endl << endl << endl;
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Exception on processing: " << endl
+ << exc.what() << endl
+ << "Aborting!" << endl
+ << "----------------------------------------------------"
+ << endl;
+ return 1;
+ }
+ catch (...)
+ {
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Unknown exception!" << endl
+ << "Aborting!" << endl
+ << "----------------------------------------------------"
+ << endl;
+ return 1;
};
-
+
return 0;
};
// $Id$
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
/**
* Generator for system matrix and right hand side.
const Function<dim>& rhs);
template<int dim, class MATRIX>
- void build_mgmatrix(MGMatrix<MATRIX>& A,
+ void build_mgmatrix(MGLevelObject<MATRIX>& A,
const MGDoFHandler<dim>& dof,
const Quadrature<dim>& quadrature);
};
template<int dim, class MATRIX>
void
-Helmholtz::build_mgmatrix(MGMatrix<MATRIX>& A,
+Helmholtz::build_mgmatrix(MGLevelObject<MATRIX>& A,
const MGDoFHandler<dim>& dof,
const Quadrature<dim>& quadrature)
{
#include <lac/precondition.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
-#include <dofs/mg_dof_handler.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_accessor.h>
-#include <dofs/mg_dof_accessor.h>
#include <grid/grid_generator.h>
#include <fe/fe_lib.lagrange.h>
#include <fe/fe_values.h>
-#include <numerics/multigrid.h>
-#include <numerics/multigrid.templates.h>
-#include <numerics/mg_smoother.h>
#include <dofs/dof_tools.h>
-#include <dofs/mg_dof_tools.h>
+
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_dof_tools.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_smoother.h>
+#include <multigrid/multigrid.h>
#include "helmholtz.h"
+#include <fstream>
+
+
+
+
template<int dim>
class RHSFunction : public Function<dim>
{
class MGSmootherLAC : public MGSmootherBase
{
private:
- SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
+ SmartPointer<MGLevelObject<SparseMatrix<double> > >matrices;
public:
- MGSmootherLAC(MGMatrix<SparseMatrix<double> >&);
+ MGSmootherLAC(MGLevelObject<SparseMatrix<double> >&);
virtual void smooth (const unsigned int level,
Vector<double> &u,
solver.solve(A,u,f,precondition);
u = 0.;
- vector<SparsityPattern> mgstruct(tr.n_levels());
- MGMatrix<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
+ MGLevelObject<SparsityPattern> mgstruct(0, tr.n_levels()-1);
+ MGLevelObject<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
for (unsigned int i=0;i<tr.n_levels();++i)
{
mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
transfer.build_matrices(mgdof);
- MG<2> multigrid(mgdof, hanging_nodes, mgA, transfer);
- PreconditionMG<MG<2> >
+ Multigrid<2> multigrid(mgdof, hanging_nodes, mgstruct, mgA, transfer);
+ PreconditionMG<Multigrid<2> >
mgprecondition(multigrid, smoother, smoother, coarse);
solver.solve(A, u, f, mgprecondition);
return 1.;
}
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<double> >& matrix)
+MGSmootherLAC::MGSmootherLAC(MGLevelObject<SparseMatrix<double> >& matrix)
:
matrices(&matrix)
{}
#include <lac/precondition.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
-#include <dofs/mg_dof_handler.h>
+#include <multigrid/mg_dof_handler.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_accessor.h>
-#include <dofs/mg_dof_accessor.h>
+#include <multigrid/mg_dof_accessor.h>
#include <grid/grid_generator.h>
#include <numerics/data_io.h>
#include <fe/fe_lib.lagrange.h>
#include <fe/fe_values.h>
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
-#include <numerics/multigrid.templates.h>
+#include <multigrid/multigrid.h>
+#include <multigrid/mg_smoother.h>
#include <dofs/dof_tools.h>
-#include <dofs/mg_dof_tools.h>
+#include <multigrid/mg_dof_tools.h>
#include <fstream>
SolverControl control(20, 1.e-12, true);
SolverRichardson<> solver(control, mem);
- vector<SparsityPattern> mgstruct(tr.n_levels());
- MGMatrix<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
+ MGLevelObject<SparsityPattern> mgstruct(0, tr.n_levels()-1);
+ MGLevelObject<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
for (unsigned int i=0;i<tr.n_levels();++i)
{
mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
transfer.build_matrices(mgdof);
- MG<2> multigrid(mgdof, hanging_nodes, mgA, transfer, tr.n_levels()-2);
- PreconditionMG<MG<2> >
+ Multigrid<2> multigrid(mgdof, hanging_nodes, mgstruct, mgA, transfer, tr.n_levels()-2);
+ PreconditionMG<Multigrid<2> >
mgprecondition(multigrid, smoother, smoother, coarse);
u = 0.;
############################################################
libs.g = $(lib-lac.g) \
- $(lib-base.g)
+ $(lib-base.g) \
+ $(lib-deal2-2d.g)
libs = $(lib-lac.o) \
- $(lib-base.o)
+ $(lib-base.o) \
+ $(lib-deal2-2d.o)
############################################################
#include <lac/solver_bicgstab.h>
#include <lac/solver_richardson.h>
#include <lac/precondition.h>
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_smoother.h>
+
#define TYPE double
#define ACCURACY 1.e-5
/**
* Pointer to the level matrices.
*/
- SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
+ SmartPointer<MGPerLevelObject<SparseMatrix<double> > >matrices;
public:
- FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<double> >& matrices,
+ FDMG(unsigned int maxlevel, MGPerLevelObject<SparseMatrix<double> >& matrices,
FDMGTransfer& transfer);
public MGSmootherBase
{
private:
- SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
+ SmartPointer<MGPerLevelObject<SparseMatrix<double> > >matrices;
public:
- MGSmootherLAC(MGMatrix<SparseMatrix<double> >&);
+ MGSmootherLAC(MGPerLevelObject<SparseMatrix<double> >&);
virtual void smooth (const unsigned int level,
Vector<double> &u,
// Make matrix
vector<SparsityPattern> structure(maxlevel+1);
- MGMatrix<SparseMatrix<double> > A(minlevel,maxlevel);
+ MGPerLevelObject<SparseMatrix<double> > A(minlevel,maxlevel);
FDMatrix testproblem(size, size);
}
}
-FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<double> >& matrices,
+FDMG::FDMG(unsigned int maxlevel, MGPerLevelObject<SparseMatrix<double> >& matrices,
FDMGTransfer& transfer)
:
MGBase(transfer, 0, maxlevel),
}
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<double> >& matrix)
+MGSmootherLAC::MGSmootherLAC(MGPerLevelObject<SparseMatrix<double> >& matrix)
:
matrices(&matrix)
{}
// deal_II_libraries.g=-ldeal_II_2d.g
// deal_II_libraries=-ldeal_II_2d
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_smoother.h>
#include <base/logstream.h>
#include <fstream>
#include <lac/forward_declarations.h>
#include <base/exceptions.h>
-#include <lac/mgbase.h>
+#include <multigrid/mg_base.h>
/**
* Finite difference matrix on uniform grid.
projects / dealii-svn.git / commitdiff