BowModel.hpp 5.17 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
#pragma once
#include "InputData.hpp"
#include "OutputData.hpp"
#include "DiscreteLimb.hpp"

#include "../fem/System.hpp"
#include "../fem/elements/BeamElement.hpp"
#include "../fem/elements/BarElement.hpp"
#include "../numerics/SecantMethod.hpp"

#include <QtCore>

class BowModel
{
public:
    static OutputData simulate(const InputData& input, bool statics, bool dynamics)
    {
        OutputData output;
        BowModel model(input, output);

        // Todo: Omit detailed static simulation if statics == false
        model.simulate_setup();
        model.simulate_statics();

        if(dynamics)
            model.simulate_dynamics();

        return output;
    }

private:
    const InputData& input;
    OutputData& output;
    DiscreteLimb limb;

    System system;
    std::vector<BeamElement> elements_limb;
    std::vector<BarElement> elements_string;
    std::vector<Node> nodes_limb;
    std::vector<Node> nodes_string;

    BowModel(const InputData& input, OutputData& output)
        : input(input),
          output(output),
          limb(input)
    {

    }

    void simulate_setup()
    {
        size_t n = input.settings.n_elements_limb;
        size_t k = input.settings.n_elements_string;

        // Create limb nodes
        for(size_t i = 0; i < n+1; ++i)
        {
            bool active = (i != 0);
            Node node = system.create_node({{limb.x[i], limb.y[i], limb.phi[i]}}, {{active, active, active}});
            nodes_limb.push_back(node);
        }

        // Create limb elements
        for(size_t i = 0; i < n; ++i)
        {
            double rhoA = 0.5*(limb.rhoA[i] + limb.rhoA[i+1]);
            double L = system.get_distance(nodes_limb[i], nodes_limb[i+1]);

            double Cee = 0.5*(limb.Cee[i] + limb.Cee[i+1]);
            double Ckk = 0.5*(limb.Ckk[i] + limb.Ckk[i+1]);
            double Cek = 0.5*(limb.Cek[i] + limb.Cek[i+1]);

            // Todo: Document this
            double phi = system.get_angle(nodes_limb[i], nodes_limb[i+1]);
            double phi0 = phi - system.get_u()(nodes_limb[i].phi);
            double phi1 = phi - system.get_u()(nodes_limb[i+1].phi);

            BeamElement element(nodes_limb[i], nodes_limb[i+1], rhoA, L);
            element.set_reference_angles(phi0, phi1);
            element.set_stiffness(Cee, Ckk, Cek);
            elements_limb.push_back(element);
        }

        // Limb tip
        double xt = system.get_u()(nodes_limb.back().x);
        double yt = system.get_u()(nodes_limb.back().y);

        // String center at brace height
        double xc = input.operation.brace_height;
        double yc = 0.0;

        // Create string nodes
        for(size_t i = 0; i < k; ++i)
        {
            double p = double(i)/double(k);
            double x = xc*(1.0 - p) + xt*p;
            double y = yc*(1.0 - p) + yt*p;

            Node node = system.create_node({{x, y, 0.0}}, {{true, (i != 0), false}});
            nodes_string.push_back(node);
        }

        nodes_string.push_back(nodes_limb.back());

        // Create string elements
        double rhoA = double(input.string.n_strands)*input.string.strand_density;
        double EA = double(input.string.n_strands)*input.string.strand_stiffness;

        for(size_t i = 0; i < k; ++i)
        {
            BarElement element(nodes_string[i], nodes_string[i+1], 0.0, EA, 0.0, rhoA); // Element lengths are reset later when string length is determined
            elements_string.push_back(element);
        }

        // Todo: Add elements inside loops, prevent iterator invalidation by using unique_ptr or shared_ptr
        // or maybe even let system own elements. Look at boost pointer container library.
        system.add_elements(elements_limb);
        system.add_elements(elements_string);

        // Takes a string length, iterates to equilibrium with the constraint of the brace height
        // and returns the angle of the string center
        auto try_string_length = [&](double string_length)
        {
            double L = 0.5*string_length/double(k);
            for(auto& element: elements_string)
            {
                element.set_length(L);
            }

            Dof dof = nodes_string[0].x;
            system.solve_statics_dc(dof, xc, 1, [](){});

            return system.get_angle(nodes_string[0], nodes_string[1]) - M_PI/2;
        };

        // Todo: Perhaps limit the step size of the root finding algorithm ti increase robustness.
        double string_length = 2.0*std::hypot(xc - xt, yc - yt);
        string_length = secant_method(try_string_length, 0.95*string_length, 0.9*string_length, 1e-8, 50);   // Todo: Magic numbers

        // Store setup results
        output.setup.limb = limb;
        output.setup.string_length = string_length;

        /*
        for(auto& node: nodes_limb)
        {
            qInfo() << system.get_u()(node.x) << ", " << system.get_u()(node.y);
        }

        qInfo() << "=============";

        for(auto& node: nodes_string)
        {
            qInfo() << system.get_u()(node.x) << ", " << system.get_u()(node.y);
        }
        */

    }

    void simulate_statics()
    {
        qInfo() << "Simulate Statics";
    }

    void simulate_dynamics()
    {
        qInfo() << "Simulate Dynamics";
    }

};