SurfaceForce.hpp

Go to the documentation of this file.

#ifndef PHYTH_SURFACE_FORCE_HPP
#define PHYTH_SURFACE_FORCE_HPP
 
#include "ForceSource.hpp"
#include "Phyth/Elements/Mechanics/IsotropicMaterial.hpp"
#include "Phyth/Physical/Constants.hpp"
 
#include <cassert>
 
namespace Phyth::Mechanics {
    class SurfaceForce : public ForceSource {
    public:
        explicit SurfaceForce(
            const Vector3<Quantity<Meter> > &center,
            const Vector3<Scalar> &u_axis,
            const Vector3<Scalar> &v_axis,
            const Quantity<Meter> width,
            const Quantity<Meter> height,
            const Vector3<Quantity<Pascal> > &traction,
            const IsotropicMaterial &material
        ) noexcept
            : center_(center),
              u_axis_(u_axis.Normalized()),
              v_axis_(v_axis.Normalized()),
              width_(width),
              height_(height),
              traction_(traction),
              material_(material) {
            assert(Utils::abs(u_axis_.Dot(v_axis_)) < Config::epsilon);
        }
 
        [[nodiscard]] Matrix3x3<Quantity<Pascal> > GetStressAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            return NumericalIntegrateStress(point, 32, 32);
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetForceAt(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            return {};
        }
 
        [[nodiscard]] Vector3<Quantity<Newton> > GetTotalForce() const override {
            return traction_ * width_ * height_;
        }
 
        [[nodiscard]] Quantity<Meter> GetMinimumDistanceTo(
            const Vector3<Quantity<Meter> > &point
        ) const override {
            const auto local = point - center_;
            auto u_proj = local.Dot(u_axis_);
            auto v_proj = local.Dot(v_axis_);
 
            u_proj = std::clamp(u_proj, -width_ / 2, width_ / 2);
            v_proj = std::clamp(v_proj, -height_ / 2, height_ / 2);
 
            const auto closest = center_ + u_proj * u_axis_ + v_proj * v_axis_;
            return (point - closest).Length();
        }
 
        [[nodiscard]] Vector3<Quantity<Pascal> > GetTraction() const {
            return traction_;
        }
 
        void SetTraction(const Vector3<Quantity<Pascal> > &traction) {
            traction_ = traction;
        }
 
    private:
        Vector3<Quantity<Meter> > center_;
        Vector3<Scalar> u_axis_, v_axis_;
        Quantity<Meter> width_, height_;
        Vector3<Quantity<Pascal> > traction_;
        IsotropicMaterial material_;
 
        [[nodiscard]] Matrix3x3<Quantity<Pascal> > NumericalIntegrateStress(
            const Vector3<Quantity<Meter> > &point,
            const int nu,
            const int nv
        ) 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 du = width_ / nu;
            const auto dv = height_ / nv;
            const auto dA = du * dv;
            const auto dF = traction_ * dA;
 
            for (int i = 0; i < nu; ++i) {
                for (int j = 0; j < nv; ++j) {
                    const auto u = (i + 0.5) / nu * width_ - width_ / 2.0;
                    const auto v = (j + 0.5) / nv * height_ - height_ / 2.0;
 
                    const auto pos = center_ + u * u_axis_ + v * v_axis_;
                    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_SURFACE_FORCE_HPP