material.jl

# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020

struct material
    E::Float64
    mass::Float64
    nu::Float64
    rho::Float64
    b::Float64
    h::Float64
    L::Float64
    area::Float64
    I::Float64
    J::Float64
    G::Float64
    a1::Float64
    a2::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    massInverse::Float64
    momentInertiaInverse::Float64
    inertia::Float64
    zeta::Float64
    zetaCollision::Float64
    muStatic::Float64
    muKinetic::Float64
    nomSize::Float64
    sqA1::Float64
    sqA2xIp::Float64 
    sqB1::Float64
    sqB2xFMp::Float64
    sqB3xIp::Float64
    _2xSqMxExS::Float64
    function material()
        E=2000
        mass=10
        nu=0.35
        rho=7.85e-9
        b=2.38
        h=2.38
        L=75
        area=b*h
        I=b*h^3/12
        J=b*h*(b*b+h*h)/12
        G = E * 1 / 3
        a1=E*b*h/L
        a2=G*J/L
        b1=12*E*I/(L^3)
        b2=6*E*I/(L^2)
        b3=2*E*I/(L)
        massInverse=1.0/mass
        momentInertiaInverse=1.92e-6
        inertia=1/momentInertiaInverse 
        zeta=1.0
        zetaCollision=0.0
        muStatic= 2.0
        muKinetic= 0.1
        nomSize=1.0
        sqA1=sqrt(a1) 
        sqA2xIp=sqrt(a2*L*L/6.0) 
        sqB1=sqrt(b1) 
        sqB2xFMp=sqrt(b2*L/2) 
        sqB3xIp=sqrt(b3*L*L/6.0)
        _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))
        new(E,mass,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,massInverse,momentInertiaInverse,inertia,zeta,zetaCollision,muStatic,muKinetic,nomSize,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,_2xSqMxExS)
    end
    function material(E,mass,nu,rho,b,h,L,momentInertiaInverse,zeta,zetaCollision,muStatic,muKinetic,nomSize)
        area=b*h
        I=b*h^3/12
        J=b*h*(b*b+h*h)/12
        G = E * 1 / 3
        a1=E*b*h/L
        a2=G*J/L
        b1=12*E*I/(L^3)
        b2=6*E*I/(L^2)
        b3=2*E*I/(L)
        massInverse=1.0/mass
        inertia=1/momentInertiaInverse
        sqA1=sqrt(a1) 
        sqA2xIp=sqrt(a2*L*L/6.0) 
        sqB1=sqrt(b1) 
        sqB2xFMp=sqrt(b2*L/2) 
        sqB3xIp=sqrt(b3*L*L/6.0)
        _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))
        new(E,mass,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,massInverse,momentInertiaInverse,inertia,zeta,zetaCollision,muStatic,muKinetic,nomSize,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,_2xSqMxExS)
    end
end

struct voxelMaterial
    E::Float64
    mass::Float64
    nu::Float64
    rho::Float64
    massInverse::Float64
    momentInertiaInverse::Float64
    inertia::Float64
    zetaInternal::Float64
    zetaGlobal::Float64
    zetaCollision::Float64
    muStatic::Float64
    muKinetic::Float64
    nomSize::Float64
    _2xSqMxExS::Float64
    eHat::Float64
    cTE::Float64
    linear::Bool
    poisson::Bool
    function voxelMaterial()
        E=2000
        mass=10
        nu=0.35
        rho=7.85e-9
        massInverse=1.0/mass
        nomSize=1.0
        inertia= mass*nomSize*nomSize / 6.0
        momentInertiaInverse=1/inertia
        zetaInternal=1.0
        zetaGlobal=0.2
        zetaCollision=0.0
        muStatic= 2.0
        muKinetic= 0.1
        _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))
        eHat=E/((1.0-2.0*nu)*(1.0+nu));
        cTE=0.0
        linear=true
        poisson=false
        new(E,mass,nu,rho,massInverse,momentInertiaInverse,inertia,zetaInternal,zetaGlobal,zetaCollision,muStatic,muKinetic,nomSize,_2xSqMxExS,eHat,cTE,linear,poisson)
    end
    function voxelMaterial(E,mass,nu,rho,zetaInternal,zetaGlobal,zetaCollision,muStatic,muKinetic,nomSize,linear,poisson,cTE)
        massInverse=1.0/mass
        inertia= mass*nomSize*nomSize / 6.0  #simple 1D approx
        momentInertiaInverse=1.0/inertia
        _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))
        eHat=E/((1.0-2.0*nu)*(1.0+nu));
        new(E,mass,nu,rho,massInverse,momentInertiaInverse,inertia,zetaInternal,zetaGlobal,zetaCollision,muStatic,muKinetic,nomSize,_2xSqMxExS,eHat,cTE,linear,poisson)
    end
end

