Method AGQIbeginGaussLoop

SRC/element/shell/ASDShellQ4.cpp:2193–2258

}

void ASDShellQ4::AGQIbeginGaussLoop(const ASDShellQ4LocalCoordinateSystem& reference_cs)
{
    // set to zero vectors and matrices for the static condensation
    // of internal DOFs before proceeding with gauss integration
    m_eas->KQU.Zero();
    m_eas->KUQ.Zero();
    m_eas->KQQ_inv.Zero();
    m_eas->Q_residual.Zero();

    // Some matrices/vectors
    auto& N = ASDShellQ4Globals::instance().N;
    auto& dN = ASDShellQ4Globals::instance().dN;

    // Jacobian
    auto& jac = ASDShellQ4Globals::instance().jac;

    // this one is already computed in the caller method
    auto& agq = ASDShellQ4Globals::instance().agq;

    // The AGQI uses incompatible modes and only the weak patch test is passed.
    // Here we computed the mean Bq matrix for the enhanced strains. It will be subtracted
    // from the gauss-wise Bq matrices to make this element pass the strict patch test:
    // BQ = BQ - 1/A*int{BQ*dV}
    auto& BQ_mean = ASDShellQ4Globals::instance().BQ_mean;
    BQ_mean.Zero();
    double Atot = 0.0;

    // Gauss loop
    std::array<double, 4> L;
    for (int igauss = 0; igauss < 4; igauss++)
    {
        // Current integration point data
        double xi = XI[igauss];
        double eta = ETA[igauss];
        double w = WTS[igauss];
        shapeFunctions(xi, eta, N);
        shapeFunctionsNaturalDerivatives(xi, eta, dN);
        jac.calculate(reference_cs, dN);
        double dA = w * jac.detJ;
        Atot += dA;

        // area coordinates of the gauss point (Eq 7)
        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));

        // strain matrix for internal dofs
        for (int i = 0; i < 2; i++)
        { 
            int j = i + 1; if (j > 3) j = 0;
            int k = j + 1; if (k > 3) k = 0;
            double NQX = (agq.b[i] * L[k] + agq.b[k] * L[i]) / agq.A / 2.0;
            double NQY = (agq.c[i] * L[k] + agq.c[k] * L[i]) / agq.A / 2.0;
            int index1 = i * 2;
            int index2 = index1 + 1;
            BQ_mean(0, index1) += NQX * dA;
            BQ_mean(1, index2) += NQY * dA;