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
GeneralModule.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 !
24 !/******************************************************************************
25 ! *
26 ! * ELMER/FEM Solver main program
27 ! *
28 ! ******************************************************************************
29 ! *
30 ! * Authors: Juha Ruokolainen
31 ! * Email: Juha.Ruokolainen@csc.fi
32 ! * Web: http://www.csc.fi/elmer
33 ! * Address: CSC - IT Center for Science Ltd.
34 ! * Keilaranta 14
35 ! * 02101 Espoo, Finland
36 ! *
37 ! * Original Date: 02 Jun 1997
38 ! *
39 ! *****************************************************************************/
40 
42 
44 
46 
48 
51 
52 !------------------------------------------------------------------------------
54 !------------------------------------------------------------------------------
55  MODULE generalmodule
56 !------------------------------------------------------------------------------
57 
58  USE mainutils
59 
60 !------------------------------------------------------------------------------
61  IMPLICIT NONE
62 !------------------------------------------------------------------------------
63 
64  INTEGER :: Initialize
65 
66 !------------------------------------------------------------------------------
67 ! Local variables
68 !------------------------------------------------------------------------------
69 
70  INTEGER :: i,j,k,n,l,t,k1,k2,iter,Ndeg,istat,nproc,tlen,nthreads
71  CHARACTER(LEN=MAX_STRING_LEN) :: threads
72 
73  REAL(KIND=dp) :: s,dt,dtfunc
74  REAL(KIND=dP), POINTER :: WorkA(:,:,:) => NULL()
75  REAL(KIND=dp), POINTER, SAVE :: sTime(:), sStep(:), sInterval(:), sSize(:), &
76  steadyIt(:),nonlinIt(:),sPrevSizes(:,:),sPeriodic(:)
77 
78  TYPE(element_t),POINTER :: CurrentElement
79 
80  LOGICAL :: GotIt,Transient,Scanning,LastSaved
81 
82  INTEGER :: TimeIntervals,interval,timestep, &
83  TotalTimesteps,SavedSteps,CoupledMaxIter,CoupledMinIter
84 
85  INTEGER, POINTER, SAVE :: Timesteps(:),OutputIntervals(:), ActiveSolvers(:)
86  REAL(KIND=dp), POINTER, SAVE :: TimestepSizes(:,:)
87 
88  INTEGER(KIND=AddrInt) :: ControlProcedure
89 
90  LOGICAL :: InitDirichlet, ExecThis
91 
92  TYPE(elementtype_t),POINTER :: elmt
93 
94  TYPE(parenv_t), POINTER :: ParallelEnv
95 
96  CHARACTER(LEN=MAX_NAME_LEN) :: ModelName, eq, ExecCommand
97  CHARACTER(LEN=MAX_STRING_LEN) :: OutputFile, PostFile, RestartFile, &
98  OutputName=' ',PostName=' ', When, OptionString
99 
100  TYPE(variable_t), POINTER :: Var
101  TYPE(mesh_t), POINTER :: Mesh
102  TYPE(solver_t), POINTER :: Solver
103 
104  REAL(KIND=dp) :: CT0,RT0,tt
105 
106  LOGICAL :: FirstLoad = .TRUE., FirstTime=.TRUE., Found
107  LOGICAL :: Silent, Version, GotModelName
108 
109  INTEGER :: NoArgs
110 
111  INTEGER :: ExtrudeLevels
112  TYPE(mesh_t), POINTER :: ExtrudedMesh
113 
114  INTEGER :: omp_get_max_threads
115  INTEGER :: MyVerbosity = 1
116 
117 
118 #ifdef HAVE_TRILINOS
119 INTERFACE
120  SUBROUTINE trilinoscleanup() BIND(C,name='TrilinosCleanup')
121  IMPLICIT NONE
122  END SUBROUTINE trilinoscleanup
123 END INTERFACE
124 #endif
125 
126 
127  CONTAINS
128 
129 
130  ! This is a dirty hack that adds an instance of ResultOutputSolver to the list of Solvers.
131  ! The idea is that it is much easier for the end user to take into use the vtu output this way.
132  ! The solver itself has limited set of parameters needed and is therefore approapriate for this
133  ! kind of hack. It can of course be also added as a regular solver also.
134  !----------------------------------------------------------------------------------------------
135  SUBROUTINE addvtuoutputsolverhack()
136  TYPE(solver_t), POINTER :: abc(:), psolver
137  CHARACTER(LEN=MAX_NAME_LEN) :: str
138  INTEGER :: i,j,j2,j3,k,n
139  TYPE(valuelist_t), POINTER :: params
140  LOGICAL :: gotit, vtuformat
141 
142  str = listgetstring( currentmodel % Simulation,'Post File',gotit)
143  IF(.NOT. gotit) RETURN
144  k = index( str,'.vtu' )
145  vtuformat = ( k /= 0 )
146 
147  IF(.NOT. vtuformat ) RETURN
148 
149  CALL info('AddVtuOutputSolverHack','Adding ResultOutputSolver to write VTU output in file: '&
150  //trim(str(1:k-1)))
151 
152  CALL listremove( currentmodel % Simulation,'Post File')
153  n = currentmodel % NumberOfSolvers+1
154  ALLOCATE( abc(n) )
155  DO i=1,n-1
156  ! Def_Dofs is the only allocatable structure within Solver_t:
157  IF( ALLOCATED( currentmodel % Solvers(i) % Def_Dofs ) ) THEN
158  j = SIZE(currentmodel % Solvers(i) % Def_Dofs,1)
159  j2 = SIZE(currentmodel % Solvers(i) % Def_Dofs,2)
160  j3 = SIZE(currentmodel % Solvers(i) % Def_Dofs,3)
161  ALLOCATE( abc(i) % Def_Dofs(j,j2,j3) )
162  END IF
163 
164  ! Copy the content of the Solver structure
165  abc(i) = currentmodel % Solvers(i)
166 
167  ! Nullify the old structure since otherwise bad things may happen at deallocation
168  nullify( currentmodel % Solvers(i) % ActiveElements )
169  nullify( currentmodel % Solvers(i) % Values )
170  nullify( currentmodel % Solvers(i) % Mesh )
171  nullify( currentmodel % Solvers(i) % BlockMatrix )
172  nullify( currentmodel % Solvers(i) % Matrix )
173  nullify( currentmodel % Solvers(i) % Variable )
174  END DO
175 
176  ! Deallocate the old structure and set the pointer to the new one
177  DEALLOCATE( currentmodel % Solvers )
178  currentmodel % Solvers => abc
179  currentmodel % NumberOfSolvers = n
180 
181  ! Now create the ResultOutputSolver instance on-the-fly
182  nullify( currentmodel % Solvers(n) % Values )
183  currentmodel % Solvers(n) % PROCEDURE = 0
184  nullify( currentmodel % Solvers(n) % Matrix )
185  nullify( currentmodel % Solvers(n) % BlockMatrix )
186  nullify( currentmodel % Solvers(n) % Values )
187  nullify( currentmodel % Solvers(n) % Variable )
188  nullify( currentmodel % Solvers(n) % ActiveElements )
189  currentmodel % Solvers(n) % NumberOfActiveElements = 0
190  nullify( currentmodel % Solvers(n) % Values )
191  j = currentmodel % NumberOfBodies
192  ALLOCATE( currentmodel % Solvers(n) % Def_Dofs(10,j,6))
193  currentmodel % Solvers(n) % Def_Dofs(:,1:j,6) = -1
194 
195  ! Add some keywords to the list
196  CALL listaddstring(currentmodel % Solvers(n) % Values,&
197  'Procedure', 'ResultOutputSolve ResultOutputSolver',.false.)
198  CALL listaddstring(currentmodel % Solvers(n) % Values,'Output Format','vtu')
199  CALL listaddstring(currentmodel % Solvers(n) % Values,'Output File Name',str(1:k-1))
200  CALL listaddstring(currentmodel % Solvers(n) % Values,'Exec Solver','after saving')
201  CALL listaddstring(currentmodel % Solvers(n) % Values,'Equation','InternalVtuOutputSolver')
202  CALL listaddlogical(currentmodel % Solvers(n) % Values,'Save Geometry IDs',.true.)
203 
204  END SUBROUTINE addvtuoutputsolverhack
205 
206 
207 !------------------------------------------------------------------------------
209 !------------------------------------------------------------------------------
210  SUBROUTINE addsolvers()
211 !------------------------------------------------------------------------------
212  INTEGER :: i,j,k,nlen
213  LOGICAL :: initsolver, found
214 !------------------------------------------------------------------------------
215 
216  ! This is a hack that sets Equation flags True for the "Active Solvers".
217  ! The Equation flag is the legacy way of setting a Solver active and is still
218  ! used internally.
219  !----------------------------------------------------------------------------
220  DO i=1,currentmodel % NumberOfSolvers
221 
222  eq = listgetstring( currentmodel % Solvers(i) % Values,'Equation', found )
223 
224  IF ( found ) THEN
225  nlen = len_trim(eq)
226  DO j=1,currentmodel % NumberOFEquations
227  activesolvers => listgetintegerarray( currentmodel % Equations(j) % Values, &
228  'Active Solvers', found )
229  IF ( found ) THEN
230  DO k=1,SIZE(activesolvers)
231  IF ( activesolvers(k) == i ) THEN
232  CALL listaddlogical( currentmodel % Equations(j) % Values, eq(1:nlen), .true. )
233  EXIT
234  END IF
235  END DO
236  END IF
237  END DO
238  END IF
239  END DO
240 
241  DO i=1,currentmodel % NumberOfSolvers
242  eq = listgetstring( currentmodel % Solvers(i) % Values,'Equation', found )
243 
244  solver => currentmodel % Solvers(i)
245  initsolver = listgetlogical( solver % Values, 'Initialize', found )
246  IF ( found .AND. initsolver ) THEN
247  CALL freematrix( solver % Matrix )
248  CALL listaddlogical( solver % Values, 'Initialize', .false. )
249  END IF
250 
251  IF ( solver % PROCEDURE == 0 .OR. initsolver ) THEN
252  IF ( .NOT. ASSOCIATED( solver % Mesh ) ) THEN
253  solver % Mesh => currentmodel % Meshes
254  END IF
255  currentmodel % Solver => solver
256  CALL addequationbasics( solver, eq, transient )
257  CALL addequationsolution( solver, transient )
258  END IF
259  END DO
260 !------------------------------------------------------------------------------
261  END SUBROUTINE addsolvers
262 !------------------------------------------------------------------------------
263 
264 
265 !------------------------------------------------------------------------------
267 !------------------------------------------------------------------------------
268  SUBROUTINE addmeshcoordinatesandtime()
269 !------------------------------------------------------------------------------
270  TYPE(variable_t), POINTER :: dtvar
271 
272  nullify( solver )
273 
274  mesh => currentmodel % Meshes
275  DO WHILE( ASSOCIATED( mesh ) )
276  CALL variableadd( mesh % Variables, mesh,solver, &
277  'Coordinate 1',1,mesh % Nodes % x )
278 
279  CALL variableadd(mesh % Variables,mesh,solver, &
280  'Coordinate 2',1,mesh % Nodes % y )
281 
282  CALL variableadd(mesh % Variables,mesh,solver, &
283  'Coordinate 3',1,mesh % Nodes % z )
284 
285  CALL variableadd( mesh % Variables, mesh, solver, 'Time', 1, stime )
286  CALL variableadd( mesh % Variables, mesh, solver, 'Periodic Time', 1, speriodic )
287  CALL variableadd( mesh % Variables, mesh, solver, 'Timestep', 1, sstep )
288  CALL variableadd( mesh % Variables, mesh, solver, 'Timestep size', 1, ssize )
289  CALL variableadd( mesh % Variables, mesh, solver, 'Timestep interval', 1, sinterval )
290 
291  ! Save some previous timesteps for variable timestep multistep methods
292  dtvar => variableget( mesh % Variables, 'Timestep size' )
293  dtvar % PrevValues => sprevsizes
294 
295  CALL variableadd( mesh % Variables, mesh, solver, &
296  'nonlin iter', 1, nonlinit )
297  CALL variableadd( mesh % Variables, mesh, solver, &
298  'coupled iter', 1, steadyit )
299  mesh => mesh % Next
300  END DO
301 !------------------------------------------------------------------------------
302  END SUBROUTINE addmeshcoordinatesandtime
303 !------------------------------------------------------------------------------
304 
305 
306 !------------------------------------------------------------------------------
308 !------------------------------------------------------------------------------
309  SUBROUTINE setinitialconditions()
310 !------------------------------------------------------------------------------
311  USE defutils
312  INTEGER :: dofs
313  CHARACTER(LEN=MAX_NAME_LEN) :: str
314  LOGICAL :: found
315  TYPE(solver_t), POINTER :: solver
316  INTEGER, ALLOCATABLE :: indexes(:)
317  REAL(KIND=dp),ALLOCATABLE :: work(:)
318 
319  INTEGER :: i,j,k,l,m,vect_dof,real_dof,dim
320 
321  REAL(KIND=dp) :: nrm(3),t1(3),t2(3),vec(3),tmp(3),udot
322  TYPE(valuelist_t), POINTER :: bc
323  TYPE(nodes_t), SAVE :: nodes
324  LOGICAL :: nt_boundary
325  TYPE(element_t), POINTER :: element
326  TYPE(variable_t), POINTER :: var, vect_var
327 
328  dim = coordinatesystemdimension()
329 
330  IF (getlogical(getsimulation(),'Restart Before Initial Conditions',found)) THEN
331  CALL restart
332  CALL initcond
333  ELSE
334  CALL initcond
335  CALL restart
336  END IF
337 
338 !------------------------------------------------------------------------------
339 ! Make sure that initial values at boundaries are set correctly.
340 ! NOTE: This overrides the initial condition setting for field variables!!!!
341 !-------------------------------------------------------------------------------
342  initdirichlet = listgetlogical( currentmodel % Simulation, &
343  'Initialize Dirichlet Conditions', gotit )
344  IF ( .NOT. gotit ) initdirichlet = .true.
345 
346  vect_var => null()
347  IF ( initdirichlet ) THEN
348  mesh => currentmodel % Meshes
349  DO WHILE( ASSOCIATED(mesh) )
350  ALLOCATE( work(mesh % MaxElementDOFs) )
351  CALL setcurrentmesh( currentmodel, mesh )
352 
353  DO t = mesh % NumberOfBulkElements + 1, &
354  mesh % NumberOfBulkElements + mesh % NumberOfBoundaryElements
355 
356  element => mesh % Elements(t)
357 
358  ! Set also the current element pointer in the model structure to
359  ! reflect the element being processed:
360  ! ---------------------------------------------------------------
361  currentmodel % CurrentElement => element
362  n = element % TYPE % NumberOfNodes
363 
364  bc => getbc()
365 
366  var => mesh % Variables
367  DO WHILE( ASSOCIATED(var) )
368  solver => var % Solver
369  IF ( .NOT. ASSOCIATED(solver) ) solver => currentmodel % Solver
370 
371  str = listgetstring( solver % Values, 'Namespace', found )
372  IF (found) CALL listsetnamespace(trim(str))
373 
374 
375  IF ( var % DOFs <= 1 ) THEN
376  work(1:n) = getreal( bc,var % Name, gotit )
377  IF ( gotit ) THEN
378 
379  nt_boundary = .false.
380  IF ( getelementfamily() /= 1 ) THEN
381  k = len_trim(var % name)
382  vect_dof = ichar(var % Name(k:k))-ichar('0');
383  IF ( vect_dof>=1 .AND. vect_dof<= 3 ) THEN
384  nt_boundary = getlogical( bc, &
385  'normal-tangential '//var % name(1:k-2), gotit)
386 
387  IF ( nt_boundary ) THEN
388  nt_boundary = .false.
389  vect_var => mesh % Variables
390  DO WHILE( ASSOCIATED(vect_var) )
391  IF ( vect_var % Dofs>=dim ) THEN
392  DO real_dof=1,vect_var % Dofs
393  nt_boundary = ASSOCIATED(var % Values, &
394  vect_var % Values(real_dof::vect_var % Dofs))
395  IF ( nt_boundary ) EXIT
396  END DO
397  IF ( nt_boundary ) EXIT
398  END IF
399  vect_var => vect_var % Next
400  END DO
401  END IF
402 
403  IF ( nt_boundary ) THEN
404  CALL getelementnodes(nodes)
405  nrm = normalvector(element,nodes,0._dp,0._dp,.true.)
406  SELECT CASE(element % TYPE % DIMENSION)
407  CASE(1)
408  t1(1) = nrm(2)
409  t1(2) = -nrm(1)
410  CASE(2)
411  CALL tangentdirections(nrm,t1,t2)
412  END SELECT
413 
414  SELECT CASE(vect_dof)
415  CASE(1)
416  vec = nrm
417  CASE(2)
418  vec = t1
419  CASE(3)
420  vec = t2
421  END SELECT
422  END IF
423  END IF
424  END IF
425 
426  DO j=1,n
427  k = element % NodeIndexes(j)
428  IF ( ASSOCIATED(var % Perm) ) k = var % Perm(k)
429  IF ( k>0 ) THEN
430  IF ( nt_boundary ) THEN
431  DO l=1,dim
432  m = l+real_dof-vect_dof
433  tmp(l)=vect_var % Values(vect_var % Dofs*(k-1)+m)
434  END DO
435  udot = sum(vec(1:dim)*tmp(1:dim))
436  tmp(1:dim)=tmp(1:dim)+(work(j)-udot)*vec(1:dim)
437  DO l=1,dim
438  m = l+real_dof-vect_dof
439  vect_var % Values(vect_var % Dofs*(k-1)+m)=tmp(l)
440  END DO
441  ELSE
442  var % Values(k) = work(j)
443  END IF
444  END IF
445  END DO
446  END IF
447 
448  IF ( transient .AND. solver % TimeOrder==2 ) THEN
449  work(1:n) = getreal( bc, trim(var % Name) // ' Velocity', gotit )
450  IF ( gotit ) THEN
451  DO j=1,n
452  k = element % NodeIndexes(j)
453  IF ( ASSOCIATED(var % Perm) ) k = var % Perm(k)
454  IF ( k>0 ) var % PrevValues(k,1) = work(j)
455  END DO
456  END IF
457  work(1:n) = getreal( bc, trim(var % Name) // ' Acceleration', gotit )
458  IF ( gotit ) THEN
459  DO j=1,n
460  k = element % NodeIndexes(j)
461  IF ( ASSOCIATED(var % Perm) ) k = var % Perm(k)
462  IF ( k>0 ) var % PrevValues(k,2) = work(j)
463  END DO
464  END IF
465  END IF
466  ELSE
467  CALL listgetrealarray( bc, &
468  var % Name, worka, n, element % NodeIndexes, gotit )
469  IF ( gotit ) THEN
470  DO j=1,n
471  k = element % NodeIndexes(j)
472  DO l=1,min(SIZE(worka,1),var % DOFs)
473  IF ( ASSOCIATED(var % Perm) ) k = var % Perm(k)
474  IF ( k>0 ) var % Values(var % DOFs*(k-1)+l) = worka(l,1,j)
475  END DO
476  END DO
477  ELSE
478  END IF
479  END IF
480 
481  CALL listsetnamespace('')
482  var => var % Next
483  END DO
484  END DO
485  DEALLOCATE( work )
486  mesh => mesh % Next
487  END DO
488  END IF
489 !------------------------------------------------------------------------------
490  END SUBROUTINE setinitialconditions
491 !------------------------------------------------------------------------------
492 
493 
494 !------------------------------------------------------------------------------
495  SUBROUTINE initcond()
496 !------------------------------------------------------------------------------
497  USE defutils
498  TYPE(element_t), POINTER :: edge
499  INTEGER :: dofs,i,j,k,l
500  CHARACTER(LEN=MAX_NAME_LEN) :: str
501  LOGICAL :: found, thingstodo
502  TYPE(solver_t), POINTER :: solver
503  INTEGER, ALLOCATABLE :: indexes(:)
504  REAL(KIND=dp) :: val
505  REAL(KIND=dp),ALLOCATABLE :: work(:)
506  TYPE(valuelist_t), POINTER :: ic
507 !------------------------------------------------------------------------------
508 
509  mesh => currentmodel % Meshes
510  DO WHILE( ASSOCIATED( mesh ) )
511  ALLOCATE( indexes(mesh % MaxElementDOFs), work(mesh % MaxElementDOFs) )
512  CALL setcurrentmesh( currentmodel, mesh )
513 
514  ! First set the global variables and check whether there is anything left to do
515  thingstodo = .false.
516  DO j=1,currentmodel % NumberOfICs
517 
518  ic => currentmodel % ICs(j) % Values
519 
520  var => mesh % Variables
521  DO WHILE( ASSOCIATED(var) )
522 
523  solver => var % Solver
524  IF ( .NOT. ASSOCIATED(solver) ) solver => currentmodel % Solver
525 
526  str = listgetstring( solver % Values, 'Namespace', found )
527  IF (found) CALL listsetnamespace(trim(str))
528 
529  ! global variable
530  IF( SIZE( var % Values ) == var % DOFs ) THEN
531  val = listgetcreal( ic, var % Name, gotit )
532  IF( gotit ) THEN
533  WRITE( message,'(A,ES12.3)') 'Initializing global variable: > '&
534  //trim(var % Name)//' < to :',val
535  CALL info('InitCond', message,level=8)
536  var % Values = val
537  END IF
538  ELSE
539  thingstodo = thingstodo .OR. &
540  listcheckpresent( ic, trim(var % Name) )
541  thingstodo = thingstodo .OR. &
542  listcheckpresent( ic, trim(var % Name)//' Velocity' )
543  thingstodo = thingstodo .OR. &
544  listcheckpresent( ic, trim(var % Name)//' Acceleration' )
545  thingstodo = thingstodo .OR. &
546  listcheckpresent( ic, trim(var % Name)//' {e}' )
547  END IF
548  var => var % Next
549  END DO
550  END DO
551 
552  ! And now do the ordinary fields
553  !--------------------------------
554  IF( thingstodo ) THEN
555  DO t=1, mesh % NumberOfBulkElements+mesh % NumberOfBoundaryElements
556 
557  currentelement => mesh % Elements(t)
558 
559  i = currentelement % BodyId
560  IF( i == 0 ) cycle
561 
562  j = listgetinteger(currentmodel % Bodies(i) % Values, &
563  'Initial Condition',gotit, 1, currentmodel % NumberOfICs )
564  IF ( .NOT. gotit ) cycle
565 
566  ic => currentmodel % ICs(j) % Values
567  currentmodel % CurrentElement => currentelement
568  n = getelementnofnodes()
569 
570  var => mesh % Variables
571  DO WHILE( ASSOCIATED(var) )
572 
573 
574  solver => var % Solver
575  IF ( .NOT. ASSOCIATED(solver) ) solver => currentmodel % Solver
576 
577  str = listgetstring( solver % Values, 'Namespace', found )
578  IF (found) CALL listsetnamespace(trim(str))
579 
580  ! global variables were already set
581  IF( SIZE( var % Values ) == var % DOFs ) THEN
582  CONTINUE
583 
584  ELSE IF ( var % DOFs <= 1 ) THEN
585 
586  work(1:n) = getreal( ic, var % Name, gotit )
587  IF ( gotit ) THEN
588  dofs = getelementdofs( indexes, usolver=var % Solver )
589  DO k=1,n
590  k1 = indexes(k)
591  IF ( ASSOCIATED(var % Perm) ) k1 = var % Perm(k1)
592  IF ( k1>0 ) var % Values(k1) = work(k)
593  END DO
594  END IF
595 
596  IF ( transient .AND. solver % TimeOrder==2 ) THEN
597  work(1:n) = getreal( ic, trim(var % Name) // ' Velocity', gotit )
598  IF ( gotit ) THEN
599  dofs = getelementdofs( indexes, usolver=var % Solver )
600  DO k=1,n
601  k1 = indexes(k)
602  IF ( ASSOCIATED(var % Perm) ) k1 = var % Perm(k1)
603  IF ( k1>0 ) var % PrevValues(k1,1) = work(k)
604  END DO
605  END IF
606  work(1:n) = getreal( ic, trim(var % Name) // ' Acceleration', gotit )
607  IF ( gotit ) THEN
608  dofs = getelementdofs( indexes, usolver=var % Solver )
609  DO k=1,n
610  k1 = indexes(k)
611  IF ( ASSOCIATED(var % Perm) ) k1 = var % Perm(k1)
612  IF ( k1>0 ) var % PrevValues(k1,2) = work(k)
613  END DO
614  END IF
615  END IF
616 
617  IF(ASSOCIATED(mesh % Edges)) THEN
618  IF ( i<=mesh % NumberOfBulkElements) THEN
619  gotit = listcheckpresent( ic, trim(var % Name)//' {e}' )
620  IF ( gotit ) THEN
621  DO k=1,currentelement % TYPE % NumberOfedges
622  edge => mesh % Edges(currentelement % EdgeIndexes(k))
623  l = var % Perm(currentelement % EdgeIndexes(k)+mesh % NumberOfNodes)
624  IF ( l>0 ) THEN
625  CALL localbcintegral( ic, &
626  edge, edge % TYPE % NumberOfNodes, currentelement, n, &
627  trim(var % Name)//' {e}', work(1) )
628  var % Values(l) = work(1)
629  END IF
630  END DO
631  END IF
632  END IF
633  END IF
634 
635  ELSE
636  CALL listgetrealarray( ic, &
637  var % Name, worka, n, currentelement % NodeIndexes, gotit )
638 
639  IF ( gotit ) THEN
640  DO k=1,n
641  k1 = indexes(k)
642  DO l=1,min(SIZE(worka,1),var % DOFs)
643  IF ( ASSOCIATED(var % Perm) ) k1 = var % Perm(k1)
644  IF ( k1>0 ) var % Values(var % DOFs*(k1-1)+l) = worka(l,1,k)
645  END DO
646  END DO
647  END IF
648  END IF
649  CALL listsetnamespace('')
650  var => var % Next
651  END DO
652  END DO
653  END IF
654 
655  DEALLOCATE( indexes, work )
656  mesh => mesh % Next
657  END DO
658 
659 !------------------------------------------------------------------------------
660  END SUBROUTINE initcond
661 !------------------------------------------------------------------------------
662 
663 !------------------------------------------------------------------------------
665 !------------------------------------------------------------------------------
666  SUBROUTINE restart()
667 !------------------------------------------------------------------------------
668  USE defutils
669  LOGICAL :: gotit
670  INTEGER :: k,starttime
671 !------------------------------------------------------------------------------
672 
673 
674  restartfile = listgetstring( currentmodel % Simulation, &
675  'Restart File', gotit )
676 
677  IF ( gotit ) THEN
678  k = listgetinteger( currentmodel % Simulation,'Restart Position',gotit, &
679  minv=0 )
680 
681  mesh => currentmodel % Meshes
682  DO WHILE( ASSOCIATED(mesh) )
683 
684  IF ( len_trim(mesh % Name) > 0 ) THEN
685  outputname = trim(mesh % Name) // '/' // trim(restartfile)
686  ELSE
687  outputname = trim(restartfile)
688  END IF
689  IF ( parenv % PEs > 1 ) &
690  outputname = trim(outputname) // '.' // trim(i2s(parenv % MyPe))
691 
692  CALL setcurrentmesh( currentmodel, mesh )
693  CALL loadrestartfile( outputname,k,mesh )
694 
695  starttime = listgetconstreal( currentmodel % Simulation,'Restart Time',gotit)
696  IF( gotit ) THEN
697  var => variableget( mesh % Variables, 'Time' )
698  IF ( ASSOCIATED( var ) ) var % Values(1) = starttime
699  END IF
700 
701  mesh => mesh % Next
702 
703  END DO
704 
705 
706  END IF
707 !------------------------------------------------------------------------------
708  END SUBROUTINE restart
709 !------------------------------------------------------------------------------
710 
711 
712 !------------------------------------------------------------------------------
714 !------------------------------------------------------------------------------
715  SUBROUTINE execsimulation(TimeIntervals, CoupledMinIter, &
716  coupledmaxiter, outputintervals, transient, scanning)
717  INTEGER :: timeintervals,coupledminiter, coupledmaxiter,outputintervals(:)
718  LOGICAL :: transient,scanning
719 !------------------------------------------------------------------------------
720  INTEGER :: interval, timestep, i, j, k, n
721  REAL(KIND=dp) :: dt, ddt, dtfunc
722  INTEGER :: timeleft,cum_timestep
723  INTEGER, SAVE :: stepcount=0, realtimestep
724  LOGICAL :: execthis,steadystatereached=.false.
725 
726  REAL(KIND=dp) :: cumtime, maxerr, adaptivelimit, &
727  adaptivemintimestep, adaptivemaxtimestep, timeperiod
728  INTEGER :: smallestcount, adaptivekeepsmallest, stepcontrol=-1
729  LOGICAL :: adaptivetime = .true.
730 
731  TYPE(solver_t), POINTER :: solver
732 
733 ! REAL(KIND=dp) :: RealTime, newtime, prevtime=0, maxtime, exitcond
734  REAL(KIND=dp) :: newtime, prevtime=0, maxtime, exitcond
735  REAL(KIND=dp), ALLOCATABLE :: xx(:,:), xxnrm(:), yynrm(:), prevxx(:,:,:)
736 
737  !Something I (Jess) am adding in order to see if I can access
738  !the displacements from this level
739  LOGICAL :: gotit
740  TYPE(valuelist_t), POINTER :: solverparams
741  TYPE(variable_t), POINTER :: stresssol
742  REAL(KIND=dp), POINTER :: displacement(:)
743  CHARACTER(LEN=MAX_NAME_LEN) :: equation
744  INTEGER :: boundaryelementcount, elementcount, t, body_id
745  INTEGER :: bc_id, nelementnodes, nt, nk
746  INTEGER, POINTER :: index(:), nodeindexes(:), myperm(:)
747  TYPE(element_t), POINTER :: element, currentelement
748  TYPE(nodes_t) :: elementnodes
749  TYPE(mesh_t), POINTER :: mesh
750  LOGICAL :: isfsi
751  !End of Jess stuff
752 
753 !$omp parallel
754 !$ IF(.NOT.GaussPointsInitialized()) CALL GaussPointsInit
755 !$omp end parallel
756 
757  IF (myverbosity > 3) WRITE(*,*) 'NumberOfSolvers = ',&
758  currentmodel % NumberOfSolvers
759 
760  DO i=1,currentmodel % NumberOfSolvers
761  solver => currentmodel % Solvers(i)
762  IF ( solver % PROCEDURE==0 ) cycle
763  IF ( solver % SolverExecWhen == solver_exec_ahead_all ) THEN
764  IF (myverbosity > 3) WRITE(*,*) &
765  'Calling SolverActivate for Solver', i
766  CALL solveractivate( currentmodel,solver,dt,transient )
767  END IF
768  END DO
769 
770  DO interval = 1, timeintervals
771  stepcount = stepcount + timesteps(interval)
772  END DO
773 
774  cum_timestep = 0
775  ddt = 0.0d0
776  DO interval = 1,timeintervals
777 
778 !------------------------------------------------------------------------------
779  IF ( transient .OR. scanning ) THEN
780  dt = timestepsizes(interval,1)
781  ELSE
782  dt = 1
783  END IF
784 !------------------------------------------------------------------------------
785 ! go trough number of timesteps within an interval
786 !------------------------------------------------------------------------------
787  timeperiod = listgetcreal(currentmodel % Simulation, 'Time Period',gotit)
788  IF(.NOT.gotit) timeperiod = huge(timeperiod)
789 
790 
791  realtimestep = 1
792  DO timestep = 1,timesteps(interval)
793 
794  cum_timestep = cum_timestep + 1
795  sstep(1) = cum_timestep
796 
797  dtfunc = listgetconstreal( currentmodel % Simulation, &
798  'Timestep Function', gotit)
799  IF(gotit) THEN
800  CALL warn('ExecSimulation','Obsolite keyword > Timestep Function < , use > Timestep Size < instead')
801  ELSE
802  dtfunc = listgetcreal( currentmodel % Simulation, &
803  'Timestep Size', gotit)
804  END IF
805  IF(gotit) dt = dtfunc
806 
807 !------------------------------------------------------------------------------
808  stime(1) = stime(1) + dt
809  speriodic(1) = stime(1)
810  DO WHILE(speriodic(1) > timeperiod)
811  speriodic(1) = speriodic(1) - timeperiod
812  END DO
813 
814  ! Move the old timesteps one step down the ladder
815  IF(timestep > 1 .OR. interval > 1) THEN
816  DO i = SIZE(sprevsizes,2),2,-1
817  sprevsizes(1,i) = sprevsizes(1,i-1)
818  END DO
819  sprevsizes(1,1) = ssize(1)
820  END IF
821  ssize(1) = dt
822 
823  sinterval(1) = interval
824  IF (.NOT. transient ) steadyit(1) = steadyit(1) + 1
825 !------------------------------------------------------------------------------
826  IF ( parenv % MyPE == 0 ) THEN
827  CALL info( 'MAIN', ' ', level=3 )
828  CALL info( 'MAIN', '-------------------------------------', level=3 )
829 
830  IF ( transient .OR. scanning ) THEN
831  WRITE( message, * ) 'Time: ',trim(i2s(cum_timestep)),'/', &
832  trim(i2s(stepcount)), stime(1)
833  CALL info( 'MAIN', message, level=3 )
834 
835  newtime= realtime()
836 
837  IF( cum_timestep > 1 ) THEN
838  maxtime = listgetconstreal( currentmodel % Simulation,'Real Time Max',gotit)
839  IF( gotit ) THEN
840  WRITE( message,'(A,F8.3)') 'Fraction of real time left: ',&
841  1.0_dp-realtime() / maxtime
842  ELSE
843  timeleft = nint((stepcount-(cum_timestep-1))*(newtime-prevtime)/60._dp);
844  IF (timeleft > 120) THEN
845  WRITE( message, *) 'Estimated time left: ', &
846  trim(i2s(timeleft/60)),' hours.'
847  ELSE IF(timeleft > 60) THEN
848  WRITE( message, *) 'Estimated time left: 1 hour ', &
849  trim(i2s(mod(timeleft,60))), ' minutes.'
850  ELSE IF(timeleft >= 1) THEN
851  WRITE( message, *) 'Estimated time left: ', &
852  trim(i2s(timeleft)),' minutes.'
853  ELSE
854  WRITE( message, *) 'Estimated time left: less than a minute.'
855  END IF
856  END IF
857  CALL info( 'MAIN', message, level=3 )
858  END IF
859  prevtime = newtime
860  ELSE
861  WRITE( message, * ) 'Steady state iteration: ',cum_timestep
862  CALL info( 'MAIN', message, level=3 )
863  END IF
864 
865  CALL info( 'MAIN', '-------------------------------------', level=3 )
866  CALL info( 'MAIN', ' ', level=3 )
867  END IF
868 
869 !------------------------------------------------------------------------------
870 ! Solve any and all governing equations in the system
871 !------------------------------------------------------------------------------
872  adaptivetime = listgetlogical( currentmodel % Simulation, &
873  'Adaptive Timestepping', gotit )
874 
875  IF ( transient .AND. adaptivetime ) THEN
876  adaptivelimit = listgetconstreal( currentmodel % Simulation, &
877  'Adaptive Time Error', gotit )
878 
879  IF ( .NOT. gotit ) THEN
880  WRITE( message, * ) 'Adaptive Time Limit must be given for' // &
881  'adaptive stepping scheme.'
882  CALL fatal( 'ElmerSolver', message )
883  END IF
884 
885  adaptivemaxtimestep = listgetconstreal( currentmodel % Simulation, &
886  'Adaptive Max Timestep', gotit )
887  IF ( .NOT. gotit ) adaptivemaxtimestep = dt
888  adaptivemaxtimestep = min(adaptivemaxtimestep, dt)
889 
890  adaptivemintimestep = listgetconstreal( currentmodel % Simulation, &
891  'Adaptive Min Timestep', gotit )
892 
893  adaptivekeepsmallest = listgetinteger( currentmodel % Simulation, &
894  'Adaptive Keep Smallest', gotit, minv=0 )
895 
896  n = currentmodel % NumberOfSolvers
897  j = 0
898  k = 0
899  DO i=1,n
900  solver => currentmodel % Solvers(i)
901  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
902  IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
903  j = max( j, SIZE( solver % Variable % PrevValues,2 ) )
904  END IF
905  k = max( k, SIZE( solver % Variable % Values ) )
906  END IF
907  END DO
908  ALLOCATE( xx(n,k), yynrm(n), xxnrm(n), prevxx( n,k,j ) )
909 
910  cumtime = 0.0d0
911  IF ( ddt == 0.0d0 .OR. ddt > adaptivemaxtimestep ) ddt = adaptivemaxtimestep
912 
913  s = stime(1) - dt
914  smallestcount = 0
915  DO WHILE( cumtime < dt-1.0d-12 )
916  ddt = min( dt - cumtime, ddt )
917 
918  DO i=1,currentmodel % NumberOFSolvers
919  solver => currentmodel % Solvers(i)
920  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
921  n = SIZE( solver % Variable % Values )
922  xx(i,1:n) = solver % Variable % Values
923  xxnrm(i) = solver % Variable % Norm
924  IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
925  DO j=1,SIZE( solver % Variable % PrevValues,2 )
926  prevxx(i,1:n,j) = solver % Variable % PrevValues(:,j)
927  END DO
928  END IF
929  END IF
930  END DO
931 
932  stime(1) = s + cumtime + ddt
933  ssize(1) = ddt
934  CALL solveequations( currentmodel, ddt, transient, &
935  coupledminiter, coupledmaxiter, steadystatereached, realtimestep )
936 
937 
938  maxerr = listgetconstreal( currentmodel % Simulation, &
939  'Adaptive Error Measure', gotit )
940 
941  DO i=1,currentmodel % NumberOFSolvers
942  solver => currentmodel % Solvers(i)
943  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
944  n = SIZE(solver % Variable % Values)
945  yynrm(i) = solver % Variable % Norm
946  solver % Variable % Values = xx(i,1:n)
947  IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
948  DO j=1,SIZE( solver % Variable % PrevValues,2 )
949  solver % Variable % PrevValues(:,j) = prevxx(i,1:n,j)
950  END DO
951  END IF
952  END IF
953  END DO
954 
955  sstep(1) = ddt / 2
956  stime(1) = s + cumtime + ddt/2
957  CALL solveequations( currentmodel, ddt/2, transient, &
958  coupledminiter, coupledmaxiter, steadystatereached, realtimestep )
959  stime(1) = s + cumtime + ddt
960  CALL solveequations( currentmodel, ddt/2, transient, &
961  coupledminiter, coupledmaxiter, steadystatereached, realtimestep )
962 
963  maxerr = abs( maxerr - listgetconstreal( currentmodel % Simulation, &
964  'Adaptive Error Measure', gotit ) )
965 
966  IF ( .NOT. gotit ) THEN
967  maxerr = 0.0d0
968  DO i=1,currentmodel % NumberOFSolvers
969  solver => currentmodel % Solvers(i)
970  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
971  IF ( yynrm(i) /= solver % Variable % Norm ) THEN
972  maxerr = max(maxerr,abs(yynrm(i)-solver % Variable % Norm)/yynrm(i))
973  END IF
974  END IF
975  END DO
976  END IF
977 
978  IF ( maxerr < adaptivelimit .OR. ddt <= adaptivemintimestep ) THEN
979  cumtime = cumtime + ddt
980  realtimestep = realtimestep+1
981  IF ( smallestcount >= adaptivekeepsmallest .OR. stepcontrol > 0 ) THEN
982  ddt = min( 2*ddt, adaptivemaxtimestep )
983  stepcontrol = 1
984  smallestcount = 0
985  ELSE
986  stepcontrol = 0
987  smallestcount = smallestcount + 1
988  END IF
989  ELSE
990  DO i=1,currentmodel % NumberOFSolvers
991  solver => currentmodel % Solvers(i)
992  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
993  n = SIZE(solver % Variable % Values)
994  solver % Variable % Norm = xxnrm(i)
995  solver % Variable % Values = xx(i,1:n)
996  IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
997  DO j=1,SIZE( solver % Variable % PrevValues,2 )
998  solver % Variable % PrevValues(:,j) = prevxx(i,1:n,j)
999  END DO
1000  END IF
1001  END IF
1002  END DO
1003  ddt = ddt / 2
1004  stepcontrol = -1
1005  END IF
1006  WRITE(*,'(a,3e20.12)') 'Adaptive(cum,ddt,err): ', cumtime, ddt, maxerr
1007  END DO
1008  ssize(1) = dt
1009  stime(1) = s + dt
1010 
1011  DEALLOCATE( xx, xxnrm, yynrm, prevxx )
1012  ELSE ! Adaptive timestepping
1013  IF( myverbosity > 3) WRITE(*,*) 'Calling SolveEquations &
1014  (no transient and no adaptive)'
1015  CALL solveequations( currentmodel, dt, transient, &
1016  coupledminiter, coupledmaxiter, steadystatereached, realtimestep )
1017  realtimestep = realtimestep+1
1018  END IF
1019 
1020  !I'm (Jess) trying to access the displacements here
1021  solver => currentmodel % Solvers(2)
1022  solverparams => getsolverparams(solver)
1023  equation = getstring(solverparams, 'Equation', gotit)
1024 
1025 
1026 !------------------------------------------------------------------------------
1027 ! Save results to disk, if requested
1028 !------------------------------------------------------------------------------
1029 
1030  lastsaved = .false.
1031  IF( outputintervals(interval) /= 0 ) THEN
1032 
1033  CALL savetopost(0)
1034  k = mod( timestep-1, outputintervals(interval) )
1035  IF ( k == 0 .OR. steadystatereached ) THEN
1036 
1037  DO i=1,currentmodel % NumberOfSolvers
1038  solver => currentmodel % Solvers(i)
1039  IF ( solver % PROCEDURE == 0 ) cycle
1040  execthis = ( solver % SolverExecWhen == solver_exec_ahead_save)
1041  when = listgetstring( solver % Values, 'Exec Solver', gotit )
1042  IF ( gotit ) execthis = ( when == 'before saving')
1043  IF( execthis ) CALL solveractivate( currentmodel,solver,dt,transient )
1044  END DO
1045 
1046  CALL savecurrent(timestep)
1047  lastsaved = .true.
1048 
1049  DO i=1,currentmodel % NumberOfSolvers
1050  solver => currentmodel % Solvers(i)
1051  IF ( solver % PROCEDURE == 0 ) cycle
1052  execthis = ( solver % SolverExecWhen == solver_exec_after_save)
1053  when = listgetstring( solver % Values, 'Exec Solver', gotit )
1054  IF ( gotit ) execthis = ( when == 'after saving')
1055  IF( execthis ) CALL solveractivate( currentmodel,solver,dt,transient )
1056  END DO
1057  END IF
1058  END IF
1059 !------------------------------------------------------------------------------
1060 
1061  maxtime = listgetcreal( currentmodel % Simulation,'Real Time Max',gotit)
1062  IF( gotit .AND. realtime() - rt0 > maxtime ) THEN
1063  CALL info('ElmerSolver','Reached allowed maximum real time, exiting...')
1064  goto 100
1065  END IF
1066 
1067  exitcond = listgetcreal( currentmodel % Simulation,'Exit Condition',gotit)
1068  IF( gotit .AND. exitcond > 0.0_dp ) THEN
1069  CALL info('ElmerSolver','Found a positive exit condition, exiting...')
1070  goto 100
1071  END IF
1072 
1073 !------------------------------------------------------------------------------
1074 
1075  IF ( steadystatereached .AND. .NOT. (transient .OR. scanning) ) THEN
1076  IF ( timestep >= coupledminiter ) EXIT
1077  END IF
1078 
1079 !------------------------------------------------------------------------------
1080  END DO ! timestep within an iterval
1081 !------------------------------------------------------------------------------
1082 
1083 !------------------------------------------------------------------------------
1084  END DO ! timestep intervals, i.e. the simulation
1085 !------------------------------------------------------------------------------
1086 
1087 100 DO i=1,currentmodel % NumberOfSolvers
1088  solver => currentmodel % Solvers(i)
1089  IF ( solver % PROCEDURE == 0 ) cycle
1090  when = listgetstring( solver % Values, 'Exec Solver', gotit )
1091  IF ( gotit ) THEN
1092  IF ( when == 'after simulation' .OR. when == 'after all' ) THEN
1093  CALL solveractivate( currentmodel,solver,dt,transient )
1094  IF (ASSOCIATED(solver % Variable % Values) ) lastsaved = .false.
1095  END IF
1096  ELSE
1097  IF ( solver % SolverExecWhen == solver_exec_after_all ) THEN
1098  CALL solveractivate( currentmodel,solver,dt,transient )
1099  IF (ASSOCIATED(solver % Variable % Values) ) lastsaved = .false.
1100  END IF
1101  END IF
1102  END DO
1103 
1104  IF( .NOT. lastsaved ) THEN
1105  DO i=1,currentmodel % NumberOfSolvers
1106  solver => currentmodel % Solvers(i)
1107  IF ( solver % PROCEDURE == 0 ) cycle
1108  execthis = ( solver % SolverExecWhen == solver_exec_ahead_save)
1109  when = listgetstring( solver % Values, 'Exec Solver', gotit )
1110  IF ( gotit ) execthis = ( when == 'before saving')
1111  IF( execthis ) CALL solveractivate( currentmodel,solver,dt,transient )
1112  END DO
1113  END IF
1114 
1115 !------------------------------------------------------------------------------
1116  END SUBROUTINE execsimulation
1117 !------------------------------------------------------------------------------
1118 
1119 
1120 !------------------------------------------------------------------------------
1122 !------------------------------------------------------------------------------
1123  SUBROUTINE savecurrent( CurrentStep )
1124 !------------------------------------------------------------------------------
1125  INTEGER :: i, j,k,l,n,q,currentstep,nlen
1126  TYPE(variable_t), POINTER :: var
1127  LOGICAL :: eiganal, gotit
1128  CHARACTER(LEN=MAX_NAME_LEN) :: simul
1129  LOGICAL :: binaryoutput, saveall
1130 
1131  simul = listgetstring( currentmodel % Simulation, 'Simulation Type' )
1132 
1133  outputfile = listgetstring(currentmodel % Simulation,'Output File',gotit)
1134  IF ( gotit ) THEN
1135  IF ( parenv % PEs > 1 ) THEN
1136  DO i=1,max_name_len
1137  IF ( outputfile(i:i) == ' ' ) EXIT
1138  END DO
1139  outputfile(i:i) = '.'
1140  WRITE( outputfile(i+1:), '(a)' ) trim(i2s(parenv % MyPE))
1141  END IF
1142 
1143  binaryoutput = listgetlogical( currentmodel % Simulation,'Binary Output',gotit )
1144  IF ( .NOT.gotit ) binaryoutput = .false.
1145 
1146  saveall = .NOT.listgetlogical( currentmodel % Simulation,&
1147  'Omit unchanged variables in output',gotit )
1148  IF ( .NOT.gotit ) saveall = .true.
1149 
1150  mesh => currentmodel % Meshes
1151  DO WHILE( ASSOCIATED( mesh ) )
1152  IF ( mesh % OutputActive ) THEN
1153  nlen = len_trim(mesh % Name )
1154  IF ( nlen > 0 ) THEN
1155  outputname = mesh % Name(1:nlen) // '/' // trim(outputfile)
1156  ELSE
1157  outputname = outputfile
1158  END IF
1159 
1160  eiganal = .false.
1161  DO i=1,currentmodel % NumberOfSolvers
1162  IF ( ASSOCIATED( currentmodel % Solvers(i) % Mesh, mesh ) ) THEN
1163  eiganal = listgetlogical( currentmodel % Solvers(i) % Values, &
1164  'Eigen Analysis', gotit )
1165  eiganal = eiganal .OR. listgetlogical( currentmodel % Solvers(i) % Values, &
1166  'Harmonic Analysis', gotit )
1167 
1168  IF ( eiganal ) THEN
1169  var => currentmodel % Solvers(i) % Variable
1170  IF ( ASSOCIATED(var % EigenValues) ) THEN
1171  IF ( totaltimesteps == 1 ) THEN
1172  DO j=1,currentmodel % Solvers(i) % NOFEigenValues
1173  IF ( currentmodel % Solvers(i) % Matrix % COMPLEX ) THEN
1174 
1175  n = SIZE(var % Values)/var % DOFs
1176  DO k=1,n
1177  DO l=1,var % DOFs/2
1178  q = var % DOFs*(k-1)
1179  var % Values(q+l) = REAL(var % eigenvectors(j,q/2+l))
1180  var % Values(q+l+var % DOFs/2) = aimag(var % EigenVectors(j,q/2+l))
1181  END DO
1182  END DO
1183  ELSE
1184  var % Values = REAL( Var % EigenVectors(j,:) )
1185  END IF
1186  savedsteps = saveresult( outputname, mesh, &
1187  j, stime(1), binaryoutput, saveall )
1188  END DO
1189  ELSE
1190  j = min( currentstep, SIZE( var % EigenVectors,1 ) )
1191  IF ( currentmodel % Solvers(i) % Matrix % COMPLEX ) THEN
1192  n = SIZE(var % Values)/var % DOFs
1193  DO k=1,n
1194  DO l=1,var % DOFs/2
1195  q = var % DOFs*(k-1)
1196  var % Values(q+l) = REAL(var % eigenvectors(j,q/2+l))
1197  var % Values(q+l+var % DOFs/2) = aimag(var % EigenVectors(j,q/2+l))
1198  END DO
1199  END DO
1200  ELSE
1201  var % Values = REAL(var % eigenvectors(j,:))
1202  END IF
1203  savedsteps = saveresult( outputname, mesh, &
1204  currentstep, stime(1), binaryoutput, saveall )
1205  END IF
1206  var % Values = 0.0d0
1207  END IF
1208  END IF
1209  END IF
1210  END DO
1211 
1212  IF ( .NOT. eiganal ) THEN
1213  savedsteps = saveresult( outputname,mesh, nint(sstep(1)), &
1214  stime(1), binaryoutput, saveall )
1215  END IF
1216  END IF
1217  mesh => mesh % Next
1218  END DO
1219  ELSE
1220  mesh => currentmodel % Meshes
1221  DO WHILE( ASSOCIATED( mesh ) )
1222  IF ( mesh % OutputActive ) &
1223  mesh % SavesDone=mesh % SavesDone+1
1224  mesh => mesh % Next
1225  END DO
1226  END IF
1227  CALL savetopost(currentstep)
1228 !------------------------------------------------------------------------------
1229  END SUBROUTINE savecurrent
1230 !------------------------------------------------------------------------------
1231 
1232 !------------------------------------------------------------------------------
1234 !------------------------------------------------------------------------------
1235  SUBROUTINE savetopost(CurrentStep)
1236 !------------------------------------------------------------------------------
1237  TYPE(variable_t), POINTER :: var
1238  LOGICAL :: eiganal = .false., found
1239  INTEGER :: i, j,k,l,n,q,currentstep,nlen,timesteps,savedeigenvalues
1240  CHARACTER(LEN=MAX_NAME_LEN) :: simul, savewhich
1241 
1242  simul = listgetstring( currentmodel % Simulation, 'Simulation Type' )
1243 
1244  outputfile = listgetstring( currentmodel % Simulation,'Output File',gotit )
1245  IF ( gotit ) THEN
1246  IF ( parenv % PEs > 1 ) THEN
1247  DO i=1,max_name_len
1248  IF ( outputfile(i:i) == ' ' ) EXIT
1249  END DO
1250  outputfile(i:i) = '.'
1251  WRITE( outputfile(i+1:), '(a)' ) trim(i2s(parenv % MyPE))
1252  END IF
1253  END IF
1254 
1255  postfile = listgetstring( currentmodel % Simulation,'Post File',gotit )
1256  IF( .NOT. gotit ) RETURN
1257 
1258  IF ( parenv % PEs > 1 ) THEN
1259  DO i=1,max_name_len
1260  IF ( postfile(i:i) == ' ' ) EXIT
1261  END DO
1262  postfile(i:i) = '.'
1263  WRITE( postfile(i+1:), '(a)' ) trim(i2s(parenv % MyPE))
1264  END IF
1265 
1266  ! Loop over all meshes
1267  !---------------------------------------------------
1268  mesh => currentmodel % Meshes
1269  DO WHILE( ASSOCIATED( mesh ) )
1270  IF ( currentstep == 0 .AND. mesh % SavesDone > 0) THEN
1271  mesh => mesh % Next
1272  cycle
1273  END IF
1274 
1275  ! Check whether this mesh is active for saving
1276  !--------------------------------------------------
1277  IF ( mesh % OutputActive ) THEN
1278  nlen = len_trim(mesh % Name)
1279  IF ( nlen==0 .OR. filenamequalified(outputfile) ) THEN
1280  outputname = outputfile
1281  ELSE
1282  outputname = mesh % Name(1:nlen)//'/'//trim(outputfile)
1283  END IF
1284 
1285  IF ( nlen==0 .OR. filenamequalified(postfile) ) THEN
1286  postname = postfile
1287  ELSE
1288  postname = mesh % Name(1:nlen)//'/'//trim(postfile)
1289  END IF
1290 
1291  IF ( listgetlogical( currentmodel % Simulation,'Filename Numbering',gotit) ) THEN
1292  IF( currentstep == 0 ) THEN
1293  postname = nextfreefilename(postname)
1294  ELSE
1295  postname = nextfreefilename(postname,lastexisting = .true. )
1296  END IF
1297  END IF
1298 
1299  ! Set the Mesh pointer in the CurrentModel
1300  CALL setcurrentmesh( currentmodel, mesh )
1301 
1302  ! Use number of timesteps or number of eigenmodes
1303  !------------------------------------------------
1304  eiganal = .false.
1305  timesteps = totaltimesteps
1306  DO i=1,currentmodel % NumberOfSolvers
1307  IF (ASSOCIATED(currentmodel % Solvers(i) % Mesh, mesh)) THEN
1308  eiganal = listgetlogical( currentmodel % &
1309  solvers(i) % Values, 'Eigen Analysis', gotit )
1310 
1311  eiganal = eiganal .OR. listgetlogical( currentmodel % &
1312  solvers(i) % Values, 'Harmonic Analysis', gotit )
1313 
1314  IF ( eiganal ) timesteps = max( timesteps, &
1315  currentmodel % Solvers(i) % NOFEigenValues )
1316  END IF
1317  END DO
1318 
1319  DO i=1,currentmodel % NumberOfSolvers
1320  IF (ASSOCIATED(currentmodel % Solvers(i) % Mesh, mesh)) THEN
1321  eiganal = listgetlogical( currentmodel % &
1322  solvers(i) % Values, 'Eigen Analysis', gotit )
1323 
1324  eiganal = eiganal .OR. listgetlogical( currentmodel % &
1325  solvers(i) % Values, 'Harmonic Analysis', gotit )
1326 
1327  IF ( eiganal ) THEN
1328  savewhich = listgetstring( currentmodel % Solvers(i) % Values, &
1329  'Eigen and Harmonic Solution Output', found )
1330 
1331  savedeigenvalues = currentmodel % Solvers(i) % NOFEigenValues
1332  IF( totaltimesteps > 1 ) THEN
1333 ! SavedEiegnValues = MIN( CurrentStep, SIZE( Var % EigenVectors,1 ) )
1334  END IF
1335 
1336  DO j=1, savedeigenvalues
1337  var => mesh % Variables
1338  DO WHILE(ASSOCIATED(var))
1339  IF ( .NOT. ASSOCIATED(var % EigenValues) ) THEN
1340  var => var % Next
1341  cycle
1342  END IF
1343 
1344  IF ( currentmodel % Solvers(i) % Matrix % COMPLEX ) THEN
1345  IF(var % DOFs==1) THEN
1346  var => var % Next
1347  cycle
1348  END IF
1349 
1350  n = SIZE(var % Values)/var % DOFs
1351  DO k=1,n
1352  DO l=1,var % DOFs/2
1353  q = var % DOFs*(k-1)
1354  var % Values(q+l) = REAL(var % eigenvectors(j,q/2+l))
1355  var % Values(q+l+var % DOFs/2) = aimag(var % EigenVectors(j,q/2+l))
1356  END DO
1357  END DO
1358 
1359  ELSE
1360  SELECT CASE( savewhich )
1361  CASE('real part')
1362  var % Values = var % EigenVectors(j,:)
1363  CASE('imag part')
1364  var % Values = aimag(var % EigenVectors(j,:))
1365  CASE('abs value')
1366  var % Values = abs(var % EigenVectors(j,:))
1367  CASE('phase angle')
1368  var % Values = atan2(aimag(var % EigenVectors(j,:)), &
1369  REAL(var % eigenvectors(j,:)))
1370  CASE default
1371  var % CValues => var % EigenVectors(j,:)
1372  END SELECT
1373  END IF
1374  var => var % Next
1375  END DO
1376 
1377  IF ( currentstep > 0 ) THEN
1378  IF ( mesh % SavesDone /= 0 ) THEN
1379  IF( totaltimesteps == 1 ) THEN
1380  mesh % SavesDone = j
1381  ELSE
1382  mesh % SavesDone = currentstep
1383  END IF
1384  END IF
1385  CALL writepostfile( postname,outputname, currentmodel, &
1386  currentmodel % Solvers(i) % NOFEigenValues, .true. )
1387  END IF
1388  END DO
1389  EXIT
1390  END IF
1391  END IF
1392  END DO
1393 
1394  ! If this mesh has not been saved, then do so
1395  !----------------------------------------------------------------------------
1396  IF ( .NOT. eiganal .OR. currentstep == 0 ) THEN
1397  WRITE(6,*) 'timesteps = ', timesteps
1398  CALL writepostfile( postname, outputname, currentmodel, timesteps, .true. )
1399  END IF
1400 
1401  var => mesh % Variables
1402  DO WHILE(ASSOCIATED(var))
1403  IF ( ASSOCIATED(var % EigenValues) ) THEN
1404  var % Values = 0._dp
1405  var % CValues => null()
1406  END IF
1407  var => var % Next
1408  END DO
1409  END IF
1410  mesh => mesh % Next
1411  END DO
1412 !------------------------------------------------------------------------------
1413  END SUBROUTINE savetopost
1414 !------------------------------------------------------------------------------
1415 
1416 
1417 !------------------------------------------------------------------------------
1419 !------------------------------------------------------------------------------
1420 ! RECURSIVE SUBROUTINE FreeMesh(Mesh)
1421 ! !------------------------------------------------------------------------------
1422 ! TYPE(Mesh_t), POINTER :: Mesh
1423 ! !------------------------------------------------------------------------------
1424 ! IF (.NOT.ASSOCIATED(Mesh)) RETURN
1425 
1426 ! CALL FreeMesh(Mesh % Next)
1427 
1428 ! Mesh % Next => NULL()
1429 ! Mesh % Child => NULL()
1430 ! Mesh % Parent => NULL()
1431 
1432 ! CALL ReleaseMesh(Mesh)
1433 ! DEALLOCATE(Mesh)
1434 ! !------------------------------------------------------------------------------
1435 ! END SUBROUTINE FreeMesh
1436 ! !------------------------------------------------------------------------------
1437 
1438 
1439 ! !------------------------------------------------------------------------------
1440 ! !> Releases structures related to the Solver.
1441 ! !------------------------------------------------------------------------------
1442 ! SUBROUTINE FreeSolver(Solver)
1443 ! !------------------------------------------------------------------------------
1444 ! TYPE(Solver_t) :: Solver
1445 ! !------------------------------------------------------------------------------
1446 ! CALL FreeValueList(Solver % Values)
1447 ! CALL FreeMatrix(Solver % Matrix)
1448 ! IF (ALLOCATED(Solver % Def_Dofs)) DEALLOCATE(Solver % Def_Dofs)
1449 ! IF (ASSOCIATED(Solver % ActiveElements)) DEALLOCATE(Solver % ActiveElements)
1450 ! !------------------------------------------------------------------------------
1451 ! END SUBROUTINE FreeSolver
1452 ! !------------------------------------------------------------------------------
1453 
1454 
1455 ! !------------------------------------------------------------------------------
1456 ! !> Releases value list which includes all the sif definitions, for example.
1457 ! !------------------------------------------------------------------------------
1458 ! SUBROUTINE FreeValueList(List)
1459 ! !------------------------------------------------------------------------------
1460 ! TYPE(ValueList_t), POINTER :: List
1461 ! !------------------------------------------------------------------------------
1462 ! TYPE(ValueList_t), POINTER :: ptr
1463 
1464 ! ptr => List
1465 ! DO WHILE(ASSOCIATED(ptr))
1466 ! IF (ASSOCIATED(ptr % TValues)) DEALLOCATE(ptr % TValues)
1467 ! IF (ASSOCIATED(ptr % FValues)) DEALLOCATE(ptr % FValues)
1468 ! IF (ASSOCIATED(ptr % IValues)) DEALLOCATE(ptr % IValues)
1469 ! ptr => ptr % Next
1470 ! END DO
1471 ! !------------------------------------------------------------------------------
1472 ! END SUBROUTINE FreeValueList
1473 ! !------------------------------------------------------------------------------
1474 
1475 
1476 ! !------------------------------------------------------------------------------
1477 ! !> Releases the whole model.
1478 ! !------------------------------------------------------------------------------
1479 ! SUBROUTINE FreeModel(Model)
1480 ! !------------------------------------------------------------------------------
1481 ! TYPE(Model_t), POINTER :: Model
1482 ! !------------------------------------------------------------------------------
1483 ! TYPE(Matrix_t), POINTER :: A,B
1484 ! INTEGER :: i
1485 ! IF (.NOT.ASSOCIATED(Model)) RETURN
1486 
1487 ! CALL FreeMesh(Model % Meshes)
1488 
1489 ! CALL FreeValueList(Model % Constants)
1490 ! CALL FreeValueList(Model % Simulation)
1491 
1492 ! IF (ASSOCIATED(Model % BCs)) THEN
1493 ! DO i=1,Model % NumberOfBCs
1494 ! #if 0
1495 ! A => Model % BCs(i) % PMatrix
1496 ! IF (ASSOCIATED(A)) THEN
1497 ! DO WHILE( ASSOCIATED(A) )
1498 ! B => A % Child
1499 ! A % Child => NULL()
1500 ! A => B
1501 ! END DO
1502 ! CALL FreeMatrix(Model % BCs(i) % PMatrix)
1503 ! END IF
1504 ! #endif
1505 ! CALL FreeValueList( Model % BCs(i) % Values)
1506 ! END DO
1507 ! DEALLOCATE(Model % BCs)
1508 ! END IF
1509 
1510 ! DO i=1,Model % NumberOfSolvers
1511 ! CALL FreeSolver(Model % Solvers(i))
1512 ! END DO
1513 ! DEALLOCATE(Model % Solvers)
1514 
1515 ! IF (ASSOCIATED(Model % ICs)) THEN
1516 ! DO i=1,Model % NumberOfICs
1517 ! CALL FreeValueList( Model % ICs(i) % Values)
1518 ! END DO
1519 ! DEALLOCATE(Model % ICs)
1520 ! END IF
1521 
1522 ! IF (ASSOCIATED(Model % Bodies)) THEN
1523 ! DO i=1,Model % NumberOfBodies
1524 ! CALL FreeValueList( Model % Bodies(i) % Values)
1525 ! END DO
1526 ! DEALLOCATE(Model % Bodies)
1527 ! END IF
1528 
1529 ! IF (ASSOCIATED(Model % Equations)) THEN
1530 ! DO i=1,Model % NumberOfEquations
1531 ! CALL FreeValueList( Model % Equations(i) % Values)
1532 ! END DO
1533 ! DEALLOCATE(Model % Equations)
1534 ! END IF
1535 
1536 ! IF (ASSOCIATED(Model % BodyForces)) THEN
1537 ! DO i=1,Model % NumberOfBodyForces
1538 ! CALL FreeValueList( Model % BodyForces(i) % Values)
1539 ! END DO
1540 ! DEALLOCATE(Model % BodyForces)
1541 ! END IF
1542 
1543 ! Model=>NULL()
1544 ! !------------------------------------------------------------------------------
1545 ! END SUBROUTINE FreeModel
1546 !------------------------------------------------------------------------------
1547 !------------------------------------------------------------------------------
1549 !------------------------------------------------------------------------------
1550  FUNCTION elmerloadmodel(ModelName,BoundariesOnly,numprocs,mype ) RESULT( Model )
1551 !------------------------------------------------------------------------------
1552  USE testobject
1553 
1554  IMPLICIT NONE
1555 
1556  TYPE(t_global), POINTER :: global
1557 
1558  CHARACTER(LEN=*) :: modelname
1559  LOGICAL :: boundariesonly
1560 
1561  INTEGER, OPTIONAL :: numprocs,mype
1562 
1563  TYPE(model_t), POINTER :: model
1564 
1565 !------------------------------------------------------------------------------
1566  TYPE(mesh_t), POINTER :: mesh,mesh1,newmesh,oldmesh
1567  INTEGER :: i,j,k,l,s,nlen,eqn,meshkeep,meshlevels
1568  LOGICAL :: gotit,gotmesh,found,onemeshname, openfile, transient
1569  LOGICAL :: stat, single, meshgrading, newloadmesh
1570  TYPE(solver_t), POINTER :: solver
1571  INTEGER(KIND=AddrInt) :: initproc
1572  INTEGER, TARGET :: def_dofs(10,6)
1573  REAL(KIND=dp) :: meshpower
1574  REAL(KIND=dp), POINTER :: h(:)
1575  CHARACTER(LEN=MAX_NAME_LEN) :: name,elementdef,elementdef0
1576  CHARACTER(LEN=MAX_NAME_LEN) :: meshdir,meshname
1577  TYPE(valuelist_t), POINTER :: lst
1578  INTEGER, ALLOCATABLE :: edgedofs(:),facedofs(:)
1579 !------------------------------------------------------------------------------
1580 
1581  ALLOCATE( model )
1582  currentmodel => model
1583  nullify( model % Variables )
1584 
1585  meshdir = ' '
1586  meshname = ' '
1587 
1588  model % DIMENSION = 0
1589  model % NumberOfBoundaries = 0
1590  model % NumberOfBodies = 0
1591  model % NumberOfICs = 0
1592  model % NumberOfBCs = 0
1593  model % NumberOfEquations = 0
1594  model % NumberOfSolvers = 0
1595  model % NumberOfMaterials = 0
1596  model % NumberOfBodyForces = 0
1597 
1598  INQUIRE( unit=infileunit, opened=openfile )
1599  IF ( .NOT. openfile ) OPEN( unit=infileunit, file=modelname, status='OLD' )
1600  CALL loadinputfile( model,infileunit,modelname,meshdir,meshname, .true., .true. )
1601  rewind( infileunit )
1602  CALL loadinputfile( model,infileunit,modelname,meshdir,meshname, .true., .false. )
1603  IF ( .NOT. openfile ) CLOSE( infileunit )
1604 
1605  CALL initializeoutputlevel()
1606 
1607  transient=listgetstring(model % Simulation, &
1608  'Simulation Type',found)=='transient'
1609 
1610  def_dofs = -1; def_dofs(:,1)=1
1611 
1612  i = 1
1613  DO WHILE(i<=model % NumberOFSolvers)
1614 
1615  solver => model % Solvers(i)
1616  model % Solver => solver
1617 
1618  name = listgetstring( solver % Values, 'Procedure', found )
1619  IF ( found ) THEN
1620  initproc = getprocaddr( trim(name)//'_Init0', abort=.false. )
1621  IF ( initproc /= 0 ) THEN
1622  CALL execsolver( initproc, model, solver, &
1623  solver % dt, transient )
1624 
1625  solver => model % Solvers(i)
1626  model % Solver => solver
1627  END IF
1628  END IF
1629 
1630  gotmesh = listcheckpresent(solver % Values, 'Mesh')
1631 
1632  IF(.NOT.ALLOCATED(solver % Def_Dofs)) THEN
1633  ALLOCATE(solver % Def_Dofs(10,model % NumberOfBodies,6))
1634  solver % Def_Dofs = -1; solver % Def_Dofs(:,:,1)=1
1635  END IF
1636 
1637  ! Define what kind of element we are working with in this solver
1638  !-----------------------------------------------------------------
1639  elementdef = listgetstring( solver % Values, 'Element', stat )
1640 
1641  IF ( .NOT. stat ) THEN
1642  IF ( listgetlogical( solver % Values, &
1643  'Discontinuous Galerkin', stat ) ) THEN
1644  solver % Def_Dofs(:,:,4) = 0
1645  IF ( .NOT. gotmesh ) def_dofs(:,4) = max(def_dofs(:,4),0 )
1646  i=i+1
1647  cycle
1648  ELSE
1649  elementdef = "n:1"
1650  END IF
1651  END IF
1652 
1653  elementdef0 = elementdef
1654  DO WHILE(.true.)
1655  j = index( elementdef0, '-' )
1656  IF (j>0) THEN
1657  elementdef = elementdef0(1:j-1)
1658  ELSE
1659  elementdef = elementdef0
1660  END IF
1661  CALL getdefs( elementdef )
1662  IF(j>0) THEN
1663  elementdef0 = elementdef0(j+1:)
1664  ELSE
1665  EXIT
1666  END IF
1667  END DO
1668 
1669  i = i + 1
1670  END DO
1671 
1672  ! Check the mesh
1673  !--------------------------------------------------------
1674  name = listgetstring( model % Simulation, 'Mesh', gotit )
1675 
1676  onemeshname = .false.
1677  IF ( gotit ) THEN
1678  k = 1
1679  i = 1
1680  nlen = len_trim(name)
1681  DO WHILE( k<=nlen .AND. name(k:k) /= ' ' )
1682  meshdir(i:i) = name(k:k)
1683  meshname(i:i) = name(k:k)
1684  k = k + 1
1685  i = i + 1
1686  END DO
1687 
1688  DO WHILE( k<=nlen .AND. name(k:k) == ' ' )
1689  k = k + 1
1690  END DO
1691 
1692  IF ( k<=nlen ) THEN
1693  meshname(i:i) = '/'
1694  i = i + 1
1695  DO WHILE( name(k:k) /= ' ' )
1696  meshname(i:i) = name(k:k)
1697  k = k + 1
1698  i = i + 1
1699  END DO
1700  ELSE
1701  onemeshname = .true.
1702  meshdir = "." // char(0)
1703  END IF
1704  meshname(i:i) = char(0)
1705  END IF
1706 
1707  nullify( model % Meshes )
1708  IF ( meshdir(1:1) /= ' ' ) THEN
1709 
1710  CALL resettimer('LoadMesh')
1711  newloadmesh = listgetlogical( model % Simulation,'New Load Mesh',gotit)
1712  IF( .NOT. gotit ) newloadmesh = .true.
1713  IF( newloadmesh ) THEN
1714  model % Meshes => loadmesh2( model, meshdir, meshname, &
1715  boundariesonly, numprocs, mype, def_dofs )
1716  ELSE
1717  model % Meshes => loadmesh( model, meshdir, meshname, &
1718  boundariesonly, numprocs, mype, def_dofs(1,:) )
1719  END IF
1720 
1721  CALL checktimer('LoadMesh',level=5,delete=.true.)
1722 
1723  CALL setcoordinatesystem( model )
1724 
1725  meshlevels = listgetinteger( model % Simulation, 'Mesh Levels', gotit )
1726  IF ( .NOT. gotit ) meshlevels=1
1727 
1728  meshkeep = listgetinteger( model % Simulation, 'Mesh keep', gotit )
1729  IF ( .NOT. gotit ) meshkeep=meshlevels
1730 
1731  meshpower = listgetconstreal( model % Simulation, 'Mesh Grading Power',gotit)
1732  meshgrading = listgetlogical( model % Simulation, 'Mesh Keep Grading', gotit)
1733 
1734  DO i=2,meshlevels
1735  oldmesh => model % Meshes
1736 
1737  IF (meshgrading) THEN
1738  ALLOCATE(h(oldmesh % NumberOfNodes))
1739  model % Mesh => oldmesh
1740  CALL getnodalelementsize(model,meshpower,.false.,h)
1741  newmesh => splitmeshequal(oldmesh,h)
1742  DEALLOCATE(h)
1743  ELSE
1744  newmesh => splitmeshequal(oldmesh)
1745  END IF
1746 
1747  IF(ASSOCIATED(oldmesh % Faces)) THEN
1748  CALL findmeshedges(newmesh)
1749 
1750  ALLOCATE( edgedofs(newmesh % NumberOfBulkElements))
1751  ALLOCATE( facedofs(newmesh % NumberOfBulkElements))
1752  edgedofs = max(0,maxval(def_dofs(:,2)))
1753  facedofs = max(0,maxval(def_dofs(:,3)))
1754  CALL setmeshedgefacedofs(newmesh,edgedofs,facedofs)
1755  DEALLOCATE(edgedofs,facedofs)
1756 
1757  CALL setmeshmaxdofs(newmesh)
1758  IF(parenv % PEs>1) CALL sparedgenumbering(newmesh)
1759  IF(parenv % PEs>1) CALL sparfacenumbering(newmesh)
1760  END IF
1761 
1762  IF ( i>meshlevels-meshkeep+1 ) THEN
1763  newmesh % Next => oldmesh
1764  newmesh % Parent => oldmesh
1765  oldmesh % Child => newmesh
1766  newmesh % OutputActive = .true.
1767  oldmesh % OutputActive = .false.
1768  ELSE
1769  CALL releasemesh(oldmesh)
1770  END IF
1771  model % Meshes => newmesh
1772  END DO
1773 
1774 
1775  IF ( onemeshname ) THEN
1776  i = 0
1777  ELSE
1778  i = len_trim(meshname)
1779  DO WHILE( i>0 )
1780  IF (meshname(i:i) /= '/' ) THEN
1781  i = i-1
1782  ELSE
1783  EXIT
1784  END IF
1785  END DO
1786  END IF
1787 
1788  i = i + 1
1789  k = 1
1790  model % Meshes % Name = ' '
1791  DO WHILE( meshname(i:i) /= char(0) )
1792  model % Meshes % Name(k:k) = meshname(i:i)
1793  k = k + 1
1794  i = i + 1
1795  END DO
1796 
1797  ! Ok, give name also to the parent meshes as they might be saved too
1798  oldmesh => model % Meshes % Parent
1799  DO WHILE( ASSOCIATED( oldmesh ) )
1800  oldmesh % Name = trim(oldmesh % Child % Name)//'p'
1801  oldmesh => oldmesh % Parent
1802  END DO
1803 
1804  DO i=1,model % NumberOfSolvers
1805  model % Solvers(i) % Mesh => model % Meshes
1806  END DO
1807  END IF
1808 
1809  DO s=1,model % NumberOfSolvers
1810  name = listgetstring( model % Solvers(s) % Values, 'Mesh', gotit )
1811 
1812  IF( gotit ) THEN
1813  WRITE(message,'(A,I0)') 'Loading solver specific mesh > '//trim(name)// ' < for solver ',s
1814  CALL info('ElmerLoadModel',message,level=7)
1815 
1816  single=.false.
1817  IF ( name(1:8) == '-single ' ) THEN
1818  single=.true.
1819  name=name(9:)
1820  END IF
1821  onemeshname = .false.
1822  k = 1
1823  i = 1
1824  nlen = len_trim(name)
1825  meshname = ' '
1826  DO WHILE( k<=nlen .AND. name(k:k) /= ' ' )
1827  meshdir(i:i) = name(k:k)
1828  meshname(i:i) = name(k:k)
1829  k = k + 1
1830  i = i + 1
1831  END DO
1832 
1833  DO WHILE( k<=nlen .AND. name(k:k) == ' ' )
1834  k = k + 1
1835  END DO
1836 
1837  IF ( k<=nlen ) THEN
1838  meshname(i:i) = '/'
1839  i = i + 1
1840  DO WHILE( name(k:k) /= ' ' )
1841  meshname(i:i) = name(k:k)
1842  k = k + 1
1843  i = i + 1
1844  END DO
1845  ELSE
1846  onemeshname = .true.
1847  meshdir = "." // char(0)
1848  END IF
1849  meshname(i:i) = char(0)
1850 
1851  IF ( onemeshname ) THEN
1852  i = 0
1853  ELSE
1854  DO WHILE( i>0 .AND. meshname(i:i) /= '/' )
1855  i = i - 1
1856  END DO
1857  END IF
1858 
1859  mesh => model % Meshes
1860  found = .false.
1861  DO WHILE( ASSOCIATED( mesh ) )
1862  found = .true.
1863  k = 1
1864  j = i+1
1865  DO WHILE( meshname(j:j) /= char(0) )
1866  IF ( mesh % Name(k:k) /= meshname(j:j) ) THEN
1867  found = .false.
1868  EXIT
1869  END IF
1870  k = k + 1
1871  j = j + 1
1872  END DO
1873  IF ( len_trim(mesh % Name) /= k-1 ) found = .false.
1874  IF ( found ) EXIT
1875  mesh => mesh % Next
1876  END DO
1877 
1878  IF ( found ) THEN
1879  model % Solvers(s) % Mesh => mesh
1880  cycle
1881  END IF
1882 
1883  DO i=1,6
1884  DO j=1,8
1885  def_dofs(j,i) = maxval(model % Solvers(s) % Def_Dofs(j,:,i))
1886  END DO
1887  END DO
1888 
1889  newloadmesh = listgetlogical( model % Simulation,'New Load Mesh',gotit)
1890  IF(.NOT. gotit) newloadmesh = .true.
1891  IF ( single ) THEN
1892  IF( newloadmesh ) THEN
1893  model % Solvers(s) % Mesh => &
1894  loadmesh2( model,meshdir,meshname,boundariesonly,1,0,def_dofs )
1895  ELSE
1896  model % Solvers(s) % Mesh => &
1897  loadmesh( model,meshdir,meshname,boundariesonly,1,0,def_dofs(1,:) )
1898  END IF
1899  ELSE
1900  IF( newloadmesh ) THEN
1901  model % Solvers(s) % Mesh => &
1902  loadmesh2( model,meshdir,meshname,boundariesonly,numprocs,mype,def_dofs )
1903  ELSE
1904  model % Solvers(s) % Mesh => &
1905  loadmesh( model,meshdir,meshname,boundariesonly,numprocs,mype,def_dofs(1,:) )
1906  END IF
1907  END IF
1908  model % Solvers(s) % Mesh % OutputActive = .true.
1909 
1910 
1911  meshlevels = listgetinteger( model % Solvers(s) % Values, 'Mesh Levels', gotit )
1912  IF ( .NOT. gotit ) meshlevels=1
1913 
1914  meshkeep = listgetinteger( model % Solvers(s) % Values, 'Mesh keep', gotit )
1915  IF ( .NOT. gotit ) meshkeep=meshlevels
1916 
1917  meshpower = listgetconstreal( model % Simulation, 'Mesh Grading Power',gotit)
1918  meshgrading = listgetlogical( model % Simulation, 'Mesh Keep Grading', gotit)
1919 
1920 
1921  DO i=2,meshlevels
1922  oldmesh => model % Solvers(s) % Mesh
1923 
1924  IF (meshgrading) THEN
1925  ALLOCATE(h(oldmesh % NumberOfNodes))
1926  model % Mesh => oldmesh
1927  CALL getnodalelementsize(model,meshpower,.false.,h)
1928  newmesh => splitmeshequal(oldmesh,h)
1929  DEALLOCATE(h)
1930  ELSE
1931  newmesh => splitmeshequal(oldmesh)
1932  END IF
1933 
1934  IF(ASSOCIATED(oldmesh % Faces)) THEN
1935  CALL findmeshedges(newmesh)
1936 
1937  ALLOCATE( edgedofs(newmesh % NumberOfBulkElements))
1938  ALLOCATE( facedofs(newmesh % NumberOfBulkElements))
1939  edgedofs = max(0,maxval(def_dofs(:,2)))
1940  facedofs = max(0,maxval(def_dofs(:,3)))
1941  CALL setmeshedgefacedofs(newmesh,edgedofs,facedofs)
1942  DEALLOCATE(edgedofs,facedofs)
1943 
1944  CALL setmeshmaxdofs(newmesh)
1945  IF(parenv % PEs>1) CALL sparedgenumbering(newmesh)
1946  IF(parenv % PEs>1) CALL sparfacenumbering(newmesh)
1947  END IF
1948 
1949  IF ( i>meshlevels-meshkeep+1 ) THEN
1950  newmesh % Next => oldmesh
1951  newmesh % Parent => oldmesh
1952  oldmesh % Child => newmesh
1953  newmesh % Name = oldmesh % Name
1954  newmesh % OutputActive = .true.
1955  oldmesh % OutputActive = .false.
1956  ELSE
1957  CALL releasemesh(oldmesh)
1958  END IF
1959  model % Solvers(s) % Mesh => newmesh
1960  END DO
1961 
1962  IF ( onemeshname ) i = 0
1963 
1964  k = 1
1965  i = i + 1
1966  model % Solvers(s) % Mesh % Name = ' '
1967  DO WHILE( meshname(i:i) /= char(0) )
1968  model % Solvers(s) % Mesh % Name(k:k) = meshname(i:i)
1969  k = k + 1
1970  i = i + 1
1971  END DO
1972 
1973  IF ( ASSOCIATED( model % Meshes ) ) THEN
1974  mesh1 => model % Meshes
1975  DO WHILE( ASSOCIATED( mesh1 % Next ) )
1976  mesh1 => mesh1 % Next
1977  END DO
1978  mesh1 % Next => model % Solvers(s) % Mesh
1979  ELSE
1980  model % Meshes => model % Solvers(s) % Mesh
1981  END IF
1982  END IF
1983  END DO
1984 
1985  CALL setcoordinatesystem( model )
1986 
1987  IF ( outputpath == ' ' ) THEN
1988  DO i=1,max_name_len
1989  IF ( meshdir(i:i) == char(0) ) EXIT
1990  outputpath(i:i) = meshdir(i:i)
1991  END DO
1992  END IF
1993 
1994  mesh => model % Meshes
1995  DO WHILE( ASSOCIATED( mesh ) )
1996  CALL meshstabparams( mesh )
1997  mesh => mesh % Next
1998  END DO
1999 !------------------------------------------------------------------------------
2000 
2001  CONTAINS
2002 
2003 !------------------------------------------------------------------------------
2004  SUBROUTINE getdefs(ElementDef)
2005 !------------------------------------------------------------------------------
2006  CHARACTER(LEN=*) :: elementdef
2007 !------------------------------------------------------------------------------
2008  INTEGER :: ind(8),i,j,n
2009  INTEGER, POINTER :: gdofs(:,:), sdofs(:,:,:)
2010 
2011  ind = [1,2,3,4,5,6,7,8]
2012 
2013  IF (elementdef(1:5) == 'point' ) ind=1
2014  IF (elementdef(1:4) == 'line' ) ind=2
2015  IF (elementdef(1:3) == 'tri' ) ind=3
2016  IF (elementdef(1:4) == 'quad' ) ind=4
2017  IF (elementdef(1:5) == 'tetra' ) ind=5
2018  IF (elementdef(1:7) == 'pyramid' ) ind=6
2019  IF (elementdef(1:5) == 'prism' ) ind=7
2020  IF (elementdef(1:5) == 'brick' ) ind=8
2021  IF (elementdef(1:8) == 'tri_face' ) ind=9
2022  IF (elementdef(1:9) == 'quad_face' ) ind=10
2023 
2024  n = index(elementdef,'-')
2025  IF (n<=0) n=len_trim(elementdef)
2026 
2027  j = index( elementdef(1:n), 'n:' )
2028  IF ( j>0 ) THEN
2029  READ( elementdef(j+2:), * ) l
2030  solver % Def_Dofs(ind,:,1) = l
2031  IF (.NOT. gotmesh ) def_dofs(ind,1) = max(def_dofs(ind,1), l)
2032  END IF
2033 
2034  j = index( elementdef(1:n), 'e:' )
2035  IF ( j>0 ) THEN
2036  READ( elementdef(j+2:), * ) l
2037  solver % Def_Dofs(ind,:,2) = l
2038  IF ( .NOT. gotmesh ) def_dofs(ind,2) = max(def_dofs(ind,2), l )
2039  END IF
2040 
2041  j = index( elementdef(1:n), 'f:' )
2042  IF ( j>0 ) THEN
2043  READ( elementdef(j+2:), * ) l
2044  solver % Def_Dofs(ind,:,3) = l
2045  IF ( .NOT. gotmesh ) def_dofs(ind,3) = max(def_dofs(ind,3), l )
2046  END IF
2047 
2048  j = index( elementdef(1:n), 'd:' )
2049  IF ( j>0 ) THEN
2050  READ( elementdef(j+2:), * ) l
2051  solver % Def_Dofs(ind,:,4) = l
2052  IF ( .NOT. gotmesh ) def_dofs(ind,4) = max(def_dofs(ind,4), l )
2053  ELSE
2054  IF ( listgetlogical( solver % Values, &
2055  'Discontinuous Galerkin', stat ) ) THEN
2056  solver % Def_Dofs(ind,:,4) = 0
2057  IF ( .NOT. gotmesh ) def_dofs(ind,4) = max(def_dofs(ind,4),0 )
2058  END IF
2059  END IF
2060 
2061  j = index( elementdef(1:n), 'b:' )
2062  IF ( j>0 ) THEN
2063  READ( elementdef(j+2:), * ) l
2064  solver % Def_Dofs(ind,:,5) = l
2065  IF ( .NOT. gotmesh ) def_dofs(ind,5) = max(def_dofs(ind,5), l )
2066  END IF
2067 
2068  j = index( elementdef(1:n), 'p:' )
2069  IF ( j>0 ) THEN
2070  IF ( elementdef(j+2:j+2)=='%' ) THEN
2071  solver % Def_Dofs(ind,:,6) = 0
2072  ELSE
2073  READ( elementdef(j+2:), * ) l
2074  solver % Def_Dofs(ind,:,6) = l
2075  IF ( .NOT. gotmesh ) def_dofs(ind,6) = max(def_dofs(ind,6), l )
2076  END IF
2077  END IF
2078 
2079 !------------------------------------------------------------------------------
2080  END SUBROUTINE getdefs
2081 !------------------------------------------------------------------------------
2082 
2083 
2084  !------------------------------------------------------------------------------
2085  ! Initialize the log file output system for Messages
2086  !------------------------------------------------------------------------------
2087  SUBROUTINE initializeoutputlevel( )
2088 
2089  INTEGER, POINTER :: outputmask(:)
2090 
2091 
2092  minoutputlevel = listgetinteger( currentmodel % Simulation, &
2093  'Min Output Level', gotit )
2094 
2095  maxoutputlevel = listgetinteger( currentmodel % Simulation, &
2096  'Max Output Level', gotit )
2097 
2098  IF ( .NOT. gotit ) maxoutputlevel = 10
2099 
2100  outputmask => listgetintegerarray( currentmodel % Simulation, &
2101  'Output Level', gotit )
2102 
2103  IF ( gotit ) THEN
2104  DO i=1,SIZE(outputmask)
2105  outputlevelmask(i-1) = outputmask(i) /= 0
2106  END DO
2107  END IF
2108 
2109  DO i=0,31
2110  outputlevelmask(i) = outputlevelmask(i) .AND. &
2111  i >= minoutputlevel .AND. i <= maxoutputlevel
2112  END DO
2113 
2114  outputprefix = listgetlogical( currentmodel % Simulation, &
2115  'Output Prefix', gotit )
2116  IF ( .NOT. gotit ) outputprefix = .false.
2117 
2118  outputcaller = listgetlogical( currentmodel % Simulation, &
2119  'Output Caller', gotit )
2120  IF ( .NOT. gotit ) outputcaller = .true.
2121 
2122  END SUBROUTINE initializeoutputlevel
2123 !------------------------------------------------------------------------------
2124  END FUNCTION elmerloadmodel
2125 !------------------------------------------------------------------------------
2126 !------------------------------------------------------------------------------
2127 
2128  END MODULE generalmodule
2129 !------------------------------------------------------------------------------
2130 
generalmodule::setinitialconditions
subroutine setinitialconditions()
Sets initial conditions for the fields.
Definition: GeneralModule.F90:309
generalmodule::savecurrent
subroutine savecurrent(CurrentStep)
Saves current timestep to external files.
Definition: GeneralModule.F90:1123
generalmodule::savetopost
subroutine savetopost(CurrentStep)
Saves results file to post processing file of ElmerPost format, if requested.
Definition: GeneralModule.F90:1235
generalmodule::elmerloadmodel
type(model_t) function, pointer elmerloadmodel(ModelName, BoundariesOnly, numprocs, mype)
Release a mesh from the list of meshes.
Definition: GeneralModule.F90:1550
generalmodule::restart
subroutine restart()
Check if we are restarting are if yes, read in field values.
Definition: GeneralModule.F90:666
generalmodule::addvtuoutputsolverhack
subroutine addvtuoutputsolverhack()
Definition: GeneralModule.F90:135
testobject::t_global
Definition: ModTestObject.F90:5
generalmodule::addsolvers
subroutine addsolvers()
Adds flags for active solvers.
Definition: GeneralModule.F90:210
generalmodule::initcond
subroutine initcond()
Definition: GeneralModule.F90:495
generalmodule::addmeshcoordinatesandtime
subroutine addmeshcoordinatesandtime()
Adds coordinate and time variables to the current mesh structure.
Definition: GeneralModule.F90:268
generalmodule::execsimulation
subroutine execsimulation(TimeIntervals, CoupledMinIter, CoupledMaxIter, OutputIntervals, Transient, Scanning)
Execute the individual solvers in defined sequence.
Definition: GeneralModule.F90:715
initializeoutputlevel
subroutine initializeoutputlevel()
Definition: GeneralModule.F90:2087
generalmodule
ElmerLib Elmer library routines
Definition: GeneralModule.F90:55
getdefs
subroutine getdefs(ElementDef)
Definition: GeneralModule.F90:2004