Permalink
Browse files

Merge pull request #2496 from fweik/sine

ElectricPlaneWave external field
  • Loading branch information...
fweik committed Feb 8, 2019
2 parents 036278e + 2f16260 commit 664307740c951a3211b41fdceea0a3b7addd4a77
@@ -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 @@ template <class ParticleRange, class Constraint> class Constraints {
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 @@ template <class ParticleRange, class Constraint> class Constraints {
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 @@ class ExternalField : public Constraint {
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 @@ class ExternalPotential : public Constraint {
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 @@ class HomogeneousMagneticField : public Constraint {

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 @@ class ShapeBasedConstraint : public Constraint {
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 @@ class ForceField : public detail::Base<Coupling, Field> {
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 @@ class PotentialField : public detail::Base<Coupling, Field> {
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 @@ template <typename T, size_t codim> class Constant {

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 @@ template <typename T, size_t codim> class Interpolated {
/*
* @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 @@ template <typename T, size_t codim> class Interpolated {
}

/*
* @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 {
Oops, something went wrong.

0 comments on commit 6643077

Please sign in to comment.