1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> \brief This module contains fstr control file data obtaining functions
6
7module fstr_ctrl_common
8  use m_fstr
9  use hecmw
10  use mContact
11  use m_timepoint
12
13  implicit none
14
15  include 'fstr_ctrl_util_f.inc'
16
17  private :: pc_strupr
18
19contains
20
21  subroutine pc_strupr( s )
22    implicit none
23    character(*) :: s
24    integer :: i, n, a, da
25
26    n = len_trim(s)
27    da = iachar('a') - iachar('A')
28    do i = 1, n
29      a = iachar(s(i:i))
30      if( a > iachar('Z')) then
31        a = a - da
32        s(i:i) = achar(a)
33      end if
34    end do
35  end subroutine pc_strupr
36
37
38  !> Read in !SOLUTION
39  function fstr_ctrl_get_SOLUTION( ctrl, type, nlgeom )
40    integer(kind=kint) :: ctrl
41    integer(kind=kint) :: type
42    logical            :: nlgeom
43    integer(kind=kint) :: fstr_ctrl_get_SOLUTION
44
45    integer(kind=kint) :: ipt
46    character(len=80) :: s
47
48    fstr_ctrl_get_SOLUTION = -1
49
50    s = 'ELEMCHECK,STATIC,EIGEN,HEAT,DYNAMIC,NLSTATIC,STATICEIGEN,NZPROF '
51    if( fstr_ctrl_get_param_ex( ctrl,      'TYPE ',     s,    1,   'P',  type )/= 0) return
52    type = type -1
53
54    ipt=0
55    if( fstr_ctrl_get_param_ex( ctrl,    'NONLINEAR ',  '# ',    0,   'E',   ipt )/= 0) return
56    if( ipt/=0 .and. ( type == kstSTATIC .or. type == kstDYNAMIC )) nlgeom = .true.
57
58    if( type == 5 ) then !if type == NLSTATIC
59      type = kstSTATIC
60      nlgeom = .true.
61    end if
62    if( type == kstSTATICEIGEN ) nlgeom = .true.
63
64    fstr_ctrl_get_SOLUTION = 0
65  end function fstr_ctrl_get_SOLUTION
66
67
68  !> Read in !SOLVER
69  function fstr_ctrl_get_SOLVER( ctrl, method, precond, nset, iterlog, timelog, steplog, nier, &
70      iterpremax, nrest, scaling, &
71      dumptype, dumpexit, usejad, ncolor_in, mpc_method, estcond, method2, recyclepre, &
72      solver_opt1, solver_opt2, solver_opt3, solver_opt4, solver_opt5, solver_opt6, &
73      resid, singma_diag, sigma, thresh, filter )
74    integer(kind=kint) :: ctrl
75    integer(kind=kint) :: method
76    integer(kind=kint) :: precond
77    integer(kind=kint) :: nset
78    integer(kind=kint) :: iterlog
79    integer(kind=kint) :: timelog
80    integer(kind=kint) :: steplog
81    integer(kind=kint) :: nier
82    integer(kind=kint) :: iterpremax
83    integer(kind=kint) :: nrest
84    integer(kind=kint) :: scaling
85    integer(kind=kint) :: dumptype
86    integer(kind=kint) :: dumpexit
87    integer(kind=kint) :: usejad
88    integer(kind=kint) :: ncolor_in
89    integer(kind=kint) :: mpc_method
90    integer(kind=kint) :: estcond
91    integer(kind=kint) :: method2
92    integer(kind=kint) :: recyclepre
93    integer(kind=kint) :: solver_opt1
94    integer(kind=kint) :: solver_opt2
95    integer(kind=kint) :: solver_opt3
96    integer(kind=kint) :: solver_opt4
97    integer(kind=kint) :: solver_opt5
98    integer(kind=kint) :: solver_opt6
99    real(kind=kreal) :: resid
100    real(kind=kreal) :: singma_diag
101    real(kind=kreal) :: sigma
102    real(kind=kreal) :: thresh
103    real(kind=kreal) :: filter
104    integer(kind=kint) :: fstr_ctrl_get_SOLVER
105
106    character(92) :: mlist = '1,2,3,4,101,CG,BiCGSTAB,GMRES,GPBiCG,DIRECT,DIRECTmkl,DIRECTlag,MUMPS,MKL '
107    character(24) :: dlist = '0,1,2,3,NONE,MM,CSR,BSR '
108
109    integer(kind=kint) :: number_number = 5
110    integer(kind=kint) :: indirect_number = 4
111    integer(kind=kint) :: iter, time, sclg, dmpt, dmpx, usjd, step
112
113    fstr_ctrl_get_SOLVER = -1
114
115    iter = iterlog+1
116    time = timelog+1
117    step = steplog+1
118    sclg = scaling+1
119    dmpt = dumptype+1
120    dmpx = dumpexit+1
121    usjd = usejad+1
122    !* parameter in header line -----------------------------------------------------------------*!
123
124    ! JP-0
125    if( fstr_ctrl_get_param_ex( ctrl, 'METHOD ',   mlist,              1,   'P',   method  ) /= 0) return
126    if( fstr_ctrl_get_param_ex( ctrl, 'PRECOND ', '1,2,3,4,5,6,7,8,9,10,11,12,20,21,30,31,32 ' ,0, 'I', precond ) /= 0) return
127    if( fstr_ctrl_get_param_ex( ctrl, 'NSET ',    '0,-1,+1 ',          0,   'I',   nset    ) /= 0) return
128    if( fstr_ctrl_get_param_ex( ctrl, 'ITERLOG ', 'NO,YES ',           0,   'P',   iter ) /= 0) return
129    if( fstr_ctrl_get_param_ex( ctrl, 'TIMELOG ', 'NO,YES,VERBOSE ',   0,   'P',   time ) /= 0) return
130    if( fstr_ctrl_get_param_ex( ctrl, 'STEPLOG ', 'NO,YES ',           0,   'P',   step ) /= 0) return
131    if( fstr_ctrl_get_param_ex( ctrl, 'SCALING ', 'NO,YES ',           0,   'P',   sclg ) /= 0) return
132    if( fstr_ctrl_get_param_ex( ctrl, 'DUMPTYPE ', dlist,              0,   'P',   dmpt ) /= 0) return
133    if( fstr_ctrl_get_param_ex( ctrl, 'DUMPEXIT ','NO,YES ',           0,   'P',   dmpx ) /= 0) return
134    if( fstr_ctrl_get_param_ex( ctrl, 'USEJAD '  ,'NO,YES ',           0,   'P',   usjd ) /= 0) return
135    if( fstr_ctrl_get_param_ex( ctrl, 'MPCMETHOD ','# ',               0, 'I',mpc_method) /= 0) return
136    if( fstr_ctrl_get_param_ex( ctrl, 'ESTCOND '  ,'# ',               0,   'I',estcond ) /= 0) return
137    if( fstr_ctrl_get_param_ex( ctrl, 'METHOD2 ',  mlist,              0,   'P',   method2 ) /= 0) return
138    ! JP-1
139    if( method > number_number ) then  ! JP-2
140      method = method - number_number
141      if( method > indirect_number ) then
142        ! JP-3
143        method = method - indirect_number + 100
144        if( method == 103 ) method = 101 ! DIRECTlag => DIRECT
145        if( method == 105 ) method = 102 ! MKL => DIRECTmkl
146      end if
147    end if
148    if( method2 > number_number ) then  ! JP-2
149      method2 = method2 - number_number
150      if( method2 > indirect_number ) then
151        ! JP-3
152        method2 = method2 - indirect_number + 100
153      end if
154    end if
155
156    dumptype = dmpt - 1
157    if( dumptype >= 4 ) then
158      dumptype = dumptype - 4
159    end if
160
161    !* data --------------------------------------------------------------------------------------- *!
162    ! JP-4
163    if( fstr_ctrl_get_data_ex( ctrl, 1,   'iiiii ', nier, iterpremax, nrest, ncolor_in, recyclepre )/= 0) return
164    if( fstr_ctrl_get_data_ex( ctrl, 2,   'rrr ', resid, singma_diag, sigma )/= 0) return
165
166    if( precond == 20 .or. precond == 21) then
167      if( fstr_ctrl_get_data_ex( ctrl, 3, 'rr ', thresh, filter)/= 0) return
168    else if( precond == 5 ) then
169      if( fstr_ctrl_get_data_ex( ctrl, 3, 'iiiiii ', solver_opt1, solver_opt2, solver_opt3, &
170        & solver_opt4, solver_opt5, solver_opt6 )/= 0) return
171    else if( method == 101 ) then
172      if( fstr_ctrl_get_data_ex( ctrl, 3, 'i ', solver_opt1 )/= 0) return
173    end if
174
175    iterlog = iter -1
176    timelog = time -1
177    steplog = step -1
178    scaling = sclg -1
179    dumpexit = dmpx -1
180    usejad = usjd -1
181
182    fstr_ctrl_get_SOLVER = 0
183
184  end function fstr_ctrl_get_SOLVER
185
186
187  !> Read in !STEP
188  function fstr_ctrl_get_STEP( ctrl, amp, iproc )
189    integer(kind=kint) :: ctrl
190    character(len=HECMW_NAME_LEN) :: amp
191    integer(kind=kint) :: iproc
192    integer(kind=kint) :: fstr_ctrl_get_STEP
193
194    integer(kind=kint) :: ipt = 0
195    integer(kind=kint) :: ip = 0
196
197    fstr_ctrl_get_STEP = -1
198
199    if( fstr_ctrl_get_param_ex( ctrl, 'AMP ',     '# ',  0, 'S', amp )/= 0) return
200    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ',   'STANDARD,NLGEOM ', 0, 'P',   ipt    )/= 0) return
201    if( fstr_ctrl_get_param_ex( ctrl, 'NLGEOM ',  '# ',           0,    'E',   ip     )/= 0) return
202
203    if( ipt == 2 .or. ip == 1 ) iproc = 1
204
205    fstr_ctrl_get_STEP = 0
206
207  end function fstr_ctrl_get_STEP
208
209  !> Read in !STEP and !ISTEP
210  logical function fstr_ctrl_get_ISTEP( ctrl, hecMESH, steps, tpname, apname )
211    use fstr_setup_util
212    use m_step
213    integer(kind=kint), intent(in)        :: ctrl      !< ctrl file
214    type (hecmwST_local_mesh), intent(in) :: hecMESH   !< mesh information
215    type(step_info), intent(out)          :: steps     !< step control info
216    character(len=*), intent(out)         :: tpname    !< name of timepoints
217    character(len=*), intent(out)         :: apname    !< name of auto increment parameter
218
219    character(len=HECMW_NAME_LEN) :: data_fmt,ss, data_fmt1
220    character(len=HECMW_NAME_LEN) :: amp
221    character(len=HECMW_NAME_LEN) :: header_name
222    integer(kind=kint) :: bcid
223    integer(kind=kint) :: i, n, sn, ierr
224    integer(kind=kint) :: bc_n, load_n, contact_n
225    real(kind=kreal) :: fn, f1, f2, f3
226
227    fstr_ctrl_get_ISTEP = .false.
228
229    write(ss,*)  HECMW_NAME_LEN
230    write( data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)), 'I '
231    write( data_fmt1, '(a,a,a)') 'S', trim(adjustl(ss)),'rrr '
232
233    call init_stepInfo(steps)
234    steps%solution = stepStatic
235    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ',   'STATIC,VISCO ', 0, 'P', steps%solution )/= 0) return
236    steps%inc_type = stepFixedInc
237    if( fstr_ctrl_get_param_ex( ctrl, 'INC_TYPE ', 'FIXED,AUTO ', 0, 'P', steps%inc_type )/= 0) return
238    if( fstr_ctrl_get_param_ex( ctrl, 'SUBSTEPS ',  '# ',  0, 'I', steps%num_substep )/= 0) return
239    steps%initdt = 1.d0/steps%num_substep
240    if( fstr_ctrl_get_param_ex( ctrl, 'ITMAX ',  '# ',  0, 'I', steps%max_iter )/= 0) return
241    if( fstr_ctrl_get_param_ex( ctrl, 'MAXITER ',  '# ',  0, 'I', steps%max_iter )/= 0) return
242    if( fstr_ctrl_get_param_ex( ctrl, 'MAXCONTITER ',  '# ',  0, 'I', steps%max_contiter )/= 0) return
243    if( fstr_ctrl_get_param_ex( ctrl, 'CONVERG ',  '# ',  0, 'R', steps%converg )/= 0) return
244    if( fstr_ctrl_get_param_ex( ctrl, 'MAXRES ',  '# ',  0, 'R', steps%maxres )/= 0) return
245    amp = ""
246    if( fstr_ctrl_get_param_ex( ctrl, 'AMP ',  '# ',  0, 'S', amp )/= 0) return
247    if( len( trim(amp) )>0 ) then
248      call amp_name_to_id( hecMESH, '!STEP', amp, steps%amp_id )
249    endif
250    tpname=""
251    if( fstr_ctrl_get_param_ex( ctrl, 'TIMEPOINTS ',  '# ',  0, 'S', tpname )/= 0) return
252    apname=""
253    if( fstr_ctrl_get_param_ex( ctrl, 'AUTOINCPARAM ',  '# ',  0, 'S', apname )/= 0) return
254
255    n = fstr_ctrl_get_data_line_n( ctrl )
256    if( n == 0 ) then
257      fstr_ctrl_get_ISTEP = .true.;  return
258    endif
259
260    f2 = steps%mindt
261    f3 = steps%maxdt
262    if( fstr_ctrl_get_data_ex( ctrl, 1, data_fmt1, ss, f1, f2, f3  )/= 0) return
263    read( ss, * , iostat=ierr ) fn
264    sn=1
265    if( ierr==0 ) then
266      steps%initdt = fn
267      steps%elapsetime = f1
268      if( steps%inc_type == stepAutoInc ) then
269        steps%mindt = min(f2,steps%initdt)
270        steps%maxdt = f3
271      endif
272      steps%num_substep = max(int((f1+0.999999999d0*fn)/fn),steps%num_substep)
273      !if( mod(f1,fn)/=0 ) steps%num_substep =steps%num_substep+1
274      sn = 2
275    endif
276
277    bc_n = 0
278    load_n = 0
279    contact_n = 0
280    do i=sn,n
281      if( fstr_ctrl_get_data_ex( ctrl, i, data_fmt, header_name, bcid  )/= 0) return
282      if( trim(header_name) == 'BOUNDARY' ) then
283        bc_n = bc_n + 1
284      else if( trim(header_name) == 'LOAD' ) then
285        load_n = load_n +1
286      else if( trim(header_name) == 'CONTACT' ) then
287        contact_n = contact_n+1
288      else if( trim(header_name) == 'TEMPERATURE' ) then
289        !   steps%Temperature = .true.
290      endif
291    end do
292
293    if( bc_n>0 ) allocate( steps%Boundary(bc_n) )
294    if( load_n>0 ) allocate( steps%Load(load_n) )
295    if( contact_n>0 ) allocate( steps%Contact(contact_n) )
296
297    bc_n = 0
298    load_n = 0
299    contact_n = 0
300    do i=sn,n
301      if( fstr_ctrl_get_data_ex( ctrl, i, data_fmt, header_name, bcid  )/= 0) return
302      if( trim(header_name) == 'BOUNDARY' ) then
303        bc_n = bc_n + 1
304        steps%Boundary(bc_n) = bcid
305      else if( trim(header_name) == 'LOAD' ) then
306        load_n = load_n +1
307        steps%Load(load_n) = bcid
308      else if( trim(header_name) == 'CONTACT' ) then
309        contact_n = contact_n+1
310        steps%Contact(contact_n) = bcid
311      endif
312    end do
313
314    fstr_ctrl_get_ISTEP = .true.
315  end function fstr_ctrl_get_ISTEP
316
317  !> Read in !SECTION
318  integer function fstr_ctrl_get_SECTION( ctrl, hecMESH, sections )
319    integer(kind=kint), intent(in)           :: ctrl
320    type (hecmwST_local_mesh), intent(inout) :: hecMESH   !< mesh information
321    type (tSection), pointer, intent(inout)  :: sections(:)
322
323    integer(kind=kint)            :: j, k, sect_id, ori_id, elemopt
324    integer(kind=kint),save       :: cache = 1
325    character(len=HECMW_NAME_LEN) :: sect_orien
326    character(16) :: form361list = 'FI,BBAR,IC,FBAR '
327
328    fstr_ctrl_get_SECTION = -1
329
330    if( fstr_ctrl_get_param_ex( ctrl, 'SECNUM ',  '# ',  1, 'I', sect_id )/= 0) return
331    if( sect_id > hecMESH%section%n_sect ) return
332
333    elemopt = 0
334    if( fstr_ctrl_get_param_ex( ctrl, 'FORM361 ',   form361list, 0, 'P', elemopt )/= 0) return
335    if( elemopt > 0 ) sections(sect_id)%elemopt361 = elemopt
336
337    ! sectional orientation ID
338    hecMESH%section%sect_orien_ID(sect_id) = -1
339    if( fstr_ctrl_get_param_ex( ctrl, 'ORIENTATION ',  '# ',  0, 'S', sect_orien )/= 0) return
340
341    if( associated(g_LocalCoordSys) ) then
342      k = size(g_LocalCoordSys)
343
344      if(cache < k)then
345        if( sect_orien == g_LocalCoordSys(cache)%sys_name ) then
346          hecMESH%section%sect_orien_ID(sect_id) = cache
347          cache = cache + 1
348          fstr_ctrl_get_SECTION = 0
349          return
350        endif
351      endif
352
353      do j=1, k
354        if( sect_orien == g_LocalCoordSys(j)%sys_name ) then
355          hecMESH%section%sect_orien_ID(sect_id) = j
356          cache = j + 1
357          exit
358        endif
359      enddo
360    endif
361
362    fstr_ctrl_get_SECTION = 0
363
364  end function fstr_ctrl_get_SECTION
365
366
367  !> Read in !WRITE
368  function fstr_ctrl_get_WRITE( ctrl, res, visual, femap )
369    integer(kind=kint) :: ctrl
370    integer(kind=kint) :: res
371    integer(kind=kint) :: visual
372    integer(kind=kint) :: femap
373    integer(kind=kint) :: fstr_ctrl_get_WRITE
374
375    fstr_ctrl_get_WRITE = -1
376
377    ! JP-6
378    if( fstr_ctrl_get_param_ex( ctrl, 'RESULT ',  '# ',    0,   'E',   res    )/= 0) return
379    if( fstr_ctrl_get_param_ex( ctrl, 'VISUAL ',  '# ',    0,   'E',   visual )/= 0) return
380    if( fstr_ctrl_get_param_ex( ctrl, 'FEMAP ',   '# ',    0,   'E',   femap  )/= 0) return
381
382    fstr_ctrl_get_WRITE = 0
383
384  end function fstr_ctrl_get_WRITE
385
386  !> Read in !ECHO
387  function fstr_ctrl_get_ECHO( ctrl, echo )
388    integer(kind=kint) :: ctrl
389    integer(kind=kint) :: echo
390    integer(kind=kint) :: fstr_ctrl_get_ECHO
391
392    echo = kON;
393
394    fstr_ctrl_get_ECHO = 0
395
396  end function fstr_ctrl_get_ECHO
397
398  !> Read in !COUPLE
399  function fstr_ctrl_get_COUPLE( ctrl, fg_type, fg_first, fg_window, surf_id, surf_id_len )
400    integer(kind=kint) :: ctrl                           !< readed data
401    integer(kind=kint) :: fg_type                        !< if type
402    integer(kind=kint) :: fg_first                       !< if first
403    integer(kind=kint) :: fg_window                      !< if window
404    character(len=HECMW_NAME_LEN),target  :: surf_id(:)  !< surface id
405    character(len=HECMW_NAME_LEN),pointer :: surf_id_p   !< surface id
406    integer(kind=kint) :: surf_id_len
407    integer(kind=kint) :: fstr_ctrl_get_COUPLE
408
409    character(len=HECMW_NAME_LEN) :: data_fmt,ss
410    write(ss,*)  surf_id_len
411    write(data_fmt,'(a,a,a)') 'S',trim(adjustl(ss)),' '
412
413    fstr_ctrl_get_COUPLE = -1
414    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ', '1,2,3,4,5,6 ', 0, 'I', fg_type )/= 0) return
415    if( fstr_ctrl_get_param_ex( ctrl, 'ISTEP ', '# ', 0, 'I', fg_first )/= 0) return
416    if( fstr_ctrl_get_param_ex( ctrl, 'WINDOW ', '# ', 0, 'I', fg_window )/= 0) return
417
418    surf_id_p => surf_id(1)
419    fstr_ctrl_get_COUPLE = &
420      fstr_ctrl_get_data_array_ex( ctrl, data_fmt, surf_id_p )
421
422  end function fstr_ctrl_get_COUPLE
423
424  !> Read in !MPC
425  function fstr_ctrl_get_MPC( ctrl, penalty )
426    integer(kind=kint), intent(in) :: ctrl      !< readed data
427    real(kind=kreal), intent(out)  :: penalty   !< penalty
428    integer(kind=kint) :: fstr_ctrl_get_MPC
429
430    fstr_ctrl_get_MPC = fstr_ctrl_get_data_ex( ctrl, 1,   'r ', penalty )
431    if( penalty <= 1.0 ) then
432      if (myrank == 0) then
433        write(IMSG,*) "Warging : !MPC : too small penalty: ", penalty
434        write(*,*) "Warging : !MPC : too small penalty: ", penalty
435      endif
436    endif
437
438  end function fstr_ctrl_get_MPC
439
440  !> Read in !OUTPUT_RES & !OUTPUT_VIS
441  logical function fstr_ctrl_get_outitem( ctrl, hecMESH, outinfo )
442    use fstr_setup_util
443    use m_out
444    integer(kind=kint), intent(in)        :: ctrl      !< readed data
445    type (hecmwST_local_mesh), intent(in) :: hecMESH   !< mesh information
446    type( output_info ), intent(out)      :: outinfo   !< output information
447
448    integer(kind=kint) :: rcode, ipos
449    integer(kind=kint) :: n, i, j
450    character(len=HECMW_NAME_LEN) :: data_fmt, ss
451    character(len=HECMW_NAME_LEN), allocatable :: header_name(:), onoff(:), vtype(:)
452
453    write( ss, * )  HECMW_NAME_LEN
454    write( data_fmt, '(a,a,a,a,a)') 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), ' '
455    !  write( data_fmt, '(a,a,a,a,a,a,a)') 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), ' '
456
457    fstr_ctrl_get_outitem = .false.
458
459    outinfo%grp_id_name = "ALL"
460    rcode = fstr_ctrl_get_param_ex( ctrl, 'GROUP ', '# ', 0, 'S', outinfo%grp_id_name )
461    ipos = 0
462    rcode = fstr_ctrl_get_param_ex( ctrl, 'ACTION ', 'SUM ', 0, 'P', ipos )
463    outinfo%actn = ipos
464
465    n = fstr_ctrl_get_data_line_n( ctrl )
466    if( n == 0 ) return
467    allocate( header_name(n), onoff(n), vtype(n) )
468    header_name(:) = ""; vtype(:) = "";  onoff(:) = ""
469    rcode = fstr_ctrl_get_data_array_ex( ctrl, data_fmt, header_name, onoff )
470    !  rcode = fstr_ctrl_get_data_array_ex( ctrl, data_fmt, header_name, onoff, vtype )
471
472    do i = 1, n
473      do j = 1, outinfo%num_items
474        if( trim(header_name(i)) == outinfo%keyWord(j) ) then
475          outinfo%on(j) = .true.
476          if( trim(onoff(i)) == 'OFF' ) outinfo%on(j) = .false.
477          if( len( trim(vtype(i)) )>0 ) then
478            if( fstr_str2index( vtype(i), ipos ) ) then
479              outinfo%vtype(j) = ipos
480            else if( trim(vtype(i)) == "SCALER" ) then
481              outinfo%vtype(j) = -1
482            else if( trim(vtype(i)) == "VECTOR" ) then
483              outinfo%vtype(j) = -2
484            else if( trim(vtype(i)) == "SYMTENSOR" ) then
485              outinfo%vtype(j) = -3
486            else if( trim(vtype(i)) == "TENSOR" ) then
487              outinfo%vtype(j) = -4
488            endif
489          endif
490        endif
491      enddo
492    enddo
493
494    deallocate( header_name, onoff, vtype )
495    fstr_ctrl_get_outitem = .true.
496
497  end function fstr_ctrl_get_outitem
498
499  !> Read in !CONTACT
500  function fstr_ctrl_get_CONTACTALGO( ctrl, algo )
501    integer(kind=kint) :: ctrl
502    integer(kind=kint) :: algo
503    integer(kind=kint) :: fstr_ctrl_get_CONTACTALGO
504
505    integer(kind=kint) :: rcode
506    character(len=80) :: s
507    algo = kcaSLagrange
508    s = 'SLAGRANGE,ALAGRANGE '
509    rcode = fstr_ctrl_get_param_ex( ctrl, 'TYPE ', s, 0, 'P', algo )
510    fstr_ctrl_get_CONTACTALGO = rcode
511  end function fstr_ctrl_get_CONTACTALGO
512
513  !>  Read in contact definition
514  logical function fstr_ctrl_get_CONTACT( ctrl, n, contact, np, tp, ntol, ttol, ctAlgo )
515    integer(kind=kint), intent(in)    :: ctrl          !< ctrl file
516    integer(kind=kint), intent(in)    :: n             !< number of item defined in this section
517    integer(kind=kint), intent(in)    :: ctAlgo        !< contact algorithm
518    type(tContact), intent(out)       :: contact(n)    !< contact definition
519    real(kind=kreal), intent(out)      :: np             !< penalty along contact nomral
520    real(kind=kreal), intent(out)      :: tp             !< penalty along contact tangent
521    real(kind=kreal), intent(out)      :: ntol           !< tolrence along contact nomral
522    real(kind=kreal), intent(out)      :: ttol           !< tolrence along contact tangent
523
524    integer           :: rcode, ipt
525    character(len=30) :: s1 = 'TIED,GLUED,SSLID,FSLID '
526    character(len=HECMW_NAME_LEN) :: data_fmt,ss
527    character(len=HECMW_NAME_LEN) :: cp_name(n)
528    real(kind=kreal)  :: fcoeff(n),tPenalty(n)
529
530    tPenalty = 1.0d6
531
532    write(ss,*)  HECMW_NAME_LEN
533    write( data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)),'Rr '
534
535    fstr_ctrl_get_CONTACT = .false.
536    contact(1)%ctype = 1   ! pure slave-master contact; default value
537    contact(1)%algtype = CONTACTSSLID ! small sliding contact; default value
538    rcode = fstr_ctrl_get_param_ex( ctrl, 'INTERACTION ', s1, 0, 'P', contact(1)%algtype )
539    if( contact(1)%algtype==CONTACTGLUED ) contact(1)%algtype=CONTACTFSLID  ! not complemented yet
540    if( fstr_ctrl_get_param_ex( ctrl, 'GRPID ', '# ', 1, 'I', contact(1)%group )/=0) return
541    do rcode=2,n
542      contact(rcode)%ctype = contact(1)%ctype
543      contact(rcode)%group = contact(1)%group
544      contact(rcode)%algtype = contact(1)%algtype
545    end do
546    if(  fstr_ctrl_get_data_array_ex( ctrl, data_fmt, cp_name, fcoeff, tPenalty ) /= 0 ) return
547    do rcode=1,n
548      contact(rcode)%pair_name = cp_name(rcode)
549      contact(rcode)%fcoeff = fcoeff(rcode)
550      contact(rcode)%tPenalty = tPenalty(rcode)
551    enddo
552
553    np = 0.d0;  tp=0.d0
554    ntol = 0.d0;  ttol=0.d0
555    if( fstr_ctrl_get_param_ex( ctrl, 'NPENALTY ',  '# ',  0, 'R', np ) /= 0 ) return
556    if( fstr_ctrl_get_param_ex( ctrl, 'TPENALTY ', '# ', 0, 'R', tp ) /= 0 ) return
557    if( fstr_ctrl_get_param_ex( ctrl, 'NTOL ',  '# ',  0, 'R', ntol ) /= 0 ) return
558    if( fstr_ctrl_get_param_ex( ctrl, 'TTOL ', '# ', 0, 'R', ttol ) /= 0 ) return
559    fstr_ctrl_get_CONTACT = .true.
560  end function fstr_ctrl_get_CONTACT
561
562  !> Read in !ELEMOPT
563  function fstr_ctrl_get_ELEMOPT( ctrl, elemopt361 )
564    integer(kind=kint) :: ctrl
565    integer(kind=kint) :: elemopt361
566    integer(kind=kint) :: fstr_ctrl_get_ELEMOPT
567
568    character(72) :: o361list = 'IC,Bbar '
569
570    integer(kind=kint) :: o361
571
572    fstr_ctrl_get_ELEMOPT = -1
573
574    o361 = elemopt361 + 1
575
576    !* parameter in header line -----------------------------------------------------------------*!
577    if( fstr_ctrl_get_param_ex( ctrl, '361 ', o361list, 0, 'P', o361 ) /= 0) return
578
579    elemopt361 = o361 - 1
580
581    fstr_ctrl_get_ELEMOPT = 0
582
583  end function fstr_ctrl_get_ELEMOPT
584
585
586  !> Read in !AUTOINC_PARAM                                                             !
587  function fstr_get_AUTOINC( ctrl, aincparam )
588    implicit none
589    integer(kind=kint)    :: ctrl
590    type( tParamAutoInc ) :: aincparam !< auto increment paramter
591    integer(kind=kint) :: fstr_get_AUTOINC
592
593    integer(kind=kint) :: rcode
594    character(len=HECMW_NAME_LEN) :: data_fmt
595    character(len=128) :: msg
596    integer(kind=kint) :: bound_s(10), bound_l(10)
597    real(kind=kreal) :: Rs, Rl
598
599    fstr_get_AUTOINC = -1
600
601    bound_s(:) = 0
602    bound_l(:) = 0
603
604    !parameters
605    aincparam%name = ''
606    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', aincparam%name ) /=0 ) return
607
608    !read first line ( decrease criteria )
609    data_fmt = 'riiii '
610    rcode = fstr_ctrl_get_data_ex( ctrl, 1, data_fmt, Rs, &
611      &  bound_s(1), bound_s(2), bound_s(3), aincparam%NRtimes_s )
612    if( rcode /= 0 ) return
613    aincparam%ainc_Rs = Rs
614    aincparam%NRbound_s(knstMAXIT) = bound_s(1)
615    aincparam%NRbound_s(knstSUMIT) = bound_s(2)
616    aincparam%NRbound_s(knstCITER) = bound_s(3)
617
618    !read second line ( increase criteria )
619    data_fmt = 'riiii '
620    rcode = fstr_ctrl_get_data_ex( ctrl, 2, data_fmt, Rl, &
621      &  bound_l(1), bound_l(2), bound_l(3), aincparam%NRtimes_l )
622    if( rcode /= 0 ) return
623    aincparam%ainc_Rl = Rl
624    aincparam%NRbound_l(knstMAXIT) = bound_l(1)
625    aincparam%NRbound_l(knstSUMIT) = bound_l(2)
626    aincparam%NRbound_l(knstCITER) = bound_l(3)
627
628    !read third line ( cutback criteria )
629    data_fmt = 'ri '
630    rcode = fstr_ctrl_get_data_ex( ctrl, 3, data_fmt, &
631      &  aincparam%ainc_Rc, aincparam%CBbound )
632    if( rcode /= 0 ) return
633
634    !input check
635    rcode = 1
636    if( Rs<0.d0 .or. Rs>1.d0 ) then
637      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : decrease ratio Rs must 0 < Rs < 1.'
638    else if( any(bound_s<0) ) then
639      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : decrease NR bound must >= 0.'
640    else if( aincparam%NRtimes_s < 1 ) then
641      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : # of times to decrease must > 0.'
642    else if( Rl<1.d0 ) then
643      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : increase ratio Rl must > 1.'
644    else if( any(bound_l<0) ) then
645      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : increase NR bound must >= 0.'
646    else if( aincparam%NRtimes_l < 1 ) then
647      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : # of times to increase must > 0.'
648    elseif( aincparam%ainc_Rc<0.d0 .or. aincparam%ainc_Rc>1.d0 ) then
649      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : cutback decrease ratio Rc must 0 < Rc < 1.'
650    else if( aincparam%CBbound < 1 ) then
651      write(msg,*) 'fstr contol file error : !AUTOINC_PARAM : maximum # of cutback times must > 0.'
652    else
653      rcode =0
654    end if
655    if( rcode /= 0 ) then
656      write(*,*) trim(msg)
657      write(ILOG,*) trim(msg)
658      return
659    endif
660
661    fstr_get_AUTOINC = 0
662  end function fstr_get_AUTOINC
663
664  !> Read in !TIME_POINTS
665  function fstr_ctrl_get_TIMEPOINTS( ctrl, tp )
666    integer(kind=kint) :: ctrl
667    type(time_points)  :: tp
668    integer(kind=kint) :: fstr_ctrl_get_TIMEPOINTS
669
670    integer(kind=kint) :: i, n, rcode
671    logical            :: generate
672    real(kind=kreal)   :: stime, etime, interval
673
674    fstr_ctrl_get_TIMEPOINTS = -1
675
676    tp%name = ''
677    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', tp%name ) /=0 ) return
678    tp%range_type = 1
679    if( fstr_ctrl_get_param_ex( ctrl, 'TIME ', 'STEP,TOTAL ', 0, 'P', tp%range_type ) /= 0 ) return
680    generate = .false.
681    if( fstr_ctrl_get_param_ex( ctrl, 'GENERATE ',  '# ', 0, 'E', generate ) /= 0) return
682
683    if( generate ) then
684      stime = 0.d0;  etime = 0.d0;  interval = 1.d0
685      if( fstr_ctrl_get_data_ex( ctrl, 1, 'rrr ', stime, etime, interval ) /= 0) return
686      tp%n_points = int((etime-stime)/interval)+1
687      allocate(tp%points(tp%n_points))
688      do i=1,tp%n_points
689        tp%points(i) = stime + dble(i-1)*interval
690      end do
691    else
692      n = fstr_ctrl_get_data_line_n( ctrl )
693      if( n == 0 ) return
694      tp%n_points = n
695      allocate(tp%points(tp%n_points))
696      if( fstr_ctrl_get_data_array_ex( ctrl, 'r ', tp%points ) /= 0 ) return
697      do i=1,tp%n_points-1
698        if( tp%points(i) < tp%points(i+1) ) cycle
699        write(*,*) 'Error in reading !TIME_POINT: time points must be given in ascending order.'
700        return
701      end do
702    end if
703
704    fstr_ctrl_get_TIMEPOINTS = 0
705  end function fstr_ctrl_get_TIMEPOINTS
706
707
708end module fstr_ctrl_common
709