MatMechLab / AsFem

Advanced Simulation kit based on Finite Element Method (AsFem)
https://matmechlab.github.io/AsFem
GNU General Public License v3.0
180 stars 53 forks source link

Accurate calc for 2d shape fun in 3d case #36

Closed yangbai90 closed 3 years ago

yangbai90 commented 3 years ago

In the 3D scenario, a correct expression for the 2D shape funs' derivative is necessary. For the curved surface in the 3D case with the gradient in the Z-direction, the current code may fail!

The following one may fail for gradu_z in 3D case:

_DetJac=(_dydxi*_dzdeta-_dydeta*_dzdxi)*(_dydxi*_dzdeta-_dydeta*_dzdxi)
           +(_dzdxi*_dxdeta-_dzdeta*_dxdxi)*(_dzdxi*_dxdeta-_dzdeta*_dxdxi)
           +(_dxdxi*_dydeta-_dxdeta*_dydxi)*(_dxdxi*_dydeta-_dxdeta*_dydxi);
    _DetJac=sqrt(_DetJac);

    if(abs(_DetJac)<1.0e-15){
        MessagePrinter::PrintErrorTxt("singular element in 2D case, this error occurs in your 2D shape function calculation");
        MessagePrinter::AsFem_Exit();
    }

    _Jac[0][0]= _dxdxi;_Jac[0][1]= _dydxi;
    _Jac[1][0]=_dxdeta;_Jac[1][1]=_dydeta;

    _XJac[0][0]= _Jac[1][1]/_DetJac;
    _XJac[0][1]=-_Jac[0][1]/_DetJac;
    _XJac[1][0]=-_Jac[1][0]/_DetJac;
    _XJac[1][1]= _Jac[0][0]/_DetJac;

    double temp;
    for(int i=1;i<=GetShapeFunNums();i++){
        if(flag){
            temp=(*this)(i,1)*_XJac[0][0]+(*this)(i,2)*_XJac[0][1];
            (*this)(i,2)=(*this)(i,1)*_XJac[1][0]+(*this)(i,2)*_XJac[1][1];
            (*this)(i,1)=temp;
        }
       _shape_value[i-1]=(*this)(i,0);
       _shape_grad[i-1].setZero();
       _shape_grad[i-1](1)=(*this)(i,1);
       _shape_grad[i-1](2)=(*this)(i,2);
    }
yangbai90 commented 3 years ago

This issue is solve in 3c9071b59623a4a9ef0b636a5d3976385a79548b