Unverified Commit 66430774 authored by Florian Weik's avatar Florian Weik Committed by GitHub

Merge pull request #2496 from fweik/sine

ElectricPlaneWave external field
parents 036278ee 2f16260e
Pipeline #5485 passed with stages
in 111 minutes and 53 seconds
......@@ -19,8 +19,6 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef CONSTRAINTS_CONSTRAINT_HPP
#define CONSTRAINTS_CONSTRAINT_HPP
#include <memory>
#include "energy.hpp"
#include "particle_data.hpp"
......@@ -29,14 +27,29 @@ class Constraint {
public:
/**
* @brief Add energy contribution of this constraints to energy.
*
* Add constraint energy for particle to observable.
*
* @param[in] p The particle to add the energy for.
* @param[in] folded_pos Folded position of the particle.
* @param[in] t The time at which the energy should be calculated.
* @param[out] Observable to add the energy to.
*/
virtual void add_energy(const Particle &p, const Vector3d &folded_pos,
Observable_stat &energy) const = 0;
double t, Observable_stat &energy) const = 0;
/**
* @brief Return constraint force on particle.
* @brief Calculate the force of the constraint on a particle.
*
* Add constraint energy for particle to observable.
*
* @param[in] p The particle to calculate the force for.
* @param[in] folded_pos Folded position of the particle.
* @param[in] t The time at which the force should be calculated.
* @return The force on the particle.
*/
virtual ParticleForce force(const Particle &p,
const Vector3d &folded_pos) = 0;
virtual ParticleForce force(const Particle &p, const Vector3d &folded_pos,
double t) = 0;
/**
* @brief Check if constraints if compatible with box size.
......
......@@ -65,7 +65,7 @@ public:
const_iterator begin() const { return m_constraints.begin(); }
const_iterator end() const { return m_constraints.end(); }
void add_forces(ParticleRange &particles) const {
void add_forces(ParticleRange &particles, double t) const {
if (m_constraints.empty())
return;
......@@ -75,19 +75,20 @@ public:
auto const pos = folded_position(p);
ParticleForce force{};
for (auto const &c : *this) {
force += c->force(p, pos);
force += c->force(p, pos, t);
}
p.f += force;
}
}
void add_energy(ParticleRange &particles, Observable_stat &energy) const {
void add_energy(ParticleRange &particles, double t,
Observable_stat &energy) const {
for (auto &p : particles) {
auto const pos = folded_position(p);
for (auto const &c : *this) {
c->add_energy(p, pos, energy);
c->add_energy(p, pos, t, energy);
}
}
}
......
......@@ -37,11 +37,12 @@ public:
const Coupling &coupling() const { return impl.coupling(); }
const Field &field() const { return impl.field(); }
void add_energy(const Particle &, const Vector3d &,
void add_energy(const Particle &, const Vector3d &, double t,
Observable_stat &) const override {}
ParticleForce force(const Particle &p, Vector3d const &folded_pos) override {
return impl.force(p, folded_pos);
ParticleForce force(const Particle &p, const Vector3d &folded_pos,
double t) override {
return impl.force(p, folded_pos, t);
}
bool fits_in_box(Vector3d const &box) const override {
......
......@@ -37,13 +37,14 @@ public:
const Coupling &coupling() const { return impl.coupling(); }
const Field &field() const { return impl.field(); }
void add_energy(const Particle &p, const Vector3d &folded_pos,
void add_energy(const Particle &p, const Vector3d &folded_pos, double t,
Observable_stat &e) const override {
e.external_fields[0] += impl.energy(p, folded_pos);
e.external_fields[0] += impl.energy(p, folded_pos, t);
}
ParticleForce force(const Particle &p, Vector3d const &folded_pos) override {
return impl.force(p, folded_pos);
ParticleForce force(const Particle &p, const Vector3d &folded_pos,
double t) override {
return impl.force(p, folded_pos, t);
}
bool fits_in_box(Vector3d const &box) const override {
......
......@@ -22,7 +22,8 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
namespace Constraints {
ParticleForce HomogeneousMagneticField::force(const Particle &p,
const Vector3d &folded_pos) {
const Vector3d &folded_pos,
double t) {
#if defined(ROTATION) && defined(DIPOLES)
return {Vector3d{}, Vector3d::cross(p.calc_dip(), m_field)};
#else
......@@ -31,7 +32,7 @@ ParticleForce HomogeneousMagneticField::force(const Particle &p,
}
void HomogeneousMagneticField::add_energy(const Particle &p,
const Vector3d &folded_pos,
const Vector3d &folded_pos, double t,
Observable_stat &energy) const {
#ifdef DIPOLES
energy.dipolar[0] += -1.0 * m_field * p.calc_dip();
......
......@@ -32,10 +32,11 @@ public:
Vector3d const &H() const { return m_field; }
void add_energy(const Particle &p, const Vector3d &folded_pos,
void add_energy(const Particle &p, const Vector3d &folded_pos, double t,
Observable_stat &energy) const override;
ParticleForce force(const Particle &p, const Vector3d &folded_pos) override;
ParticleForce force(const Particle &p, const Vector3d &folded_pos,
double t) override;
bool fits_in_box(Vector3d const &box) const override { return true; }
......
......@@ -56,7 +56,8 @@ double ShapeBasedConstraint::min_dist() {
}
ParticleForce ShapeBasedConstraint::force(const Particle &p,
const Vector3d &folded_pos) {
const Vector3d &folded_pos,
double t) {
double dist = 0.;
Vector3d dist_vec, force, torque1, torque2, outer_normal_vec;
......@@ -117,7 +118,7 @@ ParticleForce ShapeBasedConstraint::force(const Particle &p,
}
void ShapeBasedConstraint::add_energy(const Particle &p,
const Vector3d &folded_pos,
const Vector3d &folded_pos, double t,
Observable_stat &energy) const {
double dist;
IA_parameters *ia_params;
......
......@@ -37,10 +37,11 @@ public:
ShapeBasedConstraint::reset_force();
}
void add_energy(const Particle &p, const Vector3d &folded_pos,
void add_energy(const Particle &p, const Vector3d &folded_pos, double t,
Observable_stat &energy) const override;
ParticleForce force(const Particle &p, const Vector3d &folded_pos) override;
ParticleForce force(const Particle &p, const Vector3d &folded_pos,
double t) override;
bool fits_in_box(Vector3d const &) const override { return true; }
......
......@@ -152,7 +152,7 @@ void energy_calc(double *result) {
calc_long_range_energies();
auto local_parts = local_cells.particles();
Constraints::constraints.add_energy(local_parts, energy);
Constraints::constraints.add_energy(local_parts, sim_time, energy);
#ifdef CUDA
copy_energy_from_GPU();
......
......@@ -36,9 +36,10 @@ public:
using Base::Base;
template <typename Particle>
Vector3d force(const Particle &p, const Vector3d &folded_pos) const {
Vector3d force(const Particle &p, const Vector3d &folded_pos,
double t) const {
using detail::make_bind_coupling;
return m_coupling(p, m_field(folded_pos));
return m_coupling(p, m_field(folded_pos, t));
}
};
} /* namespace FieldCoupling */
......
......@@ -36,14 +36,15 @@ public:
using Base::Base;
template <typename Particle>
Vector3d force(const Particle &p, const Vector3d &folded_pos) const {
Vector3d force(const Particle &p, const Vector3d &folded_pos,
double t) const {
using detail::make_bind_coupling;
return m_coupling(p, -m_field.gradient(folded_pos));
return m_coupling(p, -m_field.jacobian(folded_pos, t));
}
template <typename Particle>
double energy(const Particle &p, const Vector3d &folded_pos) const {
return m_coupling(p, m_field(folded_pos));
double energy(const Particle &p, const Vector3d &folded_pos, double t) const {
return m_coupling(p, m_field(folded_pos, t));
}
};
} /* namespace FieldCoupling */
......
......@@ -19,7 +19,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef CORE_EXTERNAL_FIELD_FIELDS_AFFINE_MAP_HPP
#define CORE_EXTERNAL_FIELD_FIELDS_AFFINE_MAP_HPP
#include "gradient_type.hpp"
#include "jacobian_type.hpp"
#include "utils/Vector.hpp"
namespace FieldCoupling {
......@@ -54,23 +54,23 @@ template <typename T> struct matrix_vector_impl<T, 1> {
template <typename T, size_t codim> class AffineMap {
public:
using value_type = typename decay_to_scalar<Vector<codim, T>>::type;
using gradient_type = detail::gradient_type<T, codim>;
using jacobian_type = detail::jacobian_type<T, codim>;
private:
gradient_type m_A;
jacobian_type m_A;
value_type m_b;
public:
AffineMap(const gradient_type &A, const value_type &b) : m_A(A), m_b(b) {}
AffineMap(const jacobian_type &A, const value_type &b) : m_A(A), m_b(b) {}
gradient_type &A() { return m_A; }
jacobian_type &A() { return m_A; }
value_type &b() { return m_b; }
value_type operator()(const Vector3d &pos) const {
value_type operator()(const Vector3d &pos, double = {}) const {
return detail::matrix_vector_impl<T, codim>{}(m_A, pos) + m_b;
}
gradient_type gradient(const Vector3d &) const { return m_A; }
jacobian_type jacobian(const Vector3d &, double = {}) const { return m_A; }
bool fits_in_box(const Vector3d &) const { return true; }
};
......
......@@ -19,7 +19,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef CORE_EXTERNAL_FIELD_FIELDS_CONSTANT_HPP
#define CORE_EXTERNAL_FIELD_FIELDS_CONSTANT_HPP
#include "gradient_type.hpp"
#include "jacobian_type.hpp"
#include "utils/Vector.hpp"
namespace FieldCoupling {
......@@ -30,7 +30,7 @@ namespace Fields {
template <typename T, size_t codim> class Constant {
public:
using value_type = typename decay_to_scalar<Vector<codim, T>>::type;
using gradient_type = detail::gradient_type<T, codim>;
using jacobian_type = detail::jacobian_type<T, codim>;
private:
value_type m_value;
......@@ -40,9 +40,9 @@ public:
value_type &value() { return m_value; }
value_type operator()(const Vector3d &) const { return m_value; }
static constexpr gradient_type gradient(const Vector3d &) {
return gradient_type{};
value_type operator()(const Vector3d &, double = {}) const { return m_value; }
static constexpr jacobian_type jacobian(const Vector3d &) {
return jacobian_type{};
}
bool fits_in_box(const Vector3d &) const { return true; }
......
......@@ -23,7 +23,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "utils/interpolation/bspline_3d_gradient.hpp"
#include "utils/math/tensor_product.hpp"
#include "gradient_type.hpp"
#include "jacobian_type.hpp"
#include "utils/Vector.hpp"
/* Turn off range checks if release build. */
......@@ -59,7 +59,7 @@ void deep_copy(boost::multi_array<T, 3> &dst,
template <typename T, size_t codim> class Interpolated {
public:
using value_type = typename decay_to_scalar<Vector<codim, T>>::type;
using gradient_type = detail::gradient_type<T, codim>;
using jacobian_type = detail::jacobian_type<T, codim>;
using storage_type = boost::multi_array<value_type, 3>;
private:
......@@ -98,7 +98,7 @@ public:
/*
* @brief Evaluate f at pos with the field value as argument.
*/
value_type operator()(const Vector3d &pos) const {
value_type operator()(const Vector3d &pos, double = {}) const {
using Utils::Interpolation::bspline_3d_accumulate;
return bspline_3d_accumulate<2>(
pos,
......@@ -107,14 +107,14 @@ public:
}
/*
* @brief Evaluate f at pos with the gradient field value as argument.
* @brief Evaluate f at pos with the jacobian field value as argument.
*/
gradient_type gradient(const Vector3d &pos) const {
jacobian_type jacobian(const Vector3d &pos, double = {}) const {
using Utils::Interpolation::bspline_3d_gradient_accumulate;
return bspline_3d_gradient_accumulate<2>(
pos,
[this](const std::array<int, 3> &ind) { return m_global_field(ind); },
m_grid_spacing, m_origin, gradient_type{});
m_grid_spacing, m_origin, jacobian_type{});
}
bool fits_in_box(const Vector3d &box) const {
......
/*
Copyright (C) 2010-2018 The ESPResSo project
This file is part of ESPResSo.
ESPResSo is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ESPResSo is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef CORE_EXTERNAL_FIELD_FIELDS_PLAIN_WAVE_HPP
#define CORE_EXTERNAL_FIELD_FIELDS_PLAIN_WAVE_HPP
#include "jacobian_type.hpp"
#include "utils/Vector.hpp"
namespace FieldCoupling {
namespace Fields {
/**
* @brief A plane wave.
*
* A time dependent plane wave, with a certain (vector-valued)
* amplitude, wave vector frequency and phase.
*
* See @url https://en.wikipedia.org/wiki/Plane_wave
*/
template <typename T, size_t codim> class PlaneWave {
public:
using value_type = typename decay_to_scalar<Vector<codim, T>>::type;
using jacobian_type = detail::jacobian_type<T, codim>;
private:
value_type m_amplitude;
value_type m_k;
T m_omega;
T m_phase;
public:
PlaneWave(const value_type amplitude, const value_type &k, T omega, T phase)
: m_amplitude(amplitude), m_k(k), m_omega(omega), m_phase(phase) {}
value_type &amplitude() { return m_amplitude; }
value_type &k() { return m_k; }
T &omega() { return m_omega; }
T &phase() { return m_phase; }
/**
* brief Evaluate field.
*
* @param x Where?
* @param t When?
* @return Value of the field at point x and time t.
*/
value_type operator()(const Vector3d &x, T t = 0.) const {
return m_amplitude * sin(m_k * x - m_omega * t + m_phase);
}
/**
* brief Evaluate the Jacobian matrix of the field.
*
* See @url https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant
* In the special case of a scalar field, the Jacobian is the gradient of
* the field.
*
* @param x Where?
* @param t When?
* @return Jacobian matrix
*/
jacobian_type jacobian(const Vector3d &x, T t = 0.) const {
using Utils::tensor_product;
return tensor_product(m_amplitude, m_k) *
cos(m_k * x - m_omega * t + m_phase);
}
bool fits_in_box(const Vector3d &) const { return true; }
};
} // namespace Fields
} // namespace FieldCoupling
#endif
......@@ -24,23 +24,23 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
namespace FieldCoupling {
namespace Fields {
namespace detail {
template <class T, size_t codim> struct gradient_type_impl {
template <class T, size_t codim> struct jacobian_type_impl {
using type = Vector<codim, Vector<3, T>>;
};
template <class T> struct gradient_type_impl<T, 1> {
template <class T> struct jacobian_type_impl<T, 1> {
using type = Vector<3, T>;
};
/**
* @brief Deduce type for gradient from codim.
* @brief Deduce type for jacobian from codim.
*
* Small helper that returns Vector3d if codim = 1,
* and Vector<codim, Vector<3, T>> otherwise to avoid
* using Vectors of size one, where scalars would do.
*/
template <class T, size_t codim>
using gradient_type = typename gradient_type_impl<T, codim>::type;
using jacobian_type = typename jacobian_type_impl<T, codim>::type;
} // namespace detail
} // namespace Fields
} // namespace FieldCoupling
......
......@@ -140,7 +140,7 @@ void force_calc() {
}
}
auto local_parts = local_cells.particles();
Constraints::constraints.add_forces(local_parts);
Constraints::constraints.add_forces(local_parts, sim_time);
#ifdef OIF_GLOBAL_FORCES
if (max_oif_objects) {
......
......@@ -668,6 +668,9 @@ void on_parameter_change(int field) {
/* Thermalized distance bonds needs ghost velocities */
on_ghost_flags_change();
break;
case FIELD_SIMTIME:
recalc_forces = 1;
break;
}
}
......
......@@ -320,11 +320,17 @@ void integrate_vv(int n_steps, int reuse_forces) {
if (integ_switch == INTEG_METHOD_NPT_ISO) {
propagate_vel();
propagate_pos();
/* Propagate time: t = t+dt */
sim_time += time_step;
} else if (integ_switch == INTEG_METHOD_STEEPEST_DESCENT) {
if (steepest_descent_step())
break;
} else {
propagate_vel_pos();
/* Propagate time: t = t+dt */
sim_time += time_step;
}
#ifdef BOND_CONSTRAINT
......@@ -433,9 +439,6 @@ void integrate_vv(int n_steps, int reuse_forces) {
#endif
if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
/* Propagate time: t = t+dt */
sim_time += time_step;
#ifdef COLLISION_DETECTION
handle_collisions();
#endif
......
......@@ -24,7 +24,8 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "field_coupling/fields/AffineMap.hpp"
#include "field_coupling/fields/Constant.hpp"
#include "field_coupling/fields/Interpolated.hpp"
#include "field_coupling/fields/gradient_type.hpp"
#include "field_coupling/fields/PlaneWave.hpp"
#include "field_coupling/fields/jacobian_type.hpp"
#include "utils/interpolation/bspline_3d.hpp"
#include "utils/interpolation/bspline_3d_gradient.hpp"
......@@ -33,12 +34,12 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
using namespace FieldCoupling::Fields;
BOOST_AUTO_TEST_CASE(gradient_type_test) {
using FieldCoupling::Fields::detail::gradient_type;
BOOST_AUTO_TEST_CASE(jacobian_type_test) {
using FieldCoupling::Fields::detail::jacobian_type;
using std::is_same;
static_assert(is_same<gradient_type<double, 1>, Vector3d>::value, "");
static_assert(is_same<gradient_type<double, 2>, Vector<2, Vector3d>>::value,
static_assert(is_same<jacobian_type<double, 1>, Vector3d>::value, "");
static_assert(is_same<jacobian_type<double, 2>, Vector<2, Vector3d>>::value,
"");
}
......@@ -48,7 +49,7 @@ BOOST_AUTO_TEST_CASE(constant_scalar_field) {
/* Types */
{
static_assert(std::is_same<Field::value_type, double>::value, "");
static_assert(std::is_same<Field::gradient_type, Vector3d>::value, "");
static_assert(std::is_same<Field::jacobian_type, Vector3d>::value, "");
}
/* ctor */
......@@ -80,7 +81,7 @@ BOOST_AUTO_TEST_CASE(constant_scalar_field) {
{
Field field(5.);
BOOST_CHECK(Vector3d{} == field.gradient({1., 2., 3.}));
BOOST_CHECK(Vector3d{} == field.jacobian({1., 2., 3.}));
}
}
......@@ -91,7 +92,7 @@ BOOST_AUTO_TEST_CASE(constant_vector_field) {
{
static_assert(std::is_same<Field::value_type, Vector2d>::value, "");
static_assert(
std::is_same<Field::gradient_type, Vector<2, Vector3d>>::value, "");
std::is_same<Field::jacobian_type, Vector<2, Vector3d>>::value, "");
}
/* ctor */
......@@ -123,10 +124,10 @@ BOOST_AUTO_TEST_CASE(constant_vector_field) {
/* Gradient */
{
const auto zero = Field::gradient_type{};
const auto zero = Field::jacobian_type{};
Field field({5., 6.});
BOOST_CHECK(zero == field.gradient({1., 2., 3.}));
BOOST_CHECK(zero == field.jacobian({1., 2., 3.}));
}
}
......@@ -137,7 +138,7 @@ BOOST_AUTO_TEST_CASE(affine_scalar_field) {
{
static_assert(std::is_same<Field::value_type, double>::value, "");
static_assert(
std::is_same<Field::gradient_type, Vector<3, Field::value_type>>::value,
std::is_same<Field::jacobian_type, Vector<3, Field::value_type>>::value,
"");
}
......@@ -181,7 +182,7 @@ BOOST_AUTO_TEST_CASE(affine_scalar_field) {
const double b = 4.;
Field field(A, b);
BOOST_CHECK(A == field.gradient({1., 2., 3.}));
BOOST_CHECK(A == field.jacobian({1., 2., 3.}));
}
}
......@@ -195,7 +196,7 @@ BOOST_AUTO_TEST_CASE(affine_vector_field) {
{
static_assert(std::is_same<Field::value_type, Vector2d>::value, "");
static_assert(
std::is_same<Field::gradient_type, Matrix<2, 3, double>>::value, "");
std::is_same<Field::jacobian_type, Matrix<2, 3, double>>::value, "");
}
/* Field value unshifted */
......@@ -218,7 +219,7 @@ BOOST_AUTO_TEST_CASE(affine_vector_field) {
const Vector2d b = {7., 8.};
Field field(A, b);
BOOST_CHECK(A == field.gradient({1., 2., 3.}));
BOOST_CHECK(A == field.jacobian({1., 2., 3.}));
}
}
......@@ -230,7 +231,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) {
/* Types */
{
static_assert(std::is_same<Field::value_type, double>::value, "");
static_assert(std::is_same<Field::gradient_type, Vector<3, double>>::value,
static_assert(std::is_same<Field::jacobian_type, Vector<3, double>>::value,
"");
}
......@@ -276,7 +277,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) {
std::numeric_limits<double>::epsilon());
}
/* gradient value */
/* jacobian value */
{
using Utils::Interpolation::bspline_3d_gradient_accumulate;
......@@ -293,7 +294,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) {
auto const p = Vector3d{-.4, 3.14, 0.1};
auto const field_value = field.gradient(p);
auto const field_value = field.jacobian(p);
auto const interpolated_value = bspline_3d_gradient_accumulate<2>(
p, [&data](const std::array<int, 3> &ind) { return data(ind); },
......@@ -311,7 +312,7 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) {
{
static_assert(std::is_same<Field::value_type, Vector2d>::value, "");
static_assert(
std::is_same<Field::gradient_type, Vector<2, Vector3d>>::value, "");
std::is_same<Field::jacobian_type, Vector<2, Vector3d>>::value, "");
}
/* field value */
......@@ -350,7 +351,7 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) {
std::numeric_limits<double>::epsilon());
}
/* gradient value */
/* jacobian value */
{
using Utils::Interpolation::bspline_3d_gradient_accumulate;
......@@ -376,11 +377,11 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) {
auto const p = Vector3d{-.4, 3.14, 0.1};
auto const field_value = field.gradient(p);
auto const field_value = field.jacobian(p);
auto const interpolated_value = bspline_3d_gradient_accumulate<2>(
p, [&data](const std::array<int, 3> &ind) { return data(ind); },
grid_spacing, origin, Field::gradient_type{});
grid_spacing, origin, Field::jacobian_type{});
BOOST_CHECK_SMALL((interpolated_value[0] - field_value[0]).norm(),
std::numeric_limits<double>::epsilon());
......
......@@ -103,7 +103,7 @@ BOOST_AUTO_TEST_CASE(FieldBase_test) {
}
struct DummyVectorField {
Vector3d operator()(const Vector3d &x) const { return 2. * x; }
Vector3d operator()(const Vector3d &x, double t) const { return t * x; }
};
BOOST_AUTO_TEST_CASE(ForceField_test) {
......@@ -112,24 +112,23 @@ BOOST_AUTO_TEST_CASE(ForceField_test) {
const Vector3d x{1., 2., 3.};
const int p = 5;
BOOST_CHECK((2. * x) == ff.force(5, x));
BOOST_CHECK((9. * x) == ff.force(5, x, 9.));
BOOST_CHECK(1 == ff.coupling().count);
}
struct DummyScalarField {
double operator()(const Vector3d &x) const { return 2. * x.norm(); }
Vector3d gradient(const Vector3d &x) const { return 3. * x; }
double operator()(const Vector3d &x, double t) const { return t * x.norm(); }
Vector3d jacobian(const Vector3d &x, double = {}) const { return 3. * x; }
};
BOOST_AUTO_TEST_CASE(PotentialField_test) {
auto pf = PotentialField<Id<true>, DummyScalarField>(Id<true>{},
DummyScalarField{});
const Vector3d x{1., 2., 3.};
const int p = 5;
BOOST_CHECK((2. * x.norm()) == pf.energy(5, x));
BOOST_CHECK((2. * x.norm()) == pf.energy(5, x, 2.));
BOOST_CHECK(1 == pf.coupling().count);
BOOST_CHECK(-(3. * x) == pf.force(5, x));