b2beam_cross_sections.H

//---------------------------------------------------------------------------
// b2beam_cross_sections.H --
//
//      Compute pre-defined shape beam section properties
//
// Authors: S. Merazzi, templates Th. Ludwig
//
// Copyright (c) 2016, 2017, 2023
//     SMR Engineering & Development SA
//     2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//---------------------------------------------------------------------------
 
#ifndef _B2_BEAM_CROSS_SECTIONS_
#define _B2_BEAM_CROSS_SECTIONS_
 
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <tuple>
#include <vector>
 
#include "utils/b2linear_algebra.H"
#include "utils/b2rtable.H"
 
namespace {
 
template <typename T>
T pow2(const T& v) {
    return v * v;
}
 
template <typename T>
T pow3(const T& v) {
    return v * v * v;
}
 
template <typename T>
T pow4(const T& v) {
    const T vv = v * v;
    return vv * vv;
}
 
}  // namespace
 
namespace b2000 {
 
// Thin-walled sections: A point of a "thin" line strip.
template <typename T>
class Point {
public:
    T y;  // y-coordinate with respect to (0,0)
    T z;  // z-coordinate with respect to (0,0)
    T t;  // Thickness at point
    Point(const T& y_, const T& z_, const T& t_) : y(y_), z(z_), t(t_) {}
    void print(const std::string& s) { std::cout << s << y << " " << z << " " << t << std::endl; };
};
 
// // Thin-walled sections: Line strip (contains points).
template <typename T>
class LineStrip {
public:
    std::vector<Point<T>> points;  // Points of line strip
    std::string connection;        // "open", "closed"
    // add here dict for density etc...
    LineStrip() { connection = "open"; }
    LineStrip(const std::string& s) { connection = s; }
    void append(const Point<T>& p) { points.push_back(p); };
    void print(void) {
        int k = 0;
        std::cout << "LineStrip, n. of points=" << size() << ", connection=" << connection
                  << std::endl;
        for (auto i : points) {
            auto s = "Point " + std::to_string(k) + ": ";
            i.print(s);
            k++;
        }
    };
    const int size(void) { return (int)points.size(); };
};
 
// // Thin-walled sections: LineStrips contains all LineStrip
template <typename T>
class LineStrips {
public:
    std::vector<LineStrip<T>> linestrips;
    LineStrips() {}  // All LineStrips
    void append(const LineStrip<T>& ls) { linestrips.push_back(ls); };
    void print(void) {
        std::cout << "LineStrips, n. of Line Strips=" << size() << std::endl;
        for (auto i : linestrips) { i.print(); }
    };
    void print(int ind) {
        std::cout << "LineStrip " << ind << ":" << std::endl;
        linestrips[ind].print();
    };
    const int size(void) { return (int)linestrips.size(); };
};
 
template <typename T>
class BeamCrossSection {
public:
    BeamCrossSection()
        : has(false),
          swapyz(false),
          ymin(+std::numeric_limits<T>::max()),
          ymax(std::numeric_limits<T>::lowest()),
          zmin(+std::numeric_limits<T>::max()),
          zmax(std::numeric_limits<T>::lowest()),
          yc(0),
          zc(0),
          ys(0),
          zs(0),
          A(0),
          It(0),
          Iyy(0),
          Izz(0),
          Iyz(0),
          I1(0),
          I2(0),
          K1(1),
          K2(1),
          theta(0),
          Wt(+std::numeric_limits<T>::max()),
          nseg(10) {}
 
    const std::string& get_title() const { return title; }
 
    // void add_line_strip(const LineStrip& s) { segments.push_back(s); }
 
    void add_stress_point(const T& y, const T& z) {
        int i = spoints.size1();
        spoints.resize(i + 1, 2);
        spoints(i, 0) = y;
        spoints(i, 1) = z;
    }
 
