Method getStiffnessMatrix

tools/Assembler/ref/Assembler3D.cpp:168–218

}

ZSparseMatrix Assembler3D::getStiffnessMatrix() {
    // Elastic modulii
    Eigen::MatrixXd& D = m_D;
    Eigen::MatrixXd& C = m_C;

    typedef Eigen::Triplet<double> T;
    std::vector<T> triplets;

    for (size_t i=0; i<m_mesh->getNbrElements(); ++i)
    {
        VectorI idx = m_mesh->getElement(i);
        assert(idx.size() == 4);

        Eigen::MatrixXd &dN = m_DN[i];

        // Small strain-displacement matrix
        //
        Eigen::MatrixXd B(6,12);
        B << dN(0,0),0.0,0.0,dN(1,0),0.0,0.0,dN(2,0),0.0,0.0,dN(3,0),0.0,0.0,
             0.0,dN(0,1),0.0,0.0,dN(1,1),0.0,0.0,dN(2,1),0.0,0.0,dN(3,1),0.0,
             0.0,0.0,dN(0,2),0.0,0.0,dN(1,2),0.0,0.0,dN(2,2),0.0,0.0,dN(3,2),

             0.5*dN(0,1),0.5*dN(0,0),0.0,
             0.5*dN(1,1),0.5*dN(1,0),0.0,
             0.5*dN(2,1),0.5*dN(2,0),0.0,
             0.5*dN(3,1),0.5*dN(3,0),0.0,

             0.5*dN(0,2),0.0,0.5*dN(0,0),
             0.5*dN(1,2),0.0,0.5*dN(1,0),
             0.5*dN(2,2),0.0,0.5*dN(2,0),
             0.5*dN(3,2),0.0,0.5*dN(3,0),

             0.0,0.5*dN(0,2),0.5*dN(0,1),
             0.0,0.5*dN(1,2),0.5*dN(1,1),
             0.0,0.5*dN(2,2),0.5*dN(2,1),
             0.0,0.5*dN(3,2),0.5*dN(3,1);


        Eigen::MatrixXd k_el = B.transpose() * D * C * B * m_mesh->getElementVolume(i);

        for (size_t j=0; j<4; ++j)
            for (size_t k=0; k<4; ++k)
                for (size_t l=0; l<3; ++l)
                    for (size_t m=0; m<3; ++m)
                        triplets.push_back(T(3*idx[j]+l, 3*idx[k]+m, k_el(3*j+l, 3*k+m)));
    }

    Eigen::SparseMatrix<double> K = Eigen::SparseMatrix<double>(3*m_mesh->getNbrNodes(),3*m_mesh->getNbrNodes());
    K.setFromTriplets(triplets.begin(), triplets.end());
    return ZSparseMatrix(K);
}

ZSparseMatrix Assembler3D::getLaplacianMatrix() {
    typedef Eigen::Triplet<double> T;