diff --git a/NonLinearSolver/internalPoints/CMakeLists.txt b/NonLinearSolver/internalPoints/CMakeLists.txt
index 48cdd82566bd822814156f875322069bd5599ddd..a399f88b53a6b1e332014d4002b1a591ca76dfc1 100644
--- a/NonLinearSolver/internalPoints/CMakeLists.txt
+++ b/NonLinearSolver/internalPoints/CMakeLists.txt
@@ -49,6 +49,7 @@ set(SRC
ipCrystalPlasticity.cpp
ipGursonUMAT.cpp
ipVEVPUMAT.cpp
+ ipIMDEACPUMAT.cpp
ipAnIsotropicElecTherMech.cpp
ipHyperelastic.cpp
ipNonLocalDamageHyperelastic.cpp
diff --git a/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..82f24f03c1dfe00ef3393fc30156e710c83a148e
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp
@@ -0,0 +1,27 @@
+//
+// Description: storing class for cp umat interface
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "ipIMDEACPUMAT.h"
+#include "ipField.h"
+#include "restartManager.h"
+
+
+double IPIMDEACPUMAT::get(int comp) const
+{
+ // for vizualization, with tag in ipField.h (only variables)
+ //if(comp==IPField::PLASTICSTRAIN)
+ // return _statev[18];
+ //else
+ return IPUMATInterface::get(comp);
+}
+
+void IPIMDEACPUMAT::restart()
+{
+ IPUMATInterface::restart();
+ return;
+}
diff --git a/NonLinearSolver/internalPoints/ipIMDEACPUMAT.h b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.h
new file mode 100644
index 0000000000000000000000000000000000000000..9cc9d477995eeb6cf5410e9ae53739b4ee31bc18
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.h
@@ -0,0 +1,47 @@
+//
+// Description: storing class for vevp umat interface
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef IPIMDEACPUMAT_H_
+#define IPIMDEACPUMAT_H_
+#include "ipUMATInterface.h"
+#include "STensor3.h"
+class IPIMDEACPUMAT : public IPUMATInterface
+{
+ protected:
+
+ public: // All data public to avoid the creation of function to access
+
+ public:
+ IPIMDEACPUMAT(int _nsdv, double eleSize) : IPUMATInterface(_nsdv, eleSize)
+ {
+
+ // here initialised internal variables which have to be
+ //_statev[0]=1;
+ //_statev[1]=1;
+ //_statev[2]=1;
+ }
+ IPIMDEACPUMAT(const IPIMDEACPUMAT &source) : IPUMATInterface(source)
+ {
+ }
+ IPIMDEACPUMAT &operator = (const IPVariable &_source)
+ {
+ IPUMATInterface::operator=(_source);
+ return *this;
+ }
+ virtual ~IPIMDEACPUMAT()
+ {
+ }
+ virtual IPVariable* clone() const {return new IPIMDEACPUMAT(*this);};
+ virtual void restart();
+ virtual double get(const int i) const;
+
+ protected:
+ IPIMDEACPUMAT(){}
+};
+
+#endif // IPIMDEACPUMAT_H_
diff --git a/NonLinearSolver/materialLaw/CMakeLists.txt b/NonLinearSolver/materialLaw/CMakeLists.txt
index a36a3d4860982ff69989fce9cbd75758aba3e011..109a929bb9d2dd4eed2eb0dcefa06b8d0fe87e0b 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -37,6 +37,8 @@ set(SRC
mlawCrystalPlasticity.cpp
mlawGursonUMAT.cpp
mlawVEVPUMAT.cpp
+ library_crysplas_CAPSUL_v0.f
+ mlawIMDEACPUMAT.cpp
additionalFortran.f
Ksub15.f
Kevpm15.f
@@ -44,6 +46,7 @@ set(SRC
ALAMEL15.f
gurson.f
vevp.f
+ imdeacp.f
LD_utils.f
mlawAnisotropic.cpp
vumat.f
diff --git a/NonLinearSolver/materialLaw/imdeacp.f b/NonLinearSolver/materialLaw/imdeacp.f
new file mode 100755
index 0000000000000000000000000000000000000000..a051f44dfd50d14d9d13f8c54cfacf8686549854
--- /dev/null
+++ b/NonLinearSolver/materialLaw/imdeacp.f
@@ -0,0 +1,1714 @@
+C
+c CAPSUL: CrystAl Plasticity UtiLities tool set
+C Multiscale Materials Modelling Group, IMDEA Materials
+C
+C UMAT CRYSTAL PLASTICITY MODEL
+c umat_crysplas_CAPSUL_isoh.f
+C
+c Short description:
+c -Phenomenlogical hardening (Asaro-Needleman + Voce)
+c -Double residual (Lp + Dgamma)-->Quadratic convergency and robust
+c -Integration of internal variables using several subincrements (now commented)
+c -NR Residual exact, numerical Jacobian
+c -EXP MAP and derivatives not used and delted. To reincorporate, see older versions
+
+c Author: Javier Segurado
+
+c Version controls and dates
+c Update 18/01/2016: Split into library with common routines and actual umat
+c Update 22/06/2015: Cleaning, rewritting, comments, from old version 5.3.8
+c
+C Subversion 5.2- 14th June 2012
+c -- EXP MAP and DEXP MAP are suppressed (linear approach)
+c -- Includes more general symmetry of stiffness matrix
+c -- write of messages supressed
+c -- Works under parallel in all ABQ versions (6.7-6.11)
+c -- Now the former file subroutines_v4.f is included in file
+C -- orientate_tensor4 updated to use explicit rotation formulae from mathematica
+c Subversion 5.3- 18th June 2012
+c -- Voce Hardening law included
+C Subversion 5.3.1-18th July 2012
+c -- Improvement of NR jacobian
+C Subversion 5.3.2-June 2013
+c -- Pseudo-line-search to improove numerical efficiency
+C Version 5.3.8. January 2015
+c -- Implicit in Fe and tau by a New residual formulation on Lp & Delta_gamma_i
+c -- Exact Jacobian and Quadratic convergency for new residual. Robust
+C
+ SUBROUTINE umat_imdeacp(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
+ 1 RPL,DDSDDT,DRPLDE,DRPLDT,
+ 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,
+ 2 DPRED,CMNAME,
+ 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,
+ 3 DROT,PNEWDT,
+ 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,
+ 5 LAYER,KSPT,KSTEP,KINC)
+
+C INCLUDE 'ABA_PARAM.INC'
+
+ 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:
+
+C To set the maximum of slip systems
+ integer max_nsystems,nsets,nsystems
+ PARAMETER (max_nsystems=48)
+
+c File management
+ CHARACTER*1 paths
+ CHARACTER*2 FLABEL
+ CHARACTER*200 outdir,filename,patho
+ INTEGER lenoutdir,nfiles,ns,UR0
+
+c MATERIAL PARAMETERS
+
+c Elastic constants
+ REAL*8 c11,c12,c44,c13,c33,c66
+c Viscoplastic parameters
+ REAL*8 mm,gamma_0
+c Hardening law
+ real*8 tau0(max_nsystems),taus(max_nsystems)
+ real*8 h0(max_nsystems),h
+c Voce hardening law
+ real*8 h1(max_nsystems)
+c Latent hardening matrix
+ real*8 q(max_nsystems,max_nsystems)
+c Normal planes and tangential planes
+ real*8 m(max_nsystems,3),s(max_nsystems,3)
+ real*8 dnorm_m,dnorm_s
+c Schmidt matrices
+ real*8 schmidt(max_nsystems,3,3)
+
+C Stiffness matrices
+ REAL*8 STIFF(3,3,3,3),stiff_orient(3,3,3,3)
+
+c Kinematic tensors: deformation gradients and strains
+ REAL*8 DFGRD_INC(3,3)
+ REAL*8 DFGRD0_inv(3,3)
+ REAL*8 DFGRD_elas(3,3)
+ REAL*8 DFGRD_elas_t(3,3)
+ REAL*8 DFGRD_elas_pred(3,3),DFGRD_elas_pred0(3,3)
+ REAL*8 DFGRD_elas_pred_tr(3,3)
+ REAL*8 DFGRD_elas_pred_inv(3,3),DFGRD_elas_pred0_inv(3,3)
+ REAL*8 EPSILON_el(3,3),Jdet,Jdet2
+ real*8 RR_tot(3,3),RR_tot_tr(3,3)
+ real*8 RR(3,3),RR_tr(3,3),UU(3,3),VV(3,3),CC(3,3)
+
+c Stresses: skirchoff_rot: 2nd Piola-Kirchoff (intermediate configuration)
+c skirchoff: Kirchoff stress (final configuration), Cauchy
+ real*8 skirchoff(3,3)
+ REAL*8 skirchoff_rot(3,3)
+
+c Velocity gradients:
+ REAL*8 L_p(3,3),L_p_new(3,3),EXP_L_p(3,3),Lterm,term
+
+c Shear on slip systems:
+ real*8 gamma_pred(max_nsystems)
+ real*8 gamma_dot(max_nsystems),gamma_act(max_nsystems)
+ real*8 gamma_dot_old(max_nsystems),gama_inc(max_nsystems)
+ real*8 gamma_dot_res(max_nsystems),error_hard
+ real*8 gamma_TOT,gamma_tot_pred,gamma_tot_act
+
+c Tau on slip systems:
+ real*8 tau_act(max_nsystems),tau_pred(max_nsystems)
+ real*8 tau_pred0(max_nsystems)
+ real*8 tau_resolved(max_nsystems)
+
+c Definition of orientation
+ REAL*8 orient1(3),orient2(3),orient3(3)
+ REAL*8 orient_MAT(3,3),orient_MAT_tr(3,3)
+ real*8 ROTATED_ORIENT(3,3),ROTATED_ORIENT_tr(3,3)
+
+c Auxiliar (tmp) integers, scalars, vectors, tensors
+ real*8 aux_escalar
+ real*8 aux3_1(3),aux3_2(3)
+ real*8 aux33(3,3),aux33_2(3,3)
+ real*8 aux3333(3,3,3,3)
+
+
+ integer i,j,ii,jj,kk,ll,pp,qq,nm,nn,iii,jjj
+ integer ia,ib,ic,id,ii1,jj1
+ integer isystem,jsystem
+ logical noconv
+ real*8 pi
+
+c NR resolution, Jacobian, etc
+ integer index , iter, kroneker,istep
+
+ integer iflag
+
+ real*8 JACOB_NR(3,3,3,3),jacob_NR_i(3,3,3,3)
+ REAL*8 dfactor(max_nsystems),dfactor2(max_nsystems)
+ real*8 mjacob(3,3,3,3),partial_sigma(3,3,3,3)
+ real*8 partial_tau(max_nsystems,3,3)
+ real*8 partial_tau2(3,3)
+ real*8 dnorm,dnorm_NR,dnorm_old,toler,toler_jac
+ real*8 dnorm_hard
+ integer nincmax,nincmax_jac
+
+ REAL*8 JACOB3333(3,3,3,3)
+
+ real*8 BIGJAC(max_nsystems+9,max_nsystems+9)
+ real*8 BIGRES(max_nsystems+9)
+ real*8 BIGCORR(max_nsystems+9),R1_FE(9,9)
+ integer ntot,nloc
+ real*8 dgamma(max_nsystems), dgamma_new(max_nsystems)
+ real*8 dgamma_jacob(max_nsystems),signog,signoT
+
+ integer toolarge
+
+c Numerical tangent stiffness matrix (DDSDDE)
+ real*8 skirchoff2(3,3),delta_eps_jac(3,3)
+ real*8 dtime2,strain_increment
+ real*8 gamma_act_jacob(max_nsystems)
+ real*8 gamma_pred_jacob(max_nsystems)
+ real*8 tau_pred_jacob(max_nsystems),tau_act_jacob(max_nsystems)
+ real*8 dfgrd_elas_pred_jacob(3,3)
+ real*8 dfgrd_inc0(3,3)
+ real*8 tinterpol,detF
+
+
+c Common Variables saved once at the beginning
+ save mm,gamma_0,tau0,taus,h0,h1
+ save q,stiff
+ save s,m
+ save nsystems
+ save pi
+ save toler,toler_jac,nincmax,nincmax_jac,strain_increment
+ save implicit_hard
+
+c Variables for paralelization
+ integer noelprocess(10)
+ logical init(10)
+ logical init1(10)
+ save init
+ save init1
+ save noelprocess
+
+C STEP 0: READ CRYSTAL PROPERTIES AND COMPUTE AND SAVE SOME VARIABLES
+C Reads 'crystal.prop' and saves some static variables
+
+c This STEP has to be run ONLY ONCE AT BEGINNING (t=0) for the whole model (no loop on elements or GP)
+C 'init1' and 'init' variable controls if the STEP has been run, if init=.true. the block is jumped
+c To ensure compatibility for mpi and thread paralelization, paralelization subroutine is called
+
+ call paralelization(init,init1,noel,numprocess,noelprocess)
+
+ if(.not. init1(numprocess+1).and.noel.eq.
+ 1 noelprocess(numprocess+1))
+ 1 THEN
+
+c Some constants
+ pi=4d0*datan(1d0)
+
+c OPEN THE FILE WITH THE CRYSTAL DEFINITION
+c CALL GETOUTDIR( OUTDIR, LENOUTDIR )
+c FILENAME = OUTDIR(1:LENOUTDIR)//'/crystal.prop' ! Opens crystal.prop
+ paths="/"
+c CALL GETCWD(OUTDIR)
+c FILENAME = TRIM(OUTDIR)// paths // 'MAT_MODELS'//
+c 1 paths //'umat_cp_phen_iso' // paths
+c 1 //'crystal.prop' ! Opens crystal.prop
+c CALL GETCWD(OUTDIR)
+ FILENAME = 'crystal.prop' ! Opens crystal.prop
+ UR0=74
+ OPEN(UR0,FILE=FILENAME,STATUS='OLD')
+
+ CALL READPROPERTIES(UR0,max_nsystems,c11,c12,c44,c13,c33,c66, ! Subroutine readproperties read crystal.prop
+ 1 gamma_0,mm,nsets,nsystems,s,m,q,tau0,taus,h0,h1,
+ 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
+c The symmetry of the crystal taken from c33,
+c if c33=0--> cubic (3 constants) if c33!=0 ortrhotropic (6 constants)
+
+ if(abs(c33).GT.1D-10) THEN
+ CALL STIFF6(c11,c12,c44,c13,c33,c66,stiff)
+ else
+ CALL STIFF4(c11,c12,c44,stiff)
+ endif
+
+ call sleep(1)
+ init1(numprocess+1)=.true.
+ call sleep(1)
+
+ endif
+
+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
+
+ do,ii=1,nsystems
+ do,i=1,3
+ aux3_1(i)=s(ii,i)
+ aux3_2(i)=m(ii,i)
+ enddo
+
+c Orientation of schmidt tensors to material orientation, schmidt
+
+ CALL PTENS(aux33_2,aux3_1,aux3_2)
+ CALL ROTATE_TENS3(aux33,aux33_2,ORIENT_MAT_tr) ! Schmidt_cart = Orient_MAT Schmidt_Crys Orient_MAT_T
+
+ do,i=1,3
+ do,j=1,3
+ schmidt(ii,i,j)=aux33(i,j)
+ enddo
+ enddo
+
+ enddo
+
+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
+
+ call dzeros(DFGRD_ELAS,3,3)
+ call dzeros(DFGRD_ELAS_t,3,3)
+ call dzeros(L_p,3,3)
+
+ do,ii=1,3
+ DFGRD_ELAS(ii,ii)=1d0
+ DFGRD_ELAS_t(ii,ii)=1d0
+ enddo
+
+ call dzeros(gamma_act,max_nsystems,1)
+ call dzeros(tau_act,max_nsystems,1)
+
+ do,ii=1,nsystems
+ tau_act(ii)=tau0(ii)
+ enddo
+
+ else
+
+ index=0
+ do,ii=1,3
+ do,jj=1,3
+ index=index+1
+ DFGRD_ELAS_t(ii,jj)=statev(index)
+ enddo
+ enddo
+ do,i=1,nsystems
+ index=index+1
+ gamma_act(i)=statev(index)
+ enddo
+ do,i=1,nsystems
+ index=index+1
+ tau_act(i)=statev(index)
+ enddo
+ do,ii=1,3
+ do,jj=1,3
+ index=index+1
+ L_p(ii,jj)=statev(index)
+ enddo
+ enddo
+ endif
+c write(*,*)'End reading state variables'
+c PREDICTOR:
+
+c Deformation gradient increment
+
+ CALL MATINV3(DFGRD0_inv,DFGRD0,iflag)
+
+ 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)
+ cALL matinv3(DFGRD_ELAS_PRED0_inv,DFGRD_ELAS_PRED0,iflag)
+
+c Jacobian initialization using DFGRD_ELAS_PRED0_inv
+c write(*,*)'dtime ',dtime
+ CALL form_Mjacob(dtime,DFGRD_ELAS_PRED0_inv,Mjacob)
+
+c Elastic deformation gradient prediction based on last L_p
+
+ do,ii=1,3
+ do,jj=1,3
+ L_p_new(ii,jj)=L_p(ii,jj)*(-dtime)
+ enddo
+ enddo
+
+c CALL EXPMAP(EXP_L_P,noconv,L_P)
+
+ do,ii=1,3
+ do,jj=1,3
+ EXP_L_P(ii,jj)=kroneker(ii,jj)+L_p_new(ii,jj)
+ enddo
+ enddo
+
+ CALL PMAT(DFGRD_ELAS_PRED,dfgrd_elas_pred0,exp_l_p,3,3,3)
+
+
+c gamma and CRSS initialization
+
+ gamma_TOT_pred=0d0
+ call dzeros(dgamma_new,max_nsystems,1)
+ call dzeros(dgamma,max_nsystems,1)
+
+ do,isystem=1,nsystems
+
+ tau_pred(isystem)=tau_act(isystem)
+ gamma_pred(isystem)=gamma_act(isystem)
+ gamma_TOT_pred=gamma_TOT_pred+gamma_pred(isystem)
+ gamma_TOT_act=gamma_TOT_act+gamma_act(isystem)
+
+c dgamma_prediction = 0
+
+ dgamma(isystem)=0.
+
+ enddo
+
+c GLOBAL NEWTON-RAPHSON
+
+c write(*,*)'Start Newton-Raphson'
+
+ dnorm_NR=1d20
+ iter=0
+ toolarge=0
+
+ DO 100 WHILE(dnorm_NR.GT.TOLER)
+
+ iter=iter+1
+
+C 1.1: Decomposition of Fe in U and R
+c 1.3: Lagrangian deformation
+
+c Lagrangian deformation
+ CALL GREEN_LAGRANGE(DFGRD_ELAS_pred,EPSILON_EL)
+
+c 1.4: Calculation of Kirchoff stress
+
+c skirchoff_rot DEFINED ON crystal axes (undeformed configuration)
+c skirchoff DEFINED ON final configuration
+c if(noel.eq.1)then
+c do,ii=1,3
+c do,jj=1,3
+c do,kk=1,3
+c do,ll=1,3
+c write(*,*)ii,jj,kk,ll,STIFF_ORIENT(ii,jj,kk,ll)
+c enddo
+c enddo
+c enddo
+c enddo
+c endif
+ CALL PCONTRACT2(skirchoff_rot,STIFF_ORIENT,epsilon_el,3)
+
+c 2: PROYECTION OF STRESS IN SYSTEMS
+
+c 2.2: Obtention of resolved tau on all systems
+c write(*,*)skirchoff_rot
+ do,isystem=1,nsystems
+ do,ii=1,3
+ do,jj=1,3
+ AUX33_2(ii,jj)=schmidt(isystem,ii,jj)
+ enddo
+ enddo
+
+c write(*,*)AUX33_2
+ CALL PCONTRACT(tau_resolved(isystem),skirchoff_rot,AUX33_2,3)
+c write(*,*)tau_resolved(isystem)
+
+ enddo
+
+c 3: CALCULATION OF SLIP RATES
+
+
+cc NUEVO ALGORITMO PARA tau_pred
+
+ do,isystem=1,nsystems
+
+ CALL VISCOLAW(gamma_dot(isystem),gamma_0,mm,
+ 1 tau_resolved(isystem),tau_pred(isystem),noconv)
+CC write(*,*) gamma_dot(isystem),gamma_0,mm,
+CC 1 tau_resolved(isystem),tau_pred(isystem),noconv
+
+ if(noconv) THEN
+c write(*,*)'ERROR IN GAMMADOT',isystem
+ pnewdt=.75
+ return
+
+ endif
+
+ dgamma_new(isystem)=gamma_dot(isystem)*dtime
+
+ enddo
+
+
+c 4: CALCULATION OF FP
+
+C 4.1 Form L_P
+
+ call dzeros(L_p_new,3,3)
+
+ do,isystem=1,nsystems
+
+ do,ii=1,3
+ do,jj=1,3
+ L_p_new(ii,jj)=L_p_new(ii,jj)+gamma_dot(isystem)*
+ 1 schmidt(isystem,ii,jj)
+ enddo
+ enddo
+
+ enddo
+
+c 5: CALCULATION OF RESIDUAL: RES=DFGRD_ELAS_PRED-DFGRD_ELAS
+
+
+ do ii=1,3
+ do jj=1,3
+ Lterm=0d0
+ do pp=1,3
+ Lterm=Lterm+dfgrd_elas_pred0_inv(ii,pp)*
+ 1 dfgrd_elas_pred(pp,jj)
+ enddo
+ L_p(ii,jj)=(-Lterm+kroneker(ii,jj))/dtime
+ enddo
+ enddo
+
+c Forms RESIDUAL VECTOR. If double residual is used keep as it is written
+
+ call dzeros(BIGRES,9+max_nsystems,1)
+c write(*,*) L_p,L_p_new,dgamma,dgamma_new,nsystems
+ CALL form_BIGRES(L_p,L_p_new,dgamma,dgamma_new,nsystems,BIGRES)
+c write(*,*)BIGRES
+
+ CALL norm(dnorm_NR,BIGRES,9+nsystems,1)
+
+c write(*,*)'kinc,iter,dnorm',kinc,iter,dnorm_NR
+
+ if(dnorm_NR.LT.TOLER) THEN
+ dnorm=0d0
+ goto 201
+ endif
+
+ if(iter.GE.nincmax) THEN
+c write(*,*)'ERROR, too many iterations'
+ PNEWDT=0.75
+ return
+ endif
+
+ if(dnorm_NR.GT.1d10.and.toolarge.LE.1) THEN
+ toolarge=2
+ iter=0
+c write(*,*)'toolarge RES',noel,npt,kinc
+
+ dnorm=1d50
+ goto 201
+ endif
+
+
+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
+
+ CALL form_partial_sigma_F(DFGRD_ELAS_PRED,stiff_orient,
+ 1 partial_sigma)
+
+ do,isystem=1,nsystems
+ do,ii=1,3
+ do,jj=1,3
+ partial_tau(isystem,ii,jj)=0d0
+ enddo
+ enddo
+ enddo
+
+ do,isystem=1,nsystems
+
+c dfactor=\partial gamma_dot \partial tau_resolved
+
+ dfactor(isystem)=(1/mm)*gamma_0*
+ 1 (abs(tau_resolved(isystem)/tau_pred(isystem)))
+ 1 **((1/mm)-1)
+ 1 *(1d0/tau_pred(isystem))
+
+c dfactor2=\partial gamma_dot \partial tau_pred
+
+ if(tau_resolved(isystem).GE.0d0) then
+ signoT=1.
+ else
+ signoT=-1
+ endif
+
+ dfactor2(isystem)=dfactor(isystem)*
+ 1 abs(tau_resolved(isystem)/tau_pred(isystem))*
+ 1 signoT
+
+ do, nm=1,3
+ do, nn=1,3
+ do,ii=1,3
+ do,jj=1,3
+ partial_tau(isystem,nm,nn)=
+ 1 partial_tau(isystem,nm,nn)+
+ 1 dfactor(isystem)*partial_sigma(ii,jj,nm,nn)*
+ 1 schmidt(isystem,ii,jj)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+ CALL dzeros_tensor4(jacob_NR,3,3,3,3)
+
+ do,ii=1,3
+ do,jj=1,3
+ do,nm=1,3
+ do,nn=1,3
+
+ term=0d0
+
+ do, isystem=1, nsystems ! sum on isystem
+c
+ term=term+partial_tau(isystem,nm,nn)*
+ 1 schmidt(isystem,ii,jj)
+ enddo ! end isystem
+
+ jacob_NR(ii,jj,nm,nn)=
+ 1 Mjacob(ii,jj,nm,nn)-term
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call dzeros(BIGJAC,9+max_nsystems,9+max_nsystems)
+ call dzeros(BIGCORR,9+max_nsystems,1)
+
+c write(*,*)'BOX11: partial Lp / partial Fe'
+
+ CALL TENS3333(jacob_NR,R1_FE)
+ do,ii=1,9
+ do jj=1,9
+ BIGJAC(ii,jj)=R1_FE(ii,jj)
+ enddo
+ enddo
+
+c write(*,*)'BOX12: Partial Lp/ partial dgamma'
+
+ do,jsystem=1,nsystems
+c
+ call dzeros(aux33,3,3)
+ do isystem=1,nsystems
+
+ if(dgamma(jsystem).GE.0) THEN
+ signog=1.
+ else
+ signog=-1.
+ endif
+
+ aux_escalar=dfactor2(isystem)*
+ 1 q(isystem,jsystem)*
+ 1 h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*signog
+ do,ii=1,3
+ do,jj=1,3
+ aux33(ii,jj)=aux33(ii,jj)+
+ 1 aux_escalar*schmidt(isystem,ii,jj)
+ enddo
+ enddo
+ enddo
+ BIGJAC(1,9+jsystem)=aux33(1,1)
+ BIGJAC(2,9+jsystem)=aux33(2,2)
+ BIGJAC(3,9+jsystem)=aux33(3,3)
+ BIGJAC(4,9+jsystem)=aux33(1,2)
+ BIGJAC(5,9+jsystem)=aux33(1,3)
+ BIGJAC(6,9+jsystem)=aux33(2,3)
+ BIGJAC(7,9+jsystem)=aux33(2,1)
+ BIGJAC(8,9+jsystem)=aux33(3,1)
+ BIGJAC(9,9+jsystem)=aux33(3,2)
+ enddo
+
+c write(*,*)'BOX21: Partial dgamma / partial Fe'
+
+ do isystem=1,nsystems
+
+ BIGJAC(9+isystem,1)=-dtime*partial_tau(isystem,1,1)
+ BIGJAC(9+isystem,2)=-dtime*partial_tau(isystem,2,2)
+ BIGJAC(9+isystem,3)=-dtime*partial_tau(isystem,3,3)
+ BIGJAC(9+isystem,4)=-dtime*partial_tau(isystem,1,2)
+ BIGJAC(9+isystem,5)=-dtime*partial_tau(isystem,1,3)
+ BIGJAC(9+isystem,6)=-dtime*partial_tau(isystem,2,3)
+ BIGJAC(9+isystem,7)=-dtime*partial_tau(isystem,2,1)
+ BIGJAC(9+isystem,8)=-dtime*partial_tau(isystem,3,1)
+ BIGJAC(9+isystem,9)=-dtime*partial_tau(isystem,3,2)
+
+ enddo
+
+c write(*,*)'BOX22: Partial dgamma/dgamma'
+
+ do isystem=1,nsystems
+ do jsystem=1,nsystems
+
+ if(dgamma(jsystem).GE.0d0) THEN
+ signog=1.
+ else
+ signog=-1.
+ endif
+
+ BIGJAC(9+isystem,9+jsystem)=kroneker(isystem,jsystem)
+ 1 +dfactor2(isystem)*
+ 1 h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*
+ 1 signog*q(isystem,jsystem)*dtime
+ 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)
+ DFGRD_ELAS_PRED(2,2)=DFGRD_ELAS_PRED(2,2)-BIGCORR(2)
+ DFGRD_ELAS_PRED(3,3)=DFGRD_ELAS_PRED(3,3)-BIGCORR(3)
+ DFGRD_ELAS_PRED(1,2)=DFGRD_ELAS_PRED(1,2)-BIGCORR(4)
+ DFGRD_ELAS_PRED(1,3)=DFGRD_ELAS_PRED(1,3)-BIGCORR(5)
+ DFGRD_ELAS_PRED(2,3)=DFGRD_ELAS_PRED(2,3)-BIGCORR(6)
+ DFGRD_ELAS_PRED(2,1)=DFGRD_ELAS_PRED(2,1)-BIGCORR(7)
+ DFGRD_ELAS_PRED(3,1)=DFGRD_ELAS_PRED(3,1)-BIGCORR(8)
+ DFGRD_ELAS_PRED(3,2)=DFGRD_ELAS_PRED(3,2)-BIGCORR(9)
+
+ 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))
+ gamma_tot_pred=gamma_tot_pred+gamma_pred(isystem)
+ gamma_tot_act=gamma_tot_act+gamma_act(isystem)
+ enddo
+
+ do,isystem=1,nsystems
+ tau_pred(isystem)=tau_act(isystem)
+
+ do,jsystem=1,nsystems
+ tau_pred(isystem)=tau_pred(isystem)+
+ 1 q(isystem,jsystem)
+ 1 *h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))
+ enddo
+
+ enddo
+c
+c$$$ do,isystem=1,nsystems
+c$$$ tau_pred(isystem)=tau_act(isystem)
+c$$$ if(abs(h(gamma_tot_pred,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))-
+c$$$ 1 h(gamma_tot_act,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem)))/h0(isystem).GT..1) THEN
+c$$$ write(*,*)'subincrementation',kinc
+c$$$ do,jsystem=1,nsystems
+c$$$ do,ii=1,11
+c$$$ gamma_tot=gamma_tot_act+gamma_tot_pred*(ii-1.)/11.
+c$$$ tau_pred(isystem)=tau_pred(isystem)+
+c$$$ 1 q(isystem,jsystem)
+c$$$ 1 *h(gamma_tot,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))*.1
+c$$$ enddo
+c$$$ enddo
+c$$$ else
+c$$$ do,jsystem=1,nsystems
+c$$$ tau_pred(isystem)=tau_pred(isystem)+
+c$$$ 1 q(isystem,jsystem)
+c$$$ 1 *h(gamma_tot_pred,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))
+c$$$ enddo
+c$$$ endif
+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
+
+c write(*,*)'201b'
+
+ call polar_decomp(aux33,rr_tot,aux33_2,DFGRD_inc)
+ do,ii=1,3
+ do,jj=1,3
+ DFGRD_ELAS_PRED(ii,jj)=0d0
+ do,kk=1,3
+ DFGRD_ELAS_PRED(ii,jj)=DFGRD_ELAS_PRED(ii,jj)+
+ 1 rr_tot(ii,kk)*DFGRD_ELAS_t(kk,jj)
+
+ enddo
+ enddo
+ enddo
+
+ CALL GREEN_LAGRANGE(DFGRD_ELAS_pred,EPSILON_EL)
+
+ CALL PCONTRACT2(skirchoff_rot,STIFF_ORIENT,epsilon_el,3)
+ do,isystem=1,nsystems
+ do,ii=1,3
+ do,jj=1,3
+ AUX33_2(ii,jj)=schmidt(isystem,ii,jj)
+ enddo
+ enddo
+ CALL PCONTRACT(tau_resolved(isystem),skirchoff_rot,
+ 1 AUX33_2,3)
+ enddo
+
+ gamma_tot_pred=0d0
+ gamma_tot_act=0d0
+ do,isystem=1,nsystems
+
+ CALL VISCOLAW(gamma_dot(isystem),gamma_0,mm,
+ 1 tau_resolved(isystem),tau_pred(isystem),noconv)
+
+ if(noconv) THEN
+c write(*,*)'ERROR IN GAMMADOT',isystem
+ pnewdt=.75
+ return
+ endif
+
+ dgamma(isystem)=gamma_dot(isystem)*dtime
+ gamma_pred(isystem)=gamma_act(isystem)
+ 1 +abs(dgamma(isystem))
+ gamma_tot_pred=gamma_tot_pred+gamma_pred(isystem)
+ gamma_tot_act=gamma_tot_pred+gamma_act(isystem)
+
+ enddo
+
+ do,isystem=1,nsystems
+ tau_pred(isystem)=tau_act(isystem)
+
+ do,jsystem=1,nsystems
+ tau_pred(isystem)=tau_pred(isystem)+
+ 1 q(isystem,jsystem)
+ 1 *h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))
+ enddo
+
+ enddo
+c
+c$$$ do,isystem=1,nsystems
+c$$$ tau_pred(isystem)=tau_act(isystem)
+c$$$
+c$$$ if(abs(h(gamma_tot_pred,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))-
+c$$$ 1 h(gamma_tot_act,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem)))/h0(isystem).GT..1) THEN
+c$$$ write(*,*)'subincrementation'
+c$$$ do,jsystem=1,nsystems
+c$$$ do,ii=1,11
+c$$$ gamma_tot=gamma_tot_act+gamma_tot_pred*(ii-1.)/11.
+c$$$ tau_pred(isystem)=tau_pred(isystem)+
+c$$$ 1 q(isystem,jsystem)
+c$$$ 1 *h(gamma_tot,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))*
+c$$$ 1 abs(dgamma(jsystem))*.1
+c$$$ enddo
+c$$$ enddo
+c$$$ else
+c$$$ do,jsystem=1,nsystems
+c$$$ tau_pred(isystem)=tau_pred(isystem)+
+c$$$ 1 q(isystem,jsystem)
+c$$$ 1 *h(gamma_tot_pred,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))
+c$$$ enddo
+c$$$ endif
+c$$$ enddo
+ endif
+
+ 100 CONTINUE
+
+c write(*,*) 'ROTATE STRESSES TO DEFORMED CONFIGURATION'
+
+ CALL POLAR_DECOMP(UU,RR,CC,DFGRD_ELAS_pred)
+c CALL TRANSPOSE(RR_tr,RR,3,3)
+ CALL TRANSPOSE(DFGRD_ELAS_pred_tr,DFGRD_ELAS_pred,3,3)
+ CALL ROTATE_TENS3(skirchoff,skirchoff_ROT,DFGRD_ELAS_pred_tr)
+ Jdet=DFGRD_ELAS_pred(1,1)*DFGRD_ELAS_pred(2,2)*
+ 1 DFGRD_ELAS_pred(3,3)+
+ 1 DFGRD_ELAS_pred(1,2)*DFGRD_ELAS_pred(2,3)*DFGRD_ELAS_pred(3,1)+
+ 1 DFGRD_ELAS_pred(1,3)*DFGRD_ELAS_pred(2,1)*DFGRD_ELAS_pred(3,2)-
+ 1 DFGRD_ELAS_pred(1,3)*DFGRD_ELAS_pred(2,2)*DFGRD_ELAS_pred(3,1)-
+ 1 DFGRD_ELAS_pred(2,3)*DFGRD_ELAS_pred(3,2)*DFGRD_ELAS_pred(1,1)-
+ 1 DFGRD_ELAS_pred(3,3)*DFGRD_ELAS_pred(1,2)*DFGRD_ELAS_pred(2,1)
+ do,ii=1,3
+ STRESS(ii)=1.0/Jdet*skirchoff(ii,ii)
+ enddo
+ 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 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),
+c 1 skirchoff_ROT(2,3)
+c write(*,*) 'kir',skirchoff_ROT(3,1),skirchoff_ROT(3,2),
+c 1 skirchoff_ROT(3,3)
+c write(*,*) 'sig',STRESS
+c write(*,*) 'fel',DFGRD_ELAS_pred(1,1),DFGRD_ELAS_pred(1,2),
+c 1 DFGRD_ELAS_pred(1,3)
+c write(*,*) 'fel',DFGRD_ELAS_pred(2,1),DFGRD_ELAS_pred(2,2),
+c 1 DFGRD_ELAS_pred(2,3)
+c write(*,*) 'fel',DFGRD_ELAS_pred(3,1),DFGRD_ELAS_pred(3,2),
+c 1 DFGRD_ELAS_pred(3,3)
+c write(*,*) 'green',epsilon_el(1,1),epsilon_el(1,2),
+c 1 epsilon_el(1,3)
+c write(*,*) 'green',epsilon_el(2,1),epsilon_el(2,2),
+c 1 epsilon_el(2,3)
+c write(*,*) 'green',epsilon_el(3,1),epsilon_el(3,2),
+c 1 epsilon_el(3,3)
+c endif
+
+C SAVE INTERNAL VARIABLES
+C AND REDEFINE INTIAL STATE FOR JACOBIAN CALCULATION
+
+ index=0
+ do,ii=1,3
+ do,jj=1,3
+ index=index+1
+ statev(index)=DFGRD_ELAS_pred(ii,jj)
+ dfgrd_elas_pred_jacob(ii,jj)=dfgrd_elas_pred(ii,jj)
+ enddo
+ enddo
+
+ do,i=1,nsystems
+ index=index+1
+
+c Actualization of gamma_act
+ statev(index)=gamma_pred(i)
+
+c Definition of gammas for Jacobian
+ gamma_act_jacob(i)=gamma_act(i)
+ gamma_pred_jacob(i)=gamma_pred(i)
+ dgamma_jacob(i)=dgamma(i)
+
+ enddo
+
+ do,i=1,nsystems
+ index=index+1
+
+c Actualization of tau_act
+ statev(index)=tau_pred(i)
+
+c Definition of tau_act_jacob and tau_pred_jacob
+ tau_act_jacob(i)=tau_act(i)
+ tau_pred_jacob(i)=tau_pred(i)
+
+ enddo
+
+ do,ii=1,3
+ do,jj=1,3
+ index=index+1
+ statev(index)=L_p(ii,jj)
+ enddo
+ enddo
+
+C SAVE THE EULER ANGLES OF EACH IP
+C ROTATED_ORIENT=RR*ORIENT_MAT
+
+ call PMAT(ROTATED_ORIENT_tr,RR,ORIENT_MAT,3,3,3)
+ call transpose(ROTATED_ORIENT,ROTATED_ORIENT_tr,3,3)
+ call euler(1,statev(index+1),statev(index+2),statev(index+3)
+ 1 ,rotated_orient)
+
+C SAVE acumulated shear strain
+ index=index+4
+ statev(index)=gamma_tot_pred
+c write(*,*)'Save accumulated shear strain'
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+c CALCULATION OF JACOBIAN
+C numerical derivation
+c From DFGRD0 to DFGRD1+delta_epsilon
+C Prediction based on actual converged solution of NR: DFGRD_ELAS_PRED
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+c SAME TENSORS AS MAIN NR:
+c
+c DFGRD0
+c DFGRD_ELAS_t
+c DFGRD_ELAS_PRED
+
+c NEW TENSORS
+c dfgrd_inc
+c DFGRD_ELAS_PRED0
+
+
+ dtime2=dtime
+
+
+C Saving initial values before each step
+
+ do,ii=1,3
+ do,jj=1,3
+ dfgrd_inc0(ii,jj)=dfgrd_inc(ii,jj)
+ enddo
+ enddo
+
+C Initialize to cero the DDSDDE in 3x3x3x3 (aux3333j)
+
+ call dzeros(DDSDDE,6,6)
+
+c General loop on the 6 perturbations to obtain Jacobian
+C On each istep a whole NR problem is solved
+
+ do,istep=1,6
+
+c write(*,*)'istep',istep
+ call dzeros(delta_eps_jac,3,3)
+
+ if(istep.LE.3) THEN
+ iii=istep
+ jjj=istep
+ else if(istep.EQ.4) THEN
+ iii=1
+ jjj=2
+ else if(istep.EQ.5) THEN
+ iii=1
+ jjj=3
+ else if(istep.EQ.6) THEN
+ iii=2
+ jjj=3
+ endif
+
+c Definition of perturbation strain_increment
+
+ if(istep.LE.3) THEN
+ delta_eps_jac(jjj,iii)= strain_increment
+ else
+ delta_eps_jac(iii,jjj)= .5*strain_increment
+ delta_eps_jac(jjj,iii)= .5*strain_increment
+ endif
+
+c Obtention of new deformation gradient increment DFGRD_INC for this perturbation
+
+c$$$ do,ii=1,3
+c$$$ do,jj=1,3
+c$$$ DFGRD_INC(ii,jj)=DFGRD_INC0(ii,jj)+delta_eps_jac(ii,jj)
+c$$$c DFGRD_ELAS_PRED(II,JJ)=DFGRD_ELAS_PRED1(II,JJ)
+c$$$ enddo
+c$$$ enddo
+c$$$ write(*,*)'DFGRD_INC_forma1',DFGRD_INC
+
+ do,ii=1,3
+ do,jj=1,3
+
+ DFGRD_INC(ii,jj) = 0.0
+ do,pp=1,3
+ DFGRD_INC(ii,jj)=DFGRD_INC(ii,jj)+
+ 1 (kroneker(ii,pp)+delta_eps_jac(ii,pp))*
+ 1 DFGRD_INC0(pp,jj)
+ enddo
+ enddo
+ enddo
+c$$$ write(*,*)'DFGRD_INC_forma2',DFGRD_INC
+
+c NEW dfgrd_elas_pred0
+
+ CALL PMAT(DFGRD_ELAS_pred0,DFGRD_INC,DFGRD_ELAS_t,3,3,3)
+ cALL matinv3(DFGRD_ELAS_PRED0_inv,DFGRD_ELAS_PRED0,iflag)
+
+c Jacobian initialization using DFGRD_ELAS_PRED0_inv
+
+ CALL form_Mjacob(dtime,DFGRD_ELAS_PRED0_inv,Mjacob)
+
+C RELOAD INTERNAL VARIABLES
+
+ index=0
+
+ do,ii=1,3
+ do,jj=1,3
+ DFGRD_ELAS_PRED(ii,jj)=DFGRD_ELAS_PRED_JACOB(ii,jj)
+ enddo
+ enddo
+
+ do,i=1,nsystems
+ gamma_act(i)=gamma_act_jacob(i)
+ gamma_pred(i)=gamma_pred_jacob(i)
+ dgamma(i)=dgamma_jacob(i)
+ tau_act(i)=tau_act_jacob(i)
+ tau_pred(i)=tau_pred_jacob(i)
+ enddo
+
+
+c HERE ALL THE NR LOOP
+c Enter with DFGRD_INC and DFGRD_ELAS_t(ii,jj)
+
+c GLOBAL NEWTON-RAPHSON
+
+ dnorm_NR=1d0
+ iter=0
+ toolarge=0
+c
+ DO 101 WHILE(dnorm_NR.GT.TOLER_jac.and.iter.LE.nincmax_jac)
+
+ iter=iter+1
+
+
+c 1. COMPUTE LAGRANGIAN ELASTIC STRAIN
+C 1.1: Decomposition of Fe in U and R,
+
+
+ CALL green_lagrange(DFGRD_ELAS_pred,EPSILON_EL)
+
+ CALL PCONTRACT2(skirchoff_rot,STIFF_ORIENT,epsilon_el,3)
+
+c 2: PROYECTION OF STRESS IN SYSTEMS
+
+ do,isystem=1,nsystems
+ do,ii=1,3
+ do,jj=1,3
+ AUX33_2(ii,jj)=schmidt(isystem,ii,jj)
+ enddo
+ enddo
+ CALL PCONTRACT(tau_resolved(isystem),skirchoff_rot,
+ 1 AUX33_2,3)
+ enddo
+
+
+
+c 3: CALCULATION OF SLIP RATES
+
+ do,isystem=1,nsystems
+
+ CALL VISCOLAW(gamma_dot(isystem),gamma_0,mm,
+ 1 tau_resolved(isystem),tau_pred(isystem),noconv)
+ if(noconv) THEN
+c write(*,*)'ERROR IN GAMMADOT',isystem
+ pnewdt=.75
+ return
+
+ endif
+
+ dgamma_new(isystem)=gamma_dot(isystem)*dtime
+
+ enddo
+
+
+c 4: CALCULATION OF FP
+
+C 4.1 Form L_P
+
+ call dzeros(L_p_new,3,3)
+
+ do,isystem=1,nsystems
+
+ do,ii=1,3
+ do,jj=1,3
+ L_p_new(ii,jj)=L_p_new(ii,jj)+gamma_dot(isystem)*
+ 1 schmidt(isystem,ii,jj)
+ enddo
+ enddo
+
+ enddo
+
+c 5: CALCULATION OF RESIDUAL: RES=DFGRD_ELAS_PRED-DFGRD_ELAS
+
+
+
+ do ii=1,3
+ do jj=1,3
+ Lterm=0d0
+ do pp=1,3
+ Lterm=Lterm+dfgrd_elas_pred0_inv(ii,pp)*
+ 1 dfgrd_elas_pred(pp,jj)
+ enddo
+ L_p(ii,jj)=(-Lterm+kroneker(ii,jj))/dtime
+ enddo
+ enddo
+
+ call dzeros(BIGRES,9+max_nsystems,1)
+ CALL form_BIGRES(L_p,L_p_new,dgamma,dgamma_new,nsystems,
+ 1 BIGRES)
+
+ CALL norm(dnorm_NR,BIGRES,9+nsystems,1)
+
+c write(*,*)'JAC: kinc,iter,dnorm',kinc,iter,dnorm_NR
+ if(dnorm_NR.LT.TOLER_jac) then
+ dnorm=0d0
+ GOTO 202
+ endif
+
+c if(iter.GE.nincmax_jac-1) THEN
+c write(*,*)'ERROR_JAC, demasiadas iteraciones'
+c write(*,*) 'iter,dnorm_NR',iter,dnorm_NR
+c endif
+
+ if(dnorm_NR.GT.1E10.and.toolarge.LE.1) THEN
+ toolarge=2
+ iter=0
+ dnorm=1d50
+ goto 202
+
+ endif
+
+
+
+c 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
+
+ CALL form_partial_sigma_F(DFGRD_ELAS_PRED,stiff_orient,
+ 1 partial_sigma)
+
+ do,isystem=1,nsystems
+ do,ii=1,3
+ do,jj=1,3
+ partial_tau(isystem,ii,jj)=0d0
+ enddo
+ enddo
+ enddo
+
+ do,isystem=1,nsystems
+
+ dfactor(isystem)=(1/mm)*gamma_0*
+ 1 (abs(tau_resolved(isystem)/tau_pred(isystem)))
+ 1 **((1/mm)-1)*(1d0/tau_pred(isystem))
+
+c dfactor2=\partial gamma_dot \partial tau_pred
+
+ if(tau_resolved(isystem).GE.0d0) then
+ signoT=1.
+ else
+ signoT=-1
+ endif
+
+ dfactor2(isystem)=dfactor(isystem)*
+ 1 abs(tau_resolved(isystem)/tau_pred(isystem))*
+ 1 signoT
+
+ do, nm=1,3
+ do, nn=1,3
+ do,ii=1,3
+ do,jj=1,3
+ partial_tau(isystem,nm,nn)=
+ 1 partial_tau(isystem,nm,nn)+
+ 1 dfactor(isystem)*
+ 1 partial_sigma(ii,jj,nm,nn)*
+ 1 schmidt(isystem,ii,jj)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+ CALL dzeros_tensor4(jacob_NR,3,3,3,3)
+
+
+ do,ii=1,3
+ do,jj=1,3
+ do,nm=1,3
+ do,nn=1,3
+
+ term=0d0
+
+ do, isystem=1, nsystems ! sum on isystem
+c
+ term=term+partial_tau(isystem,nm,nn)*
+ 1 schmidt(isystem,ii,jj)
+ enddo ! end isystem
+
+ jacob_NR(ii,jj,nm,nn)=
+ 1 Mjacob(ii,jj,nm,nn)-term
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call dzeros(BIGJAC,9+max_nsystems,9+max_nsystems)
+ call dzeros(BIGCORR,9+max_nsystems,1)
+
+c BOX11: partial Lp / partial Fe
+
+ CALL TENS3333(jacob_NR,R1_FE)
+ do,ii=1,9
+ do jj=1,9
+ BIGJAC(ii,jj)=R1_FE(ii,jj)
+ enddo
+ enddo
+
+c BOX12: Partial Lp/ partial dgamma
+
+ do,jsystem=1,nsystems
+c
+ call dzeros(aux33,3,3)
+ do isystem=1,nsystems
+
+ if(dgamma(jsystem).GE.0) THEN
+ signog=1.
+ else
+ signog=-1.
+ endif
+
+ aux_escalar=dfactor2(isystem)*
+ 1 q(isystem,jsystem)*
+ 1 h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*signog
+
+ do,ii=1,3
+ do,jj=1,3
+ aux33(ii,jj)=aux33(ii,jj)+
+ 1 aux_escalar*schmidt(isystem,ii,jj)
+ enddo
+ enddo
+ enddo
+
+ BIGJAC(1,9+jsystem)=aux33(1,1)
+ BIGJAC(2,9+jsystem)=aux33(2,2)
+ BIGJAC(3,9+jsystem)=aux33(3,3)
+ BIGJAC(4,9+jsystem)=aux33(1,2)
+ BIGJAC(5,9+jsystem)=aux33(1,3)
+ BIGJAC(6,9+jsystem)=aux33(2,3)
+ BIGJAC(7,9+jsystem)=aux33(2,1)
+ BIGJAC(8,9+jsystem)=aux33(3,1)
+ BIGJAC(9,9+jsystem)=aux33(3,2)
+ enddo
+
+c BOX21: Partial dgamma / partial Fe
+ do isystem=1,nsystems
+ BIGJAC(9+isystem,1)=-dtime*partial_tau(isystem,1,1)
+ BIGJAC(9+isystem,2)=-dtime*partial_tau(isystem,2,2)
+ BIGJAC(9+isystem,3)=-dtime*partial_tau(isystem,3,3)
+ BIGJAC(9+isystem,4)=-dtime*partial_tau(isystem,1,2)
+ BIGJAC(9+isystem,5)=-dtime*partial_tau(isystem,1,3)
+ BIGJAC(9+isystem,6)=-dtime*partial_tau(isystem,2,3)
+ BIGJAC(9+isystem,7)=-dtime*partial_tau(isystem,2,1)
+ BIGJAC(9+isystem,8)=-dtime*partial_tau(isystem,3,1)
+ BIGJAC(9+isystem,9)=-dtime*partial_tau(isystem,3,2)
+
+ enddo
+
+c BOX22: Partial dgamma/dgamma
+ do isystem=1,nsystems
+ do jsystem=1,nsystems
+
+ if(dgamma(jsystem).GE.0) THEN
+ signog=1.
+ else
+ signog=-1.
+ endif
+ BIGJAC(9+isystem,9+jsystem)=
+ 1 kroneker(isystem,jsystem)
+ 1 +dfactor2(isystem)*
+ 1 h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*
+ 1 signog*q(isystem,jsystem)*dtime
+ enddo
+ enddo
+
+ ntot=max_nsystems+9
+ nloc=9+nsystems
+ call gauss_3(BIGJAC,BIGRES,BIGCORR,ntot,nloc,iflag)
+ if(iflag.EQ.1) THEN
+c write(*,*)'error in system of eqs'
+ PNEWDT=.5
+ return
+ endif
+
+ DFGRD_ELAS_PRED(1,1)=DFGRD_ELAS_PRED(1,1)-BIGCORR(1)
+ DFGRD_ELAS_PRED(2,2)=DFGRD_ELAS_PRED(2,2)-BIGCORR(2)
+ DFGRD_ELAS_PRED(3,3)=DFGRD_ELAS_PRED(3,3)-BIGCORR(3)
+ DFGRD_ELAS_PRED(1,2)=DFGRD_ELAS_PRED(1,2)-BIGCORR(4)
+ DFGRD_ELAS_PRED(1,3)=DFGRD_ELAS_PRED(1,3)-BIGCORR(5)
+ DFGRD_ELAS_PRED(2,3)=DFGRD_ELAS_PRED(2,3)-BIGCORR(6)
+ DFGRD_ELAS_PRED(2,1)=DFGRD_ELAS_PRED(2,1)-BIGCORR(7)
+ DFGRD_ELAS_PRED(3,1)=DFGRD_ELAS_PRED(3,1)-BIGCORR(8)
+ DFGRD_ELAS_PRED(3,2)=DFGRD_ELAS_PRED(3,2)-BIGCORR(9)
+
+ gamma_tot_pred=0d0
+ gamma_tot_act=0d0
+ do,isystem=1,nsystems
+ dgamma(isystem)=dgamma(isystem)-BIGCORR(9+isystem)
+ gamma_pred(isystem)=gamma_act(isystem)+
+ 1 abs(dgamma(isystem))
+ gamma_tot_pred=gamma_tot_pred+gamma_pred(isystem)
+ gamma_tot_act=gamma_tot_act+gamma_act(isystem)
+ enddo
+
+ do,isystem=1,nsystems
+ tau_pred(isystem)=tau_act(isystem)
+
+ do,jsystem=1,nsystems
+ tau_pred(isystem)=tau_pred(isystem)+
+ 1 q(isystem,jsystem)
+ 1 *h(gamma_tot_pred,tau0(isystem),taus(isystem),
+ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))
+ enddo
+
+ enddo
+c
+c$$$ do,isystem=1,nsystems
+c$$$ tau_pred(isystem)=tau_act(isystem)
+c$$$ if(abs(h(gamma_tot_pred,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))-
+c$$$ 1 h(gamma_tot_act,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem)))/h0(isystem).GT..1) THEN
+c$$$ write(*,*)'subincrementation'
+c$$$ do,jsystem=1,nsystems
+c$$$ do,ii=1,11
+c$$$ gamma_tot=
+c$$$ 1 gamma_tot_act+gamma_tot_pred*(ii-1.)/11.
+c$$$ tau_pred(isystem)=tau_pred(isystem)+
+c$$$ 1 q(isystem,jsystem)
+c$$$ 1 *h(gamma_tot,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))*
+c$$$ 1 abs(dgamma(jsystem))*.1
+c$$$ enddo
+c$$$ enddo
+c$$$ else
+c$$$ do,jsystem=1,nsystems
+c$$$ tau_pred(isystem)=tau_pred(isystem)+
+c$$$ 1 q(isystem,jsystem)
+c$$$ 1 *h(gamma_tot_pred,tau0(isystem),taus(isystem),
+c$$$ 1 h0(isystem),h1(isystem))*abs(dgamma(jsystem))
+c$$$ enddo
+c$$$ endif
+c$$$ enddo
+
+
+c
+c NUMERICAL PROBLEMS WITH LARGE CORRECTIONS-->PLASTIC PREDICTOR
+c
+
+ call norm(dnorm,BIGCORR,max_nsystems+9,1)
+
+ 202 IF(dnorm.GT.1d10) THEN
+
+c write(*,*)'In jacob in 202'
+
+ call polar_decomp(aux33,rr_tot,aux33_2,DFGRD_inc)
+ do,ii=1,3
+ do,jj=1,3
+ DFGRD_ELAS_PRED(ii,jj)=0d0
+ do,kk=1,3
+ DFGRD_ELAS_PRED(ii,jj)=DFGRD_ELAS_PRED(ii,jj)+
+ 1 rr_tot(ii,kk)*DFGRD_ELAS_t(kk,jj)
+
+ enddo
+ enddo
+ enddo
+
+ do,isystem=1,nsystems
+ gamma_pred(isystem)=gamma_act(isystem)
+ dgamma(isystem)=0.
+ enddo
+ do,isystem=1,nsystems
+ tau_pred(isystem)=tau_act(isystem)
+ enddo
+
+ endif
+
+
+ 101 CONTINUE
+
+c ROTATE STRESSES TO DEFORMED CONFIGURATION
+
+c CALL POLAR_DECOMP(UU,RR,CC,DFGRD_ELAS_pred)
+c CALL TRANSPOSE(RR_tr,RR,3,3)
+ CALL TRANSPOSE(DFGRD_ELAS_pred_tr,DFGRD_ELAS_pred,3,3)
+ CALL ROTATE_TENS3(skirchoff2,skirchoff_rot,DFGRD_ELAS_pred_tr)
+ Jdet2=DFGRD_ELAS_pred(1,1)*DFGRD_ELAS_pred(2,2)*
+ 1 DFGRD_ELAS_pred(3,3)+
+ 1 DFGRD_ELAS_pred(1,2)*DFGRD_ELAS_pred(2,3)*DFGRD_ELAS_pred(3,1)+
+ 1 DFGRD_ELAS_pred(1,3)*DFGRD_ELAS_pred(2,1)*DFGRD_ELAS_pred(3,2)-
+ 1 DFGRD_ELAS_pred(1,3)*DFGRD_ELAS_pred(2,2)*DFGRD_ELAS_pred(3,1)-
+ 1 DFGRD_ELAS_pred(2,3)*DFGRD_ELAS_pred(3,2)*DFGRD_ELAS_pred(1,1)-
+ 1 DFGRD_ELAS_pred(3,3)*DFGRD_ELAS_pred(1,2)*DFGRD_ELAS_pred(2,1)
+c Numerical Jacobian is J=skirchoff2-skirchoff/strain_increment
+c Directly computed in Nye notation
+
+ DDSDDE(istep,1)= (1/Jdet2*skirchoff2(1,1)-1/
+ 1 Jdet*skirchoff(1,1))/
+ 1 strain_increment
+ DDSDDE(istep,2)= (1/Jdet2*skirchoff2(2,2)-1/
+ 1 Jdet*skirchoff(2,2))/
+ 1 strain_increment
+ DDSDDE(istep,3)= (1/Jdet2*skirchoff2(3,3)-1/
+ 1 Jdet*skirchoff(3,3))/
+ 1 strain_increment
+ DDSDDE(istep,4)= (1/Jdet2*skirchoff2(1,2)-1/
+ 1 Jdet*skirchoff(1,2))/
+ 1 strain_increment
+ DDSDDE(istep,5)= (1/Jdet2*skirchoff2(1,3)-1/
+ 1 Jdet*skirchoff(1,3))/
+ 1 strain_increment
+ DDSDDE(istep,6)= (1/Jdet2*skirchoff2(2,3)-1/
+ 1 Jdet*skirchoff(2,3))/
+ 1 strain_increment
+
+ enddo ! END OF 6 perturbation steps
+c if((noel.eq.1010101).and.(npt.eq.1))then
+c write(*,*) 'DDSDDE'
+c do i=1,6
+c write(*,*) (DDSDDE(i,j), j=1,6)
+c enddo
+
+c write(*,*) 'Defo0,inc',DFGRD0,kinc
+c write(*,*) 'Defo1,inc',DFGRD1,kinc
+c endif
+c write(*,*) DFGRD0,DFGRD1
+c CALL TENS3333_sym(STIFF_orient,ddsdde)
+ 110 IF(PNEWDT.LT.1) THEN
+ CALL TENS3333_sym(STIFF_orient,ddsdde)
+ ENDIF
+
+ RETURN
+ END
+
+c The library cointatining general purpose subroutines is included
+
+c include 'library_crysplas_CAPSUL_v0.f'
+
+
+c This is the set of subroutines dependent of the constitutive equations
+
+ SUBROUTINE READPROPERTIES(UR0,max_nsystems,c11,c12,c44,c13,c33
+ 1 ,c66,gamma_0,
+ 1 mm,nsets,nsystems,s,m,q,tau0,taus,h0,h1,
+ 1 toler,toler_jac,nincmax,nincmax_jac,strain_increment,
+ 1 implicit_hard)
+
+c Reads input file with crystal properties and orientations
+c Generates data on each system
+c Normalizes the vectors
+ implicit none
+ integer ii,jj,i,j,k,max_nsystems
+ real*8 c11,c12,c44,c13,c33,c66
+ real*8 gamma_0,mm
+ integer nsets,nsystems,ur0
+ real*8 q(max_nsystems,max_nsystems)
+ real*8 qsys(6,6)
+ real*8 h0(*),h1(*),tau0(*),taus(*)
+ real*8 h0_set(6),tau0_set(6),taus_set(6),h1_set(6)
+ real*8 s(max_nsystems,3),m(max_nsystems,3)
+ integer system_set(max_nsystems)
+
+ real*8 dnorm_m,dnorm_s
+ real*8 toler,toler_jac,strain_increment
+ integer nincmax,nincmax_jac,implicit_hard
+
+ character*100 line
+c Reading all the crystal data:
+
+ READ(UR0,*)
+ READ(UR0,*)
+ READ(UR0,111) line
+ 111 FORMAT(A100)
+ c13=0d0
+ c33=0d0
+ c66=0d0
+ READ(line,*,END=112) c11,c12,c44,c13,c33,c66
+ 112 write(*,*) 'c11 a c66',c11,c12,c44,c13,c33,c66
+ READ(UR0,*)
+ READ(UR0,*) gamma_0,mm
+
+ READ(UR0,*)
+ READ(UR0,*) nsets,nsystems
+
+ READ(UR0,*)
+ do,ii=1,nsystems
+ READ(UR0,*) m(ii,1), m(ii,2), m(ii,3),s(ii,1),
+ 1 s(ii,2),s(ii,3),
+ 1 system_set(ii)
+ dnorm_m=dsqrt(m(ii,1)*m(ii,1)+m(ii,2)*m(ii,2)
+ 1 +m(ii,3)*m(ii,3))
+ dnorm_s=dsqrt(s(ii,1)*s(ii,1)+s(ii,2)*s(ii,2)
+ 1 +s(ii,3)*s(ii,3))
+ do,jj=1,3
+ m(ii,jj)=m(ii,jj)/dnorm_m
+ s(ii,jj)=s(ii,jj)/dnorm_s
+ enddo
+c write(*,*) 'system',ii
+c write(*,*) m(ii,1),m(ii,2),m(ii,3)
+c write(*,*) s(ii,1),s(ii,2),s(ii,3)
+ enddo
+
+ READ(UR0,*)
+ do,ii=1,nsets
+ READ(UR0,*) (qsys(ii,jj),jj=1,nsets)
+ enddo
+
+ READ(UR0,*)
+ do,ii=1,nsets
+ h1_set(ii)=0d0
+ READ(UR0,111) line
+ READ(line,*,END=113) tau0_set(ii),taus_set(ii),h0_set(ii)
+ 1 ,h1_set(ii)
+ 113 write(*,*) 'tau0,taus,h0,h1', tau0_set(ii),taus_set(ii),
+ 1 h0_set(ii),h1_set(ii)
+ enddo
+
+ READ(UR0,*)
+ READ(UR0,*) toler,toler_jac,nincmax,nincmax_jac,
+ 1 strain_increment,
+ 1 implicit_hard
+c Creating the latent-hardening matrix for each individual system
+
+ call dzeros(q,max_nsystems,max_nsystems)
+ call dzeros(tau0,max_nsystems,1)
+ call dzeros(taus,max_nsystems,1)
+ call dzeros(h0,max_nsystems,1)
+ call dzeros(h1,max_nsystems,1)
+
+ do,ii=1,nsystems
+ do,jj=1,nsystems
+ q(ii,jj)=qsys(system_set(ii),system_set(jj))
+ if(ii.EQ.jj) q(ii,jj)=1d0
+ enddo
+c write(*,101) (q(ii,jj),jj=1,nsystems)
+ enddo
+ 101 format(20(F2.0,' '))
+
+c Creating the properties for each individual system
+ do,ii=1,nsystems
+ tau0(ii)=tau0_set(system_set(ii))
+ taus(ii)=taus_set(system_set(ii))
+ h0(ii)=h0_set(system_set(ii))
+ h1(ii)=h1_set(system_set(ii))
+ enddo
+ write(*,*)'END READING PROPERTIES (crystal.prop)'
+ RETURN
+ END
+
+c viscolaw --> viscoplastic law defining gammadot as function of
+c resolved shear stress (tauintern), and other parameters,
+c internal variables or stress components depending on the model
+
+ SUBROUTINE VISCOLAW(gammadot,gammadot0,expo,tauintern,tau_crit,
+ 1 noconv)
+
+ implicit none
+ logical noconv
+ real*8 gammadot0,expo
+ real*8 tauintern,tau_crit,gammadot
+ real*8 logtau
+ if(abs(tauintern).GT.0d0) then
+ logtau=log(abs(tauintern/tau_crit))
+ else
+ logtau=0
+ endif
+c write(*,*)'mm=',expo
+c write(*,*)'logtau=',logtau
+ if(logtau.GT.expo*150*log(10.)) THEN
+c write(*,*)'LOGTAU>X',logtau
+ gammadot=0
+ noconv=.TRUE.
+ else
+ noconv=.FALSE.
+ if(tauintern.ge.0)then
+ gammadot=gammadot0*
+ 1 abs((tauintern/tau_crit))**(1/expo)
+ else
+ gammadot=gammadot0*-1.*
+ 1 abs((tauintern/tau_crit))**(1/expo)
+ endif
+
+ endif
+ RETURN
+ END
+
+C
+ FUNCTION h(gamma,tau0I,tausI,h0I,h1I)
+C Self Hardening law: Either Asaro-Needleman or Voce
+ implicit none
+ REAL*8 h,gamma,tau0I,tausI,h0I,secah,h1I,expo
+ REAL*8 tausVOCE
+
+ if(h1I.LE.1D-20) then
+
+ h=h0I*( secah (abs(h0I*gamma/(tausI-tau0I))) )**2
+
+ else
+
+c In this case, taus is recalculated as taus-tau0
+
+c write(*,*)'tau0,taus',tau0I,tausI
+ tausVOCE=tausI-tau0I
+ if(tausVOCE.LT.0d0) WRITE(*,*)'ERROR IN HARDENING',tausVOCE
+ expo=exp(-gamma*h0I/tausVOCE)
+ h=h1I*(1d0-expo)+(tausVOCE+h1I*gamma)*(h0I/tausVOCE)*expo
+
+ endif
+
+ RETURN
+ END
+
+ SUBROUTINE form_BIGRES(L_p,L_p_new,dgamma,dgamma_new,
+ 1 nsystems,BIGRES)
+ implicit none
+ real*8 L_p(3,3),L_p_new(3,3),dgamma(*),dgamma_new(*),BIGRES(*)
+ integer isystem,nsystems
+
+ BIGRES(1)=L_p(1,1)-L_p_new(1,1)
+ BIGRES(2)=L_p(2,2)-L_p_new(2,2)
+ BIGRES(3)=L_p(3,3)-L_p_new(3,3)
+ BIGRES(4)=L_p(1,2)-L_p_new(1,2)
+ BIGRES(5)=L_p(1,3)-L_p_new(1,3)
+ BIGRES(6)=L_p(2,3)-L_p_new(2,3)
+ BIGRES(7)=L_p(2,1)-L_p_new(2,1)
+ BIGRES(8)=L_p(3,1)-L_p_new(3,1)
+ BIGRES(9)=L_p(3,2)-L_p_new(3,2)
+
+ do,isystem=1,nsystems
+ BIGRES(9+isystem)=(dgamma(isystem)-dgamma_new(isystem))
+ enddo
+
+ return
+ end
+
diff --git a/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f b/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
new file mode 100644
index 0000000000000000000000000000000000000000..5a4d7a635994875a045ec071c7acf4409ef8059e
--- /dev/null
+++ b/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
@@ -0,0 +1,1666 @@
+C
+c CAPSUL: CrystAl Plasticity UtiLities tool set
+C Multiscale Materials Modelling Group, IMDEA Materials
+C
+C GENERAL SUBROUTINES FOR UMAT CRYSTAL PLASTICITY MODELS
+c
+c Javier Segurado, last modification 15/01/2016
+
+
+ SUBROUTINE STIFF6(c11,c12,c44,c13,c33,c66,stiff)
+
+C Generates the elastic stiffness matrix for an anisotropic crystal (i.e. HCP)
+C axis 3 is normal to isotropy plane (i.e. c-axis)
+C 11-->1, 22-->2, 33-->3, 13-->4, 23-->5, 33-->6
+ implicit none
+ REAL*8 c11,c12,c44,c13,c33,c66
+ REAL*8 stiff(3,3,3,3)
+ integer ii,jj,kk,ll
+
+ do,ii=1,3
+ do,jj=1,3
+ do,kk=1,3
+ do,ll=1,3
+ STIFF(ii,jj,kk,ll)=0D0
+ enddo
+ enddo
+ enddo
+ enddo
+
+ STIFF(1,1,1,1)=c11
+ STIFF(2,2,2,2)=c11
+ STIFF(3,3,3,3)=c33
+c
+ STIFF(1,1,2,2)=c12
+ STIFF(2,2,1,1)=c12
+ STIFF(1,1,3,3)=c13
+ STIFF(2,2,3,3)=c13
+ STIFF(3,3,1,1)=c13
+ STIFF(3,3,2,2)=c13
+c
+ STIFF(1,2,1,2)=c66
+ STIFF(2,1,2,1)=c66
+ STIFF(1,2,2,1)=c66
+ STIFF(2,1,1,2)=c66
+c
+ STIFF(1,3,1,3)=c44
+ STIFF(3,1,3,1)=c44
+ STIFF(1,3,3,1)=c44
+ STIFF(3,1,1,3)=c44
+c
+ STIFF(2,3,2,3)=c44
+ STIFF(3,2,3,2)=c44
+ STIFF(2,3,3,2)=c44
+ STIFF(3,2,2,3)=c44
+
+ RETURN
+ END
+
+ SUBROUTINE STIFF4(c11,c12,c44,stiff)
+ implicit none
+ REAL*8 c11,c12,c44,stiff(3,3,3,3)
+ integer ii,jj,kk,ll
+ do,ii=1,3
+ do,jj=1,3
+ do,kk=1,3
+ do,ll=1,3
+ STIFF(ii,jj,kk,ll)=0D0
+ enddo
+ enddo
+ enddo
+ enddo
+
+ STIFF(1,1,1,1)=c11
+ STIFF(2,2,2,2)=c11
+ STIFF(3,3,3,3)=c11
+c
+ STIFF(1,1,2,2)=c12
+ STIFF(1,1,3,3)=c12
+ STIFF(2,2,1,1)=c12
+ STIFF(2,2,3,3)=c12
+ STIFF(3,3,1,1)=c12
+ STIFF(3,3,2,2)=c12
+c
+ STIFF(1,2,1,2)=c44
+ STIFF(2,1,2,1)=c44
+ STIFF(1,2,2,1)=c44
+ STIFF(2,1,1,2)=c44
+c
+ STIFF(1,3,1,3)=c44
+ STIFF(3,1,3,1)=c44
+ STIFF(1,3,3,1)=c44
+ STIFF(3,1,1,3)=c44
+c
+ STIFF(2,3,2,3)=c44
+ STIFF(3,2,3,2)=c44
+ STIFF(2,3,3,2)=c44
+ STIFF(3,2,2,3)=c44
+ RETURN
+ END
+
+ SUBROUTINE green_lagrange(FE,epsilon_el)
+ implicit none
+ real*8 FE(3,3),epsilon_el(3,3)
+ real*8 FEt(3,3)
+ integer ii,jj
+
+ CALL TRANSPOSE(FEt,FE,3,3)
+ CALL PMAT(epsilon_el,FEt,FE,3,3,3)
+
+ do,ii=1,3
+ do,jj=1,3
+ epsilon_el(ii,jj)=epsilon_el(ii,jj)*.5d0
+ enddo
+ epsilon_el(ii,ii)=epsilon_el(ii,ii)-.5D0
+ enddo
+ return
+ end
+
+
+ SUBROUTINE lagrangian_def(CC,epsilon_el)
+ implicit none
+ real*8 cc(3,3),epsilon_el(3,3)
+ integer ii,jj
+c Lagrangian deformation
+
+ do,ii=1,3
+ do,jj=1,3
+ epsilon_el(ii,jj)=cc(ii,jj)*.5d0
+ enddo
+ epsilon_el(ii,ii)=epsilon_el(ii,ii)-.5D0
+ enddo
+ RETURN
+ END
+
+
+ FUNCTION secah(x)
+ implicit none
+ REAL*8 x,secah
+ if(abs(x).GT.50)THEN
+ secah=0d0
+ else
+ secah=2d0/(dexp(x)+dexp(-x))
+ endif
+ RETURN
+ END
+
+
+
+
+c *****************************************************************************
+ subroutine euler (iopt,ph,th,tm,a)
+c
+c CALCULATE THE EULER ANGLES ASSOCIATED WITH THE TRANSFORMATION
+c MATRIX A(I,J) IF IOPT=1 AND VICEVERSA IF IOPT=2
+c A(i,j) TRANSFORMS FROM SYSTEM sa TO SYSTEM ca.
+c ph,th,om ARE THE EULER ANGLES (in degrees) OF ca REFERRED TO sa.
+c Copied from VPSC code (Lebensohn et at 1993)
+c *****************************************************************************
+ real*8 ph,th,tm,a,pi,sph,cph,sth,cth,stm,ctm
+ integer iopt
+ dimension a(3,3)
+ pi=4.*datan(1.d0)
+
+ if(iopt.eq.1) then
+ if(abs(a(3,3)).ge.0.999999d0) then
+ if(a(3,3).GT.0) th=0d0
+ if(a(3,3).LT.0) th=180d0
+ tm=0d0
+ ph=datan2(a(1,2),a(1,1))
+ else
+ th=dacos(a(3,3))
+ sth=dsin(th)
+ tm=datan2(a(1,3)/sth,a(2,3)/sth)
+ ph=datan2(a(3,1)/sth,-a(3,2)/sth)
+ endif
+ th=th*180d0/pi
+ ph=ph*180d0/pi
+ tm=tm*180d0/pi
+ else if(iopt.eq.2) then
+ sph=dsin(ph*pi/180.)
+ cph=dcos(ph*pi/180.)
+ sth=dsin(th*pi/180.)
+ cth=dcos(th*pi/180.)
+ stm=dsin(tm*pi/180.)
+ ctm=dcos(tm*pi/180.)
+ a(1,1)=ctm*cph-sph*stm*cth
+ a(2,1)=-stm*cph-sph*ctm*cth
+ a(3,1)=sph*sth
+ a(1,2)=ctm*sph+cph*stm*cth
+ a(2,2)=-sph*stm+cph*ctm*cth
+ a(3,2)=-sth*cph
+ a(1,3)=sth*stm
+ a(2,3)=ctm*sth
+ a(3,3)=cth
+ endif
+
+ return
+ end
+
+
+
+c SUBROUTINE to orientate a 4 order rank tensor using a rotation matrix
+
+ SUBROUTINE ORIENTATE_TENSOR4(TENSOR_R,TENSOR,ROT)
+ implicit none
+ REAL*8 TENSOR_R(3,3,3,3),TENSOR(3,3,3,3),ROT(3,3)
+ integer ii,jj,kk,ll,mm,nn,pp,qq
+
+
+ do,ii=1,3
+ do,jj=1,3
+ do,kk=1,3
+ do,ll=1,3
+ TENSOR_R(ii,jj,kk,ll)=0d0
+
+ do,mm=1,3
+ do,nn=1,3
+ do,pp=1,3
+ do,qq=1,3
+ TENSOR_R(ii,jj,kk,ll)=
+ 1 TENSOR_R(ii,jj,kk,ll)+
+ 2 ROT(mm,ii)*ROT(nn,jj)*
+ 3 ROT(pp,kk)*ROT(qq,ll)*
+ 4 TENSOR(mm,nn,pp,qq)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ RETURN
+ END
+c
+
+c
+C Polar decomposition by Simo 1991
+
+ SUBROUTINE POLAR_DECOMP(UU,RR,CC,DEFGRAD)
+ implicit none
+ REAL*8 UU(3,3),UUINV(3,3),RR(3,3),DEFGRAD(3,3),DEFGRAD_T(3,3)
+ REAL*8 CC(3,3),CC2(3,3),unit(3,3)
+ REAL*8 I1,I2,I3,iu1,iu2,iu3
+ REAL*8 x(3),lambda(3),lambda2(3)
+ REAL*8 b,c,m,n,t,D,dd
+ real*8 a,p,q,phi,temp1,y1,y2,y3,temp2,u,v,y2r,y2i
+ INTEGER l,ii,jj,ndif
+
+ real*8 pi
+
+
+ PI=3.141592653589793238D0
+
+ do,ii=1,3
+ do,jj=1,3
+ unit(ii,jj)=0d0
+ enddo
+ unit(ii,ii)=1d0
+ enddo
+
+c Transpose F
+ CALL TRANSPOSE(DEFGRAD_T,DEFGRAD,3,3)
+
+c CC=Ft*F
+ CALL PMAT(CC,DEFGRAD_T,DEFGRAD,3,3,3)
+c write(*,*) 'cc',CC
+
+C Diagonalization of C (obtention of lambda(i))
+ CALL INVARIANTS(I1,I2,I3,CC)
+
+ b = I2 - I1*I1/3d0
+ c = (-2d0/27d0)*I1*I1*I1 + I1*I2/3d0 - I3
+
+ if(abs(b).LE.1D-15) then
+
+ do,l=1,3
+ x(l)=-c**(1D0/3d0)
+ enddo
+
+c write(*,*)'b<1e-12', x(1),x(2),x(3)
+
+ else
+ m = 2d0*dsqrt(-b/3D0)
+ n = 3d0*c/(m*b)
+ t = datan2(dsqrt(abs(1d0-n*n)),n)/3D0
+ do,l=1,3
+ x(l) = m*dcos(t+2d0*(l-1)*PI/3d0)
+ enddo
+ endif
+
+ do,l=1,3
+ lambda(l)=dsqrt(x(l)+I1/3D0)
+ enddo
+
+
+c Invariants of U
+
+ iu1 = lambda(1) + lambda(2) + lambda(3)
+ iu2 = lambda(1)*lambda(2) + lambda(1)*lambda(3)
+ 1 +lambda(2)*lambda(3)
+ iu3 = lambda(1) * lambda(2) * lambda(3)
+
+ D= iu1*iu2 - iu3
+ IF(ABS(D).le.1d-15) write(*,*)'ERROR D',D
+
+ CALL PMAT(CC2,CC,CC,3,3,3)
+
+ do,ii=1,3
+ do,jj=1,3
+ UU(ii,jj)=(1d0/D)*(-CC2(ii,jj)+(iu1*iu1-iu2)*CC(ii,jj)+
+ 1 iu1*iu3*unit(ii,jj))
+ UUINV(ii,jj)=(1d0/iu3)*(CC(ii,jj)-iu1*UU(ii,jj)+
+ 1 iu2*unit(ii,jj))
+ enddo
+ enddo
+
+C ROTATION
+
+ CALL PMAT(RR,DEFGRAD,UUINV,3,3,3)
+
+c write(*,*)'uu',uu
+c write(*,*)'rr',rr
+
+ RETURN
+ END
+
+
+ SUBROUTINE DECOMP(Matrix,lambda,vec1,vec2,vec3)
+ implicit none
+ REAL*8 Matrix(3,3),lambda(3),vec1(3),vec2(3),vec3(3)
+ REAL*8 I1,I2,I3
+ real*8 a,d
+ real*8 p,q,dd,phi,temp1,y1,y2,y3,temp2,u,v,y2r,y2i
+ REAL*8 Matrix2(3,3)
+ REAL*8 dnorm,det(3),detmax
+ real*8 b,c,x(3),n,t,m,pi
+ INTEGER ii,jj,ndif,nroot,i,ni,l
+
+
+
+ PI=3.141592653589793238D0
+C Diagonalization of matrix (sym) (obtention of lambda(i))
+ CALL INVARIANTS(I1,I2,I3,matrix)
+
+
+c write(*,*)'invariants',i1,i2,i3
+
+ b = I2 - I1*I1/3.d0
+ c = (-2/27.d0)*I1*I1*I1 + I1*I2/3.d0 - I3
+
+c write(*,*)'b,c',b,c
+
+ if(abs(b).LE.1D-15) then
+
+ do,l=1,3
+ x(l)=-c**(1.D0/3.D0)
+ enddo
+
+ ndif=1
+c write(*,*) x(1),x(2),x(3)
+
+ else
+
+ m = 2*dsqrt(-b/3.D0)
+ n = 3*c/(m*b)
+ if(n.GT.1d0) THEN
+c write(*,*)'RU error, SQRT() neg',1d0-n*n
+ n=0d0
+ endif
+ t = datan2(dsqrt(1d0-n*n),n)/3.D0
+
+ if(abs(t).LE.1D-15) then
+ ndif=2
+ else
+ ndif=3
+ endif
+c write(*,*) m,n,t
+
+ do,l=1,3
+ x(l) = m*dcos(t+2*(l-1)*PI/3.d0)
+ enddo
+
+c write(*,*) x(1),x(2),x(3)
+
+ endif
+
+ do,l=1,3
+ lambda(l)=dsqrt(x(l)+I1/3D0)
+ enddo
+
+ write(*,*) 'lambdas',lambda(1),lambda(2),lambda(3)
+
+
+c$$$
+c$$$
+c$$$
+c$$$
+c$$$
+c$$$ CALL INVARIANTS(I1,I2,I3,Matrix)
+c$$$
+c$$$c write(*,*)'invariants',i1,i2,i3
+c$$$
+c$$$ a=1d0
+c$$$ b=-i1
+c$$$ c=i2
+c$$$ d=-i3
+c$$$c PI=3.141592653589793238D0
+c$$$ pi=4d0*datan(1.d0)
+c$$$
+c$$$c Step 1: Calculate p and q --------------------------------------------
+c$$$ p = c/a - b*b/a/a/3D0
+c$$$ q = (2D0*b*b*b/a/a/a - 9D0*b*c/a/a + 27D0*d/a) / 27D0
+c$$$
+c$$$c Step 2: Calculate DD (discriminant) ----------------------------------
+c$$$ DD = p*p*p/27D0 + q*q/4D0
+c$$$c write(*,*)'DD',dd
+c$$$
+c$$$c Step 3: Branch to different algorithms based on DD -------------------
+c$$$ if(DD.lt.0.and.abs(dd/i3).GT.1D-15) then
+c$$$c Step 3b:
+c$$$c 3 real unequal roots -- use the trigonometric formulation
+c$$$ phi = dacos(-q/2D0/dsqrt(dabs(p*p*p)/27D0))
+c$$$ temp1=2D0*dsqrt(dabs(p)/3D0)
+c$$$ y1 = temp1*dcos(phi/3D0)
+c$$$ y2 = -temp1*dcos((phi+pi)/3D0)
+c$$$ y3 = -temp1*dcos((phi-pi)/3D0)
+c$$$ else
+c$$$ if(abs(DD/i3).LE.1D-15) DD=0d0
+c$$$c Step 3a:
+c$$$c 1 real root & 2 conjugate complex roots OR 3 real roots (some are equal)
+c$$$ temp1 = -q/2D0 + dsqrt(DD)
+c$$$ temp2 = -q/2D0 - dsqrt(DD)
+c$$$ u = dabs(temp1)**(1D0/3D0)
+c$$$ v = dabs(temp2)**(1D0/3D0)
+c$$$ if(temp1 .lt. 0D0) u=-u
+c$$$ if(temp2 .lt. 0D0) v=-v
+c$$$ y1 = u + v
+c$$$ y2r = -(u+v)/2D0
+c$$$ y2i = (u-v)*dsqrt(3D0)/2D0
+c$$$ endif
+c$$$
+c$$$c Step 4: Final transformation -----------------------------------------
+c$$$ temp1 = b/a/3D0
+c$$$ y1 = y1-temp1
+c$$$ y2 = y2-temp1
+c$$$ y3 = y3-temp1
+c$$$ y2r=y2r-temp1
+c$$$
+c$$$c Assign answers and number of lambdas differents-------------------------------------------------------
+c$$$ if(DD.lt.0.and.abs(dd/i3).GT.1D-15) then
+c$$$
+c$$$ lambda(1) = y1
+c$$$ lambda(2) = y2
+c$$$ lambda(3) = y3
+c$$$ ndif=3
+c$$$
+c$$$ elseif(abs(DD/i3).LE.1D-15) then
+c$$$
+c$$$ if(abs(y1-y2r).LE.1D-15) then
+c$$$ ndif=1
+c$$$ lambda(1)=y1
+c$$$ lambda(2)=y2r
+c$$$ lambda(3)=y2r
+c$$$ else
+c$$$ ndif=2
+c$$$ if(y1.GT.y2r) THEN
+c$$$ lambda(1)=y1
+c$$$ lambda(2)=y2r
+c$$$ lambda(3)=y2r
+c$$$ else
+c$$$ lambda(1)=y2r
+c$$$ lambda(2)=y2r
+c$$$ lambda(3)=y1
+c$$$ endif
+c$$$
+c$$$ endif
+c$$$
+c$$$ else
+c$$$
+c$$$ write(*,*)'error:COMPLEx'
+c$$$
+c$$$ endif
+c$$$
+c$$$c write(*,*)'NDIF',ndif
+c$$$c write(*,*)'LAMBDAS',lambda(1),lambda(2),lambda(3)
+c$$$
+C EIGENVECTORS
+
+C three different
+ IF(ndif.EQ.3) THEN
+
+c FIRST
+ do,ii=1,3
+ do,jj=1,3
+ Matrix2(ii,jj)=Matrix(ii,jj)
+ if(ii.EQ.jj) MATRIX2(ii,jj)=MATRIX2(ii,jj)-lambda(1)
+ enddo
+ enddo
+
+ det(1)=Matrix2(1,1)*Matrix2(2,2)-Matrix2(1,2)*Matrix2(2,1)
+ det(2)=Matrix2(2,2)*Matrix2(3,3)-Matrix2(2,3)*Matrix2(3,2)
+ det(3)=Matrix2(1,1)*Matrix2(3,3)-Matrix2(1,3)*Matrix2(3,1)
+ detmax=0d0
+ do,i=1,3
+ if(abs(det(i)).GT.detmax) THEN
+ ni=i
+ detmax=abs(det(i))
+ endif
+ enddo
+c write(*,*)'dets',det
+ if(ni.EQ.1) THEN
+c write(*,*)'supuesto1',detmax
+ a=(1/det(ni))*
+ 1 (-Matrix2(2,2)*Matrix2(1,3)+Matrix2(2,1)*Matrix2(2,3))
+ b=(1/det(ni))*
+ 1 (Matrix2(1,2)*Matrix2(1,3)-Matrix2(1,1)*Matrix2(2,3))
+ c=1d0
+ elseif(ni.EQ.2) THEN
+c write(*,*)'supuesto2',detmax
+ a=1d0
+ b=(1/det(ni))*
+ 1 (-Matrix2(3,3)*Matrix2(2,1)+Matrix2(3,2)*Matrix2(3,1))
+ c=(1/det(ni))*
+ 1 (Matrix2(2,3)*Matrix2(2,1)-Matrix2(2,2)*Matrix2(3,1))
+ else
+c write(*,*)'supuesto3',detmax
+ b=1d0
+ a=(1/det(ni))*
+ 1 (-Matrix2(3,3)*Matrix2(1,2)+Matrix2(3,1)*Matrix2(2,3))
+ c=(1/det(ni))*
+ 1 (Matrix2(1,3)*Matrix2(1,2)-Matrix2(1,1)*Matrix2(2,3))
+ endif
+
+ dnorm=dsqrt(a*a+b*b+c*c)
+ vec1(1)=a/dnorm
+ vec1(2)=b/dnorm
+ vec1(3)=c/dnorm
+c SECOND
+
+ do,ii=1,3
+ do,jj=1,3
+ Matrix2(ii,jj)=Matrix(ii,jj)
+ if(ii.EQ.jj) MATRIX2(ii,jj)=MATRIX2(ii,jj)-lambda(2)
+ enddo
+ enddo
+ detmax=0d0
+ det(1)=Matrix2(1,1)*Matrix2(2,2)-Matrix2(1,2)*Matrix2(2,1)
+ det(2)=Matrix2(2,2)*Matrix2(3,3)-Matrix2(2,3)*Matrix2(3,2)
+ det(3)=Matrix2(1,1)*Matrix2(3,3)-Matrix2(1,3)*Matrix2(3,1)
+c write(*,*)'dets',det
+ do,i=1,3
+ if(abs(det(i)).GT.detmax) THEN
+ ni=i
+ detmax=abs(det(i))
+ endif
+ enddo
+ if(ni.EQ.1) THEN
+c write(*,*)'supuesto1',detmax
+ a=(1/det(ni))*
+ 1 (-Matrix2(2,2)*Matrix2(1,3)+Matrix2(2,1)*Matrix2(2,3))
+ b=(1/det(ni))*
+ 1 (Matrix2(1,2)*Matrix2(1,3)-Matrix2(1,1)*Matrix2(2,3))
+ c=1d0
+ elseif(ni.EQ.2) THEN
+c write(*,*)'supuesto2',detmax
+ a=1d0
+ b=(1/det(ni))*
+ 1 (-Matrix2(3,3)*Matrix2(2,1)+Matrix2(3,2)*Matrix2(3,1))
+ c=(1/det(ni))*
+ 1 (Matrix2(2,3)*Matrix2(2,1)-Matrix2(2,2)*Matrix2(3,1))
+ else
+c write(*,*)'supuesto3',detmax
+ b=1d0
+ a=(1/det(ni))*
+ 1 (-Matrix2(3,3)*Matrix2(1,2)+Matrix2(3,1)*Matrix2(2,3))
+ c=(1/det(ni))*
+ 1 (Matrix2(1,3)*Matrix2(1,2)-Matrix2(1,1)*Matrix2(2,3))
+ endif
+
+ dnorm=dsqrt(a*a+b*b+c*c)
+ vec2(1)=a/dnorm
+ vec2(2)=b/dnorm
+ vec2(3)=c/dnorm
+
+c THIRD
+ CALL PVECT(vec3,vec1,vec2)
+c write(*,*)'v1:',vec1
+c write(*,*)'v2:',vec2
+c write(*,*)'v3:',vec3
+
+C TWO DIFFERENT EIGENVALUES
+ ELSE IF(ndif.EQ.2) THEN
+
+ if(y2r.GT.y1) then
+ do,ii=1,3
+ do,jj=1,3
+ Matrix2(ii,jj)=Matrix(ii,jj)
+ if(ii.EQ.jj) MATRIX2(ii,jj)=MATRIX2(ii,jj)-lambda(3)
+ enddo
+ enddo
+
+ else
+
+ do,ii=1,3
+ do,jj=1,3
+ Matrix2(ii,jj)=Matrix(ii,jj)
+ if(ii.EQ.jj) MATRIX2(ii,jj)=MATRIX2(ii,jj)-lambda(1)
+ enddo
+ enddo
+
+ endif
+
+ detmax=0d0
+ det(1)=Matrix2(1,1)*Matrix2(2,2)-Matrix2(1,2)*Matrix2(2,1)
+ det(2)=Matrix2(2,2)*Matrix2(3,3)-Matrix2(2,3)*Matrix2(3,2)
+ det(3)=Matrix2(1,1)*Matrix2(3,3)-Matrix2(1,3)*Matrix2(3,1)
+ do,i=1,3
+ if(abs(det(i)).GT.detmax) THEN
+ ni=i
+ detmax=abs(det(i))
+ endif
+ enddo
+ if(ni.EQ.1) THEN
+c write(*,*)'supuesto1',detmax
+ a=(1/det(ni))*
+ 1 (-Matrix2(2,2)*Matrix2(1,3)+Matrix2(2,1)*Matrix2(2,3))
+ b=(1/det(ni))*
+ 1 (Matrix2(1,2)*Matrix2(1,3)-Matrix2(1,1)*Matrix2(2,3))
+ c=1d0
+ elseif(ni.EQ.2) THEN
+c write(*,*)'supuesto2',detmax
+ a=1d0
+ b=(1/det(ni))*
+ 1 (-Matrix2(3,3)*Matrix2(2,1)+Matrix2(3,2)*Matrix2(3,1))
+ c=(1/det(ni))*
+ 1 (Matrix2(2,3)*Matrix2(2,1)-Matrix2(2,2)*Matrix2(3,1))
+ else
+c write(*,*)'supuesto3',detmax
+ b=1d0
+ a=(1/det(ni))*
+ 1 (-Matrix2(3,3)*Matrix2(1,2)+Matrix2(3,1)*Matrix2(2,3))
+ c=(1/det(ni))*
+ 1 (Matrix2(1,3)*Matrix2(1,2)-Matrix2(1,1)*Matrix2(2,3))
+ endif
+
+ dnorm=dsqrt(a*a+b*b+c*c)
+
+ if(y2r.GT.y1) then
+ vec3(1)=a/dnorm
+ vec3(2)=b/dnorm
+ vec3(3)=c/dnorm
+ if(abs(vec3(3)).GT.1D-15) THEN
+
+ vec1(1)=-vec3(3)/dsqrt(vec3(1)**2+vec3(3)**2)
+ vec1(2)=0d0
+ vec1(3)=vec3(1)/dsqrt(vec3(1)**2+vec3(3)**2)
+
+ else
+ vec1(1)=0d0
+ vec1(2)=0d0
+ vec1(3)=1d0
+ endif
+ call pvect(vec2,vec3,vec1)
+
+
+ else
+
+ vec1(1)=a/dnorm
+ vec1(2)=b/dnorm
+ vec1(3)=c/dnorm
+ if(abs(vec1(3)).GT.1D-15) THEN
+ vec2(1)=-vec1(3)/dsqrt(vec1(1)**2+vec1(3)**2)
+ vec2(2)=0d0
+ vec2(3)=vec1(1)/dsqrt(vec1(1)**2+vec1(3)**2)
+ else
+ vec2(1)=0d0
+ vec2(2)=0d0
+ vec2(3)=1d0
+ endif
+ call pvect(vec3,vec1,vec2)
+
+
+
+ endif
+
+c lambda1=lambda2=lambda3
+ ELSE
+ vec1(1)=1d0
+ vec1(2)=0d0
+ vec1(3)=0d0
+ vec2(1)=0d0
+ vec2(2)=1d0
+ vec2(3)=0d0
+ vec3(1)=0d0
+ vec3(2)=0d0
+ vec3(3)=1d0
+
+ ENDIF
+
+ RETURN
+ END
+
+
+ SUBROUTINE LOGtens(matrix,lambda,vec1,vec2,vec3)
+ REAL*8 matrix(3,3),lambda(3),vec1(3),vec2(3),vec3(3)
+ REAL*8 TENS1(3,3),TENS2(3,3),TENS3(3,3)
+ do,ii=1,3
+ lambda(ii)=log(lambda(ii))
+ enddo
+ CALL PTENS(TENS1,vec1,vec1)
+ CALL PTENS(TENS2,vec2,vec2)
+ CALL PTENS(TENS3,vec3,vec3)
+ do,ii=1,3
+ do,jj=1,3
+ matrix(ii,jj)=lambda(1)*TENS1(ii,jj)+
+ 1 lambda(2)*TENS2(ii,jj)+
+ 1 lambda(3)*TENS3(ii,jj)
+ enddo
+ enddo
+ RETURN
+ END
+
+ SUBROUTINE INVARIANTS(I1,I2,I3,CC)
+c Calculate invariants of a symmetric matrix
+ implicit none
+ REAL*8 I1,I2,I3
+ REAL*8 CC(3,3),CC2(3,3)
+ real*8 trac2
+
+
+ I1=CC(1,1)+CC(2,2)+CC(3,3)
+
+ CALL PMAT(CC2,CC,CC,3,3,3)
+
+ trac2=CC2(1,1)+CC2(2,2)+CC2(3,3)
+
+ I2=0.5d0*(I1*I1-trac2)
+
+ I3=CC(1,1)*CC(2,2)*CC(3,3)+CC(1,2)*CC(2,3)*CC(3,1)
+ 1 +CC(1,3)*CC(2,1)*CC(3,2)
+ 2 -CC(1,3)*CC(2,2)*CC(3,1)-CC(2,3)*CC(3,2)*CC(1,1)
+ 3 -CC(3,3)*CC(1,2)*CC(2,1)
+
+ RETURN
+ END
+
+ SUBROUTINE PMAT(matout,mat1,mat2,id,jd,kd)
+C Matrix product
+ REAL*8 MAT1,MAT2,MATOUT
+ DIMENSION mat1(id,kd),mat2(kd,jd),matout(id,jd)
+ do,i=1,id
+ Do,j=1,jd
+ matout(i,j)=0D0
+ do,k=1,kd
+ matout(i,j)=matout(i,j)+mat1(i,k)*mat2(k,j)
+ enddo
+ enddo
+ enddo
+ RETURN
+ END
+c
+ SUBROUTINE PTENS(matout,vectors1,vectors2)
+c Tensorial product of 2 vectors
+ integer i,j
+ real*8 vectors1(3),vectors2(3),matout(3,3)
+ do,i=1,3
+ do,j=1,3
+ matout(i,j)=vectors1(i)*vectors2(j)
+ enddo
+ enddo
+ RETURN
+ END
+
+ SUBROUTINE PTENS2(matout4,tens1,tens2)
+ integer i,j,k,l
+ real*8 tens1(3,3),tens2(3,3),matout4(3,3,3,3)
+ do,i=1,3
+ do,j=1,3
+ do,k=1,3
+ do,l=1,3
+ matout4(i,j,k,l)=tens1(i,j)*tens2(k,l)
+ enddo
+ enddo
+ enddo
+ enddo
+ return
+ end
+
+ SUBROUTINE TRANSPOSE(Mat2,Mat1,ifil,icol)
+ REAL*8 MAT1,MAT2
+ DIMENSION MAT1(ifil,icol),MAT2(icol,ifil)
+ DO,i=1,ifil
+ do,j=1,icol
+ MAT2(j,i)=0D0
+ MAT2(j,i)=MAT1(i,j)
+C write(8,*) j,i,'-->',i,j,MAT2(i,j)
+ enddo
+ enddo
+ RETURN
+ END
+
+ SUBROUTINE PCONTRACT(pcont,MAT1,MAT2,n)
+ REAL*8 MAT1(n,n),MAT2(n,n),pcont
+ INTEGER n,ii,jj
+ pcont=0d0
+ do,ii=1,n
+ do,jj=1,n
+ pcont=pcont+MAT1(ii,jj)*MAT2(ii,jj)
+ enddo
+ enddo
+ RETURN
+ END
+
+ SUBROUTINE PCONTRACT2(pcont,MAT1,MAT2,n)
+ REAL*8 MAT1(n,n,n,n),MAT2(n,n),pcont(n,n)
+ INTEGER n,ii,jj,kk,ll
+ do,ii=1,n
+ do,jj=1,n
+ pcont(ii,jj)=0d0
+ do,kk=1,n
+ do,ll=1,n
+ pcont(ii,jj)=pcont(ii,jj)+MAT1(ii,jj,kk,ll)*
+ 1 MAT2(kk,ll)
+ enddo
+ enddo
+ enddo
+ enddo
+
+
+ RETURN
+ END
+
+
+ SUBROUTINE PVECT(vecout,vectors1,vectors2)
+C Vectorial product
+ REAL*8 vectors1,vectors2,vecout
+ DIMENSION vectors1(3),vectors2(3),vecout(3)
+ Dimension ii(3,2)
+ ii(1,1)=2
+ ii(1,2)=3
+ ii(2,1)=1
+ ii(2,2)=3
+ ii(3,1)=1
+ ii(3,2)=2
+ Do,i=1,3
+ vecout(i)=0d0
+ vecout(i)=(-1)**(i+1)*( vectors1(ii(i,1))*vectors2(ii(i,2))-
+ 1 vectors2(ii(i,1))*vectors1(ii(i,2)) )
+ ENDDO
+
+ RETURN
+ END
+
+c
+ SUBROUTINE MATINV3(CC2,CC,iflag)
+ REAL*8 CC(3,3),CC2(3,3)
+ REAL*8 I3
+ integer iflag
+
+ I3=CC(1,1)*CC(2,2)*CC(3,3)+CC(1,2)*CC(2,3)*CC(3,1)
+ 1 +CC(1,3)*CC(2,1)*CC(3,2)
+ 2 -CC(1,3)*CC(2,2)*CC(3,1)-CC(2,3)*CC(3,2)*CC(1,1)
+ 3 -CC(3,3)*CC(1,2)*CC(2,1)
+
+ if(abs(I3).GT.1D-50) THEN
+ iflag=0
+ I3=1D0/I3
+ else
+ iflag=1
+ return
+ endif
+
+ CC2(1,1)=I3*(cc(2,2)*cc(3,3)-cc(2,3)*cc(3,2))
+ CC2(2,2)=I3*(cc(1,1)*cc(3,3)-cc(1,3)*cc(3,1))
+ CC2(3,3)=I3*(cc(1,1)*cc(2,2)-cc(1,2)*cc(2,1))
+
+ CC2(1,2)=-I3*(cc(1,2)*cc(3,3)-cc(3,2)*cc(1,3))
+ CC2(2,1)=-I3*(cc(2,1)*cc(3,3)-cc(2,3)*cc(3,1))
+
+ CC2(1,3)=I3*(cc(1,2)*cc(2,3)-cc(1,3)*cc(2,2))
+ CC2(3,1)=I3*(cc(2,1)*cc(3,2)-cc(2,2)*cc(3,1))
+
+ CC2(2,3)=-I3*(cc(1,1)*cc(2,3)-cc(1,3)*cc(2,1))
+ CC2(3,2)=-I3*(cc(1,1)*cc(3,2)-cc(1,2)*cc(3,1))
+
+ RETURN
+ END
+
+
+ SUBROUTINE EXPONENTIAL3(EXPO,MATRIX)
+ REAL*8 EXPO(3,3),MATRIX(3,3)
+ REAL*8 AUX33(3,3),MATRIXn(3,3)
+ REAL*8 dnorm,toler
+ toler=1D-5
+ do,ii=1,3
+ do,jj=1,3
+ EXPO(ii,jj)=0d0
+ MATRIXn(ii,jj)=0d0
+ enddo
+ EXPO(ii,ii)=1d0
+ MATRIXn(ii,ii)=1d0
+ enddo
+ n=0
+ nfact=1
+ dnorm=1D0
+ DO 100 WHILE(dnorm.GE.toler)
+ n=n+1
+ nfact=nfact*n
+ CALL PMAT(AUX33,MATRIX,MATRIXn,3,3,3)
+ do,ii=1,3
+ do,jj=1,3
+ MATRIXn(ii,jj)=AUX33(ii,jj)
+ EXPO(ii,jj)=EXPO(ii,jj)+(1D0/nfact)*MATRIXn(ii,jj)
+ enddo
+ enddo
+
+ dnorm=0d0
+ do,ii=1,3
+ do,jj=1,3
+ dnorm=dnorm+MATRIXn(ii,jj)*MATRIXn(ii,jj)
+ enddo
+ enddo
+ dnorm=dsqrt(dnorm)/nfact
+c write(*,*) dnorm
+ 100 CONTINUE
+ RETURN
+ 102 FORMAT(3(3(E16.4,' '),/))
+
+ END
+
+ SUBROUTINE EXPMAP
+ 1 ( EXPX ,NOCONV ,X )
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ PARAMETER
+ 1 ( NDIM=3 ,NDIM2=9 )
+C Arguments
+ LOGICAL NOCONV
+ DIMENSION
+ 1 EXPX(NDIM,NDIM) ,X(NDIM,NDIM)
+C Local arrays and variables
+ DIMENSION
+ 1 XN(NDIM,NDIM) ,XNM1(NDIM,NDIM) ,XNM2(NDIM,NDIM) ,
+ 2 XNM3(NDIM,NDIM) ,X2(NDIM,NDIM)
+ DATA
+ 1 R0 ,RP5 ,R1 ,R2 ,TOL ,OVER ,UNDER /
+ 2 0.0D0,0.5D0,1.0D0,2.0D0,1.0D-5,1.0D+100,1.0D-100/
+ DATA
+ 1 NMAX / 100 /
+C***********************************************************************
+C COMPUTES THE EXPONENTIAL OF A (GENERALLY UNSYMMETRIC) 3-D TENSOR. USES
+C THE SERIES REPRESENTATION OF THE TENSOR EXPONENTIAL
+C
+C REFERENCE: Section B.1
+C Box B.1
+C***********************************************************************
+C Initialise series convergence flag
+ NOCONV=.FALSE.
+C Compute X square
+ CALL RVZERO(X2,NDIM2)
+ DO 30 I=1,NDIM
+ DO 20 J=1,NDIM
+ DO 10 K=1,NDIM
+ X2(I,J)= X2(I,J)+X(I,K)*X(K,J)
+ 10 CONTINUE
+ 20 CONTINUE
+ 30 CONTINUE
+C Compute principal invariants of X
+ C1=X(1,1)+X(2,2)+X(3,3)
+ C2=RP5*(C1*C1-(X2(1,1)+X2(2,2)+X2(3,3)))
+ C3=X(1,1)*X(2,2)*X(3,3)+X(1,2)*X(2,3)*X(3,1)+
+ 1 X(1,3)*X(2,1)*X(3,2)-X(1,2)*X(2,1)*X(3,3)-
+ 2 X(1,1)*X(2,3)*X(3,2)-X(1,3)*X(2,2)*X(3,1)
+C Start computation of exponential using its series definition
+C ============================================================
+ DO 50 I=1,NDIM
+ DO 40 J=1,NDIM
+ XNM1(I,J)=X2(I,J)
+ XNM2(I,J)=X(I,J)
+ 40 CONTINUE
+ 50 CONTINUE
+ XNM3(1,1)=R1
+ XNM3(1,2)=R0
+ XNM3(1,3)=R0
+ XNM3(2,1)=R0
+ XNM3(2,2)=R1
+ XNM3(2,3)=R0
+ XNM3(3,1)=R0
+ XNM3(3,2)=R0
+ XNM3(3,3)=R1
+C Add first three terms of series
+C -------------------------------
+ DO 70 I=1,NDIM
+ DO 60 J=1,NDIM
+ EXPX(I,J)=RP5*XNM1(I,J)+XNM2(I,J)+XNM3(I,J)
+ 60 CONTINUE
+ 70 CONTINUE
+C Add remaining terms (with X to the powers 3 to NMAX)
+C ----------------------------------------------------
+ FACTOR=R2
+ DO 140 N=3,NMAX
+C Use recursive formula to obtain X to the power N
+ DO 90 I=1,NDIM
+ DO 80 J=1,NDIM
+ XN(I,J)=C1*XNM1(I,J)-C2*XNM2(I,J)+C3*XNM3(I,J)
+ 80 CONTINUE
+ 90 CONTINUE
+C Update factorial
+ FACTOR=DBLE(N)*FACTOR
+ R1DFAC=R1/FACTOR
+C Add Nth term of the series
+ DO 110 I=1,NDIM
+ DO 100 J=1,NDIM
+ EXPX(I,J)=EXPX(I,J)+R1DFAC*XN(I,J)
+ 100 CONTINUE
+ 110 CONTINUE
+C Check convergence of series
+ XNNORM=SQRT(SCAPRD(XN(1,1),XN(1,1),NDIM2))
+ IF(XNNORM.GT.OVER.OR.(XNNORM.LT.UNDER.AND.XNNORM.GT.R0)
+ 1 .OR.R1DFAC.LT.UNDER)THEN
+C... first check possibility of overflow or underflow.
+C... numbers are to small or too big: Break (unconverged) loop and exit
+ NOCONV=.TRUE.
+ GOTO 999
+ ELSEIF(XNNORM*R1DFAC.LT.TOL)THEN
+C... converged: Break series summation loop and exit with success
+ GOTO 999
+ ENDIF
+ DO 130 I=1,NDIM
+ DO 120 J=1,NDIM
+ XNM3(I,J)=XNM2(I,J)
+ XNM2(I,J)=XNM1(I,J)
+ XNM1(I,J)=XN(I,J)
+ 120 CONTINUE
+ 130 CONTINUE
+ 140 CONTINUE
+C Re-set convergence flag if series did not converge
+ NOCONV=.TRUE.
+C
+ 999 CONTINUE
+ RETURN
+ END
+
+
+
+ SUBROUTINE DEXPMP
+ 1( DEXPX ,NOCONV ,X )
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ PARAMETER
+ 1( NDIM=3 ,NDIM2=9 ,NDIM4=81 ,MAXN=100 )
+C Arguments
+ LOGICAL NOCONV
+ DIMENSION
+ 1 DEXPX(NDIM,NDIM,NDIM,NDIM),X(NDIM,NDIM)
+C Local arrays and variables
+C...matrix of powers of X
+ DIMENSION
+ 1 R1DFAC(MAXN) ,XMATX(NDIM,NDIM,0:MAXN)
+C...initialise identity matrix: X to the power 0
+ DATA
+ 1 XMATX(1,1,0) ,XMATX(1,2,0) ,XMATX(1,3,0) /
+ 2 1.D0 ,0.D0 ,0.D0 /
+ 3 XMATX(2,1,0) ,XMATX(2,2,0) ,XMATX(2,3,0) /
+ 4 0.D0 ,1.D0 ,0.D0 /
+ 5 XMATX(3,1,0) ,XMATX(3,2,0) ,XMATX(3,3,0) /
+ 6 0.D0 ,0.D0 ,1.D0 /
+ DATA
+ 1 R0 ,RP5 ,R1 ,TOL ,OVER ,UNDER /
+ 2 0.0D0,0.5D0,1.0D0,1.0D-10,1.0D+100,1.0D-100/
+C***********************************************************************
+C COMPUTES THE DERIVATIVE OF THE EXPONENTIAL OF A (GENERALLY
+C UNSYMMETRIC) 3-D TENSOR. USES THE SERIES REPRESENTATION OF THE TENSOR
+C EXPONENTIAL.
+C
+C REFERENCE: Section B.2
+C Box B.2
+C***********************************************************************
+C Initialise convergence flag
+ NOCONV=.FALSE.
+C X to the power 1
+ DO 20 I=1,NDIM
+ DO 10 J=1,NDIM
+ XMATX(I,J,1)=X(I,J)
+ 10 CONTINUE
+ 20 CONTINUE
+C Zero remaining powers of X
+ CALL RVZERO(XMATX(1,1,2),NDIM*NDIM*(MAXN-1))
+C Compute X square
+ DO 50 I=1,NDIM
+ DO 40 J=1,NDIM
+ DO 30 K=1,NDIM
+ XMATX(I,J,2)=XMATX(I,J,2)+X(I,K)*X(K,J)
+ 30 CONTINUE
+ 40 CONTINUE
+ 50 CONTINUE
+C Compute principal invariants of X
+ C1=X(1,1)+X(2,2)+X(3,3)
+ C2=RP5*(C1*C1-(XMATX(1,1,2)+XMATX(2,2,2)+XMATX(3,3,2)))
+ C3=X(1,1)*X(2,2)*X(3,3)+X(1,2)*X(2,3)*X(3,1)+
+ 1 X(1,3)*X(2,1)*X(3,2)-X(1,2)*X(2,1)*X(3,3)-
+ 2 X(1,1)*X(2,3)*X(3,2)-X(1,3)*X(2,2)*X(3,1)
+C Compute X to the powers 3,4,...,NMAX using recursive formula
+ R1DFAC(1)=R1
+ R1DFAC(2)=RP5
+ DO 80 N=3,MAXN
+ R1DFAC(N)=R1DFAC(N-1)/DBLE(N)
+ DO 70 I=1,NDIM
+ DO 60 J=1,NDIM
+ XMATX(I,J,N)=C1*XMATX(I,J,N-1)-C2*XMATX(I,J,N-2)+
+ 1 C3*XMATX(I,J,N-3)
+ 60 CONTINUE
+ 70 CONTINUE
+ XNNORM=SQRT(SCAPRD(XMATX(1,1,N),XMATX(1,1,N),NDIM2))
+C...check number of terms required for series convergence
+ IF(XNNORM.GT.OVER.OR.(XNNORM.LT.UNDER.AND.XNNORM.GT.R0)
+ 1 .OR.R1DFAC(N).LT.UNDER)THEN
+C...numbers are to small or too big: Exit without computing derivative
+ NOCONV=.TRUE.
+ GOTO 999
+ ELSEIF(XNNORM*R1DFAC(N).LT.TOL)THEN
+C...series will converge with NMAX terms:
+C Carry on to derivative computation
+ NMAX=N
+ GOTO 90
+ ENDIF
+ 80 CONTINUE
+C...series will not converge for the currently prescribed tolerance
+C with the currently prescribed maximum number of terms MAXN:
+C Exit without computing derivative
+ NOCONV=.TRUE.
+ GOTO 999
+ 90 CONTINUE
+C Compute the derivative of exponential map
+ CALL RVZERO(DEXPX,NDIM4)
+ DO 150 I=1,NDIM
+ DO 140 J=1,NDIM
+ DO 130 K=1,NDIM
+ DO 120 L=1,NDIM
+ DO 110 N=1,NMAX
+ DO 100 M=1,N
+ DEXPX(I,J,K,L)=DEXPX(I,J,K,L)+
+ 1 R1DFAC(N)*XMATX(I,K,M-1)*XMATX(L,J,N-M)
+ 100 CONTINUE
+ 110 CONTINUE
+ 120 CONTINUE
+ 130 CONTINUE
+ 140 CONTINUE
+ 150 CONTINUE
+C
+ 999 CONTINUE
+ RETURN
+ END
+
+
+
+
+ SUBROUTINE RVZERO
+ 1 ( V ,N )
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ DIMENSION V(N)
+ DATA R0/0.0D0/
+C***********************************************************************
+C INITIALISES TO ZERO A DOUBLE PRECISION ARRAY OF DIMENSION N
+C***********************************************************************
+ DO 10 I=1,N
+ V(I)=R0
+ 10 CONTINUE
+ RETURN
+ END
+
+
+ DOUBLE PRECISION FUNCTION SCAPRD(U ,V ,N )
+C
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ DIMENSION U(N), V(N)
+ DATA R0 / 0.0D0 /
+C***********************************************************************
+C SCALAR PRODUCT OF DOUBLE PRECISION VECTORS U AND V OF DIMENSION N
+C***********************************************************************
+ SCAPRD=R0
+ DO 10 I=1,N
+ SCAPRD=SCAPRD+U(I)*V(I)
+ 10 CONTINUE
+ RETURN
+ END
+
+
+
+
+ SUBROUTINE norm(dnorm,MAT,m,n)
+ REAL*8 MAT(m,n),dnorm
+ integer ii,jj,m,n
+
+ dnorm=0d0
+ DO,ii=1,m
+ do,jj=1,n
+ dnorm=dnorm+MAT(ii,jj)*MAT(ii,jj)
+ enddo
+ enddo
+ dnorm=dsqrt(dnorm)
+ RETURN
+ END
+
+ subroutine TENS3333_sym(A,A66)
+ real*8 A(3,3,3,3)
+ real*8 a66(6,6)
+ integer indexi(6),indexj(6),iflag
+ integer indexii(3,3)
+ integer ii,jj
+
+
+ do,ii=1,3
+ indexi(ii)=ii
+ indexj(ii)=ii
+ enddo
+
+ indexi(4)=1
+ indexj(4)=2
+ indexi(5)=1
+ indexj(5)=3
+ indexi(6)=2
+ indexj(6)=3
+ indexii(1,1)=1
+ indexii(2,2)=2
+ indexii(3,3)=3
+ indexii(1,2)=4
+ indexii(1,3)=5
+ indexii(2,3)=6
+ indexii(2,1)=4
+ indexii(3,1)=5
+ indexii(3,2)=6
+
+ do,ii=1,6
+ do,jj=1,6
+
+ a66(ii,jj) = A( indexi(ii),indexj(ii),indexi(jj),indexj(jj))
+
+ enddo
+ enddo
+
+ RETURN
+ 100 FORMAT( 6(E16.8,' '))
+ END
+
+
+
+ subroutine TENS3333(A,A99)
+ real*8 A(3,3,3,3)
+ real*8 a99(9,9)
+ integer indexi(9),indexj(9),iflag
+ integer indexii(3,3)
+
+ do,ii=1,3
+ indexi(ii)=ii
+ indexj(ii)=ii
+ enddo
+
+ indexi(4)=1
+ indexj(4)=2
+ indexi(5)=1
+ indexj(5)=3
+ indexi(6)=2
+ indexj(6)=3
+ indexi(7)=2
+ indexj(7)=1
+ indexi(8)=3
+ indexj(8)=1
+ indexi(9)=3
+ indexj(9)=2
+
+ indexii(1,1)=1
+ indexii(2,2)=2
+ indexii(3,3)=3
+ indexii(1,2)=4
+ indexii(2,1)=7
+ indexii(1,3)=5
+ indexii(3,1)=8
+ indexii(2,3)=6
+ indexii(3,2)=9
+ do,ii=1,9
+ do,jj=1,9
+ a99(ii,jj)=A(indexi(ii),indexj(ii),indexi(jj),indexj(jj))
+ enddo
+ enddo
+
+ return
+ 100 FORMAT( 9(E16.8,' '))
+ end
+
+ integer function kroneker(i,j)
+ integer i,j
+ kroneker=0
+ if(i.EQ.j) kroneker=1
+ return
+ end
+
+
+ SUBROUTINE ROTATE_TENS3(AROT,A,R)
+c INPUTS: A (second order tensor, in conf 2, i.e. crystal axes), R (rot matrix from 1 to 2)
+C OUTPUT: AROT (second order tensor, in conf 1, i.e. XYZ)
+C AROT=RT A R !!!
+C
+c A given in system ei
+C ei*=Aei
+C Arot=R^T A R
+ REAL*8 A(3,3),R(3,3),AROT(3,3)
+ REAL*8 aux33(3,3)
+
+ do,i=1,3
+ do,j=1,3
+ aux33(i,j)=0d0
+ do,k=1,3
+ aux33(i,j)=R(k,i)*A(k,j)+aux33(i,j)
+ enddo
+ enddo
+ enddo
+ do,i=1,3
+ do,j=1,3
+ arot(i,j)=0d0
+ do,k=1,3
+ arot(i,j)=aux33(i,k)*R(k,j)+arot(i,j)
+ enddo
+ enddo
+ enddo
+ return
+ end
+
+ Subroutine LUDCMP_IMDEACP(A,N,INDX,D,CODE)
+ PARAMETER(NMAX=100,TINY=1d-16)
+ REAL*8 AMAX,DUM, SUM, A(N,N),VV(NMAX)
+ INTEGER CODE, D, INDX(N)
+
+ D=1; CODE=0
+
+ DO I=1,N
+ AMAX=0.d0
+ DO J=1,N
+ IF (DABS(A(I,J)).GT.AMAX) AMAX=DABS(A(I,J))
+ END DO ! j loop
+ IF(AMAX.LT.TINY) THEN
+ CODE = 1
+ RETURN
+ END IF
+ VV(I) = 1.d0 / AMAX
+ END DO ! i loop
+
+ DO J=1,N
+ DO I=1,J-1
+ SUM = A(I,J)
+ DO K=1,I-1
+ SUM = SUM - A(I,K)*A(K,J)
+ END DO ! k loop
+ A(I,J) = SUM
+ END DO ! i loop
+ AMAX = 0.d0
+ DO I=J,N
+ SUM = A(I,J)
+ DO K=1,J-1
+ SUM = SUM - A(I,K)*A(K,J)
+ END DO ! k loop
+ A(I,J) = SUM
+ DUM = VV(I)*DABS(SUM)
+ IF(DUM.GE.AMAX) THEN
+ IMAX = I
+ AMAX = DUM
+ END IF
+ END DO ! i loop
+
+ IF(J.NE.IMAX) THEN
+ DO K=1,N
+ DUM = A(IMAX,K)
+ A(IMAX,K) = A(J,K)
+ A(J,K) = DUM
+ END DO ! k loop
+ D = -D
+ VV(IMAX) = VV(J)
+ END IF
+
+ INDX(J) = IMAX
+ IF(DABS(A(J,J)) < TINY) A(J,J) = TINY
+
+ IF(J.NE.N) THEN
+ DUM = 1.d0 / A(J,J)
+ DO I=J+1,N
+ A(I,J) = A(I,J)*DUM
+ END DO ! i loop
+ END IF
+ END DO ! j loop
+
+ RETURN
+ END subroutine LUDCMP_IMDEACP
+
+c
+ SUBROUTINE gauss_3(AA,BB,xx,nmax,n,iflag)
+ integer iflag,d
+ real*8 AA(nmax,nmax),bb(nmax),xx(nmax)
+ real*8 A(n,n),b(n)
+ integer n,nmax,i,j,indx(n)
+ det=0.
+
+ do,i=1,n
+ b(i)=bb(i)
+ do,j=1,n
+ A(i,j)=AA(i,j)
+ enddo
+ enddo
+ call ludcmp_IMDEACP(a,n,indx,d,iflag)
+ call LUBKSB_IMDEACP(A,N,INDX,B)
+ do,i=1,n
+ xx(i)=b(i)
+ enddo
+ return
+ end
+
+ Subroutine LUBKSB_IMDEACP(A,N,INDX,B)
+ REAL*8 SUM, A(N,N),B(N)
+ INTEGER INDX(N)
+
+ II = 0
+
+ DO I=1,N
+ LL = INDX(I)
+ SUM = B(LL)
+ B(LL) = B(I)
+ IF(II.NE.0) THEN
+ DO J=II,I-1
+ SUM = SUM - A(I,J)*B(J)
+ END DO ! j loop
+ ELSE IF(SUM.NE.0.d0) THEN
+ II = I
+ END IF
+ B(I) = SUM
+ END DO ! i loop
+
+ DO I=N,1,-1
+ SUM = B(I)
+ IF(I < N) THEN
+ DO J=I+1,N
+ SUM = SUM - A(I,J)*B(J)
+ END DO ! j loop
+ END IF
+ B(I) = SUM / A(I,I)
+ END DO ! i loop
+
+
+ RETURN
+ END subroutine LUBKSB_IMDEACP
+
+
+! end of file lu.f90
+
+
+ subroutine dzeros_vec(A,i)
+ real*8 A(i)
+ integer ii,i
+ do,ii=1,i
+ A(ii)=0d0
+ enddo
+ return
+ end
+
+ subroutine dzeros(A,i,j)
+ real*8 A(i,j)
+ integer i,j,ii,jj
+ do,ii=1,i
+ do,jj=1,j
+ A(ii,jj)=0d0
+ enddo
+ enddo
+ return
+ end
+
+ subroutine dzeros_tensor4(A,i,j,k,l)
+ real*8 A(i,j,k,l)
+ integer i,j,ii,jj
+ do,ii=1,i
+ do,jj=1,j
+ do,kk=1,k
+ do,ll=1,l
+ A(ii,jj,kk,ll)=0d0
+ enddo
+ enddo
+ enddo
+ enddo
+ return
+ end
+
+ subroutine dunit(A,i,j)
+ real*8 A(i,j)
+ integer i,j,ii,jj
+ do,ii=1,i
+ do,jj=1,j
+ A(ii,jj)=0d0
+ enddo
+ A(ii,ii)=1d0
+ enddo
+ return
+ end
+
+
+ subroutine min_location(min_value,min_loc,input_array,n)
+ implicit none
+
+ integer ii, n, min_loc
+ real*8 min_value, input_array(n)
+
+ min_value=input_array(1)
+ min_loc=1
+
+ do, ii=2,n
+ if(min_value.gt.input_array(ii)) then
+ min_value=input_array(ii)
+ min_loc=ii
+ endif
+ enddo
+
+ return
+ end
+
+ SUBROUTINE paralelization(init,init1,noel,numprocess,noelprocess)
+ implicit NONE
+ logical init(10),init1(10)
+ INTEGER NUMPROCESS
+ integer i,noel,noelprocess(10)
+c save init,init1
+ numprocess=1
+c write(*,*)' numprocess,',numprocess
+ DO,I=0,7
+ if(.not. init(numprocess+1).and.numprocess==I) THEN
+ noelprocess(numprocess+1)=noel
+ write(*,*)'noel cluster node number and element',I,noel
+ init(numprocess+1)=.true.
+ endif
+ enddo
+
+ do,i=0,7
+ if(.not. init1(numprocess+1).and.numprocess==I) THEN
+ if((noel.NE.noelprocess(numprocess+1))) then
+ write(*,*)'noel waiting node number and element ',noel
+ call sleep(1)
+ endif
+ endif
+ enddo
+ return
+ end
+
+ SUBROUTINE form_orient_MAT(props,orient_MAT,orient_MAT_tr)
+ implicit none
+ integer i
+ real*8 props(*)
+ real*8 orient1(3),orient2(3),orient3(3),dnorm
+ real*8 orient_MAT(3,3),orient_mat_tr(3,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
+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
+ enddo
+ CALL pvect(orient3,orient1,orient2)
+ call norm(dnorm,orient3,3,1)
+ do,i=1,3
+ orient3(i)=orient3(i)/dnorm
+ enddo
+ do,i=1,3
+ orient_MAT(i,1)=orient1(i) ! Orient_MAT is the rotation matrix from crystal axis to reference configuration
+ orient_MAT(i,2)=orient2(i) !
+ orient_MAT(i,3)=orient3(i) ! Orient_MAT u_crys= u_cart
+ enddo
+
+ call transpose(orient_MAT_tr,orient_mat,3,3)
+c write(*,*)'I leave lybrary_crysplas_CAPSUL cm3Libraries--------'
+ return
+ end
+
+ SUBROUTINE form_Mjacob(dtime,DFGRD_ELAS_PRED0_inv,Mjacob)
+ implicit none
+ real*8 dtime,DFGRD_ELAS_PRED0_inv(3,3),Mjacob(3,3,3,3)
+ real*8 term
+ integer ii,jj,kk,ll,kroneker
+
+ CALL dzeros_tensor4(Mjacob,3,3,3,3)
+
+ do ii=1,3
+ do jj=1,3
+ do kk=1,3
+ do ll=1,3
+ term=DFGRD_ELAS_PRED0_inv(ii,kk)*
+ 1 kroneker(jj,ll)
+ Mjacob(ii,jj,kk,ll)=-(1/dtime)*term
+ enddo
+ enddo
+ enddo
+ enddo
+ return
+ end
+
+ SUBROUTINE form_partial_sigma_F(DFGRD_ELAS_PRED,stiff_orient,
+ 1 partial_sigma)
+ implicit none
+ REAL*8 DFGRD_ELAS_PRED(3,3),stiff_orient(3,3,3,3)
+ REAL*8 partial_sigma(3,3,3,3),aux3333(3,3,3,3)
+ integer ii,jj,nm,nn,kk,ll,kroneker
+
+ CALL dzeros_tensor4(partial_sigma,3,3,3,3)
+
+c partial Epsilon/Fe (rs,ss,nm,nn)
+ do,ii=1,3
+ do,jj=1,3
+ do,nm=1,3
+ do,nn=1,3
+ aux3333(ii,jj,nm,nn)=
+ 1 .5d0*(kroneker(ii,nn)*
+ 1 DFGRD_ELAS_PRED(nm,jj)+
+ 1 kroneker(jj,nn)*
+ 1 DFGRD_ELAS_PRED(nm,ii))
+ enddo
+ enddo
+ enddo
+ enddo
+
+
+
+c partial sigmaij/Fmn
+
+ do,ii=1,3
+ do,jj=1,3
+ do,nm=1,3
+ do,nn=1,3
+ do,kk=1,3
+ do,ll=1,3
+ partial_sigma(ii,jj,nm,nn)=
+ 1 partial_sigma(ii,jj,nm,nn)+
+ 1 stiff_orient(ii,jj,kk,ll)*
+ 1 aux3333(kk,ll,nm,nn)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ return
+ end
+
+
+
diff --git a/NonLinearSolver/materialLaw/mlaw.h b/NonLinearSolver/materialLaw/mlaw.h
index a682feea1f180ae1d371e81b9a86f0bda27d0ff7..27db4d202b2490936c134ab715658241fd6738f6 100644
--- a/NonLinearSolver/materialLaw/mlaw.h
+++ b/NonLinearSolver/materialLaw/mlaw.h
@@ -35,7 +35,7 @@ class materialLaw{
ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic,nonLocalDamageGursonThermoMechanics,
localDamageJ2SmallStrain,nonLocalDamageJ2SmallStrain,cluster,tfa,ANN, DMN, torchANN, NonlocalDamageTorchANN, LinearElecMagTherMech, LinearElecMagInductor, hyperviscoelastic,
GenericThermoMechanics, ElecMagGenericThermoMechanics, ElecMagInductor,
- nonlineartvm,nonlinearTVE,nonlinearTVP,vevpmfh, VEVPUMAT, Hill48,nonlinearTVEnonlinearTVP};
+ nonlineartvm,nonlinearTVE,nonlinearTVP,vevpmfh, VEVPUMAT, IMDEACPUMAT, Hill48,nonlinearTVEnonlinearTVP};
protected :
diff --git a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d33bd712b31eacd069ee85c0e8725bc49ec9c2bc
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
@@ -0,0 +1,126 @@
+//
+// C++ Interface: material law
+//
+// Description: UMAT interface for IMDEACP model
+//
+//
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include <math.h>
+#include "MInterfaceElement.h"
+#include <fstream>
+#include "mlawIMDEACPUMAT.h"
+#include "STensorOperations.h"
+
+#if !defined(F77NAME)
+#define F77NAME(x) (x##_)
+#endif
+
+extern "C" {
+ void F77NAME(umat_imdeacp)(double *STRESS, double *STATEV, double *DDSDDE, double *SSE, double *SPD, double *SCD,
+ double *RPL, double *DDSDDT, double *DRPLDE, double *DRPLDT,
+ double *STRAN, double *DSTRAN, double *TIM, double *DTIME, double *TEMP, double *DTEMP, double *PREDEF, double *DPRED, const char *CMNAME,
+ int *NDI, int*NSHR, int* NTENS, int*NSTATV, double *PROPS, int *NPROPS, double *COORDS, double *DROT, double *PNEWDT,
+ double *CELENT, double *DFGRD0, double *DFGRD1, int *NOEL, int *NPT, int *LAYER, int *KSPT, int *KSTEP, int *KINC);
+}
+
+
+
+mlawIMDEACPUMAT::mlawIMDEACPUMAT(const int num, const char *propName) : mlawUMATInterface(num,0.,propName)
+{
+ //should depend on the umat: here we read the grains.inp file with 3 euler angle, 2 parameters equal to 1, and the number of internal variables (23)
+ std::ifstream in(propName);
+ if(!in) Msg::Error("Cannot open the file %s! Maybe is missing ",propName);
+ std::string line;
+ std::getline(in, line,'\n');
+ std::getline(in, line,'\n');
+ 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,',') and i<3)
+ {
+ 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 C1111 = props[0];
+ double C1122 = props[1];
+ double C1212 = props[2];
+ _rho = 0.;
+
+ 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)
+{
+
+}
+
+mlawIMDEACPUMAT& mlawIMDEACPUMAT::operator=(const materialLaw &source)
+{
+ mlawUMATInterface::operator=(source);
+ const mlawIMDEACPUMAT* src =static_cast<const mlawIMDEACPUMAT*>(&source);
+ return *this;
+};
+
+mlawIMDEACPUMAT::~mlawIMDEACPUMAT()
+{
+
+}
+void mlawIMDEACPUMAT::createIPState(IPIMDEACPUMAT *ivi, IPIMDEACPUMAT *iv1, IPIMDEACPUMAT *iv2) const
+{
+ mlawUMATInterface::createIPState(ivi, iv1, iv2);
+}
+
+void mlawIMDEACPUMAT::createIPVariable(IPIMDEACPUMAT *ipv, bool hasBodyForce, const MElement *ele,const int nbFF,const IntPt *GP, const int gpt) const
+{
+ mlawUMATInterface::createIPVariable(ipv, hasBodyForce, ele, nbFF,GP, gpt);
+}
+
+
+void mlawIMDEACPUMAT::callUMAT(double *stress, double *statev, double **ddsdde, double &sse, double &spd, double &scd, double &rpl,
+ double *ddsddt, double *drplde, double &drpldt, double *stran, double *dtsran,
+ double *tim, double timestep, double temperature, double deltaTemperature, double *predef, double *dpred,
+ const char *CMNAME, int &ndi, int &nshr, int tensorSize,
+ int statevSize, double *prs, int matPropSize, double *coords, double **dRot,
+ double &pnewdt, double &celent, double **F0, double **F1,
+ int &noel, int &npt, int &layer, int &kspt, int &kstep, int &kinc) const
+{
+
+
+ double dRtmp[9];
+ double F0tmp[9];
+ double F1tmp[9];
+
+ double ddsddetmp[36];
+
+ convert3x3To9((const double**)dRot, dRtmp);
+ convert3x3To9((const double**)F0, F0tmp);
+ convert3x3To9((const double**)F1, F1tmp);
+
+ convert6x6To36((const double**)ddsdde,ddsddetmp);
+
+ F77NAME(umat_imdeacp)(stress, statev, ddsddetmp, &sse, &spd, &scd, &rpl,
+ ddsddt, drplde, &drpldt, stran, dtsran,
+ tim, ×tep, &temperature, &deltaTemperature, predef, dpred,
+ CMNAME, &ndi, &nshr, &tensorSize, &statevSize, prs, &matPropSize, coords, dRtmp,
+ &pnewdt, &celent, F0tmp, F1tmp,
+ &noel, &npt, &layer, &kspt, &kstep, &kinc);
+
+ convert9To3x3((const double*)dRtmp,dRot);
+ convert36To6x6((const double*)ddsddetmp,ddsdde);
+}
+
diff --git a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h
new file mode 100644
index 0000000000000000000000000000000000000000..2346dd18a392ce0da9945f908e4baf83756c2c0d
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h
@@ -0,0 +1,109 @@
+//
+// C++ Interface: material law
+//
+// Description: UMAT interface for VEVP model
+//
+//
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef MLAWVEVPUMAT_H_
+#define MLAWVEVPUMAT_H_
+#include "mlawUMATInterface.h"
+#include "ipVEVPUMAT.h"
+
+class mlawVEVPUMAT : public mlawUMATInterface
+{
+
+ public:
+ mlawVEVPUMAT(const int num, const char *propName);
+
+ #ifndef SWIG
+ mlawVEVPUMAT(const mlawVEVPUMAT &source);
+ mlawVEVPUMAT& operator=(const materialLaw &source);
+ virtual ~mlawVEVPUMAT();
+ virtual materialLaw* clone() const {return new mlawVEVPUMAT(*this);};
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{}; // do nothing
+ // function of materialLaw
+ virtual matname getType() const{return materialLaw::VEVPUMAT;}
+
+ virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const
+ {
+ Msg::Error("Cannot be called");
+ }
+
+ virtual void createIPState(IPVEVPUMAT *ivi, IPVEVPUMAT *iv1, IPVEVPUMAT *iv2) const;
+ virtual void createIPVariable(IPVEVPUMAT *ipv,bool hasBodyForce, const MElement *ele,const int nbFF, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
+
+ virtual void callUMAT(double *stress, double *statev, double **ddsdde, double &sse, double &spd, double &scd, double &rpl,
+ double *ddsddt, double *drplde, double &drpldt, double *stran, double *dtsran,
+ double *tim, double timestep, double temperature, double deltaTemperature, double *predef, double *dpred,
+ const char *CMNAME, int &ndi, int &nshr, int tensorSize,
+ int statevSize, double *prs, int matPropSize, double *coords, double **dRot,
+ double &pnewdt, double &celent, double **F0, double **F1,
+ int &noel, int &npt, int &layer, int &kspt, int &kstep, int &kinc) const;
+
+ virtual const char* getCMNAME() const {return "VEVP";}
+ virtual double getDensity() const {return _rho;}
+
+ #endif // SWIG
+};
+
+#endif // MLAWVEVPUMAT_H_//
+// C++ Interface: material law
+//
+// Description: UMAT interface for IMDEACP model
+//
+//
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef MLAWIMDEACPUMAT_H_
+#define MLAWIMDEACPUMAT_H_
+#include "mlawUMATInterface.h"
+#include "ipIMDEACPUMAT.h"
+
+class mlawIMDEACPUMAT : public mlawUMATInterface
+{
+
+ public:
+ mlawIMDEACPUMAT(const int num, const char *propName);
+
+ #ifndef SWIG
+ mlawIMDEACPUMAT(const mlawIMDEACPUMAT &source);
+ mlawIMDEACPUMAT& operator=(const materialLaw &source);
+ virtual ~mlawIMDEACPUMAT();
+ virtual materialLaw* clone() const {return new mlawIMDEACPUMAT(*this);};
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{}; // do nothing
+ // function of materialLaw
+ virtual matname getType() const{return materialLaw::IMDEACPUMAT;}
+
+ virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const
+ {
+ Msg::Error("Cannot be called");
+ }
+
+ virtual void createIPState(IPIMDEACPUMAT *ivi, IPIMDEACPUMAT *iv1, IPIMDEACPUMAT *iv2) const;
+ virtual void createIPVariable(IPIMDEACPUMAT *ipv,bool hasBodyForce, const MElement *ele,const int nbFF, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
+
+ virtual void callUMAT(double *stress, double *statev, double **ddsdde, double &sse, double &spd, double &scd, double &rpl,
+ double *ddsddt, double *drplde, double &drpldt, double *stran, double *dtsran,
+ double *tim, double timestep, double temperature, double deltaTemperature, double *predef, double *dpred,
+ const char *CMNAME, int &ndi, int &nshr, int tensorSize,
+ int statevSize, double *prs, int matPropSize, double *coords, double **dRot,
+ double &pnewdt, double &celent, double **F0, double **F1,
+ int &noel, int &npt, int &layer, int &kspt, int &kstep, int &kinc) const;
+
+ virtual const char* getCMNAME() const {return "IMDEACP";}
+ virtual double getDensity() const {return _rho;}
+
+ #endif // SWIG
+};
+
+#endif // MLAWIMDEACPUMAT_H_
diff --git a/NonLinearSolver/materialLaw/mlawUMATInterface.cpp b/NonLinearSolver/materialLaw/mlawUMATInterface.cpp
index 048d7991db600f78f72c15a6e20d4750ab2debbd..3de5f5bbaa84e1b6ac37f35cc723da3aee838c89 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/CMakeLists.txt b/dG3D/benchmarks/imdeaCP/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3667d9747b7c516dd7833cc4e59989e55b267c01
--- /dev/null
+++ b/dG3D/benchmarks/imdeaCP/CMakeLists.txt
@@ -0,0 +1,11 @@
+# test file
+
+set(PYFILE cylinder.py)
+
+set(FILES2DELETE
+ disp*.msh
+ stress*.msh
+ *.csv
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")
diff --git a/dG3D/benchmarks/imdeaCP/crystal.prop b/dG3D/benchmarks/imdeaCP/crystal.prop
new file mode 100644
index 0000000000000000000000000000000000000000..e27c272d1aa9c3f6e914e24029d2aea6d33010fb
--- /dev/null
+++ b/dG3D/benchmarks/imdeaCP/crystal.prop
@@ -0,0 +1,29 @@
+# Name: INC718
+# C11, C12,C44
+259.6D9,179D9,109.6D9
+#Viscoplastic law: gamma_0, exponent
+2.42D-3,0.017
+# Set of slip systems, Total number of systems
+1,12
+# 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
+ 1 , -1 , 1 , 1 , 1 , 0, 1
+ 1 , -1 , 1 , 0 , 1 , 1 , 1
+ 1 , -1 , 1 , -1 , 0 , 1 , 1
+ -1 ,-1 , 1 , 1 , -1 , 0 , 1
+ -1 ,-1 , 1 , 0 , 1 , 1 , 1
+ -1 ,-1 , 1 , 1 , 0 , 1 , 1
+ -1 ,1 , 1 , 1, 1 , 0 , 1
+ -1 ,1 , 1 , 0 , -1 ,1 , 1
+ -1 ,1 , 1 , 1 , 0 , 1 , 1
+# q(alfa,alfa), q(alfa,beta), q (alfa, gamma), q(alfa, eta), q (beta, beta), q(beta, gamma), q(beta,eta), q(gamma, gamma), q (gamma, eta), q (eta, eta)
+1.
+# Single crystal behavior based on Voce: tau0,taus,h0 (.9tau0,CRSS,fixed) x0.25
+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.geo b/dG3D/benchmarks/imdeaCP/cube.geo
new file mode 100644
index 0000000000000000000000000000000000000000..9654d90c5497bf8a1ff446e39d1b9ad26a25086c
--- /dev/null
+++ b/dG3D/benchmarks/imdeaCP/cube.geo
@@ -0,0 +1,51 @@
+// Cube.geo
+mm=1.0; // Units
+n=1; // Number of layers (in thickness / along z)
+m = 1; // Number of element + 1 along y
+p = 3; // Number of element + 1 along x (min 3 to have a crack)
+L=2.*mm; // Cube lenght
+sl1=L/n;
+
+// Face z = 0 (point in anti-clockwise order)
+Point(1)={0,0,0,sl1};
+Point(2)={L,0,0,sl1};
+Point(3)={L,L,0,sl1};
+Point(4)={0,L,0,sl1};
+Line(1)={1,2};
+Line(2)={2,3};
+Line(3)={3,4};
+Line(4)={4,1};
+Line Loop(5) = {3, 4, 1, 2};
+Plane Surface(6) = {5};
+
+// Cube extrusion
+Extrude {0, 0, L} {
+ Surface{6}; Layers{n}; Recombine;
+}
+
+// Physical entities
+ // Volume
+Physical Volume(29) = {1};
+ // Surface
+Physical Surface(30) = {19}; // Left face x = 0
+Physical Surface(31) = {27}; // Right face x = L
+Physical Surface(32) = {6}; // Back face z = 0
+Physical Surface(33) = {28}; // Front face z = L
+Physical Surface(34) = {15}; // Down face y = 0
+Physical Surface(35) = {23}; // Up face y = L
+
+// Mesh
+Transfinite Line {2, 4} = m Using Progression 1;
+Transfinite Line {1, 3} = p Using Progression 1;
+Transfinite Surface {6};
+Recombine Surface{6};
+
+
+// To save cube elongation
+Physical Point(41) = {1}; // Point on left face (0,0,0)
+Physical Point(42) = {4}; // Point on left face (0,L,0)
+Physical Point(43) = {2}; // Point on right face (L,0,0)
+Physical Point(44) = {3}; // Point on right face (L,L,0)
+
+
+
diff --git a/dG3D/benchmarks/imdeaCP/cube.msh b/dG3D/benchmarks/imdeaCP/cube.msh
new file mode 100644
index 0000000000000000000000000000000000000000..a7382df16647c759d6e8e08f9a57c54f24e473bd
--- /dev/null
+++ b/dG3D/benchmarks/imdeaCP/cube.msh
@@ -0,0 +1,109 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+8 12 6 1
+1 0 0 0 1 41
+2 2 0 0 1 43
+3 2 2 0 1 44
+4 0 2 0 1 42
+5 2 2 2 0
+6 0 2 2 0
+10 0 0 2 0
+14 2 0 2 0
+1 0 0 0 2 0 0 0 2 1 -2
+2 2 0 0 2 2 0 0 2 2 -3
+3 0 2 0 2 2 0 0 2 3 -4
+4 0 0 0 0 2 0 0 2 4 -1
+8 0 2 2 2 2 2 0 2 5 -6
+9 0 0 2 0 2 2 0 2 6 -10
+10 0 0 2 2 0 2 0 2 10 -14
+11 2 0 2 2 2 2 0 2 14 -5
+13 2 2 0 2 2 2 0 2 3 -5
+14 0 2 0 0 2 2 0 2 4 -6
+18 0 0 0 0 0 2 0 2 1 -10
+22 2 0 0 2 0 2 0 2 2 -14
+6 0 0 0 2 2 0 1 32 4 3 4 1 2
+15 0 2 0 2 2 2 1 34 4 3 14 -8 -13
+19 0 0 0 0 2 2 1 30 4 4 18 -9 -14
+23 0 0 0 2 0 2 1 35 4 1 22 -10 -18
+27 2 0 0 2 2 2 1 31 4 2 13 -11 -22
+28 0 0 2 2 2 2 1 33 4 8 9 10 11
+1 0 0 0 2 2 2 1 29 6 -6 28 15 19 23 27
+$EndEntities
+$Nodes
+19 12 1 12
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+2 0 0
+0 3 0 1
+3
+2 2 0
+0 4 0 1
+4
+0 2 0
+0 5 0 1
+5
+2 2 2
+0 6 0 1
+6
+0 2 2
+0 10 0 1
+7
+0 0 2
+0 14 0 1
+8
+2 0 2
+1 1 0 1
+9
+0.9999999999973842 0 0
+1 3 0 1
+10
+1.000000000004119 2 0
+1 8 0 1
+11
+1.000000000004119 2 2
+1 10 0 1
+12
+0.9999999999973842 0 2
+2 6 0 0
+2 15 0 0
+2 19 0 0
+2 23 0 0
+2 27 0 0
+2 28 0 0
+3 1 0 0
+$EndNodes
+$Elements
+11 16 1 16
+0 1 15 1
+1 1
+0 2 15 1
+2 2
+0 3 15 1
+3 3
+0 4 15 1
+4 4
+2 6 3 2
+5 3 10 9 2
+6 10 4 1 9
+2 15 3 2
+7 3 10 11 5
+8 10 4 6 11
+2 19 3 1
+9 4 1 7 6
+2 23 3 2
+10 1 9 12 7
+11 9 2 8 12
+2 27 3 1
+12 2 3 5 8
+2 28 3 2
+13 5 11 12 8
+14 11 6 7 12
+3 1 5 2
+15 3 10 9 2 5 11 12 8
+16 10 4 1 9 11 6 7 12
+$EndElements
diff --git a/dG3D/benchmarks/imdeaCP/cube.py b/dG3D/benchmarks/imdeaCP/cube.py
new file mode 100644
index 0000000000000000000000000000000000000000..3193d7dbeb9c6312b81b19832925461c6b4631fa
--- /dev/null
+++ b/dG3D/benchmarks/imdeaCP/cube.py
@@ -0,0 +1,88 @@
+#coding-Utf-8-*-
+
+from gmshpy import *
+from dG3Dpy import*
+from math import*
+
+#script to launch PBC problem with a python script
+
+# material law
+lawnum = 11 # unique number of law
+
+# material law
+lawnum = 11 # unique number of law
+law1 = IMDEACPUMATDG3DMaterialLaw(lawnum,'crystal.prop');
+#lawPert1 = dG3DMaterialLawWithTangentByPerturbation(law1, 1e-8)
+
+
+# geometry
+meshfile="cube.msh" # name of mesh file
+
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1 = 100
+# creation of part Domain
+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)
+
+# solver
+sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 100 # number of step (used only if soltype=1)
+ftime =1. # Final time (used only if soltype=1)
+tol=1.e-8 # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+
+
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+#mysolver.addMaterialLaw(lawPert1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+
+mysolver.displacementBC('Face',32,2,0.)
+mysolver.displacementBC('Face',34,1,0.)
+mysolver.displacementBC('Face',30,0,0.)
+
+#mysolver.displacementBC('Face',1,0,0)
+#mysolver.displacementBC('Face',1,1,0)
+mysolver.displacementBC('Face',31,0,0.2)
+
+# build view
+mysolver.internalPointBuildView("Strain_xx",IPField.STRAIN_XX, 1, 1);
+mysolver.internalPointBuildView("Strain_yy",IPField.STRAIN_YY, 1, 1);
+mysolver.internalPointBuildView("Strain_zz",IPField.STRAIN_ZZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xy",IPField.STRAIN_XY, 1, 1);
+mysolver.internalPointBuildView("Strain_yz",IPField.STRAIN_YZ, 1, 1);
+mysolver.internalPointBuildView("Strain_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("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
+mysolver.internalPointBuildView("pression",IPField.PRESSION,1,1)
+mysolver.internalPointBuildView("plastic possion ratio",IPField.PLASTIC_POISSON_RATIO,1,1)
+
+#archiving
+mysolver.archivingForceOnPhysicalGroup('Face',30,0)
+mysolver.archivingForceOnPhysicalGroup('Face',30,1)
+mysolver.archivingForceOnPhysicalGroup('Face',30,2)
+mysolver.archivingNodeDisplacement(43,0,1)
+mysolver.archivingNodeDisplacement(43,1,1)
+mysolver.archivingNodeDisplacement(43,2,1)
+# solve
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-5.574939e+09,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-3)
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index be4664fe98b08b3c10d8ed08f9fb9dd11b2ae6d3..cd7f56808858eb9ad7782cf44aef4ebeea4c415d 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -4296,6 +4296,71 @@ void VEVPUMATDG3DIPVariable::restart()
return;
}
+//
+IMDEACPUMATDG3DIPVariable::IMDEACPUMATDG3DIPVariable(int _nsdv, double size, const bool createBodyForceHO, const bool oninter):
+ dG3DIPVariable(createBodyForceHO, oninter,0)
+{
+ _ipv = new IPIMDEACPUMAT(_nsdv,size);
+}
+
+IMDEACPUMATDG3DIPVariable::IMDEACPUMATDG3DIPVariable(const IMDEACPUMATDG3DIPVariable &source) :
+ dG3DIPVariable(source)
+{
+ _ipv = NULL;
+ if (source.getIPIMDEACPUMAT() != NULL){
+ _ipv = new IPIMDEACPUMAT((int)source.getIPIMDEACPUMAT()->getNsdv(),(double)source.getIPIMDEACPUMAT()->elementSize());
+ _ipv->operator= (*dynamic_cast<const IPVariable*>(source.getIPIMDEACPUMAT()));
+ }
+
+}
+
+IMDEACPUMATDG3DIPVariable& IMDEACPUMATDG3DIPVariable::operator=(const IPVariable &source)
+{
+ dG3DIPVariable::operator=(source);
+ const IMDEACPUMATDG3DIPVariable* src = dynamic_cast<const IMDEACPUMATDG3DIPVariable*>(&source);
+ if(src != NULL)
+ {
+ if (src->getIPIMDEACPUMAT() != NULL){
+ if (_ipv != NULL){
+ _ipv->operator=(*dynamic_cast<const IPVariable*>(src->getIPIMDEACPUMAT()));
+ }
+ else{
+ _ipv = dynamic_cast<IPIMDEACPUMAT*>(src->getIPIMDEACPUMAT()->clone());
+ }
+ }
+ else{
+ if (_ipv != NULL){
+ delete _ipv;
+ _ipv = NULL;
+ }
+ }
+ }
+ return *this;
+}
+double IMDEACPUMATDG3DIPVariable::get(const int i) const
+{
+ double val = dG3DIPVariable::get(i);
+ if (val == 0)
+ val = getIPIMDEACPUMAT()->get(i);
+ return val;
+}
+double IMDEACPUMATDG3DIPVariable::defoEnergy() const
+{
+ return getIPIMDEACPUMAT()->defoEnergy();
+}
+double IMDEACPUMATDG3DIPVariable::plasticEnergy() const
+{
+ return getIPIMDEACPUMAT()->plasticEnergy();
+}
+
+
+void IMDEACPUMATDG3DIPVariable::restart()
+{
+ dG3DIPVariable::restart();
+ if (_ipv != NULL)
+ restartManager::restart(_ipv);
+ return;
+}
//
gursonUMATDG3DIPVariable::gursonUMATDG3DIPVariable(int _nsdv, double size, const bool createBodyForceHO, const bool oninter):
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 12912f1ca4d5cc9eabd4d72d3708da39e5e26d32..44eef4a5f6e7f5a33d4b73ee034a99d02ff1d824 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -35,6 +35,7 @@
#include "ipCrystalPlasticity.h"
#include "ipGursonUMAT.h"
#include "ipVEVPUMAT.h"
+#include "ipIMDEACPUMAT.h"
#include "ipViscoelastic.h"
#include "ipAnIsotropicElecTherMech.h"
#include "ipLinearElecTherMech.h"
@@ -3604,6 +3605,51 @@ class VEVPUMATDG3DIPVariable : public dG3DIPVariable
virtual void restart();
};
+class IMDEACPUMATDG3DIPVariable : public dG3DIPVariable
+{
+
+ protected:
+ IPIMDEACPUMAT *_ipv;
+
+ public:
+ IMDEACPUMATDG3DIPVariable(int _nsdv, double size, const bool createBodyForceHO, const bool oninter);
+ IMDEACPUMATDG3DIPVariable(const IMDEACPUMATDG3DIPVariable &source);
+ virtual IMDEACPUMATDG3DIPVariable& operator=(const IPVariable &source);
+
+ /* specific function */
+ IPIMDEACPUMAT* getIPIMDEACPUMAT(){return _ipv;}
+ const IPIMDEACPUMAT* getIPIMDEACPUMAT() const{return _ipv;}
+ virtual ~IMDEACPUMATDG3DIPVariable()
+ {
+ if (_ipv) delete _ipv;
+ _ipv = NULL;
+ }
+
+ virtual void setLocation(const IPVariable::LOCATION loc) {
+ dG3DIPVariable::setLocation(loc);
+ if (_ipv)
+ _ipv->setLocation(loc);
+ };
+
+ virtual IPVariable* getInternalData() {return _ipv;};
+ virtual const IPVariable* getInternalData() const {return _ipv;};
+
+ virtual bool dissipationIsActive() const {return _ipv->dissipationIsActive();};
+
+ virtual void blockDissipation(const bool fl){if (_ipv) _ipv->blockDissipation(fl);};
+ virtual bool dissipationIsBlocked() const {
+ if (_ipv) return _ipv->dissipationIsBlocked();
+ else return false;
+ }
+
+ virtual double get(const int i) const;
+ virtual double defoEnergy() const;
+ virtual double plasticEnergy() const;
+
+ virtual IPVariable* clone() const {return new IMDEACPUMATDG3DIPVariable(*this);};
+ virtual void restart();
+};
+
class gursonUMATDG3DIPVariable : public dG3DIPVariable
{
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index dcdbf9b1ae8108935b7657fb44943e85edc6aee1..bb1f96e5a58056f33e28d35c213f9be63d353855 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -7774,6 +7774,154 @@ void VEVPUMATDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
_vevplaw->setMacroSolver(sv);
}
};
+//
+
+//
+IMDEACPUMATDG3DMaterialLaw::IMDEACPUMATDG3DMaterialLaw(const int num, const char *propName) :
+ dG3DMaterialLaw(num,0.,true)
+{
+ _vevplaw = new mlawIMDEACPUMAT(num,propName);
+ double nu = _vevplaw->poissonRatio();
+ double mu = _vevplaw->shearModulus();
+ double E = 2.*mu*(1.+nu);
+ _rho = _vevplaw->getDensity();
+ fillElasticStiffness(E, nu, elasticStiffness);
+
+
+}
+
+IMDEACPUMATDG3DMaterialLaw::~IMDEACPUMATDG3DMaterialLaw()
+{
+ if (_vevplaw !=NULL) delete _vevplaw;
+ _vevplaw = NULL;
+};
+
+
+IMDEACPUMATDG3DMaterialLaw::IMDEACPUMATDG3DMaterialLaw(const IMDEACPUMATDG3DMaterialLaw &source) :
+ dG3DMaterialLaw(source)
+{
+ _vevplaw = NULL;
+ if (source._vevplaw != NULL){
+ _vevplaw = dynamic_cast<mlawIMDEACPUMAT*>(source._vevplaw->clone());
+ }
+
+}
+
+void IMDEACPUMATDG3DMaterialLaw::setTime(const double t,const double dtime){
+ dG3DMaterialLaw::setTime(t,dtime);
+ _vevplaw->setTime(t,dtime);
+ return;
+}
+
+materialLaw::matname IMDEACPUMATDG3DMaterialLaw::getType() const {return _vevplaw->getType();}
+
+bool IMDEACPUMATDG3DMaterialLaw::withEnergyDissipation() const {return _vevplaw->withEnergyDissipation();};
+
+void IMDEACPUMATDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ // check interface element
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele==NULL) inter=false;
+
+ int nsdv = _vevplaw->getNsdv();
+ double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
+
+ IPVariable* ipvi = new IMDEACPUMATDG3DIPVariable(nsdv,size,hasBodyForce,inter);
+ IPVariable* ipv1 = new IMDEACPUMATDG3DIPVariable(nsdv,size,hasBodyForce,inter);
+ IPVariable* ipv2 = new IMDEACPUMATDG3DIPVariable(nsdv,size,hasBodyForce,inter);
+
+ if(ips != NULL) delete ips;
+ ips = new IP3State(state_,ipvi,ipv1,ipv2);
+ _vevplaw->createIPState((static_cast <IMDEACPUMATDG3DIPVariable*> (ipvi))->getIPIMDEACPUMAT(),
+ (static_cast <IMDEACPUMATDG3DIPVariable*> (ipv1))->getIPIMDEACPUMAT(),
+ (static_cast <IMDEACPUMATDG3DIPVariable*> (ipv2))->getIPIMDEACPUMAT());
+
+}
+
+void IMDEACPUMATDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ if(ipv !=NULL) delete ipv;
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele == NULL) inter=false;
+ double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
+
+ ipv = new IMDEACPUMATDG3DIPVariable(_vevplaw->getNsdv(),size,hasBodyForce,inter);
+
+ IPIMDEACPUMAT * ipvnl = static_cast <IMDEACPUMATDG3DIPVariable*>(ipv)->getIPIMDEACPUMAT();
+ _vevplaw->createIPVariable(ipvnl,hasBodyForce, ele,nbFF_,GP,gpt);
+}
+
+
+double IMDEACPUMATDG3DMaterialLaw::soundSpeed() const{return _vevplaw->soundSpeed();}
+
+void IMDEACPUMATDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{
+ /* get ipvariable */
+ IMDEACPUMATDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+ const IMDEACPUMATDG3DIPVariable* ipvprev; //= static_cast <const nonLocalDamageDG3DIPVariable *> (ipvp);
+
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+
+ ipvcur = static_cast<IMDEACPUMATDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const IMDEACPUMATDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<IMDEACPUMATDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const IMDEACPUMATDG3DIPVariable*>(ipvp);
+ }
+ _vevplaw->checkInternalState(ipvcur->getInternalData(),ipvprev->getInternalData());
+
+};
+
+void IMDEACPUMATDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent)
+{
+ /* get ipvariable */
+ IMDEACPUMATDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+ const IMDEACPUMATDG3DIPVariable* ipvprev; //= static_cast <const nonLocalDamageDG3DIPVariable *> (ipvp);
+
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+
+ ipvcur = static_cast<IMDEACPUMATDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const IMDEACPUMATDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<IMDEACPUMATDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const IMDEACPUMATDG3DIPVariable*>(ipvp);
+ }
+
+ /* strain xyz */
+ const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
+ const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
+
+ /* data for J2 law */
+ IPIMDEACPUMAT* q1 = ipvcur->getIPIMDEACPUMAT();
+ const IPIMDEACPUMAT* q0 = ipvprev->getIPIMDEACPUMAT();
+
+ /* compute stress */
+ _vevplaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
+ stiff, NULL,dTangent,ipvcur->getPtrTodCmdFm());
+
+ ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+
+}
+
+double IMDEACPUMATDG3DMaterialLaw::scaleFactor() const{return _vevplaw->scaleFactor();};
+
+void IMDEACPUMATDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
+ dG3DMaterialLaw::setMacroSolver(sv);
+ if (_vevplaw != NULL){
+ _vevplaw->setMacroSolver(sv);
+ }
+};
//
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index a0efb8fee0d4c88b6baf66ee1b009f15c8fc8558..20aaacf2d10764416e83b10c9744b6b24dd45424 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -30,6 +30,7 @@
#include "mlawCrystalPlasticity.h"
#include "mlawGursonUMAT.h"
#include "mlawVEVPUMAT.h"
+#include "mlawIMDEACPUMAT.h"
#include "mlawViscoelastic.h"
#include "mlawLinearElecTherMech.h"
#include "mlawAnIsotropicElecTherMech.h"
@@ -2604,6 +2605,42 @@ class VEVPUMATDG3DMaterialLaw : public dG3DMaterialLaw{
#endif
};
+class IMDEACPUMATDG3DMaterialLaw : public dG3DMaterialLaw{
+ protected:
+ #ifndef SWIG
+ mlawIMDEACPUMAT *_vevplaw;
+ #endif //SWIG
+
+ public:
+ IMDEACPUMATDG3DMaterialLaw(const int num, const char *propName);
+ #ifndef SWIG
+ IMDEACPUMATDG3DMaterialLaw(const IMDEACPUMATDG3DMaterialLaw &source);
+ virtual ~IMDEACPUMATDG3DMaterialLaw();
+
+ // set the time of _cplaw
+ virtual void setTime(const double t,const double dtime);
+ virtual materialLaw::matname getType() const;
+ virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
+ virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}
+ virtual double soundSpeed() const;// or change value ??
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const;
+ virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
+ virtual double scaleFactor() const;
+ virtual materialLaw* clone() const{return new IMDEACPUMATDG3DMaterialLaw(*this);};
+ virtual bool withEnergyDissipation() const;
+ virtual void setMacroSolver(const nonLinearMechSolver* sv);
+ virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
+ {
+ return _vevplaw;
+ }
+ virtual materialLaw *getNonLinearSolverMaterialLaw()
+ {
+ return _vevplaw;
+ }
+ #endif
+};
+
class gursonUMATDG3DMaterialLaw : public dG3DMaterialLaw{
protected:
#ifndef SWIG