    // Computes area moments and centroid of arbitrary thin-walled section
    // composed of segments (y, z, dy, dz, angle)
    void compute_area_moments() {
        // Compute area and centroid coordinates.
        A = yc = zc = T(0);
        has = false;
        for (auto linestrip : linestrips.linestrips) {
            for (int k = 0; k != linestrip.size() - 1; ++k) {
                // mean thickness
                const T t = T(0.5) * (linestrip.points[k].t + linestrip.points[k + 1].t);
                // std::cerr << "Linestrip " << k <<": t=" << t << std::endl;
                // Compute delta-centroid
                T c[2] = {
                      linestrip.points[k + 1].y - linestrip.points[k].y,
                      linestrip.points[k + 1].z - linestrip.points[k].z};  // direction
                T ab = (std::sqrt(c[0] * c[0] + c[1] * c[1]));             // length
                // std::cerr << "ab=" << ab << std::endl;
                T yyc = T(0.5) * (linestrip.points[k + 1].y + linestrip.points[k].y);
                T zzc = T(0.5) * (linestrip.points[k + 1].z + linestrip.points[k].z);
                // std::cerr << "yyc=" << yyc << " zzc=" << zzc << std::endl;
                T darea = ab * t;
                // std::cerr << "darea=" << darea << std::endl;
                yc += yyc * darea;
                zc += zzc * darea;
                A += darea;
            }
        }
        // if (std::abs(A) < std::numeric_limits<T>::epsilon())
        //     Exception() << "Cross section area (" << A << ") is zero or "nearly zero."  << THROW;
        yc /= A;
        zc /= A;
 
        // Compute inertia moments at centroid
        ymin = +std::numeric_limits<T>::max();
        ymax = std::numeric_limits<T>::lowest();
        zmin = +std::numeric_limits<T>::max();
        zmax = std::numeric_limits<T>::lowest();
        It = Iyy = Izz = Iyz = K1 = K2 = T(0);
        T tmax = std::numeric_limits<T>::lowest();
        for (auto linestrip : linestrips.linestrips) {
            // std::cerr << "---Linestrip " << std::endl;
            for (int k = 0; k != linestrip.size() - 1; ++k) {
                // mean thickness
                const T t = T(0.5) * (linestrip.points[k].t + linestrip.points[k + 1].t);
                // std::cerr << "t=" << t << std::endl;
                // Compute delta-centroid
                T c[2] = {
                      linestrip.points[k + 1].y - linestrip.points[k].y,
                      linestrip.points[k + 1].z - linestrip.points[k].z};
                // std::cerr << "c[0]=" << c[0] << " c[1]=" << c[1] << std::endl;
                T seglength = (std::sqrt(c[0] * c[0] + c[1] * c[1]));  // length
                // std::cerr << "seglength=" << seglength << std::endl;
                T yyc = T(0.5) * (linestrip.points[k + 1].y + linestrip.points[k].y) - yc;
                T zzc = T(0.5) * (linestrip.points[k + 1].z + linestrip.points[k].z) - zc;
                // std::cerr << "yyc=" << yyc << " zzc=" << zzc << std::endl;
                T darea = seglength * t;
                // // Inertia moments of rectangle about its centroid.
                T Iyyp = seglength * pow3(t) / T(12);
                T Izzp = t * pow3(seglength) / T(12);
                // std::cerr << "Iyyp=" << Iyyp << " Izzp=" << Izzp << std::endl;
                T alfa = std::atan2(c[1], c[0]);
                // // Rotate
                T c2 = std::cos(2. * alfa);
                T s2 = std::sin(2. * alfa);
                // std::cerr << "alfa=" << alfa<< " cos=" << c2<< " sin=" << s2 << std::endl;
                T p = T(0.5) * (Iyyp + Izzp);
                T m = T(0.5) * (Iyyp - Izzp);
                T Iyypp = p + m * c2;
                T Izzpp = p - m * c2;
                T Iyzpp = -m * s2;
                // std::cerr << "Iyypp=" << Iyypp+ zzc * zzc * darea << " Izzpp=" << Izzpp + yyc *
                // yyc * darea << " Iyzpp=" << Iyzpp + yyc * zzc * darea << std::endl;
                // // Compute at centroid
                Iyy += Iyypp + zzc * zzc * darea;
                Izz += Izzpp + yyc * yyc * darea;
                Iyz += Iyzpp + yyc * zzc * darea;
                if (linestrip.connection == "open") {
                    // Torsional moment: Bredts' formula
                    It += seglength * pow3(t) * T(0.3333333333333333);
                    tmax = std::max(tmax, t);
                    Wt = It / tmax;
                }
                ymin = std::min(ymin, linestrip.points[k].y - yc);
                ymin = std::min(ymin, linestrip.points[k + 1].y - yc);
                zmin = std::min(zmin, linestrip.points[k].z - zc);
                zmin = std::min(zmin, linestrip.points[k + 1].z - zc);
 
                ymax = std::max(ymax, linestrip.points[k].y - yc);
                ymax = std::max(ymax, linestrip.points[k + 1].y - yc);
                zmax = std::max(zmax, linestrip.points[k].z - zc);
                zmax = std::max(zmax, linestrip.points[k + 1].z - zc);
            }
        }
        if (std::abs(Iyz) > std::numeric_limits<T>::epsilon() * std::abs(Iyy)
            && std::abs(Iyz) > std::numeric_limits<T>::epsilon() * std::abs(Izz)) {
            const T pi = std::acos(T(-1));
            theta = T(90) * std::atan(2. * Iyz / (-Iyy + Izz)) / pi;
        } else {
            theta = T(0);
        }
        // principal moments
        T a = T(0.5) * (Iyy + Izz);
        T b = std::sqrt(pow2(T(0.5) * (Iyy - Izz)) + pow2(Iyz));
        I1 = a + b;
        I2 = a - b;
 
        // Shear center = centroid by default (not (yet) computed
        // here. Modified in shape dependent code.
        ys = yc;
        zs = zc;
        has = true;
    }
 
