#include <grid/tria_accessor.h>
-#include <lac/forward_declarations.h>
#include <vector>
+template <typename number> class FullMatrix;
+template <typename number> class SparseMatrix;
+template <typename number> class Vector;
+
+template <int dim> class DoFHandler;
// note: in non-debug mode, i.e. with optimizations, the file
// dof_accessor.templates.h is included at the end of this file.
#include <utility>
#include <base/exceptions.h>
#include <base/subscriptor.h>
-#include <lac/forward_declarations.h>
+class SparsityPattern;
+template <int rows, int columns> class BlockSparsityPattern;
+template <typename number> class SparseMatrix;
+template <typename number, int rows, int columns> class BlockSparseMatrix;
+template <int n_blocks> class BlockIndices;
/**
* This class represents the matrix denoting the distribution of the degrees
#include <vector>
#include <map>
#include <base/exceptions.h>
-#include <base/forward_declarations.h>
#include <base/smartpointer.h>
-#include <grid/forward_declarations.h>
-#include <lac/forward_declarations.h>
+template <int dim> class Function;
+
+template <int dim> class DoFCellAccessor;
+template <int dim> class DoFLevel;
+template<int celldim, int dim> class DoFObjectAccessor;
+template <int dim> class FiniteElement;
+template <int dim, typename Accessor> class TriaRawIterator;
+template <int dim, typename Accessor> class TriaIterator;
+template <int dim, typename Accessor> class TriaActiveIterator;
+template <int dim> class Triangulation;
/**
* Define some types which differ between the dimensions. This class
#include <base/exceptions.h>
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
#include <vector>
+template <int dim> class MGDoFHandler;
+template <template <int> class GridClass, int dim> class InterGridMap;
/**
* This is a collection of functions operating on, and manipulating
#include <grid/tria.h>
#include <fe/fe_update_flags.h>
+template <int dim> class Quadrature;
+
/**
* This class offers a multitude of arrays and other fields which are used by
#define __deal2__grid_generator_h
-#include <base/forward_declarations.h>
-#include <grid/forward_declarations.h>
#include <base/exceptions.h>
+template <int dim> class Point;
+template <int dim> class Triangulation;
+
/**
* This class offers triangulations of some standard domains such as hypercubes,
#include <base/exceptions.h>
#include <base/smartpointer.h>
-#include <grid/forward_declarations.h>
#include <iostream>
#include <vector>
#include <string>
+template <int dim> class Triangulation;
/**
* This class implements an input mechanism for grid data. It allows to
#define __deal2__grid_out_h
-#include <grid/forward_declarations.h>
#include <base/exceptions.h>
#include <string>
+template <int dim> class Triangulation;
/**
* This class provides a means to output a triangulation to a file in different
#ifndef __deal2__intergrid_map_h
#define __deal2__intergrid_map_h
+template <class T> class SmartPointer;
/**
* This class provides a map between two grids which are derived from
#include <vector>
-#include <grid/forward_declarations.h>
-#include <lac/forward_declarations.h>
-#include <base/forward_declarations.h>
#include <base/point.h>
#include <grid/geometry_info.h>
#include <base/subscriptor.h>
+template <int dim> class Boundary;
+template <int dim> class StraightBoundary;
+template <int dim> class CellAccessor;
+template<int celldim, int dim> class TriaObjectAccessor;
+template <int dim, typename Accessor> class TriaRawIterator;
+template <int dim, typename Accessor> class TriaIterator;
+template <int dim, typename Accessor> class TriaActiveIterator;
+template <int dim> class TriangulationLevel;
+template <int dim> class DoFHandler;
+template <int dim> class MGDoFHandler;
+
/*------------------------------------------------------------------------*/
#include <base/exceptions.h>
-#include <base/forward_declarations.h>
-#include <grid/forward_declarations.h>
#include <grid/tria_iterator_base.h>
+template <int dim> class Point;
+
+template <int dim> class Triangulation;
+template <int dim, typename Accessor> class TriaRawIterator;
+template <int dim, typename Accessor> class TriaIterator;
+template <int dim, typename Accessor> class TriaActiveIterator;
+class Line;
+class Quad;
+class Hexahedron;
+
// note: in non-debug mode, i.e. with optimizations, the file
// tria_accessor.templates.h is included at the end of this file.
present_index (index),
tria (parent) {};
-
/**
* Copy operator. Since this is only
* called from iterators, do not
#include <base/point.h>
#include <base/subscriptor.h>
#include <grid/geometry_info.h>
-#include <grid/forward_declarations.h>
+template <int dim> class Triangulation;
/**
* This class is used to represent a boundary to a triangulation.
#include <base/exceptions.h>
#include <base/point.h>
#include <grid/tria_iterator_base.h>
-#include <grid/forward_declarations.h>
+template <int dim> class Triangulation;
// note: in non-debug mode, i.e. with optimizations, the file
// tria_iterator.templates.h is included at the end of this file.
#include <base/smartpointer.h>
#include <lac/vector.h>
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
-
#include <vector>
+class MGSmootherBase;
+
/**
* An array with an object for each level. The purpose of this class
#include <dofs/dof_accessor.h>
+template <int dim> class MGDoFHandler;
/**
* Define the basis for accessors to the degrees of freedom for
#include <dofs/dof_handler.h>
+template <int dim> class MGDoFCellAccessor;
+template <int celldim, int dim> class MGDoFObjectAccessor;
/**
* Define some types which differ between the dimensions. This class
#define __deal2__mg_dof_tools_h
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
-
//TODO: Consider incorporating these functions in DoFTools
/**
#include <lac/forward_declarations.h>
#include <multigrid/mg_base.h>
-#include <grid/forward_declarations.h>
#include <base/smartpointer.h>
#include <vector>
+template <int dim> class MGDoFHandler;
/**
* Abstract base class for multigrid smoothers.
#include <map>
#include <string>
+template <int dim> class DataOut;
+
/**
* The use of this class is now deprecated!
#include <base/exceptions.h>
-#include <grid/forward_declarations.h>
-#include <lac/forward_declarations.h>
#include <vector>
#include <string>
+template <typename number> class Vector;
+template <int dim> class DoFHandler;
+class EpsOutputData;
/**
* This class is deprecated. Use the #DataOut# class instead.
#define __deal2__data_out_h
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
#include <base/data_out_base.h>
#include <base/multithread_info.h>
+template <int dim> class DoFHandler;
/**
* This is an abstract class which provides the functionality to generate
#include <base/data_out_base.h>
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
#include <base/smartpointer.h>
#include <string>
#include <vector>
+template <int dim> class DoFHandler;
/**
* This class is used to stack the output from several computations
#include <base/exceptions.h>
-#include <grid/forward_declarations.h>
/**
#include <base/exceptions.h>
#include <base/function.h>
#include <base/multithread_info.h>
-#include <grid/forward_declarations.h>
#include <dofs/dof_handler.h>
#include <map>
#include <base/exceptions.h>
#include <map>
+template <int dim> class Quadrature;
+
+template<typename number> class Vector;
+template<typename number> class FullMatrix;
+template<typename number> class SparseMatrix;
+
+template <typename number, int rows, int columns > class BlockSparseMatrix;
+template <int n_blocks, typename Number> class BlockVector;
+template <int dim> class FEValues;
+template <int dim> class Equation;
/**
* Provide a class which assembles certain standard matrices for a given
#define __deal2__vectors_h
-#include <base/forward_declarations.h>
#include <base/exceptions.h>
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
#include <map>
#include <vector>
+template <int dim> class Function;
+template <int dim> class Quadrature;
+template <int dim> class QGauss2;
+
+template <typename number> class Vector;
+template <typename number> class FullMatrix;
+template <int dim> class DoFHandler;
+class ConstraintMatrix;
+
/**
* Denote which norm/integral is to be computed. The following possibilities
const MGSmootherBase &post_smooth,
const MGCoarseGridSolver &coarse_grid_solver)
{
+#ifdef MG_DEBUG
char *name = new char[100];
-
sprintf(name, "MG%d-defect",level);
print_vector(level, defect[level], name);
-
+#endif
+
solution[level] = 0.;
if (level == minlevel)
{
coarse_grid_solver(level, solution[level], defect[level]);
+#ifdef MG_DEBUG
sprintf(name, "MG%d-solution",level);
print_vector(level, solution[level], name);
+#endif
return;
}
pre_smooth.smooth(level, solution[level], defect[level]);
// t = d-As
+#ifdef MG_DEBUG
sprintf(name, "MG%d-pre",level);
print_vector(level, solution[level], name);
+#endif
t.reinit(solution[level].size());
level_vmult(level, t, solution[level], defect[level]);
transfer->prolongate(level, t, solution[level-1]);
+#ifdef MG_DEBUG
sprintf(name, "MG%d-cgc",level);
print_vector(level, t, name);
+#endif
solution[level] += t;
post_smooth.smooth(level, solution[level], defect[level]);
+#ifdef MG_DEBUG
sprintf(name, "MG%d-post",level);
print_vector(level, solution[level], name);
delete[] name;
+#endif
}
#include <lac/sparse_matrix.h>
#include <lac/block_sparsity_pattern.h>
-
+template <int n_blocks, typename Number> class BlockVector;
/**
* Blocked sparse matrix. The behaviour of objects of this type is
#include <base/exceptions.h>
#include <base/subscriptor.h>
-#include <lac/forward_declarations.h>
#include <lac/sparsity_pattern.h>
#include <lac/block_indices.h>
#define __deal2__eigen_h
-#include <lac/forward_declarations.h>
#include <lac/solver.h>
#include <lac/solver_control.h>
#include <lac/solver_cg.h>
#define __deal2__precondition_block_h
-#include <lac/forward_declarations.h>
#include <base/exceptions.h>
#include <base/subscriptor.h>
#include <base/smartpointer.h>
#include <vector>
+template <typename number> class FullMatrix;
+template <typename number> class SparseMatrix;
+template <typename number> class Vector;
+
/**
* Base class for #PreconditionBlockJacobi#, #PreconditionBlockSOR#, ...
};
+/**
+ * Block Jacobi preconditioning.
+ */
+template<typename number,
+ typename inverse_type = number>
+class PreconditionBlockJacobi : public PreconditionBlock<number,inverse_type>
+{
+ public:
+ /**
+ * Execute block Jacobi preconditioning.
+ * Make sure that the right
+ * block size
+ * of the matrix is set by
+ * #set_block_size#
+ * before calling this function.
+ *
+ * This function will automatically use the
+ * inverse matrices if they exist, if not
+ * then BlockJacobi will need much time
+ * inverting the diagonal block
+ * matrices in each preconditioning step.
+ */
+ template <typename number2>
+ void operator() (Vector<number2>&, const Vector<number2>&) const;
+};
+
+
/**
* Block SOR preconditioning.
*
}
+/*--------------------- PreconditionBlockJacobi -----------------------*/
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockJacobi<number,inverse_type>
+::operator() (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ Assert(A!=0, ExcNoMatrixGivenToUse());
+ const SparseMatrix<number> &M=*A;
+ Assert (M.m() == M.n(), ExcMatrixNotSquare());
+ Assert (blocksize!=0, ExcBlockSizeNotSet());
+ Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
+ const unsigned int n_cells=M.m()/blocksize;
+ Assert (inverse.size()==0 || inverse.size()==n_cells,
+ ExcWrongNumberOfInverses(inverse.size(), n_cells));
+
+ Vector<number2> b_cell(blocksize), x_cell(blocksize);
+
+ // cell_row, cell_column are the
+ // numbering of the blocks (cells).
+ // row_cell, column_cell are the local
+ // numbering of the unknowns in the
+ // blocks.
+ // row, column are the global numbering
+ // of the unkowns.
+ unsigned int row, row_cell, begin_diag_block=0;
+
+ if (inverse.size()==0)
+ {
+ FullMatrix<number> M_cell(blocksize);
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ {
+ b_cell(row_cell)=src(row);
+ for (unsigned int column_cell=0, column=cell*blocksize;
+ column_cell<blocksize; ++column_cell, ++column)
+ M_cell(row_cell,column_cell)=M(row,column);
+ }
+ M_cell.householder(b_cell);
+ M_cell.backward(x_cell,b_cell);
+ // distribute x_cell to dst
+ for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ dst(row)=x_cell(row_cell);
+
+ begin_diag_block+=blocksize;
+ }
+ }
+ else
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ {
+ b_cell(row_cell)=src(row);
+ }
+ inverse[cell].vmult(x_cell, b_cell);
+ // distribute x_cell to dst
+ for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ dst(row)=x_cell(row_cell);
+
+ begin_diag_block+=blocksize;
+ }
+}
/*--------------------- PreconditionBlockSOR -----------------------*/
template<typename number, typename inverse_type>
#ifndef __deal2__solver_h
#define __deal2__solver_h
-
-#include <lac/forward_declarations.h>
+template<typename number> class Vector;
+template<class VECTOR> class VectorMemory;
+class SolverControl;
/**
#include <lac/vector_memory.h>
#include <lac/solver_richardson.h>
#include <lac/precondition.h>
-#include <lac/forward_declarations.h>
/**
#include <base/smartpointer.h>
#include <lac/sparsity_pattern.h>
+template<typename number> class Vector;
/**
* Sparse matrix.
#include <base/smartpointer.h>
#include <base/multithread_info.h>
-#include <lac/forward_declarations.h>
#include <vector>
#include <map>
+template <typename number> class FullMatrix;
+template <typename number> class SparseMatrix;
+template <typename number> class Vector;
/**
* Point-wise Vanka preconditioning.
#include <base/exceptions.h>
#include <base/subscriptor.h>
-#include <lac/forward_declarations.h>
+
+template <typename number> class SparseMatrix;
#include <vector>
template class PreconditionBlock<double, double>;
+/*--------------------- PreconditionBlockJacobi -----------------------*/
+
+
+// explicit instantiations for "float" PreconditionBlock
+template class PreconditionBlockJacobi<float, float>;
+
+template void PreconditionBlockJacobi<float, float>::operator() (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::operator() (
+ Vector<double> &, const Vector<double> &) const;
+
+
+template class PreconditionBlockJacobi<double, float>;
+
+template void PreconditionBlockJacobi<double, float>::operator() (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::operator() (
+ Vector<double> &, const Vector<double> &) const;
+
+template class PreconditionBlockJacobi<double, double>;
+
+template void PreconditionBlockJacobi<double, double>::operator() (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::operator() (
+ Vector<double> &, const Vector<double> &) const;
+
/*--------------------- PreconditionBlockGaussSeidel -----------------------*/
// logfile.setf(ios::fixed);
// logfile.precision (3);
deallog.attach(logfile);
+ deallog.log_execution_time(true);
// deallog.depth_console(0);
Helmholtz equation;
FEQ2<2> fe2;
FEQ3<2> fe3;
FEQ4<2> fe4;
- for (unsigned int degree=1;degree<=4;degree++)
+ for (unsigned int degree=1;degree<=3;degree++)
{
Triangulation<2> tr;
MGDoFHandler<2> mgdof(tr);
cell->set_refine_flag();
tr.execute_coarsening_and_refinement();
- tr.refine_global(2);
+ tr.refine_global(1);
dof.distribute_dofs(*fe);
const unsigned int size = dof.n_dofs();
deallog << "DoFs " << size << endl;
deallog << "Levels: " << tr.n_levels() << endl;
- for (unsigned int step=14;step < 15; ++step)
+ for (unsigned int step=1;step < 3; ++step)
{
deallog << "smoothing-steps" << step << endl;
SparsityPattern structure(size, dof.max_couplings_between_dofs());
Vector<double> u;
u.reinit(f);
PrimitiveVectorMemory<> mem;
- SolverControl control(20, 1.e-8, true);
+
+ SolverControl control(20, 1.e-12, true);
SolverCG<> solver(control, mem);
MGLevelObject<SparsityPattern> mgstruct(0, tr.n_levels()-1);