!***********************************************************************!
! TITLE: CFVR3901 - 4 node 2D laplace element test
! SUBTITLE: Test 4 node 2D Laplace elements (PLANE55)
!
! DESCRIPTION:  This example test the thermal properties of a 4 node 2D PLANE55
! DESCRIPTION:  model.
! DESCRIPTION:
! DESCRIPTION:  The model consists on an irregular hepthagon, same as the one on
! DESCRIPTION:  the following picture.
! DESCRIPTION:
! DESCRIPTION:  The example is launched with ANSYS and with CivilFEM separately;
! DESCRIPTION:  it compares ANSYS results with the ~LPSOLVE command.
! DESCRIPTION:
! DESCRIPTION:  ANSYS uses a thermal laplace ecuation model, and CivilFEM solves
! DESCRIPTION:  it as a Seepage analysis.
!
! ELEMENT TYPE: PLANE55
! MODULES:
! UNITS: SI
! KEYWORD1: Material Models
! KEYWORD2: Materials
!
!***********************************************************************!
  FINISH
  ~CFCLEAR,,1
  NomFile='CFVR3901'
  /TITLE, %NomFile%, 4 node 2D Laplace element test

!--------------------------------------------------------------------------------------
! MODEL WITH ANSYS (thermal)
!--------------------------------------------------------------------------------------
/prep7
et,1,plane55

mp,kxx,1,10
mp,kyy,1,20
mp,kxx,2,2
mp,kyy,2,6

! Area 1

x_off = 5
y_off = 5

L1 = x_off+5
L2 = L1+3
L3 = L2+2
H1 = y_off+4
H2 = H1+4

k,,x_off,y_off
k,,L1,y_off
k,,L2,H1
k,,L3,H1
k,,L3,H2
k,,x_off,H2

l,1,2
*repeat,5,1,1
l,6,1

al,all

! Area 2

k,,x_off+2.5,y_off+1
k,,x_off+4.5,y_off+4
k,,x_off+3.5,y_off+5
k,,x_off+1.5,y_off+4

l,7,8
l,8,9
l,9,10
l,10,7

al,7,8,9,10

asba,1,2,,,keep

H3 = y_off+12
H4 = y_off+20
L4 = x_off+15
L5 = x_off+20

k,,x_off,H3
k,,L4,H3
k,,L4,H1

k,,x_off,H4
k,,L5,H4
k,,L5,H3

A,6,11,12,5
A,5,4,13,12
A,11,14,15,16,12,11
A,13,16,12,13

lesize,13,,,5,0.5
lesize,4,,,5,0.5

! Area 2 con cuadr
esize,0.8
mat,1
amesh,2
mshape,1
esize,1
mat,1
amesh,3
mat,2
amesh,1
mat,1
amesh,4
amesh,6
esize,2
mat,2
amesh,5
/solu

dl,16,,temp,120
dl,20,,temp,30
dl,6,,temp,50
dl,2,,temp,80

SOLVE
/post1
SET,1

! Results

NNUM  = NDINQR(0,12)
ELNUM = ELMIQR(0,12)
*DIM,H_ANSYS,,NNUM
*DIM,TGX_ANSYS,,NNUM
*DIM,TGY_ANSYS,,NNUM
*DIM,TFX_ANSYS,,NNUM
*DIM,TFY_ANSYS,,NNUM
*VGET,H_ANSYS,NODE,1,TEMP
*VGET,TFX_ANSYS,NODE,1,TF,X
*VGET,TFY_ANSYS,NODE,1,TF,Y
*VGET,TGX_ANSYS,NODE,1,TG,X
*VGET,TGY_ANSYS,NODE,1,TG,Y

PARSAV,ALL,E4_ANSYS


!--------------------------------------------------------------------------------------
! MODEL WITH CIVILFEM
!--------------------------------------------------------------------------------------

finish
~cfclear,,1

~CFACTIV,CIVILFEM,Y
~CFACTIV,GETC,Y

FINISH
~SEEPAGE,2D

/prep7

~cfmp,1,lib,soil,,sc
~cfmp,1,soil,kx,,10
~cfmp,1,soil,ky,,20
~cfmp,2,lib,soil,,sc
~cfmp,2,soil,kx,,2
~cfmp,2,soil,ky,,6

! Area 1

x_off = 5
y_off = 5

L1 = x_off+5
L2 = L1+3
L3 = L2+2
H1 = y_off+4
H2 = H1+4

k,,x_off,y_off
k,,L1,y_off
k,,L2,H1
k,,L3,H1
k,,L3,H2
k,,x_off,H2

l,1,2
*repeat,5,1,1
l,6,1

al,all

! Area 2

k,,x_off+2.5,y_off+1
k,,x_off+4.5,y_off+4
k,,x_off+3.5,y_off+5
k,,x_off+1.5,y_off+4

l,7,8
l,8,9
l,9,10
l,10,7

al,7,8,9,10

asba,1,2,,,keep

