#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
+#include <deal.II/base/smartpointer.h>
#include <deal.II/base/table.h>
#include <array>
const bool integrate_to_one = false,
const double unitary_integral_value = 1.0);
+ /**
+ * Virtual destructor.
+ */
+ virtual ~CutOffFunctionBase() = default;
+
/**
* Move the center of the ball to new point <tt>p</tt>.
*
/**
* Set the center of the ball to the point @p p.
*/
- void
+ virtual void
set_center(const Point<dim> &p);
/**
*
* @deprecated Use set_radius() instead.
*/
- void
+ virtual void
set_radius(const double r);
/**
double
get_radius() const;
+ /**
+ * Return a boolean indicating if this function integrates to one.
+ */
+ bool
+ integrates_to_one() const;
+
protected:
/**
* Center of the integration ball.
/**
* Flag that controls wether we rescale the value when the radius changes.
*/
- const bool integrate_to_one;
+ bool integrate_to_one;
/**
* The reference integral value. Derived classes should specify what their
};
+ /**
+ * Tensor product of CutOffFunctionBase objects.
+ *
+ * Instead of using the distance to compute the cut-off function, this class
+ * performs a tensor product of the same CutOffFunctionBase object in each
+ * coordinate direction.
+ *
+ * @ingroup functions
+ * @author Luca Heltai, 2019.
+ */
+ template <int dim>
+ class CutOffFunctionTensorProduct : public CutOffFunctionBase<dim>
+ {
+ public:
+ /**
+ * Construct an empty CutOffFunctionTensorProduct object.
+ *
+ * Before you can use this class, you have to call the set_base() method
+ * with a class derived from the CutOffFunctionBase object.
+ *
+ * If you try to use this class before you call the set_base() method,
+ * and exception will be triggered.
+ */
+ CutOffFunctionTensorProduct(
+ double radius = 1.0,
+ const Point<dim> & center = Point<dim>(),
+ const unsigned int n_components = 1,
+ const unsigned int select = CutOffFunctionBase<dim>::no_component,
+ const bool integrate_to_one = false);
+
+ /**
+ * Initialize the class with an objet of type
+ * @tparam CutOffFunctionBaseType<1>.
+ */
+ template <template <int> class CutOffFunctionBaseType>
+ void
+ set_base();
+
+ /**
+ * Set the new center.
+ */
+ virtual void
+ set_center(const Point<dim> ¢er) override;
+
+ /**
+ * Set the new radius.
+ */
+ virtual void
+ set_radius(const double radius) override;
+
+ /**
+ * Function value at one point.
+ */
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+ /**
+ * Function gradient at one point.
+ */
+ virtual Tensor<1, dim>
+ gradient(const Point<dim> & p,
+ const unsigned int component = 0) const override;
+
+ private:
+ std::array<std::unique_ptr<CutOffFunctionBase<1>>, dim> base;
+
+ bool initialized;
+ };
+
+
/**
* Cut-off function in L-infinity for an arbitrary ball. This function is
const std::vector<double> coefficients;
};
+#ifndef DOXYGEN
+
+ // Template definitions
+ template <int dim>
+ template <template <int> class CutOffFunctionBaseType>
+ void
+ CutOffFunctionTensorProduct<dim>::set_base()
+ {
+ initialized = true;
+ for (unsigned int i = 0; i < dim; ++i)
+ base[i].reset(new CutOffFunctionBaseType<1>(this->radius,
+ Point<1>(this->center[i]),
+ this->n_components,
+ this->selected,
+ this->integrate_to_one));
+ }
+
+
+
+#endif
+
} // namespace Functions
DEAL_II_NAMESPACE_CLOSE
}
+
+ template <int dim>
+ bool
+ CutOffFunctionBase<dim>::integrates_to_one() const
+ {
+ return integrate_to_one;
+ }
+
+
+
+ template <int dim>
+ CutOffFunctionTensorProduct<dim>::CutOffFunctionTensorProduct(
+ double radius,
+ const Point<dim> & center,
+ const unsigned int n_components,
+ const unsigned int select,
+ const bool integrate_to_one)
+ : CutOffFunctionBase<dim>(radius,
+ center,
+ n_components,
+ select,
+ integrate_to_one)
+ , initialized(false)
+ {}
+
+
+
+ template <int dim>
+ void
+ CutOffFunctionTensorProduct<dim>::set_center(const Point<dim> &p)
+ {
+ Assert(initialized, ExcNotInitialized());
+ for (unsigned int i = 0; i < dim; ++i)
+ base[i]->set_center(Point<1>(p[i]));
+ CutOffFunctionBase<dim>::set_center(p);
+ }
+
+
+
+ template <int dim>
+ void
+ CutOffFunctionTensorProduct<dim>::set_radius(const double r)
+ {
+ Assert(initialized, ExcNotInitialized());
+ for (unsigned int i = 0; i < dim; ++i)
+ base[i]->set_radius(r);
+ CutOffFunctionBase<dim>::set_radius(r);
+ }
+
+
+
+ template <int dim>
+ double
+ CutOffFunctionTensorProduct<dim>::value(const Point<dim> & p,
+ const unsigned int component) const
+ {
+ Assert(initialized, ExcNotInitialized());
+ double ret = 1.0;
+ for (unsigned int i = 0; i < dim; ++i)
+ ret *= base[i]->value(Point<1>(p[i]), component);
+ return ret;
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ CutOffFunctionTensorProduct<dim>::gradient(const Point<dim> & p,
+ const unsigned int component) const
+ {
+ Assert(initialized, ExcNotInitialized());
+ Tensor<1, dim> ret;
+ for (unsigned int i = 0; i < dim; ++i)
+ ret[i] = base[i]->gradient(Point<1>(p[i]), component)[0];
+ return ret;
+ }
+
+
+
//////////////////////////////////////////////////////////////////////
namespace
{
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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 the functionality of the CutOffTensorProductFunction
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <map>
+
+#include "../tests.h"
+
+using namespace Functions;
+
+template <int dim>
+struct Domain
+{
+ Domain()
+ : fe(1)
+ , dh(tria)
+ , quad(5)
+ {
+ deallog << "Testing dim = " << dim << std::endl;
+ GridGenerator::hyper_cube(tria, -2, 2);
+ tria.refine_global(5 - dim);
+ dh.distribute_dofs(fe);
+ zero.reinit(dh.n_dofs());
+ cell_integral.reinit(tria.n_active_cells());
+ }
+
+ template <class Function>
+ void
+ integrate(const Function &fun)
+ {
+ deallog << "Integrating " << Utilities::type_to_string(fun)
+ << " -2,2 cube: ";
+
+ VectorTools::integrate_difference(
+ dh, zero, fun, cell_integral, quad, VectorTools::L1_norm);
+ const double integral =
+ VectorTools::compute_global_error(tria,
+ cell_integral,
+ VectorTools::L1_norm);
+
+ deallog << integral << std::endl;
+ }
+
+ Triangulation<dim> tria;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dh;
+ Vector<double> zero;
+ Vector<double> cell_integral;
+ QGauss<dim> quad;
+};
+
+template <int dim, template <int> class TestCutOffFunction>
+void
+test()
+{
+ static Domain<dim> domain;
+
+ CutOffFunctionTensorProduct<dim> fun(
+ 1, Point<dim>(), 1, CutOffFunctionBase<dim>::no_component, true);
+ fun.template set_base<TestCutOffFunction>();
+
+ deallog << "Center: " << fun.get_center() << std::endl
+ << "Radius: " << fun.get_radius() << std::endl;
+
+ domain.integrate(fun);
+
+ Point<dim> new_center;
+ for (unsigned int i = 0; i < dim; ++i)
+ new_center[i] = .5;
+
+ fun.set_center(new_center);
+ fun.set_radius(.5);
+
+ deallog << "Center: " << fun.get_center() << std::endl
+ << "Radius: " << fun.get_radius() << std::endl;
+
+ domain.integrate(fun);
+}
+
+int
+main()
+{
+ initlog(1);
+
+ test<1, Functions::CutOffFunctionLinfty>();
+ test<2, Functions::CutOffFunctionLinfty>();
+ test<3, Functions::CutOffFunctionLinfty>();
+
+ test<1, Functions::CutOffFunctionW1>();
+ test<2, Functions::CutOffFunctionW1>();
+ test<3, Functions::CutOffFunctionW1>();
+
+ test<1, Functions::CutOffFunctionCinfty>();
+ test<2, Functions::CutOffFunctionCinfty>();
+ test<3, Functions::CutOffFunctionCinfty>();
+}
--- /dev/null
+
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.17474
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00002
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 0.999735
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 0.999471
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00322
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00484
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00779
projects / dealii.git / commitdiff