PointForce.hpp

Go to the documentation of this file.

#ifndef PHYTH_POINT_FORCE_HPP
#define PHYTH_POINT_FORCE_HPP
 
#include "ForceSource.hpp"
#include "Phyth/Elements/Mechanics/IsotropicMaterial.hpp"
#include "Phyth/Physical/Constants.hpp"
 
namespace Phyth::Mechanics {
    class PointForce : public ForceSource {
    public:
        explicit PointForce(const Vector3<Quantity<Meter> > &position,
                            const Vector3<Quantity<Newton> > &force,
                            const IsotropicMaterial &material) noexcept
            : position_(position), force_(force), material_(material) {
        }
 
        [[nodiscard]] Matrix3x3<Quantity<Pascal> > GetStressAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            const auto r = point - position_;
            const auto r_mag = r.Length();
 
            if (r_mag < Quantity<Meter>(Config::epsilon)) {
                throw std::runtime_error("PointForce: stress diverges at force location");
            }
 
            const auto poisson_ratio = material_.GetPoissonRatio();
            const Quantity<Joule> F_dot_r = force_.Dot(r);
            const Scalar kelvin_prefactor = -1.0 / (8.0 * Consts::pi * (1.0 - poisson_ratio));
 
            const auto r3 = r_mag * r_mag * r_mag;
            const auto r5 = r3 * r_mag * r_mag;
 
            const Scalar one_minus_2nu = 1.0 - 2.0 * poisson_ratio;
            const Scalar diagonal_coeff = 2.0 * one_minus_2nu;
            const Quantity<Joule> dipole_coeff = 3.0 * F_dot_r;
 
            Matrix3x3<Quantity<Pascal> > stress;
 
            // Diagonal components
            stress(0, 0) = kelvin_prefactor * (diagonal_coeff * force_.x * r.x / r3 - one_minus_2nu * F_dot_r / r3
                                               -
                                               dipole_coeff * r.x * r.x / r5);
            stress(1, 1) = kelvin_prefactor * (diagonal_coeff * force_.y * r.y / r3 - one_minus_2nu * F_dot_r / r3
                                               -
                                               dipole_coeff * r.y * r.y / r5);
            stress(2, 2) = kelvin_prefactor * (diagonal_coeff * force_.z * r.z / r3 - one_minus_2nu * F_dot_r / r3
                                               -
                                               dipole_coeff * r.z * r.z / r5);
 
            // Off-diagonal components (symmetric)
            stress(0, 1) = stress(1, 0) = kelvin_prefactor * (
                                              one_minus_2nu * (force_.x * r.y + force_.y * r.x) / r3 - dipole_coeff * r.
                                              x * r.y / r5);
            stress(0, 2) = stress(2, 0) = kelvin_prefactor * (
                                              one_minus_2nu * (force_.x * r.z + force_.z * r.x) / r3 - dipole_coeff * r.
                                              x * r.z / r5);
            stress(1, 2) = stress(2, 1) = kelvin_prefactor * (
                                              one_minus_2nu * (force_.y * r.z + force_.z * r.y) / r3 - dipole_coeff * r.
                                              y * r.z / r5);
 
            return stress;
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetForceAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            if ((point - position_).Length() < Quantity<Meter>(Config::epsilon)) {
                return force_;
            }
            return {};
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetTotalForce() const override {
            return force_;
        }
 
        [[nodiscard]] Quantity<Meter> GetMinimumDistanceTo(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            return (point - position_).Length();
        }
 
        [[nodiscard]] Vector3<Quantity<Meter> > GetPosition() const {
            return position_;
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetForce() const {
            return force_;
        }
 
    private:
        Vector3<Quantity<Meter> > position_;
        Vector3<Quantity<Newton> > force_;
        IsotropicMaterial material_;
    };
}
 
#endif // PHYTH_POINT_FORCE_HPP