HOB 코드 전환 저장용

// matlab

function [Ke, Me, TLG] = MITC4(Young, Poisson, thick, density, node, Vn)
    
    nNPE = 4; % # of nodes per element
    nDPN = 5; % # of DOFs per node
    nDIM = 3; % Problem dimension
    nGP  = 2;  % # of Gauss points
    
    % Vn = Make_Vn(node);
    [V1,V2] = Director(Vn(1,:)); % norn vec 구하기 
    TLG = [V1' V2' Vn(1,:)'];
    
    Ke = Shell_Stiffness(Young, Poisson, thick, node, Vn);
    
    Me = Shell_Mass(thick, density, node, Vn);
    
    % nodewise2dofwise = reshape(1:nDPN*nNPE,nDPN,nNPE)';
    % nodewise2dofwise = nodewise2dofwise(:)';
    dofwise2nodewise = reshape(1:nDPN*nNPE,nNPE,nDPN)';
    dofwise2nodewise = dofwise2nodewise(:)';
    
    Ke = Ke(dofwise2nodewise,dofwise2nodewise);
    Me = Me(dofwise2nodewise,dofwise2nodewise);
    
    
    
        function Me = Shell_Mass(thick, density, node, Vn)
    
            Me = zeros(nDPN*nNPE);
    
            for i = 1:nGP
                for j = 1:nGP
                    for k = 1:2
    
                        % ! Gauss integration point and weight factor
                        [r, s, t, weight] = Gauss22_Point(i, j, k);
    
                        % ! Jacobian
                        Jacob = Jacobian(r, s, t, thick, node, Vn);
    
                        % ! uvw: shape function matrix for displacment
                        uvw = Displacement(r, s, t, thick, Vn);
    
                        %       Me = Me + weight * det(Jacob) * thick * transpose(uvw) * uvw * density;
                        Me = Me + weight * det(Jacob) * transpose(uvw) * uvw * density;
                    end
                end
            end
    
        end

C++

#include <Eigen/Dense>
#include <vector>

using namespace Eigen;

void MITC4(double Young, double Poisson, double thick, double density, const MatrixXd& node, const MatrixXd& Vn, MatrixXd& Ke, MatrixXd& Me, MatrixXd& TLG) {
    
    const int nNPE = 4; // # of nodes per element
    const int nDPN = 5; // # of DOFs per node
    const int nDIM = 3; // Problem dimension
    const int nGP  = 2; // # of Gauss points
    
    Vector3d V1, V2;
    std::tie(V1, V2) = Director(Vn.row(0)); // norn vec 구하기 
    TLG = (MatrixXd(3, 1) << V1, V2, Vn(0, 0)).finished();
    
    Ke = Shell_Stiffness(Young, Poisson, thick, node, Vn);
    
    Me = Shell_Mass(thick, density, node, Vn);
    
    VectorXi dofwise2nodewise = VectorXi::LinSpaced(nDPN * nNPE, 1, nDPN * nNPE);
    dofwise2nodewise = dofwise2nodewise.array() - 1; // Convert to 0-based indexing
    
    Ke = Ke(dofwise2nodewise, dofwise2nodewise);
    Me = Me(dofwise2nodewise, dofwise2nodewise);
}

MatrixXd Shell_Mass(double thick, double density, const MatrixXd& node, const MatrixXd& Vn) {
    
    const int nDPN = 5; // # of DOFs per node
    const int nNPE = 4; // # of nodes per element
    const int nGP = 2;  // # of Gauss points
    
    MatrixXd Me = MatrixXd::Zero(nDPN * nNPE, nDPN * nNPE);
    
    for (int i = 0; i < nGP; ++i) {
        for (int j = 0; j < nGP; ++j) {
            for (int k = 0; k < 2; ++k) {
                
                // ! Gauss integration point and weight factor
                double r, s, t, weight;
                std::tie(r, s, t, weight) = Gauss22_Point(i, j, k);
                
                // ! Jacobian
                MatrixXd Jacob = Jacobian(r, s, t, thick, node, Vn);
                
                // ! uvw: shape function matrix for displacement
                MatrixXd uvw = Displacement(r, s, t, thick, Vn);
                
                // Me = Me + weight * det(Jacob) * thick * transpose(uvw) * uvw * density;
                Me += weight * det(Jacob) * uvw.transpose() * uvw * density;
            }
        }
    }
    
    return Me;
}

 

 

matlab

function Vn = Make_Avg_Vn(Node,etable)
    NN = size(Node,1);       
    EN = size(etable,1);     
    nNPE = size(etable,2);   
    vn = zeros(nNPE,3);      
    Vn = zeros(NN,3);        
    
    for i = 1:EN
        enode = Node(etable(i,:),:);    
    
        idx1 = [2:nNPE 1];             
        idx2 = [nNPE 1:nNPE-1];        
        for j = 1:nNPE
            v1 = enode(idx1(j),:) - enode(j,:);  
            v2 = enode(idx2(j),:) - enode(j,:);   
            vc = cross(v1,v2);                    
            vn(j,:) = vc/norm(vc);                
        end
        Vn(etable(i,:),:) = Vn(etable(i,:),:) + vn;   
    end
    Vn = Vn./vecnorm(Vn,2,2);      
end

function [V1,V2] = Director(Vn)
            
    e2 = [0 1 0];
    
    if (Equal_Vector(Vn,e2))||(Equal_Vector(Vn,-e2))
        V1 = [0 0 1];
    else % for the special case Vn = e2
        V1 = cross(e2,Vn);
        V1 = Normalize(V1);
    end
    
    V2 = cross(Vn,V1);
    
end
    
function equal = Equal_Vector(Va, Vb)
     

    cond = abs(Va(1) - Vb(1)) < eps & ...
        abs(Va(2) - Vb(2)) < eps & ...
        abs(Va(3) - Vb(3)) < eps;
    
    if cond
        equal = true;
    else
        equal = false;
    end
    
end
    
function v1 = Normalize(v)
        
    v1 = v / sqrt((v(1)^2 + v(2)^2 + v(3)^2));
    
end

function C = Material_Law(Young, Poisson)

    C = zeros(6,6);
    k = 1;
    C(1,1) = Young / (1 - Poisson^2);
    C(2,2) = C(1,1);
    C(1,2) = Poisson * C(1,1);
    C(2,1) = C(1,2);
    C(4,4) = Young / (1 + Poisson) /2;
    C(5,5) = k * C(4,4);
    C(6,6) = C(5,5);
end

function [r, s, t, w] = Gauss22_Point(i, j, k)

    gauss_point2 = [-0.577350269189626, 0.577350269189626];
    
    w = 1.0;
    r = gauss_point2(i);
    s = gauss_point2(j);
    t = gauss_point2(k);
    
end

function dHrs = dHrs_Matrix(r,s)
    
    dHrs = 0.25*[1+s, -(1+s), -(1-s), 1-s; 1+r, 1-r, -(1-r), -(1+r)];
    
end

function H = Shape_Function(r,s)
    
    H = 0.25*[(1+r)*(1+s) (1-r)*(1+s) (1-r)*(1-s) (1+r)*(1-s)];
    
end

function [e1,e2,e3] = Local_Basis(Jacob)

    g2 = Jacob(2,:);
    g3 = Jacob(3,:);
    
    e3 = g3;
    e3 = Normalize(e3);
    
    e1 = cross(g2,g3);
    e1 = Normalize(e1);

    e2 = cross(g3,e1);
    e2 = Normalize(e2);
    
end

function [cg1,cg2,cg3] = Contravariant_Basis(Jacob)
    
    g1 = Jacob(1,:);
    g2 = Jacob(2,:);
    g3 = Jacob(3,:);
    
    detJ = det(Jacob);
    
    cg1 = cross(g2,g3) / detJ;
    cg2 = cross(g3,g1) / detJ;
    cg3 = cross(g1,g2) / detJ;
    
end

function [g1,g2,g3] = Covariant_Basis(Jacob)
    
    g1 = Jacob(1,:);
    g2 = Jacob(2,:);
    g3 = Jacob(3,:);
    
end

function Jacob = Jacobian(r,s,t,thick,node,Vn)
    
    H = Shape_Function(r,s);
    dHrs = dHrs_Matrix(r,s);
    
    nDIM = 3;
    Jacob = zeros(nDIM,nDIM);
    for i = 1:nDIM
        Jacob(1,i) = dot(dHrs(1,:),node(:,i)) + 0.5*t*dot(dHrs(1,:),thick*Vn(:,i));
        Jacob(2,i) = dot(dHrs(2,:),node(:,i)) + 0.5*t*dot(dHrs(2,:),thick*Vn(:,i));
        Jacob(3,i) = 0.5*dot(thick*H,Vn(:,i));
    end
end

 

C++

#include <Eigen/Dense>
#include <vector>
#include <cmath>

using namespace Eigen;
using namespace std;

MatrixXd Make_Avg_Vn(const MatrixXd& Node, const MatrixXi& etable) {
    int NN = Node.rows();
    int EN = etable.rows();
    int nNPE = etable.cols();
    MatrixXd vn = MatrixXd::Zero(nNPE, 3);
    MatrixXd Vn = MatrixXd::Zero(NN, 3);

    for (int i = 0; i < EN; i++) {
        MatrixXd enode = Node.row(etable(i, 0)).transpose();
        for (int j = 1; j < nNPE; j++) {
            enode.conservativeResize(3, enode.cols() + 1);
            enode.col(j) = Node.row(etable(i, j)).transpose();
        }

        VectorXi idx1(nNPE);
        VectorXi idx2(nNPE);
        for (int j = 0; j < nNPE; j++) {
            idx1(j) = (j + 1) % nNPE;
            idx2(j) = (j - 1 + nNPE) % nNPE;
        }

        for (int j = 0; j < nNPE; j++) {
            Vector3d v1 = enode.row(idx1(j)) - enode.row(j);
            Vector3d v2 = enode.row(idx2(j)) - enode.row(j);
            Vector3d vc = v1.cross(v2);
            vn.row(j) = vc.normalized();
        }
        Vn.row(etable.row(i)) += vn;
    }
    for (int i = 0; i < NN; i++) {
        Vn.row(i) = Vn.row(i).normalized();
    }
    return Vn;
}

void Director(const Vector3d& Vn, Vector3d& V1, Vector3d& V2) {
    Vector3d e2(0, 1, 0);

    if (Equal_Vector(Vn, e2) || Equal_Vector(Vn, -e2)) {
        V1 = Vector3d(0, 0, 1);
    } else {
        V1 = Vn.cross(e2).normalized();
    }

    V2 = Vn.cross(V1);
}

bool Equal_Vector(const Vector3d& Va, const Vector3d& Vb) {
    return (abs(Va(0) - Vb(0)) < numeric_limits<double>::epsilon() &&
            abs(Va(1) - Vb(1)) < numeric_limits<double>::epsilon() &&
            abs(Va(2) - Vb(2)) < numeric_limits<double>::epsilon());
}

Vector3d Normalize(const Vector3d& v) {
    return v.normalized();
}

MatrixXd Material_Law(double Young, double Poisson) {
    MatrixXd C = MatrixXd::Zero(6, 6);
    C(0, 0) = Young / (1 - Poisson * Poisson);
    C(1, 1) = C(0, 0);
    C(0, 1) = Poisson * C(0, 0);
    C(1, 0) = C(0, 1);
    C(3, 3) = Young / (1 + Poisson) / 2;
    C(4, 4) = C(3, 3);
    C(5, 5) = C(4, 4);
    return C;
}

void Gauss22_Point(int i, int j, int k, double& r, double& s, double& t, double& w) {
    double gauss_point2[2] = {-0.577350269189626, 0.577350269189626};

    w = 1.0;
    r = gauss_point2[i];
    s = gauss_point2[j];
    t = gauss_point2[k];
}

MatrixXd dHrs_Matrix(double r, double s) {
    return 0.25 * (MatrixXd(2, 4) << 1 + s, -(1 + s), -(1 - s), 1 - s,
                   1 + r, 1 - r, -(1 - r), -(1 + r)).finished();
}

MatrixXd Shape_Function(double r, double s) {
    return 0.25 * (MatrixXd(1, 4) << (1 + r) * (1 + s), (1 - r) * (1 + s), (1 - r) * (1 - s), (1 + r) * (1 - s)).finished();
}

void Local_Basis(const MatrixXd& Jacob, Vector3d& e1, Vector3d& e2, Vector3d& e3) {
    Vector3d g2 = Jacob.row(1);
    Vector3d g3 = Jacob.row(2);

    e3 = Normalize(g3);
    e1 = Normalize(g2.cross(g3));
    e2 = Normalize(g3.cross(e1));
}

void Contravariant_Basis(const MatrixXd& Jacob, Vector3d& cg1, Vector3d& cg2, Vector3d& cg3) {
    Vector3d g1 = Jacob.row(0);
    Vector3d g2 = Jacob.row(1);
    Vector3d g3 = Jacob.row(2);

    double detJ = Jacob.determinant();

    cg1 = g2.cross(g3) / detJ;
    cg2 = g3.cross(g1) / detJ;
    cg3 = g1.cross(g2) / detJ;
}

void Covariant_Basis(const MatrixXd& Jacob, Vector3d& g1, Vector3d& g2, Vector3d& g3) {
    g1 = Jacob.row(0);
    g2 = Jacob.row(1);
    g3 = Jacob.row(2);
}

MatrixXd Jacobian(double r, double s, double t, const MatrixXd& thick, const MatrixXd& node, const MatrixXd& Vn) {
    MatrixXd H = Shape_Function(r, s);
    MatrixXd dHrs = dHrs_Matrix(r, s);

    int nDIM = 3;
    MatrixXd Jacob = MatrixXd::Zero(nDIM, nDIM);
    for (int i = 0; i < nDIM; i++) {
        Jacob(0, i) = dHrs.row(0).dot(node.col(i)) + 0.5 * t * dHrs.row(0).dot(thick * Vn.col(i));
        Jacob(1, i) = dHrs.row(1).dot(node.col(i)) + 0.5 * t * dHrs.row(1).dot(thick * Vn.col(i));
        Jacob(2, i) = 0.5 * (thick * H).dot(Vn.col(i));
    }
    return Jacob;
}

 

matlab

function uvw = Displacement(r,s,t,thick,Vn)
    
    V1 = zeros(nNPE,nDIM);
    V2 = zeros(nNPE,nDIM);
    
    for i = 1:nNPE
        [V1(i,:), V2(i,:)] = Director(Vn(i,:));
    end
    
    H = Shape_Function(r,s);
    uvw = zeros(nDIM, nNPE*nDPN);
    
    for j = 1:nNPE
        uvw(1,nNPE*0+j) = H(j);
        uvw(2,nNPE*1+j) = H(j);
        uvw(3,nNPE*2+j) = H(j);
    end
    
    for i = 1:nDIM
        for j = 1:nNPE
            uvw(i,nNPE*3+j) = -0.5*t*thick*H(j)*V2(j,i);
            uvw(i,nNPE*4+j) =  0.5*t*thick*H(j)*V1(j,i);
        end
    end
end
        
        
function [duvwdr, duvwds, duvwdt] = Derivative_Displacement(r,s,t,thick,node,Vn)
    
    V1 = zeros(nNPE,nDIM);
    V2 = zeros(nNPE,nDIM);
    
    for i = 1:nNPE
        [V1(i,:), V2(i,:)] = Director(Vn(i,:));
    end
    
    H = Shape_Function(r,s);
    dHrs = dHrs_Matrix(r,s);
    
    duvwdr = zeros(nDIM, nNPE*nDPN);
    duvwds = zeros(nDIM, nNPE*nDPN);
    duvwdt = zeros(nDIM, nNPE*nDPN);
    
    for j = 1:nNPE
        duvwdr(1,nNPE*0+j) = dHrs(1,j);
        duvwdr(2,nNPE*1+j) = dHrs(1,j);
        duvwdr(3,nNPE*2+j) = dHrs(1,j);
    
        duvwds(1,nNPE*0+j) = dHrs(2,j);
        duvwds(2,nNPE*1+j) = dHrs(2,j);
        duvwds(3,nNPE*2+j) = dHrs(2,j);
    end
    
    for i = 1:nDIM
        for j = 1:nNPE
            duvwdr(i,nNPE*3+j) = -0.5*t*thick*dHrs(1,j)*V2(j,i);
            duvwdr(i,nNPE*4+j) =  0.5*t*thick*dHrs(1,j)*V1(j,i);

            duvwds(i,nNPE*3+j) = -0.5*t*thick*dHrs(2,j)*V2(j,i);
            duvwds(i,nNPE*4+j) =  0.5*t*thick*dHrs(2,j)*V1(j,i);
    
                   
            duvwdt(i,nNPE*3+j) = -0.5*thick*H(j)*V2(j,i);
            duvwdt(i,nNPE*4+j) =  0.5*thick*H(j)*V1(j,i);
        end
    end
    
end

function cov_B = Covariant_Strain(r,s,t,thick,node,Vn)
    
    [duvwdr, duvwds, duvwdt] = Derivative_Displacement(r,s,t,thick,node,Vn);
    

    Jacob = Jacobian(r, s, t, thick, node, Vn);
    

    [g1, g2, g3] = Covariant_Basis(Jacob);
    

    cov_B(1,:) = 0.5* (...
        g1(1)*duvwdr(1,:) + g1(2)*duvwdr(2,:) + g1(3)*duvwdr(3,:) + ...
        g1(1)*duvwdr(1,:) + g1(2)*duvwdr(2,:) + g1(3)*duvwdr(3,:) );
    cov_B(2,:) = 0.5* (...
        g2(1)*duvwds(1,:) + g2(2)*duvwds(2,:) + g2(3)*duvwds(3,:) + ...
        g2(1)*duvwds(1,:) + g2(2)*duvwds(2,:) + g2(3)*duvwds(3,:) );
    cov_B(3,:) = 0;

    cov_B(4,:) = 0.5* (...
        g2(1)*duvwdr(1,:) + g2(2)*duvwdr(2,:) + g2(3)*duvwdr(3,:) + ...
        g1(1)*duvwds(1,:) + g1(2)*duvwds(2,:) + g1(3)*duvwds(3,:) );
    cov_B(5,:) = 0.5* (...
        g3(1)*duvwds(1,:) + g3(2)*duvwds(2,:) + g3(3)*duvwds(3,:) + ...
        g2(1)*duvwdt(1,:) + g2(2)*duvwdt(2,:) + g2(3)*duvwdt(3,:) );
    cov_B(6,:) = 0.5* (...
        g1(1)*duvwdt(1,:) + g1(2)*duvwdt(2,:) + g1(3)*duvwdt(3,:) + ...
        g3(1)*duvwdr(1,:) + g3(2)*duvwdr(2,:) + g3(3)*duvwdr(3,:) );
end

function ass_B = Assumed_Covariant_Strain(r, s, t, thick, node, Vn)
    

    cov_B_st1 = Covariant_Strain( 1.0, 0.0, t, thick, node, Vn);
    cov_B_st2 = Covariant_Strain(-1.0, 0.0, t, thick, node, Vn);
    
    cov_B_tr1 = Covariant_Strain( 0.0, 1.0, t, thick, node, Vn);
    cov_B_tr2 = Covariant_Strain( 0.0,-1.0, t, thick, node, Vn);
    
    ass_B = Covariant_Strain(r, s, t, thick, node, Vn);

    ass_B(3,:) = 0.0;
    ass_B(5,:) = 0.5 * (1.0+r) * cov_B_st1(5,:) + 0.5*(1.0 - r) * cov_B_st2(5,:);
    ass_B(6,:) = 0.5 * (1.0+s) * cov_B_tr1(6,:) + 0.5*(1.0 - s) * cov_B_tr2(6,:);
    
end

C++

#include <Eigen/Dense>
#include <vector>

using namespace Eigen;
using namespace std;

MatrixXd Displacement(double r, double s, double t, double thick, const MatrixXd& Vn) {
    MatrixXd V1 = MatrixXd::Zero(nNPE, nDIM);
    MatrixXd V2 = MatrixXd::Zero(nNPE, nDIM);

    for (int i = 0; i < nNPE; i++) {
        tie(V1.row(i), V2.row(i)) = Director(Vn.row(i));
    }

    VectorXd H = Shape_Function(r, s);
    MatrixXd uvw = MatrixXd::Zero(nDIM, nNPE * nDPN);

    for (int j = 0; j < nNPE; j++) {
        uvw(0, nNPE * 0 + j) = H(j);
        uvw(1, nNPE * 1 + j) = H(j);
        uvw(2, nNPE * 2 + j) = H(j);
    }

    for (int i = 0; i < nDIM; i++) {
        for (int j = 0; j < nNPE; j++) {
            uvw(i, nNPE * 3 + j) = -0.5 * t * thick * H(j) * V2(j, i);
            uvw(i, nNPE * 4 + j) = 0.5 * t * thick * H(j) * V1(j, i);
        }
    }

    return uvw;
}

void Derivative_Displacement(double r, double s, double t, double thick, const MatrixXd& node, const MatrixXd& Vn, 
                             MatrixXd& duvwdr, MatrixXd& duvwds, MatrixXd& duvwdt) {
    MatrixXd V1 = MatrixXd::Zero(nNPE, nDIM);
    MatrixXd V2 = MatrixXd::Zero(nNPE, nDIM);

    for (int i = 0; i < nNPE; i++) {
        tie(V1.row(i), V2.row(i)) = Director(Vn.row(i));
    }

    VectorXd H = Shape_Function(r, s);
    MatrixXd dHrs = dHrs_Matrix(r, s);

    duvwdr = MatrixXd::Zero(nDIM, nNPE * nDPN);
    duvwds = MatrixXd::Zero(nDIM, nNPE * nDPN);
    duvwdt = MatrixXd::Zero(nDIM, nNPE * nDPN);

    for (int j = 0; j < nNPE; j++) {
        duvwdr(0, nNPE * 0 + j) = dHrs(0, j);
        duvwdr(1, nNPE * 1 + j) = dHrs(0, j);
        duvwdr(2, nNPE * 2 + j) = dHrs(0, j);

        duvwds(0, nNPE * 0 + j) = dHrs(1, j);
        duvwds(1, nNPE * 1 + j) = dHrs(1, j);
        duvwds(2, nNPE * 2 + j) = dHrs(1, j);
    }

    for (int i = 0; i < nDIM; i++) {
        for (int j = 0; j < nNPE; j++) {
            duvwdr(i, nNPE * 3 + j) = -0.5 * t * thick * dHrs(0, j) * V2(j, i);
            duvwdr(i, nNPE * 4 + j) = 0.5 * t * thick * dHrs(0, j) * V1(j, i);

            duvwds(i, nNPE * 3 + j) = -0.5 * t * thick * dHrs(1, j) * V2(j, i);
            duvwds(i, nNPE * 4 + j) = 0.5 * t * thick * dHrs(1, j) * V1(j, i);

            duvwdt(i, nNPE * 3 + j) = -0.5 * thick * H(j) * V2(j, i);
            duvwdt(i, nNPE * 4 + j) = 0.5 * thick * H(j) * V1(j, i);
        }
    }
}

MatrixXd Covariant_Strain(double r, double s, double t, double thick, const MatrixXd& node, const MatrixXd& Vn) {
    MatrixXd duvwdr, duvwds, duvwdt;
    Derivative_Displacement(r, s, t, thick, node, Vn, duvwdr, duvwds, duvwdt);

    MatrixXd Jacob = Jacobian(r, s, t, thick, node, Vn);

    VectorXd g1, g2, g3;
    tie(g1, g2, g3) = Covariant_Basis(Jacob);

    MatrixXd cov_B = MatrixXd::Zero(6, duvwdr.cols());

    cov_B.row(0) = 0.5 * (g1(0) * duvwdr.row(0) + g1(1) * duvwdr.row(1) + g1(2) * duvwdr.row(2));
    cov_B.row(1) = 0.5 * (g2(0) * duvwds.row(0) + g2(1) * duvwds.row(1) + g2(2) * duvwds.row(2));
    cov_B.row(2) = VectorXd::Zero(duvwdr.cols());

    cov_B.row(3) = 0.5 * (g2(0) * duvwdr.row(0) + g2(1) * duvwdr.row(1) + g2(2) * duvwdr.row(2) +
                          g1(0) * duvwds.row(0) + g1(1) * duvwds.row(1) + g1(2) * duvwds.row(2));
    cov_B.row(4) = 0.5 * (g3(0) * duvwds.row(0) + g3(1) * duvwds.row(1) + g3(2) * duvwds.row(2) +
                          g2(0) * duvwdt.row(0) + g2(1) * duvwdt.row(1) + g2(2) * duvwdt.row(2));
    cov_B.row(5) = 0.5 * (g1(0) * duvwdt.row(0) + g1(1) * duvwdt.row(1) + g1(2) * duvwdt.row(2) +
                          g3(0) * duvwdr.row(0) + g3(1) * duvwdr.row(1) + g3(2) * duvwdr.row(2));

    return cov_B;
}

MatrixXd Assumed_Covariant_Strain(double r, double s, double t, double thick, const MatrixXd& node, const MatrixXd& Vn) {
    MatrixXd cov_B_st1 = Covariant_Strain(1.0, 0.0, t, thick, node, Vn);
    MatrixXd cov_B_st2 = Covariant_Strain(-1.0, 0.0, t, thick, node, Vn);

    MatrixXd cov_B_tr1 = Covariant_Strain(0.0, 1.0, t, thick, node, Vn);
    MatrixXd cov_B_tr2 = Covariant_Strain(0.0, -1.0, t, thick, node, Vn);

    MatrixXd ass_B = Covariant_Strain(r, s, t, thick, node, Vn);

    ass_B.row(2) = VectorXd::Zero(ass_B.cols());
    ass_B.row(4) = 0.5 * (1.0 + r) * cov_B_st1.row(4) + 0.5 * (1.0 - r) * cov_B_st2.row(4);
    ass_B.row(5) = 0.5 * (1.0 + s) * cov_B_tr1.row(5) + 0.5 * (1.0 - s) * cov_B_tr2.row(5);

    return ass_B;
}

 

matlab

function B = Strain_Displacement(r, s, t, thick, node, Vn)
    
    MITC = true;
    
    if MITC == true
        cov_B = Assumed_Covariant_Strain(r, s, t, thick, node, Vn);
    else
        cov_B = Covariant_Strain(r, s, t, thick, node, Vn);
    end
    

    Jacob = Jacobian(r, s, t, thick, node, Vn);
    

    [cg1, cg2, cg3] = Contravariant_Basis(Jacob);
    
    [e1, e2, e3] = Local_Basis(Jacob);
    
 
    B(1,:) = ...
        (dot(cg1,e1)*dot(cg1,e1)) * cov_B(1,:) + ...
        (dot(cg2,e1)*dot(cg2,e1)) * cov_B(2,:) + ...
        (dot(cg1,e1)*dot(cg2,e1) + dot(cg2,e1)*dot(cg1,e1)) * cov_B(4,:) + ...
        (dot(cg2,e1)*dot(cg3,e1) + dot(cg3,e1)*dot(cg2,e1)) * cov_B(5,:) + ...
        (dot(cg3,e1)*dot(cg1,e1) + dot(cg1,e1)*dot(cg3,e1)) * cov_B(6,:);
    
    B(2,:) = ...
        (dot(cg1,e2)*dot(cg1,e2)) * cov_B(1,:) + ...
        (dot(cg2,e2)*dot(cg2,e2)) * cov_B(2,:) + ...
        (dot(cg1,e2)*dot(cg2,e2) + dot(cg2,e2)*dot(cg1,e2)) * cov_B(4,:) + ...
        (dot(cg2,e2)*dot(cg3,e2) + dot(cg3,e2)*dot(cg2,e2)) * cov_B(5,:) + ...
        (dot(cg3,e2)*dot(cg1,e2) + dot(cg1,e2)*dot(cg3,e2)) * cov_B(6,:);
    
    B(3,:) = 0.0;
    B(4,:) = ...
        (dot(cg1,e1)*dot(cg1,e2)) * cov_B(1,:) + ...
        (dot(cg2,e1)*dot(cg2,e2)) * cov_B(2,:) + ...
        (dot(cg1,e1)*dot(cg2,e2) + dot(cg2,e1)*dot(cg1,e2)) * cov_B(4,:) + ...
        (dot(cg2,e1)*dot(cg3,e2) + dot(cg3,e1)*dot(cg2,e2)) * cov_B(5,:) + ...
        (dot(cg3,e1)*dot(cg1,e2) + dot(cg1,e1)*dot(cg3,e2)) * cov_B(6,:);
    
    B(5,:) = ...
        (dot(cg1,e2)*dot(cg1,e3)) * cov_B(1,:) + ...
        (dot(cg2,e2)*dot(cg2,e3)) * cov_B(2,:) + ...
        (dot(cg1,e2)*dot(cg2,e3) + dot(cg2,e2)*dot(cg1,e3)) * cov_B(4,:) + ...
        (dot(cg2,e2)*dot(cg3,e3) + dot(cg3,e2)*dot(cg2,e3)) * cov_B(5,:) + ...
        (dot(cg3,e2)*dot(cg1,e3) + dot(cg1,e2)*dot(cg3,e3)) * cov_B(6,:);

    B(6,:) = ...
        (dot(cg1,e3)*dot(cg1,e1)) * cov_B(1,:) + ...
        (dot(cg2,e3)*dot(cg2,e1)) * cov_B(2,:) + ...
        (dot(cg1,e3)*dot(cg2,e1) + dot(cg2,e3)*dot(cg1,e1)) * cov_B(4,:) + ...
        (dot(cg2,e3)*dot(cg3,e1) + dot(cg3,e3)*dot(cg2,e1)) * cov_B(5,:) + ...
        (dot(cg3,e3)*dot(cg1,e1) + dot(cg1,e3)*dot(cg3,e1)) * cov_B(6,:);
    
    B(4:6,:) = 2.0 * B(4:6,:);
    
end

function Ke = Shell_Stiffness(Young, Poisson, thick, node, Vn)

    C = Material_Law(Young, Poisson);
    
    Ke = zeros(nNPE*nDPN, nNPE*nDPN);
    
    for i = 1:nGP
        for j = 1:nGP
            for k = 1:2
    
                [r, s, t, weight] = Gauss22_Point(i, j, k);
    

                B = Strain_Displacement(r, s, t, thick, node, Vn);
    

                Jacob = Jacobian(r, s, t, thick, node, Vn);
    

                Ke = Ke + weight * abs(det(Jacob)) * B'*C*B;
            end
        end
    end
end

C++

#include <Eigen/Dense>
#include <vector>

using namespace Eigen;
using namespace std;

MatrixXd Strain_Displacement(double r, double s, double t, double thick, const MatrixXd& node, const VectorXd& Vn) {
    bool MITC = true;
    MatrixXd cov_B;

    if (MITC) {
        cov_B = Assumed_Covariant_Strain(r, s, t, thick, node, Vn);
    } else {
        cov_B = Covariant_Strain(r, s, t, thick, node, Vn);
    }

    MatrixXd Jacob = Jacobian(r, s, t, thick, node, Vn);
    VectorXd cg1, cg2, cg3;
    tie(cg1, cg2, cg3) = Contravariant_Basis(Jacob);
    
    VectorXd e1, e2, e3;
    tie(e1, e2, e3) = Local_Basis(Jacob);
    
    MatrixXd B(6, cov_B.cols());
    
    B.row(0) = (cg1.dot(e1) * cg1.dot(e1)) * cov_B.row(0) +
               (cg2.dot(e1) * cg2.dot(e1)) * cov_B.row(1) +
               (cg1.dot(e1) * cg2.dot(e1) + cg2.dot(e1) * cg1.dot(e1)) * cov_B.row(3) +
               (cg2.dot(e1) * cg3.dot(e1) + cg3.dot(e1) * cg2.dot(e1)) * cov_B.row(4) +
               (cg3.dot(e1) * cg1.dot(e1) + cg1.dot(e1) * cg3.dot(e1)) * cov_B.row(5);
    
    B.row(1) = (cg1.dot(e2) * cg1.dot(e2)) * cov_B.row(0) +
               (cg2.dot(e2) * cg2.dot(e2)) * cov_B.row(1) +
               (cg1.dot(e2) * cg2.dot(e2) + cg2.dot(e2) * cg1.dot(e2)) * cov_B.row(3) +
               (cg2.dot(e2) * cg3.dot(e2) + cg3.dot(e2) * cg2.dot(e2)) * cov_B.row(4) +
               (cg3.dot(e2) * cg1.dot(e2) + cg1.dot(e2) * cg3.dot(e2)) * cov_B.row(5);
    
    B.row(2) = VectorXd::Zero(cov_B.cols());
    
    B.row(3) = (cg1.dot(e1) * cg1.dot(e2)) * cov_B.row(0) +
               (cg2.dot(e1) * cg2.dot(e2)) * cov_B.row(1) +
               (cg1.dot(e1) * cg2.dot(e2) + cg2.dot(e1) * cg1.dot(e2)) * cov_B.row(3) +
               (cg2.dot(e1) * cg3.dot(e2) + cg3.dot(e1) * cg2.dot(e2)) * cov_B.row(4) +
               (cg3.dot(e1) * cg1.dot(e2) + cg1.dot(e1) * cg3.dot(e2)) * cov_B.row(5);
    
    B.row(4) = (cg1.dot(e2) * cg1.dot(e3)) * cov_B.row(0) +
               (cg2.dot(e2) * cg2.dot(e3)) * cov_B.row(1) +
               (cg1.dot(e2) * cg2.dot(e3) + cg2.dot(e2) * cg1.dot(e3)) * cov_B.row(3) +
               (cg2.dot(e2) * cg3.dot(e3) + cg3.dot(e2) * cg2.dot(e3)) * cov_B.row(4) +
               (cg3.dot(e2) * cg1.dot(e3) + cg1.dot(e2) * cg3.dot(e3)) * cov_B.row(5);
    
    B.row(5) = (cg1.dot(e3) * cg1.dot(e1)) * cov_B.row(0) +
               (cg2.dot(e3) * cg2.dot(e1)) * cov_B.row(1) +
               (cg1.dot(e3) * cg2.dot(e1) + cg2.dot(e3) * cg1.dot(e1)) * cov_B.row(3) +
               (cg2.dot(e3) * cg3.dot(e1) + cg3.dot(e3) * cg2.dot(e1)) * cov_B.row(4) +
               (cg3.dot(e3) * cg1.dot(e1) + cg1.dot(e3) * cg3.dot(e1)) * cov_B.row(5);
    
    B.row(3) *= 2.0;
    B.row(4) *= 2.0;
    B.row(5) *= 2.0;

    return B;
}

MatrixXd Shell_Stiffness(double Young, double Poisson, double thick, const MatrixXd& node, const VectorXd& Vn) {
    MatrixXd C = Material_Law(Young, Poisson);
    int nNPE = node.rows(); // Assuming nNPE is the number of nodes
    int nDPN = node.cols(); // Assuming nDPN is the number of degrees of freedom per node
    MatrixXd Ke = MatrixXd::Zero(nNPE * nDPN, nNPE * nDPN);
    
    for (int i = 0; i < nGP; ++i) {
        for (int j = 0; j < nGP; ++j) {
            for (int k = 0; k < 2; ++k) {
                double r, s, t, weight;
                tie(r, s, t, weight) = Gauss22_Point(i, j, k);
                
                MatrixXd B = Strain_Displacement(r, s, t, thick, node, Vn);
                MatrixXd Jacob = Jacobian(r, s, t, thick, node, Vn);
                
                Ke += weight * abs(Jacob.determinant()) * B.transpose() * C * B;
            }
        }
    }
    
    return Ke;
}

 

matlab

Nel = size(etable_4,1); % shell 수 
K = sparse(tdof,tdof);

for i = 1:Nel
    if nnz(etable_4(i,:))==4 % 사각형 이니까 
        nNPE = 4;
    else
        nNPE = 3;
    end

    eNode = node(etable_4(i,1:nNPE),:); % 노드 위치 정보 뽑기 
    eVn = Vn(etable_4(i,1:nNPE),:); % 노말 벡터 
    
    if nNPE==4
        [ke, me] = MITC4(E, nu, thk, rho, eNode, eVn);
    else
        [ke, me] = MITC3(E, nu, thk, rho, eNode, eVn);
    end

    edofm = zeros(Ndof,nNPE);
    for j = 1:Ndof % 5 자유도까지 있다.
        edofm(j,1:nNPE) = Ndof*(etable_4(i,1:nNPE)-1)+j;
    end
    edof = edofm(:); % 위치 지정해주기 
    K(edof,edof) = K(edof,edof) + ke;% k 메트릭스 채우기 
end

c++

#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <vector>

using namespace Eigen;
using namespace std;

int Nel = etable_4.rows(); // shell 수 
SparseMatrix<double> K(tdof, tdof);

for (int i = 0; i < Nel; i++) {
    int nNPE;
    if (etable_4.row(i).nonZeros() == 4) { // 사각형 이니까 
        nNPE = 4;
    } else {
        nNPE = 3;
    }

    MatrixXd eNode = node.row(etable_4(i, Eigen::seq(0, nNPE - 1))); // 노드 위치 정보 뽑기 
    MatrixXd eVn = Vn.row(etable_4(i, Eigen::seq(0, nNPE - 1))); // 노말 벡터 

    MatrixXd ke, me;
    if (nNPE == 4) {
        tie(ke, me) = MITC4(E, nu, thk, rho, eNode, eVn);
    } else {
        tie(ke, me) = MITC3(E, nu, thk, rho, eNode, eVn);
    }

    MatrixXd edofm = MatrixXd::Zero(Ndof, nNPE);
    for (int j = 0; j < Ndof; j++) { // 5 자유도까지 있다.
        edofm.row(j).head(nNPE) = Ndof * (etable_4.row(i).head(nNPE).array() - 1) + j + 1;
    }
    VectorXd edof = Map<VectorXd>(edofm.data(), Ndof * nNPE); // 위치 지정해주기 
    K.coeffRef(edof.cast<int>(), edof.cast<int>()) += ke; // k 메트릭스 채우기 
}

 

MATLAB

clear, clc, close all

filename = 'ex1_MITC4_static.hob';

[MATERIAL, PSHELL, NODE, ~, SHELL_4, BC, FORCE] = InputReader(filename);

node = NODE(:, 2:4); 
etable_4 = SHELL_4(:, 3:6); 
fnode = FORCE(:,1); 
fix_node = BC(:,1); 

E = MATERIAL(1, 2);
nu = MATERIAL(1, 3);
rho = MATERIAL(1, 4);
thk = PSHELL(1, 3);

Nn = size(node,1);
Ndof = 5;
tdof = Ndof*Nn; 

Vn = Make_Avg_Vn(node,etable_4);      

load ex1_static.csv
load ex1_static_etable.txt

Nn_abaqus = size(ex1_static,1);
node_abaqus = ex1_static(:,2:4);
mag_abaqus = ex1_static(:,5);
sol_abaqus = ex1_static(:,6:8);
etable_abaqus = ex1_static_etable(:,2:5);



Nel = size(etable_4,1);
K = sparse(tdof,tdof);

for i = 1:Nel
    if nnz(etable_4(i,:))==4  
        nNPE = 4;
    else
        nNPE = 3;
    end

    eNode = node(etable_4(i,1:nNPE),:);
    eVn = Vn(etable_4(i,1:nNPE),:); 
    
    if nNPE==4
        [ke, me] = MITC4(E, nu, thk, rho, eNode, eVn);
    else
        [ke, me] = MITC3(E, nu, thk, rho, eNode, eVn);
    end

    edofm = zeros(Ndof,nNPE);
    for j = 1:Ndof
        edofm(j,1:nNPE) = Ndof*(etable_4(i,1:nNPE)-1)+j;
    end
    edof = edofm(:); 
    K(edof,edof) = K(edof,edof) + ke;
end

%% Apply Boundary Condition
F = zeros(tdof,1);
F(Ndof*(fnode-1)+FORCE(1,2)) = FORCE(1,3);
d = zeros(tdof,1);

fix_node = fix_node';
fix_dofm = zeros(Ndof,length(fix_node));
for i = 1:Ndof
    fix_dofm(i,:) = Ndof*(fix_node-1)+i;
end
fix_dof = fix_dofm(:);

%% Static Analysis
dof_a = fix_dof;
dof_b = setdiff(1:tdof,dof_a);

Ub = K(dof_b,dof_b)\F(dof_b);
d(dof_b) = Ub;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Vn = Make_Avg_Vn(Node,etable)
    NN = size(Node,1);       % NN: 전체 노드 수
    EN = size(etable,1);     % EN: 전체 요소 수
    nNPE = size(etable,2);   % nNPE: 각 요소의 노드 수
    vn = zeros(nNPE,3);      % 각 요소의 법선 벡터를 저장하는 행렬 (nNPE x 3 크기)
    Vn = zeros(NN,3);        % 각 노드의 평균 법선 벡터를 저장하는 행렬 (NN x 3 크기)
    
    for i = 1:EN
        enode = Node(etable(i,:),:);    % 요소의 노드 좌표 가져오기
    
        idx1 = [2:nNPE 1];             % 순환 인덱스를 위한 인덱스 배열 1
        idx2 = [nNPE 1:nNPE-1];        % 순환 인덱스를 위한 인덱스 배열 2
        for j = 1:nNPE
            v1 = enode(idx1(j),:) - enode(j,:);   % 첫 번째 벡터
            v2 = enode(idx2(j),:) - enode(j,:);   % 두 번째 벡터
            vc = cross(v1,v2);                    % 두 벡터의 외적(법선 벡터)
            vn(j,:) = vc/norm(vc);                % 법선 벡터를 단위 벡터로 정규화하여 저장
        end
        Vn(etable(i,:),:) = Vn(etable(i,:),:) + vn;   % 노드의 법선 벡터에 요소의 법선 벡터를 누적
    end
    Vn = Vn./vecnorm(Vn,2,2);       % 모든 법선 벡터를 정규화

end

 

C++

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <vector>
#include <fstream>
#include <sstream>

using namespace Eigen;
using namespace std;

void Make_Avg_Vn(const MatrixXd& Node, const MatrixXi& etable, MatrixXd& Vn) {
    int NN = Node.rows();       // NN: total number of nodes
    int EN = etable.rows();     // EN: total number of elements
    int nNPE = etable.cols();   // nNPE: number of nodes per element
    MatrixXd vn(nNPE, 3);      // matrix to store normal vectors for each element (nNPE x 3)

    for (int i = 0; i < EN; ++i) {
        MatrixXd enode = Node.row(etable.row(i).cast<int>()).eval(); // get node coordinates of the element

        VectorXi idx1(nNPE);
        VectorXi idx2(nNPE);
        for (int j = 0; j < nNPE; ++j) {
            idx1(j) = (j + 1) % nNPE; // cyclic index array 1
            idx2(j) = (j == 0) ? nNPE - 1 : j - 1; // cyclic index array 2
        }

        for (int j = 0; j < nNPE; ++j) {
            Vector3d v1 = enode.row(idx1(j)) - enode.row(j);   // first vector
            Vector3d v2 = enode.row(idx2(j)) - enode.row(j);   // second vector
            Vector3d vc = v1.cross(v2);                          // cross product (normal vector)
            vn.row(j) = vc.normalized();                         // normalize and store the normal vector
        }
        Vn.row(etable.row(i).cast<int>()) += vn;   // accumulate the normal vectors to the nodes
    }
    Vn = Vn.array().rowwise().normalized(); // normalize all normal vectors
}

int main() {
    string filename = "ex1_MITC4_static.hob";
    // Assuming InputReader is defined elsewhere and returns the required matrices
    MatrixXd MATERIAL, PSHELL, NODE, SHELL_4, BC, FORCE;
    tie(MATERIAL, PSHELL, NODE, ignore, SHELL_4, BC, FORCE) = InputReader(filename);

    MatrixXd node = NODE.rightCols(3); 
    MatrixXd etable_4 = SHELL_4.rightCols(4); 
    VectorXi fnode = FORCE.col(0).cast<int>(); 
    VectorXi fix_node = BC.col(0).cast<int>(); 

    double E = MATERIAL(0, 1);
    double nu = MATERIAL(0, 2);
    double rho = MATERIAL(0, 3);
    double thk = PSHELL(0, 2);

    int Nn = node.rows();
    int Ndof = 5;
    int tdof = Ndof * Nn; 

    MatrixXd Vn(Nn, 3);
    Make_Avg_Vn(node, etable_4.cast<int>(), Vn);      

    MatrixXd ex1_static;
    MatrixXd ex1_static_etable;
    // Load ex1_static.csv and ex1_static_etable.txt into ex1_static and ex1_static_etable respectively

    int Nn_abaqus = ex1_static.rows();
    MatrixXd node_abaqus = ex1_static.rightCols(3);
    VectorXd mag_abaqus = ex1_static.col(4);
    MatrixXd sol_abaqus = ex1_static.rightCols(3);
    MatrixXd etable_abaqus = ex1_static_etable.rightCols(4);

    int Nel = etable_4.rows();
    SparseMatrix<double> K(tdof, tdof);

    for (int i = 0; i < Nel; ++i) {
        int nNPE = (etable_4.row(i).array() != 0).count() == 4 ? 4 : 3;

        MatrixXd eNode = node.row(etable_4.row(i).head(nNPE).cast<int>());
        MatrixXd eVn = Vn.row(etable_4.row(i).head(nNPE).cast<int>());

        MatrixXd ke, me;
        if (nNPE == 4) {
            tie(ke, me) = MITC4(E, nu, thk, rho, eNode, eVn);
        } else {
            tie(ke, me) = MITC3(E, nu, thk, rho, eNode, eVn);
        }

        MatrixXi edofm(Ndof, nNPE);
        for (int j = 0; j < Ndof; ++j) {
            edofm.row(j) = Ndof * (etable_4.row(i).head(nNPE).cast<int>() - 1).transpose() + j + 1;
        }
        VectorXi edof = Map<VectorXi>(edofm.data(), Ndof * nNPE);
        for (int j = 0; j < edof.size(); ++j) {
            K.coeffRef(edof(j), edof(j)) += ke(j);
        }
    }

    // Apply Boundary Condition
    VectorXd F = VectorXd::Zero(tdof);
    F(Ndof * (fnode - 1).array() + FORCE(0, 1)) = FORCE(0, 2);
    VectorXd d = VectorXd::Zero(tdof);

    fix_node = fix_node.transpose();
    MatrixXi fix_dofm(Ndof, fix_node.size());
    for (int i = 0; i < Ndof; ++i) {
        fix_dofm.row(i) = Ndof * (fix_node - 1).array() + i + 1;
    }
    VectorXi fix_dof = Map<VectorXi>(fix_dofm.data(), Ndof * fix_node.size());

    // Static Analysis
    VectorXi dof_a = fix_dof;
    VectorXi dof_b = VectorXi::LinSpaced(tdof, 1, tdof).array().set_difference(dof_a).eval();

    VectorXd Ub = K.block(dof_b, dof_b).ldlt().solve(F(dof_b));
    d(dof_b) = Ub;

    return 0;
}

 

 

Pshell=
    1     1 0.001

Shetable=
  0   1   1  41  42  21  20
  1   2   1  62  63  42  41
398 399   1 400 401 380 379
399 400   1 421 422 401 400

bbns=


Materials=
    1 2e+11   0.3  7850

shellcoord[0]=
-0.45 -0.45  -0.5  -0.5
 0.45   0.5   0.5  0.45
    0     0     0     0


shellcoord[399]=
  0.5   0.5  0.45  0.45
 -0.5 -0.45 -0.45  -0.5
    0     0     0     0

scoord=
       1     -0.5     -0.5        0
       2     -0.5    -0.45        0
       3     -0.5     -0.4        0
       4     -0.5    -0.35        0
       5     -0.5     -0.3        0
       6     -0.5    -0.25        0
       7     -0.5     -0.2        0
       8     -0.5    -0.15        0
     432      0.5     0.05        0
     433      0.5      0.1        0
     434      0.5     0.15        0
     435      0.5      0.2        0
     436      0.5     0.25        0
     437      0.5      0.3        0
     438      0.5     0.35        0
     439      0.5      0.4        0
     440      0.5     0.45        0
     441      0.5      0.5        0