    // Generates a LineSegment of a straight line from a to b with
    // thickness t. Return area of segment.
    T add_line_straight(const T ya, const T za, const T yb, const T zb, const T& t) {
        LineStrip<T> ls;
        ls.append(Point<T>(ya, za, t));
        ls.append(Point<T>(yb, zb, t));
        linestrips.append(ls);
        T c[2] = {yb - ya, zb - za};                    // direction
        T ab = (std::sqrt(c[0] * c[0] + c[1] * c[1]));  // length
        return t * std::sqrt(c[0] * c[0] + c[1] * c[1]);
    }
 
    void print(const std::string& text) {
        std::cout << "*****************************************************************************"
                  << std::endl
                  << "Cross section properties ";
        if (text.size() > 0) { std::cout << "\"" << text << "\""; }
        std::cout << std::endl
                  << "*****************************************************************************"
                  << std::endl
                  << std::endl;
        std::cout << "Shape=\"" << shape << "\"" << std::endl
                  << "length unit=\"" << length_unit << "\"" << std::endl
                  << std::setprecision(std::numeric_limits<double>::digits10 + 1) << std::scientific
                  << "yc=" << yc << " zc=" << zc << std::endl
                  << "ys=" << ys << " zs=" << zs << std::endl
                  << "A=" << A << std::endl
                  << "Iyy=" << Iyy << std::endl
                  << "Izz=" << Izz << std::endl
                  << "Iyz=" << Iyz << std::endl
                  << "It=" << It << std::endl
                  << "Wt=" << Wt << std::endl
                  << "theta=" << theta << std::endl
                  << "I1=" << I1 << std::endl
                  << "I2=" << I2 << std::endl
                  << "K1=" << K1 << std::endl
                  << "K2=" << K2 << std::endl
                  << "ymin=" << ymin << " ymax=" << ymax << std::endl
                  << "zmin=" << zmin << " zmax=" << zmax << std::endl
                  << std::endl;
 
        if (spoints.size1() > 0) {
            std::cout << "Stress points y, z:" << std::endl;
            for (int row = 0; row < spoints.size1(); ++row) {
                for (int col = 0; col < 2; ++col) { std::cout << spoints(row, col) << " "; }
                std::cout << std::endl;
            }
        }
    }
 