H3 = y_off+12
H4 = y_off+20
L4 = x_off+15
L5 = x_off+20

k,,x_off,H3
k,,L4,H3
k,,L4,H1

k,,x_off,H4
k,,L5,H4
k,,L5,H3

A,6,11,12,5
A,5,4,13,12
A,11,14,15,16,12,11
A,13,16,12,13

lesize,13,,,5,0.5
lesize,4,,,5,0.5

! Area 2 con cuadr
esize,0.8
mat,1
amesh,2
mshape,1
esize,1
mat,1
amesh,3
mat,2
amesh,1
mat,1
amesh,4
amesh,6
esize,2
mat,2
amesh,5
/solu

~DLHEAD,ADD,16,120
~DLHEAD,ADD,20,30
~DLHEAD,ADD, 6,50
~DLHEAD,ADD, 2,80

~LPSOLVE
/post1
file,file,seep
SET,1

PARRES,new,E4_ANSYS

*DIM,H_CF,,NNUM
*DIM,TGX_CF,,NNUM
*DIM,TGY_CF,,NNUM
*DIM,TFX_CF,,NNUM
*DIM,TFY_CF,,NNUM
*VGET,H_CF,NODE,1,HEAD
*VGET,TFX_CF,NODE,1,S,X
*VGET,TFY_CF,NODE,1,S,Y
*VGET,TGX_CF,NODE,1,EPEL,X
*VGET,TGY_CF,NODE,1,EPEL,Y
! Tolerance
TOL = 1e-4
*DIM,VTOL,,NNUM
*VFILL,VTOL,RAMP,TOL

*DIM,H_COMP,,NNUM
*DIM,H_OK,,NNUM
*DIM,TGX_COMP,,NNUM
*DIM,TGX_OK,,NNUM
*DIM,TGY_COMP,,NNUM
*DIM,TGY_OK,,NNUM
*DIM,TFX_COMP,,NNUM
*DIM,TFX_OK,,NNUM
*DIM,TFY_COMP,,NNUM
*DIM,TFY_OK,,NNUM

! a-b<Tolerance

*VOPER, H_COMP(1),H_ANSYS(1),SUB,H_CF(1)
*VOPER, H_OK(1),H_COMP(1),LE,VTOL(1)
*VSCFUN,H_BOOL,SUM,H_OK ! HBOOL debe ser = NNUM

*VOPER, TGX_COMP(1),TGX_ANSYS(1),SUB,TGX_CF(1)
*VOPER, TGX_OK(1),TGX_COMP(1),LE,VTOL(1)
*VSCFUN,TGX_BOOL,SUM,TGX_OK

*VOPER, TGY_COMP(1),TGY_ANSYS(1),SUB,TGY_CF(1)
*VOPER, TGY_OK(1),TGY_COMP(1),LE,VTOL(1)
*VSCFUN,TGY_BOOL,SUM,TGY_OK

*VOPER, TFX_COMP(1),TFX_ANSYS(1),SUB,TFX_CF(1)
*VOPER, TFX_OK(1),TFX_COMP(1),LE,VTOL(1)
*VSCFUN,TFX_BOOL,SUM,TFX_OK

*VOPER, TFY_COMP(1),TFY_ANSYS(1),SUB,TFY_CF(1)
*VOPER, TFY_OK(1),TFY_COMP(1),LE,VTOL(1)
*VSCFUN,TFY_BOOL,SUM,TFY_OK
/IMAGE,SAVE,%NomFile%,BMP	!Save Image to File
!--------------------------------------------------------------------------------------
! DATA CHECK
!--------------------------------------------------------------------------------------
! Data comparison number
  NComp = 5
  NComp_ch = 0

! Marix dim.
  *DIM,LABEL,CHAR,Ncomp,1
  *DIM,LABEL_CH,CHAR,Ncomp_ch,1
  *DIM,VALUE,,Ncomp,3
  *DIM,VALUE_CH,CHAR,Ncomp_ch,3
  *DIM,TOLER,,Ncomp,2

! Labels
  LABEL(1,1) = 'H'
  LABEL(2,1) = 'TGX'
  LABEL(3,1) = 'TGY'
  LABEL(4,1) = 'TFX'
  LABEL(5,1) = 'TFY'

! Correct values
  VALUE(1,1)=NNUM
  VALUE(2,1)=NNUM
  VALUE(3,1)=NNUM
  VALUE(4,1)=NNUM
  VALUE(5,1)=NNUM

! Obtained values
  VALUE(1,2)= H_BOOL
  VALUE(2,2)= TGX_BOOL
  VALUE(3,2)= TGY_BOOL
  VALUE(4,2)= TFX_BOOL
  VALUE(5,2)= TFY_BOOL

! Warning and error tolerances
  *DO,II,1,5
    TOLER(II, 1)= 1E-9
    TOLER(II, 2)= 1E-9
  *ENDDO

!--------------------------------------------------------------------------------------
! Results comparison
!--------------------------------------------------------------------------------------
  COMPARA.MAC