#coding-Utf-8-*-

from gmshpy import *
from dG3Dpy import*

#script to launch PBC problem with a python script

# material law
E = 3.2E3
nu = 0.3
K = E/3./(1.-2.*nu)	# Bulk mudulus
mu =E/2./(1.+nu)	  # Shear mudulus
rho = 2.7e-9 # Bulk mass

# creation of material law

law1 = dG3DLinearElasticMaterialLaw(11,rho,E,nu)

# geometry
meshfile="rve.msh" # name of mesh file

# creation of part Domain
myfield1 = dG3DDomain(10,11,0,11,0,3)


# solver
sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 1  # number of step (used only if soltype=1)
ftime =1.   # Final time (used only if soltype=1)
tol=1.e-6  # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
system = 1 # Displacement elimination =0 Multiplier elimination = 1 Displacement+ multiplier = 2
control = 0 # load control = 0 arc length control euler = 1

# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.loadModel(meshfile)
mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(law1)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.setSystemType(system)
mysolver.setControlType(control)
mysolver.stiffnessModification(bool(1))
mysolver.iterativeProcedure(bool(1))
mysolver.setMessageView(bool(1))



#boundary condition
microBC = nonLinearDisplacementBC(1000,3)
microBC.setOrder(2)
microBC.setBCPhysical(1,3,5,2,4,6) 

 # Deformation gradient
microBC.setDeformationGradient(1.0,0.02,0.02,0.0,0.97,0,0,0,1.02)

 # Gradient of deformation gradient
microBC.setGradientOfDeformationGradient(1,1,1,0.1)
microBC.setGradientOfDeformationGradient(0,1,1,1) 
microBC.setGradientOfDeformationGradient(1,0,1,1) 
microBC.setGradientOfDeformationGradient(1,1,0,1) 
microBC.setGradientOfDeformationGradient(0,0,0,1)

mysolver.addMicroBC(microBC)

#stress averaging flag and averaging method 0- VOLUME, 1- SURFACE
mysolver.stressAveragingFlag(bool(1)) # set stress averaging ON- 0 , OFF-1
mysolver.setStressAveragingMethod(1) # 0 -volume 1- surface
#tangent averaging flag
mysolver.tangentAveragingFlag(bool(1)) # set tangent averaging ON -0, OFF -1
mysolver.setTangentAveragingMethod(1,1e-6) # 0- perturbation 1- condensation

#mysolver.setExtractPerturbationToFileFlag(0)	

# build view
mysolver.internalPointBuildView("Green-Lagrange_xx",IPField.STRAIN_XX, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_yy",IPField.STRAIN_YY, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_zz",IPField.STRAIN_ZZ, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_xy",IPField.STRAIN_XY, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_yz",IPField.STRAIN_YZ, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_xz",IPField.STRAIN_XZ, 1, 1);
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1);
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1);
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1);
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1);
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1);
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1);
mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange equivalent strain",IPField.GL_EQUIVALENT_STRAIN, 1, 1);
mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);


# solve
mysolver.solve()

# test check
check = TestCheck()
check.equal(-1.144128e+01,mysolver.getHomogenizedStress(0,0),1.e-4)
check.equal(1.840176e+01,mysolver.getHomogenizedStress(1,0),1.e-4)
check.equal(2.605796e+01,mysolver.getHomogenizedStress(2,2),1.e-4)

check.equal(3.019089e+03,mysolver.getHomogenizedTangent(0,0,0,0),1.e-4)
check.equal(1.144128e+03,mysolver.getHomogenizedTangent(0,0,1,1),1.e-4)
check.equal(9.200879e+02,mysolver.getHomogenizedTangent(0,1,0,1),1.e-4)

check.equal(1.720444e+02,mysolver.getHomogenizedSecondTangent(0,0,0,0,0,0),1.e-4)
check.equal(8.453953e+01,mysolver.getHomogenizedSecondTangent(1,0,1,1,0,1),1.e-4)
check.equal(7.866244e+01,mysolver.getHomogenizedSecondTangent(1,0,0,1,0,0),1.e-4)
check.equal(28.491995815,mysolver.getHomogenizedSecondTangent(0,1,0,1,0,0),1.e-4)