struct edgeMaterial
    E::Float64
    nu::Float64
    rho::Float64
    b::Float64
    h::Float64
    L::Float64
    area::Float64
    I::Float64
    J::Float64
    G::Float64
    a1::Float64
    a2::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    sqA1::Float64
    sqA2xIp::Float64 
    sqB1::Float64
    sqB2xFMp::Float64
    sqB3xIp::Float64
    dampingM::Float64
    loaded::Float64
    strainRatio::Float64
    eHat::Float64
    cTE::Float64
    linear::Bool
    poisson::Bool
    sigmaYield::Float64
    sigmaFail::Float64
    epsilonYield::Float64
    epsilonFail::Float64
    strainData::SVector{3, Float64}
    stressData::SVector{3, Float64}
    function edgeMaterial()
        E=2000
        mass=10.0
        nu=0.35
        rho=7.85e-9
        b=2.38
        h=2.38
        L=75
        area=b*h
        I=b*h^3/12
        J=b*h*(b*b+h*h)/12
        G = E * 1 / 3
        a1=E*b*h/L
        a2=G*J/L
        b1=12*E*I/(L^3)
        b2=6*E*I/(L^2)
        b3=2*E*I/(L)
        sqA1=sqrt(a1) 
        sqA2xIp=sqrt(a2*L*L/6.0) 
        sqB1=sqrt(b1) 
        sqB2xFMp=sqrt(b2*L/2) 
        sqB3xIp=sqrt(b3*L*L/6.0)
        zeta=1.0
        dampingM=2.0*sqrt(mass)*zeta
        loaded=0.0
        strainRatio=1.0
        eHat=E/((1.0-2.0*nu)*(1.0+nu)); 
        linear=true
        poisson=false
        sigmaYield=-1.0
        sigmaFail=-1.0
        epsilonYield=-1.0
        epsilonFail=10000*E
        strainData=SVector(0.0, -1.0, -1.0)
        stressData=SVector(0.0, -1.0, -1.0)
        cTE=0.0
        new(E,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,dampingM,loaded,strainRatio,eHat,cTE,linear,poisson,sigmaYield,sigmaFail,epsilonYield,epsilonFail,strainData,stressData)
    end
    function edgeMaterial(mat,E,strainData,stressData,sigmaYield ,sigmaFail,epsilonYield,epsilonFail)
        linear=false
        new(E,mat.nu,mat.rho,mat.b,mat.h,mat.L,mat.area,mat.I,mat.J,mat.G,mat.a1,mat.a2,mat.b1,mat.b2,mat.b3,mat.sqA1,mat.sqA2xIp,mat.sqB1,mat.sqB2xFMp,mat.sqB3xIp,mat.dampingM,mat.loaded,mat.strainRatio,mat.eHat,mat.cTE,
            linear,mat.poisson,sigmaYield ,sigmaFail,epsilonYield,epsilonFail,strainData,stressData)
    end
    function edgeMaterial(E,mass,nu,rho,b,h,L,loaded,strainRatio,linear,poisson,cTE)
        area=b*h
        I=b*h^3/12
        J=b*h*(b*b+h*h)/12
        G = E  / (2.0*(1.0+nu))
        a1=E*b*h/L
        a2=G*J/L
        b1=12*E*I/(L^3)
        b2=6*E*I/(L^2)
        b3=2*E*I/(L)
        sqA1=sqrt(a1) 
        sqA2xIp=sqrt(a2*L*L/6.0) 
        sqB1=sqrt(b1) 
        sqB2xFMp=sqrt(b2*L/2) 
        sqB3xIp=sqrt(b3*L*L/6.0)
        zeta=1.0
        dampingM=2.0*sqrt(mass)*zeta
        eHat=E/((1.0-2.0*nu)*(1.0+nu))
        sigmaYield=-1.0
        sigmaFail=-1.0
        epsilonYield=-1.0
        epsilonFail=10000*E
        strainData=SVector(0.0,-1.0, -1.0)
        stressData=SVector(0.0,-1.0, -1.0)
        new(E,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,dampingM,loaded,strainRatio,eHat,cTE,linear,poisson,sigmaYield,sigmaFail,epsilonYield,epsilonFail,strainData,stressData)
    end
end

function setModelBilinear(mat, plasticModulus, yieldStress, failureStress=-1.0)
	# if (youngsModulus<=0){
	# 	error = "Young's modulus must be positive";
	# 	return false;
	# }
	# if (plasticModulus<=0 || plasticModulus>=youngsModulus){
	# 	error = "Plastic modulus must be positive but less than Young's modulus";
	# 	return false;
	# }
	# if (yieldStress<=0){
	# 	error = "Yield stress must be positive";
	# 	return false;
	# }
	# if (failureStress != -1.0f && failureStress <= yieldStress){
	# 	error = "Failure stress must be positive and greater than the yield stress";
	# 	return false;
	# }
    youngsModulus=mat.E
	yieldStrain = yieldStress / youngsModulus;
	tmpfailureStress = failureStress; #create a dummy failure stress if none was provided
    if (tmpfailureStress == -1) 
        tmpfailureStress = 3.0 * yieldStress;
    end

	tM = plasticModulus;
	tB = yieldStress-tM*yieldStrain; #y-mx=b
	tmpfailStrain = (tmpfailureStress-tB)/tM; # (y-b)/m = x


    strainData=SVector(0.0, yieldStrain, tmpfailStrain);#add in the zero data point (required always)

    stressData=SVector(0.0, yieldStress, tmpfailureStress);

	# linear = false;
	# E=youngsModulus;
	sigmaYield = yieldStress;
	sigmaFail = failureStress;
	epsilonYield = yieldStrain;
	epsilonFail = (failureStress == -1.0) ? -1.0 : tmpfailStrain
    # return updateDerived();
    mat=edgeMaterial(mat,mat.E,strainData,stressData,sigmaYield ,sigmaFail,epsilonYield,epsilonFail)

	return mat

