Function computeBdrilling

SRC/element/shell/ASDShellQ4.cpp:422–494

computes only the drilling B matrix

    // computes only the drilling B matrix
    void computeBdrilling(
        const ASDShellQ4LocalCoordinateSystem& LCS,
        double xi, double eta,
        const JacobianOperator& Jac,
        const AGQIParams& agq,
        const Vector& N,
        const Matrix& dN,
        Vector& Bd,
        bool use_eas
    )
    {
        // cartesian derivatives of standard shape function
        auto& dNdX = ASDShellQ4Globals::instance().dNdX;
        dNdX.addMatrixProduct(0.0, dN, Jac.invJ, 1.0);

        // initialize
        Bd.Zero();

        // AGQI proc ***********************************************************************************************

        if (use_eas) {

            // area coordinates of the gauss point (Eq 7)
            std::array<double, 4> L;
            L[0] = 0.25 * (1.0 - xi) * (agq.g[1] * (1.0 - eta) + agq.g[2] * (1.0 + eta));
            L[1] = 0.25 * (1.0 - eta) * (agq.g[3] * (1.0 - xi) + agq.g[2] * (1.0 + xi));
            L[2] = 0.25 * (1.0 + xi) * (agq.g[0] * (1.0 - eta) + agq.g[3] * (1.0 + eta));
            L[3] = 0.25 * (1.0 + eta) * (agq.g[0] * (1.0 - xi) + agq.g[1] * (1.0 + xi));

            // computed modified shape function gradients for the strain matrix for external dofs
            // using the area coordinate method as written in the reference paper of the AGQI element
            static constexpr std::array<double, 4> NXI = { -1.0, 1.0, 1.0, -1.0 };
            static constexpr std::array<double, 4> NETA = { -1.0, -1.0, 1.0, 1.0 };
            for (int i = 0; i < 4; i++)
            {
                int j = i + 1; if (j > 3) j = 0;
                int k = j + 1; if (k > 3) k = 0;
                double SX = 0.0;
                double SY = 0.0;
                for (int ii = 0; ii < 4; ii++)
                {
                    int jj = ii + 1; if (jj > 3) jj = 0;
                    int kk = jj + 1; if (kk > 3) kk = 0;
                    int mm = kk + 1; if (mm > 3) mm = 0;
                    SX = SX + agq.b[ii] * NXI[ii] * NETA[ii] * (3.0 * (L[jj] - L[mm]) + (agq.g[jj] - agq.g[kk]));
                    SY = SY + agq.c[ii] * NXI[ii] * NETA[ii] * (3.0 * (L[jj] - L[mm]) + (agq.g[jj] - agq.g[kk]));
                }
                dNdX(i, 0) = (agq.b[i] + agq.b[j]) / agq.A / 2.0 +
                    NXI[i] * NETA[i] * agq.g[k] * SX / 2.0 / agq.A / (1.0 + agq.g[0] * agq.g[2] + agq.g[1] * agq.g[3]);
                dNdX(i, 1) = (agq.c[i] + agq.c[j]) / agq.A / 2.0 +
                    NXI[i] * NETA[i] * agq.g[k] * SY / 2.0 / agq.A / (1.0 + agq.g[0] * agq.g[2] + agq.g[1] * agq.g[3]);
            }
        }

        // We use the drilling penalty as defined by hughes and brezzi,
        // where we link the independent rotation to the skew symmetric part of the in-plane displacement field.

        Bd(0) = -0.5 * dNdX(0, 1);