VolumeForce.hpp

Go to the documentation of this file.

#ifndef PHYTH_VOLUME_FORCE_HPP
#define PHYTH_VOLUME_FORCE_HPP
 
#include "ForceSource.hpp"
#include "Phyth/Elements/Mechanics/IsotropicMaterial.hpp"
#include "Phyth/Physical/Constants.hpp"
 
namespace Phyth::Mechanics {
    class VolumeForce : public ForceSource {
    public:
        explicit VolumeForce(
            const Vector3<Quantity<Meter> > &center,
            const Quantity<Meter> width,
            const Quantity<Meter> height,
            const Quantity<Meter> depth,
            const Vector3<Quantity<NewtonPerMeterCubed> > &force_density,
            const IsotropicMaterial &material
        ) noexcept
            : center_(center),
              width_(width),
              height_(height),
              depth_(depth),
              force_density_(force_density),
              material_(material) {
        }
 
        [[nodiscard]] Matrix3x3<Quantity<Pascal> > GetStressAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            return NumericalIntegrateStress(point, 16, 16, 16);
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetForceAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            return {};
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetTotalForce() const override {
            return force_density_ * width_ * height_ * depth_;
        }
 
        [[nodiscard]] Vector3<Quantity<NewtonPerMeterCubed> > GetBodyForceDensityAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            const auto half_w = width_ / 2.0;
            const auto half_h = height_ / 2.0;
            const auto half_d = depth_ / 2.0;
 
            if (const auto local = point - center_; Utils::abs(local.x) <= half_w &&
                                                    Utils::abs(local.y) <= half_h &&
                                                    Utils::abs(local.z) <= half_d) {
                return force_density_;
            }
            return {};
        }
 
        [[nodiscard]] Quantity<Meter> GetMinimumDistanceTo(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            const auto half_w = width_ / 2.0;
            const auto half_h = height_ / 2.0;
            const auto half_d = depth_ / 2.0;
 
            const auto local = point - center_;
            const auto dx = std::clamp(local.x, -half_w, half_w);
            const auto dy = std::clamp(local.y, -half_h, half_h);
            const auto dz = std::clamp(local.z, -half_d, half_d);
 
            const auto closest = center_ + Vector3<Quantity<Meter>>(dx, dy, dz);
            return (point - closest).Length();
        }
 
        [[nodiscard]] Vector3<Quantity<NewtonPerMeterCubed> > GetForceDensity() const {
            return force_density_;
        }
 
        void SetForceDensity(const Vector3<Quantity<NewtonPerMeterCubed> > &force_density) {
            force_density_ = force_density;
        }
 
    private:
        Vector3<Quantity<Meter> > center_;
        Quantity<Meter> width_, height_, depth_;
        Vector3<Quantity<NewtonPerMeterCubed> > force_density_;
        IsotropicMaterial material_;
 
        [[nodiscard]] Matrix3x3<Quantity<Pascal> > NumericalIntegrateStress(
            const Vector3<Quantity<Meter> > &point,
            const int nx,
            const int ny,
            const int nz
        ) const {
            Matrix3x3<Quantity<Pascal> > stress;
 
            const auto poisson_ratio = material_.GetPoissonRatio();
            const auto prefactor = -1.0 / (8.0 * Consts::pi * (1.0 - poisson_ratio));
 
            const auto dx = width_ / nx;
            const auto dy = height_ / ny;
            const auto dz = depth_ / nz;
            const auto dV = dx * dy * dz;
            const auto dF = force_density_ * dV;
 
            const auto half_w = width_ / 2.0;
            const auto half_h = height_ / 2.0;
            const auto half_d = depth_ / 2.0;
 
            for (int i = 0; i < nx; ++i) {
                for (int j = 0; j < ny; ++j) {
                    for (int k = 0; k < nz; ++k) {
                        const auto x = (i + 0.5) / nx * width_ - half_w;
                        const auto y = (j + 0.5) / ny * height_ - half_h;
                        const auto z = (k + 0.5) / nz * depth_ - half_d;
 
                        const auto pos = center_ + Vector3<Quantity<Meter>>(x, y, z);
                        const auto R = point - pos;
                        const auto R_mag = R.Length();
 
                        if (R_mag < Quantity<Meter>(Config::epsilon)) {
                            continue;
                        }
 
                        const auto R3 = R_mag * R_mag * R_mag;
                        const auto R5 = R3 * R_mag * R_mag;
 
                        const auto one_minus_2nu = 1.0 - 2.0 * poisson_ratio;
                        const auto dF_dot_R = dF.Dot(R);
 
                        for (int m = 0; m < 3; ++m) {
                            for (int n = 0; n < 3; ++n) {
                                const auto term1 = one_minus_2nu * (
                                                       dF[m] * R[n] + dF[n] * R[m]
                                                       - (m == n ? 1.0 : 0.0) * dF_dot_R
                                                   ) / R3;
 
                                const auto term2 = 3.0 * dF_dot_R * R[m] * R[n] / R5;
 
                                stress(m, n) += prefactor * (term1 - term2);
                            }
                        }
                    }
                }
            }
 
            return stress;
        }
    };
}
 
#endif // PHYTH_VOLUME_FORCE_HPP