    void add_properties(RTable& keys) {
        keys.set("ymin", ymin);
        keys.set("ymax", ymax);
        keys.set("zmin", zmin);
        keys.set("zmax", zmax);
        keys.set("yc", yc);
        keys.set("zc", zc);
        keys.set("ys", ys);
        keys.set("zs", zs);
        keys.set("A", A);
        keys.set("It", It);
        keys.set("Iyy", Iyy);
        keys.set("Izz", Izz);
        keys.set("Iyz", Iyz);
        keys.set("I1", I1);
        keys.set("I2", I2);
        keys.set("K1", K1);
        keys.set("K2", K2);
        keys.set("theta", theta);
        keys.set("Wt", Wt);
        keys.set("length_unit", length_unit);
        if (spoints.size1() > 0) {
            std::vector<T> points;
            points.reserve(spoints.size1() * spoints.size2());
            for (int row = 0; row < spoints.size1(); ++row) {
                points.push_back(spoints(row, 0));
                points.push_back(spoints(row, 1));
            }
            keys.set("SPOINTS", points);
        }
    }
 
    void doswap() {
        std::swap(Iyy, Izz);
        std::swap(I1, I2);
        std::swap(K1, K2);
        std::swap(yc, zc);
        std::swap(ys, zs);
        std::swap(ymin, zmin);
        std::swap(ymax, zmax);
    }
 
    // Update stress y, z point coordinates with respect to centroid.
    void update_stress_points() {
        for (int row = 0; row < spoints.size1(); ++row) {
            spoints(row, 0) -= yc;
            spoints(row, 1) -= zc;
        }
    }
 
protected:
    std::string title;
    std::string shape;
    std::string geometry;  // Type of section: "open", "closed"
    LineStrips<T> linestrips;
    bool has;     // Set to true if section proerties are defined
    bool swapyz;  // Set to 1 if y and z axes parameters have to be swapped
 
public:
    T ymin;  // minmax with respect to centroid if applicable
    T ymax;
    T zmin;
    T zmax;
    T yc;  // centroid
    T zc;
    T ys;  // Shear center
    T zs;
    T A;    // Area
    T It;   // Torsion moment
    T Iyy;  // Inertia (area) moments
    T Izz;
    T Iyz;
    T I1;  // Major pricipal moment
    T I2;  // Minor pricipal moment
    T K1;
    T K2;     // Shear correction factors
    T theta;  // Angle to principal y axis
    T Wt;     // Torsional resistance moment
    std::string length_unit;
    // Stress evaluation points (y,z)
    // b2linalg::Matrix<T, b2linalg::Mrectangle> spoints;
    b2000::b2linalg::Matrix<T, b2000::b2linalg::Mrectangle> spoints;
    int nseg = 10;
};
 
// Computes bar (rectangle) section properties with respect to centroid.
template <typename T>
class BeamCrossSectionBar : public BeamCrossSection<T> {
private:
    T width_;
    T height_;
 
public:
    void compute_properties(const T& width, const T& height, const bool swap) {
        width_ = width;
        height_ = height;
        this->nseg = 0;  // Not used (analytical)
        this->shape = "Bar";
        this->geometry = "closed";
        this->swapyz = swap;
        this->A = height_ * width_;
        this->Iyy = width_ * pow3(height_) / T(12);
        this->Izz = height_ * pow3(width_) / T(12);
        this->ymax = T(0.5) * width_;
        this->ymin = -this->ymax;
        this->zmax = T(0.5) * height_;
        this->zmin = -this->zmax;
        // const T K = 10  *  (1 + nu) / (12 + 1 * nu);
        this->K1 = this->K2 = T(5) / T(6);
        // Roark
        T a = width_;
        T b = height_;
        if (width_ >= height_) {
            this->I1 = this->Izz;
            this->I2 = this->Iyy;
            this->theta = std::acos(T(-1)) * T(.5);
        } else {
            this->I1 = this->Iyy;
            this->I2 = this->Izz;
            this->theta = 0;
            std::swap(a, b);
        }
        this->It =
              a * pow3(b) * (T(1) / T(3) - T(0.21) * b / a * (T(1) - pow4(b) / (T(12) * pow4(a))));
        // From Roark and Young
        // if (dstd::swap(1 < d2)
        //    std::swap(d1, d2);
        // const T J = d1 * std::pow(d2, 3) * (1.0 / 3 - 0.21 * d2 / d1 * (1 - std::pow(d2 / d1, 4)
        // / 12));
        // self.Wt = ?? # must be calculated by pp
 
        // Add stress points
        this->add_stress_point(this->ymax, this->zmax);
        this->add_stress_point(this->ymax, this->zmin);
        this->add_stress_point(this->ymin, this->zmin);
        this->add_stress_point(this->ymin, this->zmax);
 
        if (this->swapyz) { this->doswap(); }
    }
 
