From: Denis Davydov Date: Tue, 23 Oct 2018 11:05:07 +0000 (+0200) Subject: add IncrementalFunction X-Git-Tag: v9.1.0-rc1~588^2~8 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=575ddd9e13399f856d5884bf81abfb651fb4bf88;p=dealii.git add IncrementalFunction --- diff --git a/doc/news/changes/minor/20181023DenisDavydov b/doc/news/changes/minor/20181023DenisDavydov new file mode 100644 index 0000000000..e96dcd0221 --- /dev/null +++ b/doc/news/changes/minor/20181023DenisDavydov @@ -0,0 +1,4 @@ +New: Add Functions::IncrementalFunction which returns the difference in values of +a given function at two time steps. +
+(Denis Davydov, Jean-Paul Pelteret, 2018/10/05) diff --git a/include/deal.II/base/incremental_function.h b/include/deal.II/base/incremental_function.h new file mode 100644 index 0000000000..00e30a6eec --- /dev/null +++ b/include/deal.II/base/incremental_function.h @@ -0,0 +1,119 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_incremental_function_h +#define dealii_incremental_function_h + + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template +class Vector; + +namespace Functions +{ + /** + * This class represents an incremental function. That is, given arbitrary + * function func, this class will return + * f(t)-f(t-delta). The decrement is set by the method + * set_decrement(). The main application is to transform a given Dirichlet + * boundary condition function into the incremental form. + * + * @ingroup functions + * @author Denis Davydov, Jean-Paul Pelteret, 2018 + */ + template + class IncrementalFunction : public Function + { + public: + /** + * Export the value of the template parameter as a static member constant. + * Sometimes useful for some expression template programming. + */ + static const unsigned int dimension = dim; + + /** + * The scalar-valued real type used for representing time. + */ + using time_type = typename Function::time_type; + + /** + * Constructor which wraps a given function @p base. + * + * @note this class stores a non-constant reference to @p base + * and will call base.set_time() during evaluation. + */ + IncrementalFunction(Function &base); + + /** + * Virtual destructor + */ + virtual ~IncrementalFunction() = default; + + /** + * Return the value of the function at the given point. Unless there is only + * one component (i.e. the function is scalar), you should state the + * component you want to have evaluated; it defaults to zero, i.e. the first + * component. + */ + virtual RangeNumberType + value(const Point &p, const unsigned int component = 0) const override; + + /** + * Return all components of a vector-valued function at a given point. + * + * values shall have the right size beforehand. + */ + virtual void + vector_value(const Point & p, + Vector &values) const override; + + /** + * Set (positive) time decrement. + */ + void + set_decrement(const time_type delta_t); + + private: + /** + * Reference to the function being wrapped. + */ + Function &base; + + /** + * Time decrement. + */ + time_type delta_t; + + /** + * Auxiliary vector to store values + */ + mutable Vector values_old; + + /** + * Thread mutex + */ + mutable Threads::Mutex mutex; + }; + +} // namespace Functions + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 8764dc271c..8e03e013c5 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -39,6 +39,7 @@ SET(_unity_include_src geometry_info.cc geometric_utilities.cc graph_coloring.cc + incremental_function.cc index_set.cc job_identifier.cc logstream.cc @@ -112,6 +113,7 @@ SET(_inst data_out_base.inst.in function.inst.in function_time.inst.in + incremental_function.inst.in geometric_utilities.inst.in mpi.inst.in partitioner.inst.in diff --git a/source/base/incremental_function.cc b/source/base/incremental_function.cc new file mode 100644 index 0000000000..d65d479e9e --- /dev/null +++ b/source/base/incremental_function.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2007 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions +{ + template + IncrementalFunction::IncrementalFunction( + Function &base) + : base(base) + , delta_t(time_type()) + , values_old(base.n_components) + {} + + + + template + void + IncrementalFunction::set_decrement( + const time_type delta) + { + delta_t = delta; + } + + + + template + void + IncrementalFunction::vector_value( + const Point & p, + Vector &values) const + { + // since we modify a mutable member variable, lock the + // the data via a mutex + std::lock_guard lock(mutex); + + base.set_time(this->get_time()); + base.vector_value(p, values); + + base.set_time(this->get_time() - delta_t); + base.vector_value(p, values_old); + + values -= values_old; + } + + + + template + RangeNumberType + IncrementalFunction::value(const Point &p, + unsigned int comp) const + { + base.set_time(this->get_time()); + const RangeNumberType current = base.value(p, comp); + + base.set_time(this->get_time() - delta_t); + const RangeNumberType old = base.value(p, comp); + + return current - old; + } + + +// Explicit instantiations +#include "incremental_function.inst" +} // namespace Functions +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/incremental_function.inst.in b/source/base/incremental_function.inst.in new file mode 100644 index 0000000000..e02fa9553c --- /dev/null +++ b/source/base/incremental_function.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2017 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +for (S : REAL_SCALARS; dim : SPACE_DIMENSIONS) + { + template class Functions::IncrementalFunction; + } + +for (S : COMPLEX_SCALARS; dim : SPACE_DIMENSIONS) + { + template class Functions::IncrementalFunction; + } diff --git a/tests/base/functions_incremental.cc b/tests/base/functions_incremental.cc new file mode 100644 index 0000000000..099a2f521b --- /dev/null +++ b/tests/base/functions_incremental.cc @@ -0,0 +1,103 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test IncrementalFunction + +#include + +#include "../tests.h" + +// homogeneous linear in time function +template +class MyFunc : public Function +{ +public: + MyFunc(const double a, const double b) + : Function(dim) + , a(a) + , b(b) + {} + + virtual double + value(const Point /*p*/, const unsigned int component = 0) const + { + (void)component; + return a * this->get_time() + b; + } + + virtual void + vector_value(const Point & /*p*/, Vector &values) const override + { + for (unsigned int i = 0; i < dim; ++i) + values[i] = a * this->get_time() + b; + } + +private: + const double a; + const double b; +}; + +template +void +check() +{ + Vector v1(dim); + Vector v2(dim); + Vector v3(dim); + + const Point point; + + MyFunc func(0.5, 10); + Functions::IncrementalFunction inc(func); + + // default delta is zero: + inc.set_time(1.0); + + std::vector> time_delta = {{1.0, 0.}, {1.0, 0.2}}; + + for (auto &el : time_delta) + { + deallog << "time = " << el.first << " delta = " << el.second << std::endl; + func.set_time(el.first); + func.vector_value(point, v1); + v1.print(deallog); + deallog << "-" << std::endl; + + func.set_time(el.first - el.second); + func.vector_value(point, v2); + v2.print(deallog); + deallog << "=" << std::endl; + + inc.set_time(el.first); + inc.set_decrement(el.second); + inc.vector_value(point, v3); + v3.print(deallog); + + v1.add(-1, v2); + v1.add(-1, v3); + AssertThrow(v1.l2_norm() == 0., ExcInternalError()); + } + + deallog << "Ok" << std::endl; +} + +int +main() +{ + initlog(); + + check<2>(); +} diff --git a/tests/base/functions_incremental.output b/tests/base/functions_incremental.output new file mode 100644 index 0000000000..7956c1eb52 --- /dev/null +++ b/tests/base/functions_incremental.output @@ -0,0 +1,14 @@ + +DEAL::time = 1.00000 delta = 0.00000 +DEAL::10.5000 10.5000 +DEAL::- +DEAL::10.5000 10.5000 +DEAL::= +DEAL::0.00000 0.00000 +DEAL::time = 1.00000 delta = 0.200000 +DEAL::10.5000 10.5000 +DEAL::- +DEAL::10.4000 10.4000 +DEAL::= +DEAL::0.100000 0.100000 +DEAL::Ok