* communicate the current time step information through AnyData as
* well. Therefore, the AnyData objects handed as input to
* #op_explicit and #op_implicit contain two entries of type
- * <tt>const double</tt> reference named "Time" and "Timestep". Note
+ * <tt>const double*</tt> named "Time" and "Timestep". Note
* that "Time" refers to the time at the beginning of the current
* step for #op_explicit and at the end for #op_implicit,
* respectively.
// timestep
AnyData src1;
src1.add<const VECTOR*>(&solution, "Previous iterate");
- src1.add<const double&>(d_explicit.time, "Time");
- src1.add<const double&>(d_explicit.step, "Timestep");
+ src1.add<const double*>(&d_explicit.time, "Time");
+ src1.add<const double*>(&d_explicit.step, "Timestep");
src1.merge(in);
AnyData src2;
out1.add<VECTOR*>(aux, "Solution");
// The data provided to the inner solver
src2.add<const VECTOR*>(aux, "Previous time");
- src2.add<const double&>(d_implicit.time, "Time");
- src2.add<const double&>(d_explicit.step, "Timestep");
+ src2.add<const double*>(&d_implicit.time, "Time");
+ src2.add<const double*>(&d_explicit.step, "Timestep");
src2.merge(in);
if (output != 0)
#define __deal2__mesh_worker_functional_h
#include <deal.II/base/named_data.h>
+#include <deal.II/algorithms/any_data.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/lac/block_vector.h>
{
public:
/**
- * The data type for
- * communicating the cell and
- * face vectors.
- */
- typedef NamedData<BlockVector<number>*> DataVectors;
-
- /**
- * The initialization
- * function, specifying the
- * @p results vectors and
- * whether face data should
- * be collected separately.
+ * The initialization function, specifying the @p results
+ * vectors and whether face data should be collected separately.
*
- * @p results should contain
- * two block vectors named
- * "cells" and "faces" (the
- * latter only if
- * @p separate_faces is
- * true). In each of the two,
- * each block should have
- * equal size and be large
- * enough to accommodate all
- * user indices set in the
- * cells and faces covered by
- * the loop it is used
- * in. Typically, for
- * estimators, this is
- * Triangulation::n_active_cells()
- * and
- * Triangulation::n_faces(),
- * respectively.
+ * @p results should contain two block vectors named "cells" and
+ * "faces" (the latter only if @p separate_faces is true). In
+ * each of the two, each block should have equal size and be
+ * large enough to accommodate all user indices set in the cells
+ * and faces covered by the loop it is used in. Typically, for
+ * estimators, this is Triangulation::n_active_cells() and
+ * Triangulation::n_faces(), respectively.
*
- * The use of BlockVector may
- * seem cumbersome, but it
- * allows us to assemble
- * several functionals at the
- * same time, one in each
- * block. The typical
- * situation for error
- * estimate is just having a
- * single block in each vector.
+ * The use of BlockVector may seem cumbersome, but it allows us
+ * to assemble several functionals at the same time, one in each
+ * block. The typical situation for error estimate is just
+ * having a single block in each vector.
*/
- void initialize(DataVectors &results,
+ void initialize(AnyData &results, bool separate_faces = true);
+
+ /**
+ * @deprecated
+ */
+ void initialize(NamedData<BlockVector<number>*> &results,
bool separate_faces = true);
/**
* Initialize the local data
*/
number operator() (const unsigned int i) const;
private:
- DataVectors results;
+ AnyData results;
bool separate_faces;
};
//----------------------------------------------------------------------//
template <typename number>
inline void
- CellsAndFaces<number>::initialize(DataVectors &r, bool sep)
+ CellsAndFaces<number>::initialize(AnyData &r, bool sep)
+ {
+ Assert(r.name(0) == "cells", AnyData::ExcNameMismatch(0, "cells"));
+ if (sep)
+ {
+ Assert(r.name(1) == "faces", AnyData::ExcNameMismatch(1, "faces"));
+ AssertDimension(r.entry<BlockVector<double>*>(0)->n_blocks(),
+ r.entry<BlockVector<double>*>(1)->n_blocks());
+ }
+
+ results = r;
+ separate_faces = sep;
+ }
+
+ template <typename number>
+ inline void
+ CellsAndFaces<number>::initialize(NamedData<BlockVector<number>*> &r, bool sep)
{
- Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells"));
+ Assert(r.name(0) == "cells", AnyData::ExcNameMismatch(0, "cells"));
if (sep)
{
- Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces"));
- AssertDimension(r(0)->n_blocks(), r(1)->n_blocks());
+ Assert(r.name(1) == "faces", AnyData::ExcNameMismatch(1, "faces"));
+ AssertDimension(r(0)->n_blocks(),r(1)->n_blocks());
}
results = r;
inline void
CellsAndFaces<number>::initialize_info(DOFINFO &info, bool) const
{
- info.initialize_numbers(results(0)->n_blocks());
+ info.initialize_numbers(results.entry<BlockVector<double>*>(0)->n_blocks());
}
inline void
CellsAndFaces<number>::assemble(const DOFINFO &info)
{
+ BlockVector<double>* v;
+ if (separate_faces &&
+ info.face_number != deal_II_numbers::invalid_unsigned_int)
+ v = results.entry<BlockVector<double>*>(1);
+ else
+ v = results.entry<BlockVector<double>*>(0);
+
for (unsigned int i=0; i<info.n_values(); ++i)
- {
- if (separate_faces &&
- info.face_number != deal_II_numbers::invalid_unsigned_int)
- results(1)->block(i)(info.face->user_index()) += info.value(i);
- else
- results(0)->block(i)(info.cell->user_index()) += info.value(i);
- }
+ v->block(i)(info.cell->user_index()) += info.value(i);
}
{
if (separate_faces)
{
+ BlockVector<double>* v1 = results.entry<BlockVector<double>*>(1);
const double J = info1.value(i) + info2.value(i);
- results(1)->block(i)(info1.face->user_index()) += J;
+ v1->block(i)(info1.face->user_index()) += J;
if (info2.face != info1.face)
- results(1)->block(i)(info2.face->user_index()) += J;
+ v1->block(i)(info2.face->user_index()) += J;
}
else
{
- results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i);
- results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i);
+ BlockVector<double>* v0 = results.entry<BlockVector<double>*>(0);
+ v0->block(i)(info1.cell->user_index()) += .5*info1.value(i);
+ v0->block(i)(info2.cell->user_index()) += .5*info2.value(i);
}
}
}
* undertaken to use this class.
*
* First, you should consider if you need values from any vectors in a
- * NamedData object. If so, fill the VectorSelector objects
+ * AnyData object. If so, fill the VectorSelector objects
* #cell_selector, #boundary_selector and #face_selector with their names
* and the data type (value, gradient, Hessian) to be extracted.
*
* formulas, which integrate polynomial bilinear forms exactly.
*/
template <typename VECTOR>
+ void initialize(const FiniteElement<dim, spacedim> &el,
+ const Mapping<dim, spacedim> &mapping,
+ const AnyData &data,
+ const VECTOR& dummy,
+ const BlockInfo *block_info = 0);
+ /**
+ * @deprecated Use AnyData instead of NamedData.
+ */
+ template <typename VECTOR>
void initialize(const FiniteElement<dim, spacedim> &el,
const Mapping<dim, spacedim> &mapping,
const NamedData<VECTOR *> &data,
const BlockInfo *block_info = 0);
/**
- * Initialize the IntegrationInfo objects contained.
- *
- * Before doing so, add update flags necessary to produce the data
- * needed and also set uninitialized quadrature rules to Gauss
- * formulas, which integrate polynomial bilinear forms exactly.
+ * @deprecated Use AnyData instead of NamedData.
*/
template <typename VECTOR>
void initialize(const FiniteElement<dim, spacedim> &el,
}
+ template <int dim, int sdim>
+ template <typename VECTOR>
+ void
+ IntegrationInfoBox<dim,sdim>::initialize(
+ const FiniteElement<dim,sdim> &el,
+ const Mapping<dim,sdim> &mapping,
+ const AnyData &data,
+ const VECTOR& dummy,
+ const BlockInfo *block_info)
+ {
+ initialize(el, mapping, block_info);
+ std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> > p;
+ VectorDataBase<dim,sdim>* pp;
+
+ p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (cell_selector));
+ // Public member function of parent class was not found without
+ // explicit cast
+ pp = &*p;
+ pp->initialize(data);
+ cell_data = p;
+ cell.initialize_data(p);
+
+ p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (boundary_selector));
+ pp = &*p;
+ pp->initialize(data);
+ boundary_data = p;
+ boundary.initialize_data(p);
+
+ p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (face_selector));
+ pp = &*p;
+ pp->initialize(data);
+ face_data = p;
+ face.initialize_data(p);
+ subface.initialize_data(p);
+ neighbor.initialize_data(p);
+ }
+
+
template <int dim, int sdim>
template <typename VECTOR>
void
#define __deal2__mesh_worker_simple_h
#include <deal.II/base/named_data.h>
+#include <deal.II/algorithms/any_data.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/lac/block_vector.h>
class ResidualSimple
{
public:
- void initialize(NamedData<VECTOR *> &results);
+ /**
+ * Initialize with an AnyData object holding the result of
+ * assembling.
+ *
+ * Assembling currently writes into the first vector of <tt>results</tt>.
+ */
+ void initialize(AnyData& results);
+ /**
+ * @deprecated Use initialize(AnyData&) instead.
+ */
+ void initialize(NamedData<VECTOR *> &results);
/**
* Initialize the constraints.
*/
/**
* The global residal vectors filled by assemble().
*/
- NamedData<SmartPointer<VECTOR,ResidualSimple<VECTOR> > > residuals;
+ AnyData residuals;
/**
* A pointer to the object containing constraints.
*/
//----------------------------------------------------------------------//
+ template <class VECTOR>
+ inline void
+ ResidualSimple<VECTOR>::initialize(AnyData &results)
+ {
+ residuals = results;
+ }
+
template <class VECTOR>
inline void
ResidualSimple<VECTOR>::initialize(NamedData<VECTOR *> &results)
{
for (unsigned int k=0; k<residuals.size(); ++k)
{
+ VECTOR* v = residuals.entry<VECTOR*>(k);
if (constraints == 0)
{
for (unsigned int i=0; i<info.vector(k).block(0).size(); ++i)
- (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
+ (*v)(info.indices[i]) += info.vector(k).block(0)(i);
}
else
{
if (info.indices_by_block.size() == 0)
- constraints->distribute_local_to_global(info.vector(k).block(0), info.indices, (*residuals(k)));
+ constraints->distribute_local_to_global(info.vector(k).block(0), info.indices, *v);
else
for (unsigned int i=0; i != info.vector(k).n_blocks(); ++i)
- constraints->distribute_local_to_global(info.vector(k).block(i), info.indices_by_block[i], (*residuals(k)));
+ constraints->distribute_local_to_global(info.vector(k).block(i), info.indices_by_block[i], *v);
}
}
}
{
for (unsigned int k=0; k<residuals.size(); ++k)
{
+ VECTOR* v = residuals.entry<VECTOR*>(k);
if (constraints == 0)
{
for (unsigned int i=0; i<info1.vector(k).block(0).size(); ++i)
- (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
+ (*v)(info1.indices[i]) += info1.vector(k).block(0)(i);
for (unsigned int i=0; i<info2.vector(k).block(0).size(); ++i)
- (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
+ (*v)(info2.indices[i]) += info2.vector(k).block(0)(i);
}
else
{
if (info1.indices_by_block.size() == 0 && info2.indices_by_block.size() == 0)
{
constraints->distribute_local_to_global
- (info1.vector(k).block(0), info1.indices, (*residuals(k)));
+ (info1.vector(k).block(0), info1.indices, *v);
constraints->distribute_local_to_global
- (info2.vector(k).block(0), info2.indices, (*residuals(k)));
+ (info2.vector(k).block(0), info2.indices, *v);
}
else if (info1.indices_by_block.size() != 0 && info2.indices_by_block.size() != 0)
{
for (unsigned int i=0; i<info1.vector(k).n_blocks(); ++i)
{
constraints->distribute_local_to_global
- (info1.vector(k).block(i), info1.indices_by_block[i], (*residuals(k)));
+ (info1.vector(k).block(i), info1.indices_by_block[i], *v);
constraints->distribute_local_to_global
- (info2.vector(k).block(i), info2.indices_by_block[i], (*residuals(k)));
+ (info2.vector(k).block(i), info2.indices_by_block[i], *v);
}
}
else
inline void
SystemSimple<MATRIX,VECTOR>::initialize(MATRIX &m, VECTOR &rhs)
{
- NamedData<VECTOR *> data;
+ AnyData data;
VECTOR *p = &rhs;
data.add(p, "right hand side");
#define __deal2__mesh_worker_vector_selector_h
#include <deal.II/base/named_data.h>
+#include <deal.II/algorithms/any_data.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/mg_level_object.h>
/**
* A class that selects vectors from a list of named vectors.
*
- * Since the number of vectors in a NamedData object may grow with every
+ * Since the number of vectors in an AnyData object may grow with every
* nesting of applications or loops, it is important to be able to
* select those, which are actually used in computing residuals etc.
* This class organizes the selection.
{
public:
/**
- * Add a vector to the
- * selection of finite element
- * functions. The arguments are
- * the name of the vector and
- * indicators, which
- * information is to be
- * extracted from the
- * vector. The name refers to
- * an entry in a NamedData
- * object, which will be
- * identified by initialize().
- * The three bool parameters
- * indicate, whether values,
- * greadients and Hessians of
- * the finite element function
- * are to be computed on each
- * cell or face.
+ * Add a vector to the selection of finite element functions. The
+ * arguments are the name of the vector and indicators, which
+ * information is to be extracted from the vector. The name refers
+ * to an entry in a AnyData object, which will be identified by
+ * initialize(). The three bool parameters indicate, whether
+ * values, greadients and Hessians of the finite element function
+ * are to be computed on each cell or face.
*/
void add(const std::string &name,
const bool values = true,
const bool hessians = false);
/**
- * Does the same as the function
- * above but it is possible to
- * select a block of the global
- * vector.
+ * Does the same as the function above but it is possible to
+ * select a block of the global vector.
*/
// void add(const std::string& name,
// const unsigned int selected_block,
// bool hessians = false);
/**
- * Initialize the selection
- * field with a data
- * vector. While add() only
- * enters the names of vectors,
- * which will be used in the integration
- * loop over cells and faces,
- * this function links the
- * names to actual vectos in a
- * NamedData object.
+ * Initialize the selection field with a data vector. While add()
+ * only enters the names of vectors, which will be used in the
+ * integration loop over cells and faces, this function links the
+ * names to actual vectos in a AnyData object.
*
- * @note This function caches
- * the index associated with a
- * name. Therefore, it must be
- * called every time after the
- * NamedData object has changed.
+ * @note This function caches the index associated with a
+ * name. Therefore, it must be called every time after the AnyData
+ * object has changed.
+ */
+ void initialize(const AnyData &);
+
+ /**
+ * @deprecated Use AnyData instead of NamedData.
*/
template <class DATA>
void initialize(const NamedData<DATA> &);
bool empty () const;
/**
- * Returns true if values are
- * selected for any vector.
+ * Returns true if values are selected for any vector.
*/
bool has_values () const;
/**
- * Returns true if gradients are
- * selected for any vector.
+ * Returns true if gradients are selected for any vector.
*/
bool has_gradients () const;
/**
- * Returns true if hessians are
- * selected for any vector.
+ * Returns true if hessians are selected for any vector.
*/
bool has_hessians () const;
* selection to the stream.
*/
template <class STREAM, typename DATA>
+ void print (STREAM &s, const AnyData &v) const;
+
+ /**
+ * @deprecated Use AnyData instead of NamedData.
+ */
+ template <class STREAM, typename DATA>
void print (STREAM &s, const NamedData<DATA> &v) const;
/**
NamedSelection value_selection;
/**
- * Selection of the vectors
- * used to compute gradients.
+ * Selection of the vectors used to compute gradients.
*/
NamedSelection gradient_selection;
/**
- * Selection of the vectors
- * used to compute hessians.
+ * Selection of the vectors used to compute hessians.
*/
NamedSelection hessian_selection;
};
VectorDataBase();
/**
- * Constructor from a base
- * class object
+ * Constructor from a base class object
*/
VectorDataBase(const VectorSelector &);
/**
- * Virtual, but empty
- * destructor.
+ * Initialize with a AnyData object and cache the indices in the
+ * VectorSelector base class.
+ *
+ * @note Make sure the VectorSelector base class was filled with
+ * reasonable data before calling this function.
+ */
+ void initialize(const AnyData&);
+
+ /**
+ * Virtual, but empty destructor.
*/
virtual ~VectorDataBase();
/**
- * The only function added to
- * VectorSelector is an
- * abstract virtual function
- * implemented in the derived
- * class template and called by
- * IntegrationInfo.
+ * The only function added to VectorSelector is an abstract
+ * virtual function implemented in the derived class template and
+ * called by IntegrationInfo.
*
- * Depending on the selections
- * made in our base class, this
- * fills the first three
- * arguments with the local
- * data of the finite element
- * functions. It is usually
- * called either for the whole
- * FESystem, or for each base
- * element separately.
+ * Depending on the selections made in our base class, this fills
+ * the first three arguments with the local data of the finite
+ * element functions. It is usually called either for the whole
+ * FESystem, or for each base element separately.
*
- * @param values is the vector
- * filled with the values of
- * the finite element function
- * in the quadrature points.
+ * @param values is the vector filled with the values of the
+ * finite element function in the quadrature points.
*
- * @param gradients is the vector
- * filled with the derivatives of
- * the finite element function
- * in the quadrature points.
+ * @param gradients is the vector filled with the derivatives of
+ * the finite element function in the quadrature points.
*
- * @param hessians is the
- * vector filled with the
- * second derivatives of the
- * finite element function in
- * the quadrature points.
+ * @param hessians is the vector filled with the second
+ * derivatives of the finite element function in the quadrature
+ * points.
*
- * @param fe is the
- * FEValuesBase object which is
- * used to compute the function
- * values. Its UpdateFlags must
- * have been set appropriately.
+ * @param fe is the FEValuesBase object which is used to compute
+ * the function values. Its UpdateFlags must have been set
+ * appropriately.
*
- * @param index is the local
- * index vector. If @p fe
- * refers to base elements of
- * the system, this vector
- * should be sorted by block
- * and the arguments @p start
- * and @p size below specify
- * the subset of @p indices
- * used.
+ * @param index is the local index vector. If @p fe refers to base
+ * elements of the system, this vector should be sorted by block
+ * and the arguments @p start and @p size below specify the subset
+ * of @p indices used.
*
- * @param component is the first
- * index in @p values, @p
- * gradients and @p hessians
- * entered in this function.
+ * @param component is the first index in @p values, @p gradients
+ * and @p hessians entered in this function.
*
- * @param n_comp is the number
- * of components to be filled.
+ * @param n_comp is the number of components to be filled.
*
- * @param start is the first
- * index of this block in @p
- * indices, or zero if
- * no base elements are used.
+ * @param start is the first index of this block in @p indices, or
+ * zero if no base elements are used.
*
- * @param size is the number of
- * dofs per cell of the current
+ * @param size is the number of dofs per cell of the current
* element or base element.
*/
virtual void fill(
const unsigned int size) const;
/**
- * Fill the local data vector
- * from level vectors. Performs
- * exactly what the other
- * fill() does, but uses the
- * cell level to access a
- * single level out of a
- * hierarchy of level vectors,
- * instead of a global data
- * vector on the active cells.
+ * Fill the local data vector from level vectors. Performs exactly
+ * what the other fill() does, but uses the cell level to access a
+ * single level out of a hierarchy of level vectors, instead of a
+ * global data vector on the active cells.
*/
virtual void mg_fill(
std::vector<std::vector<std::vector<double> > > &values,
* The memory used by this object.
*/
std::size_t memory_consumption () const;
+ protected:
+ AnyData data;
};
/**
* Based on VectorSelector, this is the class that implements the
* function VectorDataBase::fill() for a certain type of vector, using
- * NamedData to identify vectors by name.
+ * AnyData to identify vectors by name.
*
* @ingroup MeshWorker
* @author Guido Kanschat, 2009
VectorData(const VectorSelector &);
/**
- * Initialize with a NamedData
- * object and cache the indices
- * in the VectorSelector base
- * class.
- *
- * @note Make sure the
- * VectorSelector base class
- * was filled with reasonable
- * data before calling this
- * function.
+ * @deprecated Use AnyData instead of NamedData.
*/
void initialize(const NamedData<VECTOR *> &);
/**
- * Initialize with a single
- * vector and cache the indices
- * in the VectorSelector base
- * class.
+ * Initialize with a single vector and cache the indices in the
+ * VectorSelector base class.
*
- * @note Make sure the
- * VectorSelector base class
- * was filled with reasonable
- * data before calling this
- * function.
+ * @note Make sure the VectorSelector base class was filled with
+ * reasonable data before calling this function.
*/
void initialize(const VECTOR *, const std::string &name);
const unsigned int start,
const unsigned int size) const;
- /**
+ virtual void mg_fill(
+ std::vector<std::vector<std::vector<double> > > &values,
+ std::vector<std::vector<std::vector<Tensor<1,dim> > > > &gradients,
+ std::vector<std::vector<std::vector<Tensor<2,dim> > > > &hessians,
+ const FEValuesBase<dim,spacedim> &fe,
+ const unsigned int level,
+ const std::vector<types::global_dof_index> &index,
+ const unsigned int component,
+ const unsigned int n_comp,
+ const unsigned int start,
+ const unsigned int size) const;
+
+ /**
* The memory used by this object.
*/
std::size_t memory_consumption () const;
-
- private:
- NamedData<SmartPointer<const VECTOR,VectorData<VECTOR,dim,spacedim> > > data;
};
/**
* Based on VectorSelector, this is the class that implements the
* function VectorDataBase::fill() for a certain type of multilevel vectors, using
- * NamedData to identify vectors by name.
+ * AnyData to identify vectors by name.
*
* @ingroup MeshWorker
* @author Guido Kanschat, 2010
*/
MGVectorData();
/**
- * Constructor using a
- * prefilled VectorSelector
+ * Constructor using a prefilled VectorSelector
*/
MGVectorData(const VectorSelector &);
/**
- * Initialize with a NamedData
- * object and cache the indices
- * in the VectorSelector base
- * class.
- *
- * @note Make sure the
- * VectorSelector base class
- * was filled with reasonable
- * data before calling this
- * function.
+ * @deprecated Use AnyData instead of NamedData.
*/
void initialize(const NamedData<MGLevelObject<VECTOR>*> &);
/**
- * Initialize with a single
- * vector and cache the indices
- * in the VectorSelector base
- * class.
+ * Initialize with a single vector and cache the indices in the
+ * VectorSelector base class.
*
- * @note Make sure the
- * VectorSelector base class
- * was filled with reasonable
- * data before calling this
- * function.
+ * @note Make sure the VectorSelector base class was filled with
+ * reasonable data before calling this function.
*/
void initialize(const MGLevelObject<VECTOR> *, const std::string &name);
- virtual void mg_fill(
- std::vector<std::vector<std::vector<double> > > &values,
- std::vector<std::vector<std::vector<Tensor<1,dim> > > > &gradients,
- std::vector<std::vector<std::vector<Tensor<2,dim> > > > &hessians,
- const FEValuesBase<dim,spacedim> &fe,
- const unsigned int level,
- const std::vector<types::global_dof_index> &index,
- const unsigned int component,
- const unsigned int n_comp,
- const unsigned int start,
- const unsigned int size) const;
/**
* The memory used by this object.
*/
std::size_t memory_consumption () const;
-
- private:
- NamedData<SmartPointer<const MGLevelObject<VECTOR>,MGVectorData<VECTOR,dim,spacedim> > > data;
};
//}
+ inline void
+ VectorSelector::initialize(const AnyData &src)
+ {
+ value_selection.initialize(src);
+ gradient_selection.initialize(src);
+ hessian_selection.initialize(src);
+ }
+
template <typename DATA>
inline void
VectorSelector::initialize(const NamedData<DATA> &src)
}
+ template <class STREAM, typename DATA>
+ inline void
+ VectorSelector::print(STREAM &s, const AnyData &v) const
+ {
+ s << "values: ";
+ for (unsigned int i=0; i<n_values(); ++i)
+ s << " '" << v.name(value_selection(i)) << '\'';
+ s << std::endl << "gradients:";
+ for (unsigned int i=0; i<n_gradients(); ++i)
+ s << " '" << v.name(gradient_selection(i)) << '\'';
+ s << std::endl << "hessians: ";
+ for (unsigned int i=0; i<n_hessians(); ++i)
+ s << " '" << v.name(hessian_selection(i)) << '\'';
+ s << std::endl;
+ }
+
+
template <class STREAM, typename DATA>
inline void
VectorSelector::print(STREAM &s, const NamedData<DATA> &v) const
{}
+ template <int dim, int spacedim>
+ void
+ VectorDataBase<dim, spacedim>::initialize(const AnyData &d)
+ {
+ this->data = d;
+ VectorSelector::initialize(d);
+ }
+
+
template <int dim, int spacedim>
void
VectorDataBase<dim, spacedim>::fill(
void
VectorData<VECTOR, dim, spacedim>::initialize(const NamedData<VECTOR *> &d)
{
- data = d;
+ this->data = d;
VectorSelector::initialize(d);
}
VectorData<VECTOR, dim, spacedim>::initialize(const VECTOR *v, const std::string &name)
{
SmartPointer<const VECTOR,VectorData<VECTOR, dim, spacedim> > p = v;
- data.add(p, name);
- VectorSelector::initialize(data);
+ this->data.add(p, name);
+ VectorSelector::initialize(this->data);
}
AssertDimension(gradients.size(), this->n_gradients());
AssertDimension(hessians.size(), this->n_hessians());
+ const AnyData& data = this->data;
for (unsigned int i=0; i<this->n_values(); ++i)
{
- const VECTOR &src = *data(this->value_index(i));
+ const VECTOR* src = data.read_ptr<VECTOR>(this->value_index(i));
VectorSlice<std::vector<std::vector<double> > > dst(values[i], component, n_comp);
- fe.get_function_values(src, make_slice(index, start, size), dst, true);
+ fe.get_function_values(*src, make_slice(index, start, size), dst, true);
}
for (unsigned int i=0; i<this->n_gradients(); ++i)
{
- const VECTOR &src = *data(this->gradient_index(i));
+ const VECTOR* src = data.read_ptr<VECTOR>(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);
+ 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->hessian_index(i));
+ const VECTOR* src = data.read_ptr<VECTOR>(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);
+ fe.get_function_hessians(*src, make_slice(index, start, size), dst, true);
}
}
VectorData<VECTOR, dim, spacedim>::memory_consumption () const
{
std::size_t mem = VectorSelector::memory_consumption();
- mem += sizeof (data);
+ mem += sizeof (this->data);
return mem;
}
void
MGVectorData<VECTOR, dim, spacedim>::initialize(const NamedData<MGLevelObject<VECTOR>*> &d)
{
- data = d;
+ this->data = d;
VectorSelector::initialize(d);
}
{
SmartPointer<const MGLevelObject<VECTOR>, MGVectorData<VECTOR, dim, spacedim> >
p = v;
- data.add(p, name);
- VectorSelector::initialize(data);
+ this->data.add(p, name);
+ VectorSelector::initialize(this->data);
}
template <class VECTOR, int dim, int spacedim>
void
- MGVectorData<VECTOR, dim, spacedim>::mg_fill(
+ VectorData<VECTOR, dim, spacedim>::mg_fill(
std::vector<std::vector<std::vector<double> > > &values,
std::vector<std::vector<std::vector<Tensor<1,dim> > > > &gradients,
std::vector<std::vector<std::vector<Tensor<2,dim> > > > &hessians,
AssertDimension(gradients.size(), this->n_gradients());
AssertDimension(hessians.size(), this->n_hessians());
+ const AnyData& data = this->data;
for (unsigned int i=0; i<this->n_values(); ++i)
{
- const VECTOR &src = (*data(this->value_index(i)))[level];
+ const MGLevelObject<VECTOR>* src = data.read_ptr<MGLevelObject<VECTOR> >(this->value_index(i));
VectorSlice<std::vector<std::vector<double> > > dst(values[i], component, n_comp);
- fe.get_function_values(src, make_slice(index, start, size), dst, true);
+ fe.get_function_values((*src)[level], make_slice(index, start, size), dst, true);
}
for (unsigned int i=0; i<this->n_gradients(); ++i)
{
- const VECTOR &src = (*data(this->value_index(i)))[level];
+ const MGLevelObject<VECTOR>* src = data.read_ptr<MGLevelObject<VECTOR> >(this->value_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);
+ fe.get_function_gradients((*src)[level], 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)))[level];
+ const MGLevelObject<VECTOR>* src = data.read_ptr<MGLevelObject<VECTOR> >(this->value_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);
+ fe.get_function_hessians((*src)[level], make_slice(index, start, size), dst, true);
}
}
MGVectorData<VECTOR, dim, spacedim>::memory_consumption () const
{
std::size_t mem = VectorSelector::memory_consumption();
- mem += sizeof (data);
+ mem += sizeof (this->data);
return mem;
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: theta_timestepping.cc 30050 2013-07-18 20:03:31Z maier $
+//
+// Copyright (C) 2005 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// See documentation of ThetaTimestepping for documentation of this example
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/algorithms/operator.h>
+#include <deal.II/algorithms/theta_timestepping.h>
+
+#include <iostream>
+
+using namespace dealii;
+using namespace Algorithms;
+
+
+class Explicit
+ : public Operator<Vector<double> >
+{
+public:
+ Explicit(const FullMatrix<double> &matrix);
+ void operator() (AnyData &out, const AnyData &in);
+private:
+ SmartPointer<const FullMatrix<double>, Explicit> matrix;
+ FullMatrix<double> m;
+};
+
+
+class Implicit
+ : public Operator<Vector<double> >
+{
+public:
+ Implicit(const FullMatrix<double> &matrix);
+ void operator() (AnyData &out, const AnyData &in);
+private:
+ SmartPointer<const FullMatrix<double>, Implicit> matrix;
+ FullMatrix<double> m;
+};
+
+// End of declarations
+
+int main()
+{
+ std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ FullMatrix<double> matrix(2);
+ matrix(0,0) = 0.;
+ matrix(1,1) = 0.;
+ matrix(0,1) = numbers::PI;
+ matrix(1,0) = -numbers::PI;
+
+ OutputOperator<Vector<double> > out;
+ out.initialize_stream(logfile);
+
+ Explicit op_explicit(matrix);
+ Implicit op_implicit(matrix);
+ ThetaTimestepping<Vector<double> > solver(op_explicit, op_implicit);
+ solver.timestep_control().start_step(0.1);
+ //solver.set_output(out);
+
+ Vector<double> value(2);
+ value(0) = 1.;
+ AnyData indata;
+ AnyData outdata;
+ Vector<double> *p = &value;
+ outdata.add(p, "value");
+ deallog << "Initial: " << value(0) << ' ' << value(1) << std::endl;
+ solver.notify(Events::initial);
+ solver(outdata, indata);
+ deallog << "Result: " << value(0) << ' ' << value(1)
+ << " Norm " << value.l2_norm() << std::endl;
+}
+
+
+Explicit::Explicit(const FullMatrix<double> &M)
+ :
+ matrix(&M)
+{
+ m.reinit(M.m(), M.n());
+}
+
+
+void
+Explicit::operator() (AnyData &out, const AnyData &in)
+{
+ const double* step = in.read_ptr<double>("Timestep");
+
+ if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size))
+ {
+ m.equ(-*step, *matrix);
+ for (unsigned int i=0; i<m.m(); ++i)
+ m(i,i) += 1.;
+ }
+ this->notifications.clear();
+ unsigned int i = in.find("Previous iterate");
+ m.vmult(*out.entry<Vector<double>*>(0), *in.read_ptr<Vector<double> >(i));
+}
+
+
+Implicit::Implicit(const FullMatrix<double> &M)
+ :
+ matrix(&M)
+{
+ m.reinit(M.m(), M.n());
+}
+
+
+void
+Implicit::operator() (AnyData &out, const AnyData &in)
+{
+ const double * step = in.read_ptr<double>("Timestep");
+
+ if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size))
+ {
+ m.equ(*step, *matrix);
+ for (unsigned int i=0; i<m.m(); ++i)
+ m(i,i) += 1.;
+ m.gauss_jordan();
+ }
+ this->notifications.clear();
+
+ unsigned int i = in.find("Previous time");
+ m.vmult(*out.entry<Vector<double>*>(0), *in.read_ptr<Vector<double> >(i));
+}
+
+
--- /dev/null
+
+DEAL::Initial: 1.00000 0
+DEAL:Theta::Time step:0.100000
+DEAL:Theta::Time step:0.200000
+DEAL:Theta::Time step:0.300000
+DEAL:Theta::Time step:0.400000
+DEAL:Theta::Time step:0.500000
+DEAL:Theta::Time step:0.600000
+DEAL:Theta::Time step:0.700000
+DEAL:Theta::Time step:0.800000
+DEAL:Theta::Time step:0.900000
+DEAL:Theta::Time step:1.00000
+DEAL::Result: -0.999676 0.0254599 Norm 1.00000
projects / dealii-svn.git / commitdiff