    void load(const RTable& keys) {
        width_ = keys.get<T>("D1");
        height_ = keys.get<T>("D2");
        const bool s = keys.get_bool("SWAP", false);
        compute_properties(width_, height_, s);
    }
};
 
// Computes box cross section properties with respect to centroid.
template <typename T>
class BeamCrossSectionBox : public BeamCrossSection<T> {
private:
    T width_;    // D1 (width of box)
    T height_;   // D2(height of box)
    T tflange_;  // D3 (thickness of horizontal walls along y)
    T tweb_;     // D4 (thickness of vertical walls along z)
 
public:
    void compute_properties(
          const T& width, const T& height, const T& tflange, const T& tweb, const bool swap) {
        width_ = width;
        height_ = height;
        tflange_ = tflange;
        tweb_ = tweb;
        this->shape = "Box";
        this->geometry = "closed";
        this->swapyz = swap;
        this->nseg = 1;
        this->A = height_ * width_ - (height_ - T(2) * tflange) * (width_ - T(2) * tflange);
        this->Iyy = width_ * pow3(height_) / T(12)
                    - (width_ - T(2) * tweb_) * pow3(height_ - T(2) * tflange_) / T(12);
        this->Izz = height_ * pow3(width_) / T(12)
                    - (height_ - T(2) * tflange_) * pow3(width_ - T(2) * tweb_) / T(12);
        //    It=2.*tflange*tweb*(width_-tweb)**2*(h-tflange)**2/(width_*tweb+h*tflange-tweb**2-tflange**2);
        T A1 = T(2) * height_ * tweb_;
        T A2 = T(2) * width_ * tflange_;
        if (width_ >= height_) {
            this->I1 = this->Izz;
            this->I2 = this->Iyy;
            this->theta = T(.5) * std::acos(T(-1));
        } else {
            this->I1 = this->Iyy;
            this->I2 = this->Izz;
        }
        this->K1 = A2 / this->A;
        this->K2 = A1 / this->A;
        // Torsion: Zweite Bredt'sche Formel
        T ww = width_ - tweb_;
        T hh = height_ - tflange_;
        T Am = ww * hh;
        this->It = 4. * Am * Am / (2. * (ww / tflange_ + hh / tweb_));
        this->Wt = 2. * Am * tweb_;
 
        this->add_stress_point(this->ymax, this->zmax);
        this->add_stress_point(this->ymax, this->zmin);
        this->add_stress_point(this->ymin, this->zmin);
        this->add_stress_point(this->ymin, this->zmax);
 
        if (this->swapyz) { this->doswap(); }
    }
 
