- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
I'm finishing my PhD in Geotechincal Engineering, in the Polytechnic University of Catalonia.
I've implemented Fortran Code within Plaxis environtment (User Soil Defined Model). To create the necessaries dynamic libraries .dll, I've been using Intel Fortran Compiler (32/64) and LF95.
Intel Fortran (32/64) compiles the principal part of the code, and Lahey compiles an additional part of the code that let me use Plaxis Interface introducing input parameters.
I've using next make orders (which I've adapted from the original ones) for compiling the code and create the .dll:
IFORT32 (ifcomp.bat)
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
ifort @comp_rel UsrMod.f90 UsrAddDF.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf
IFORT64 (ifcomp64.bat)
REM compile 64-bit DLL: make sure you're in a 64-bit development environment
copy impivf impexp
ifort /static UsrMod.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf64
LF95 (lfcomp.bat)
echo DLL_Export User_Mod
lf95 -c usrlib.f90
lib /out:lfusrlib.lib usrlib.ob
lf95 -c -ml winapi okmsglah.f90
lf95 -win -dll UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll
I've been able to compile the code of the constitutive code that I've developed in PLAXIS, introduce the model parameters, calculate and analyze the results. But now, LF95 trial license has expired (>30 days), and Lahey doesn't offer me any Student Version at a reasonable price.
That's why I'd like to ask you how to change LF95 make orders to Intel Fortran. I've tried all it was in my hand, but I'm not an expert developer (can manage well numerical models, but get lost with compiler orders -except basic ones-).
Thank you very much.
Best regards,
Link Copied
11 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The requirements for the Fortran DLL are on page 174 of the Plaxis documentation http://www.plaxis.nl/files/files/2D-3-Material-Models.pdf .
You need to add the DEC$ directive given on the same page to your Fortran source code, and use the compiler command listed on the same page.
Mentioning the Lahey compiler in this context is an unnecessary distraction.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thak you for your answer.
The structure of the compiler orders that I've written in the previous message is the same as the one you could find in Plaxis USRMOD example in their web (they use both compilers):"example_ivf.dll" and "example_ivf64.dll" are created with Intel Compiler and "example_lah.dll" with Lahey Compiler. "example_ivf.dll" and "example_ivf64.dll" are correct.
I've tried many times to change the directives of lfcomp.bat (the one that creates "example_lah.dll") to Intel directives, but I failed (Plaxis didn't work). That's why I asked.
echo DLL_Export User_Mod ==> echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
lf95 -c usrlib.f90 ==> ¿?¿? (userlib.f90 have no DEC of DLL_Export directive)
lib /out:lfusrlib.lib usrlib.obj ==> ¿?¿?
lf95 -c -ml winapi okmsglah.f90 ==> ¿?¿? (there is an equivalent code with DEC directives, okmsgivf.f90)
lf95 -win -dll UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll ==> If I write any of the following orders, it doesn'st work:
ifort /winapp UsrMod.f90 dfusrlib.lib /dll (not work)
ifort /winapp UsrMod.f90 usrlib.lib /dll (not work)
ifort /winapp UsrMod.f90 usrlib.90 /dll (not work)
ifort /winapp UsrMod.f90 UsrAddDF.f90 usrlib.lib /dll (not work)
ifort /winapp UsrMod.f90 UsrAddDF.f90 okmsgivf.f90 usrlib.lib /dll (not work)
ifort /winapp UsrMod.f90 UsrAddDF.f90 okmsgivf.f90 usrlib.f90 /dll (not work)
ifort /winapp UsrMod.f90 UsrAddDF.f90 usrlib.f90 /dll /exe:example_lah (not work)
ifort @comp_rel UsrMod.f90 UsrAddDF.f90 usrlib.f90 /dll /exe:example_lah (not work)
Comp_rel include: /real_size:64
/fpconstant
/optimize:4
/static
/assume:noaccuracy_sensitive
/assume:byterecl
/assume:buffered_io
/notraceback
/debug:none
/check:none
/heap-arrays:1
/fpp
/recursive
ifort @comp_rel UsrMod.f90 UTsrAddDF.f90 okmsgivf.f90 usrlib.lib /dll /exe:example_lah (not work)
Etc...
Thank You
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I don't understand what the problem is.
Apparently, you have Intel Fortran and you are able to produce with it DLLs that Plaxis can use correctly? And, you used to have Lahey F95 and you were able to use it (until the trial license expired) to produce DLLs which also worked correctly with Plaxis? If so, what is holding you back from creating IFort-compatible versions of the sources that you used with LF95?
What do you mean by "not work" and "it doesn't work"? Failure to compile, failure to link, rejection of DLL by Plaxis? In any case, you should show us the details of the attempted commands and the error messages that you received.
Note that you have not shown any source code. Plaxis is a big commercial package and I am not going to be able to do anything that requires having Plaxis or Delphi. If you post the files UsrMod.f90, UTsrAddDF.f90, okmsgivf.f90 and usrlib.f90, we can help with building an IFort DLL with them.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The problem is that I don't know how to "translate" the make file (and some lines in the code) from LF95 to Intel.
I'm sure that it's possible to work only with Intel, but the files that Plaxis provides as an example of User Defined Soil Model (http://kb.plaxis.nl/models/plaxis-user-defined-soil-models) uses both compilers, and, as I'm not a professional developer, I've followed this structure.
Example Plaxis file uses:
- INTEL to create two .dll: example_ivf.dll (compiled with IA-32) and example_ivf64.dll (compiled with IA-64)
- LF95: to create one .dll: example_lah.dll
I changed ".for" => ".f90", and adapted makefiles. With that, I implemented a very simple model, and after finishing compilation, I put this 3 files in Plaxis folder and all worked well.
In order to use LF95 or Intel Compiler, it's necessary to change one line of the code of USRMOD.f90 (see below).
After that, I started the implementation of an elastoplastic model which already exist in Plaxis (hardening soil). I'm copying below the code of this one.
When I finish its validation (making comparations with Plaxis Hardening Soil Model), I will implement my model code (which tries to improve the elastic part of that one).
I don't have any problem to copy here the code, but it has more than 15.000 lines, so I'm going to copy the main files and some subroutines (if you need more information, please, tell me).
Thank you very much.
The code I'm attaching works perfectly once compiled with INTEL and LF95 (example_ivf.dll with IA-32, example_ivf64.dll with IA-64, and example_lah.dll with LF95).
It uses the libraries provided by Plaxis, and Plaxis structure for subrotuines. As well, it has some errors from a physical point of view which I've already correct but not implemented yet (in any case, they're not important for our propose here).
I compiled this code with next make files (which I executed from the respectives compiler consoles):
IFCOMP.BAT:
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
ifort @comp_rel UsrMod.f90 UsrAddDF.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf
COMP_REL:
/real_size:64
/fpconstant
/optimize:4
/static
/assume:noaccuracy_sensitive
/assume:byterecl
/assume:buffered_io
/notraceback
/debug:none
/check:none
/heap-arrays:1
/fpp
/recursive
IFCOMP64.BAT:
REM compile 64-bit DLL: make sure you're in a 64-bit development environment
copy impivf impexp
ifort /static UsrMod.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf64
IMPIVF:
: User_Mod
LFCOMP.BAT:
echo DLL_Export User_Mod
lf95 -c usrlib.f90
lib /out:lfusrlib.lib usrlib.obj
lf95 -c -ml winapi okmsglah.f90
lf95 -win -dll UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll
CODE:
USRMOD.F90 (TO USE WITH INTEL OR LF95)
_____________________________________________________________________________________________________
Subroutine User_Mod ( IDTask, iMod, IsUndr, &
iStep, iTer, iEl, Int, &
X, Y, Z, &
Time0, dTime, &
Props, Sig0, Swp0, StVar0, &
dEps, D, BulkW, &
Sig, Swp, StVar, ipl, &
nStat, NonSym, iStrsDep, iTimeDep,iTang, &
iPrjDir, iPrjLen, iAbort )
!
Implicit Double Precision (A-H, O-Z)
!
Dimension Props(*), Sig0(*), StVar0(*), dEps(*), D(6,6), &
Sig(*), StVar(*), iPrjDir(*)
Character*255 PrjDir, Dbg_Name
Data iounit / 0 /
Save iounit
! La siguiente linea debe utilizarse para la compilación del DLL
!NEXT TWO LINES MUST BE ACTIVATE OR DESACTIVATE DEPENDING ON THE COMPILER USED:
!DEC$ ATTRIBUTES DLLExport, StdCall, Reference :: User_Mod ! activate with Intel IA-32/64
!DLL_Export User_Mod !activate with LF95
! Posibilidad de debugging
If (iounit.Eq.0) Then
PrjDir=' '
Do i=1,iPrjLen
PrjDir(i:i) = Char(iPrjDir(i))
End Do
Dbg_Name=PrjDir(:iPrjLen)//'udsmex'
nErr=0
1 Continue
Open(Unit= 1, File= Dbg_Name, iostat=ios)
If (ios.Eq.0) Close(Unit=1, Status='delete', iostat=ios)
If (ios.Ne.0) Then
!
! in case of error try ...udsmex1 or udsmex2 or ..
!
nErr=nErr+1
dbg_name=PrjDir(:iPrjLen)//'udsmex'//char(48+nErr)
If (nErr.Lt.10) Goto 1
End If
Open(Unit= 1, File=Dbg_Name, blocksize=4096)
Write(1,*)'File 1 opened'
iounit = 1
! Call WriVec(-1,'Props',Props,14)
End If
! Call WriIvl(-1, 'iounit',iounit )
! Call WriIvl(-1, 'IDTask',IDTask )
Select Case (iMod)
Case (1)
Call MyModel1( IDTask, iMod, IsUndr, iStep, iTer, iEl, Int, &
X, Y, Z, Time0, dTime, &
Props, Sig0, Swp0, StVar0, &
dEps, D, BulkW, Sig, Swp, StVar, ipl, &
nStat, NonSym, iStrsDep, iTimeDep, iTang, &
iAbort )
! Case (2)
! Call MyModel3( IDTask, ....
Case Default
Write(1,*) 'invalid model number in UsrMod', iMod
Write(1,*) 'IDTask: ',IDTask
Stop 'invalid model number in UsrMod'
iAbort=1
Return
End Select ! iMod
If (IDTask .EQ. 5 .AND. iel+int .EQ. 2) Then
Write(1,*)'nStat : ',nStat
Write(1,*)'NonSym : ',NonSym
Write(1,*)'StrsDep : ',iStrsDep
Write(1,*)'TimeDep : ',iTimeDep
Write(1,*)'Tangent : ',iTang
End If
If (IDTask .EQ. 3 .AND. iel+int .EQ. 2) Then
Write(1,*)'IDTask: ', IDTask,'iStep: ', iStep, 'iTer :', iTer
Call Flush(1)
End If
! Call WriIvl(-1, 'IDTask end',IDTask )
Return
End ! User_Mod
include 'Matriz6.f90'
include 'Matriz3.f90'
include 'mymodel1.f90'
include 'HS_plast.f90'
include 'MC_plast.f90'
include 'TRAC_plast.f90'
include 'CAP_plast.f90'
include 'HS_CAP_plast.f90'
include 'MC_CAP_plast.f90'
include 'HS_TRAC_plast.f90'
include 'MC_TRAC_plast.f90'
_____________________________________________________________________________________________________
USR_ADD.F90 (TO USE WITH LF95)
_____________________________________________________________________________________________________
! Subroutines in this file:
! Exported:
! Subroutine GetModelCount( nMod )
! Subroutine GetModelName ( iMod , ModelName )
! Subroutine GetParamCount( iMod , nParam )
! Subroutine GetParamName ( iMod , iParam, ParamName )
! Subroutine GetParamUnit ( iMod , iParam, Units )
! Subroutine GetStateVarCount( iMod , nVar )
! Subroutine GetStateVarName ( iMod , iVar, Name )
! Subroutine GetStateVarUnit ( iMod , iVar, Unit )
!
! Local:
! Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
! Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
!
Subroutine GetModelCount(nMod)
!
! Return the maximum model number (iMod) in this DLL
!
Integer (Kind=4) nMod
Dll_Export GetModelCount
nMod = 1 ! Maximum model number (iMod) in current DLL
Return
End ! GetModelCount
Subroutine GetModelName( iMod , ModelName )
!
! Return the name of the different models
!
Integer iMod
Character (Len= * ) ModelName
Character (Len=255) tName
DLL_Export GetModelName
Select Case (iMod)
Case (1)
tName = 'HS'
Case Default
tName = 'not in DLL'
End Select
LT = Len_Trim(tName)
ModelName= Char(lt) // tName(1:Lt)
Return
End ! GetModelName
Subroutine GetParamCount( iMod , nParam )
!
! Return the number of parameters of the different models
!
DLL_Export GetParamCount
Select Case (iMod)
Case ( 1 )
nParam = 14
Case Default
nParam = 0
End Select
Return
End ! GetParamCount
Subroutine GetParamName( iMod , iParam, ParamName )
!
! Return the parameters name of the different models
!
Character (Len=255) ParamName, Units
DLL_Export GetParamName
Call GetParamAndUnit(iMod,iParam,ParamName,Units)
Return
End
Subroutine GetParamUnit( iMod , iParam, Units )
!
! Return the units of the different parameters of the different models
!
Character (Len=255) ParamName, Units
DLL_Export GetParamUnit
Call GetParamAndUnit(iMod,iParam,ParamName,Units)
Return
End
Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
!
! Return the parameters name and units of the different models
!
! Units: use F for force unit
! L for length unit
! T for time unit
!
Character (Len=255) ParamName, Units, tName
Select Case (iMod)
Case (1)
! ModName = 'EL'
Select Case (iParam)
Case (1)
ParamName = 'G' ; Units = 'F/L^2#'
Case (2)
ParamName = 'nu' ; Units = '-'
Case (3)
ParamName = 'E50_ref' ; Units = 'F/L^2#'
Case (4)
ParamName = 'Eoed_ref' ; Units = 'F/L^2#'
Case (5)
ParamName = 'Eur_ref' ; Units = 'F/L^2#'
Case (6)
ParamName = 'm' ; Units = '-'
Case (7)
ParamName = 'c' ; Units = 'F/L^2#'
Case (8)
ParamName = 'phi' ; Units = 'deg'
Case (9)
ParamName = 'psi' ; Units = 'deg'
Case (10)
ParamName = 'nu_ur' ; Units = '-'
Case (11)
ParamName = 'p_ref' ; Units = 'F/L^2#'
Case (12)
ParamName = 'K0_nc' ; Units = '-'
Case (13)
ParamName = 'R_f' ; Units = '-'
Case (14)
ParamName = 'sitma_tracc'; Units = 'F/L^2#'
Case Default
ParamName = '???' ; Units = '???'
End Select
Case Default
! model not in DLL
ParamName = ' N/A ' ; Units = ' N/A '
End Select
tName = ParamName
LT = Len_Trim(tName)
ParamName= Char(lt) // tName(1:Lt)
tName = Units
LT = Len_Trim(tName)
Units = Char(lt) // tName(1:Lt)
Return
End ! GetParamAndUnit
Subroutine GetStateVarCount( iMod , nVar )
!
! Return the number of state variables of the different models
!
Dll_export GetStateVarCount
Select Case (iMod)
Case (1)
nVar = 2
Case Default
nVar = 0
End Select
Return
End
Subroutine GetStateVarName( iMod , iVar, Name )
!
! Return the name of the different state variables
! of the different models
!
Character (Len=255) Name, Unit
Dll_export GetStateVarName
Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
End
Subroutine GetStateVarUnit( iMod , iVar, Unit )
!
! Return the units of the different state variables of the different models
!
Character (Len=255) Name, Unit
Dll_export GetStateVarUnit
Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
End
Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
!
! Return the name and unit of the different state variables of the different models
!
Character (Len=255) Name, Unit, tName
Select Case (iMod)
Case (1)
Select Case (iVar)
Case (1)
Name = 'gamma_p' ; Unit = '-'
Case (2)
Name = 'p_p' ; Unit = 'F/L^2#'
Case Default
Name='N/A' ; Unit = '?'
End Select
Case Default
Name='N/A' ; Unit = '?'
End Select
tName = Name
LT = Len_Trim(tName)
Name = Char(lt) // tName(1:Lt)
tName = Unit
LT = Len_Trim(tName)
Unit = Char(lt) // tName(1:Lt)
Return
End
_____________________________________________________________________________________________________
USRADDDF.F90 (TO USE WITH INTEL)
_____________________________________________________________________________________________________
! Subroutines in this file:
! Exported:
! Subroutine GetModelCount( nMod )
! Subroutine GetModelName ( iMod , ModelName )
! Subroutine GetParamCount( iMod , nParam )
! Subroutine GetParamName ( iMod , iParam, ParamName )
! Subroutine GetParamUnit ( iMod , iParam, Units )
! Subroutine GetStateVarCount( iMod , nVar )
! Subroutine GetStateVarName ( iMod , iVar, Name )
! Subroutine GetStateVarUnit ( iMod , iVar, Unit )
!
! Local:
! Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
! Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
Subroutine GetModelCount(nMod)
!
! Return the maximum model number (iMod) in this DLL
!
Integer (Kind=4) nMod
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetModelCount
nMod = 1 ! Maximum model number (iMod) in current DLL
Return
End ! GetModelCount
Subroutine GetModelName( iMod , ModelName )
!
! Return the name of the different models
!
Integer iMod
Character (Len= * ) ModelName
Character (Len=255) tName
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetModelName
Select Case (iMod)
Case (1)
tName = 'HS'
Case Default
tName = 'not in DLL'
End Select
LT = Len_Trim(tName)
ModelName= Char(lt) // tName(1:Lt)
Return
End ! GetModelName
Subroutine GetParamCount( iMod , nParam )
!
! Return the number of parameters of the different models
!
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetParamCount
Select Case (iMod)
Case ( 1 )
nParam = 14
Case Default
nParam = 0
End Select
Return
End ! GetParamCount
Subroutine GetParamName( iMod , iParam, ParamName )
!
! Return the parameters name of the different models
!
Character (Len=255) ParamName, Units
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetParamName
Call GetParamAndUnit(iMod,iParam,ParamName,Units)
Return
End
Subroutine GetParamUnit( iMod , iParam, Units )
!
! Return the units of the different parameters of the different models
!
Character (Len=255) ParamName, Units
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetParamUnit
Call GetParamAndUnit(iMod,iParam,ParamName,Units)
Return
End
Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
!
! Return the parameters name and units of the different models
!
! Units: use F for force unit
! L for length unit
! T for time unit
!
Character (Len=255) ParamName, Units, tName
Select Case (iMod)
Case (1)
! ModName = 'HS'
Select Case (iParam)
Case (1)
ParamName = 'G' ; Units = 'F/L^2#'
Case (2)
ParamName = 'nu' ; Units = '-'
Case (3)
ParamName = 'E50_ref' ; Units = 'F/L^2#'
Case (4)
ParamName = 'Eoed_ref' ; Units = 'F/L^2#'
Case (5)
ParamName = 'Eur_ref' ; Units = 'F/L^2#'
Case (6)
ParamName = 'm' ; Units = '-'
Case (7)
ParamName = 'c' ; Units = 'F/L^2#'
Case (8)
ParamName = 'phi' ; Units = 'deg'
Case (9)
ParamName = 'psi' ; Units = 'deg'
Case (10)
ParamName = 'nu_ur' ; Units = '-'
Case (11)
ParamName = 'p_ref' ; Units = 'F/L^2#'
Case (12)
ParamName = 'K0_nc' ; Units = '-'
Case (13)
ParamName = 'R_f' ; Units = '-'
Case (14)
ParamName = 'sitma_tracc'; Units = 'F/L^2#'
Case Default
ParamName = '???' ; Units = '???'
End Select
Case Default
! model not in DLL
ParamName = ' N/A ' ; Units = ' N/A '
End Select
tName = ParamName
LT = Len_Trim(tName)
ParamName= Char(lt) // tName(1:Lt)
tName = Units
LT = Len_Trim(tName)
Units = Char(lt) // tName(1:Lt)
Return
End ! GetParamAndUnit
Subroutine GetStateVarCount( iMod , nVar )
!
! Return the number of state variables of the different models
!
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarCount
Select Case (iMod)
Case (1)
nVar = 2
Case Default
nVar = 0
End Select
Return
End
Subroutine GetStateVarName( iMod , iVar, Name )
!
! Return the name of the different state variables of the different models
!
Character (Len=255) Name, Unit
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarName
Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
End
Subroutine GetStateVarUnit( iMod , iVar, Unit )
!
! Return the units of the different state variables of the different models
!
Character (Len=255) Name, Unit
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarUnit
Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
End
Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
!
! Return the name and unit of the different state variables of the different models
!
Character (Len=255) Name, Unit, tName
Select Case (iMod)
Case (1)
Select Case (iVar)
Case (1)
Name = 'gamma_p' ; Unit = '-'
Case (2)
Name = 'p_p' ; Unit = 'F/L^2#'
Case Default
Name='N/A' ; Unit = '?'
End Select
Case Default
Name='N/A' ; Unit = '?'
End Select
tName = Name
LT = Len_Trim(tName)
Name = Char(lt) // tName(1:Lt)
tName = Unit
LT = Len_Trim(tName)
Unit = Char(lt) // tName(1:Lt)
Return
End
_____________________________________________________________________________________________________
OKMSGIVF.F90 (TO USE WITH INTEL)
_____________________________________________________________________________________________________
Subroutine OK_MessageBox(t)
use dfwin
use dfwinty
integer (kind=Int_Ptr_Kind()) hWnd
character*(*) t
character(len=256) :: msg, title
! Display a messagebox with an OK button
! Note that all strings must be null terminated for C's sake
msg = Trim(t) // char(0)
title = 'UDSM_Example' // char(0)
hWnd = 0
iret = MessageBox( hWnd, &
msg, &
title, &
MB_OK )
Return
End
_____________________________________________________________________________________________________
OKMSGLAH.F90 (TO USE WITH LF95)
_____________________________________________________________________________________________________
Subroutine OK_MessageBox(t)
dll_import MessageBoxA
integer :: MessageBoxA
Integer is3DF
! Some (not all) Message box flags. The C header file winuser.h
! contains a complete list of flags. Winuser.h is not provided
! in the LF95 distribution.
integer, parameter :: MB_OK = 0
integer :: iret
character(len=*) :: t
character(len=256) :: msg, title
! Display a messagebox with an OK button
! Note that all strings must be null terminated for C's sake
msg = Trim(t) // char(0)
title = 'UDSM_example' // char(0)
iret = MessageBoxA(val(0),&
val(pointer(msg)), &
val(pointer(title)), &
val(MB_OK))
Return
End
_____________________________________________________________________________________________________
USRLIB.F90 (SUBROUTINES; LIBRARY)
_____________________________________________________________________________________________________
Subroutine MZEROR(R,K)
!C
!C***********************************************************************
!C
!C Function: To make a real array R with dimension K to zero
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension R(*)
Do J=1,K
R(J) = 0.0D0
End Do
Return
End
Subroutine MZEROI(I,K)
!C
!C***********************************************************************
!C
!C Function: To make an integre array I with Dimension K to zero
!C
!C***********************************************************************
!C
Dimension I(*)
Do J=1,K
I(J)=0
End Do
Return
End
Subroutine SETRVAL(R,K,V)
!C
!C***********************************************************************
!C
!C Function: To fill a real array R with Dimension K with value V
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension R(*)
Do J=1,K
R(J)=V
End Do
Return
End
Subroutine SETIVAL(I,K,IV)
!C
!C***********************************************************************
!C
!C Function: To fill an integer array I with Dimension K with value IV
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension I(*)
Do J=1,K
I(J)=IV
End Do
Return
End
Subroutine COPYIVEC(I1,I2,K)
!C
!C***********************************************************************
!C
!C Function: To copy an integer array I1 with Dimension K to I2
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension I1(*),I2(*)
Do J=1,K
I2(J)=I1(J)
End Do
Return
End
Subroutine COPYRVEC(R1,R2,K)
!C
!C***********************************************************************
!C
!C Function: To copy a Double array R1 with Dimension K to R2
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension R1(*),R2(*)
Do J=1,K
R2(J)=R1(J)
End Do
Return
End
Logical Function IS0ARR(A,N)
!C
!C***********************************************************************
!C Function : To check whether a real array contains only zero values.
!C When an array contains only zero's is might not need to be
!C written to the XXX file.
!C exit Function when first non-zero value occured or when
!C all elements are checked and are zero.
!C
!C Input: A : array to be checked
!C N : number of elements in array that should be checked
!C
!!C Output : .TRUE. when all elements are 0
!C .FALSE. when at least one element is not zero
!C
!C Called by : Subroutine TOBXX
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension A(*)
Is0Arr=.False.
Do I=1,N
If ( A(I) .Ne. 0 ) Return
End Do
Is0Arr=.True.
Return
End
Logical Function IS0IARR(IARR,N)
!C
!C***********************************************************************
!C Function : To check whether a integer array contains only zero values.
!C Similar to IS0ARR
!C
!C Input: IARR : array to be checked
!C N : number of elements in array that should be checked
!C
!C Output : .TRUE. when all elements are 0
!C .FALSE. when at least one element is not zero
!C
!C Called by : Subroutine TOBXX
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Dimension IARR(*)
Is0IArr=.False.
Do I=1,N
If ( IARR(I) .Ne. 0 ) Return
End Do
Is0IArr=.True.
Return
End
!C***********************************************************************
Subroutine MulVec(V,F,K)
!C***********************************************************************
!C
!C Function: To multiply a real vector V with dimension K by F
!!C
!C***********************************************************************
!C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION V(*)
Do J=1,K
V(J)=F*V(J)
End Do
Return
End ! Subroutine Mulvec
!C***********************************************************************
Subroutine MatVec(xMat,IM,Vec,N,VecR)
!!C***********************************************************************
!C
!C Calculate VecR = xMat*Vec
!!C
!C I xMat : (Square) Matrix (IM,*)
!C I Vec : Vector
!C I N : Number of rows/colums
!C O VecR : Resulting vector
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension xMat(IM,*),Vec(*),VecR(*)
!C***********************************************************************
Do I=1,N
X=0
Do J=1,N
X=X+xMat(I,J)*Vec(J)
End Do
VecR(I)=X
End Do
Return
End ! Subroutine MatVec
!C***********************************************************************
Subroutine AddVec(Vec1,Vec2,R1,R2,N,VecR)
!C***********************************************************************
!C
!C Calculate VecR() = R1*Vec1()+R2*Vec2()
!C
!C I Vec1,
!C I Vec2 : Vectors
!C I R1,R2 : Multipliers
!C I N : Number of rows
!C O VecR : Resulting vector
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension Vec1(*),Vec2(*),VecR(*)
!C***********************************************************************
Do I=1,N
X=R1*Vec1(I)+R2*Vec2(I)
VecR(I)=X
End Do
Return
End ! Subroutine AddVec
!C
!C***********************************************************************
Subroutine DInProd(A,B,N,RES)
!C***********************************************************************
!C
!C Returns the Inproduct of two vectors
!C
!C I A,B : Two vectors
!C I N : Used length of vectors
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension A(*),B(*)
!C***********************************************************************
X = 0
Do I=1,N
X = X + A(I)*B(I)
End Do
RES = X
Return
End ! Function DInProd
!C
!C***********************************************************************
Subroutine MatMat(xMat1,Id1,xMat2,Id2,nR1,nC2,nC1,xMatR,IdR)
!C***********************************************************************
!!C
!C Calculate xMatR = xMat1*xMat2
!C
!C I xMat1 : Matrix (Id1,*)
!C I xMat2 : Matrix (Id2,*)
!C I nR1 : Number of rows in resulting matrix (= No rows in xMat1)
!C I nC2 : Number of columns in resulting matrix (= No cols in xMat2)
!C I nC1 : Number of columns in matrix xMat1
!C = Number rows in matrix xMat2
!C O xMatR : Resulting matrix (IdR,*)
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension xMat1(Id1,*),xMat2(Id2,*),xMatR(IdR,*)
!C**********************************************************************
Do I=1,nR1
Do J=1,nC2
X=0
Do K=1,nC1
X=X+xMat1(I,K)*xMat2(K,J)
End Do
xMatR(I,J)=X
End Do
End Do
Return
End ! Subroutine MatMat
!C***********************************************************************
Subroutine MatMatSq(n, xMat1, xMat2, xMatR)
!C***********************************************************************
!C
!C Calculate xMatR = xMat1*xMat2 for square matrices, size n
!C
!C I n : Dimension of matrices
!C I xMat1 : Matrix (n,*)
!C I xMat2 : Matrix (n,*)
!C O xMatR : Resulting matrix (n,*)
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension xMat1(n,*),xMat2(n,*),xMatR(n,*)
!C**********************************************************************
Do I=1,n
Do J=1,n
X=0
Do K=1,n
X=X+xMat1(I,K)*xMat2(K,J)
End Do
xMatR(I,J)=X
End Do
End Do
Return
End ! Subroutine MatMatSq
!C***********************************************************************
Subroutine WriVal ( io, C , V )
!C***********************************************************************
!C
!!C Write (Double) value to file unit io (when io>0)
!C
!C***********************************************************************
!C
Implicit Double Precision (A-H,O-Z)
Character C*(*)
If (io.Le.0) Return
Write(io,*) C,V
1 Format( A,3x, 1x,1p,e12.5)
Return
End
!C***********************************************************************
Subroutine WriIVl ( io, C , I )
!C***********************************************************************
!C
!C Write (integer) value to file unit io (when io>0)
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Character C*(*)
If (io.Le.0) Return
Write(io,*) C,I
1 Format( A,3x, 1x,I6)
Return
End
!C***********************************************************************
Subroutine WriIVc ( io, C , iV , n )
!C***********************************************************************
!C
!C Write (integer) vector to file unit io (when io>0)
!C
!C***********************************************************************
Character C*(*)
Dimension iV(*)
If (io.Le.0) Return
Write(io,*) C
Write(io,1) (iv(i),i=1,n)
1 Format( ( 2(3x,5i4) ) )
Return
End
!C***********************************************************************
Subroutine WriVec ( io, C , V , n )
!C***********************************************************************
!C
!C Write (Double) vector to file unit io (when io>0)
!C 6 values per line
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Character C*(*)
Dimension V(*)
If (io.Le.0) Return
If (Len_Trim(C).Le.6) Then
Write(io,2) C,( V(i),i=1,n)
Else
Write(io,*) C
Write(io,1) ( V(i),i=1,n)
End If
1 Format( ( 2(1x, 3(1x,1p,e10.3) ) ) )
2 Format( A, ( T7, 2(1x, 3(1x,1p,e10.3) ) ) )
Return
End
!C***********************************************************************
Subroutine WriVec5( io, C , V , n )
!C***********************************************************************
!C
!C Write (Double) vector to file unit io (when io>0)
!C 5 values per line
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Character C*(*)
Dimension V(*)
If (io.Le.0) Return
Write(io,*) C
Write(io,1) ( V(i),i=1,n)
1 Format( 5(1x,1p,e12.5) )
Return
End
!C***********************************************************************
Subroutine WriMat ( io, C , V , nd, nr, nc )
!C***********************************************************************
!C
!C Write (Double) matrix to file unit io (when io>0)
!C 6 values per line
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Character C*(*)
Dimension V(nd,*)
If (io.Le.0) Return
Write(io,*) C
Do j=1,nr
Write(io,1) j,( V(j,i),i=1,nc)
End Do
1 Format(i4, ( T7,2(1x, 3(1x,1p,e10.3) ) ) )
Return
End
!C***********************************************************************
Subroutine MatInvPiv(Aorig,B,N)
Implicit Double Precision (A-H,O-Z)
Dimension Aorig(n,*), B(n,*),A(:,:)
Allocatable :: A
Allocate ( A(n,n) ) ! No error checking !!
Call CopyRVec(AOrig, A, n*n )
Call MZeroR(B,n*n)
Do i=1,n
B(i,i) = 1d0
End Do
Do I=1,n
T=A(I,I)
iPiv=i
Do j=i+1,n
If ( Abs(A(j,i)) .Gt. Abs(A(iPiv,i)) ) iPiv=j
End Do
If (iPiv.Ne.i) Then
Do j=1,n
x = A( i ,j)
A( i ,j) = A(iPiv,j)
A(iPiv,j) = x
x = B( i ,j)
B( i ,j) = B(iPiv,j)
B(iPiv,j) = x
End Do
T=A(I,I)
End If
Do J=1,n
A(I,J)=A(I,J)/T
B(I,J)=B(I,J)/T
End Do
Do K=1,n
If (K.Ne.I) Then
T=A(K,I)
Do J=1,n
A(K,J)=A(K,J)-T*A(I,J)
B(K,J)=B(K,J)-T*B(I,J)
End Do
End If
End Do
End Do
DeAllocate ( A )
Return
End ! MatinvPiv
!C***********************************************************************
Subroutine PrnSig(IOpt,S,xN1,xN2,xN3,S1,S2,S3,P,Q)
Implicit Double Precision (A-H,O-Z)
Dimension S(*),xN1(*),xN2(*),xN3(*)
If (iOpt.Eq.1) Then
Call Eig_3(0,S,xN1,xN2,xN3,S1,S2,S3,P,Q) ! with Eigenvectors
Else
Call Eig_3a(0,S,S1,S2,S3,P,Q) ! no Eigenvectors
End If
Return
End
!C***********************************************************************
Subroutine Eig_3(iOpt,St,xN1,xN2,xN3,S1,S2,S3,P,Q)
Implicit Double Precision (A-H,O-Z)
Dimension St(6),A(3,3),V(3,3), &
xN1(3),xN2(3),xN3(3)
!
! Get Eigenvalues/Eigenvectors for 3*3 matrix
! Wim Bomhof 15/11/'01
! PGB : adaption to Principal stress calculation
!
! Applied on principal stresses, directions
! Stress vector St(): XX, YY, ZZ, XY, YZ, ZX
!
A(1,1) = St(1) ! xx
A(1,2) = St(4) ! xy = yx
A(1,3) = St(6) ! zx = xz
A(2,1) = St(4) ! xy = yx
A(2,2) = St(2) ! yy
A(2,3) = St(5) ! zy = yz
A(3,1) = St(6) ! zx = xz
A(3,2) = St(5) ! zy = yz
A(3,3) = St(3) ! zz
! Set V to unity matrix
V(1,1) = 1
V(2,1) = 0
V(3,1) = 0
V(1,2) = 0
V(2,2) = 1
V(3,2) = 0
V(1,3) = 0
V(2,3) = 0
V(3,3) = 1
abs_max_s=0.0
Do i=1,3
Do j=1,3
if (abs(a(i,j)) .Gt. abs_max_s) abs_max_s=abs(a(i,j))
End Do
End Do
Tol = 1d-20 * abs_max_s
it = 0
itmax = 50
Do While ( it.Lt.itMax .And. &
abs(a(1,2))+abs(a(2,3))+abs(a(1,3)) .Gt. Tol )
it=it+1
Do k=1,3
If (k .Eq. 1) Then
ip=1
iq=2
Else If (k .Eq.2) Then
ip=2
iq=3
Else
ip=1
iq=3
End If
If (a(ip,iq) .Ne. 0.0) Then
tau=(a(iq,iq)-a(ip,ip))/(2.0*a(ip,iq))
If (tau .Ge.0.0) Then
sign_tau=1.0
Else
sign_tau=-1.0
End If
t=sign_tau/(abs(tau)+sqrt(1.0+tau*tau))
c=1.0/sqrt(1.0+t*t)
s=t*c
a1p=c*a(1,ip)-s*a(1,iq)
a2p=c*a(2,ip)-s*a(2,iq)
a3p=c*a(3,ip)-s*a(3,iq)
a(1,iq)=s*a(1,ip)+c*a(1,iq)
a(2,iq)=s*a(2,ip)+c*a(2,iq)
a(3,iq)=s*a(3,ip)+c*a(3,iq)
a(1,ip)=a1p
a(2,ip)=a2p
a(3,ip)=a3p
v1p=c*v(1,ip)-s*v(1,iq)
v2p=c*v(2,ip)-s*v(2,iq)
v3p=c*v(3,ip)-s*v(3,iq)
v(1,iq)=s*v(1,ip)+c*v(1,iq)
v(2,iq)=s*v(2,ip)+c*v(2,iq)
v(3,iq)=s*v(3,ip)+c*v(3,iq)
v(1,ip)=v1p
v(2,ip)=v2p
v(3,ip)=v3p
ap1=c*a(ip,1)-s*a(iq,1)
ap2=c*a(ip,2)-s*a(iq,2)
ap3=c*a(ip,3)-s*a(iq,3)
a(iq,1)=s*a(ip,1)+c*a(iq,1)
a(iq,2)=s*a(ip,2)+c*a(iq,2)
a(iq,3)=s*a(ip,3)+c*a(iq,3)
a(ip,1)=ap1
a(ip,2)=ap2
a(ip,3)=ap3
End If ! a(ip,iq)<>0
End Do ! k
End Do ! While
! principal values on diagonal of a
S1 = a(1,1)
S2 = a(2,2)
S3 = a(3,3)
! Derived invariants
P = (S1+S2+S3)/3
Q = Sqrt( ( (S1-S2)**2 + (S2-S3)**2 + (S3-S1)**2 ) / 2 )
! Sort eigenvalues S1 <= S2 <= S3
is1 = 1
is2 = 2
is3 = 3
if (s1.Gt.s2) Then
t = s2
s2 = s1
s1 = t
it = is2
is2 = is1
is1 = it
End If
if (s2.Gt.s3) Then
t = s3
s3 = s2
s2 = t
it = is3
is3 = is2
is2 = it
End If
if (s1.Gt.s2) Then
t = s2
s2 = s1
s1 = t
it = is2
is2 = is1
is1 = it
End If
Do i=1,3
xN1(i) = v(i,is1) ! first column
xN2(i) = v(i,is2) ! second column
xN3(i) = v(i,is3) ! third column
End Do
Return
End ! Eig_3
Subroutine Eig_3a(iOpt,St,S1,S2,S3,P,Q) ! xN1,xN2,xN3,
Implicit Double Precision (A-H,O-Z)
Dimension St(6),A(3,3) ! V(3,3),xN1(3),xN2(3),xN3(3)
!
! Get Eigenvalues ( no Eigenvectors) for 3*3 matrix
! Wim Bomhof 15/11/'01
!
! Applied on principal stresses, directions
! Stress vector XX, YY, ZZ, XY, YZ, ZX
!
A(1,1) = St(1) ! xx
A(1,2) = St(4) ! xy = yx
A(1,3) = St(6) ! zx = xz
A(2,1) = St(4) ! xy = yx
A(2,2) = St(2) ! yy
A(2,3) = St(5) ! zy = yz
A(3,1) = St(6) ! zx = xz
A(3,2) = St(5) ! zy = yz
A(3,3) = St(3) ! zz
abs_max_s=0.0
Do i=1,3
Do j=1,3
if (abs(a(i,j)) .Gt. abs_max_s) abs_max_s=abs(a(i,j))
End Do
End Do
Tol = 1d-20 * abs_max_s
If (iOpt.Eq.1) Tol = 1d-50*abs_max_s
it=0
itmax = 50
Do While ( it.lt.itmax .And.&
abs(a(1,2))+abs(a(2,3))+abs(a(1,3)) .Gt. Tol )
it=it+1
Do k=1,3
If (k .Eq. 1) Then
ip=1
iq=2
Else If (k .Eq.2) Then
ip=2
iq=3
Else
ip=1
iq=3
End If
If (a(ip,iq) .Ne. 0.0) Then ! ongelijk nul ?
tau=(a(iq,iq)-a(ip,ip))/(2.0*a(ip,iq))
If (tau .Ge.0.0) Then
sign_tau=1.0
Else
sign_tau=-1.0
End If
t=sign_tau/(abs(tau)+sqrt(1.0+tau*tau))
c=1.0/sqrt(1.0+t*t)
s=t*c
a1p=c*a(1,ip)-s*a(1,iq)
a2p=c*a(2,ip)-s*a(2,iq)
a3p=c*a(3,ip)-s*a(3,iq)
a(1,iq)=s*a(1,ip)+c*a(1,iq)
a(2,iq)=s*a(2,ip)+c*a(2,iq)
a(3,iq)=s*a(3,ip)+c*a(3,iq)
a(1,ip)=a1p
a(2,ip)=a2p
a(3,ip)=a3p
ap1=c*a(ip,1)-s*a(iq,1)
ap2=c*a(ip,2)-s*a(iq,2)
ap3=c*a(ip,3)-s*a(iq,3)
a(iq,1)=s*a(ip,1)+c*a(iq,1)
a(iq,2)=s*a(ip,2)+c*a(iq,2)
a(iq,3)=s*a(ip,3)+c*a(iq,3)
a(ip,1)=ap1
a(ip,2)=ap2
a(ip,3)=ap3
End If ! a(ip,iq)<>0
End Do ! k
End Do ! While
! principal values on diagonal of a
S1 = a(1,1)
S2 = a(2,2)
S3 = a(3,3)
! Derived invariants
P = (S1+S2+S3)/3
Q = Sqrt( ( (S1-S2)**2 + (S2-S3)**2 + (S3-S1)**2 ) / 2 )
if (s1.Gt.s2) Then
t = s2
s2 = s1
s1 = t
End If
if (s2.Gt.s3) Then
t = s3
s3 = s2
s2 = t
End If
if (s1.Gt.s2) Then
t = s2
s2 = s1
s1 = t
End If
Return
End ! Eig_3a
!C
!C***********************************************************************
Logical Function LEqual(A,B,Eps)
!C***********************************************************************
!C
!C Returns .TRUE. when two real values are (almost) equal,
!C .FALSE. otherwise
!!C
!C I A,B : Two real values to be compared
!C I Eps : Toleration (Magnitude ~= 1E-5)
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
!C***********************************************************************
LEqual =.True.
If (A .Eq. B) Return
If (DAbs(A-B) .LT. 0.5D0*Eps*( DAbs(A) + DAbs(B) + Eps ) )Return
LEqual =.False.
Return
End ! function LEqual
!C
!C***********************************************************************
Subroutine CrossProd(xN1,xN2,xN3)
!C***********************************************************************
!!C
!C Returns cross product of xN1 and xN2
!C
!C I xN1,xN2 : Two basic vectors
!C O xN3 : Resulting vector
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension xN1(*),xN2(*),xN3(*)
!C***********************************************************************
xN3(1) = xN1(2)*xN2(3) - xN1(3)*xN2(2)
xN3(2) = xN1(3)*xN2(1) - xN1(1)*xN2(3)
xN3(3) = xN1(1)*xN2(2) - xN1(2)*xN2(1)
Return
End ! Subroutine CrossProd
!C
!C***********************************************************************
Double Precision Function ArcSin(X,ie)
!C***********************************************************************
!C
!C Returns the Arc Sine of X
!C
!C I X : Input value
!C
!C Note : In stead of using default routine DASIN we use this one
!C because ³X³ can be slightly beyond 1 and this will give
!C a RTE using DASIN(X)
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
!C***********************************************************************
Ie=0
S = (1-X*X)
! If (S .Lt. -1E-10) Ie=1
! If (S .Lt. -1E-10) Write(*,1) X,S
! If (S .Lt. -1E-10) Write(2,1) X,S
1 Format(' ArcSin(',1x,1p,e13.5e3,') , S =',1x,1p,e13.5e3)
If (S.LT.0) S = 0
S = DSQRT(S)
ArcSin = DATan2(X,S)
Return
End ! function ArcSin
!C
!C***********************************************************************
Subroutine CarSig(S1,S2,S3,xN1,xN2,xN3,SNew)
!C***********************************************************************
!C
!C Returns the Cartesian stresses using the principal stresses S1..S3
!C and the principal directions
!C
!C I S1..S3 : Principal stresses
!C I xN1..xN3 : Principal directions (xNi for Si)
!C
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension xN1(*),xN2(*),xN3(*),SNew(*)
Dimension SM(3,3),T(3,3),TT(3,3),STT(3,3)
!C***********************************************************************
!C
!C**** Fill transformation (rotation) matrix
!C
Do I=1,3
T(I,1) = xN1(I)
T(I,2) = xN2(I)
T(I,3) = xN3(I)
TT(1,I) = T(I,1)
TT(2,I) = T(I,2)
TT(3,I) = T(I,3)
End Do
! Call MatTranspose(T,3,TT,3,3,3)
Call MZeroR(SM,9)
SM(1,1) = S1
SM(2,2) = S2
SM(3,3) = S3
!C
!C**** SMnew = T*SM*TT
!C
Call MatMat(SM ,3, TT,3 , 3,3,3 ,STT,3)
Call MatMat( T ,3, STT,3 , 3,3,3 ,SM ,3)
! Call MatMatSq(3, SM, TT, STT ) ! STT = SM*TT
! Call MatMatSq(3, T, STT, SM ) ! SM = T*STT
!C
!C**** Extract cartesian stress vector from stress matrix
!C
Do I=1,3
SNew(I) = SM(I,I)
End Do
SNew(4) = SM(2,1)
SNew(5) = SM(3,2)
SNew(6) = SM(3,1)
Return
End ! Subroutine CarSig
!C**********************************************************************
subroutine setveclen(xn,n,xl)
!C**********************************************************************
Implicit Double Precision (A-H,O-Z)
Dimension xN(*)
x=0
do i=1,n
x=x+xn(i)**2
end do
if (x.Ne.0) Then
f=xl/sqrt(x)
do i=1,3
xn(i)=xn(i)*f
end do
end if
return
end ! setveclen
!C**********************************************************************
!C End Of file
!C**********************************************************************
_____________________________________________________________________________________________________
MYMODEL1.F90 (SUBROUTINE)
_____________________________________________________________________________________________________
Subroutine MyModel1 ( IDTask, iMod, IsUndr, &
iStep, iTer, iEl, Int, &
X, Y, Z, &
Time0, dTime, &
Props, Sig0, Swp0, StVar0, &
dEps, D, BulkW, &
Sig, Swp, StVar, ipl, &
nStat, &
NonSym, iStrsDep, iTimeDep, iTang, &
iAbort )
!
! Número de tarea IDTask 1 : Inicializar variables
! 2 : Calcular tensiones
! 3 : Calcular la matriz de rigidez
! 4 : Número de variables de estado
! 5 : Preguntar propiedades de la matriz
! 6 : Calcular la matriz de rigidez elástica del material
!
! Argumentos:
! I/O Type
! IDTask I I : Ver arriba
! iMod I I : Número de modelo (1..10)
! IsUndr I I : =1 para no drenado, =0 en otro caso
! iStep I I : Número de paso global
! iter I I : Número de iteración global
! iel I I : Número de elemento global
! Int I I : Número de punto de integración global
! X I R : Posición X del punto de integración
! Y I R : Posición Y del punto de integración
! Z I R : Posición Z del punto de integración
! Time0 I R : Tiempo al inicio del paso de cálculo
! dTime I R : Incremento de tiempo
! Props I R() : Lista con los parámetros del modelo
! Sig0 I R() : Tensiones al inicio del paso de cálculo
! Swp0 I R : Exceso d mme presión de poros al inicio del paso de cálculo
! StVar0 I R() : Variable de estado al inicio del paso de cálculo
! dEps I R() : Incremento de deformación
! D I/O R(,) : Matriz de rigidez del material
! BulkW I/O R : Módulo volumétrico del agua (sólo caso no drenado)
! Sig O R() : Tensión resultante
! Swp O R : Exceso de presión de agua resultante
! StVar O R() : Resulting values state variables
! ipl O I : Indicador plástico
! nStat O I : Número de variables de estado
! NonSym O I : Matriz D no simétrica
! iStrsDep O I : =1 para matriz D dependiente de las tesniones
! iTimeDep O I : =1 para matriz D dependiente del tiempo
! iAbort O I : =1 para detenr proceso de cálculo
!
Implicit Double Precision (A-H, O-Z)
Dimension Props(*), Sig0(*), StVar0(*), dEps(*), D(6,6), &
Sig(*), StVar(*)
!--- Variables locales
!Prs_E : Tensiones principales trial (compresión positiva)
!Prs : Tensiones principales al final del paso de cálculo (compresión positiva)
!Si : Tensiones principales trial (compresión negativa)
!xN(i) : Autocvectores de SigE (xx-yy-zz) = Prs_E (1-2-3)
Dimension dSig(6), xN1(3), xN2(3), xN3(3), xN10(3), xN20(3), xN30(3), &
Prs_E(3), Prs(3), Prs_0(3)
pi = 3.1415926535d0
Rad = pi/180d0
io = 0
ttol = 1d-9
nStatV = 2
ideb = 1 !para debugging
!Contenido de Props()
! I/O Type
! 1: G I R Módulo de corte
! 2: nu I R Coeficiente de Poisson
! 3: E50_ref I R Rigidez triaxial secante (superficie cono)
! 4: Eoed_ref I R Rigidez edométrica tangente (superficie cap)
! 5: Eur_ref I R Rigidez en descarga/recarga
! 6: xm I R Potencia de la dependencia tensional
! 7: c I R Cohesión (efectva)
! 8: phi I R Ángulo de fricción (efectiva)
! 9: psi I R Ángulo de dilatancia
! 10: xnu_ur I R Coeficiente de poisson en descarga/recarga
! 11: p_ref I R Presión de refencia para la rigidez
! 12: xK0_nc I R K0 (normalmente consolidado)
! 13: R_f I R Ratio de rotura
! 14: Sigma_trac I R Máxima tensión de tracción (aplica sobre Prs, tracción +, compresión -)
! Prs_E I R(3) : Tensiones principales del predictor elástico S1 >= S2 >= S3
! Prs O R(3) : Tensiones principales resultantes
! ipl O I : Indicador de punto plástico
! 1 : Superfície límite de Coulomb (MC_plast)
! 2 : Superficie de tracción (vertice)
! 3 : Punto plástico en CAP (CAP_plast)
! 4 : Punto plástico en HS_CAP (HS_CAP_plast)
! 5 : Punto plástico en HS (HS_plast)
!Contenido de StVar0
! I/O Type
! 1: gamma_p I/O R Variable de endurecimiento HS
! 2: p_p I/O R Variable de endurecimiento CAP
If (IDTask .EQ. 1) Then !Inicializo variables
!Parámetros y variables de estado
G = Props(1)
xNu = Props(2)
E50_ref = Props(3)
Eoed_ref = Props(4)
Eur_ref = Props(5)
xm = Props(6)
c = Props(7)
phi = Props(8)*Rad
psi = Props(9)*Rad
xNu_ur = Props(10)
p_ref = Props(11)
xK0_nc = Props(12)
R_f = Props(13)
sigma_trac = Props(14)
If (iStep .EQ. 0) Then
!Cálculo variables endurecimiento en primer passo de cálculo
iOpt = 0
Call PrnSig(iOpt, Sig0, xN10, xN20, xN30, S10, S20, S30, P, Q)
Prs_0(1) = S10*(-1d0)
Prs_0(2) = S20*(-1d0)
Prs_0(3) = S30*(-1d0)
!Superficie HS (StVar0(1))
qf1 = max(c*cos(phi)/sin(phi)+Prs_0(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs_0(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf) Then
f13_hs_fac1 = ((Prs_0(1)-Prs_0(3))/(qa-(Prs_0(1)-Prs_0(3))))
f13_hs_fac2 = 2d0*(Prs_0(1)-Prs_0(3))/Eur
gamma_p_0 = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2
Else
f13_hs_fac1 = qf/(qa-qf)
f13_hs_fac2 = 2d0*qf/Eur
gamma_p_0 = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2
End If
!Superficie CAP (StVar0(2))
xlam_ast = p_ref/Eoed_ref
xkap_ast = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
alfa1 = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
alfa21 = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
alfa221 = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
alfa222 = (1d0-xK0_nc)*(1d0+xNu_ur)
alfa = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0
delta = (3d0+sin(phi))/(3d0-sin(phi))
q_m_0 = Prs_0(1)+(delta-1d0)*Prs_0(2)-delta*Prs_0(3)
p_0 = (Prs_0(1)+Prs_0(2)+Prs_0(3))/3d0
p_p_0 = ((q_m_0**2d0)/(alfa**2d0)+p_0**2d0)**0.5d0
!Cálculo variables en primer paso de cálculo
StVar0(1) = max(max(StVar0(1), gamma_p_0), 1d-10)
StVar0(2) = min(StVar0(2), min((-1d0)*p_p_0,(-1d0)*c*cos(phi)/sin(phi)))
End If
!Corrección valor xNu
If (xNu .GT. 0.5d0-1d-6 .OR. xNu .LT. -1d-6) Then
Call OK_MessageBox('Invalid value for "nu" (param 2)')
iAbort = 1
Return
End If
If (xNu .LT. 0d0) Then
xNu = 0.0d0
Props(2) = xNu
End If
If (xNu .GT. 0.495d0) Then
xNu = 0.495d0
Props(2) = xNu
End If
Call CopyRVec(StVar0, StVar, nStatV) !Copio StVar0 en StVar
End If ! IDTask = 1
If (IDTask .EQ. 2) Then !Cálculo de tensiones
Write(1,*) 'IDTask:', IDTask, 'iEl:', iEl, 'Int:', Int, 'iStep:', iStep, 'iTer:', iTer
Call WriIvl(ideb, 'IsUndr', IsUndr)
!Parámetros y variables de estado
Call CopyRVec(StVar0, StVar, nStatV) !Copio StVar0 en StVar
G = Props(1)
xNu = Props(2)
E50_ref = Props(3)
Eoed_ref = Props(4)
Eur_ref = Props(5)
xm = Props(6)
c = Props(7)
phi = Props(8)*Rad
psi = Props(9)*Rad
xNu_ur = Props(10)
p_ref = Props(11)
xK0_nc = Props(12)
R_f = Props(13)
sigma_trac = Props(14)
gamma_p_E = StVar0(1)
p_p_E = (-1d0)*StVar0(2)
Call WriVal(ideb, 'gamma_p_E', gamma_p_E)
Call WriVal(ideb, 'p_p_E', p_p_E)
!Construcción matriz D6
iOpt = 0
Call PrnSig(iOpt, Sig0, xN10, xN20, xN30, S10, S20, S30, P0, Q0)
Prs_0(1) = S10*(-1d0)
Prs_0(2) = S20*(-1d0)
Prs_0(3) = S30*(-1d0)
Call WriVec(ideb, 'Prs_0', Prs_0, 3)
Call Matriz6(Props, Prs_0, D, IsUndr)
Call WriMat(ideb, 'D6' , D , 6, 6, 6)
!Paso de cálculo elástico
Call MatVec(D, 6, dEps, 6, dSig)
Call AddVec(Sig0, dSig, 1d0, 1d0, 6, Sig )
Call WriVec(ideb,'dEps', dEps, 6)
Call WriVec(ideb,'Sig0', Sig0,6)
Call WriVec(ideb,'Sig', Sig,6)
iOpt = 1
Call PrnSig(iOpt, Sig, xN1, xN2, xN3, S1, S2, S3, P, Q)
!Sig : compresión negativa (S1<S2<S3)
!Prs, Prs_E: compresión positiva (S1>S2>S3)
Prs_E(1) = S1*(-1d0)
Prs_E(2) = S2*(-1d0)
Prs_E(3) = S3*(-1d0)
p_E = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
Call WriVec(ideb, 'Prs_E', Prs_E, 3)
!Cálculo superficies de fluencia f13_hs, f13_mc y f_cap en Prs_E
!Cálculo f13_hs en Prs_E
qf1 = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf_E = qf1*qf2
qa = qf_E/R_f
E_fac1 = max(c*cos(phi)+Prs_E(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then
f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
Else
f13_hs = 1d0
End If
Call WriVal(ideb, 'Prs_E(1)-Prs_E(3)', Prs_E(1)-Prs_E(3))
Call WriVal(ideb, 'qa', qa)
Call WriVal(ideb, 'qa-(Prs_E(1)-Prs_E(3))', qa-(Prs_E(1)-Prs_E(3)))
Call WriVal(ideb, 'f13_hs', f13_hs)
!Cálculo f13_mc en PrS_E
f13_mc_fac1 = (1d0/2d0)*(Prs_E(1)-Prs_E(3))
f13_mc_fac2 = (1d0/2d0)*(Prs_E(1)+Prs_E(3))*sin(phi)
f13_mc = f13_mc_fac1-f13_mc_fac2-c*cos(phi)
Call WriVal(ideb, 'f13_mc', f13_mc)
!Cálculo f_cap en Prs_E
xlam_ast = p_ref/Eoed_ref
xkap_ast = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
alfa1 = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
alfa21 = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
alfa221 = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
alfa222 = (1d0-xK0_nc)*(1d0+xNu_ur)
alfa = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0
delta = (3d0+sin(phi))/(3d0-sin(phi))
q_m = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
p = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+(p**2d0)-(p_p_E**2d0)
f_cap_E = f_cap
Call WriVal(ideb, 'f_cap', f_cap)
! Call WriVal(ideb, 'xK0_nc', xK0_nc)
! Call WriVal(ideb, 'alfa', alfa)
! Call WriVal(ideb, 'xH', xH)
!Cálculo f_trac en Prs_E
f13_trac = sigma_trac-Prs_E(3)
Call WriVal(ideb, 'f13_trac', f13_trac)
!Cálculo tolerancias
f13_hs_TOL = ttol*abs(f13_hs)+ttol
f13_mc_TOL = ttol*abs(f13_mc)+ttol
f_cap_TOL = ttol*abs(f_cap)+ttol
f13_trac_TOL = ttol*abs(f13_trac)+ttol
!Call WriVal(ideb, 'f13_hs_TOL', f13_hs_TOL)
!Call WriVal(ideb, 'f13_mc_TOL', f13_mc_TOL)
!Call WriVal(ideb, 'f_cap_TOL', f_cap_TOL)
!Call WriVal(ideb, 'f13_trac_TOL', f13_trac_TOL)
!Comprobación si el paso de scálculo es elástico
If (f13_hs .LT. f13_hs_TOL .AND. &
f13_mc .LT. f13_mc_TOL .AND. &
f13_trac .LT. f13_trac_TOL .AND. &
f_cap .LT. f_cap_TOL) Then
write(1,*) "PASO DE CÁLCULO ELÁSTICO, DENTRO DE f13_hs, f13_mc, f_cap y f13_trac"
Prs(1) = Prs_E(1)
Prs(2) = Prs_E(2)
Prs(3) = Prs_E(3)
ipl = 0
GOTO 800
End If
!Comprobación si Prs_E está fuera de f13_hs
If (f13_hs .GE. f13_hs_TOL) Then
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then
write(1,*) "ESTOY FUERA DE f13_hs"
!Cálculo Prs en superficies de HS y endurecimiento gamma_p (StVar)
iArea=1 !Calculo en plano f13_hs
Call HS_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_hs, AREA=2"
iArea = 2 !Tr. Compresion (estoy en eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_hs, AREA=3"
iArea = 3 !Tr. Extension (estoy en eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call HS_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS, AREA 2,3', Prs, 3)
End If
Else
Prs(1) = Prs_E(1)
Prs(2) = Prs_E(2)
Prs(3) = Prs_E(3)
End If
!Comprobación si nuevas tensiones Prs no superan MC
f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
f13_mc = f13_mc_fac1-f13_mc_fac2-c*cos(phi)
Call WriVal(ideb, 'f13_mc', f13_mc)
If (f13_mc .GE. f13_mc_TOL) Then
!No sirve Prs anterior
!Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Cálculo Prs en superficies MC (no tienen endurecimiento)
iArea = 1
write(1,*) "ESTOY FUERA DE f13_mc"
Call MC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_mc, AREA=2"
iArea = 2 !Tr. Compresion (estoy en eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_mc, AREA=3"
iArea = 3 !Tr. Extension (estoy en eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call MC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC, AREA 2,3', Prs, 3)
End If
End If
!Comprobación si nuevas tensiones Prs no superan TRAC
f13_trac = sigma_trac-Prs(3)
Call WriVal(ideb, 'f13_trac', f13_trac)
If (f13_trac .GE. f13_trac_TOL) Then
!No sirve Prs anterior
!Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Cálculo Prs en superficie TRAC (no tienen endurecimiento)
iArea = 1
write(1,*) "ESTOY FUERA DE f13_trac"
Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_trac, AREA=2"
iArea = 2 !Tr. Compresion (eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f12, AREA=3"
iArea = 3 !Tr. Extension (eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC, AREA 2,3', Prs, 3)
End If
!Comprobación si nuevas tensiones Prs no superan HS
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
Eur = Eur_ref*(E_fac1/E_fac2)**xm
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
If (abs(Prs(1)-Prs(3)) .LT. qf) Then
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
Else
f13_hs = 1d0
End If
If (f13_hs .GE. f13_hs_TOL) Then
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then
!No sirve Prs anterior
!Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Cálculo Prs en superficies HS y TRAC (no tienen endurecimiento)
iArea = 1
write(1,*) "ESTOY FUERA DE f13_hs y f13_trac"
Call HS_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS Y TRAC', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_hs y f13_trac, AREA=2"
iArea = 2 !Tr. Compresion (estoy en eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_hs y f13_trac, AREA=3"
iArea = 3 !Tr. Extension (estoy en eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call HS_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS Y TRAC, AREA 2,3', Prs, 3)
End If
End If
!Comprobación si nuevas tensiones Prs no superan MC
f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
f13_mc = f13_mc_fac1-f13_mc_fac2-c*cos(phi)
Call WriVal(ideb, 'f13_mc', f13_mc)
If (f13_mc .GE. f13_mc_TOL) Then
!No sirve Prs anterior
!Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Cálculo Prs en superficies MC y TRAC (no tienen endurecimiento)
iArea = 1
write(1,*) "ESTOY FUERA DE f13_mc y f13_trac"
Call MC_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y TRAC', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
iArea = 2 !Tr. Compresion (estoy en eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_mc y f13_trac, AREA=3"
iArea = 3 !Tr. Extension (estoy en eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call MC_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y TRAC, AREA 2,3', Prs, 3)
End If
End If
End If
End If
!Compruebo que las nuevas tensiones Prs (de HS o MC) no superan f_cap (alfa calculado antes)
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0
Call WriVal(ideb, 'f_cap', f_cap)
If (f_cap .GE. f_cap_TOL) Then
!No sirve Prs anterior
!Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Calculo Prs en superficie CAP y endurezco p_p
write(1,*) "ESTOY FUERA DE f_cap"
If (f_cap_E .GT. f_cap_TOL) Then
Call CAP_plast(Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
End If
!Compruebo que las nuevas tensiones Prs no superan f13_hs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
Eur = Eur_ref*(E_fac1/E_fac2)**xm
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
If (abs(Prs(1)-Prs(3)) .LT. qf) Then
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
Else
f13_hs = 1d0
End If
Call WriVal(ideb, 'f13_hs', f13_hs)
If (f13_hs .GE. f13_hs_TOL) Then
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then
!Mixed Return HS-CAP
!No sirve Prs anterior
!Desendurecimiento CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Calculo Prs en superficie HS y CAP y endurecimiento gamma_p y p_p
iArea = 1
write(1,*) "ESTOY FUERA DE f13_hs y f_cap"
Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS y CAP', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=2"
iArea = 2 !Tr. Compresion (eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=3"
iArea = 3 !Tr. Extension (eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS y CAP, AREA 2,3', Prs, 3)
End If
End If
!Compruebo que las nuevas tensiones Prs no superan MC
f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
f13_mc = f13_mc_fac1-f13_mc_fac2-c*cos(phi)
Call WriVal(ideb, 'f13_mc', f13_mc)
If (f13_mc .GE. f13_mc_TOL) Then
!Returning Map mixto MC-CAP
!No sirve Prs anterior
!Desendurecimiento HS y CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Calculo Prs en superficies CAP y MC y endurecimiento gamma_p
iArea = 1
write(1,*) "ESTOY FUERA DE f13_mc y f_cap"
Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y CAP', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=2"
iArea = 2 !Tr. Compresion (eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=3"
iArea = 3 !Tr. Extension (eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y CAP, AREA 2,3', Prs, 3)
End If
End If
End If
End If
Else
!Compruebo si Prs_E está fuera de f_cap (alfa calculado antes)
q_m = Prs_E(1)+(delta-1)*Prs_E(2)-delta*Prs_E(3)
p = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0
If (f_cap .GE. f_cap_TOL) Then
!Calculo Prs en superficie CAP y endurezco p_p
write(1,*) "ESTOY FUERA DE f_cap"
Call CAP_plast(Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Compruebo que las nuevas tensiones Prs no superan f13_hs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
Eur = Eur_ref*(E_fac1/E_fac2)**xm
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
If (abs(Prs(1)-Prs(3)) .LT. qf) Then
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
Else
f13_hs = 1d0
End If
Call WriVal(ideb, 'f13_hs', f13_hs)
If (f13_hs .GE. f13_hs_TOL) Then
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then
!Returning Map mixto HS-CAP
!No sirve Prs anterior
!Desendurecimiento CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Calculo Prs en superficie HS y CAP y endurecimiento gamma_p y p_p
iArea = 1
write(1,*) "ESTOY FUERA DE f13_hs y f_cap"
Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=2"
iArea = 2 !Tr. Compresion (eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=3"
iArea = 3 !Tr. Extension (eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS, AREA 2,3', Prs, 3)
End If
End If
End If
!Compruebo que las nuevas tensiones Prs no superan MC
f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
f13_mc = f13_mc_fac1-f13_mc_fac2-c*cos(phi)
Call WriVal(ideb, 'f13_mc', f13_mc)
If (f13_mc .GE. f13_mc_TOL) Then
!Returning Map mixto MC-CAP
!No sirve Prs anterior
!Desendurecimiento HS y CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
!Calculo Prs en superficies CAP y MC y endurecimiento gamma_p
iArea = 1
write(1,*) "ESTOY FUERA DE f13_mc y f_cap"
Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=2"
iArea = 2 !Tr. Compresion (eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=3"
iArea = 3 !Tr. Extension (eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC, AREA 2,3', Prs, 3)
End If
End If
End If
End If
!Cálculo f13_trac, f13_hs, f_cap y f13_mc con Prs_E
qf1 = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf_E = qf1*qf2
qa = qf_E/R_f
E_fac1 = max(c*cos(phi)+Prs_E(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
Eur = Eur_ref*(E_fac1/E_fac2)**xm
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then
f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
Else
f13_hs = 1d0
End If
f13_trac = sigma_trac-Prs_E(3)
If (f13_trac .GE. f13_trac_TOL .AND. &
f13_hs .LT. f13_hs_TOL) Then
!Calculo Prs en superficies TRAC
iArea = 1
write(1,*) "ESTOY FUERA DE f13_trac"
Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!Comprobación cambios en el orden de las tensiones principales
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC', Prs, 3)
If (Prs(2) .LT. Prs(3)) Then
write(1,*) "ESTOY FUERA DE f13_trac, AREA=2"
iArea = 2 !Tr. Compresion (eje S2=S3)
End If
If (Prs(1) .LT. Prs(2)) Then
write(1,*) "ESTOY FUERA DE f12, AREA=3"
iArea = 3 !Tr. Extension (eje S1=S2)
End If
If (iArea .NE. 1) Then !Returning Map en eje
Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC, AREA 2,3', Prs, 3)
End If
End If
p = (Prs(1)+Prs(2)+Prs(3))/3d0
If (p .LT. (-1d0)*sigma_trac) Then !Si p<-sigma_trac quiere decir que se ha hecho
!returning map en cono complementario
Prs(1) = (-1d0)*sigma_trac
Prs(2) = (-1d0)*sigma_trac
Prs(3) = (-1d0)*sigma_trac
ipl=2
End If
800 Continue
Write(1,*) "Las tensiones ppales al final de este paso de cálculo son Prs="
Call WriVec(1, 'Prs', Prs, 3)
Call WriVal(ideb, 'gamma_p_E', StVar0(1))
Call WriVal(ideb, 'gamma_p', StVar(1))
Call WriVal(ideb, 'p_p_E', (-1d0)*StVar0(2))
Call WriVal(ideb, 'p_p', (-1d0)*StVar(2))
!Se vuelve a las tensiones en coordenadas cartesianas
!Prs: tracción negativa, compresión positiva
S1 = (-1d0)*Prs(1)
S2 = (-1d0)*Prs(2)
S3 = (-1d0)*Prs(3)
Call CarSig(S1, S2, S3, xN1, xN2, xN3, Sig)
!Sig: tracción positiva, compresión negativa
! Write(1,*) "Las tensiones al final de este paso de cálculo son Sig"
! Call WriVec(ideb, 'Sig', Sig, 6)
E_fac1 = max(c*cos(phi)+Prs_0(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
Eur = Eur_ref*(E_fac1/E_fac2)**xm
G = Eur/(2d0*(1d0+xNu_ur))
xNu_U = 0.495d0
Fac = (1d0+xNu_U)/(1d0-2d0*xNu_U)-(1d0+xNu_ur)/(1d0-2d0*xNu_ur)
If (IsUndr .EQ. 1) Then
BulkW = 2d0*G/3d0*Fac
Else
BulkW = 0d0
End If
!Caso no drenado
If (IsUndr .EQ. 1) Then
dEpsV = dEps(1)+dEps(2)+dEps(3)
dSwp = BulkW*dEpsV
Swp = Swp0+dSwp
Else
Swp = Swp0
End If
End If ! IDTask = 2
If (IDTask .EQ. 3 .OR. &
IDTask .EQ. 6 ) Then ! Cálculo de la matriz D
!Construcción matriz D
iOpt = 1
Call PrnSig(iOpt, Sig0, xN1, xN2, xN3, S10, S20, S30, P0, Q0)
Prs_0(1) = S10*(-1d0)
Prs_0(2) = S20*(-1d0)
Prs_0(3) = S30*(-1d0)
Call Matriz6(Props, Prs_0, D, IsUndr)
BulkW = 0
If (IsUndr .EQ. 1) Then
! BulkW = ...
E_fac1 = max(c*cos(phi)+Prs_0(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
Eur = Eur_ref*(E_fac1/E_fac2)**xm
G = Eur/(2d0*(1d0+xNu_ur))
xNu_U = 0.495d0
Fac = (1d0+xNu_U)/(1d0-2d0*xNu_U)-(1d0+xNu_ur)/(1d0-2d0*xNu_ur)
Fac = 2d0*G/3d0*Fac
BulkW = Fac
End If
End If ! IDTask = 3, 6
If (IDTask .EQ. 4) Then !Número de variables de estado
nStat = nStatV
End If ! IDTask = 4
If (IDTask .Eq. 5) Then !Tipo de matriz
NonSym = 0 ! 1 para matrices D no simétricas
iStrsDep = 1 ! 1 para matrices D que dependen de las tensiones
iTang = 0 ! 1 para matrices D tangentes
iTimeDep = 0 ! 1 para matrices D que dependen del tiempo
End If ! IDTask = 5
End ! MyModel1
____________________________________________________________________________________________________
HS_CAP.F90 (SUBROUTINE)
_____________________________________________________________________________________________________
!C***********************************************************************
Subroutine HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!C***********************************************************************
Implicit Double Precision (A-H,O-Z)
!
! Rutina para el Returning Map en superficie HS
! Tensión de compresión positiva!
!
! I/O Tipo
! iArea I I 1 : Superficie f13_hs Prs_E(1) > Prs_E(2) > Prs_E(3)
! 2 : Eje entre f13_hs y f12_hs Prs_E(1) > Prs_E(2) = Prs_E(3)
! 3 : Eje entre f13_hs y f23_hs Prs_E(1) = Prs_E(2) > Prs_E(3)
!
! Local:
! f13_hs : superficie HS para Prs_E(1) > Prs_E(2) > Prs_E(3)
! f12_hs : superficie HS para Prs_E(1) > Prs_E(2) = Prs_E(3)
! f23_hs : superficie HS para Prs_E(1) = Prs_E(2) > Prs_E(3)
! xLam13_hs : multiplicador plástico f13_hs
! xLam12_hs : multiplicador plástico f12_hs
! xLam23_hs : multiplicador plástico f23_hs
! xLam_cap : multiplicador plástico f_cap
! D : matriz elástica
! df13_hsmdS(3) : Primera derivada función de fluencia f13_hs/dS
! df12_hsmdS(3) : Primera derivada función de fluencia f12_hs/dS
! df23_hsmdS(3) : Primera derivada función de fluencia f23_hs/dS
! df_capmdS(3) : Primera derivada función de fluencia f_cap/dS
! df13_hsmdq : Primera derivada función de fluencia f13_hs/dq
! df12_hsmdq : Primera derivada función de fluencia f12_hs/dq
! df23_hsmdq : Primera derivada función de fluencia f23_hs/dq
! df_capmdq : Primera derivada función de fluencia f_cap/dq
! dg13_hsmdS(3) : Primera derivada potencial plástico g13_hs/dS
! dg12_hsmdS(3) : Primera derivada potencial plástico g12_hs/dS
! dg23_hsmdS(3) : Primera derivada potencial plástico g23_hs/dS
! dg_capmdS(3) : Primera derivada potencial plástico g_cap/dS
! d2g13_hsmdS2(3,3) : Segunda derivada potencial plástico g13_hs/dS2
! d2g12_hsmdS2(3,3) : Segunda derivada potencial plástico g12_hs/dS2
! d2g23_hsmdS2(3,3) : Segunda derivada potencial plástico g23_hs/dS2
! d2g_capmdS2(3,3) : Segunda derivada potencial plástico g_cap/dS2
! d2g13_hsmdSdq(3) : Derivada dS y dq del potencial plástico g13_hs
! d2g12_hsmdSdq(3) : Derivada dS y dq del potencial plástico g12_hs
! d2g23_hsmdSdq(3) : Derivada dS y dq del potencial plástico g23_hs
! d2g_capmdSdq(3) : Derivada dS y dq del potencial plástico g_cap
! D_dgij : Producto de D por dgij_hs/dS
! D_dgij_2 : Producto de D por d2gij_hs/dS2
! HH_dfij : Producto de D por dfij_hs/dS
! thetaij : Operación dgij_hs/dS+d2gij_hs/dS2*xLamij_hs*xh_hs
! D_thetaij : Productor de D por thetaij
! HH : Inversa de la matriz Unity+xLamij_hs*D_gij_2
Dimension Prs_E(3), Prs(3), Prs_res1(3), Prs_res2(3), Prs_res3(3), Prs_0(3), xN10(3), xN20(3), xN30(3)
Dimension Prs_2sup(3), Prs_3sup(3), r(3)
Dimension Props(*)
Dimension StVar0(*), StVar(*)
Dimension df13_hsmdS(3), df12_hsmdS(3), df23_hsmdS(3), df_capmdS(3)
Dimension dg13_hsmdS(3), dg12_hsmdS(3), dg23_hsmdS(3), dg_capmdS(3)
Dimension d2g13_hsmdS2(3,3), d2g12_hsmdS2(3,3), d2g23_hsmdS2(3,3), d2g_capmdS2(3,3)
Dimension d2g13_hsmdSdq(3), d2g12_hsmdSdq(3), d2g23_hsmdSdq(3), d2g_capmdSdq(3)
Dimension D_dg13(3), D_dg12(3), D_dg23(3), D_dgcap(3)
Dimension D_dg13_2(3,3), D_dg12_2(3,3), D_dg23_2(3,3), D_dgcap_2(3,3)
Dimension HH_df13(3), HH_df12(3), HH_df23(3), HH_dfcap(3)
Dimension theta13(3), theta12(3), theta23(3), thetacap(3)
Dimension D_theta13(3), D_theta12(3), D_theta23(3), D_thetacap(3)
Dimension Unity(3,3), D(3,3), HH(3,3), HH_inv(3,3), xMat(3,3), xMat_inv(3,3), f_b(3), d_lam(3)
Call MZeroR(Prs_res1, 3)
Call MZeroR(Prs_res2, 3)
Call MZeroR(Prs_res3, 3)
Call MZeroR(Prs_2sup, 3)
Call MZeroR(Prs_3sup, 3)
Call MZeroR(df13_hsmdS, 3)
Call MZeroR(dg13_hsmdS, 3)
Call MZeroR(d2g13_hsmdS2, 9)
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(df12_hsmdS, 3)
Call MZeroR(dg12_hsmdS, 3)
Call MZeroR(d2g12_hsmdS2, 9)
Call MZeroR(d2g12_hsmdSdq, 3)
Call MZeroR(df23_hsmdS, 3)
Call MZeroR(dg23_hsmdS, 3)
Call MZeroR(d2g23_hsmdS2, 9)
Call MZeroR(d2g23_hsmdSdq, 3)
Call MZeroR(df_capmdS, 3)
Call MZeroR(dg_capmdS, 3)
Call MZeroR(d2g_capmdS2, 9)
Call MZeroR(d2g_capmdSdq, 3)
Call MZeroR(D_dg13_2, 9)
Call MZeroR(D_dg12_2, 9)
Call MZeroR(D_dg23_2, 9)
Call MZeroR(D_dgcap_2, 9)
Call MZeroR(D_dg13, 3)
Call MZeroR(D_dg12, 3)
Call MZeroR(D_dg23, 3)
Call MZeroR(D_dgcap, 3)
Call MZeroR(HH_df13, 3)
Call MZeroR(HH_df12, 3)
Call MZeroR(HH_df23, 3)
Call MZeroR(HH_dfcap, 3)
Call MZeroR(theta13, 3)
Call MZeroR(theta12, 3)
Call MZeroR(theta23, 3)
Call MZeroR(thetacap, 3)
Call MZeroR(D_theta13, 3)
Call MZeroR(D_theta12, 3)
Call MZeroR(D_theta23, 3)
Call MZeroR(D_thetacap, 3)
Call MZeroR(HH, 3)
Call MZeroR(HH_inv, 3)
Call MZeroR(xMat, 9)
Call MZeroR(xMat_inv, 9)
Call MZeroR(f_b, 3)
Call MZeroR(d_lam, 3)
ipl = 5
it = 0
itmax = 20
ttol = 1d-9
pi = 3.1415926535d0
Rad = pi/180d0
ideb = -1 !para debugging
Call MZeroR(Unity, 9)
Unity(1,1) = 1d0
Unity(2,2) = 1d0
Unity(3,3) = 1d0
G = Props(1)
xNu = Props(2)
E50_ref = Props(3)
Eoed_ref = Props(4)
Eur_ref = Props(5)
xm = Props(6)
c = Props(7)
phi = Props(8)*Rad
psi = Props(9)*Rad
xNu_ur = Props(10)
p_ref = Props(11)
xK0_nc = Props(12)
R_f = Props(13)
sigma_trac= Props(14)
xKs = Eur_ref/((3d0*(1d0-2d0*xNu_ur)))
xKsmKc = (Eur_ref/Eoed_ref)*(xK0_nc/((1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)))
xH = xKs/(xKsmKc-1d0)
gamma_p_E = StVar0(1)
p_p_E = (-1d0)*StVar0(2)
dxLam13_hs = 0d0
dxLam12_hs = 0d0
dxLam23_hs = 0d0
dxLam_cap = 0d0
xLam13_hs = 0d0
xLam12_hs = 0d0
xLam23_hs = 0d0
xLam_cap = 0d0
!Construcción matriz D
iOpt = 0
Call PrnSig(iOpt, Sig0, xN10, xN20, xN30, S10, S20, S30, P0, Q0)
Prs_0(1) = S10*(-1d0)
Prs_0(2) = S20*(-1d0)
Prs_0(3) = S30*(-1d0)
Call Matriz3(Props, Prs_0, D, IsUndr)
Select Case (IArea)
Case (1) ! 1 : Superficie f13_hs y f_cap (Sig1 > Sig2 > Sig3)
!STARTING POINT
!Cálculo valor de f13_hs y f_cap en Prs_E
qf1 = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs_E(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
xlam_ast = p_ref/Eoed_ref
xkap_ast = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
alfa1 = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
alfa21 = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
alfa221 = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
alfa222 = (1d0-xK0_nc)*(1d0+xNu_ur)
alfa = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0
delta = (3d0+sin(phi))/(3d0-sin(phi))
q_m = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
p = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0
Call WriVal(ideb, 'qa', qa)
Call WriVal(ideb, 'qf', qf)
Call WriVal(ideb, 'Prs_E(1)-Prs_E(3)', Prs_E(1)-Prs_E(3))
Call WriVal(ideb, 'qa-(Prs_E(1)-Prs_E(3))', qa-(Prs_E(1)-Prs_E(3)))
!Cálculo tolerancia de f13_hs y f_cap
f13_hs_TOL = ttol*abs(f13_hs)+ttol
f_cap_TOL = ttol*abs(f_cap)+ttol
Call WriVal(ideb, 'f13_hs', f13_hs)
Call WriVal(ideb, 'f_cap', f_cap)
Call WriVal(ideb, 'f13_hs_TOL', f13_hs_TOL)
Call WriVal(ideb, 'f_cap_TOL', f_cap_TOL)
!Calculo df13_hs/dS y df_cap/dS en Prs_E
df13_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
df13_hsmdS(2) = 0d0
df13_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur
df_capmdS(1) = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(2) = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(3) = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)
Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)
!Cálculo df13_hs/dq en Prs_E (q1=gamma_p, q2=p_p)
df13_hsmdq = -1d0
df_capmdq = -2d0*p_p_E
!Cálculo h en Prs_E (dq1=dgamma_p=dlam1*h1, dq2=dp_p=dlam2*h2)
xh_hs = 1d0
xh_cap = 2d0*xH*p*(max(Prs_E(3)+c*cos(phi)/sin(phi),c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
!Cálculo dg13_hs/dS y dg_cap/dS en Prs_E
sin_phi_m = (Prs_E(1)-Prs_E(3))/(Prs_E(1)+Prs_E(3)-2d0*c*cos(phi)/sin(phi))
sin_phi_cv = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
sin_psi_m = 0d0
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
psi_m = psi
sin_psi_m = sin(psi_m)
ElseIf (phi .EQ. 0d0) Then
sin_psi_m = 0d0
End If
dg13_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg13_hsmdS(2) = 0d0
dg13_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg_capmdS(1) = df_capmdS(1)
dg_capmdS(2) = df_capmdS(2)
dg_capmdS(3) = df_capmdS(3)
Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)
Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)
!Cálculo d2g13_hs/dS2 y d2gcap/dS2 en Prs_E
Call MZeroR(d2g13_hsmdS2, 9)
d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(2,1) = d2g_capmdS2(1,2)
d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(3,1) = d2g_capmdS2(1,3)
d2g_capmdS2(3,2) = d2g_capmdS2(2,3)
d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)
!Cálculo d2g13_hs/dSdq en Prs_E
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(d2g_capmdSdq, 3)
!Cálculo de los multiplicadores plásticos
Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)
Call DInProd(df13_hsmdS, D_dg13, 3, df13_D_dg13)
xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs
Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)
Call DInProd(df13_hsmdS, D_dgcap, 3, df13_D_dgcap)
xOmega12 = df13_D_dgcap
Call DInProd(df_capmdS, D_dg13,3, dfcap_D_dg13)
xOmega21 = dfcap_D_dg13
Call DInProd(df_capmdS, D_dgcap,3, dfcap_D_dgcap)
xOmega22 = dfcap_D_dgcap-df_capmdq*xh_cap
! Call WriVal(ideb, 'df13_D_dg13', df13_D_dg13)
! Call WriVal(ideb, 'df13_hsmdq', df13_hsmdq)
! Call WriVal(ideb, 'xh_hs', xh_hs)
! Call WriVal(ideb, 'xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs', df13_D_dg13-df13_hsmdq*xh_hs)
! Call WriVal(ideb, 'xOmega12 = df13_D_dgcap', df13_D_dgcap)
! Call WriVal(ideb, 'xOmega21 = dfcap_D_dg13p', dfcap_D_dg13)
! Call WriVal(ideb, 'dfcap_D_dgcap', dfcap_D_dgcap)
! Call WriVal(ideb, 'df_capmdq', df_capmdq)
! Call WriVal(ideb, 'xh_cap', xh_cap)
! Call WriVal(ideb, 'xOmega22 = dfcap_D_dgcap-df_capmdq*xh_cap ', dfcap_D_dgcap-df_capmdq*xh_cap)
!Resuolución sistema y cálculo dxLam13_hs y dxLam_cap
dxLam13_hs = (f13_hs*xOmega22-f_cap*xOmega12)/(xOmega11*xOmega22-xOmega12*xOmega21)
dxLam_cap = (f_cap*xOmega11-f13_hs*xOmega21)/(xOmega11*xOmega22-xOmega12*xOmega21)
xLam13_hs = dxLam13_hs
xLam_cap = dxLam_cap
Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)
Call WriVal(ideb, 'dxLam_cap', dxLam_cap)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
!Actualizo variables
Call AddVec(D_dg13, D_dgcap, -xLam13_hs, -xLam_cap, 3, Prs_2sup)
Call AddVec(Prs_E, Prs_2sup, 1d0, 1d0, 3, Prs)
gamma_p = gamma_p_E + xLam13_hs*xh_hs
p_p = p_p_E + xLam_cap*xh_cap
!Cálculo valor de f13_hs y f_cap con la nueva Prs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0)) !OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0
!EULER BACKWARD
Do while ((abs(f13_hs) .GT. f13_hs_TOL .OR. &
abs(f_cap) .GT. f_cap_TOL) .AND. &
it .LE. itmax)
!Construcción matriz D
!Call Matriz3(Props, Prs, D, IsUndr)
!Calculo df13_hs/dS y df_cap/dS en Prs
Call WriVal(ideb, 'qa', qa)
Call WriVal(ideb, 'qf', qf)
Call WriVal(ideb, 'Prs(1)-Prs(3)', Prs(1)-Prs(3))
Call WriVal(ideb, 'qa-(Prs(1)-Prs(3))', qa-(Prs(1)-Prs(3)))
df13_hsmdS(1) = ((2d0*(qa**2d0))/Ei)*(1/((qa-(Prs(1)-Prs(3)))**2d0))-2d0/Eur
df13_hsmdS(2) = 0d0
df13_hsmdS(3) = -((2d0*(qa**2d0))/Ei)*(1/((qa-(Prs(1)-Prs(3)))**2d0))+2d0/Eur
df_capmdS(1) = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(2) = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(3) = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)
Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)
!Cálculo df13_hs/dq en Prs (q1=gamma_p, q2=p_p)
df13_hsmdq = -1d0
df_capmdq = -2d0*p_p
Call WriVal(ideb, 'df_capmdq', df_capmdq)
!Cálculo h en Prs (dq=dlam1*h1 y dq2=dlam2*h2)
xh_hs = 1d0
xh_cap = 2d0*xH*p*(max(Prs(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
Call WriVal(ideb, 'xh_cap', xh_cap)
!Cálculo dg13_hs/dS y dg_cap/dS en Prs
sin_phi_m = (Prs(1)-Prs(3))/(Prs(1)+Prs(3)-2d0*c*cos(phi)/sin(phi))
sin_phi_cv = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
sin_psi_m = 0d0
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
psi_m = psi
sin_psi_m = sin(psi_m)
ElseIf (phi .EQ. 0d0) Then
sin_psi_m = 0d0
End If
dg13_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg13_hsmdS(2) = 0d0
dg13_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg_capmdS(1) = df_capmdS(1)
dg_capmdS(2) = df_capmdS(2)
dg_capmdS(3) = df_capmdS(3)
Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)
Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)
!Cálculo d2g13_hs/dS2 y d2g_cap/dS2 en Prs
Call MZeroR(d2g13_hsmdS2, 9)
d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(2,1) = d2g_capmdS2(1,2)
d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(3,1) = d2g_capmdS2(1,3)
d2g_capmdS2(3,2) = d2g_capmdS2(2,3)
d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)
!Cálculo d2g13_hs/dSdq en Prs
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(d2g_capmdSdq, 3)
!Cálculo del multiplicador plástico
Call MatMat(D, 3, d2g13_hsmdS2, 3, 3, 3, 3, D_dg13_2, 3)
Call MatMat(D, 3, d2g_capmdS2, 3, 3, 3, 3, D_dgcap_2, 3)
HH_inv = Unity + xLam13_hs*D_dg13_2 + xLam_cap*D_dgcap_2
Call MatInvPiv(HH_inv, HH, 3)
Call AddVec(dg13_hsmdS, d2g13_hsmdSdq, 1d0, xLam13_hs*xh_hs, 3, theta13)
Call AddVec(dg_capmdS, d2g_capmdSdq, 1d0, xLam_cap*xh_cap, 3, thetacap)
Call MatMat(HH, 3, df13_hsmdS, 3, 3, 1, 3, HH_df13,3)
Call MatMat(HH, 3, df_capmdS, 3, 3, 1, 3, HH_dfcap,3)
Call MatMat(D, 3, theta13, 3, 3, 1, 3, D_theta13,3)
Call MatMat(D, 3, thetacap, 3, 3, 1, 3, D_thetacap,3)
Call DInProd(D_theta13, HH_df13, 3, D_theta13_HH_df13)
Call DInProd(D_thetacap, HH_df13, 3, D_thetacap_HH_df13)
Call DInProd(D_theta13, HH_dfcap, 3, D_theta13_HH_dfcap)
Call DInProd(D_thetacap, HH_dfcap, 3, D_thetacap_HH_dfcap)
xOmega11 = D_theta13_HH_df13-df13_hsmdq*xh_hs
xOmega12 = D_thetacap_HH_df13
xOmega21 = D_theta13_HH_dfcap
xOmega22 = D_thetacap_HH_dfcap-df_capmdq*xh_cap
!Call WriMat(1, 'HH', HH, 3, 3, 3)
Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)
Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)
Call AddVec(D_dg13, D_dgcap, -xLam13_hs, -xLam_cap, 3, Prs_res1)
Call AddVec(Prs_E, Prs_res1, 1d0, 1d0, 3, Prs_res2)
Call AddVec(Prs, Prs_res2, 1d0, -1d0, 3, r)
Call DInProd(HH_df13, r, 3, r_HH_df13)
Call DInProd(HH_dfcap, r, 3, r_HH_dfcap)
xOmegaf13 = f13_hs - r_HH_df13
xOmegafcap = f_cap - r_HH_dfcap
!Resolución sistema y cálculo dxLam13_hs y dxLam_cap
dxLam13_hs = (xOmegaf13*xOmega22-xOmegafcap*xOmega12)/(xOmega11*xOmega22-xOmega12*xOmega21)
dxLam_cap = (xOmegafcap*xOmega11-xOmegaf13*xOmega21)/(xOmega11*xOmega22-xOmega12*xOmega21)
xLam13_hs = xLam13_hs + dxLam13_hs
xLam_cap = xLam_cap + dxLam_cap
Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)
Call WriVal(ideb, 'dxLam_cap', dxLam_cap)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
!Actualización variables
Call AddVec(D_dg13, D_dgcap, -xLam13_hs, -xLam_cap, 3, Prs_2sup)
Call AddVec(Prs_E, Prs_2sup, 1d0, 1d0, 3, Prs)
gamma_p = gamma_p_E + xLam13_hs*xh_hs
p_p = p_p_E + xLam_cap*xh_cap
Call WriVal(ideb, 'gamma_p', gamma_p)
Call WriVal(ideb, 'p_p', p_p)
Call WriVal(ideb, 'qa', qa)
Call WriVal(ideb, 'Prs(1)-Prs(3)', Prs(1)-Prs(3))
Call WriVal(ideb, 'qa-(Prs(1)-Prs(3))', qa-(Prs(1)-Prs(3)))
Call WriVal(ideb, 'Ei', Ei)
Call WriVal(ideb, 'Eur', Eur)
Call WriVal(ideb, '(qa-(Prs(1)-Prs(3)))**2d0', (qa-(Prs(1)-Prs(3)))**2d0)
Call WriVec(ideb, 'Prs durante iteración', Prs, 3)
!Cálculo valor de f13_hs en (k+1)
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0
Call WriIvl(ideb, 'it', it)
!Actualizo iteración
it = it+1
End do
!Actualizo variables de estado
StVar(1) = gamma_p
StVar(2) = min((-1d0)*p_p, (-1d0)*p_p_E)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
If (xLam13_hs .GT. (-1d0)*ttol .AND. &
xLam_cap .GT. (-1d0)*ttol .AND. &
it .LE. itmax) Then
write(1,*)'Return mapping Ok!'
Else
write(1,*)'Numerical problem'
End If
Case (2) ! 2 : Eje entre f13_hs, f12_hs y f_cap (Prs(1) > Prs(2) = Prs(3))
!STARTING POINT
!Cálculo el valor de f13_hs, f12_hs y f_cap en Prs_E
qf1 = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs_E(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
f12_hs_fac1 = ((Prs_E(1)-Prs_E(2))/(qa-(Prs_E(1)-Prs_E(2))))
f12_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(2))/Eur
f12_hs = (2d0*qa/Ei)*f12_hs_fac1-f12_hs_fac2-gamma_p_E
xlam_ast = p_ref/Eoed_ref
xkap_ast = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
alfa1 = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
alfa21 = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
alfa221 = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
alfa222 = (1d0-xK0_nc)*(1d0+xNu_ur)
alfa = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0
delta = (3d0+sin(phi))/(3d0-sin(phi))
q_m = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
p = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0
Call WriVal(ideb, 'qa', qa)
Call WriVal(ideb, 'qf', qf)
Call WriVal(ideb, 'Prs_E(1)-Prs_E(3)', Prs_E(1)-Prs_E(3))
Call WriVal(ideb, 'qa-(Prs_E(1)-Prs_E(3))', qa-(Prs_E(1)-Prs_E(3)))
!Cálculo tolerancia f13_hs, f12_hs y f_cap en Prs_E
f13_hs_TOL = ttol*abs(f13_hs)+ttol
f12_hs_TOL = ttol*abs(f12_hs)+ttol
f_cap_TOL = ttol*abs(f_cap)+ttol
!Calculo df13_hs/dS, df12_hs/dS y df_cap/dS en Prs_E
df13_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
df13_hsmdS(2) = 0d0
df13_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur
df12_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(2)))**2d0)-2d0/Eur
df12_hsmdS(2) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(2)))**2d0)+2d0/Eur
df12_hsmdS(3) = 0d0
df_capmdS(1) = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(2) = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(3) = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)
Call WriVec(ideb, 'df12_hsmdS', df12_hsmdS, 3)
Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)
!Cálculo df13_hs/dq, df12_hs/dq y df_cap/dS en Prs_E (q1=gamma_p -para f13_hs y f12_hs-, q2=p_p)
df13_hsmdq = -1d0
df12_hsmdq = -1d0
df_capmdq = -2d0*p_p_E
Call WriVal(ideb, 'df13_hsmdq', df13_hsmdq)
Call WriVal(ideb, 'df12_hsmdq', df12_hsmdq)
Call WriVal(ideb, 'df_capmdq', df_capmdq)
!Cálculo h en Prs_E (dq=dlam1*h1 -para f13_hs y f12_hs- y dq2=dlam2*h2)
xh_hs = 1d0
xh_cap = 2d0*xH*p*(max(Prs_E(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
Call WriVal(ideb, 'xh_hs', xh_hs)
Call WriVal(ideb, 'xh_cap', xh_cap)
!Cálculo dg13_hs/dS, dg12_hs/dS y dg_cap/dS en Prs_E
sin_phi_m = (Prs_E(1)-Prs_E(3))/(Prs_E(1)+Prs_E(3)-2d0*c*cos(phi)/sin(phi))
sin_phi_cv = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
sin_psi_m = 0d0
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
psi_m = psi
sin_psi_m = sin(psi_m)
ElseIf (phi .EQ. 0d0) Then
sin_psi_m = 0d0
End If
dg13_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg13_hsmdS(2) = 0d0
dg13_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg12_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg12_hsmdS(2) = -(1d0/2d0)*(1d0+sin_psi_m)
dg12_hsmdS(3) = 0d0
dg_capmdS(1) = df_capmdS(1)
dg_capmdS(2) = df_capmdS(2)
dg_capmdS(3) = df_capmdS(3)
Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)
Call WriVec(ideb, 'dg12_hsmdS', dg12_hsmdS, 3)
Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)
!Cálculo d2g13_hs/dS2, d2g12_hs/dS2 y d2g_cap/dS2 en Prs_E
Call MZeroR(d2g13_hsmdS2, 9)
Call MZeroR(d2g12_hsmdS2, 9)
d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(2,1) = d2g_capmdS2(1,2)
d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(3,1) = d2g_capmdS2(1,3)
d2g_capmdS2(3,2) = d2g_capmdS2(2,3)
d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)
!Cálculo d2g13_hs/dSdq, d2g12_hs/dSdq y d2g_cap/dSdq en Prs_E
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(d2g12_hsmdSdq, 3)
Call MZeroR(d2g_capmdSdq, 3)
!Cálculo de los multiplicadores plásticos
Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)
Call DInProd(df13_hsmdS, D_dg13,3, df13_D_dg13)
xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs
Call MatMat(D,3,dg12_hsmdS,3,3,1,3,D_dg12,3)
Call DInProd(df13_hsmdS, D_dg12,3, df13_D_dg12)
xOmega12 = df13_D_dg12
Call MatMat(D,3,dg_capmdS,3,3,1,3,D_dgcap,3)
Call DInProd(df13_hsmdS, D_dgcap,3, df13_D_dgcap)
xOmega13 = df13_D_dgcap
Call DInProd(df12_hsmdS, D_dg13,3, df12_D_dg13)
xOmega21 = df12_D_dg13
Call DInProd(df12_hsmdS, D_dg12,3, df12_D_dg12)
xOmega22 = df12_D_dg12-df12_hsmdq*xh_hs
Call DInProd(df12_hsmdS, D_dgcap,3, df12_D_dgcap)
xOmega23 = df12_D_dgcap
Call DInProd(df_capmdS, D_dg13,3, dfcap_D_dg13)
xOmega31 = dfcap_D_dg13
Call DInProd(df_capmdS, D_dg12,3, dfcap_D_dg12)
xOmega32 = dfcap_D_dg12
Call DInProd(df_capmdS, D_dgcap,3, dfcap_D_dgcap)
xOmega33 = dfcap_D_dgcap-df_capmdq*xh_cap
!Resolución sistema y calculo dxLam13_hs, dxLam12_hs y dxLam_cap
xMat(1,1) = xOmega11
xMat(1,2) = xOmega12
xMat(1,3) = xOmega13
xMat(2,1) = xOmega21
xMat(2,2) = xOmega22
xMat(2,3) = xOmega23
xMat(3,1) = xOmega31
xMat(3,2) = xOmega32
xMat(3,3) = xOmega33
f_b(1) = f13_hs
f_b(2) = f12_hs
f_b(3) = f_cap
Call MatInvPiv(xMat, xMat_inv, 3)
Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3)
dxLam13_hs = d_lam(1)
dxLam12_hs = d_lam(2)
dxLam_cap = d_lam(3)
xLam13_hs = dxLam13_hs
xLam12_hs = dxLam12_hs
xLam_cap = dxLam_cap
Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)
Call WriVal(ideb, 'dxLam12_hs', dxLam12_hs)
Call WriVal(ideb, 'dxLam_cap', dxLam_cap)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam12_hs', xLam12_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
!Actualización variables
Call AddVec(D_dg13, D_dg12, -xLam13_hs, -xLam12_hs, 3, Prs_2sup)
Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)
Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
gamma_p = gamma_p_E + xLam13_hs*xh_hs
p_p = p_p_E + xLam_cap*xh_cap
Call WriVal(ideb, 'gamma_p', gamma_p)
Call WriVal(ideb, 'p_p', p_p)
!Cálculo valor de f13_hs, f12_hs y f_cap con la nueva Prs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = (c*cos(phi)+Prs(3)*sin(phi))
E_fac2 = max(c*cos(phi)+p_ref*sin(phi), c*cos(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p
f12_hs_fac1 = ((Prs(1)-Prs(2))/(qa-(Prs(1)-Prs(2))))
f12_hs_fac2 = 2d0*(Prs(1)-Prs(2))/Eur
f12_hs = (2d0*qa/Ei)*f12_hs_fac1-f12_hs_fac2-gamma_p
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0
!EULER BACKWARD
Do while ((abs(f13_hs) .GT. f13_hs_TOL .OR. &
abs(f12_hs) .GT. f12_hs_TOL .OR. &
abs(f_cap) .GT. f_cap_TOL) .AND. &
it .LE. itmax)
!Construcción matriz D
!Call Matriz3(Props, Prs, D, IsUndr)
Call WriVal(ideb, 'qa', qa)
Call WriVal(ideb, 'qf', qf)
Call WriVal(ideb, 'Prs(1)-Prs(3)', Prs(1)-Prs(3))
Call WriVal(ideb, 'qa-(Prs(1)-Prs(3))', qa-(Prs(1)-Prs(3)))
!Cálculo df13_hs/dS, df12_hs/dS y df_cap/dS en Prs
df13_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)-2d0/Eur
df13_hsmdS(2) = 0d0
df13_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)+2d0/Eur
df12_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(2)))**2d0)-2d0/Eur
df12_hsmdS(2) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(2)))**2d0)+2d0/Eur
df12_hsmdS(3) = 0d0
df_capmdS(1) = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(2) = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(3) = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)
Call WriVec(ideb, 'df12_hsmdS', df12_hsmdS, 3)
Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)
!Cálculo df13_hs/dq, df12_hs/dq y df_cap/dq en Prs (q1=gamma_p -para f13_hs y f12_hs-, q2=p_p)
df13_hsmdq = -1d0
df12_hsmdq = -1d0
df_capmdq = -2d0*p_p
Call WriVal(ideb, 'df13_hsmdq', df13_hsmdq)
Call WriVal(ideb, 'df12_hsmdq', df12_hsmdq)
Call WriVal(ideb, 'df_capmdq', df_capmdq)
!Cálculo h en Prs (dq=dlam1*h1 -para f13_hs y f12_hs- y dq2=dlam2*h2)
xh_hs = 1d0
xh_cap = 2d0*xH*p*(max(Prs(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
Call WriVal(ideb, 'xh_hs', xh_hs)
Call WriVal(ideb, 'xh_cap', xh_cap)
!Cálculo dg13_hs/dS, dg12_hs/dS y dg_cap/dS en Prs
sin_phi_m = (Prs(1)-Prs(3))/(Prs(1)+Prs(3)-2d0*c*cos(phi)/sin(phi))
sin_phi_cv = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
sin_psi_m = 0d0
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
psi_m = psi
sin_psi_m = sin(psi_m)
ElseIf (phi .EQ. 0d0) Then
sin_psi_m = 0d0
End If
dg13_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg13_hsmdS(2) = 0d0
dg13_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg12_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg12_hsmdS(2) = -(1d0/2d0)*(1d0+sin_psi_m)
dg12_hsmdS(3) = 0d0
dg_capmdS(1) = df_capmdS(1)
dg_capmdS(2) = df_capmdS(2)
dg_capmdS(3) = df_capmdS(3)
Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)
Call WriVec(ideb, 'dg12_hsmdS', dg12_hsmdS, 3)
Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)
!Cálculo d2g13_hs/dS2, d2g12_hs/dS2 y d2g_cap/dS2 en Prs
Call MZeroR(d2g13_hsmdS2, 9)
Call MZeroR(d2g12_hsmdS2, 9)
d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(2,1) = d2g_capmdS2(1,2)
d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(3,1) = d2g_capmdS2(1,3)
d2g_capmdS2(3,2) = d2g_capmdS2(2,3)
d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)
!Cálculo d2g13_hs/dSdq, d2g12_hs/dSdq y d2g_cap/dSdq en Prs
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(d2g12_hsmdSdq, 3)
Call MZeroR(d2g_capmdSdq, 3)
!Cálculo del multiplicador plástico
Call MatMat(D, 3, d2g13_hsmdS2, 3, 3, 3, 3, D_dg13_2, 3)
Call MatMat(D, 3, d2g12_hsmdS2, 3, 3, 3, 3, D_dg12_2, 3)
Call MatMat(D, 3, d2g_capmdS2, 3, 3, 3, 3, D_dgcap_2, 3)
HH_inv = Unity + xLam13_hs*D_dg13_2 + xLam12_hs*D_dg12_2 + xLam_cap*D_dgcap_2
Call MatInvPiv(HH_inv, HH, 3)
Call AddVec(dg13_hsmdS, d2g13_hsmdSdq, 1d0, xLam13_hs*xh_hs, 3, theta13)
Call AddVec(dg12_hsmdS, d2g12_hsmdSdq, 1d0, xLam12_hs*xh_hs, 3, theta12)
Call AddVec(dg_capmdS, d2g_capmdSdq, 1d0, xLam_cap*xh_cap, 3, thetacap)
Call MatMat(HH, 3, df13_hsmdS, 3, 3, 1, 3, HH_df13,3)
Call MatMat(HH, 3, df12_hsmdS, 3, 3, 1, 3, HH_df12,3)
Call MatMat(HH, 3, df_capmdS, 3, 3, 1, 3, HH_dfcap,3)
Call MatMat(D, 3, theta13, 3, 3, 1, 3, D_theta13,3)
Call MatMat(D, 3, theta12, 3, 3, 1, 3, D_theta12,3)
Call MatMat(D, 3, thetacap, 3, 3, 1, 3, D_thetacap,3)
Call DInProd(D_theta13, HH_df13, 3, D_theta13_HH_df13)
Call DInProd(D_theta12, HH_df13, 3, D_theta12_HH_df13)
Call DInProd(D_thetacap, HH_df13, 3, D_thetacap_HH_df13)
Call DInProd(D_theta13, HH_df12, 3, D_theta13_HH_df12)
Call DInProd(D_theta12, HH_df12, 3, D_theta12_HH_df12)
Call DInProd(D_thetacap, HH_df12, 3, D_thetacap_HH_df12)
Call DInProd(D_theta13, HH_dfcap, 3, D_theta13_HH_dfcap)
Call DInProd(D_theta12, HH_dfcap, 3, D_theta12_HH_dfcap)
Call DInProd(D_thetacap, HH_dfcap, 3, D_thetacap_HH_dfcap)
xOmega11 = D_theta13_HH_df13-df13_hsmdq*xh_hs
xOmega12 = D_theta12_HH_df13
xOmega13 = D_thetacap_HH_df13
xOmega21 = D_theta13_HH_df12
xOmega22 = D_theta12_HH_df12-df12_hsmdq*xh_hs
xOmega23 = D_thetacap_HH_df12
xOmega31 = D_theta13_HH_dfcap
xOmega32 = D_theta12_HH_dfcap
xOmega33 = D_thetacap_HH_dfcap-df_capmdq*xh_cap
Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)
Call MatMat(D, 3, dg12_hsmdS, 3, 3, 1, 3, D_dg12, 3)
Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)
Call AddVec(D_dg13, D_dg12, -xLam13_hs, -xLam12_hs, 3, Prs_res1)
Call AddVec(Prs_res1, D_dgcap, 1d0, -xLam_cap, 3, Prs_res2)
Call AddVec(Prs_E, Prs_res2, 1d0, 1d0, 3, Prs_res3)
Call AddVec(Prs, Prs_res3, 1d0,-1d0, 3, r)
Call DInProd(HH_df13, r, 3, r_HH_df13)
Call DInProd(HH_df12, r, 3, r_HH_df12)
Call DInProd(HH_dfcap, r, 3, r_HH_dfcap)
xOmegaf13 = f13_hs - r_HH_df13
xOmegaf12 = f12_hs - r_HH_df12
xOmegafcap = f_cap - r_HH_dfcap
!Resolución sistema y calculo dxLam13_hs y dxLam12_hs
xMat(1,1) = xOmega11
xMat(1,2) = xOmega12
xMat(1,3) = xOmega13
xMat(2,1) = xOmega21
xMat(2,2) = xOmega22
xMat(2,3) = xOmega23
xMat(3,1) = xOmega31
xMat(3,2) = xOmega32
xMat(3,3) = xOmega33
f_b(1) = xOmegaf13
f_b(2) = xOmegaf12
f_b(3) = xOmegafcap
Call MatInvPiv(xMat, xMat_inv, 3)
Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3)
dxLam13_hs = d_lam(1)
dxLam12_hs = d_lam(2)
dxLam_cap = d_lam(3)
xLam13_hs = xLam13_hs + dxLam13_hs
xLam12_hs = xLam12_hs + dxLam12_hs
xLam_cap = xLam_cap + dxLam_cap
Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)
Call WriVal(ideb, 'dxLam12_hs', dxLam12_hs)
Call WriVal(ideb, 'dxLam_cap', dxLam_cap)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam12_hs', xLam12_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
!Actualizo variables
Call AddVec(D_dg13, D_dg12, -xLam13_hs, -xLam12_hs, 3, Prs_2sup)
Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)
Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
gamma_p = gamma_p_E + xLam13_hs*xh_hs
p_p = p_p_E + xLam_cap*xh_cap
Call WriVal(ideb, 'gamma_p', gamma_p)
Call WriVal(ideb, 'p_p', p_p)
!Calculo el valor de f13_hs, f12_hs y f_cap en (k+1) en Prs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p
f12_hs_fac1 = ((Prs(1)-Prs(2))/(qa-(Prs(1)-Prs(2))))
f12_hs_fac2 = 2d0*(Prs(1)-Prs(2))/Eur
f12_hs = (2d0*qa/Ei)*f12_hs_fac1-f12_hs_fac2-gamma_p
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0
Call WriIvl(ideb, 'it', it)
!Actualización iteración
it = it+1
End do
!Actualización variables de estado
StVar(1) = gamma_p
StVar(2) = min((-1d0)*p_p, (-1d0)*p_p_E)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam12_hs', xLam12_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
If (xLam13_hs .GT. (-1d0)*ttol .AND. &
xLam12_hs .GT. (-1d0)*ttol .AND. &
xLam_cap .GT. (-1d0)*ttol .AND. &
it .LE. itmax) Then
write(1,*)'Return mapping Ok!'
Else
write(1,*)'Numerical problem'
End If
Case (3) ! 3 : Eje entre f13_hs y f23_hs (Prs(1) = Prs(2) > Prs(3))
!STARTING POINT
!Cálculo valor de f13_hs, f23_hs y f_cap en Prs_E
qf1 = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs_E(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
f23_hs_fac1 = ((Prs_E(2)-Prs_E(3))/(qa-(Prs_E(2)-Prs_E(3))))
f23_hs_fac2 = 2d0*(Prs_E(2)-Prs_E(3))/Eur
f23_hs = (2d0*qa/Ei)*f23_hs_fac1-f23_hs_fac2-gamma_p_E
xlam_ast = p_ref/Eoed_ref
xkap_ast = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
alfa1 = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
alfa21 = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
alfa221 = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
alfa222 = (1d0-xK0_nc)*(1d0+xNu_ur)
alfa = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0
delta = (3d0+sin(phi))/(3d0-sin(phi))
q_m = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
p = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0
!Cálculo tolerancia f13_hs, f23_hs y f_cap en Prs_E
f13_hs_TOL = ttol*abs(f13_hs)+ttol
f23_hs_TOL = ttol*abs(f23_hs)+ttol
f_cap_TOL = ttol*abs(f_cap)+ttol
!Cálculo df13_hs/dS, df23_hs/dS y df_cap/dS en Prs_E
df13_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
df13_hsmdS(2) = 0d0
df13_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur
df23_hsmdS(1) = 0d0
df23_hsmdS(2) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
df23_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur
df_capmdS(1) = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(2) = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(3) = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
!Cálculo df13_hs/dq, df23_hs/dq y df_cap/dq en Prs_E (q1=gamma_p -para f13_hs y f23_hs-, q2=p_p)
df13_hsmdq = -1d0
df23_hsmdq = -1d0
df_capmdq = -2d0*p_p_E
!Cálculo h en Prs_E (dq=dlam1*h1 -para f13_hs y f23_hs- y dq2=dlam2*h2)
xh_hs = 1d0
xh_cap = 2d0*xH*p*(max(Prs_E(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
!Cálculo dg13_hs/dS, dg23_hs/dS y dg_cap/dS en Prs_E
sin_phi_m = (Prs_E(1)-Prs_E(3))/(Prs_E(1)+Prs_E(3)-2d0*c*cos(phi)/sin(phi))
sin_phi_cv = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
sin_psi_m = 0d0
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
psi_m = psi
sin_psi_m = sin(psi_m)
ElseIf (phi .EQ. 0d0) Then
sin_psi_m = 0d0
End If
dg13_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg13_hsmdS(2) = 0d0
dg13_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg23_hsmdS(1) = 0d0
dg23_hsmdS(2) = (1d0/2d0)*(1d0-sin_psi_m)
dg23_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg_capmdS(1) = df_capmdS(1)
dg_capmdS(2) = df_capmdS(2)
dg_capmdS(3) = df_capmdS(3)
!Cálculo d2g13_hs/dS2, d2g23_hs/dS2 y d2g_cap/dS2 en Prs_E
Call MZeroR(d2g13_hsmdS2, 9)
Call MZeroR(d2g23_hsmdS2, 9)
d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(2,1) = d2g_capmdS2(1,2)
d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(3,1) = d2g_capmdS2(1,3)
d2g_capmdS2(3,2) = d2g_capmdS2(2,3)
d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)
!Cálculo d2g13_hs/dSdq, d2g23_hs/dSdq y d2g_cap/dSdq en Prs_E
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(d2g23_hsmdSdq, 3)
Call MZeroR(d2g_capmdSdq, 3)
!Cálculo de los multiplicadores plásticos
Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)
Call DInProd(df13_hsmdS, D_dg13,3, df13_D_dg13)
xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs
Call MatMat(D,3,dg23_hsmdS,3,3,1,3,D_dg23,3)
Call DInProd(df13_hsmdS, D_dg23,3, df13_D_dg23)
xOmega12 = df13_D_dg23
Call MatMat(D,3,dg_capmdS,3,3,1,3,D_dgcap,3)
Call DInProd(df13_hsmdS, D_dgcap,3, df13_D_dgcap)
xOmega13 = df13_D_dgcap
Call DInProd(df23_hsmdS, D_dg13,3, df23_D_dg13)
xOmega21 = df23_D_dg13
Call DInProd(df23_hsmdS, D_dg23,3, df23_D_dg23)
xOmega22 = df23_D_dg23-df23_hsmdq*xh_hs
Call DInProd(df23_hsmdS, D_dgcap,3, df23_D_dgcap)
xOmega23 = df23_D_dgcap
Call DInProd(df_capmdS, D_dg13,3, dfcap_D_dg13)
xOmega31 = dfcap_D_dg13
Call DInProd(df_capmdS, D_dg23,3, dfcap_D_dg23)
xOmega32 = dfcap_D_dg23
Call DInProd(df_capmdS, D_dgcap,3, dfcap_D_dgcap)
xOmega33 = dfcap_D_dgcap-df_capmdq*xh_cap
!Resolución sistema y cálculo dxLam13_hs, dxLam23_hs y dxLam_cap
xMat(1,1) = xOmega11
xMat(1,2) = xOmega12
xMat(1,3) = xOmega13
xMat(2,1) = xOmega21
xMat(2,2) = xOmega22
xMat(2,3) = xOmega23
xMat(3,1) = xOmega31
xMat(3,2) = xOmega32
xMat(3,3) = xOmega33
f_b(1) = f13_hs
f_b(2) = f23_hs
f_b(3) = f_cap
Call MatInvPiv(xMat, xMat_inv, 3)
Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3)
dxLam13_hs = d_lam(1)
dxLam23_hs = d_lam(2)
dxLam_cap = d_lam(3)
xLam13_hs = dxLam13_hs
xLam23_hs = dxLam23_hs
xLam_cap = dxLam_cap
!Actualización variables
Call AddVec(D_dg13, D_dg23, -xLam13_hs, -xLam23_hs, 3, Prs_2sup)
Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)
Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
gamma_p = gamma_p_E + xLam13_hs*xh_hs
p_p = p_p_E + xLam_cap*xh_cap
!Cálculo valor de f13_hs, f23_hs y f_cap con la nueva Prs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p
f23_hs_fac1 = ((Prs(2)-Prs(3))/(qa-(Prs(2)-Prs(3))))
f23_hs_fac2 = 2d0*(Prs(2)-Prs(3))/Eur
f23_hs = (2d0*qa/Ei)*f23_hs_fac1-f23_hs_fac2-gamma_p
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0
!EULER BACKWARD
Do while ((abs(f13_hs) .GT. f13_hs_TOL .OR. &
abs(f23_hs) .GT. f23_hs_TOL .OR. &
abs(f_cap) .GT. f_cap_TOL) .AND. &
it .LE. itmax)
!Construcción matriz D
!Call Matriz3(Props, Prs, D, IsUndr)
!Cálculo df13_hs/dS, df23_hs/dS y df_cap/dS en Prs
df13_hsmdS(1) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)-2d0/Eur
df13_hsmdS(2) = 0d0
df13_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)+2d0/Eur
df23_hsmdS(1) = 0d0
df23_hsmdS(2) = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(2)))**2d0)-2d0/Eur
df23_hsmdS(3) = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)+2d0/Eur
df_capmdS(1) = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(2) = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
df_capmdS(3) = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
!Cálculo df13_hs/dq, df23_hs/dq y df_cap/dq en Prs (q1=gamma_p -para f13_hs y f23_hs-, q2=p_p)
df13_hsmdq = -1d0
df23_hsmdq = -1d0
df_capmdq = -2d0*p_p
!Cálculo h en Prs (dq=dlam1*h1 -para f13_hs y f23_hs- y dq2=dlam2*h2)
xh_hs = 1d0
xh_cap = 2d0*xH*p*(max(Prs(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
!Cálculo dg13_hs/dS, dg23_hs/dS y dg_cap/dS en Prs
sin_phi_m = (Prs(1)-Prs(3))/(Prs(1)+Prs(3)-2d0*c*cos(phi)/sin(phi))
sin_phi_cv = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
sin_psi_m = 0d0
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
psi_m = psi
sin_psi_m = sin(psi_m)
ElseIf (phi .EQ. 0d0) Then
sin_psi_m = 0d0
End If
dg13_hsmdS(1) = (1d0/2d0)*(1d0-sin_psi_m)
dg13_hsmdS(2) = 0d0
dg13_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg23_hsmdS(1) = 0d0
dg23_hsmdS(2) = (1d0/2d0)*(1d0-sin_psi_m)
dg23_hsmdS(3) = -(1d0/2d0)*(1d0+sin_psi_m)
dg_capmdS(1) = df_capmdS(1)
dg_capmdS(2) = df_capmdS(2)
dg_capmdS(3) = df_capmdS(3)
!Cálculo d2g13_hs/dS2, d2g23_hs/dS2 y d2g_cap/dS2 en Prs
Call MZeroR(d2g13_hsmdS2, 9)
Call MZeroR(d2g23_hsmdS2, 9)
d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(2,1) = d2g_capmdS2(1,2)
d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)
d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)
d2g_capmdS2(3,1) = d2g_capmdS2(1,3)
d2g_capmdS2(3,2) = d2g_capmdS2(2,3)
d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)
!Cálculo d2g13_hs/dSdq, d2g23_hs/dSdq y d2g_cap/dSdq en Prs
Call MZeroR(d2g13_hsmdSdq, 3)
Call MZeroR(d2g23_hsmdSdq, 3)
Call MZeroR(d2g_capmdSdq, 3)
!Cálculo del multiplicador plástico
Call MatMat(D, 3, d2g13_hsmdS2, 3, 3, 3, 3, D_dg13_2, 3)
Call MatMat(D, 3, d2g23_hsmdS2, 3, 3, 3, 3, D_dg23_2, 3)
Call MatMat(D, 3, d2g_capmdS2, 3, 3, 3, 3, D_dgcap_2, 3)
HH_inv = Unity + xLam13_hs*D_dg13_2 + xLam23_hs*D_dg23_2 + xLam_cap*D_dgcap_2
Call MatInvPiv(HH_inv, HH, 3)
Call AddVec(dg13_hsmdS, d2g13_hsmdSdq, 1d0, xLam13_hs*xh_hs, 3, theta13)
Call AddVec(dg23_hsmdS, d2g23_hsmdSdq, 1d0, xLam23_hs*xh_hs, 3, theta23)
Call AddVec(dg_capmdS, d2g_capmdSdq, 1d0, xLam_cap*xh_cap, 3, thetacap)
Call MatMat(HH, 3, df13_hsmdS, 3, 3, 1, 3, HH_df13,3)
Call MatMat(HH, 3, df23_hsmdS, 3, 3, 1, 3, HH_df23,3)
Call MatMat(HH, 3, df_capmdS, 3, 3, 1, 3, HH_dfcap,3)
Call MatMat(D, 3, theta13, 3, 3, 1, 3, D_theta13,3)
Call MatMat(D, 3, theta23, 3, 3, 1, 3, D_theta23,3)
Call MatMat(D, 3, thetacap, 3, 3, 1, 3, D_thetacap,3)
Call DInProd(D_theta13, HH_df13, 3, D_theta13_HH_df13)
Call DInProd(D_theta23, HH_df13, 3, D_theta23_HH_df13)
Call DInProd(D_thetacap, HH_df13, 3, D_thetacap_HH_df13)
Call DInProd(D_theta13, HH_df23, 3, D_theta13_HH_df23)
Call DInProd(D_theta23, HH_df23, 3, D_theta23_HH_df23)
Call DInProd(D_thetacap, HH_df23, 3, D_thetacap_HH_df23)
Call DInProd(D_theta13, HH_dfcap, 3, D_theta13_HH_dfcap)
Call DInProd(D_theta23, HH_dfcap, 3, D_theta23_HH_dfcap)
Call DInProd(D_thetacap, HH_dfcap, 3, D_thetacap_HH_dfcap)
xOmega11 = D_theta13_HH_df13-df13_hsmdq*xh_hs
xOmega12 = D_theta23_HH_df13
xOmega13 = D_thetacap_HH_df13
xOmega21 = D_theta13_HH_df23
xOmega22 = D_theta23_HH_df23-df23_hsmdq*xh_hs
xOmega23 = D_thetacap_HH_df23
xOmega31 = D_theta13_HH_dfcap
xOmega32 = D_theta23_HH_dfcap
xOmega33 = D_thetacap_HH_dfcap-df_capmdq*xh_cap
Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13,3)
Call MatMat(D, 3, dg23_hsmdS, 3, 3, 1, 3, D_dg23, 3)
Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)
Call AddVec(D_dg13, D_dg23, -xLam13_hs, -xLam23_hs, 3, Prs_res1)
Call AddVec(Prs_res1, D_dgcap, 1d0, -xLam_cap, 3, Prs_res2)
Call AddVec(Prs_E, Prs_res2, 1d0, 1d0, 3, Prs_res3)
Call AddVec(Prs, Prs_res3, 1d0, -1d0, 3, r)
Call DInProd(HH_df13, r, 3, r_HH_df13)
Call DInProd(HH_df23, r, 3, r_HH_df23)
Call DInProd(HH_dfcap, r, 3, r_HH_dfcap)
xOmegaf13 = f13_hs - r_HH_df13
xOmegaf23 = f23_hs - r_HH_df23
xOmegafcap = f_cap - r_HH_dfcap
!Resolución sistema y cálculo dxLam13_hs y dxLam23_hs
xMat(1,1) = xOmega11
xMat(1,2) = xOmega12
xMat(1,3) = xOmega13
xMat(2,1) = xOmega21
xMat(2,2) = xOmega22
xMat(2,3) = xOmega23
xMat(3,1) = xOmega31
xMat(3,2) = xOmega32
xMat(3,3) = xOmega33
f_b(1) = xOmegaf13
f_b(2) = xOmegaf23
f_b(3) = xOmegafcap
Call MatInvPiv(xMat, xMat_inv, 3)
Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3)
dxLam13_hs = d_lam(1)
dxLam23_hs = d_lam(2)
dxLam_cap = d_lam(3)
xLam13_hs = xLam13_hs + dxLam13_hs
xLam23_hs = xLam23_hs + dxLam23_hs
xLam_cap = xLam_cap + dxLam_cap
!Actualización variables
Call AddVec(D_dg13, D_dg23, -xLam13_hs, -xLam23_hs, 3, Prs_2sup)
Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)
Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
gamma_p = gamma_p_E + xLam13_hs*xh_hs
p_p = p_p_E + xLam_cap*xh_cap
!Cálculo valor de f13_hs, f23_hs y f_cap en (k+1) en Prs
qf1 = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
qf2 = (2d0*sin(phi)/(1d0-sin(phi)))
qf = qf1*qf2
qa = qf/R_f
E_fac1 = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
E_fac2 = (c*cos(phi)+p_ref*sin(phi))
E50 = E50_ref*(E_fac1/E_fac2)**xm
Ei = 2d0*E50/(2d0-R_f)
Eur = Eur_ref*(E_fac1/E_fac2)**xm
f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
f13_hs = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p
f23_hs_fac1 = ((Prs(2)-Prs(3))/(qa-(Prs(2)-Prs(3))))
f23_hs_fac2 = 2d0*(Prs(2)-Prs(3))/Eur
f23_hs = (2d0*qa/Ei)*f23_hs_fac1-f23_hs_fac2-gamma_p
q_m = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
p = (Prs(1)+Prs(2)+Prs(3))/3d0
f_cap = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0
Call WriIvl(ideb, 'it', it)
!Actualizo iteración
it = it+1
End do
!Actualizo variables de estado
StVar(1) = gamma_p
StVar(2) = min((-1d0)*p_p, (-1d0)*p_p_E)
Call WriVal(ideb, 'xLam13_hs', xLam13_hs)
Call WriVal(ideb, 'xLam23_hs', xLam23_hs)
Call WriVal(ideb, 'xLam_cap', xLam_cap)
If (xLam13_hs .GT. (-1d0)*ttol .AND. &
xLam23_hs .GT. (-1d0)*ttol .AND. &
xLam_cap .GT. (-1d0)*ttol .AND. &
it .LE. itmax) Then
write(1,*)'Return mapping Ok!'
Else
write(1,*)'Numerical problem'
End If
Case Default
write(1,*) ' incorrect iArea in HS_plast'
Stop
End Select
End ! HS_CAP_plast
____________________________________________________________________________________________________
CAP_PLAST.F90 , HS_PLAST.F90, HS_TRAC_PLAST.F90, MATRIZ6.f90, ETC... (OTHER SUBROUTINES)
_____________________________________________________________________________________________________
ETC...
ETC...
ETC...
_____________________________________________________________________________________________________
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
OK, using the link that you gave in the second paragraph of #5, I downloaded the Zip file http://kb.plaxis.nl/sites/kb.plaxis.nl/files/kb-downloads/udsm_examples_2012.zip. I built the 32-bit DLL using Intel Fortran 16 (32-bit) after fixing a couple of errors in the batch and options files that were in the Zip:
- In file ifcomp.bat, surround the replace > by ">", i.e., surround the > character with double quotes.
- In file comp_rel, change /optimize:4 to /O2
What they have done is to have the batch file first put the compiler-specific directive/declarations for the DLL entry point user_mod into a one-line file (impexp)using the echo command in the batch file. This line is then included in the source file usrmod.for. You do not have this INCLUDE statement in your source file in #5. Without it, the exported name will not be correctly seen when Plaxis attempts to load the DLL.
Everything that I have seen about Plaxis today tells me that they expect users to use one compiler to build the DLL. All that they care about is that the DLL routine be callable by Delphi calling protocol. Therefore, you can use the batch file ifcomp.bat to build the DLL (assuming that your copy of Plaxis wants to call a 32-bit DLL), and ignore the other batch files.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
As I'm using Plaxis in a 64 bit environtment, it's necessary to compile the corresponding 64-bit .dll file too (using ifcomp64.bat).
I've done what you pointed in the previous post, but Plaxis still not worrking (I mean Plaxis can't read all necessary information from the .dll's, because it doesn't let me to introduce model parameters in its input page).
In 64 bit dll compilation appeared next message:
">copy impivf impexp
System cannot find the file specified"
(however, this message has appeared always since I started compiling weeks ago)
Yes, I erased "include 'impexp'" statement, as you can see in my Usemod.f90 (copy in #5). What I do to compile using Intel is to activate the line: !DEC$ ATTRIBUTES DLLExport, StdCall, Reference :: User_Mod in Usermod.f90.
I think that until now, there were some necessary information for Plaxis missed in the dynamic libraries example_ivf.dll and example_ivf64.dll (due an error in compilation? maybe the error I pointed out before in 64bit dll compilation?).
And maybe Plaxis found this information in the dynamic library example_lah.dll (compiled with LF95).
I guess this because once add example_lah.dll to Plaxis folder, it work well (I mean that let me introduce the model parameters; this is the only thing that change adding example_lah.dll, but i'ts primordial). That's why I initially thoght that both compilers were necessary.
This supposed missed information in example_ivf.dll and example_ivf64.dll due to a possible error in compilation, might be related with user_add.f90 / useradddf.f90 , which are the responsible of parameter input in Plaxis. I've copy user_add.f90 and useradddf.f90 code (and others .f90) in my last post.
Now, I have the following orders in make files, but still not working:
IFCOMP.BAT
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod">"impexp
ifort @comp_rel UsrMod.f90 UsrAddDF.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf
IMPEXP
: User_Mod
COMP_REL
/real_size:64
/fpconstant
/O2
/static
/assume:noaccuracy_sensitive
/assume:byterecl
/assume:buffered_io
/notraceback
/debug:none
/check:none
/heap-arrays:1
/fpp
/recursive
IFCOMP64.BAT
REM compile 64-bit DLL: make sure you're in a 64-bit development environment
copy impivf impexp
ifort /static UsrMod.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf64
And, as I said before, there must be some error/missed orders in the previous files (IFCOMP.BAT, IMPEXP, COMP_REL, IFCOMP64.BAT) that is contained in the next one:
LFCOMP.BAT
echo DLL_Export User_Mod
lf95 -c usrlib.f90
lib /out:lfusrlib.lib usrlib.obj
lf95 -c -ml winapi okmsglah.f90
lf95 -win -dll UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll
Thank you so much for your support.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Page 175 of http://www.plaxis.nl/files/files/2D-3-Material-Models.pdf says:
In order to compile as 64-bit, you need both a 32-bit compiled 'USRMOD.DLL' and a
64-bit compiled ' USRMOD64.DLL' file in the usdm folder under the PLAXIS program
directory. The last one only needs to contain the 'User_Mod' subroutine.
I take this to imply that you should produce two versions of your DLL, one 32-bit and the other 64-bit. For the example code provided in http://kb.plaxis.nl/sites/kb.plaxis.nl/files/kb-downloads/udsm_examples_2012.zip, all that I had to do was
Open a 32-bit IFort command window, create a new directory and extract the examples zip file and correct comp_rel as in #6. Delete the extracted DLLs or put them into a different directory. Issue the commands
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod > impexp ifort @comp_rel UsrMod.for UsrAddDF.for usrlib.for okmsgivf.for /dll /exe:example_ivf
Open a 64-bit IFort command window, change to the same working directory as in the previous step, and issue the commands
echo !DEC$ ATTRIBUTES DLLExport :: User_Mod > impexp ifort @comp_rel UsrMod.for usrlib.for okmsgivf.for /dll /exe:example_ivf64
Note that I did not use the troublesome batch files with their redirected echo commands, nor did I use the one-line file impivf that is in the Plaxis examples zip but you do not have (or have erased).
The two DLLs were produced. They are not the same as the DLLs in the Plaxis examples zip, probably because they used an older version of the compiler.
To avoid confusion, I suggest that you re-do what I have described in #8 in a separate directory into which you have extracted the Zip file from Plaxis. Once this works on your system, you can move to the directory with your own sources and repeat the procedure but with your sources instead.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I did what you post on #8 with Plaxis Zip (erasing the line with " include 'impexp' " in usrmod.for, changing optimize:4 to O2 in comp_rel, and issuing the comands you've written in #8 in my IA-32 and IA-64 command windows)
Dll's are created, but Plaxis still not working, because it doesn't let me introduce model parameters once opened Plaxis Console and selected "Used-defined model" in Material Model type.
As it's explained in the manual (pp. 172-174), Plaxis uses the procedures in UsrAddDF.for to activate user-interface functions (this is need for model parameters input and its recognition). I think that something is going wrong related with this.
Each procedure has it's own directive (I guess all is correct ¿?):
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetModelCount
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetModelName
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetParamCount
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetParamName
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: GetParamUnit
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarCount
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarName
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarUnit
I've tried also some variations of what you've told me to do:
- Activate the following line in UsrMod.for !DEC$ ATTRIBUTES DLLExport, StdCall, Reference :: User_Mod . Try compilation.
- Activate the following line in UsrMod.for: include 'usr_add.for' ! give parameter names and units and erase UsrAddDF.for from issue commands. Try compilation.
But none of them fix the problem.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Javier C. wrote:
I did what you post on #8 with Plaxis Zip (erasing the line with " include 'impexp' " in usrmod.for, changing optimize:4 to O2 in comp_rel, and issuing the comands you've written in #8 in my IA-32 and IA-64 command windows)
You should not do that. The include 'impexp' line is essential for making the Plaxis-DLL interface work. What gets included (just one line) has to be different for 32-bit and 64-bit, because the calling conventions are different.
If you prefer, you could make two copies of usrmod.for, say usrmod32.for and usrmod64.for, in which you replace the include 'impexp' line by
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
in usrmod32.for and by
!DEC$ ATTRIBUTES DLLExport :: User_Mod
in usrmod64.for. When you build the DLLs, use usrmod32.for for the 32-bit DLL, and usrmod64.for for the 64-bit DLL.
Beyond that, I cannot comment about the !DEC$ lines that you showed in #9 for the other routines because I do not have Plaxis and I do not know what its requirements are. The documentation seems to indicate that those routines are not needed in the 64-bit DLL, when it says "The last one (USRMOD64.DLL) only needs to contain the 'User_Mod' subroutine.".
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Oh, ok, I'm Sorry.
Now I have two copies of UsrMod (UsrMod32.for and UsrMod64.for) with its corresponding lines (post #10) containing !DEC$ directives in substitution of include 'impexp' (I understand now which was the function of line include 'ìmpexp' UsrMod.for, thank you).
Then, I write issues commands:
ifort @comp_rel UsrMod32.for UsrAddDF.for usrlib.for okmsgivf.for /dll /exe:example_ivf
for UsrMod32.for for 32-bit DLL and
ifort @comp_rel UsrMod64.for usrlib.for okmsgivf.for /dll /exe:example_ivf64
for UsrMod64.for for 64-bit DLL
Then two .dll are created (example_ivf.dll and example_ivf64.dll). I have copy them in Plaxis Folder.
After that, I open Plaxis, but when I'm going to introduce Model Input Parameters, Plaxis doesn't shown anything.
If all was correct, in the table of Input Parameters of the Model (while executing Plaxis) it should appear the name of the .dll that contains this information, the name of the model, the name of the parameters, their units, etc. All this information is provided by UsrAddDF.for.
UsrAddDF.for contains the code and directives that let Plaxis "read" how many parameters have the model/s that we are defining, what are their names and units, etc. (I don't know much more about UsrAddDF.for; in my code, I just adapted UsrAddDF from the one in Plaxis ZIP).
How I said before, I solved this using LF95 (now I can't, since License has expired). To compile with LF95, I only changed nex line
!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
in UsrMod by
DLL_Export User_Mod
Then, after execute lfcomp.bat with LF95, the file example_lah.dll was created and once added together with the files example_ivf.dll and example_ivf64.dll (compiled with Intel) in Plaxis Folder, all works well (Plaxis could read parameters, etc.)
Maybe some clue is in lfcomp.bat (used in LF95 compilation).
Thank you for all your time and comments, and sorry for my ignorance in this issue.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I think I was wrong in one thing. What made Plaxis run until now when using the option of User Soil Defined Model was example_lah.dll compiled with LF95. I've tried to put and old version of this dll alone in Plaxis Folder (without example_ivf.dll and example_ivf64.dll), and it works well.
In any case, I'm sure that two dll compiled with Intel will finally works.
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page