From 456d0fefecd5611d7ba49aee23f104a112bdd05c Mon Sep 17 00:00:00 2001
From: gabriel <gabriel.zarzoso@imdea.org>
Date: Tue, 20 Feb 2024 12:21:23 +0100
Subject: [PATCH] running
---
NonLinearSolver/materialLaw/imdeacp.f | 66 ++++++++++++++-----
.../materialLaw/library_crysplas_CAPSUL_v0.f | 25 +++++--
.../materialLaw/mlawIMDEACPUMAT.cpp | 28 ++++----
.../materialLaw/mlawUMATInterface.cpp | 4 ++
dG3D/benchmarks/imdeaCP/crystal.prop | 4 +-
dG3D/benchmarks/imdeaCP/cube.py | 8 +--
6 files changed, 94 insertions(+), 41 deletions(-)
diff --git a/NonLinearSolver/materialLaw/imdeacp.f b/NonLinearSolver/materialLaw/imdeacp.f
index 6f49aec18..a051f44df 100755
--- a/NonLinearSolver/materialLaw/imdeacp.f
+++ b/NonLinearSolver/materialLaw/imdeacp.f
@@ -46,12 +46,15 @@ C
C INCLUDE 'ABA_PARAM.INC'
- DIMENSiON STRESS(NTENS),STATEV(NSTATV),
+ REAL*8 STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),DFGRD0(3,3),COORDS(3),DROT(3,3),
4 DFGRD1(3,3)
+ REAL*8 SSE, SPD, SCD, RPL, DTIME, TEMP, CELENT, DRPLDT
+ REAL*8 PNEWDT, DTEMP, R, dR, ddR
+
CHARACTER*8 CMNAME
C List of internal variables:
@@ -225,6 +228,8 @@ c CALL GETCWD(OUTDIR)
1 toler,toler_jac,nincmax,nincmax_jac,strain_increment,
1 implicit_hard)
+c write(*,*)'properties read in main umat subroutine'
+
CLOSE(unit=UR0)
c Generation of a fourth order stiffness matrix STIFF
@@ -247,7 +252,8 @@ C STEP 1: READ INTEGRATION POINT PROPERTIES (here only orientation)
C DONE FOR EACH INTEGRATION POINT AT EVERY TIME (using static variables could reduce operation but will result in large data)
CALL form_orient_MAT(props,orient_MAT,orient_MAT_tr)
-
+c write(*,*)'-',orient_MAT
+c write(*,*)'I have orientated the crystal'
c DEFINITION OF THE SCHMDIT MATRICES
@@ -273,10 +279,14 @@ c Orientation of schmidt tensors to material orientation, schmidt
c Orientation of stiffness matrix to material orientation, STIFF_ORIENT
CALL ORIENTATE_TENSOR4(STIFF_ORIENT,STIFF,orient_MAT_tr) ! STIFF_ORIENT is STIFF in cartesians
-
+
+c write(*,*)'STIFF',STIFF_ORIENT
+
c READ STATE VARIABLES:
-
+c write(*,*)'Start read state variables'
c If time=0 initialize variables
+c write(*,*)'TIME',time(1),time(2)
+c write(*,*)'props',props(1)
if(time(2).EQ.0d0) THEN
@@ -320,7 +330,7 @@ c If time=0 initialize variables
enddo
enddo
endif
-
+c write(*,*)'End reading state variables'
c PREDICTOR:
c Deformation gradient increment
@@ -330,7 +340,9 @@ c Deformation gradient increment
if(iflag.NE.0) WRITE(*,*)'ERROR INVERTIGN DFGRD0_inv'
CALL PMAT(DFGRD_INC,DFGRD1,DFGRD0_inv,3,3,3)
-
+
+c write(*,*)'Def grad inc',DFGRD_INC
+
c Elastic deformation gradient predictor, DGGRD_ELAS_pred0
CALL PMAT(DFGRD_ELAS_pred0,DFGRD_INC,DFGRD_ELAS_t,3,3,3)
@@ -380,6 +392,8 @@ c dgamma_prediction = 0
c GLOBAL NEWTON-RAPHSON
+c write(*,*)'Start Newton-Raphson'
+
dnorm_NR=1d20
iter=0
toolarge=0
@@ -515,7 +529,7 @@ c write(*,*)'toolarge RES',noel,npt,kinc
endif
-c 6: NON LINEAR SYSTEM RESOLUTION
+c write(*,*) '6: NON LINEAR SYSTEM RESOLUTION'
c 6.1 Form Jacobian analytic J=\partial(RESIDUAL \partial(DFGRD_ELAS_PRED)
c partial_sigma = \partial skirchoff / \partial F
@@ -592,7 +606,7 @@ c
call dzeros(BIGJAC,9+max_nsystems,9+max_nsystems)
call dzeros(BIGCORR,9+max_nsystems,1)
-c BOX11: partial Lp / partial Fe
+c write(*,*)'BOX11: partial Lp / partial Fe'
CALL TENS3333(jacob_NR,R1_FE)
do,ii=1,9
@@ -601,7 +615,7 @@ c BOX11: partial Lp / partial Fe
enddo
enddo
-c BOX12: Partial Lp/ partial dgamma
+c write(*,*)'BOX12: Partial Lp/ partial dgamma'
do,jsystem=1,nsystems
c
@@ -636,7 +650,7 @@ c
BIGJAC(9,9+jsystem)=aux33(3,2)
enddo
-c BOX21: Partial dgamma / partial Fe
+c write(*,*)'BOX21: Partial dgamma / partial Fe'
do isystem=1,nsystems
@@ -652,7 +666,7 @@ c BOX21: Partial dgamma / partial Fe
enddo
-c BOX22: Partial dgamma/dgamma
+c write(*,*)'BOX22: Partial dgamma/dgamma'
do isystem=1,nsystems
do jsystem=1,nsystems
@@ -671,17 +685,27 @@ c BOX22: Partial dgamma/dgamma
enddo
enddo
+
ntot=max_nsystems+9
nloc=9+nsystems
+c write(*,*) 'ntot,nloc',ntot,nloc
+c write(*,*) 'gauss3',ntot,nloc,BIGRES
+c write(*,*) 'ntot,nloc',ntot,nloc
+
+c write(*,*) 'going tp call gauss_3'
call gauss_3(BIGJAC,BIGRES,BIGCORR,ntot,nloc,iflag)
+
+c write(*,*) 'gauss3 called'
+
if(iflag.EQ.1) THEN
c write(*,*)'error in system of eqs'
PNEWDT=.5
return
endif
-
+
+
104 FORMAT(21(F10.2,' '))
DFGRD_ELAS_PRED(1,1)=DFGRD_ELAS_PRED(1,1)-BIGCORR(1)
@@ -696,6 +720,8 @@ c write(*,*)'error in system of eqs'
gamma_tot_pred=0d0
gamma_tot_act=0d0
+
+c write(*,*) 'dgamma'
do,isystem=1,nsystems
dgamma(isystem)=dgamma(isystem)-BIGCORR(9+isystem)
gamma_pred(isystem)=gamma_act(isystem)+abs(dgamma(isystem))
@@ -743,7 +769,8 @@ c$$$ enddo
c
c NUMERICAL PROBLEMS WITH LARGE CORRECTIONS-->PLASTIC PREDICTOR
c
-
+
+c write(*,*) 'large deformations go to plastic predictor'
call norm(dnorm,BIGCORR,nsystems+9,1)
201 IF(dnorm.GT.1d10) THEN
@@ -839,7 +866,7 @@ c$$$ enddo
100 CONTINUE
-c ROTATE STRESSES TO DEFORMED CONFIGURATION
+c write(*,*) 'ROTATE STRESSES TO DEFORMED CONFIGURATION'
CALL POLAR_DECOMP(UU,RR,CC,DFGRD_ELAS_pred)
c CALL TRANSPOSE(RR_tr,RR,3,3)
@@ -858,8 +885,11 @@ c CALL TRANSPOSE(RR_tr,RR,3,3)
STRESS(4)=1.0/Jdet*skirchoff(1,2)
STRESS(5)=1.0/Jdet*skirchoff(1,3)
STRESS(6)=1.0/Jdet*skirchoff(2,3)
-c if(noel.eq.1)then
-c
+
+c write(*,*) '=====IM NOEL===',noel
+c write(*,*)stress
+
+c if(noel.ne.35)then
c write(*,*) 'kir',skirchoff_ROT(1,1),skirchoff_ROT(1,2),
c 1 skirchoff_ROT(1,3)
c write(*,*) 'kir',skirchoff_ROT(2,1),skirchoff_ROT(2,2),
@@ -936,7 +966,7 @@ C ROTATED_ORIENT=RR*ORIENT_MAT
C SAVE acumulated shear strain
index=index+4
statev(index)=gamma_tot_pred
-
+c write(*,*)'Save accumulated shear strain'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
@@ -1592,7 +1622,7 @@ c Creating the properties for each individual system
h0(ii)=h0_set(system_set(ii))
h1(ii)=h1_set(system_set(ii))
enddo
-
+ write(*,*)'END READING PROPERTIES (crystal.prop)'
RETURN
END
diff --git a/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f b/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
index 55e6e9aa6..5a4d7a635 100644
--- a/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
+++ b/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
@@ -1552,16 +1552,27 @@ c write(*,*)' numprocess,',numprocess
real*8 props(*)
real*8 orient1(3),orient2(3),orient3(3),dnorm
real*8 orient_MAT(3,3),orient_mat_tr(3,3)
- orient1(1)=props(1)
- orient1(2)=props(2)
- orient1(3)=props(3)
+
+c write(*,*)'I am in lybrary_crysplas_CAPSUL cm3Libraries--------'
+c orient1(1)=props(1)
+c orient1(2)=props(2)
+c orient1(3)=props(3)
+ orient1(1)=0.494651
+ orient1(2)=0.869092
+ orient1(3)=0.000000
+c write(*,*)'orient',orient1
call norm(dnorm,orient1,3,1)
+c write(*,*)'dnorm',dnorm
do,i=1,3
orient1(i)=orient1(i)/dnorm
enddo
- orient2(1)=props(4)
- orient2(2)=props(5)
- orient2(3)=props(6)
+c write(*,*)'orient after dnorm',orient1
+c orient2(1)=props(4)
+c orient2(2)=props(5)
+c orient2(3)=props(6)
+ orient2(1)=-0.869092
+ orient2(2)=0.494651
+ orient2(3)=0.000000
call norm(dnorm,orient2,3,1)
do,i=1,3
orient2(i)=orient2(i)/dnorm
@@ -1578,7 +1589,7 @@ c write(*,*)' numprocess,',numprocess
enddo
call transpose(orient_MAT_tr,orient_mat,3,3)
-
+c write(*,*)'I leave lybrary_crysplas_CAPSUL cm3Libraries--------'
return
end
diff --git a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
index 7e4ab5291..d33bd712b 100644
--- a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
+++ b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
@@ -36,26 +36,32 @@ mlawIMDEACPUMAT::mlawIMDEACPUMAT(const int num, const char *propName) : mlawUMA
if(!in) Msg::Error("Cannot open the file %s! Maybe is missing ",propName);
std::string line;
std::getline(in, line,'\n');
- nsdv=atoi(line.c_str());
std::getline(in, line,'\n');
- nprops1=atoi(line.c_str());
+ nsdv=76;//nsdv=atoi(line.c_str());
+ nprops1=109;//nprops1=atoi(line.c_str());
std::cout << "number of internal variables: "<< nsdv << "; number of properties: "<< nprops1 << std::endl;
props=new double[nprops1];
int i=0;
- while (std::getline(in, line,'\n') and i<nprops1)
+ while (std::getline(in, line,',') and i<3)
{
- //std::cout << line << std::endl;
- props[i]=atof(line.c_str());
- std::cout << "property "<< i << ": "<< props[i] << std::endl;
- i++;
+ props[i]=atof(line.c_str());
+ std::cout << "elastic property "<< i << ": "<< props[i] << std::endl;
+ i++;
}
+
+
+
+
in.close();
// adapt propertie ratio and shear modulus and density
- double Kinfty = props[1];
- double Ginfty = props[2];
+ double C1111 = props[0];
+ double C1122 = props[1];
+ double C1212 = props[2];
_rho = 0.;
- _nu0 = (3*Kinfty - 2*Ginfty)/2./(3*Kinfty+Ginfty);
- _mu0 = Ginfty;
+
+ std::cout << "C1111 "<< C1111 << ", C1122 "<< C1122<< ", C1212 "<< C1212 << std::endl;
+ _mu0 = C1212/2.;
+ _nu0 = C1111/2./(_mu0+C1111);
}
mlawIMDEACPUMAT::mlawIMDEACPUMAT(const mlawIMDEACPUMAT &source) :
mlawUMATInterface(source)
diff --git a/NonLinearSolver/materialLaw/mlawUMATInterface.cpp b/NonLinearSolver/materialLaw/mlawUMATInterface.cpp
index 048d7991d..3de5f5bba 100644
--- a/NonLinearSolver/materialLaw/mlawUMATInterface.cpp
+++ b/NonLinearSolver/materialLaw/mlawUMATInterface.cpp
@@ -438,7 +438,11 @@ void mlawUMATInterface::constitutive(const STensor3& F0,
if(timestep<1.e-16) timestep=1.e-7;
double tim[2];
tim[0]=this->getTime()-timestep;
+ if(tim[0]<0.)
+ tim[0]=0.;
tim[1]=this->getTime()-timestep; //not sure, to check
+ if(tim[1]<0.)
+ tim[1]=0.;
double predef[1], dpred[1]; //array of intrpolated values of field variables
predef[0]=0.; dpred[0]=0.;
int ndi=0, nshr=0; //ndi: number of direct stress components, nshr number of engineeing shear stress component
diff --git a/dG3D/benchmarks/imdeaCP/crystal.prop b/dG3D/benchmarks/imdeaCP/crystal.prop
index e62cbeb81..e27c272d1 100644
--- a/dG3D/benchmarks/imdeaCP/crystal.prop
+++ b/dG3D/benchmarks/imdeaCP/crystal.prop
@@ -5,7 +5,7 @@
2.42D-3,0.017
# Set of slip systems, Total number of systems
1,12
-# normalv, tangent, set
+# normalv, tangent, set 84
1 , 1 , 1 , 1 , -1 , 0, 1
1 , 1 , 1 , 0 , -1 , 1, 1
1 , 1 , 1 , -1 , 0 , 1 , 1
@@ -24,4 +24,6 @@
465.5D6,598.5D6,6000D6,300D6
# Control of subroutine: TOLER, TOLER_JAC, Nmax iter, Nmax iter JAC,strain_inc, IMPLICIT HARD (yes=1)
5D-10,5D-10,50,20,1D-5,0
+# euler angles cos alpha -sin alpha 0 sin alpha cos alpha 0 0
+1. 0. 0. 0. 1. 0. 0.
diff --git a/dG3D/benchmarks/imdeaCP/cube.py b/dG3D/benchmarks/imdeaCP/cube.py
index eec64109e..3193d7dbe 100644
--- a/dG3D/benchmarks/imdeaCP/cube.py
+++ b/dG3D/benchmarks/imdeaCP/cube.py
@@ -11,7 +11,7 @@ lawnum = 11 # unique number of law
# material law
lawnum = 11 # unique number of law
-law1 = IMDEACPUMATDG3DMaterialLaw(lawnum,'material.i02');
+law1 = IMDEACPUMATDG3DMaterialLaw(lawnum,'crystal.prop');
#lawPert1 = dG3DMaterialLawWithTangentByPerturbation(law1, 1e-8)
@@ -26,7 +26,7 @@ nfield = 29 # number of the field (physical number of surface)
myfield1 = dG3DDomain(1000,nfield,space1,lawnum,fullDg,3)
#myfield1.matrixByPerturbation(1,1,1,1e-8)
-myfield1.stabilityParameters(beta1)
+#myfield1.stabilityParameters(beta1)
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
@@ -54,7 +54,7 @@ mysolver.displacementBC('Face',30,0,0.)
#mysolver.displacementBC('Face',1,0,0)
#mysolver.displacementBC('Face',1,1,0)
-mysolver.displacementBC('Face',31,0,1)
+mysolver.displacementBC('Face',31,0,0.2)
# build view
mysolver.internalPointBuildView("Strain_xx",IPField.STRAIN_XX, 1, 1);
@@ -85,4 +85,4 @@ mysolver.archivingNodeDisplacement(43,2,1)
mysolver.solve()
check = TestCheck()
-check.equal(-5.878380e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-3)
+check.equal(-5.574939e+09,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-3)
--
GitLab