    void load(const RTable& keys) {
        width_ = keys.get<T>("D1");
        height_ = keys.get<T>("D2");
        tflange_ = keys.get<T>("D3");
        tweb_ = keys.get<T>("D4");
        const bool s = keys.get_bool("SWAP", false);
        compute_properties(width_, height_, tflange_, tweb_, s);
    }
};
 
// Computes C cross section properties with respect to centroid.
template <typename T>
class BeamCrossSectionC : public BeamCrossSection<T> {
private:
    T width_;    // D1 (profile flange width)
    T height_;   // D2 (profile web height)
    T tweb_;     // D3 (web thickness)
    T tflange_;  // D4 (thickness of flanges)
 
public:
    void compute_properties(
          const T& width, const T& height, const T& tweb, const T& tflange, const bool swap) {
        width_ = width;      // d1
        height_ = height;    // d2
        tweb_ = tweb;        // d3
        tflange_ = tflange;  // d4;
        this->shape = "C";
        this->geometry = "open";
        this->swapyz = swap;
        this->nseg = 1;
        T a[2];
        T b[2];
        T z1 = T(0.5) * tflange_;
        T z2 = height_ - z1;
        // Lower flange
        T A1 = this->add_line_straight(tweb_, z1, width_, z1, tflange);
        // Upper flange
        T A2 = this->add_line_straight(tweb_, z2, width_, z2, tflange_);
        // Web
        T tweb2 = T(0.5) * tweb_;
        T A3 = this->add_line_straight(tweb2, T(0), tweb2, height_, tweb_);
        this->compute_area_moments();
 
        // Update numerically (thin walled) computed values
 
        // Add shear center (thin-walled) analytically!
        this->ys = -pow2((width_ - T(0.5) * tweb_)) * tflange_
                   / (T(2.) * (width_ - T(0.5) * tweb_) * tflange_
                      + (height_ - tflange_) * tweb_ / T(3));
        this->ys -= this->yc;  // update for centroid origin
        this->zs = T(0);
 
        this->K2 = height_ * tweb_ / this->A;
        this->K1 = 2. * width_ * tflange_ / this->A;
 
        // Add stress points
        this->add_stress_point(this->ymax, this->zmax);
        this->add_stress_point(this->ymax, this->zmin);
        this->add_stress_point(this->ymin, this->zmin);
        this->add_stress_point(this->ymin, this->zmax);
 
        if (this->swapyz) { this->doswap(); }
    }
 
    void load(const RTable& keys) {
        width_ = keys.get<T>("D1");
        height_ = keys.get<T>("D2");
        tweb_ = keys.get<T>("D3");
        tflange_ = keys.get<T>("D4");
        const bool s = keys.get_bool("SWAP", false);
        compute_properties(width_, height_, tweb_, tflange_, s);
    }
};
 
// Computes I cross section properties with respect to centroid.
template <typename T>
class BeamCrossSectionI : public BeamCrossSection<T> {
private:
    T height_;
    T width1_;
    T width2_;
    T tweb_;
    T tflange1_;
    T tflange2_;
 
public:
    void compute_properties(
          const T& height, const T& width1, const T& width2, const T& tweb, const T& tflange1,
          const T& tflange2, const bool swap) {
        height_ = height;      // Total height of profile (=web height)
        width1_ = width1;      // Width of lower flange
        width2_ = width2;      // Width of upper flange
        tweb_ = tweb;          // Thickness of web
        tflange1_ = tflange1;  // Thickness of lower flange
        tflange2_ = tflange2;  // Thickness of upper flange
        this->swapyz = swap;
        this->shape = "I";
        this->geometry = "open";
        this->nseg = 1;
        // Lower flange
        T y = T(0.5) * width1_;
        T z = T(0.5) * (height_ - tflange1_);
        T A1 = this->add_line_straight(-y, -z, y, -z, tflange1_);
        // Upper flange
        y = T(0.5) * width2_;
        z = T(0.5) * (height_ - tflange2_);
        T A2 = this->add_line_straight(-y, z, y, z, tflange2_);
        // web
        T A3 = this->add_line_straight(
              T(0), -T(0.5) * height_ + tflange1_, T(0), T(0.5) * height_ - tflange2_, tweb_);
        this->compute_area_moments();
 
        this->K1 = (A1 + A2) / this->A;
        this->K2 = A3 / this->A;
 
        // Add stress points
        this->add_stress_point(this->ymax, this->zmax);
        this->add_stress_point(this->ymin, this->zmax);
        this->add_stress_point(this->ymin, this->zmin);
        this->add_stress_point(this->ymax, this->zmin);
 
        if (this->swapyz) { this->doswap(); }
    }
 
