ElmerFoamFSI  2.0
ElmerFoamFSI is fluid-solid interaction simulation application built up from OpenFOAM CFD and Elmer CSM coupled through the IMPACT multiphysics software integration infrastructure.
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Groups Pages
ElmerCSC.F90
Go to the documentation of this file.
1 !/*****************************************************************************/
2 ! *
3 ! * Elmer, A Finite Element Software for Multiphysical Problems
4 ! *
5 ! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 ! *
7 ! * This library is free software; you can redistribute it and/or
8 ! * modify it under the terms of the GNU Lesser General Public
9 ! * License as published by the Free Software Foundation; either
10 ! * version 2.1 of the License, or (at your option) any later version.
11 ! *
12 ! * This library is distributed in the hope that it will be useful,
13 ! * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ! * Lesser General Public License for more details.
16 ! *
17 ! * You should have received a copy of the GNU Lesser General Public
18 ! * License along with this library (in file ../LGPL-2.1); if not, write
19 ! * to the Free Software Foundation, Inc., 51 Franklin Street,
20 ! * Fifth Floor, Boston, MA 02110-1301 USA
21 ! *
22 ! *****************************************************************************/
23 
25 
26 !------------------------------------------------------------------------------
28 !------------------------------------------------------------------------------
29 
31 
32 !------------------------------------------------------------------------------
34 !------------------------------------------------------------------------------
35 
36 
37 #include "../config.h"
38 
39 ! ******************************************************************************
40 
41 SUBROUTINE initializer(global, runs, verbIn)
42  USE testobject
43  USE types
44  USE generalutils
45  USE parallelutils
46  USE timemodule
47 
48  IMPLICIT NONE
49 
50  include 'comf90.h'
51 
52  TYPE(t_global), POINTER :: global
53  INTEGER :: runs, counter, nindex, verbin
54 
55  REAL(KIND=dp) :: ct, rt
56 ! INTEGER, PARAMETER :: Init=0
57  CHARACTER(LEN=MAX_NAME_LEN) :: datestr, toutput
58  INTEGER, POINTER :: tempnodelist(:)
59  LOGICAL :: nodepresent
60  INTEGER :: lastfsielement
61 ! INTEGER, POINTER :: Conn2(:)
62 #if defined(MINGW32)
63  CALL set_stdio_bufs()
64 #endif
65 
66  WRITE(6,*) 'ElmerCSC:Initializer: Starting ...'
67 
68  !Allocate space for storing the verbosity
69  ALLOCATE(global%verbosity(1))
70  !Initialize verbosity to 1
71  global%verbosity = verbin
72 
73  WRITE(*,*) 'verbosity = ', global%verbosity
74 
75  CALL envir( 'ELMERSOLVER_OUTPUT_TOTAL'//char(0), toutput, tlen )
76  silent = toutput(1:1)=='0' .OR. toutput(1:5)=='false'
77 
78  ct = cputime()
79  rt = realtime()
80 
81  IF ( .NOT. silent ) THEN
82  datestr = formatdate()
83  WRITE( *,'(A,A)' ) "ELMER SOLVER (v " // version // ") STARTED AT: ", trim(datestr)
84  CALL flush(6)
85  END IF
86 
87  IF( myverbosity > 3) WRITE( *, * ) 'Calling ElmerInitialize'
88 
89  CALL elmerinitialize(runs)
90 
91  global%MyModel => currentmodel
92 
93  !Assigning stuff & printing out to check
94  global%SolverId = -1
95  DO i=1, currentmodel % NumberOfSolvers
96  solver => global % MyModel % Solvers(i)
97  solverparams => getsolverparams(solver)
98  equation = getstring(solverparams, 'Equation', gotit)
99  WRITE(*,*) 'Equation = ',equation
100  IF( equation .eq. 'nonlinear elasticity' ) THEN
101  global%SolverId = i
102  EXIT
103  ELSE IF( equation .eq. 'linear elasticity' ) THEN
104  global%SolverId = i
105  EXIT
106  END IF
107  END DO
108 
109  !Find out which BC tag is the FSI BC
110  global%FSIbcId = -1
111  DO t=1,currentmodel % NumberOfBCs
112  isfsi = listgetlogical(currentmodel % BCs(t) % Values,&
113  'FSI BC', gotit )
114  IF ( isfsi ) THEN
115  IF (myverbosity > 3) WRITE(*,*) 'BC', t, 'is FSI'
116  global%FSIbcId = t
117  END IF
118  END DO
119 
120  IF (global%FSIbcId == -1) THEN
121  WRITE(*,*) 'Not registering data no FSI BC!'
122  RETURN
123  END IF
124 
125  !Determine the number of elements & nodes on the
126  !FSI BC so we can allocate the arrays
127  mymesh => currentmodel % Meshes
128  global%nNodes = 0
129  global%nElem = 0
130 
131  !Allocate a temporay array to avoid repeating nodes
132  t = mymesh % NumberOfBulkElements
133  mycurrentelement => mymesh % Elements(t+1)
134  t = mycurrentelement % TYPE % DIMENSION
135  counter = (mymesh % NumberOfBoundaryElements)*(t + 1)
136  ALLOCATE(tempnodelist(counter))
137  tempnodelist = -1
138 
139  !Loop over the boundary elements
140  counter = 0
141  DO t = mymesh % NumberOfBulkElements+1, &
142  mymesh % NumberOfBulkElements + &
143  mymesh % NumberOfBoundaryElements
144 
145  mycurrentelement => mymesh % Elements(t)
146  bc_id = getbcid(mycurrentelement)
147 
148  IF ( bc_id == global%FSIbcId ) THEN
149  lastfsielement = t
150  IF( myverbosity > 3) THEN
151  WRITE(*,*) 'Initializer element = ',t,' bc_id = ',bc_id
152  END IF
153  global%nElem = global%nElem + 1
154  DO nt = 1,mycurrentelement % TYPE % NumberOfNodes
155  mynodeindexes => mycurrentelement % NodeIndexes
156  IF( myverbosity > 3) THEN
157  WRITE(*,*) 'Node index = ', mynodeindexes(nt)
158  END IF
159  nodepresent = .false.
160  DO nk=1,counter
161  IF (tempnodelist(nk) == mynodeindexes(nt)) THEN
162  nodepresent = .true.
163  EXIT
164  END IF
165  END DO
166  IF ( nodepresent .eqv. .false. ) THEN
167  counter = counter + 1
168  IF( myverbosity > 4) THEN
169  WRITE(*,*) 'counter = ', counter
170  END IF
171  tempnodelist(counter) = mynodeindexes(nt)
172  END IF
173  END DO
174  END IF
175  END DO
176 
177  IF( myverbosity > 5) THEN
178  WRITE(*,*) 'SIZE(TempNodeList) =', SIZE(tempnodelist)
179  DO t = 1,SIZE(tempnodelist)
180  WRITE(*,*) tempnodelist(t)
181  END DO
182  END IF
183 
184  global%nNodes = counter
185 
186  IF( myverbosity > 4) THEN
187  WRITE(*,*) 'global%nElem = ',global%nElem
188  WRITE(*,*) 'global%nNodes = ',global%nNodes
189  END IF
190 
191  !Allocate space for storing the NodeDisplacements
192  !This location will be registered with IMPACT
193  IF( myverbosity > 3) WRITE(*,*) 'global%nNodes =', global%nNodes
194  ALLOCATE(global%NodeDisplacements(3*global%nNodes))
195  !Masoud
196  ALLOCATE(global%PreviousNodeDisplacements(3*global%nNodes))
197  global%PreviousNodeDisplacements = 0.d0
198  ALLOCATE(global%IsLoadFreeNode(global%nNodes))
199  global%IsLoadFreeNode = .false.
200  !Masoud End
201  ALLOCATE(global%MyToElmerNodes(global%nNodes))
202  IF( myverbosity > 3) THEN
203  WRITE(*,*) 'SIZE(NodeDisplacements)=',SIZE(global%NodeDisplacements)
204  END IF
205  ALLOCATE(currentmodel%ElmerToMyNodes(currentmodel % NumberOfNodes))
206  global%MyToElmerNodes = -1
207  currentmodel%ElmerToMyNodes = -1
208 
209  !Allocate space for storing the FaceLoads
210  !This location will be registered with IMPACT
211  ALLOCATE(global%FaceLoads(3*global%nElem))
212  !Allocate space for storing the FacePressures
213  !This location will be registered with IMPACT
214  ALLOCATE(global%FacePressures(global%nElem))
215  !Allocate space for storing the NodePressures
216  !This location will be registered with IMPACT
217  ALLOCATE(global%NodePressures(global%nNodes))
218  !Allocate space for storing the NodeLoads
219  !This location will be registered with IMPACT
220  ALLOCATE(global%NodeLoads(3*global%nNodes))
221  !Allocate space for storing the previous loads
222  !so linear interpolation can be done
223  ALLOCATE(global%PreviousLoads(3,global%nNodes))
224  !Allocate the NodeLoads that gets passed to the user defined function
225  ALLOCATE(currentmodel%NodeLoadsPass(3,global%nNodes))
226  IF( myverbosity > 3) WRITE(*,*) 'SIZE(NodeLoads)=',SIZE(global%NodeLoads)
227  !Initialize the loads to 0.0
228  DO t = 1, global%nNodes
229  DO j =1,3
230 ! CurrentModel%NodeLoadsPass(j,t) = (t-1)*3.0 + j*1.0
231  currentmodel%NodeLoadsPass(j,t) = 0.0
232  global%PreviousLoads(j,t) = 0.0
233  END DO
234  END DO
235 
236  !Allocate the mesh arrays now that their sizes are known
237  mycurrentelement => mymesh % Elements(lastfsielement)
238  IF (mycurrentelement % TYPE % DIMENSION == 1) THEN
239  IF (mycurrentelement % TYPE % NumberOfNodes == 2) THEN
240  global%MeshType = ':b2:'
241  global%nConn = 2
242  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 3) THEN
243  global%MeshType = ':b3:'
244  global%nConn = 3
245  ELSE
246  CALL fatal( ' ', 'FSI elements of 1 dimension must have 2 or 3 nodes!' )
247  END IF
248  ELSE IF (mycurrentelement % TYPE % DIMENSION == 2) THEN
249  IF( myverbosity > 3) THEN
250  WRITE(*,*) 'FSI BC element NumberOfNodes =',&
251  mycurrentelement % TYPE % NumberOfNodes
252  END IF
253  IF (mycurrentelement % TYPE % NumberOfNodes == 3) THEN
254  global%MeshType = ':t3:'
255  global%nConn = 3
256  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 6) THEN
257  global%MeshType = ':t6:'
258  global%nConn = 6
259  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 4) THEN
260  global%MeshType = ':q4:'
261  global%nConn = 4
262  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 8) THEN
263  global%MeshType = ':q8:'
264  global%nConn = 8
265  ELSE
266  CALL fatal( ' ', 'FSI elements of 2 dimensions must have 3 or 6 node&
267  triangles or 4 or 8 node quadrilaterals!' )
268  END IF
269  ELSE
270  CALL fatal( ' ', 'FSI elements cannot have dimension greater than 2!' )
271  END IF
272 
273  IF( myverbosity > 3) THEN
274  WRITE(*,*) 'global%MeshType = ',global%MeshType
275  WRITE(*,*) 'global%nConn =',global%nConn
276  END IF
277  ALLOCATE(global%Conn(global%nConn*global%nElem))
278  ALLOCATE(global%Coords(3*global%nNodes))
279 
280  !Populate the mesh arrays
281  !Loop over the boundary elements
282  counter = 0
283  DO t = mymesh % NumberOfBulkElements+1, &
284  mymesh % NumberOfBulkElements + &
285  mymesh % NumberOfBoundaryElements
286 
287  !WRITE(*,*) "ElmerCSC element = ", t
288 
289  mycurrentelement => mymesh % Elements(t)
290  bc_id = getbcid(mycurrentelement)
291 
292  !For debugging - you can remove later
293  CALL getelementnodes(elementnodes,mycurrentelement)
294  !WRITE(*,*) 'Elmer elements:'
295  !DO nt = 1,MyCurrentElement % TYPE % NumberOfNodes
296  ! WRITE(*,*) nt,':',ElementNodes % x(nt),' ',ElementNodes % y(nt),' ',&
297  ! ElementNodes % z(nt)
298  !END DO
299  !end of debugging
300 
301  IF ( bc_id == global%FSIbcId ) THEN
302  counter = counter + 1
303  IF( myverbosity > 3) THEN
304  WRITE(*,*) 'element = ',t
305  END IF
306  CALL getelementnodes(elementnodes,mycurrentelement)
307 
308  global%FacePressures(counter) = 0.0
309  global%FaceLoads(3*(counter-1) + 1) = 0.0
310  global%FaceLoads(3*(counter-1) + 1) = 0.0
311  global%FaceLoads(3*(counter-1) + 1) = 0.0
312 
313  DO nt = 1,mycurrentelement % TYPE % NumberOfNodes
314  mynodeindexes => mycurrentelement % NodeIndexes
315 
316  IF( myverbosity > 3) THEN
317  WRITE(*,*) 'Node index = ', mynodeindexes(nt)
318  WRITE(*,*) elementnodes % x(nt), elementnodes % y(nt), &
319  elementnodes % z(nt)
320  END IF
321 
322  DO nk=1,global%nNodes
323  IF (tempnodelist(nk) == mynodeindexes(nt)) THEN
324  nindex = nk
325  EXIT
326  END IF
327  END DO
328 
329  !Conn2((counter-1)*2+nt) = 10
330  global%MyToElmerNodes(nindex) = mynodeindexes(nt)
331  currentmodel%ElmerToMyNodes(mynodeindexes(nt)) = nindex
332  global%Conn(global%nConn*(counter-1) + nt) = nindex
333  global%Coords(3*(nindex-1) + 1) = elementnodes % x(nt)
334  global%Coords(3*(nindex-1) + 2) = elementnodes % y(nt)
335  global%Coords(3*(nindex-1) + 3) = elementnodes % z(nt)
336  global%NodeLoads(3*(nindex-1) + 1) = 0.0
337  global%NodeLoads(3*(nindex-1) + 2) = 0.0
338  global%NodeLoads(3*(nindex-1) + 3) = 0.0
339  global%NodePressures(nindex) = 0.0
340  !Masoud: Marking LoadLessNodes
341  IF (abs(elementnodes % x(nt) - 0.6) < 10.0*epsilon(1.0)) THEN
342  !global%IsLoadFreeNode(nindex-1) = .TRUE.
343  END IF
344  !Masoud End
345  END DO
346  END IF
347  END DO
348 
349  !Masoud : printing IsLoadFreeNode
350  !WRITE(*,*) 'global%IsLoadFreeNode = ', global%IsLoadFreeNode
351  !Masoud End
352 
353  IF( myverbosity >= 5) THEN
354  WRITE(6,*) 'global%MyToElmerNodes = '
355  DO i=1,global%nNodes
356  WRITE(6,*) i,' = ',global%MyToElmerNodes(i)
357  END DO
358  WRITE(6,*) 'CurrentModel%ElmerToMyNodes = '
359  DO i=1,currentmodel % NumberOfNodes
360  WRITE(6,*) i,' = ',currentmodel%ElmerToMyNodes(i)
361  END DO
362 
363  WRITE(6,*)'-------------------------------------------'
364  WRITE(*,*) 'global%Conn = '
365  WRITE(6,*)'-------------------------------------------'
366  DO t = 1, global%nElem
367  WRITE(*,*) t,':',global%Conn(global%nConn*(t-1) + 1),&
368  global%Conn(global%nConn*(t-1) + 2)
369  END DO
370 
371  WRITE(6,*)
372 
373  WRITE(6,*)'-------------------------------------------'
374  WRITE(*,*) 'global%Coords = '
375  WRITE(6,*)'-------------------------------------------'
376  DO t = 1, global%nNodes
377  WRITE(*,*) t,':',global%Coords(3*(t-1)+1),global%Coords(3*(t-1)+2),&
378  global%Coords(3*(t-1)+3)
379  END DO
380  END IF
381 
382  WRITE(6,*)
383 
384  !Register the mesh
385  IF( myverbosity > 3) THEN
386  WRITE(*,*) 'window.mesh = ', &
387  trim(global%window_name)//'.'//trim(global%MeshType)
388  END IF
389  CALL com_set_size(trim(global%window_name)//'.nc',11,global%nNodes)
390  CALL com_set_array(trim(global%window_name)//'.nc',11,global%Coords,3)
391  CALL com_set_size(trim(global%window_name)//'.'//trim(global%MeshType),11,&
392  global%nElem)
393  CALL com_set_array(trim(global%window_name)//'.'//trim(global%MeshType),11,&
394  global%Conn, global%nConn)
395 
396  !CALL COM_SET_SIZE('Window1.:b2:real', 11,global%nElem)
397  !CALL COM_RESIZE_ARRAY('Window1.:b2:real',11)
398  !CALL COM_RESIZE_ARRAY(TRIM(global%window_name)//'.'//TRIM(global%MeshType),11)
399  !CALL COM_SET_ARRAY('Window1.:b2:real',11,Conn,2)
400 
401  !Set the displacments array with COM now that the mesh is registered
402  CALL com_new_dataitem(trim(global%window_name)//'.Displacements', 'n', com_double_precision, 3, 'm')
403  CALL com_set_array(trim(global%window_name)//'.Displacements',11,&
404  global%NodeDisplacements,3)
405  !Masoud
406  ! n: means node quantity, 3: number of columns, 11: window number (currently
407  ! default used everywhere), 'm': is units of the quanitity
408  CALL com_new_dataitem(trim(global%window_name)//'.PreviousDisplacements', 'n', com_double_precision, 3, 'm')
409  CALL com_set_array(trim(global%window_name)//'.PreviousDisplacements',11,&
410  global%PreviousNodeDisplacements,3)
411  !Masoud End
412 
413  !Set the loads array with COM now that the mesh is registered
414  CALL com_new_dataitem(trim(global%window_name)//'.Loads', 'n', com_double_precision, 3, '')
415  CALL com_set_array(trim(global%window_name)//'.Loads',11,&
416  global%NodeLoads,3)
417 
418  !Set the loads array with COM now that the mesh is registered
419  CALL com_new_dataitem(trim(global%window_name)//'.NodePressures',&
420  'n', com_double_precision, 1, '')
421  CALL com_set_array(trim(global%window_name)//'.NodePressures',11,&
422  global%NodePressures,1)
423 
424  !Set the pressures array with COM now that the mesh is registered
425  CALL com_new_dataitem(trim(global%window_name)//'.Pressures', 'e',&
426  com_double_precision, 1, '')
427  CALL com_set_array(trim(global%window_name)//'.Pressures',11,&
428  global%FacePressures,1)
429 
430  !Set the face loads array with COM now that the mesh is registered
431  CALL com_new_dataitem(trim(global%window_name)//'.FaceLoads', 'e',&
432  com_double_precision, 3, '')
433  CALL com_set_array(trim(global%window_name)//'.FaceLoads',11,&
434  global%FaceLoads,3)
435 
436  !Set the verbosity with COM
437  !CALL COM_NEW_DATAITEM(TRIM(global%window_name)//'.verbosity',&
438  ! 'w', COM_INTEGER, 1, '')
439  !CALL COM_SET_ARRAY(TRIM(global%window_name)//'.verbosity',11,&
440  ! global%verbosity,1)
441 
442  previoustime = 0.0
443 
444  IF( myverbosity > 4) THEN
445  WRITE(6,*)'-------------------------------------------'
446  WRITE(*,*) 'global%NodeDisplacements = '
447  WRITE(6,*)'-------------------------------------------'
448  DO i = 1,global%nNodes
449  DO j=1,3
450  global%NodeDisplacements(3*(i-1) + j) = 0.0
451  ENDDO
452  WRITE(6,*) (global%NodeDisplacements(3*(i-1) + j),j=1,3)
453  ENDDO
454 
455  WRITE(6,*)
456  END IF
457 
458 
459  WRITE(6,*) 'ElmerCSC:Initializer: Done.....'
460 
461  nullify(tempnodelist)
462 
463 END SUBROUTINE initializer
464 
465 ! *****************************************************************************
466 
467 SUBROUTINE run(global, runs, tFinal)
468 
469  USE testobject
470  USE generalmodule
471  USE timemodule
472 
473  IMPLICIT NONE
474 
475  include 'comf90.h'
476 
477  TYPE(t_global), POINTER :: global
478  INTEGER :: runs
479  DOUBLE PRECISION, INTENT(IN) :: tfinal
480  DOUBLE PRECISION :: standardtimestep
481  DOUBLE PRECISION :: var1, var2, var3
482  DOUBLE PRECISION :: deltatime
483  !INTEGER, PARAMETER :: Initialize=0
484  !INTEGER :: TimeIntervals, CoupledMinIter, CoupledMaxIter, Timestep
485  !INTEGER, POINTER, SAVE :: OutputIntervals(:)
486  !LOGICAL :: LastSaved, Transient, Scanning
487 
488  INTERFACE
489  SUBROUTINE updatedisplacements(global,runs,tFinal)
490  USE testobject
491  USE generalmodule
492  USE timemodule
493 
494  TYPE(t_global), POINTER :: global
495  INTEGER :: runs
496  DOUBLE PRECISION :: tfinal
497 
498  END SUBROUTINE updatedisplacements
499  END INTERFACE
500 
501  currentmodel => global%MyModel
502 
503  finaltime = tfinal
504 
505  IF (previoustime >= finaltime) THEN
506  CALL fatal( ' ', 'Next time step passed to RUN is greater than &
507  previous time step!')
508  END IF
509 
510  IF (myverbosity > 3) WRITE(*,*) 'In RUN function'
511 
512  standardtimestep = timestepsizes(1,1)
513  deltatime = tfinal - stime(1)
514  timesteps(1) = floor(deltatime/standardtimestep)
515 
516  IF (myverbosity > 3) THEN
517  WRITE(*,*) 'tFinal = ', tfinal
518  WRITE(*,*) 'sTime(1) = ', stime(1)
519  WRITE(*,*) 'standardTimestep = ', standardtimestep
520  END IF
521 
522  var3 = deltatime - int(deltatime/standardtimestep)*standardtimestep
523 
524  IF ( abs(var3) < 1.0d-12 ) THEN
525  timeintervals = 1
526  IF (myverbosity > 3) WRITE(*,*) 'standardTimestep evenly divides time'
527  ELSE
528  timestepsizes(2,1) = deltatime - timesteps(1)*standardtimestep
529  timeintervals = 2
530  timesteps(2) = 1
531  IF (outputintervals(1) == 1) THEN
532  outputintervals(2) = 1
533  END IF
534  IF (myverbosity > 3) THEN
535  WRITE(*,*) 'Remainder for time'
536  WRITE(*,*) 'TimestepSizes(2,1) = ', timestepsizes(2,1)
537  END IF
538  END IF
539 
540  IF (myverbosity > 3) THEN
541  WRITE(*,*) 'In Run function'
542  WRITE(*,*) 'tFinal = ', tfinal
543  WRITE(*,*) 'sTime(1) = ', stime(1)
544  WRITE(*,*) 'Timesteps(1) = ', timesteps(1)
545  WRITE(*,*) 'standardTimestep = ', standardtimestep
546  !Printing out the Pressures as a check
547  WRITE(*,*) 'Pressures '
548  DO t = 1, global%nElem
549  WRITE(*,*) global%FacePressures(t)
550  END DO
551  !Printing out the FaceLoads as a check
552  WRITE(*,*) 'FaceLoads '
553  DO t = 1, global%nElem
554  WRITE(*,*) global%FaceLoads(3*(t-1) + 1), global%FaceLoads(3*(t-1) + 2),&
555  global%FaceLoads(3*(t-1) + 3)
556  END DO
557  !Print out the NodePressures as a check
558  WRITE(*,*) 'NodePressures '
559  DO t = 1, global%nNodes
560  WRITE(*,*) global%NodePressures(t)
561  END DO
562 
563  END IF
564 
565 !------------------------------------------------------------------------------
566 ! Here we actually start the simulation ....
567 ! First go trough timeintervals
568 !------------------------------------------------------------------------------
569  !WRITE(*,*) 'global%nNodes = ',global%nNodes
570  IF (myverbosity > 3) WRITE(*,*) 'Calling TimeStepper'
571  CALL timestepper(global,runs)
572  IF (myverbosity > 3) WRITE(*,*) 'Done with TimeStepper'
573 
574  IF (global%FSIbcId /= -1) THEN
575  IF (myverbosity > 3) WRITE(*,*) 'Calling UpdateDisplacements'
576  CALL updatedisplacements(global,runs, tfinal)
577  END IF
578 
579  !Setting PreviousLoads values to NodeLoads now that
580  !Run is finished
581  DO t = 1, global%nNodes
582  DO j =1,3
583  global%PreviousLoads(j,t) = global%NodeLoads(3*(t-1) + j)
584  END DO
585  END DO
586 
587  !Setting PreviousTime to FinalTime now that
588  !Run is finished
589  previoustime = finaltime
590 
591  WRITE(*,*) 'ElmerCSC:Run: runs = ',runs
592  WRITE(*,*) 'ElmerCSC:Run: CurrentModel%GetTestLoads = ',&
593  currentmodel%GetTestLoads
594  WRITE(*,*) 'ElmerCSC:Run: CurrentModel%UDFUsed = ',&
595  currentmodel%UDFUsed
596  IF (runs == 1 .AND. currentmodel%GetTestLoads .eqv. .true. &
597  .AND. currentmodel%UDFUsed .eqv. .true.) THEN
598  runs = 2
599  END IF
600 
601 END SUBROUTINE run
602 
603 !**************************************************************************
604 
605 SUBROUTINE updatedisplacements(global,runs, tFinal)
606  USE testobject
607  USE generalmodule
608  USE timemodule
609 
610  IMPLICIT NONE
611 
612  include 'comf90.h'
613 
614  TYPE(t_global), POINTER :: global
615  INTEGER :: runs
616  DOUBLE PRECISION :: tfinal
617  INTEGER :: ncount, counter, ncountmax
618 
619  WRITE(*,*) 'ElmerCSC:UpdateDisplacements:',&
620  ' Reporting displacements to the fluid solver'
621 
622  !access the displacements here
623  solver => currentmodel % Solvers(global%SolverId)
624 
625  !IF (TRIM(Equation) == "nonlinear elasticity") THEN
626  stresssol => solver % Variable
627  displacement => stresssol % Values
628  myperm => stresssol % Perm
629 
630  IF( myverbosity > 3) THEN
631  WRITE(*,*) 'The SIZE(Displacement) = ', SIZE(displacement)
632  WRITE(*,*) 'After SolverActivate Call, Displacement(1) = ',&
633  displacement(1)
634  DO t=1,SIZE(displacement)
635  WRITE(*,*) 'Displacement',t,'=',displacement(t)
636  END DO
637  END IF
638 
639  !Jess adding stuff to try accessing wetted surface
640  elementcount = getnofactive()
641  boundaryelementcount = getnofboundaryelements()
642 
643  mymesh => currentmodel % Meshes
644 
645  !Write out nodes, coordinates, and connectivities
646  !Store at registered address
647  ncount = 0
648  ncountmax = 0
649  counter = 0
650  IF( myverbosity > 3) THEN
651  WRITE(*,*) 'Number of bulk elements = ', mymesh % NumberOfBulkElements
652  WRITE(*,*) 'Number of boundary elements = ', mymesh % NumberOfBoundaryElements
653  END IF
654  DO t = mymesh % NumberOfBulkElements+1, &
655  mymesh % NumberOfBulkElements + &
656  mymesh % NumberOfBoundaryElements
657  mycurrentelement => mymesh % Elements(t)
658  bc_id = getbcid(mycurrentelement)
659  IF ( bc_id == global%FSIbcId ) THEN
660  counter = counter+1
661  IF( myverbosity > 3) THEN
662  WRITE(*,*) '*************************************'
663  WRITE(*,*) 'element ',counter,'on FSI boundary'
664  WRITE(*,*) 'Update t = ',t,' bc_id = ',bc_id
665  WRITE(*,*)
666  WRITE(*,*) 'Element no. of nodes:', &
667  mycurrentelement % TYPE % NumberOfNodes
668  END IF
669 
670  CALL getelementnodes(elementnodes,mycurrentelement)
671 
672  DO nt = 1,mycurrentelement % TYPE % NumberOfNodes
673  mynodeindexes => mycurrentelement % NodeIndexes
674 
675  IF( myverbosity > 3) THEN
676  WRITE(*,*) '*********************'
677  WRITE(*,*) 'Node index = ', mynodeindexes(nt)
678  WRITE(*,*) elementnodes % x(nt), elementnodes % y(nt), &
679  elementnodes % z(nt)
680  WRITE(*,*) 'MyPerm = ', myperm(mynodeindexes(nt))
681  END IF
682 
683  IF ( myperm(mynodeindexes(nt)) > 0 ) THEN
684  nk = (stresssol % DOFs)*(myperm(mynodeindexes(nt)) - 1)
685  IF( myverbosity > 3) WRITE(*,*) 'DOFs = ', stresssol % DOFs
686 
687  ncount = currentmodel%ElmerToMyNodes(mynodeindexes(nt))
688 
689  ! Masoud : Avoding displacement update duplicates and
690  ! implementing a correct displacement update
691 
692  IF ( ncount > ncountmax ) THEN
693  IF( myverbosity > 3) THEN
694  WRITE(6,*) 'Updating NodeDisplacement(',ncount,')'
695  END IF
696 
697  IF ( stresssol % DOFs == 1 ) THEN
698  IF( myverbosity > 3) WRITE(*,*) displacement(nk+1), 0.0, 0.0
699  global%NodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)
700  global%NodeDisplacements(3*(ncount-1) + 2) = 0.0d0
701  global%NodeDisplacements(3*(ncount-1) + 3) = 0.0d0
702  ELSE IF ( stresssol % DOFs == 2 ) THEN
703  IF( myverbosity > 3) WRITE(*,*) displacement(nk+1), displacement(nk+2)
704  global%NodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)
705  global%NodeDisplacements(3*(ncount-1) + 2) = displacement(nk+2)
706  global%NodeDisplacements(3*(ncount-1) + 3) = 0.0d0
707  ELSE IF ( stresssol % DOFs == 3 ) THEN
708  ! Original
709  !global%NodeDisplacements(3*(nCount-1) + 1) = Displacement(nk+1)
710  !global%NodeDisplacements(3*(nCount-1) + 2) = Displacement(nk+2)
711  !global%NodeDisplacements(3*(nCount-1) + 3) = Displacement(nk+3)
712  ! Original: End
713 
714  ! Masoud
715  !Printing old and new displacements for the user
716  !WRITE(*,*) '----------------------'
717  !WRITE(*,*) ' new disps= ', Displacement(nk+1),&
718  ! Displacement(nk+2), Displacement(nk+3)
719  !WRITE(*,*) ' old disps= ',&
720  !global%PreviousNodeDisplacements(3*(nCount-1) + 1),&
721  !global%PreviousNodeDisplacements(3*(nCount-1) + 2), &
722  !global%PreviousNodeDisplacements(3*(nCount-1) + 3)
723  !Calculating displacement differences to pass to fluid solver
724  global%NodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)&
725  - global%PreviousNodeDisplacements(3*(ncount-1) + 1)
726  global%NodeDisplacements(3*(ncount-1) + 2) = displacement(nk+2)&
727  - global%PreviousNodeDisplacements(3*(ncount-1) + 2)
728  global%NodeDisplacements(3*(ncount-1) + 3) = displacement(nk+3)&
729  - global%PreviousNodeDisplacements(3*(ncount-1) + 3)
730  !Updating previous displacements to the current values for the
731  !next timestep
732  global%PreviousNodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)
733  global%PreviousNodeDisplacements(3*(ncount-1) + 2) = displacement(nk+2)
734  global%PreviousNodeDisplacements(3*(ncount-1) + 3) = displacement(nk+3)
735  !Priniting increament displacement values for the user
736  !WRITE(*,*) ' increament disps= ', global%NodeDisplacements(3*(nCount-1) + 1),&
737  !global%NodeDisplacements(3*(nCount-1) + 2), &
738  !global%NodeDisplacements(3*(nCount-1) + 3)
739  !WRITE(*,*) '----------------------'
740  ! Masoud : End
741  ELSE
742  WRITE(*,*) 'StressSol % DOFs = ', stresssol % DOFs
743  WRITE(*,*) 'DOFs are assumed to be <= 3'
744  CALL fatal( ' ', 'StressSol DOFs are greater than 3 ' )
745  END IF
746  ncountmax = max(ncount, ncountmax)
747  END IF
748  ! Masoud: End
749 
750  ! Original: uses a funny update order with dubplicates
751  !!IF( MyVerbosity > 3) THEN
752  ! WRITE(6,*) 'Updating NodeDisplacement(',nCount,')'
753  !!END IF
754 
755  !IF ( StressSol % DOFs == 1 ) THEN
756  ! IF( MyVerbosity > 3) WRITE(*,*) Displacement(nk+1), 0.0, 0.0
757  ! global%NodeDisplacements(3*(nCount-1) + 1) = Displacement(nk+1)
758  ! global%NodeDisplacements(3*(nCount-1) + 2) = 0.0d0
759  ! global%NodeDisplacements(3*(nCount-1) + 3) = 0.0d0
760  !ELSE IF ( StressSol % DOFs == 2 ) THEN
761  ! IF( MyVerbosity > 3) WRITE(*,*) Displacement(nk+1), Displacement(nk+2)
762  ! global%NodeDisplacements(3*(nCount-1) + 1) = Displacement(nk+1)
763  ! global%NodeDisplacements(3*(nCount-1) + 2) = Displacement(nk+2)
764  ! global%NodeDisplacements(3*(nCount-1) + 3) = 0.0d0
765  !ELSE IF ( StressSol % DOFs == 3 ) THEN
766  ! ! Masoud
767  ! !WRITE(*,*) '----------------------'
768  ! WRITE(*,*) ' new = ', Displacement(nk+1),&
769  ! Displacement(nk+2), Displacement(nk+3)
770  ! WRITE(*,*) ' old = ', global%NodeDisplacements(3*(nCount-1) + 1),&
771  ! global%NodeDisplacements(3*(nCount-1) + 2), &
772  ! global%NodeDisplacements(3*(nCount-1) + 3)
773  ! ! Masoud : end
774 
775  ! ! Original
776  ! global%NodeDisplacements(3*(nCount-1) + 1) = Displacement(nk+1)
777  ! global%NodeDisplacements(3*(nCount-1) + 2) = Displacement(nk+2)
778  ! global%NodeDisplacements(3*(nCount-1) + 3) = Displacement(nk+3)
779  ! ! Original: End
780 
781  ! !WRITE(*,*) '----------------------'
782  ! !global%NodeDisplacements(3*(nCount-1) + 1) = Displacement(nk+1) - global%NodeDisplacements(3*(nCount-1) + 1)
783  ! !global%NodeDisplacements(3*(nCount-1) + 2) = Displacement(nk+2) - global%NodeDisplacements(3*(nCount-1) + 2)
784  ! !global%NodeDisplacements(3*(nCount-1) + 3) = Displacement(nk+3) - global%NodeDisplacements(3*(nCount-1) + 3)
785  ! ! Masoud : End
786  !ELSE
787  ! WRITE(*,*) 'StressSol % DOFs = ', StressSol % DOFs
788  ! WRITE(*,*) 'DOFs are assumed to be <= 3'
789  ! CALL Fatal( ' ', 'StressSol DOFs are greater than 3 ' )
790  !END IF
791  ! Original : End
792  END IF
793 
794  END DO
795  END IF
796  END DO
797 
798  IF( myverbosity > 50) THEN
799  WRITE(6,*) 'tFinal = ', tfinal
800  WRITE(6,*) 'NodeDisplacements'
801  DO i = 1,global%nNodes
802  WRITE(6,*) (global%NodeDisplacements(3*(i-1) + j),j=1,3)
803  ENDDO
804  WRITE(6,*) 'Exiting UpdateDisplacements'
805  END IF
806 
807 END SUBROUTINE updatedisplacements
808 
809 ! ********************************************************************
810 
811 SUBROUTINE finalize(global, runs)
812  USE testobject
813  USE types
814  USE generalutils
815  USE parallelutils
816  USE generalmodule, only : firsttime
817 
818  IMPLICIT NONE
819 
820  include 'comf90.h'
821 
822  TYPE(t_global), POINTER :: global
823  INTEGER :: runs
824 
825  !INTEGER, PARAMETER :: Initialize=0
826  INTEGER :: tlen
827  LOGICAL :: silent
828  CHARACTER(LEN=MAX_NAME_LEN) :: datestr, toutput
829 
830  currentmodel => global%MyModel
831 
832  WRITE(*,*) 'ElmerCSC:Finalize: Finishing simulation'
833 
834  CALL elmerfinalize(runs)
835 
836  CALL envir( 'ELMERSOLVER_OUTPUT_TOTAL'//char(0), toutput, tlen )
837  silent = toutput(1:1)=='0' .OR. toutput(1:5)=='false'
838 
839  IF ( .NOT. silent ) THEN
840  IF ( parenv % myPE == 0 ) THEN
841  !WRITE( *,'(a,F12.2,F12.2)' ) 'SOLVER TOTAL TIME(CPU,REAL): ', &
842  ! CPUTime()-CT, RealTime()-RT
843  datestr = formatdate()
844  WRITE( *,'(A,A)' ) 'ELMER SOLVER FINISHED AT: ', trim(datestr)
845  END IF
846  END IF
847 
848  WRITE(6,*) 'ElmerCSC:Finalize: End Finalize. FirstTime = ',firsttime
849 
850 END SUBROUTINE finalize
851 
852 ! ********************************************************************
853 
854 SUBROUTINE elmercsc_load_module(name)
855 
856  USE testobject
857 
858  IMPLICIT NONE
859 
860  include "comf90.h"
861 
862  INTERFACE
863  SUBROUTINE com_set_pointer(attr,ptr,asso)
864  USE testobject
865  CHARACTER(*), INTENT(IN) :: attr
866  TYPE(t_global), POINTER :: ptr
867  EXTERNAL asso
868  END SUBROUTINE com_set_pointer
869 
870  SUBROUTINE initializer(global, runs, verbIn)
871  USE testobject
872  TYPE(t_global), POINTER :: global
873  INTEGER :: runs, verbin
874  END SUBROUTINE initializer
875 
876  SUBROUTINE run(global, runs, tFinal)
877  USE testobject
878  TYPE(t_global), POINTER :: global
879  INTEGER :: runs
880  DOUBLE PRECISION, INTENT(IN) :: tfinal
881  END SUBROUTINE run
882 
883  SUBROUTINE updatedisplacements(global, runs, tFinal)
884  USE testobject
885  TYPE(t_global), POINTER :: global
886  INTEGER :: runs
887  DOUBLE PRECISION :: tfinal
888  END SUBROUTINE updatedisplacements
889 
890  SUBROUTINE finalize(global, runs)
891  USE testobject
892  TYPE(t_global), POINTER :: global
893  INTEGER :: runs
894  END SUBROUTINE finalize
895 
896  END INTERFACE
897 
898  CHARACTER(*),intent(in) :: name
899  INTEGER :: com_types(7)
900  TYPE(t_global), POINTER :: glb
901 
902 
903  WRITE(*,'(A)') "Loading ElmerCSC: "//trim(name)
904 
905 
906  ALLOCATE(glb)
907  glb%window_name = trim(name)
908  glb%other_window_handle = -1
909  glb%c_window_handle = -1
910  CALL com_new_window(trim(name))
911 
912  CALL com_new_dataitem(trim(name)//'.global','w',com_f90pointer,1,'')
913  CALL com_allocate_array(trim(name)//'.global')
914 
915  com_types(1) = com_f90pointer
916  com_types(2) = com_integer
917  com_types(3) = com_integer
918 
919  CALL com_set_member_function(trim(name)//'.Initialize',initializer, &
920  trim(name)//'.global','bbi',com_types)
921 
922  com_types(3) = com_double_precision
923 
924  CALL com_set_member_function(trim(name)//'.Run',run, &
925  trim(name)//'.global','bbi',com_types)
926 
927  CALL com_set_member_function(trim(name)//'.Finalize',finalize, &
928  trim(name)//'.global','bb',com_types)
929 
930  CALL com_window_init_done(name)
931 
932  CALL com_set_pointer(name//'.global',glb,associate_pointer )
933 
934 END SUBROUTINE elmercsc_load_module
935 
936 
937 SUBROUTINE elmercsc_unload_module(name)
938  USE testobject
939  IMPLICIT NONE
940  include "comf90.h"
941  INTERFACE
942  SUBROUTINE com_get_pointer(attr,ptr,asso)
943  USE testobject
944  CHARACTER(*), INTENT(IN) :: attr
945  TYPE(t_global), POINTER :: ptr
946  EXTERNAL asso
947  END SUBROUTINE com_get_pointer
948  END INTERFACE
949  character(*),intent(in) :: name
950  TYPE(t_global), POINTER :: glb
951  INTEGER :: window_handle,other_window_handle,c_window_handle, owlen
952 
953  WRITE(*,'(A)') "Unloading ElmerCSC: "//trim(name)
954  nullify(glb)
955  window_handle = com_get_window_handle(trim(name))
956  if(window_handle .gt. 0) then
957  CALL com_get_pointer(trim(name)//'.global',glb,associate_pointer)
958  IF(ASSOCIATED(glb).eqv..true.) THEN
959  WRITE(*,'(A)') 'Fortran module '//trim(glb%window_name)//' unloading name '//trim(name)
960  if(glb%other_window_handle .gt. 0) then
961  WRITE(*,*) 'Fortran module '//trim(glb%window_name)//&
962  ' unloading external Fortran module '//trim(glb%other_window_name)//'.'
963  other_window_handle = com_get_window_handle(trim(glb%other_window_name))
964  IF(other_window_handle .gt. 0) THEN
965  CALL com_unload_module("ElmerCSC",trim(glb%other_window_name))
966  ENDIF
967  endif
968  if(glb%c_window_handle .gt. 0) then
969  WRITE(*,*) 'Fortran module '//trim(glb%window_name)//&
970  ' unloading external C module '//trim(glb%c_window_name)//'.'
971  c_window_handle = com_get_window_handle(trim(glb%c_window_name))
972  IF(c_window_handle .gt. 0) THEN
973  CALL com_unload_module("ElmerCSC",trim(glb%c_window_name))
974  ENDIF
975  endif
976 ! DEALLOCATE(glb)
977  ENDIF
978  CALL com_delete_window(trim(name))
979  endif
980 END SUBROUTINE elmercsc_unload_module
981 
982 ! ******************************************************************************
983 
timemodule
ElmerLib Elmer library routines
Definition: TimeModule.F90:55
timemodule::timestepper
subroutine timestepper(global, runs)
Definition: TimeModule.F90:184
elmerfinalize
subroutine elmerfinalize(runs)
ElmerLib Elmer library routines
Definition: Finalize.F90:55
finalize
subroutine finalize(global, runs)
Definition: ElmerCSC.F90:811
testobject::t_global
Definition: ModTestObject.F90:5
elmerinitialize
subroutine elmerinitialize(runs)
ElmerLib Elmer library routines
Definition: Initialize.F90:55
elmercsc_load_module
subroutine elmercsc_load_module(name)
Definition: ElmerCSC.F90:854
updatedisplacements
subroutine updatedisplacements(global, runs, tFinal)
Definition: ElmerCSC.F90:605
testobject::associate_pointer
subroutine associate_pointer(attr, ptr)
Definition: ModTestObject.F90:30
generalmodule
ElmerLib Elmer library routines
Definition: GeneralModule.F90:55
initializer
subroutine initializer(global, runs, verbIn)
A caller for the Elmer main program.
Definition: ElmerCSC.F90:41
elmercsc_unload_module
subroutine elmercsc_unload_module(name)
Definition: ElmerCSC.F90:937
run
subroutine run(global, runs, tFinal)
Definition: ElmerCSC.F90:467