#include <stdio.h>

double caldet(double f[][3]);
double keisuu(double xy[][2],double abc[][3]);
void makeN(double abc[][3],double a2det,double n[][6],double x,double y);
void makeB(double abc[][3],double a2det,double b[][6]);
void makeD(double young,double nyu,double d[][3]);
void calsigma(double sigma[3],double d[][3],double ep[]);


int main(void)
{
    double abc[3][3];
    double xy[3][2];
    double u1,v1,u2,v2,u3,v3;
    double a2det;
    double a[3][3];
    double n[2][6];
    double x,y;
    double ux,uy;
    double b[3][6];
    double ep[3];
    double young,nyu;
    double d[3][3];
    double sigma[3];

    //三角形の各点の座標値ここから
    xy[0][0]=0.0;
    xy[0][1]=0.0;
    xy[1][0]=10.0;
    xy[1][1]=0.0;
    xy[2][0]=0.0;
    xy[2][1]=5.0;
    //三角形の各点の座標値ここまで

    a2det=keisuu(xy,abc);//keisuuをつくる

    x=10.0;//求めたいx座標
    y=0.0;//求めたいy座標
    
    makeN(abc,a2det,n,x,y);//Nマトリックスをつくる

    //検証するための変位ここから
    u1=0.0;
    v1=0.0;
    u2=4.0;
    v2=0.0;
    u3=0.0;
    v3=0.0;
    //検証するための変位ここまで
    
    ux=n[0][0]*u1+n[0][2]*u2+n[0][4]*u3;
    uy=n[1][1]*v1+n[1][3]*v2+n[1][5]*v3;
    
    printf("ux=%lf uy=%lf\n",ux,uy);

    makeB(abc,a2det,b);//Bマトリックスをつくる

    ep[0]=b[0][0]*u1+b[0][1]*v1+b[0][2]*u2+b[0][3]*v2+b[0][4]*u3+b[0][5]*v3;
    ep[1]=b[1][0]*u1+b[1][1]*v1+b[1][2]*u2+b[1][3]*v2+b[1][4]*u3+b[1][5]*v3;
    ep[2]=b[2][0]*u1+b[2][1]*v1+b[2][2]*u2+b[2][3]*v2+b[2][4]*u3+b[2][5]*v3;

    printf("exx=%10.2lf eyy=%10.2lf exy=%10.2lf\n",ep[0],ep[1],ep[2]);

    young=21000;
    nyu=0.3;

    makeD(young,nyu,d);//Dマトリックスをつくる

    calsigma(sigma,d,ep);//応力計算(calculate stress)

    return 0;
}


double caldet(double f[][3])
{
    double det;
    
det=f[0][0]*f[1][1]*f[2][2]+
    f[1][0]*f[2][1]*f[0][2]+
    f[2][0]*f[1][2]*f[0][1]-
    f[0][2]*f[1][1]*f[2][0]-
    f[0][0]*f[2][1]*f[1][2]-
    f[0][1]*f[1][0]*f[2][2];

    printf("%10.2lf\n",det);    

    return det;
}


double keisuu(double xy[][2],double abc[][3])
{
    double a2det;
    double a[3][3];

    abc[0][0]=xy[1][0]*xy[2][1]-xy[2][0]*xy[1][1];
    abc[0][1]=xy[1][1]-xy[2][1];
    abc[0][2]=xy[2][0]-xy[1][0];
    abc[1][0]=xy[2][0]*xy[0][1]-xy[0][0]*xy[2][1];
    abc[1][1]=xy[2][1]-xy[0][1];
    abc[1][2]=xy[0][0]-xy[2][0];
    abc[2][0]=xy[0][0]*xy[1][1]-xy[1][0]*xy[0][1];
    abc[2][1]=xy[0][1]-xy[1][1];
    abc[2][2]=xy[1][0]-xy[0][0];

    a[0][0]=1.0;
    a[1][0]=1.0;
    a[2][0]=1.0;
    a[0][1]=xy[0][0];
    a[1][1]=xy[1][0];
    a[2][1]=xy[2][0];
    a[0][2]=xy[0][1];
    a[1][2]=xy[1][1];
    a[2][2]=xy[2][1];

    a2det=caldet(a);

    return a2det;
}


void makeN(double abc[][3],double a2det,double n[][6],double x,double y)
{
    double n1,n2,n3;

    n1=(abc[0][0]+abc[0][1]*x+abc[0][2]*y)/a2det;
    n2=(abc[1][0]+abc[1][1]*x+abc[1][2]*y)/a2det;
    n3=(abc[2][0]+abc[2][1]*x+abc[2][2]*y)/a2det;

    n[0][0]=n1;
    n[0][1]=0.0;
    n[0][2]=n2;
    n[0][3]=0.0;
    n[0][4]=n3;
    n[0][5]=0.0;
    n[1][0]=0.0;
    n[1][1]=n1;
    n[1][2]=0.0;
    n[1][3]=n2;
    n[1][4]=0.0;
    n[1][5]=n3;

}


void makeB(double abc[][3],double a2det,double b[][6])
{
    double dn1dx,dn2dx,dn3dx;
    double dn1dy,dn2dy,dn3dy;

    dn1dx=abc[0][1]/a2det;
    dn2dx=abc[1][1]/a2det;
    dn3dx=abc[2][1]/a2det;
    dn1dy=abc[0][2]/a2det;
    dn2dy=abc[1][2]/a2det;
    dn3dy=abc[2][2]/a2det;
    
    b[0][0]=dn1dx;
    b[0][1]=0.0;
    b[0][2]=dn2dx;
    b[0][3]=0.0;
    b[0][4]=dn3dx;
    b[0][5]=0.0;
    b[1][0]=0.0;
    b[1][1]=dn1dy;
    b[1][2]=0.0;
    b[1][3]=dn2dy;
    b[1][4]=0.0;
    b[1][5]=dn3dy;
    b[2][0]=dn1dy;
    b[2][1]=dn1dx;
    b[2][2]=dn2dy;
    b[2][3]=dn2dx;
    b[2][4]=dn3dy;
    b[2][5]=dn3dx;

}


void makeD(double young,double nyu,double d[][3])
{
    double c;
    
    c=young/(1-nyu*nyu);

    d[0][0]=c*1.0;
    d[0][1]=c*nyu;
    d[0][2]=0.0;
    d[1][0]=c*nyu;
    d[1][1]=c*1.0;
    d[1][2]=0.0;
    d[2][0]=0.0;
    d[2][1]=0.0;
    d[2][2]=c*(1-nyu)/2;
}


void calsigma(double sigma[],double d[][3],double ep[])
{
    sigma[0]=d[0][0]*ep[0]+d[0][1]*ep[1]+d[0][2]*ep[2];
    sigma[1]=d[1][0]*ep[0]+d[1][1]*ep[1]+d[1][2]*ep[2];
    sigma[2]=d[2][0]*ep[0]+d[2][1]*ep[1]+d[2][2]*ep[2];

    printf("sxx=%10.2lf syy=%10.2lf sxy=%10.2lf\n",sigma[0],sigma[1],sigma[2]);
}