    void load(const RTable& keys) {
        height_ = keys.get<T>("D1");
        width1_ = keys.get<T>("D2");
        width2_ = keys.get<T>("D3");
        tweb_ = keys.get<T>("D4");
        tflange1_ = keys.get<T>("D5");
        tflange2_ = keys.get<T>("D6");
        const bool s = keys.get_bool("SWAP", false);
        compute_properties(height_, width1_, width2_, tweb_, tflange1_, tflange2_, s);
    }
};
 
// Computes L cross section properties with respect to centroid.
template <typename T>
class BeamCrossSectionL : public BeamCrossSection<T> {
private:
    T width_;
    T height_;
    T tflange_;
    T tweb_;
 
public:
    void compute_properties(
          const T& width, const T& height, const T& tflange, const T& tweb, const bool swap) {
        width_ = width;      // D1
        height_ = height;    // D2
        tflange_ = tflange;  // D3
        tweb_ = tweb;        // D4
        this->shape = "L";
        this->geometry = "open";
        this->swapyz = swap;
        this->nseg = 1;
        T a[2];
        T b[2];
 
        // Web
        T tweb2 = T(0.5) * tweb_;
        T A1 = this->add_line_straight(tweb2, T(0), tweb2, height_, tweb_);
        // Flange
        T tflange2 = T(0.5) * tflange_;
        T A2 = this->add_line_straight(tweb_, tflange2, width_, tflange2, tflange_);
        this->compute_area_moments();
 
        // shear center from Roark table 8.12: Thin walled ey=ez=0
        this->ys = tweb2 - this->yc;
        this->zs = tflange2 - this->zc;
 
        this->K1 = A2 / this->A;
        this->K2 = A1 / this->A;
        // Stress points in definition system: (0,0) in lower left corner of L)
        this->add_stress_point(tweb_, height_);
        this->add_stress_point(width_, T(0));
        this->add_stress_point(T(0), T(0));
        this->add_stress_point(T(0), height_);
        this->update_stress_points();  // Correct for Centroid
        if (this->swapyz) { this->doswap(); }
    }
 
    void load(const RTable& keys) {
        width_ = keys.get<T>("D1");
        height_ = keys.get<T>("D2");
        tflange_ = keys.get<T>("D3");
        tweb_ = keys.get<T>("D4");
        const bool s = keys.get_bool("SWAP", false);
        compute_properties(width_, height_, tflange_, tweb_, s);
    }
};
// Computes upside down T shape cross section properties with respect to centroid.
template <typename T>
class BeamCrossSectionT2 : public BeamCrossSection<T> {
private:
    T width_;
    T height_;
    T tflange_;
    T tweb_;
 
public:
    void compute_properties(
          const T& width, const T& height, const T& tflange, const T& tweb, const bool swap) {
        width_ = width;      // D1
        height_ = height;    // D2
        tflange_ = tflange;  // D3
        tweb_ = tweb;        // D4
        this->shape = "T2";
        this->geometry = "open";
        this->swapyz = swap;
        this->nseg = 1;
        T a[2];
        T b[2];
 
        // Web
        T tweb2 = T(0.5) * tweb_;
        T A1 = this->add_line_straight(T(0), tflange_, T(0), height_, tweb_);
        // Flange
        T w2 = T(0.5) * width_;
        T tflange2 = T(0.5) * tflange_;
        T A2 = this->add_line_straight(-w2, tflange2, w2, tflange2, tflange_);
        this->compute_area_moments();
 
        // shear center from Roark table 8.12: Thin walled ey=ez=0
        T t2 = height_ - tweb_;
        this->ys = T(0);
        this->zs = -this->zc + tflange2
                   + 0.5 * (tflange_ + t2) / (T(1) + pow3(width_) * tflange_ / (pow3(tweb_) * t2));
 
        this->K1 = A2 / this->A;
        this->K2 = A1 / this->A;
        // Stress points in definition system: (0,0) in lower left corner of L
        this->add_stress_point(tflange2, height_);
        this->add_stress_point(w2, T(0));
        this->add_stress_point(-w2, T(0));
        this->add_stress_point(-tflange2, height_);
        this->update_stress_points();  // Correct for Centroid
        if (this->swapyz) { this->doswap(); }
    }
 