end

function setYieldFromData(mat,strainData,stressData,percentStrainOffset,E,sigmaFail,epsilonFail)
	sigmaYield = -1.0; #assume we fail until we succeed.
	epsilonYield = -1.0; #assume we fail until we succeed.

	oM = mat.E; #the offset line slope (y=Mx+B)
	oB = (-percentStrainOffset/100.0*oM); #offset line intercept (100 factor turns percent into absolute

	# assert(strainData.size() == stressData.size());
	# assert(strainData.size() > 2); # more than 2 data points (more than bilinear)
    
    dataPoints = length(strainData)-1;
    for i = 2:dataPoints-1
		x1=strainData[i];
		x2=strainData[i+1];
		y1=stressData[i];
		y2=stressData[i+1];

		tM = (y2-y1)/(x2-x1); #temporary slope
		tB = y1-tM*x1; #temporary intercept

		if (oM!=tM) #if not parallel lines...
			xIntersect = (tB-oB)/(oM-tM);
			if (xIntersect>x1 && xIntersect<x2) #if intersects at this segment...
				percentBetweenPoints = (xIntersect-x1)/(x2-x1);
				sigmaYield = y1+percentBetweenPoints*(y2-y1);
                epsilonYield = xIntersect;
                mat=edgeMaterial(mat,E,strainData,stressData,sigmaYield ,sigmaFail,epsilonYield,epsilonFail)
				return true,mat;
            end
		end
    end
    
	sigmaYield = sigmaFail
    epsilonYield = epsilonFail
    mat=edgeMaterial(mat,E,strainData,stressData,sigmaYield ,sigmaFail,epsilonYield,epsilonFail)
	return false,mat;
end

function setModel(mat,dataPointCount, strainData, stressData)
	#Pre-checks
	# if (pStrainValues==0 && pStressValues==0)  #if first data point is 0,0, ignore it
	# 	pStrainValues++; #advance the pointers...
	# 	pStressValues++;
	# 	dataPointCount-=1; #decrement the count
    # end
	# if (dataPointCount<=0)
	# 	error = "Not enough data points";
	# 	return false;
    # end
	# if (pStrainValues<=0 || pStressValues<=0)
	# 	error = "First stress and strain data points negative or zero";
	# 	return false; 
    # end

	#Copy the data into something more usable (and check for monotonically increasing)
	# tmpStrainData=[0]
	# tmpStressData=[0]
	# tmpStrainData.push_back(0); #add in the zero data point (required always)
	# tmpStressData.push_back(0);
    # sweepStrain = 0.0
    # sweepStress = 0.0;
	# for i = 1:dataPointCount
	# 	thisStrain = *(pStrainValues+i); #grab the values
	# 	thisStress = *(pStressValues+i);

	# 	if (thisStrain <= sweepStrain)
	# 		error = "Out of order strain data";
	# 		return false;
    #     end

	# 	if (thisStress <= sweepStress)
	# 		error = "Stress data is not monotonically increasing";
    #     end

	# 	if (i>0 && (thisStress-sweepStress)/(thisStrain-sweepStrain) > tmpStressData[0]/tmpStrainData[0])
	# 		error = "Slope of stress/strain curve should never exceed that of the first line segment (youngs modulus)";
	# 		return false;
    #     end

	# 	sweepStrain = thisStrain;
	# 	sweepStress = thisStress;

	# 	tmpStrainData.push_back(thisStrain); #add to the temporary vector
	# 	tmpStressData.push_back(thisStress);
	# end


	#at this point, we know we have valid data and will return true
	# strainData = tmpStrainData;
    # stressData = tmpStressData;
    
	E=stressData[2]/strainData[2]; #youngs modulus is the inital slope
	sigmaFail = stressData[dataPointCount]; #failure stress is the highest stress data point
    epsilonFail = strainData[dataPointCount]; #failure strain is the highest strain data point
    
	linear = (dataPointCount==1);

	if (dataPointCount == 1 || dataPointCount == 2) #linear or bilinear
		sigmaYield = stressData[2];
		epsilonYield = strainData[2];
	else  #.2% (0.002) strain offset to find a good yield point...
        percentStrainOffset=0.2;
        linera,mat=setYieldFromData(mat,strainData,stressData,percentStrainOffset,E,sigmaFail,epsilonFail);
    end

	return mat
    
end