#include <base/config.h>
#include <base/std_cxx1x/function.h>
#include <base/geometry_info.h>
-#include <base/named_data.h>
-#include <lac/block_indices.h>
-#include <base/mg_level_object.h>
+#include <lac/matrix_block.h>
+#include <lac/block_vector.h>
#include <numerics/mesh_worker_vector_selector.h>
-#include <numerics/mesh_worker_info.h>
-#include <numerics/mesh_worker_assembler.h>
-
DEAL_II_NAMESPACE_OPEN
+class BlockIndices;
template<int,int> class DoFHandler;
template<int,int> class MGDoFHandler;
*/
namespace MeshWorker
{
-}
+/**
+ * The class providing the scrapbook to fill with local integration
+ * results. These can be values, local contributions to forms or cell
+ * and face matrices.
+ *
+ * The local matrices initialized by reinit() of the info object and
+ * then assembled into the global system by Assembler classes.
+ *
+ * @ingroup MeshWorker
+ * @author Guido Kanschat, 2009
+ */
+ template <typename number>
+ class LocalResults
+ {
+ private:
+ /**
+ * Initialize a single local
+ * matrix block. A helper
+ * function for initialize()
+ */
+ void initialize_local(MatrixBlock<FullMatrix<number> >& M,
+ const unsigned int row,
+ const unsigned int col);
+
+ public:
+ void initialize_numbers(const unsigned int n);
+ void initialize_vectors(const unsigned int n);
+ /**
+ * Allocate @p n local
+ * matrices. Additionally,
+ * set their block row and
+ * column coordinates to
+ * zero. The matrices
+ * themselves are resized by
+ * reinit().
+ *
+ * The template parameter @p
+ * MatrixPtr should point to
+ * a MatrixBlock
+ * instantiation in order to
+ * provide row and column info.
+ */
+ void initialize_matrices(unsigned int n, bool both);
+
+ /**
+ * Allocate a local matrix
+ * for each of the global
+ * ones in @p
+ * matrices. Additionally,
+ * set their block row and
+ * column coordinates. The
+ * matrices themselves are
+ * resized by reinit().
+ *
+ * The template parameter @p
+ * MatrixPtr should point to
+ * a MatrixBlock
+ * instantiation in order to
+ * provide row and column info.
+ */
+ template <class MATRIX>
+ void initialize_matrices(const MatrixBlockVector<MATRIX>& matrices,
+ bool both);
+
+ /**
+ * Initialize quadrature values
+ * to <tt>nv</tt> values in
+ * <tt>np</tt> quadrature points.
+ */
+ void initialize_quadrature(unsigned int np, unsigned int nv);
+
+ /**
+ * Reinitialize matrices for
+ * new cell. Resizes the
+ * matrices for hp and sets
+ * them to zero.
+ */
+ void reinit(const BlockIndices& local_sizes);
+
+ /**
+ * The number of scalar values.
+ */
+ unsigned int n_values () const;
+
+ /**
+ * The number of vectors.
+ */
+ unsigned int n_vectors () const;
+
+ /**
+ * The number of matrices.
+ */
+ unsigned int n_matrices () const;
+
+ /**
+ * The number of quadrature
+ * points in #quadrature_values.
+ */
+ unsigned int n_quadrature_points() const;
+
+ /**
+ * The number of values in each
+ * quadrature point in
+ * #quadrature_values.
+ */
+ unsigned int n_quadrature_values() const;
+
+ /**
+ * Access scalar value at index
+ * @p i.
+ */
+ number& value(unsigned int i);
+
+ /**
+ * Read scalar value at index
+ * @p i.
+ */
+ number value(unsigned int i) const;
+
+ /**
+ * Access vector at index @p i.
+ */
+ BlockVector<number>& vector(unsigned int i);
+
+ /**
+ * Read vector at index @p i.
+ */
+ const BlockVector<number>& vector(unsigned int i) const;
+
+ /**
+ * Access matrix at index @p
+ * i. For results on internal
+ * faces, a true value for @p
+ * external refers to the flux
+ * between cells, while false
+ * refers to entries coupling
+ * inside the cell.
+ */
+ MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false);
+
+ /**
+ * Read matrix at index @p
+ * i. For results on internal
+ * faces, a true value for @p
+ * external refers to the flux
+ * between cells, while false
+ * refers to entries coupling
+ * inside the cell.
+ */
+ const MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false) const;
+
+ /**
+ * Access the <i>i</i>th value
+ * at quadrature point <i>k</i>
+ */
+ number& quadrature_value(unsigned int k, unsigned int i);
+
+ /**
+ * Read the <i>i</i>th value
+ * at quadrature point <i>k</i>
+ */
+ number quadrature_value(unsigned int k, unsigned int i) const;
+
+ private:
+ /**
+ * The local numbers,
+ * computed on a cell or on a
+ * face.
+ */
+ std::vector<number> J;
+
+ /**
+ * The local vectors. This
+ * field is public, so that
+ * local integrators can
+ * write to it.
+ */
+ std::vector<BlockVector<number> > R;
+
+ /**
+ * The local matrices
+ * coupling degrees of
+ * freedom in the cell
+ * itself or within the
+ * first cell on a face.
+ */
+ std::vector<MatrixBlock<FullMatrix<number> > > M1;
+
+ /**
+ * The local matrices
+ * coupling test functions on
+ * the cell with trial
+ * functions on the other
+ * cell.
+ *
+ * Only used on interior
+ * faces.
+ */
+ std::vector<MatrixBlock<FullMatrix<number> > > M2;
+ /**
+ * Values in quadrature points.
+ */
+ std::vector<std::vector<number> > quadrature_values;
+ };
+
+//----------------------------------------------------------------------//
+
+ template <typename number>
+ inline void
+ LocalResults<number>::initialize_numbers(unsigned int n)
+ {
+ J.resize(n);
+ }
+
+
+ template <typename number>
+ inline void
+ LocalResults<number>::initialize_vectors(const unsigned int n)
+ {
+ R.resize(n);
+ }
+
+
+ template <typename number>
+ template <class MATRIX>
+ inline void
+ LocalResults<number>::initialize_matrices(
+ const MatrixBlockVector<MATRIX>& matrices,
+ bool both)
+ {
+ M1.resize(matrices.size());
+ if (both)
+ M2.resize(matrices.size());
+ for (unsigned int i=0;i<matrices.size();++i)
+ {
+ const unsigned int row = matrices.block(i).row;
+ const unsigned int col = matrices.block(i).column;
+
+ M1[i].row = row;
+ M1[i].column = col;
+ if (both)
+ {
+ M2[i].row = row;
+ M2[i].column = col;
+ }
+ }
+ }
+
+
+ template <typename number>
+ inline void
+ LocalResults<number>::initialize_matrices(const unsigned int n,
+ const bool both)
+ {
+ M1.resize(n);
+ if (both)
+ M2.resize(n);
+ for (unsigned int i=0;i<n;++i)
+ {
+ M1[i].row = 0;
+ M1[i].column = 0;
+ if (both)
+ {
+ M2[i].row = 0;
+ M2[i].column = 0;
+ }
+ }
+ }
+
+
+ template <typename number>
+ inline void
+ LocalResults<number>::initialize_quadrature(unsigned int np, unsigned int nv)
+ {
+ quadrature_values.resize(np, std::vector<number>(nv));
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_values() const
+ {
+ return J.size();
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_vectors() const
+ {
+ return R.size();
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_matrices() const
+ {
+ return M1.size();
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_quadrature_points() const
+ {
+ return quadrature_values.size();
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_quadrature_values() const
+ {
+ Assert(quadrature_values.size() != 0, ExcNotInitialized());
+ return quadrature_values[0].size();
+ }
+
+
+ template <typename number>
+ inline
+ number&
+ LocalResults<number>::value(unsigned int i)
+ {
+ AssertIndexRange(i,J.size());
+ return J[i];
+ }
+
+
+ template <typename number>
+ inline
+ BlockVector<number>&
+ LocalResults<number>::vector(unsigned int i)
+ {
+ AssertIndexRange(i,R.size());
+ return R[i];
+ }
+
+
+ template <typename number>
+ inline
+ MatrixBlock<FullMatrix<number> >&
+ LocalResults<number>::matrix(unsigned int i, bool external)
+ {
+ if (external)
+ {
+ AssertIndexRange(i,M2.size());
+ return M2[i];
+ }
+ AssertIndexRange(i,M1.size());
+ return M1[i];
+ }
+
+
+ template <typename number>
+ inline
+ number&
+ LocalResults<number>::quadrature_value(unsigned int k, unsigned int i)
+ {
+ AssertIndexRange(k,quadrature_values.size());
+ AssertIndexRange(i,quadrature_values[0].size());
+ return quadrature_values[k][i];
+ }
+
+
+ template <typename number>
+ inline
+ number
+ LocalResults<number>::value(unsigned int i) const
+ {
+ AssertIndexRange(i,J.size());
+ return J[i];
+ }
+
+
+ template <typename number>
+ inline
+ const BlockVector<number>&
+ LocalResults<number>::vector(unsigned int i) const
+ {
+ AssertIndexRange(i,R.size());
+ return R[i];
+ }
+
+
+ template <typename number>
+ inline
+ const MatrixBlock<FullMatrix<number> >&
+ LocalResults<number>::matrix(unsigned int i, bool external) const
+ {
+ if (external)
+ {
+ AssertIndexRange(i,M2.size());
+ return M2[i];
+ }
+ AssertIndexRange(i,M1.size());
+ return M1[i];
+ }
+
+
+ template <typename number>
+ inline
+ number
+ LocalResults<number>::quadrature_value(unsigned int k, unsigned int i) const
+ {
+ AssertIndexRange(k,quadrature_values.size());
+ AssertIndexRange(i,quadrature_values[0].size());
+ return quadrature_values[k][i];
+ }
+}
DEAL_II_NAMESPACE_CLOSE
#define __deal2__mesh_worker_info_h
#include <base/config.h>
-#include <lac/matrix_block.h>
#include <boost/shared_ptr.hpp>
#include <dofs/block_info.h>
#include <fe/fe_values.h>
+#include <numerics/mesh_worker.h>
#include <numerics/mesh_worker_vector_selector.h>
DEAL_II_NAMESPACE_OPEN
namespace MeshWorker
{
template <int dim, class DOFINFO> class DoFInfoBox;
-
-/**
- * The class providing the scrapbook to fill with local integration
- * results. These can be values, local contributions to forms or cell
- * and face matrices.
- *
- * The local matrices initialized by reinit() of the info object and
- * then assembled into the global system by Assembler classes.
- *
- * @ingroup MeshWorker
- * @author Guido Kanschat, 2009
- */
- template <typename number>
- class LocalResults
- {
- private:
- /**
- * Initialize a single local
- * matrix block. A helper
- * function for initialize()
- */
- void initialize_local(MatrixBlock<FullMatrix<number> >& M,
- const unsigned int row,
- const unsigned int col);
-
- public:
- void initialize_numbers(const unsigned int n);
- void initialize_vectors(const unsigned int n);
- /**
- * Allocate @p n local
- * matrices. Additionally,
- * set their block row and
- * column coordinates to
- * zero. The matrices
- * themselves are resized by
- * reinit().
- *
- * The template parameter @p
- * MatrixPtr should point to
- * a MatrixBlock
- * instantiation in order to
- * provide row and column info.
- */
- void initialize_matrices(unsigned int n, bool both);
-
- /**
- * Allocate a local matrix
- * for each of the global
- * ones in @p
- * matrices. Additionally,
- * set their block row and
- * column coordinates. The
- * matrices themselves are
- * resized by reinit().
- *
- * The template parameter @p
- * MatrixPtr should point to
- * a MatrixBlock
- * instantiation in order to
- * provide row and column info.
- */
- template <class MATRIX>
- void initialize_matrices(const MatrixBlockVector<MATRIX>& matrices,
- bool both);
-
- /**
- * Initialize quadrature values
- * to <tt>nv</tt> values in
- * <tt>np</tt> quadrature points.
- */
- void initialize_quadrature(unsigned int np, unsigned int nv);
-
- /**
- * Reinitialize matrices for
- * new cell. Resizes the
- * matrices for hp and sets
- * them to zero.
- */
- void reinit(const BlockIndices& local_sizes);
-
- /**
- * The number of scalar values.
- */
- unsigned int n_values () const;
-
- /**
- * The number of vectors.
- */
- unsigned int n_vectors () const;
-
- /**
- * The number of matrices.
- */
- unsigned int n_matrices () const;
-
- /**
- * The number of quadrature
- * points in #quadrature_values.
- */
- unsigned int n_quadrature_points() const;
-
- /**
- * The number of values in each
- * quadrature point in
- * #quadrature_values.
- */
- unsigned int n_quadrature_values() const;
-
- /**
- * Access scalar value at index
- * @p i.
- */
- number& value(unsigned int i);
-
- /**
- * Read scalar value at index
- * @p i.
- */
- number value(unsigned int i) const;
-
- /**
- * Access vector at index @p i.
- */
- BlockVector<number>& vector(unsigned int i);
-
- /**
- * Read vector at index @p i.
- */
- const BlockVector<number>& vector(unsigned int i) const;
-
- /**
- * Access matrix at index @p
- * i. For results on internal
- * faces, a true value for @p
- * external refers to the flux
- * between cells, while false
- * refers to entries coupling
- * inside the cell.
- */
- MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false);
-
- /**
- * Read matrix at index @p
- * i. For results on internal
- * faces, a true value for @p
- * external refers to the flux
- * between cells, while false
- * refers to entries coupling
- * inside the cell.
- */
- const MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false) const;
-
- /**
- * Access the <i>i</i>th value
- * at quadrature point <i>k</i>
- */
- number& quadrature_value(unsigned int k, unsigned int i);
-
- /**
- * Read the <i>i</i>th value
- * at quadrature point <i>k</i>
- */
- number quadrature_value(unsigned int k, unsigned int i) const;
-
- private:
- /**
- * The local numbers,
- * computed on a cell or on a
- * face.
- */
- std::vector<number> J;
-
- /**
- * The local vectors. This
- * field is public, so that
- * local integrators can
- * write to it.
- */
- std::vector<BlockVector<number> > R;
-
- /**
- * The local matrices
- * coupling degrees of
- * freedom in the cell
- * itself or within the
- * first cell on a face.
- */
- std::vector<MatrixBlock<FullMatrix<number> > > M1;
-
- /**
- * The local matrices
- * coupling test functions on
- * the cell with trial
- * functions on the other
- * cell.
- *
- * Only used on interior
- * faces.
- */
- std::vector<MatrixBlock<FullMatrix<number> > > M2;
-
- /**
- * Values in quadrature points.
- */
- std::vector<std::vector<number> > quadrature_values;
- };
-
- template <int dim, class DOFINFO> class DoFInfoBox;
-
/**
* A class containing information on geometry and degrees of freedom
* of a mesh object.
};
-//----------------------------------------------------------------------//
-
- template <typename number>
- inline void
- LocalResults<number>::initialize_numbers(unsigned int n)
- {
- J.resize(n);
- }
-
-
- template <typename number>
- inline void
- LocalResults<number>::initialize_vectors(const unsigned int n)
- {
- R.resize(n);
- }
-
-
- template <typename number>
- template <class MATRIX>
- inline void
- LocalResults<number>::initialize_matrices(
- const MatrixBlockVector<MATRIX>& matrices,
- bool both)
- {
- M1.resize(matrices.size());
- if (both)
- M2.resize(matrices.size());
- for (unsigned int i=0;i<matrices.size();++i)
- {
- const unsigned int row = matrices.block(i).row;
- const unsigned int col = matrices.block(i).column;
-
- M1[i].row = row;
- M1[i].column = col;
- if (both)
- {
- M2[i].row = row;
- M2[i].column = col;
- }
- }
- }
-
-
- template <typename number>
- inline void
- LocalResults<number>::initialize_matrices(const unsigned int n,
- const bool both)
- {
- M1.resize(n);
- if (both)
- M2.resize(n);
- for (unsigned int i=0;i<n;++i)
- {
- M1[i].row = 0;
- M1[i].column = 0;
- if (both)
- {
- M2[i].row = 0;
- M2[i].column = 0;
- }
- }
- }
-
-
- template <typename number>
- inline void
- LocalResults<number>::initialize_quadrature(unsigned int np, unsigned int nv)
- {
- quadrature_values.resize(np, std::vector<number>(nv));
- }
-
-
- template <typename number>
- inline void
- LocalResults<number>::reinit(const BlockIndices& bi)
- {
- for (unsigned int i=0;i<J.size();++i)
- J[i] = 0.;
- for (unsigned int i=0;i<R.size();++i)
- R[i].reinit(bi);
- for (unsigned int i=0;i<M1.size();++i)
- M1[i].matrix.reinit(bi.block_size(M1[i].row),
- bi.block_size(M1[i].column));
- for (unsigned int i=0;i<M2.size();++i)
- M2[i].matrix.reinit(bi.block_size(M2[i].row),
- bi.block_size(M2[i].column));
- for (unsigned int k=0;k<quadrature_values.size();++k)
- for (unsigned int i=0;i<quadrature_values[k].size();++i)
- quadrature_values[k][i] = 0.;
- }
-
-
- template <typename number>
- inline
- unsigned int
- LocalResults<number>::n_values() const
- {
- return J.size();
- }
-
-
- template <typename number>
- inline
- unsigned int
- LocalResults<number>::n_vectors() const
- {
- return R.size();
- }
-
-
- template <typename number>
- inline
- unsigned int
- LocalResults<number>::n_matrices() const
- {
- return M1.size();
- }
-
-
- template <typename number>
- inline
- unsigned int
- LocalResults<number>::n_quadrature_points() const
- {
- return quadrature_values.size();
- }
-
-
- template <typename number>
- inline
- unsigned int
- LocalResults<number>::n_quadrature_values() const
- {
- Assert(quadrature_values.size() != 0, ExcNotInitialized());
- return quadrature_values[0].size();
- }
-
-
- template <typename number>
- inline
- number&
- LocalResults<number>::value(unsigned int i)
- {
- AssertIndexRange(i,J.size());
- return J[i];
- }
-
-
- template <typename number>
- inline
- BlockVector<number>&
- LocalResults<number>::vector(unsigned int i)
- {
- AssertIndexRange(i,R.size());
- return R[i];
- }
-
-
- template <typename number>
- inline
- MatrixBlock<FullMatrix<number> >&
- LocalResults<number>::matrix(unsigned int i, bool external)
- {
- if (external)
- {
- AssertIndexRange(i,M2.size());
- return M2[i];
- }
- AssertIndexRange(i,M1.size());
- return M1[i];
- }
-
-
- template <typename number>
- inline
- number&
- LocalResults<number>::quadrature_value(unsigned int k, unsigned int i)
- {
- AssertIndexRange(k,quadrature_values.size());
- AssertIndexRange(i,quadrature_values[0].size());
- return quadrature_values[k][i];
- }
-
-
- template <typename number>
- inline
- number
- LocalResults<number>::value(unsigned int i) const
- {
- AssertIndexRange(i,J.size());
- return J[i];
- }
-
-
- template <typename number>
- inline
- const BlockVector<number>&
- LocalResults<number>::vector(unsigned int i) const
- {
- AssertIndexRange(i,R.size());
- return R[i];
- }
-
-
- template <typename number>
- inline
- const MatrixBlock<FullMatrix<number> >&
- LocalResults<number>::matrix(unsigned int i, bool external) const
- {
- if (external)
- {
- AssertIndexRange(i,M2.size());
- return M2[i];
- }
- AssertIndexRange(i,M1.size());
- return M1[i];
- }
-
-
- template <typename number>
- inline
- number
- LocalResults<number>::quadrature_value(unsigned int k, unsigned int i) const
- {
- AssertIndexRange(k,quadrature_values.size());
- AssertIndexRange(i,quadrature_values[0].size());
- return quadrature_values[k][i];
- }
-
-
//----------------------------------------------------------------------//
template <int dim, int spacedim, typename number>
p->initialize(data);
cell_data = p;
cell.initialize_data(p);
-
+
p = boost::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (boundary_selector));
p->initialize(data);
boundary_data = p;
for (unsigned int i=0;i<this->n_gradients();++i)
{
- const VECTOR& src = *data(this->value_index(i));
+ const VECTOR& src = *data(this->gradient_index(i));
VectorSlice<std::vector<std::vector<Tensor<1,dim> > > > dst(gradients[i], component, n_comp);
fe.get_function_gradients(src, make_slice(index, start, size), dst, true);
}
for (unsigned int i=0;i<this->n_hessians();++i)
{
- const VECTOR& src = *data(this->value_index(i));
+ const VECTOR& src = *data(this->hessian_index(i));
VectorSlice<std::vector<std::vector<Tensor<2,dim> > > > dst(hessians[i], component, n_comp);
fe.get_function_hessians(src, make_slice(index, start, size), dst, true);
}
//
//---------------------------------------------------------------------------
-#include <multigrid/mg_tools.h>
-#include <fe/fe.h>
-#include <fe/fe_tools.h>
#include <numerics/mesh_worker.h>
-#include <base/quadrature_lib.h>
+#include <lac/block_indices.h>
DEAL_II_NAMESPACE_OPEN
+using namespace MeshWorker;
+
+template <typename number>
+void
+LocalResults<number>::reinit(const BlockIndices& bi)
+{
+ for (unsigned int i=0;i<J.size();++i)
+ J[i] = 0.;
+ for (unsigned int i=0;i<R.size();++i)
+ R[i].reinit(bi);
+ for (unsigned int i=0;i<M1.size();++i)
+ M1[i].matrix.reinit(bi.block_size(M1[i].row),
+ bi.block_size(M1[i].column));
+ for (unsigned int i=0;i<M2.size();++i)
+ M2[i].matrix.reinit(bi.block_size(M2[i].row),
+ bi.block_size(M2[i].column));
+ for (unsigned int k=0;k<quadrature_values.size();++k)
+ for (unsigned int i=0;i<quadrature_values[k].size();++i)
+ quadrature_values[k][i] = 0.;
+}
+
+template class LocalResults<float>;
+template class LocalResults<double>;
+
DEAL_II_NAMESPACE_CLOSE
// for using the MeshWorker framework:
#include <numerics/mesh_worker.h>
#include <numerics/mesh_worker_info.h>
+#include <numerics/mesh_worker_assembler.h>
#include <numerics/mesh_worker_loop.h>
// Like in all programs, we finish
// The include files for using the
// MeshWorker framework
#include <numerics/mesh_worker.h>
+#include <numerics/mesh_worker_info.h>
+#include <numerics/mesh_worker_assembler.h>
#include <numerics/mesh_worker_loop.h>
// Support for multigrid methods
// solution we just computed.
NamedData<Vector<double>* > solution_data;
solution_data.add(&solution, "solution");
+
// Then, we tell the Meshworker::VectorSelector
// for cells, that we need the
// second derivatives of this
projects / dealii-svn.git / commitdiff