    void load(const RTable& keys) {
        width_ = keys.get<T>("D1");
        height_ = keys.get<T>("D2");
        tflange_ = keys.get<T>("D3");
        tweb_ = keys.get<T>("D4");
        const bool s = keys.get_bool("SWAP", false);
        compute_properties(width_, height_, tflange_, tweb_, s);
    }
};
 
// Computes circular tube section properties with respect to centroid.
template <typename T>
class BeamCrossSectionTube : public BeamCrossSection<T> {
private:
    T r_;   // D1 outer radius
    T ri_;  // D2 inner radius
 
public:
    void compute_properties(const T& r, const T& ri, const bool swap) {
        r_ = r;
        ri_ = ri;
        if (real(r_) < std::numeric_limits<T>::epsilon()) {
            Exception() << "Tube cross section: radius (" << r_ << ") <= 0" << THROW;
        }
        if (ri_ >= r_) { Exception() << "Tube cross section: inner radius is >= radius" << THROW; }
        this->nseg = 0;  // Not used (analytical)
        this->shape = "CircularTube";
        this->geometry = "closed";
        this->swapyz = swap;
        if (ri_ < std::numeric_limits<T>::epsilon()) { ri_ = 0; }
        const T pi = T(M_PI);  // std::acos(T(-1));
        this->A = pi * (r_ * r_ - ri_ * ri_);
        this->Iyy = this->Izz = this->I1 = this->I2 = T(0.25) * pi * (pow4(r_) - pow4(ri_));
        this->ymax = this->zmax = r;
        this->ymin = this->zmin = -r;
        this->K1 = this->K2 = 0.9;
        this->It = T(2.0) * this->Iyy;
        this->Wt = this->It / r_;
        // Add stress points
        T delta = T(0.125) * pi;
        T phi = 0.0;
        for (int i = 0; i < 16; i++) {
            this->add_stress_point(r_ * std::cos(phi), r_ * std::sin(phi));
            phi += delta;
        }
        // not needed if (swapyz) doswap();
    }
 
    void load(const RTable& keys) {
        r_ = keys.get<T>("D1");
        ri_ = keys.get<T>("D2", T(0));
        compute_properties(r_, ri_, false);
    }
};
 
// Computes rod section properties with respect to centroid. Is not
// used, because properties (B2000++ MDL property option) are not
// (yet) supported by Rx.S (rod) elements. Rx.S elements take mid and
// area instead.
 
template <typename T>
class BeamCrossSectionRod : public BeamCrossSection<T> {
public:
    void compute_properties(const T& r) {
        if (real(r) < std::numeric_limits<T>::epsilon()) {
            Exception() << "Rod cross section radius (" << r << ") is <= 0" << THROW;
        }
        this->shape = "Rod";
        this->geometry = "closed";
        this->A = r * r * T(M_PI);
        // All other properties are undefined (rod element!)
    }
 
    void load(const RTable& keys) {
        const double area = keys.get<T>("D1");
        compute_properties(area);
    }
};
 
}  // namespace b2000
 
#endif  // _B2_BEAM_CROSS_SECTIONS_