From c242b71a31a7cf8623a3152c01df7233261ab773 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 20 May 2014 19:05:42 +0000 Subject: [PATCH] merge any_data branch git-svn-id: https://svn.dealii.org/trunk@32947 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/algorithms/theta_timestepping.h | 2 +- .../algorithms/theta_timestepping.templates.h | 8 +- .../include/deal.II/meshworker/functional.h | 113 ++++--- .../deal.II/meshworker/integration_info.h | 55 +++- deal.II/include/deal.II/meshworker/simple.h | 44 ++- .../deal.II/meshworker/vector_selector.h | 302 ++++++++---------- .../meshworker/vector_selector.templates.h | 53 +-- tests/algorithms/theta_01.cc | 146 +++++++++ tests/algorithms/theta_01.output | 13 + 9 files changed, 462 insertions(+), 274 deletions(-) create mode 100644 tests/algorithms/theta_01.cc create mode 100644 tests/algorithms/theta_01.output diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.h b/deal.II/include/deal.II/algorithms/theta_timestepping.h index e6bb52410a..0fc7df05cf 100644 --- a/deal.II/include/deal.II/algorithms/theta_timestepping.h +++ b/deal.II/include/deal.II/algorithms/theta_timestepping.h @@ -114,7 +114,7 @@ namespace Algorithms * 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 - * const double reference named "Time" and "Timestep". Note + * const double* 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. diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h index 4072b37b4e..ab0d36a243 100644 --- a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h +++ b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h @@ -91,8 +91,8 @@ namespace Algorithms // timestep AnyData src1; src1.add(&solution, "Previous iterate"); - src1.add(d_explicit.time, "Time"); - src1.add(d_explicit.step, "Timestep"); + src1.add(&d_explicit.time, "Time"); + src1.add(&d_explicit.step, "Timestep"); src1.merge(in); AnyData src2; @@ -102,8 +102,8 @@ namespace Algorithms out1.add(aux, "Solution"); // The data provided to the inner solver src2.add(aux, "Previous time"); - src2.add(d_implicit.time, "Time"); - src2.add(d_explicit.step, "Timestep"); + src2.add(&d_implicit.time, "Time"); + src2.add(&d_explicit.step, "Timestep"); src2.merge(in); if (output != 0) diff --git a/deal.II/include/deal.II/meshworker/functional.h b/deal.II/include/deal.II/meshworker/functional.h index a651ba1acc..0205951f21 100644 --- a/deal.II/include/deal.II/meshworker/functional.h +++ b/deal.II/include/deal.II/meshworker/functional.h @@ -19,6 +19,7 @@ #define __deal2__mesh_worker_functional_h #include +#include #include #include #include @@ -107,49 +108,28 @@ namespace MeshWorker { public: /** - * The data type for - * communicating the cell and - * face vectors. - */ - typedef NamedData*> 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*> &results, bool separate_faces = true); /** * Initialize the local data @@ -188,7 +168,7 @@ namespace MeshWorker */ number operator() (const unsigned int i) const; private: - DataVectors results; + AnyData results; bool separate_faces; }; //----------------------------------------------------------------------// @@ -247,13 +227,29 @@ namespace MeshWorker template inline void - CellsAndFaces::initialize(DataVectors &r, bool sep) + CellsAndFaces::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*>(0)->n_blocks(), + r.entry*>(1)->n_blocks()); + } + + results = r; + separate_faces = sep; + } + + template + inline void + CellsAndFaces::initialize(NamedData*> &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; @@ -266,7 +262,7 @@ namespace MeshWorker inline void CellsAndFaces::initialize_info(DOFINFO &info, bool) const { - info.initialize_numbers(results(0)->n_blocks()); + info.initialize_numbers(results.entry*>(0)->n_blocks()); } @@ -275,14 +271,15 @@ namespace MeshWorker inline void CellsAndFaces::assemble(const DOFINFO &info) { + BlockVector* v; + if (separate_faces && + info.face_number != deal_II_numbers::invalid_unsigned_int) + v = results.entry*>(1); + else + v = results.entry*>(0); + for (unsigned int i=0; iblock(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); } @@ -296,15 +293,17 @@ namespace MeshWorker { if (separate_faces) { + BlockVector* v1 = results.entry*>(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* v0 = results.entry*>(0); + v0->block(i)(info1.cell->user_index()) += .5*info1.value(i); + v0->block(i)(info2.cell->user_index()) += .5*info2.value(i); } } } diff --git a/deal.II/include/deal.II/meshworker/integration_info.h b/deal.II/include/deal.II/meshworker/integration_info.h index 4eb299711e..a120eced0f 100644 --- a/deal.II/include/deal.II/meshworker/integration_info.h +++ b/deal.II/include/deal.II/meshworker/integration_info.h @@ -248,7 +248,7 @@ namespace MeshWorker * 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. * @@ -313,17 +313,22 @@ namespace MeshWorker * formulas, which integrate polynomial bilinear forms exactly. */ template + void initialize(const FiniteElement &el, + const Mapping &mapping, + const AnyData &data, + const VECTOR& dummy, + const BlockInfo *block_info = 0); + /** + * @deprecated Use AnyData instead of NamedData. + */ + template void initialize(const FiniteElement &el, const Mapping &mapping, const NamedData &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 void initialize(const FiniteElement &el, @@ -767,6 +772,44 @@ namespace MeshWorker } + template + template + void + IntegrationInfoBox::initialize( + const FiniteElement &el, + const Mapping &mapping, + const AnyData &data, + const VECTOR& dummy, + const BlockInfo *block_info) + { + initialize(el, mapping, block_info); + std_cxx1x::shared_ptr > p; + VectorDataBase* pp; + + p = std_cxx1x::shared_ptr >(new VectorData (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 >(new VectorData (boundary_selector)); + pp = &*p; + pp->initialize(data); + boundary_data = p; + boundary.initialize_data(p); + + p = std_cxx1x::shared_ptr >(new VectorData (face_selector)); + pp = &*p; + pp->initialize(data); + face_data = p; + face.initialize_data(p); + subface.initialize_data(p); + neighbor.initialize_data(p); + } + + template template void diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h index ed692c57f3..297e2421a6 100644 --- a/deal.II/include/deal.II/meshworker/simple.h +++ b/deal.II/include/deal.II/meshworker/simple.h @@ -19,6 +19,7 @@ #define __deal2__mesh_worker_simple_h #include +#include #include #include #include @@ -60,8 +61,18 @@ namespace MeshWorker class ResidualSimple { public: - void initialize(NamedData &results); + /** + * Initialize with an AnyData object holding the result of + * assembling. + * + * Assembling currently writes into the first vector of results. + */ + void initialize(AnyData& results); + /** + * @deprecated Use initialize(AnyData&) instead. + */ + void initialize(NamedData &results); /** * Initialize the constraints. */ @@ -114,7 +125,7 @@ namespace MeshWorker /** * The global residal vectors filled by assemble(). */ - NamedData > > residuals; + AnyData residuals; /** * A pointer to the object containing constraints. */ @@ -500,6 +511,13 @@ namespace MeshWorker //----------------------------------------------------------------------// + template + inline void + ResidualSimple::initialize(AnyData &results) + { + residuals = results; + } + template inline void ResidualSimple::initialize(NamedData &results) @@ -537,18 +555,19 @@ namespace MeshWorker { for (unsigned int k=0; k(k); if (constraints == 0) { for (unsigned int i=0; idistribute_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); } } } @@ -561,30 +580,31 @@ namespace MeshWorker { for (unsigned int k=0; k(k); if (constraints == 0) { for (unsigned int i=0; idistribute_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; idistribute_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 @@ -1132,7 +1152,7 @@ namespace MeshWorker inline void SystemSimple::initialize(MATRIX &m, VECTOR &rhs) { - NamedData data; + AnyData data; VECTOR *p = &rhs; data.add(p, "right hand side"); diff --git a/deal.II/include/deal.II/meshworker/vector_selector.h b/deal.II/include/deal.II/meshworker/vector_selector.h index 5976df4dbd..ff4d2e2fd4 100644 --- a/deal.II/include/deal.II/meshworker/vector_selector.h +++ b/deal.II/include/deal.II/meshworker/vector_selector.h @@ -18,6 +18,7 @@ #define __deal2__mesh_worker_vector_selector_h #include +#include #include #include #include @@ -32,7 +33,7 @@ namespace MeshWorker /** * 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. @@ -49,23 +50,13 @@ namespace MeshWorker { 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, @@ -73,10 +64,8 @@ namespace MeshWorker 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, @@ -85,21 +74,19 @@ namespace MeshWorker // 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 void initialize(const NamedData &); @@ -110,20 +97,17 @@ namespace MeshWorker 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; @@ -162,6 +146,12 @@ namespace MeshWorker * selection to the stream. */ template + void print (STREAM &s, const AnyData &v) const; + + /** + * @deprecated Use AnyData instead of NamedData. + */ + template void print (STREAM &s, const NamedData &v) const; /** @@ -183,14 +173,12 @@ namespace MeshWorker 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; }; @@ -213,81 +201,61 @@ namespace MeshWorker 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( @@ -302,15 +270,10 @@ namespace MeshWorker 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 > > &values, @@ -328,13 +291,15 @@ namespace MeshWorker * 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 @@ -355,30 +320,16 @@ namespace MeshWorker 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 &); /** - * 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); @@ -393,20 +344,29 @@ namespace MeshWorker const unsigned int start, const unsigned int size) const; - /** + virtual void mg_fill( + std::vector > > &values, + std::vector > > > &gradients, + std::vector > > > &hessians, + const FEValuesBase &fe, + const unsigned int level, + const std::vector &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 > > 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 @@ -421,58 +381,29 @@ namespace MeshWorker */ 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*> &); /** - * 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 *, const std::string &name); - virtual void mg_fill( - std::vector > > &values, - std::vector > > > &gradients, - std::vector > > > &hessians, - const FEValuesBase &fe, - const unsigned int level, - const std::vector &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,MGVectorData > > data; }; @@ -504,6 +435,14 @@ namespace MeshWorker //} + inline void + VectorSelector::initialize(const AnyData &src) + { + value_selection.initialize(src); + gradient_selection.initialize(src); + hessian_selection.initialize(src); + } + template inline void VectorSelector::initialize(const NamedData &src) @@ -596,6 +535,23 @@ namespace MeshWorker } + template + inline void + VectorSelector::print(STREAM &s, const AnyData &v) const + { + s << "values: "; + for (unsigned int i=0; i inline void VectorSelector::print(STREAM &s, const NamedData &v) const diff --git a/deal.II/include/deal.II/meshworker/vector_selector.templates.h b/deal.II/include/deal.II/meshworker/vector_selector.templates.h index dc0e6fbbff..c48e19e25d 100644 --- a/deal.II/include/deal.II/meshworker/vector_selector.templates.h +++ b/deal.II/include/deal.II/meshworker/vector_selector.templates.h @@ -40,6 +40,15 @@ namespace MeshWorker {} + template + void + VectorDataBase::initialize(const AnyData &d) + { + this->data = d; + VectorSelector::initialize(d); + } + + template void VectorDataBase::fill( @@ -93,7 +102,7 @@ namespace MeshWorker void VectorData::initialize(const NamedData &d) { - data = d; + this->data = d; VectorSelector::initialize(d); } @@ -103,8 +112,8 @@ namespace MeshWorker VectorData::initialize(const VECTOR *v, const std::string &name) { SmartPointer > p = v; - data.add(p, name); - VectorSelector::initialize(data); + this->data.add(p, name); + VectorSelector::initialize(this->data); } @@ -125,25 +134,26 @@ namespace MeshWorker AssertDimension(gradients.size(), this->n_gradients()); AssertDimension(hessians.size(), this->n_hessians()); + const AnyData& data = this->data; for (unsigned int i=0; in_values(); ++i) { - const VECTOR &src = *data(this->value_index(i)); + const VECTOR* src = data.read_ptr(this->value_index(i)); VectorSlice > > 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; in_gradients(); ++i) { - const VECTOR &src = *data(this->gradient_index(i)); + const VECTOR* src = data.read_ptr(this->gradient_index(i)); VectorSlice > > > 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; in_hessians(); ++i) { - const VECTOR &src = *data(this->hessian_index(i)); + const VECTOR* src = data.read_ptr(this->hessian_index(i)); VectorSlice > > > 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); } } @@ -153,7 +163,7 @@ namespace MeshWorker VectorData::memory_consumption () const { std::size_t mem = VectorSelector::memory_consumption(); - mem += sizeof (data); + mem += sizeof (this->data); return mem; } @@ -175,7 +185,7 @@ namespace MeshWorker void MGVectorData::initialize(const NamedData*> &d) { - data = d; + this->data = d; VectorSelector::initialize(d); } @@ -186,14 +196,14 @@ namespace MeshWorker { SmartPointer, MGVectorData > p = v; - data.add(p, name); - VectorSelector::initialize(data); + this->data.add(p, name); + VectorSelector::initialize(this->data); } template void - MGVectorData::mg_fill( + VectorData::mg_fill( std::vector > > &values, std::vector > > > &gradients, std::vector > > > &hessians, @@ -209,25 +219,26 @@ namespace MeshWorker AssertDimension(gradients.size(), this->n_gradients()); AssertDimension(hessians.size(), this->n_hessians()); + const AnyData& data = this->data; for (unsigned int i=0; in_values(); ++i) { - const VECTOR &src = (*data(this->value_index(i)))[level]; + const MGLevelObject* src = data.read_ptr >(this->value_index(i)); VectorSlice > > 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; in_gradients(); ++i) { - const VECTOR &src = (*data(this->value_index(i)))[level]; + const MGLevelObject* src = data.read_ptr >(this->value_index(i)); VectorSlice > > > 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; in_hessians(); ++i) { - const VECTOR &src = (*data(this->value_index(i)))[level]; + const MGLevelObject* src = data.read_ptr >(this->value_index(i)); VectorSlice > > > 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); } } @@ -237,7 +248,7 @@ namespace MeshWorker MGVectorData::memory_consumption () const { std::size_t mem = VectorSelector::memory_consumption(); - mem += sizeof (data); + mem += sizeof (this->data); return mem; } } diff --git a/tests/algorithms/theta_01.cc b/tests/algorithms/theta_01.cc new file mode 100644 index 0000000000..bb18059d1b --- /dev/null +++ b/tests/algorithms/theta_01.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// $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 +#include +#include + +#include +#include + +#include + +using namespace dealii; +using namespace Algorithms; + + +class Explicit + : public Operator > +{ +public: + Explicit(const FullMatrix &matrix); + void operator() (AnyData &out, const AnyData &in); +private: + SmartPointer, Explicit> matrix; + FullMatrix m; +}; + + +class Implicit + : public Operator > +{ +public: + Implicit(const FullMatrix &matrix); + void operator() (AnyData &out, const AnyData &in); +private: + SmartPointer, Implicit> matrix; + FullMatrix 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 matrix(2); + matrix(0,0) = 0.; + matrix(1,1) = 0.; + matrix(0,1) = numbers::PI; + matrix(1,0) = -numbers::PI; + + OutputOperator > out; + out.initialize_stream(logfile); + + Explicit op_explicit(matrix); + Implicit op_implicit(matrix); + ThetaTimestepping > solver(op_explicit, op_implicit); + solver.timestep_control().start_step(0.1); + //solver.set_output(out); + + Vector value(2); + value(0) = 1.; + AnyData indata; + AnyData outdata; + Vector *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 &M) + : + matrix(&M) +{ + m.reinit(M.m(), M.n()); +} + + +void +Explicit::operator() (AnyData &out, const AnyData &in) +{ + const double* step = in.read_ptr("Timestep"); + + if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size)) + { + m.equ(-*step, *matrix); + for (unsigned int i=0; inotifications.clear(); + unsigned int i = in.find("Previous iterate"); + m.vmult(*out.entry*>(0), *in.read_ptr >(i)); +} + + +Implicit::Implicit(const FullMatrix &M) + : + matrix(&M) +{ + m.reinit(M.m(), M.n()); +} + + +void +Implicit::operator() (AnyData &out, const AnyData &in) +{ + const double * step = in.read_ptr("Timestep"); + + if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size)) + { + m.equ(*step, *matrix); + for (unsigned int i=0; inotifications.clear(); + + unsigned int i = in.find("Previous time"); + m.vmult(*out.entry*>(0), *in.read_ptr >(i)); +} + + diff --git a/tests/algorithms/theta_01.output b/tests/algorithms/theta_01.output new file mode 100644 index 0000000000..7c4d9ddfc7 --- /dev/null +++ b/tests/algorithms/theta_01.output @@ -0,0 +1,13 @@ + +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 -- 2.39.5