OpenCMISS-Iron Internal API Documentation
fluid_mechanics_IO_routines.f90
Go to the documentation of this file.
1 
44 
46 
48 
49  USE base_routines
51  USE field_routines
52  USE types
53  USE input_output
54  USE kinds
55 
56 #include "macros.h"
57 #include "dllexport.h"
58 
59  IMPLICIT NONE
60 
61  PRIVATE
62 
63  !Module parameters
64 
65  !Module types
66 
67  !1=M, 2=V, 3=P !
69  INTEGER(INTG) lx ! External Structure: Number of Tessellation Coefficients
70  INTEGER(INTG) lt ! External Structure: Number of Tessellation elements
71  INTEGER(INTG) id ! Internal Structure: Basis number to call letter map
72  INTEGER(INTG) lb
73  REAL(DP), pointer :: x(:,:) ! External Structure: Tessellation coefficients
74  INTEGER(INTG), pointer :: t(:,:) ! External Structure: Tessellation map
75  INTEGER(INTG), pointer :: nbt(:,:) ! Internal Structure: Neighboring element map
76  INTEGER(INTG), pointer :: b(:,:)
77  END TYPE array_mesh
78 
79  TYPE array_int
80  REAL(DP), pointer :: y(:,:) ! Internal Structure: Basis @ given volume quadrature points
81  REAL(DP), pointer :: y_f(:,:,:) ! Internal Structure: Basis @ given facet quadrature points
82  REAL(DP), pointer :: dy(:,:,:), tdy(:,:,:), dy_f(:,:,:,:)
83  END TYPE array_int
84 
86  INTEGER(INTG) n ! External Structure: Function Space dimension on T_e
87  INTEGER(INTG) nl ! Internal Structure: Breakdown for tensor input
88  INTEGER(INTG) discont ! External Structure: Continuity constraint
89  INTEGER(INTG) dm ! External Structure: Field dimension
90  TYPE(array_int) i ! External Structure: Field integrations
91  REAL(DP), pointer :: q(:,:) ! External Structure: Basis coefficient tensor
92  REAL(DP), pointer :: p(:,:) ! External Structure: Basis power tensor
93  REAL(DP), pointer :: xi(:,:) ! External Structure: Basis xi-point coordinates
94  INTEGER(INTG), pointer :: b_id(:,:) ! External Structure: Basis ordering tensor
95  CHARACTER*2 cl ! External Structure: Basis call letter
96  END TYPE array_base
97 
99  TYPE(array_base), pointer :: b(:) ! Internal Structure: Basis
100  INTEGER(INTG) n_pts ! External Structure: number of volume quadrature points
101  INTEGER(INTG) n_pts_f ! External Structure: number of facet quadrature points
102  INTEGER(INTG) n_ptsl ! Internal Structure: volume quadrature for tensor input
103  INTEGER(INTG) n_pts_fl ! Internal Structure: facet quadrature for tensor input
104  INTEGER(INTG) hexa ! External Structure: tensor input parameter
105  INTEGER(INTG) faces ! External Structure: number of master element facets
106  INTEGER(INTG) fnodes ! External Structure: number of nodes per facet (spatial map)
107  INTEGER(INTG) n_b ! External Structure: number of basis
108  INTEGER(INTG) dm ! External Structure: spatial dimension of master element
109  INTEGER(INTG) spl, vel, prs
111  REAL(DP) vl ! External Structure: volume of master element
112  REAL(DP) vl_f(6) ! External Structure: facet volume of master element
113  REAL(DP) nrm(3,6) ! Internal Structure: facet normals of master element
114  REAL(DP), pointer :: gpt(:,:) ! External Structure: volume quadrature points
115  REAL(DP), pointer :: gw(:) ! External Structure: volume quadrature weights
116  REAL(DP), pointer :: gpt_f(:,:,:) ! External Structure: facet quadrature points
117  REAL(DP), pointer :: gw_f(:) ! External Structure: facet quadrature weights
118  END TYPE array_problem_base
119 
121  REAL(DP), POINTER::n(:,:)
122  INTEGER(INTG), POINTER::m(:,:),v(:,:),p(:,:)
123  INTEGER(INTG):: d,f,id_m,id_v,id_p,it_m,it_v,it_p,it_t
124  INTEGER(INTG):: e_m,e_p,e_v,e_t,en_m,en_p,en_v,en_t,n_m,n_p,n_v,n_t
125  END TYPE export_container
126 
128  !chrm, 20.08.09: For passing data of Darcy problems
129  INTEGER:: testcase
130  INTEGER:: bc_number_of_wall_nodes, number_of_bcs
131  REAL(DP):: perm, vis, perm_over_vis, p_sink
132  REAL(DP):: x1, x2, y1, y2, z1, z2
133  REAL(DP):: length, geom_tol
134  REAL(DP):: max_node_spacing
135  LOGICAL :: stab, debug_mode, analytic
136  END TYPE darcy_parameters
137 
139  INTEGER:: number_of_couplings
140  INTEGER, POINTER:: interface_element_number(:)
141  INTEGER, POINTER:: interface_element_local_node(:)
142  INTEGER:: mesh1_id
143  INTEGER, POINTER:: mesh1_element_number(:)
144  REAL(DP), POINTER:: mesh1_element_xi(:,:)
145  INTEGER:: mesh2_id
146  INTEGER, POINTER:: mesh2_element_number(:)
147  REAL(DP), POINTER:: mesh2_element_xi(:,:)
148  END TYPE coupling_parameters
149 
151  INTEGER:: number_of_bc
152  INTEGER:: number_of_bc1
153  INTEGER:: number_of_bc2
154  INTEGER:: number_of_bc3
155  INTEGER:: number_of_bc4
156  INTEGER:: number_of_bc5
157  INTEGER, POINTER:: bc_node_number1(:)
158  INTEGER, POINTER:: bc_node_number2(:)
159  INTEGER, POINTER:: bc_node_number3(:)
160  INTEGER, POINTER:: bc_node_number4(:)
161  INTEGER, POINTER:: bc_node_number5(:)
162  END TYPE boundary_parameters
163 
164 
165  !Module variables
166 
169  type(export_container), SAVE :: tmp, tmp1
170  !DLLEXPORT(TMP)
177 
180 
181 
182  INTEGER(INTG), DIMENSION(:), ALLOCATABLE:: nodesperelement
183  INTEGER(INTG), DIMENSION(:,:), ALLOCATABLE::elementnodes
184  INTEGER(INTG):: numberoffields
185  INTEGER(INTG):: numberofdimensions
186  INTEGER(INTG):: valueindex
188  INTEGER(INTG):: numberofmeshcomponents
190  INTEGER(INTG):: numberofsourcecomponents
191  INTEGER(INTG):: numberofnodesdefined
192  INTEGER(INTG):: numberoffieldcomponent(4) !(3)
193  INTEGER(INTG):: numberofelements
194  INTEGER(INTG):: globalelementnumber(10)
195  INTEGER(INTG):: maxnodesperelement
196  INTEGER(INTG):: maxnodespermeshcomponent
197  INTEGER(INTG):: element_number
198  INTEGER(INTG):: lagrange_simplex
199  INTEGER(INTG), DIMENSION(:), ALLOCATABLE:: nodespermeshcomponent
200  INTEGER(INTG), DIMENSION(:), ALLOCATABLE:: dofspermeshcomponent
201  INTEGER(INTG), DIMENSION(:), ALLOCATABLE::simplexoutputhelp,hexoutputhelp
202  INTEGER(INTG) fld, dimen, opencmiss_interpolation(3),a,b
204  INTEGER(INTG), DIMENSION(:,:), ALLOCATABLE::opencmiss_elem_m,opencmiss_elem_v,opencmiss_elem_p
206  INTEGER(INTG):: alloc_error
207  INTEGER(INTG):: field_var_type, var_idx
208 
209  INTEGER(INTG):: parameter_set_idx
210 
211  LOGICAL :: dn
212 
214 
215  REAL(DP), DIMENSION(:,:), ALLOCATABLE::elementnodesscales
216  REAL(DP), DIMENSION(:), ALLOCATABLE:: xi_coordinates,coordinates
217 ! REAL(DP):: test
218  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodexvalue
219  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeyvalue
220  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodezvalue
221  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeuvalue
222  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodevvalue
223  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodewvalue
224  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeuvalueorg
225  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodevvalueorg
226  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodewvalueorg
227  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodepvalue
228  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodepvalue2
229  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodemivalue ! Mass increase for coupled elasticity Darcy INRIA model
230  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodemuvalue
231  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodelabelvalue
232  REAL(DP), DIMENSION(:), ALLOCATABLE:: noderhovalue
233  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodea0value
234  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodea1value
235  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodea2value
236  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeh0value
237  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeevalue
238  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodekappavalue
239 
240  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeuvalue_analytic
241  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodevvalue_analytic
242  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodewvalue_analytic
243  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodepvalue_analytic
244 
245  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeuvalue_error
246  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodevvalue_error
247  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodewvalue_error
248  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodepvalue_error
249 
250  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeperm2value
251  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeperm3value
252  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeperm4value
253  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeperm5value
254  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodeperm6value
255 
256  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodesourcevalue1
257  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodesourcevalue2
258  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodesourcevalue3
259  REAL(DP), DIMENSION(:), ALLOCATABLE:: nodesourcevalue4
260 
261  REAL(DP):: scalefactorsperelementnodes(10,10)
262  REAL(DP), DIMENSION(:,:), ALLOCATABLE::opencmiss_node_coord
263 
264  CHARACTER*2 nms(99),knot
265  CHARACTER*60 in_char
266  CHARACTER*90 nimz
267 ! CHARACTER*30 NAMz
268  CHARACTER*90 namz
269 
270 
271  !Interfaces
273  MODULE PROCEDURE fluid_mechanics_io_read_cmheart1
274  MODULE PROCEDURE fluid_mechanics_io_read_cmheart2
275  MODULE PROCEDURE fluid_mechanics_io_read_cmheart3
276  END INTERFACE !FLUID_MECHANICS_IO_READ_CMHEART
277 
278 
283 
289 
290 
293  PUBLIC darcy
294 
295 CONTAINS
296 
297  ! OK
298  !================================================================================================================================
299  !
300 
302  SUBROUTINE fluid_mechanics_io_write_cmgui(REGION, EQUATIONS_SET_GLOBAL_NUMBER, NAME, ERR, ERROR,*)
304  !Argument variables
305  TYPE(region_type), POINTER :: REGION
306 ! TYPE(VARYING_STRING), INTENT(IN) :: NAME !<the prefix name of file.
307  CHARACTER(14) :: NAME
308  INTEGER(INTG) :: ERR
309  INTEGER(INTG) :: EQUATIONS_SET_GLOBAL_NUMBER
310  TYPE(varying_string):: ERROR
311  !Local Variables
312  INTEGER(INTG):: I,J,K,icompartment,FIELD_COMPONENT
313  INTEGER(INTG):: MATERIAL_INTERPOLATION_TYPE
314  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
315  TYPE(field_type), POINTER :: EQUATIONS_SET_FIELD_FIELD
316  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
317 
318  REAL(DP), POINTER :: ANALYTIC_VALUES(:)
319  LOGICAL:: ANALYTIC
320 
321  enters("FLUID_MECHANICS_IO_WRITE_CMGUI",err,error,*999)
322 
323  NULLIFY(equations_set)
324  NULLIFY(equations_set_field_field)
325  NULLIFY(equations_set_field_data)
326 
327  IF (ALLOCATED(nodesperelement)) DEALLOCATE(nodesperelement)
328  IF (ALLOCATED(nodespermeshcomponent)) DEALLOCATE(nodespermeshcomponent)
329  IF (ALLOCATED(nodespermeshcomponent)) DEALLOCATE(dofspermeshcomponent)
330  IF (ALLOCATED(xi_coordinates)) DEALLOCATE(xi_coordinates)
331  IF (ALLOCATED(coordinates)) DEALLOCATE(coordinates)
332  IF (ALLOCATED(nodexvalue)) DEALLOCATE(nodexvalue)
333  IF (ALLOCATED(nodeyvalue)) DEALLOCATE(nodeyvalue)
334  IF (ALLOCATED(nodezvalue)) DEALLOCATE(nodezvalue)
335  IF (ALLOCATED(nodeuvalue)) DEALLOCATE(nodeuvalue)
336  IF (ALLOCATED(nodevvalue)) DEALLOCATE(nodevvalue)
337  IF (ALLOCATED(nodewvalue)) DEALLOCATE(nodewvalue)
338  IF (ALLOCATED(nodeuvalueorg)) DEALLOCATE(nodeuvalueorg)
339  IF (ALLOCATED(nodevvalueorg)) DEALLOCATE(nodevvalueorg)
340  IF (ALLOCATED(nodewvalueorg)) DEALLOCATE(nodewvalueorg)
341  IF (ALLOCATED(nodepvalue)) DEALLOCATE(nodepvalue)
342  IF (ALLOCATED(nodepvalue2)) DEALLOCATE(nodepvalue2)
343  IF (ALLOCATED(nodemivalue)) DEALLOCATE(nodemivalue)
344  IF (ALLOCATED(nodemuvalue)) DEALLOCATE(nodemuvalue)
345  IF (ALLOCATED(nodelabelvalue)) DEALLOCATE(nodelabelvalue)
346  IF (ALLOCATED(noderhovalue)) DEALLOCATE(noderhovalue)
347  IF (ALLOCATED(nodea0value)) DEALLOCATE(nodea0value)
348  IF (ALLOCATED(nodea1value)) DEALLOCATE(nodea1value)
349  IF (ALLOCATED(nodea2value)) DEALLOCATE(nodea2value)
350  IF (ALLOCATED(nodeh0value)) DEALLOCATE(nodeh0value)
351  IF (ALLOCATED(nodeevalue)) DEALLOCATE(nodeevalue)
352  IF (ALLOCATED(nodekappavalue)) DEALLOCATE(nodekappavalue)
353  IF (ALLOCATED(elementnodesscales)) DEALLOCATE(elementnodesscales)
354  IF (ALLOCATED(elementnodes)) DEALLOCATE(elementnodes)
355 
356  IF (ALLOCATED(nodeperm2value)) DEALLOCATE(nodeperm2value)
357  IF (ALLOCATED(nodeperm3value)) DEALLOCATE(nodeperm3value)
358  IF (ALLOCATED(nodeperm4value)) DEALLOCATE(nodeperm4value)
359  IF (ALLOCATED(nodeperm5value)) DEALLOCATE(nodeperm5value)
360  IF (ALLOCATED(nodeperm6value)) DEALLOCATE(nodeperm6value)
361 
362  IF (ALLOCATED(nodesourcevalue1)) DEALLOCATE(nodesourcevalue1)
363  IF (ALLOCATED(nodesourcevalue2)) DEALLOCATE(nodesourcevalue2)
364  IF (ALLOCATED(nodesourcevalue3)) DEALLOCATE(nodesourcevalue3)
365  IF (ALLOCATED(nodesourcevalue4)) DEALLOCATE(nodesourcevalue4)
366 
367  !chrm, 20.08.09
368  IF (ALLOCATED(nodeuvalue_analytic)) DEALLOCATE(nodeuvalue_analytic)
369  IF (ALLOCATED(nodevvalue_analytic)) DEALLOCATE(nodevvalue_analytic)
370  IF (ALLOCATED(nodewvalue_analytic)) DEALLOCATE(nodewvalue_analytic)
371  IF (ALLOCATED(nodepvalue_analytic)) DEALLOCATE(nodepvalue_analytic)
372 
373  IF (ALLOCATED(nodeuvalue_error)) DEALLOCATE(nodeuvalue_error)
374  IF (ALLOCATED(nodevvalue_error)) DEALLOCATE(nodevvalue_error)
375  IF (ALLOCATED(nodewvalue_error)) DEALLOCATE(nodewvalue_error)
376  IF (ALLOCATED(nodepvalue_error)) DEALLOCATE(nodepvalue_error)
377 
378  knot = '0'
379  nms(1) = '1'
380  nms(2) = '2'
381  nms(3) = '3'
382  nms(4) = '4'
383  nms(5) = '5'
384  nms(6) = '6'
385  nms(7) = '7'
386  nms(8) = '8'
387  nms(9) = '9'
388 
389  k = 9
390  DO i = 1,9
391  k = k + 1
392  nms(k) = trim(nms(i))//trim(knot)
393  DO j = 1,9
394  k = k + 1
395  nms(k) = trim(nms(i))//trim(nms(j))
396  END DO
397  END DO
398 
399  equations_set => region%equations_sets%equations_sets(equations_set_global_number)%ptr
400 
401 !---tob
402 ! FIELD_VAR_TYPE=EQUATIONS_SET%EQUATIONS%EQUATIONS_MAPPING%LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE%VARIABLE_TYPE
403 ! ! '1' associated with linear matrix
404 
405  var_idx = 1
406  field_var_type = field_u_variable_type
407  SELECT CASE(equations_set%SPECIFICATION(1))
409  SELECT CASE(equations_set%SPECIFICATION(2))
411  SELECT CASE(equations_set%SPECIFICATION(3))
414  var_idx = 3
415  field_var_type = field_v_variable_type
417  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
418  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
419  & field_values_set_type,equations_set_field_data,err,error,*999)
420  icompartment=equations_set_field_data(1)
421  var_idx = 3+2*(icompartment-1)
422  field_var_type=field_v_variable_type+(field_number_of_variable_subtypes*(icompartment-1))
423  END SELECT
424  END SELECT
425  END SELECT
426 
427 
428 
429 
431  SELECT CASE(equations_set%SPECIFICATION(1))
433  SELECT CASE(equations_set%SPECIFICATION(2))
435  SELECT CASE(equations_set%SPECIFICATION(3))
438 ! parameter_set_idx = 3 ! of shared dependent variable field
439  END SELECT
440  END SELECT
441  END SELECT
442 !---toe
443 
444 ! NumberOfFields=REGION%fields%number_of_fields
445 ! Hack for ALE... to be removed later
447  numberofdimensions=region%coordinate_system%number_of_dimensions
448 ! NumberOfVariableComponents=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%dependent%dependent_field% &
449 ! & variables(1)%number_of_components
450  numberofvariablecomponents=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
451  & variables(var_idx)%number_of_components
452  NULLIFY(material_field)
454  IF(ASSOCIATED(region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials)) THEN
455  material_field=>region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field
456  IF(ASSOCIATED(material_field)) THEN
457  numberofmaterialcomponents=material_field%variables(1)%number_of_components
459  ENDIF
460  ENDIF
461  numberofelements=region%meshes%meshes(1)%ptr%number_of_elements
462  numberofmeshcomponents=region%meshes%meshes(1)%ptr%number_of_components
463 ! IF(.NOT.ALLOCATED(NodesPerElement)) ALLOCATE(NodesPerElement(NumberOfMeshComponents))
464 
465  IF(.NOT.ALLOCATED(nodesperelement)) ALLOCATE(nodesperelement(max(numberofmeshcomponents,numberofelements)))
466 
468 
471 
473  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
474  & %ptr%topology%elements%elements(1)%basis%number_of_element_parameters
475  nodespermeshcomponent(i)=region%meshes%meshes(1)%ptr%topology(i)%ptr%nodes%numberOfNodes
476  dofspermeshcomponent(i)=region%meshes%meshes(1)%ptr%topology(i)%ptr%dofs%numberOfDofs
477  END DO
478 
479 ! MaxNodesPerElement=NodesPerElement(1)
481 
482 
483  DO i=1,numberofelements
484  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
485  & %ptr%topology%elements%elements(i)%basis%number_of_element_parameters
487  END DO
488 
489 
490  IF(.NOT.ALLOCATED(xi_coordinates)) ALLOCATE(xi_coordinates(numberofdimensions))
491  IF(.NOT.ALLOCATED(coordinates)) ALLOCATE(coordinates(numberofdimensions))
492  IF(.NOT.ALLOCATED(nodexvalue)) ALLOCATE(nodexvalue(nodespermeshcomponent(1)))
493  IF(.NOT.ALLOCATED(nodeyvalue)) ALLOCATE(nodeyvalue(nodespermeshcomponent(1)))
494  IF(.NOT.ALLOCATED(nodezvalue)) ALLOCATE(nodezvalue(nodespermeshcomponent(1)))
495  IF(.NOT.ALLOCATED(nodeuvalue)) ALLOCATE(nodeuvalue(nodespermeshcomponent(1)))
496  IF(.NOT.ALLOCATED(nodevvalue)) ALLOCATE(nodevvalue(nodespermeshcomponent(1)))
497  IF(.NOT.ALLOCATED(nodewvalue)) ALLOCATE(nodewvalue(nodespermeshcomponent(1)))
498  IF(.NOT.ALLOCATED(nodeuvalueorg)) ALLOCATE(nodeuvalueorg(nodespermeshcomponent(1)))
499  IF(.NOT.ALLOCATED(nodevvalueorg)) ALLOCATE(nodevvalueorg(nodespermeshcomponent(1)))
500  IF(.NOT.ALLOCATED(nodewvalueorg)) ALLOCATE(nodewvalueorg(nodespermeshcomponent(1)))
501  IF(.NOT.ALLOCATED(nodepvalue)) ALLOCATE(nodepvalue(nodespermeshcomponent(1)))
502  IF(.NOT.ALLOCATED(nodepvalue2)) ALLOCATE(nodepvalue2(nodespermeshcomponent(1)))
503  IF(.NOT.ALLOCATED(nodemivalue)) ALLOCATE(nodemivalue(nodespermeshcomponent(1)))
504  IF(.NOT.ALLOCATED(nodemuvalue)) ALLOCATE(nodemuvalue(nodespermeshcomponent(1)))
505  IF(.NOT.ALLOCATED(nodelabelvalue)) ALLOCATE(nodelabelvalue(nodespermeshcomponent(1)))
506  IF(.NOT.ALLOCATED(noderhovalue)) ALLOCATE(noderhovalue(nodespermeshcomponent(1)))
507  IF(.NOT.ALLOCATED(nodea0value)) ALLOCATE(nodea0value(nodespermeshcomponent(1)))
508  IF(.NOT.ALLOCATED(nodea1value)) ALLOCATE(nodea1value(nodespermeshcomponent(1)))
509  IF(.NOT.ALLOCATED(nodea2value)) ALLOCATE(nodea2value(nodespermeshcomponent(1)))
510  IF(.NOT.ALLOCATED(nodeh0value)) ALLOCATE(nodeh0value(nodespermeshcomponent(1)))
511  IF(.NOT.ALLOCATED(nodeevalue)) ALLOCATE(nodeevalue(nodespermeshcomponent(1)))
512  IF(.NOT.ALLOCATED(nodekappavalue)) ALLOCATE(nodekappavalue(nodespermeshcomponent(1)))
513 ! IF(.NOT.ALLOCATED(ElementNodesScales)) ALLOCATE(ElementNodesScales(NumberOfElements,NodesPerElement(1)))
515 ! IF(.NOT.ALLOCATED(ElementNodes)) ALLOCATE(ElementNodes(NumberOfElements,NodesPerElement(1)))
516  IF(.NOT.ALLOCATED(elementnodes)) ALLOCATE(elementnodes(numberofelements,maxnodesperelement))
517 
518  IF(.NOT.ALLOCATED(nodeperm2value)) ALLOCATE(nodeperm2value(nodespermeshcomponent(1)))
519  IF(.NOT.ALLOCATED(nodeperm3value)) ALLOCATE(nodeperm3value(nodespermeshcomponent(1)))
520  IF(.NOT.ALLOCATED(nodeperm4value)) ALLOCATE(nodeperm4value(nodespermeshcomponent(1)))
521  IF(.NOT.ALLOCATED(nodeperm5value)) ALLOCATE(nodeperm5value(nodespermeshcomponent(1)))
522  IF(.NOT.ALLOCATED(nodeperm6value)) ALLOCATE(nodeperm6value(nodespermeshcomponent(1)))
523 
524  IF(.NOT.ALLOCATED(nodesourcevalue1)) ALLOCATE(nodesourcevalue1(nodespermeshcomponent(1)))
525  IF(.NOT.ALLOCATED(nodesourcevalue2)) ALLOCATE(nodesourcevalue2(nodespermeshcomponent(1)))
526  IF(.NOT.ALLOCATED(nodesourcevalue3)) ALLOCATE(nodesourcevalue3(nodespermeshcomponent(1)))
527  IF(.NOT.ALLOCATED(nodesourcevalue4)) ALLOCATE(nodesourcevalue4(nodespermeshcomponent(1)))
528 
529  !chrm, 20.08.09
530  IF(.NOT.ALLOCATED(nodeuvalue_analytic)) ALLOCATE(nodeuvalue_analytic(nodespermeshcomponent(1)))
531  IF(.NOT.ALLOCATED(nodevvalue_analytic)) ALLOCATE(nodevvalue_analytic(nodespermeshcomponent(1)))
532  IF(.NOT.ALLOCATED(nodewvalue_analytic)) ALLOCATE(nodewvalue_analytic(nodespermeshcomponent(1)))
533  IF(.NOT.ALLOCATED(nodepvalue_analytic)) ALLOCATE(nodepvalue_analytic(nodespermeshcomponent(1)))
534 
535  IF(.NOT.ALLOCATED(nodeuvalue_error)) ALLOCATE(nodeuvalue_error(nodespermeshcomponent(1)))
536  IF(.NOT.ALLOCATED(nodevvalue_error)) ALLOCATE(nodevvalue_error(nodespermeshcomponent(1)))
537  IF(.NOT.ALLOCATED(nodewvalue_error)) ALLOCATE(nodewvalue_error(nodespermeshcomponent(1)))
538  IF(.NOT.ALLOCATED(nodepvalue_error)) ALLOCATE(nodepvalue_error(nodespermeshcomponent(1)))
539 
540  enters("CMGUI OUTPUT",err,error,*999)
541 
542  NULLIFY(source_field)
543 
544  field=>region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field
545 
546  !material field
547  output_material = ASSOCIATED(material_field)
548 
549  !source field
550  output_source = .false.
551  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
552  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) &
553  & .AND.(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) )THEN
554  source_field=>region%equations_sets%equations_sets(equations_set_global_number)%ptr%source%source_field
555  IF( ASSOCIATED(source_field) ) output_source = .true.
556  END IF
557 
561 
562  CALL field_interpolation_parameters_initialise(field,interpolation_parameters,err,error,*999)
563  CALL field_interpolated_points_initialise(interpolation_parameters,interpolated_point,err,error,*999)
564 
565  IF(output_material) THEN
566  !material
567  CALL field_interpolation_parameters_initialise(material_field,material_interpolation_parameters,err,error,*999)
568  CALL field_interpolated_points_initialise(material_interpolation_parameters,material_interpolated_point,err,error,*999)
569  ENDIF
570 
571  !source field
572  IF( output_source ) THEN
573  CALL field_interpolation_parameters_initialise(source_field,source_interpolation_parameters,err,error,*999)
574  CALL field_interpolated_points_initialise(source_interpolation_parameters,source_interpolated_point,err,error,*999)
575  numberofsourcecomponents=region%equations_sets%equations_sets(equations_set_global_number)%ptr%source%source_field% &
576  & variables(1)%number_of_components
577 
579  END IF
580 
581  !analytic field
582  analytic=.false.
583  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
584  analytic=.true.
585  CALL field_interpolation_parameters_initialise(field,analytic_interpolation_parameters,err,error,*999)
586  CALL field_interpolated_points_initialise(analytic_interpolation_parameters,analytic_interpolated_point,err,error,*999)
587  CALL field_parameter_set_data_get(field,field_u_variable_type,field_analytic_values_set_type,analytic_values, &
588  & err,error,*999)
589  ENDIF
590 
591  DO i=1,numberofelements
592  numberofdimensions=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
593  & %ptr%topology%elements%elements(i)%basis%number_of_xi
594  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
595  & %ptr%topology%elements%elements(i)%basis%number_of_element_parameters
596 
597 ! DO J=1,NodesPerElement(1)
598  DO j=1,nodesperelement(i)
599 
601 ! XI_COORDINATES(1)=(REGION%equations_sets%equations_sets(1)%ptr%equations%interpolation% &
602 ! & geometric_interp_parameters%bases(1)%ptr%node_position_index(J,1)-1.0)/(REGION%equations_sets% &
603 ! & equations_sets(1)%ptr%equations%interpolation%geometric_interp_parameters%bases(1) &
604 ! & %ptr%number_of_nodes_xi(1)-1.0)
605 ! XI_COORDINATES(2)=(REGION%equations_sets%equations_sets(1)%ptr%equations%interpolation% &
606 ! & geometric_interp_parameters%bases(1)%ptr%node_position_index(J,2)-1.0)/(REGION%equations_sets% &
607 ! & equations_sets(1)%ptr%equations%interpolation%geometric_interp_parameters%bases(1) &
608 ! & %ptr%number_of_nodes_xi(2)-1.0)
609  IF(equations_set%SPECIFICATION(3)==equations_set_transient1d_navier_stokes_subtype.OR. &
610  & equations_set%SPECIFICATION(3)==equations_set_coupled1d0d_navier_stokes_subtype)THEN
611  xi_coordinates(1)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
612  & geometric_field%variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(j))
613  ELSE
614  xi_coordinates(1)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
615  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%node_position_index(j,1)-1.0)/(region% &
616  & equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
617  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%number_of_nodes_xic(1)-1.0)
618  ENDIF
619  IF(numberofdimensions==2 .OR. numberofdimensions==3)THEN
620  xi_coordinates(2)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
621  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%node_position_index(j,2)-1.0)/(region% &
622  & equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
623  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%number_of_nodes_xic(2)-1.0)
624  ! XI_COORDINATES(2)=(REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
625  ! & geometric_field%variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(J))
626  END IF
627  IF(numberofdimensions==3)THEN
628  xi_coordinates(3)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
629  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%node_position_index(j,3)-1.0)/(region% &
630  & equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
631  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%number_of_nodes_xic(3)-1.0)
632  END IF
633 
634 !Start: This is a hack for 3D simplex elements
635  IF(numberofdimensions==2)THEN
636  IF (nodesperelement(1)==3) THEN
637  IF(j==1) xi_coordinates(1:2)=[0.0_dp,1.0_dp]
638  IF(j==2) xi_coordinates(1:2)=[1.0_dp,0.0_dp]
639  IF(j==3) xi_coordinates(1:2)=[1.0_dp,1.0_dp]
640  ELSE IF (nodesperelement(1)==6) THEN
641  IF(j==1) xi_coordinates(1:2)=[0.0_dp,1.0_dp]
642  IF(j==2) xi_coordinates(1:2)=[1.0_dp,0.0_dp]
643  IF(j==3) xi_coordinates(1:2)=[1.0_dp,1.0_dp]
644  IF(j==4) xi_coordinates(1:2)=[0.5_dp,0.5_dp]
645  IF(j==5) xi_coordinates(1:2)=[1.0_dp,0.5_dp]
646  IF(j==6) xi_coordinates(1:2)=[0.5_dp,1.0_dp]
647  ELSE IF (nodesperelement(1)==10) THEN
648  IF(j==1) xi_coordinates(1:2)=[0.0_dp,1.0_dp]
649  IF(j==2) xi_coordinates(1:2)=[1.0_dp,0.0_dp]
650  IF(j==3) xi_coordinates(1:2)=[1.0_dp,1.0_dp]
651  IF(j==4) xi_coordinates(1:2)=[1.0_dp/3.0_dp,2.0_dp/3.0_dp]
652  IF(j==5) xi_coordinates(1:2)=[2.0_dp/3.0_dp,1.0_dp/3.0_dp]
653  IF(j==6) xi_coordinates(1:2)=[1.0_dp,1.0_dp/3.0_dp]
654  IF(j==7) xi_coordinates(1:2)=[1.0_dp,2.0_dp/3.0_dp]
655  IF(j==8) xi_coordinates(1:2)=[2.0_dp/3.0_dp,1.0_dp]
656  IF(j==9) xi_coordinates(1:2)=[1.0_dp/3.0_dp,1.0_dp]
657  IF(j==10) xi_coordinates(1:2)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp]
658  ENDIF
659  ELSE IF(numberofdimensions==3)THEN
660  IF (nodesperelement(1)==4) THEN
661  IF(j==1) xi_coordinates(1:3)=[0.0_dp,1.0_dp,1.0_dp]
662  IF(j==2) xi_coordinates(1:3)=[1.0_dp,0.0_dp,1.0_dp]
663  IF(j==3) xi_coordinates(1:3)=[1.0_dp,1.0_dp,0.0_dp]
664  IF(j==4) xi_coordinates(1:3)=[1.0_dp,1.0_dp,1.0_dp]
665  ELSE IF (nodesperelement(1)==10) THEN
666  IF(j==1) xi_coordinates(1:3)=[0.0_dp,1.0_dp,1.0_dp]
667  IF(j==2) xi_coordinates(1:3)=[1.0_dp,0.0_dp,1.0_dp]
668  IF(j==3) xi_coordinates(1:3)=[1.0_dp,1.0_dp,0.0_dp]
669  IF(j==4) xi_coordinates(1:3)=[1.0_dp,1.0_dp,1.0_dp]
670  IF(j==5) xi_coordinates(1:3)=[0.5_dp,0.5_dp,1.0_dp]
671  IF(j==6) xi_coordinates(1:3)=[0.5_dp,1.0_dp,0.5_dp]
672  IF(j==7) xi_coordinates(1:3)=[0.5_dp,1.0_dp,1.0_dp]
673  IF(j==8) xi_coordinates(1:3)=[1.0_dp,0.5_dp,0.5_dp]
674  IF(j==9) xi_coordinates(1:3)=[1.0_dp,1.0_dp,0.5_dp]
675  IF(j==10) xi_coordinates(1:3)=[1.0_dp,0.5_dp,1.0_dp]
676  ELSE IF (nodesperelement(1)==20) THEN
677  IF(j==1) xi_coordinates(1:3)=[0.0_dp,1.0_dp,1.0_dp]
678  IF(j==2) xi_coordinates(1:3)=[1.0_dp,0.0_dp,1.0_dp]
679  IF(j==3) xi_coordinates(1:3)=[1.0_dp,1.0_dp,0.0_dp]
680  IF(j==4) xi_coordinates(1:3)=[1.0_dp,1.0_dp,1.0_dp]
681  IF(j==5) xi_coordinates(1:3)=[1.0_dp/3.0_dp,2.0_dp/3.0_dp,1.0_dp]
682  IF(j==6) xi_coordinates(1:3)=[2.0_dp/3.0_dp,1.0_dp/3.0_dp,1.0_dp]
683  IF(j==7) xi_coordinates(1:3)=[1.0_dp/3.0_dp,1.0_dp,2.0_dp/3.0_dp]
684  IF(j==8) xi_coordinates(1:3)=[2.0_dp/3.0_dp,1.0_dp,1.0_dp/3.0_dp]
685  IF(j==9) xi_coordinates(1:3)=[1.0_dp/3.0_dp,1.0_dp,1.0_dp]
686  IF(j==10) xi_coordinates(1:3)=[2.0_dp/3.0_dp,1.0_dp,1.0_dp]
687  IF(j==11) xi_coordinates(1:3)=[1.0_dp,1.0_dp/3.0_dp,2.0_dp/3.0_dp]
688  IF(j==12) xi_coordinates(1:3)=[1.0_dp,2.0_dp/3.0_dp,1.0_dp/3.0_dp]
689  IF(j==13) xi_coordinates(1:3)=[1.0_dp,1.0_dp,1.0_dp/3.0_dp]
690  IF(j==14) xi_coordinates(1:3)=[1.0_dp,1.0_dp,2.0_dp/3.0_dp]
691  IF(j==15) xi_coordinates(1:3)=[1.0_dp,1.0_dp/3.0_dp,1.0_dp]
692  IF(j==16) xi_coordinates(1:3)=[1.0_dp,2.0_dp/3.0_dp,1.0_dp]
693  IF(j==17) xi_coordinates(1:3)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp,2.0_dp/3.0_dp]
694  IF(j==18) xi_coordinates(1:3)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp,1.0_dp]
695  IF(j==19) xi_coordinates(1:3)=[2.0_dp/3.0_dp,1.0_dp,2.0_dp/3.0_dp]
696  IF(j==20) xi_coordinates(1:3)=[1.0_dp,2.0_dp/3.0_dp,2.0_dp/3.0_dp]
697  ENDIF
698  ENDIF
699 
700 !End: This is a hack for 3D simplex elements
701 
702  !K is global node number
703  k=region%meshes%meshes(1)%ptr%topology(1)%ptr%elements%elements(i)%global_element_nodes(j)
704 
705  IF(numberofdimensions==3)THEN
706  coordinates(1:3)=[1,1,1]
707  ELSE IF(numberofdimensions==2)THEN
708  coordinates(1:2)=[1,1]
709  END IF
710 
711  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
712  & interpolation_parameters(field_var_type)%ptr,err,error,*999)
713  CALL field_interpolate_xi(no_part_deriv,xi_coordinates,interpolated_point(field_var_type)%ptr,err,error,*999)
714 
715  !analytic field
716  IF( analytic ) THEN
717  CALL field_interpolation_parameters_element_get(field_analytic_values_set_type,element_number, &
718  & analytic_interpolation_parameters(field_u_variable_type)%ptr,err,error,*999)
719  CALL field_interpolate_xi(no_part_deriv,xi_coordinates,analytic_interpolated_point(field_u_variable_type)%ptr, &
720  & err,error,*999)
721  END IF
722 
723  !material
724  IF( output_material ) THEN
725  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
726  & material_interpolation_parameters(field_u_variable_type)%ptr,err,error,*999)
727  CALL field_interpolate_xi(no_part_deriv,xi_coordinates,material_interpolated_point(field_u_variable_type)%ptr, &
728  & err,error,*999)
729  ENDIF
730 
731  !source field
732  IF( output_source ) THEN
733  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
734  & source_interpolation_parameters(field_u_variable_type)%ptr,err,error,*999)
735  CALL field_interpolate_xi(no_part_deriv,xi_coordinates,source_interpolated_point(field_u_variable_type)%ptr, &
736  & err,error,*999)
737  END IF
738 
739 
740  nodexvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field%variables(1) &
741  & %parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k)
742 ! NodeYValue(K)=REGION%equations_sets%equations_sets(1)%ptr%geometry%geometric_field%variables(1) &
743 ! & %parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(K+NodesPerMeshComponent(1))
744  IF(numberofdimensions==2 .OR. numberofdimensions==3)THEN
745  nodeyvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field% &
746  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+nodespermeshcomponent(1))
747  END IF
748  IF(numberofdimensions==3)THEN
749  nodezvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field% &
750  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+2*nodespermeshcomponent(1))
751  END IF
752 
753  !chrm, 19.12.2010
754  !If geometry uses a quadratic mesh, but Darcy velocity and mass increase only a linear one,
755  ! then we do need to interpolate the velocities and mass increase on the geometric mesh
756  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
757  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) ) THEN
758  IF ((equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) .OR. &
759  & (equations_set%SPECIFICATION(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype)) THEN
760  nodeuvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(1,1)
761  IF(numberofdimensions==2 .OR. numberofdimensions==3)THEN
762  nodevvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(2,1)
763  ENDIF
764  IF(numberofdimensions==3)THEN
765  nodewvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(3,1)
766  ENDIF
767  ENDIF
768  ELSE
769  nodeuvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
770  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
771  IF(equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class .OR. &
772  & equations_set%SPECIFICATION(1)==equations_set_elasticity_class) THEN
773  nodevvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
774  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters% &
775  & cmiss%data_dp(k+nodespermeshcomponent(1))
776 
777  IF(numberofdimensions==3)THEN
778  nodewvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
779  & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters% &
780  & cmiss%data_dp(k+2*nodespermeshcomponent(1))
781  END IF
782  END IF
783  END IF
784 
785 ! ! ! NodeUValue(K)=INTERPOLATED_POINT(FIELD_VAR_TYPE)%ptr%VALUES(1,1)
786 ! ! ! NodeVValue(K)=INTERPOLATED_POINT(FIELD_VAR_TYPE)%ptr%VALUES(2,1)
787 ! ! ! NodeWValue(K)=INTERPOLATED_POINT(FIELD_VAR_TYPE)%ptr%VALUES(3,1)
788 
789  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
790  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
791  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
792  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
793  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_constitutive_and_growth_law_in_cellml_subtype) &
794  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype) )THEN
795  IF(numberofdimensions==3)THEN
796  nodepvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(4,1)
797  ELSE IF(numberofdimensions==2)THEN
798  nodepvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(3,1)
799  END IF
800  END IF
801 
802  IF(equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) THEN
803  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
804  nodeuvalue_analytic(k)=analytic_values(k)
805  IF(numberofdimensions==2 .OR. numberofdimensions==3)THEN
806  nodevvalue_analytic(k)=analytic_values(k+nodespermeshcomponent(1))
807  IF(numberofdimensions==3)THEN
808  nodewvalue_analytic(k)=analytic_values(k+2*nodespermeshcomponent(1))
809  ENDIF
810  ENDIF
811  SELECT CASE(numberofdimensions)
812  CASE(1)
814  CASE(2)
816  CASE(3)
818  END SELECT
819  ENDIF
820  ENDIF
821 
822 
823 !---tob: Mass increase for coupled elasticity Darcy INRIA model
824  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
825  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) &
826  & .AND.(equations_set%SPECIFICATION(3)==equations_set_elasticity_darcy_inria_model_subtype) )THEN
827  IF(numberofdimensions==3)THEN
828  nodemivalue(k)=interpolated_point(field_var_type)%ptr%VALUES(5,1)
829  END IF
830  END IF
831 !---toe
832 
833  IF(output_material) THEN
834  material_interpolation_type=material_field%variables(1)%COMPONENTS(1)%INTERPOLATION_TYPE
835 
836  IF(material_interpolation_type==field_node_based_interpolation)THEN
837 
838  IF( (equations_set%specification(1)==equations_set_fluid_mechanics_class) &
839  & .AND.(equations_set%specification(2)==equations_set_darcy_equation_type) &
840  & .AND.(equations_set%specification(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) )THEN
841  nodemuvalue(k) =material_interpolated_point(field_u_variable_type)%ptr%VALUES(1,1)
842  noderhovalue(k) =material_interpolated_point(field_u_variable_type)%ptr%VALUES(2,1)
843  ! !--- Remaining tensor material data of permeability tensor
844  nodeperm2value(k)=material_interpolated_point(field_u_variable_type)%ptr%VALUES(3,1)
845  nodeperm3value(k)=material_interpolated_point(field_u_variable_type)%ptr%VALUES(4,1)
846  nodeperm4value(k)=material_interpolated_point(field_u_variable_type)%ptr%VALUES(5,1)
847  nodeperm5value(k)=material_interpolated_point(field_u_variable_type)%ptr%VALUES(6,1)
848  nodeperm6value(k)=material_interpolated_point(field_u_variable_type)%ptr%VALUES(7,1)
849  ELSE
850  nodemuvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
851  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k)
852  noderhovalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
853  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+nodespermeshcomponent(1))
854  END IF
855 
856  IF(equations_set%specification(1)==equations_set_elasticity_class)THEN
857  IF(equations_set%specification(2)==equations_set_finite_elasticity_type)THEN
858  IF( (equations_set%specification(3)==equations_set_compressible_finite_elasticity_subtype) &
859  & .OR.(equations_set%specification(3)==equations_set_elasticity_darcy_inria_model_subtype) &
860  & .OR.(equations_set%specification(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) &
861  & .OR. (equations_set%specification(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) )THEN
862  nodekappavalue(k)=material_field%variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters% &
863  & cmiss%data_dp(k+2*nodespermeshcomponent(1))
864  END IF
865  END IF
866  END IF
867 
868  ELSE !default to FIELD_CONSTANT_INTERPOLATION
869  nodemuvalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
870  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(1)
871  noderhovalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
872  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(2)
873 
874  IF(equations_set%specification(3)==equations_set_transient1d_navier_stokes_subtype) THEN
875  nodeevalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
876  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(3)
877  nodeh0value=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
878  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(6)
879  nodea0value=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
880  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(9)
881  nodea1value=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
882  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(12)!39
883  nodea2value=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
884  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(15)!51
885  END IF
886 
887  IF(equations_set%specification(1)==equations_set_elasticity_class)THEN
888  IF(equations_set%specification(2)==equations_set_finite_elasticity_type)THEN
889  IF( (equations_set%specification(3)==equations_set_compressible_finite_elasticity_subtype) &
890  & .OR.(equations_set%specification(3)==equations_set_elasticity_darcy_inria_model_subtype) &
891  & .OR.(equations_set%specification(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) &
892  & .OR. (equations_set%specification(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) )THEN
893  nodekappavalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
894  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(3)
895  END IF
896  END IF
897  END IF
898  END IF
899  ENDIF !oUTPUT MATERIAL
900  !source field
901  IF( output_source ) THEN
902  IF( (equations_set%specification(1)==equations_set_fluid_mechanics_class) &
903  & .AND.(equations_set%specification(2)==equations_set_darcy_equation_type) &
904  & .AND.(equations_set%specification(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) )THEN
905  IF(source_field%variables(1)%COMPONENTS(1)%INTERPOLATION_TYPE==field_node_based_interpolation)THEN
906  nodesourcevalue1(k)=source_interpolated_point(field_u_variable_type)%ptr%VALUES(1,1)
907  nodesourcevalue2(k)=source_interpolated_point(field_u_variable_type)%ptr%VALUES(2,1)
908  nodesourcevalue3(k)=source_interpolated_point(field_u_variable_type)%ptr%VALUES(3,1)
909  nodesourcevalue4(k)=source_interpolated_point(field_u_variable_type)%ptr%VALUES(4,1)
910  END IF
911  END IF
912  ENDIF
913  END DO
914  END DO
915 
916  dn=.true.
917  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
918  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) &
919  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) )THEN
920  dn=.false.
921  END IF
922  IF(equations_set%SPECIFICATION(1)==equations_set_classical_field_class)THEN
923  dn=.false.
924  END IF
925 
926  IF(numberofdimensions==2 .OR. numberofdimensions==3)THEN
927  IF(dn) THEN
928  ! output for DN only
929  DO k=1,nodespermeshcomponent(2)
930  IF(numberofdimensions==2)THEN
931  nodepvalue2(k)=region%equations_sets%equations_sets(1)%ptr%dependent%dependent_field%variables(1) &
932  & %parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(2*nodespermeshcomponent(1)+k)
933  ELSE IF(numberofdimensions==3)THEN
934  nodepvalue2(k)=region%equations_sets%equations_sets(1)%ptr%dependent%dependent_field%variables(1) &
935  & %parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(3*nodespermeshcomponent(1)+k)
936  ENDIF
937  ENDDO
938  ENDIF
939  ENDIF
940 
941  IF( numberofdimensions==3 )THEN
942  !For 3D, the following call works ...
943  lagrange_simplex=region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations% &
944  & interpolation%geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%type
945 ! lagrange_simplex=2
946 
947  ELSE IF (numberofdimensions==1) THEN
949  ELSE IF (numberofdimensions==2) THEN
950  !chrm, 20.08.09:
951  ! ... but the above call does not work for 2D.
952  !Thus, for 2D, we hard-wire it to 'quad':
955  ELSE IF(maxnodesperelement==3.OR.maxnodesperelement==6.OR.maxnodesperelement==10) THEN
957  ENDIF
958  END IF
959 
962  field_component=2
963  IF(output_material) THEN
964  field_component=field_component+1
966  ENDIF
967 
968  IF( output_source ) THEN
969  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
970  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) &
971  & .AND.(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) )THEN
972  field_component=field_component+1
974  END IF
975  END IF
976 
977  DO i=1,numberofelements
978 ! DO J=1,NodesPerElement(1)
979  DO j=1,nodesperelement(i)
980  elementnodes(i,j)=region%meshes%meshes(1)%ptr%topology(1)% &
981  & ptr%elements%elements(i)%global_element_nodes(j)
982  elementnodesscales(i,j)=1.0000000000000000e+00
983  END DO
984  END DO
985 
986  CALL fluid_mechanics_io_write_nodes_cmgui(name,equations_set)
988 
989  exits("FLUID_MECHANICS_IO_WRITE_CMGUI")
990  RETURN
991 999 errorsexits("FLUID_MECHANICS_IO_WRITE_CMGUI",err,error)
992  RETURN
993 
994  END SUBROUTINE fluid_mechanics_io_write_cmgui
995 
996 
997  ! OK
998  !================================================================================================================================
999  !
1000 
1002 ! SUBROUTINE FLUID_MECHANICS_IO_WRITE_ENCAS(REGION, NAME, ERR, ERROR,*)
1003  SUBROUTINE fluid_mechanics_io_write_encas(REGION, EQUATIONS_SET_GLOBAL_NUMBER, NAME, ERR, ERROR,*)
1005  !Argument variables
1006  TYPE(region_type), POINTER :: REGION
1007 ! TYPE(VARYING_STRING), INTENT(IN) :: NAME !<the prefix name of file.
1008  CHARACTER(14) :: NAME
1009  INTEGER(INTG) :: EQUATIONS_SET_GLOBAL_NUMBER
1010  INTEGER(INTG) :: ERR
1011  TYPE(varying_string):: ERROR
1012  !Local Variables
1013  INTEGER(INTG):: I,J,K,icompartment
1014  INTEGER(INTG):: MATERIAL_INTERPOLATION_TYPE
1015  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1016  TYPE(field_type), POINTER :: EQUATIONS_SET_FIELD_FIELD
1017  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
1018 
1019  enters("FLUID_MECHANICS_IO_WRITE_ENCAS",err,error,*999)
1020 
1021  IF (ALLOCATED(nodesperelement)) DEALLOCATE(nodesperelement)
1022  IF (ALLOCATED(nodespermeshcomponent)) DEALLOCATE(nodespermeshcomponent)
1023  IF (ALLOCATED(xi_coordinates)) DEALLOCATE(xi_coordinates)
1024  IF (ALLOCATED(coordinates)) DEALLOCATE(coordinates)
1025  IF (ALLOCATED(nodexvalue)) DEALLOCATE(nodexvalue)
1026  IF (ALLOCATED(nodeyvalue)) DEALLOCATE(nodeyvalue)
1027  IF (ALLOCATED(nodezvalue)) DEALLOCATE(nodezvalue)
1028  IF (ALLOCATED(nodeuvalue)) DEALLOCATE(nodeuvalue)
1029  IF (ALLOCATED(nodevvalue)) DEALLOCATE(nodevvalue)
1030  IF (ALLOCATED(nodewvalue)) DEALLOCATE(nodewvalue)
1031  IF (ALLOCATED(nodeuvalueorg)) DEALLOCATE(nodeuvalueorg)
1032  IF (ALLOCATED(nodevvalueorg)) DEALLOCATE(nodevvalueorg)
1033  IF (ALLOCATED(nodewvalueorg)) DEALLOCATE(nodewvalueorg)
1034  IF (ALLOCATED(nodepvalue)) DEALLOCATE(nodepvalue)
1035  IF (ALLOCATED(nodepvalue2)) DEALLOCATE(nodepvalue2)
1036  IF (ALLOCATED(nodemivalue)) DEALLOCATE(nodemivalue)
1037  IF (ALLOCATED(nodemuvalue)) DEALLOCATE(nodemuvalue)
1038  IF (ALLOCATED(nodelabelvalue)) DEALLOCATE(nodelabelvalue)
1039  IF (ALLOCATED(noderhovalue)) DEALLOCATE(noderhovalue)
1040  IF (ALLOCATED(nodea0value)) DEALLOCATE(nodea0value)
1041  IF (ALLOCATED(nodea1value)) DEALLOCATE(nodea1value)
1042  IF (ALLOCATED(nodea2value)) DEALLOCATE(nodea2value)
1043  IF (ALLOCATED(nodeh0value)) DEALLOCATE(nodeh0value)
1044  IF (ALLOCATED(nodeevalue)) DEALLOCATE(nodeevalue)
1045  IF (ALLOCATED(nodekappavalue)) DEALLOCATE(nodekappavalue)
1046  IF (ALLOCATED(elementnodesscales)) DEALLOCATE(elementnodesscales)
1047  IF (ALLOCATED(elementnodes)) DEALLOCATE(elementnodes)
1048 
1049  !chrm, 20.08.09
1050  IF (ALLOCATED(nodeuvalue_analytic)) DEALLOCATE(nodeuvalue_analytic)
1051  IF (ALLOCATED(nodevvalue_analytic)) DEALLOCATE(nodevvalue_analytic)
1052  IF (ALLOCATED(nodewvalue_analytic)) DEALLOCATE(nodewvalue_analytic)
1053  IF (ALLOCATED(nodepvalue_analytic)) DEALLOCATE(nodepvalue_analytic)
1054 
1055  IF (ALLOCATED(nodeuvalue_error)) DEALLOCATE(nodeuvalue_error)
1056  IF (ALLOCATED(nodevvalue_error)) DEALLOCATE(nodevvalue_error)
1057  IF (ALLOCATED(nodewvalue_error)) DEALLOCATE(nodewvalue_error)
1058  IF (ALLOCATED(nodepvalue_error)) DEALLOCATE(nodepvalue_error)
1059 
1060  NULLIFY(equations_set_field_data)
1061 
1062  knot = '0'
1063  nms(1) = '1'
1064  nms(2) = '2'
1065  nms(3) = '3'
1066  nms(4) = '4'
1067  nms(5) = '5'
1068  nms(6) = '6'
1069  nms(7) = '7'
1070  nms(8) = '8'
1071  nms(9) = '9'
1072 
1073  k = 9
1074  DO i = 1,9
1075  k = k + 1
1076  nms(k) = trim(nms(i))//trim(knot)
1077  DO j = 1,9
1078  k = k + 1
1079  nms(k) = trim(nms(i))//trim(nms(j))
1080  END DO
1081  END DO
1082 
1083  equations_set => region%equations_sets%equations_sets(equations_set_global_number)%ptr
1084 
1085 !---tob
1086 ! FIELD_VAR_TYPE=EQUATIONS_SET%EQUATIONS%EQUATIONS_MAPPING%LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE%VARIABLE_TYPE
1087 ! ! '1' associated with linear matrix
1088 
1089  var_idx = 1
1090  field_var_type = field_u_variable_type
1091  SELECT CASE(equations_set%SPECIFICATION(1))
1093  SELECT CASE(equations_set%SPECIFICATION(2))
1095  SELECT CASE(equations_set%SPECIFICATION(3))
1098  var_idx = 3
1099  field_var_type = field_v_variable_type
1101  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1102  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1103  & field_values_set_type,equations_set_field_data,err,error,*999)
1104  icompartment=equations_set_field_data(1)
1105  var_idx = 3+2*(icompartment-1)
1106  field_var_type=field_v_variable_type+(field_number_of_variable_subtypes*(icompartment-1))
1107  END SELECT
1108  END SELECT
1109  END SELECT
1110 
1111  parameter_set_idx = 1
1112  SELECT CASE(equations_set%SPECIFICATION(1))
1114  SELECT CASE(equations_set%SPECIFICATION(2))
1116  SELECT CASE(equations_set%SPECIFICATION(3))
1119 ! parameter_set_idx = 3 ! of shared dependent variable field
1120  END SELECT
1121  END SELECT
1122  END SELECT
1123 !---toe
1124 
1125 ! NumberOfFields=REGION%fields%number_of_fields
1126 ! Hack for ALE... to be removed later
1127  numberoffields=3
1128  numberofdimensions=region%coordinate_system%number_of_dimensions
1129 ! NumberOfVariableComponents=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%dependent%dependent_field% &
1130 ! & variables(1)%number_of_components
1131  numberofvariablecomponents=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1132  & variables(var_idx)%number_of_components
1133 
1134  numberofmaterialcomponents=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1135  & variables(1)%number_of_components
1136  numberofelements=region%meshes%meshes(1)%ptr%number_of_elements
1137  numberofmeshcomponents=region%meshes%meshes(1)%ptr%number_of_components
1138 ! IF(.NOT.ALLOCATED(NodesPerElement)) ALLOCATE(NodesPerElement(NumberOfMeshComponents))
1139 
1140  IF(.NOT.ALLOCATED(nodesperelement)) ALLOCATE(nodesperelement(max(numberofmeshcomponents,numberofelements)))
1141 
1144 
1145  DO i=1,numberofmeshcomponents
1146  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
1147  & %ptr%topology%elements%elements(1)%basis%number_of_element_parameters
1148  nodespermeshcomponent(i)=region%meshes%meshes(1)%ptr%topology(i)%ptr%nodes%numberOfNodes
1149  END DO
1150 
1151 
1152 ! MaxNodesPerElement=NodesPerElement(1)
1154 
1155 
1156  DO i=1,numberofelements
1157  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
1158  & %ptr%topology%elements%elements(i)%basis%number_of_element_parameters
1160  END DO
1161 
1162 
1163  IF(.NOT.ALLOCATED(xi_coordinates)) ALLOCATE(xi_coordinates(numberofdimensions))
1164  IF(.NOT.ALLOCATED(coordinates)) ALLOCATE(coordinates(numberofdimensions))
1165  IF(.NOT.ALLOCATED(nodexvalue)) ALLOCATE(nodexvalue(nodespermeshcomponent(1)))
1166  IF(.NOT.ALLOCATED(nodeyvalue)) ALLOCATE(nodeyvalue(nodespermeshcomponent(1)))
1167  IF(.NOT.ALLOCATED(nodezvalue)) ALLOCATE(nodezvalue(nodespermeshcomponent(1)))
1168  IF(.NOT.ALLOCATED(nodeuvalue)) ALLOCATE(nodeuvalue(nodespermeshcomponent(1)))
1169  IF(.NOT.ALLOCATED(nodevvalue)) ALLOCATE(nodevvalue(nodespermeshcomponent(1)))
1170  IF(.NOT.ALLOCATED(nodewvalue)) ALLOCATE(nodewvalue(nodespermeshcomponent(1)))
1171  IF(.NOT.ALLOCATED(nodeuvalueorg)) ALLOCATE(nodeuvalueorg(nodespermeshcomponent(1)))
1172  IF(.NOT.ALLOCATED(nodevvalueorg)) ALLOCATE(nodevvalueorg(nodespermeshcomponent(1)))
1173  IF(.NOT.ALLOCATED(nodewvalueorg)) ALLOCATE(nodewvalueorg(nodespermeshcomponent(1)))
1174  IF(.NOT.ALLOCATED(nodepvalue)) ALLOCATE(nodepvalue(nodespermeshcomponent(1)))
1175  IF(.NOT.ALLOCATED(nodepvalue2)) ALLOCATE(nodepvalue2(nodespermeshcomponent(1)))
1176  IF(.NOT.ALLOCATED(nodemivalue)) ALLOCATE(nodemivalue(nodespermeshcomponent(1)))
1177  IF(.NOT.ALLOCATED(nodemuvalue)) ALLOCATE(nodemuvalue(nodespermeshcomponent(1)))
1178  IF(.NOT.ALLOCATED(nodelabelvalue)) ALLOCATE(nodelabelvalue(nodespermeshcomponent(1)))
1179  IF(.NOT.ALLOCATED(noderhovalue)) ALLOCATE(noderhovalue(nodespermeshcomponent(1)))
1180  IF(.NOT.ALLOCATED(nodea0value)) ALLOCATE(nodea0value(nodespermeshcomponent(1)))
1181  IF(.NOT.ALLOCATED(nodea1value)) ALLOCATE(nodea1value(nodespermeshcomponent(1)))
1182  IF(.NOT.ALLOCATED(nodea2value)) ALLOCATE(nodea2value(nodespermeshcomponent(1)))
1183  IF(.NOT.ALLOCATED(nodeh0value)) ALLOCATE(nodeh0value(nodespermeshcomponent(1)))
1184  IF(.NOT.ALLOCATED(nodeevalue)) ALLOCATE(nodeevalue(nodespermeshcomponent(1)))
1185  IF(.NOT.ALLOCATED(nodekappavalue)) ALLOCATE(nodekappavalue(nodespermeshcomponent(1)))
1186 ! IF(.NOT.ALLOCATED(ElementNodesScales)) ALLOCATE(ElementNodesScales(NumberOfElements,NodesPerElement(1)))
1188 ! IF(.NOT.ALLOCATED(ElementNodes)) ALLOCATE(ElementNodes(NumberOfElements,NodesPerElement(1)))
1189  IF(.NOT.ALLOCATED(elementnodes)) ALLOCATE(elementnodes(numberofelements,maxnodesperelement))
1190 
1191  !chrm, 20.08.09
1192  IF(.NOT.ALLOCATED(nodeuvalue_analytic)) ALLOCATE(nodeuvalue_analytic(nodespermeshcomponent(1)))
1193  IF(.NOT.ALLOCATED(nodevvalue_analytic)) ALLOCATE(nodevvalue_analytic(nodespermeshcomponent(1)))
1194  IF(.NOT.ALLOCATED(nodewvalue_analytic)) ALLOCATE(nodewvalue_analytic(nodespermeshcomponent(1)))
1195  IF(.NOT.ALLOCATED(nodepvalue_analytic)) ALLOCATE(nodepvalue_analytic(nodespermeshcomponent(1)))
1196 
1197  IF(.NOT.ALLOCATED(nodeuvalue_error)) ALLOCATE(nodeuvalue_error(nodespermeshcomponent(1)))
1198  IF(.NOT.ALLOCATED(nodevvalue_error)) ALLOCATE(nodevvalue_error(nodespermeshcomponent(1)))
1199  IF(.NOT.ALLOCATED(nodewvalue_error)) ALLOCATE(nodewvalue_error(nodespermeshcomponent(1)))
1200  IF(.NOT.ALLOCATED(nodepvalue_error)) ALLOCATE(nodepvalue_error(nodespermeshcomponent(1)))
1201 
1202  enters("CMGUI OUTPUT",err,error,*999)
1203 
1204  field=>region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field
1205 
1206  NULLIFY(interpolation_parameters)
1207  NULLIFY(interpolated_point)
1208 
1209  CALL field_interpolation_parameters_initialise(field,interpolation_parameters &
1210  & ,err,error,*999)
1211  CALL field_interpolated_points_initialise(interpolation_parameters,interpolated_point,err,error,*999)
1212 
1213  DO i=1,numberofelements
1214 
1215  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
1216  & %ptr%topology%elements%elements(i)%basis%number_of_element_parameters
1217 
1218 ! DO J=1,NodesPerElement(1)
1219  DO j=1,nodesperelement(i)
1220 
1221  element_number=i
1222  xi_coordinates(1)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
1223  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%node_position_index(j,1)-1.0)/(region% &
1224  & equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
1225  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%number_of_nodes_xic(1)-1.0)
1226  xi_coordinates(2)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
1227  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%node_position_index(j,2)-1.0)/(region% &
1228  & equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
1229  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%number_of_nodes_xic(2)-1.0)
1230  xi_coordinates(3)=(region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
1231  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%node_position_index(j,3)-1.0)/(region% &
1232  & equations_sets%equations_sets(equations_set_global_number)%ptr%equations%interpolation% &
1233  & geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%number_of_nodes_xic(3)-1.0)
1234  IF(numberofdimensions==2)THEN
1235  stop 'Encas format only available for 3D hex and tets'
1236  END IF
1237 
1238  !K is global node number
1239  k=region%meshes%meshes(1)%ptr%topology(1)%ptr%elements%elements(i)%global_element_nodes(j)
1240 
1241  coordinates(1:3)=[1,1,1]
1242 
1243  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1244  & interpolation_parameters(field_var_type)%ptr,err,error,*999)
1245  CALL field_interpolate_xi(no_part_deriv,xi_coordinates,interpolated_point(field_var_type)%ptr,err,error,*999)
1246  nodexvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field%variables(1) &
1247  & %parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k)
1248  nodeyvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field%variables(1) &
1249  & %parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+nodespermeshcomponent(1))
1250 
1251  IF(numberofdimensions==3)THEN
1252  nodezvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field% &
1253  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+2*nodespermeshcomponent(1))
1254  END IF
1255 
1256 ! NodeUValue(K)=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%dependent%dependent_field%variables(1) &
1257  nodeuvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1258  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
1259 ! NodeVValue(K)=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%dependent%dependent_field%variables(1) &
1260  nodevvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1261  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters% &
1262  & cmiss%data_dp(k+nodespermeshcomponent(1))
1263 
1264  IF(numberofdimensions==3)THEN
1265  nodewvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1266 ! & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(K+2*NodesPerMeshComponent(1))
1267  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters% &
1268  & cmiss%data_dp(k+2*nodespermeshcomponent(1))
1269  END IF
1270 
1271 ! ! ! NodeUValue(K)=INTERPOLATED_POINT%VALUES(1,1)
1272 ! ! ! NodeVValue(K)=INTERPOLATED_POINT%VALUES(2,1)
1273 ! ! ! NodeWValue(K)=INTERPOLATED_POINT%VALUES(3,1)
1274 
1275  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
1276  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
1277  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
1278  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
1279  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype) )THEN
1280  IF(numberofdimensions==3)THEN
1281  nodepvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(4,1)
1282  ELSE IF(numberofdimensions==2)THEN
1283  nodepvalue(k)=interpolated_point(field_var_type)%ptr%VALUES(3,1)
1284  END IF
1285  END IF
1286 
1287  material_interpolation_type=region%equations_sets%equations_sets(equations_set_global_number)%ptr% &
1288  & materials%materials_field%variables(1)%COMPONENTS(1)%INTERPOLATION_TYPE
1289 
1290  IF(material_interpolation_type==field_node_based_interpolation)THEN
1291  nodemuvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1292  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k)
1293  noderhovalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1294  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+nodespermeshcomponent(1))
1295 
1296  IF(equations_set%SPECIFICATION(1)==equations_set_elasticity_class)THEN
1297  IF(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type)THEN
1298  IF( (equations_set%SPECIFICATION(3)==equations_set_compressible_finite_elasticity_subtype) &
1299  & .OR.(equations_set%SPECIFICATION(3)==equations_set_elasticity_darcy_inria_model_subtype) &
1300  & .OR.(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) &
1301  & .OR. (equations_set%SPECIFICATION(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) )THEN
1302  nodekappavalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1303  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+2*nodespermeshcomponent(1))
1304  END IF
1305  END IF
1306  END IF
1307  ELSE !default to FIELD_CONSTANT_INTERPOLATION
1308  nodemuvalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1309  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(1)
1310  noderhovalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1311  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(2)
1312  IF(equations_set%SPECIFICATION(1)==equations_set_elasticity_class)THEN
1313  IF(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type)THEN
1314  IF( (equations_set%SPECIFICATION(3)==equations_set_compressible_finite_elasticity_subtype) &
1315  & .OR.(equations_set%SPECIFICATION(3)==equations_set_elasticity_darcy_inria_model_subtype) &
1316  & .OR.(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) )THEN
1317  nodekappavalue=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1318  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(3)
1319  END IF
1320  END IF
1321  END IF
1322  END IF
1323 
1324  END DO
1325  END DO
1326 
1327  IF( numberofdimensions==3 )THEN
1328  !For 3D, the following call works ...
1329  lagrange_simplex=region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations% &
1330  & interpolation%geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%type
1331 ! lagrange_simplex=2
1332  ELSE IF ( numberofdimensions==1 )THEN
1334  ELSE IF ( numberofdimensions==2 )THEN
1335  !chrm, 20.08.09:
1336  ! ... but the above call does not work for 2D.
1337  !Thus, for 2D, we hard-wire it to 'quad':
1340  ELSE IF(maxnodesperelement==3.OR.maxnodesperelement==6.OR.maxnodesperelement==10) THEN
1342  ENDIF
1343  END IF
1344 
1348 
1349  DO i=1,numberofelements
1350 ! DO J=1,NodesPerElement(1)
1351  DO j=1,nodesperelement(i)
1352  elementnodes(i,j)=region%meshes%meshes(1)%ptr%topology(1)% &
1353  & ptr%elements%elements(i)%global_element_nodes(j)
1354  elementnodesscales(i,j)=1.0000000000000000e+00
1355  END DO
1356  END DO
1357 
1358  CALL fluid_mechanics_io_write_data_encas(name,equations_set)
1359 
1360 
1361  exits("FLUID_MECHANICS_IO_WRITE_ENCAS")
1362  RETURN
1363 999 errorsexits("FLUID_MECHANICS_IO_WRITE_ENCAS",err,error)
1364  RETURN
1365 
1366  END SUBROUTINE fluid_mechanics_io_write_encas
1367 
1368 
1369 
1370  ! OK
1371  !================================================================================================================================
1372  !
1373 
1375  SUBROUTINE fluid_mechanics_io_write_data_encas(NAME,EQUATIONS_SET)
1377  IMPLICIT NONE
1378  CHARACTER(14), INTENT(IN) :: NAME
1379  TYPE(varying_string) :: FILENAME
1380 ! CHARACTER :: FILENAME !<the prefix name of file.
1381  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1382  INTEGER:: I,K
1383  INTEGER(INTG) :: ERR
1384  TYPE(varying_string):: ERROR
1385  DOUBLE PRECISION:: velocity_magnitude
1386 ! CHARACTER(80):: OUTPUT_FILE
1387 
1388  !==================
1389  filename="./output/"//name//".geo"
1390  OPEN(unit=80, file=char(filename),status='unknown')
1391  WRITE(80,*)'OpenCMISS Exported Encas Model Geometry File'
1392  WRITE(80,*)''
1393  WRITE(80,*)'node id assign'
1394  WRITE(80,*)'element id assign'
1395  WRITE(80,*)'part'
1396  WRITE(80,*)' 1'
1397  WRITE(80,*)'OpenCMISS'
1398  WRITE(80,*)'coordinates'
1399  WRITE(80,'(i10)')nodespermeshcomponent(1)
1400 
1401  DO i = 1,nodespermeshcomponent(1)
1402  WRITE(80,'(e12.5)')nodexvalue(i)
1403  ENDDO
1404  DO i = 1,nodespermeshcomponent(1)
1405  WRITE(80,'(e12.5)')nodeyvalue(i)
1406  ENDDO
1407  DO i = 1,nodespermeshcomponent(1)
1408  WRITE(80,'(e12.5)')nodezvalue(i)
1409  ENDDO
1410  IF(lagrange_simplex==2) THEN !simplex
1411  IF(.NOT.ALLOCATED(simplexoutputhelp)) ALLOCATE(simplexoutputhelp(nodesperelement(1)))
1412  IF(maxnodesperelement==4) THEN
1413  WRITE(80,*)'tetra4'
1414  WRITE(80,'(i10)')numberofelements
1415  DO k = 1,numberofelements
1420  WRITE(80,'(4(i10))')simplexoutputhelp
1421  END DO
1422  ELSE IF(maxnodesperelement==10) THEN
1423  WRITE(80,*)'tetra10'
1424  WRITE(80,'(i10)')numberofelements
1425  DO k = 1,numberofelements
1436  WRITE(80,'(10(i10))')simplexoutputhelp
1437  END DO
1438  ELSE
1439  stop 'Encas format only available for 3D hex and tets'
1440  ENDIF
1441  ELSE IF (lagrange_simplex==1) THEN !hexa
1442  IF(.NOT.ALLOCATED(hexoutputhelp)) ALLOCATE(hexoutputhelp(nodesperelement(1)))
1443  IF(maxnodesperelement==8) THEN
1444  WRITE(80,*)'hexa8'
1445  WRITE(80,'(i10)')numberofelements
1446  DO k = 1,numberofelements
1447  hexoutputhelp(1)=elementnodes(k,1)
1448  hexoutputhelp(2)=elementnodes(k,2)
1449  hexoutputhelp(3)=elementnodes(k,4)
1450  hexoutputhelp(4)=elementnodes(k,3)
1451  hexoutputhelp(5)=elementnodes(k,5)
1452  hexoutputhelp(6)=elementnodes(k,6)
1453  hexoutputhelp(7)=elementnodes(k,8)
1454  hexoutputhelp(8)=elementnodes(k,7)
1455  WRITE(80,'(8(i10))')hexoutputhelp
1456  END DO
1457  ELSE IF(maxnodesperelement==27) THEN
1458  WRITE(80,*)'hexa8'
1459  WRITE(80,'(i10)')8*numberofelements
1460  DO k = 1,numberofelements
1461  hexoutputhelp(1)=elementnodes(k,1)
1462  hexoutputhelp(2)=elementnodes(k,2)
1463  hexoutputhelp(3)=elementnodes(k,5)
1464  hexoutputhelp(4)=elementnodes(k,4)
1465  hexoutputhelp(5)=elementnodes(k,10)
1466  hexoutputhelp(6)=elementnodes(k,11)
1467  hexoutputhelp(7)=elementnodes(k,14)
1468  hexoutputhelp(8)=elementnodes(k,13)
1469  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1470  hexoutputhelp(1)=elementnodes(k,2)
1471  hexoutputhelp(2)=elementnodes(k,3)
1472  hexoutputhelp(3)=elementnodes(k,6)
1473  hexoutputhelp(4)=elementnodes(k,5)
1474  hexoutputhelp(5)=elementnodes(k,11)
1475  hexoutputhelp(6)=elementnodes(k,12)
1476  hexoutputhelp(7)=elementnodes(k,15)
1477  hexoutputhelp(8)=elementnodes(k,14)
1478  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1479  hexoutputhelp(1)=elementnodes(k,4)
1480  hexoutputhelp(2)=elementnodes(k,5)
1481  hexoutputhelp(3)=elementnodes(k,8)
1482  hexoutputhelp(4)=elementnodes(k,7)
1483  hexoutputhelp(5)=elementnodes(k,13)
1484  hexoutputhelp(6)=elementnodes(k,14)
1485  hexoutputhelp(7)=elementnodes(k,17)
1486  hexoutputhelp(8)=elementnodes(k,16)
1487  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1488  hexoutputhelp(1)=elementnodes(k,5)
1489  hexoutputhelp(2)=elementnodes(k,6)
1490  hexoutputhelp(3)=elementnodes(k,9)
1491  hexoutputhelp(4)=elementnodes(k,8)
1492  hexoutputhelp(5)=elementnodes(k,14)
1493  hexoutputhelp(6)=elementnodes(k,15)
1494  hexoutputhelp(7)=elementnodes(k,18)
1495  hexoutputhelp(8)=elementnodes(k,17)
1496  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1497  hexoutputhelp(1)=elementnodes(k,10)
1498  hexoutputhelp(2)=elementnodes(k,11)
1499  hexoutputhelp(3)=elementnodes(k,14)
1500  hexoutputhelp(4)=elementnodes(k,13)
1501  hexoutputhelp(5)=elementnodes(k,19)
1502  hexoutputhelp(6)=elementnodes(k,20)
1503  hexoutputhelp(7)=elementnodes(k,23)
1504  hexoutputhelp(8)=elementnodes(k,22)
1505  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1506  hexoutputhelp(1)=elementnodes(k,11)
1507  hexoutputhelp(2)=elementnodes(k,12)
1508  hexoutputhelp(3)=elementnodes(k,15)
1509  hexoutputhelp(4)=elementnodes(k,14)
1510  hexoutputhelp(5)=elementnodes(k,20)
1511  hexoutputhelp(6)=elementnodes(k,21)
1512  hexoutputhelp(7)=elementnodes(k,24)
1513  hexoutputhelp(8)=elementnodes(k,23)
1514  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1515  hexoutputhelp(1)=elementnodes(k,13)
1516  hexoutputhelp(2)=elementnodes(k,14)
1517  hexoutputhelp(3)=elementnodes(k,17)
1518  hexoutputhelp(4)=elementnodes(k,16)
1519  hexoutputhelp(5)=elementnodes(k,22)
1520  hexoutputhelp(6)=elementnodes(k,23)
1521  hexoutputhelp(7)=elementnodes(k,26)
1522  hexoutputhelp(8)=elementnodes(k,25)
1523  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1524  hexoutputhelp(1)=elementnodes(k,14)
1525  hexoutputhelp(2)=elementnodes(k,15)
1526  hexoutputhelp(3)=elementnodes(k,18)
1527  hexoutputhelp(4)=elementnodes(k,17)
1528  hexoutputhelp(5)=elementnodes(k,23)
1529  hexoutputhelp(6)=elementnodes(k,24)
1530  hexoutputhelp(7)=elementnodes(k,27)
1531  hexoutputhelp(8)=elementnodes(k,26)
1532  WRITE(80,'(8(i10))')hexoutputhelp(1:8)
1533  END DO
1534  ELSE
1535  stop 'Encas format only available for 3D hex and tets'
1536  ENDIF
1537  ENDIF
1538  CLOSE(80)
1539  !==================
1540  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
1541  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
1542  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
1543  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
1544  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype) )THEN
1545  filename="./output/"//name//".scl1"
1546  OPEN(unit=81, file=char(filename),status='unknown')
1547  WRITE(81,*)'Absolute Pressure'
1548  WRITE(81,*)'part'
1549  WRITE(81,*)' 1'
1550  WRITE(81,*)'coordinates'
1551  DO i = 1,nodespermeshcomponent(1)
1552  WRITE(81,'(e12.5)')nodepvalue(i)
1553  ENDDO
1554  CLOSE(81)
1555  ENDIF
1556 
1557  !==================
1558  filename="./output/"//name//".scl2"
1559  OPEN(unit=82, file=char(filename),status='unknown')
1560  WRITE(82,*)'Velocity Magnitude'
1561  WRITE(82,*)'part'
1562  WRITE(82,*)' 1'
1563  WRITE(82,*)'coordinates'
1564  DO i = 1,nodespermeshcomponent(1)
1565  velocity_magnitude=sqrt(nodeuvalue(i)*nodeuvalue(i)+ &
1567  WRITE(82,'(e12.5)')velocity_magnitude
1568  ENDDO
1569  CLOSE(82)
1570 
1571  !==================
1572  filename="./output/"//name//".scl3"
1573  OPEN(unit=83, file=char(filename),status='unknown')
1574  WRITE(83,*)'X Velocity'
1575  WRITE(83,*)'part'
1576  WRITE(83,*)' 1'
1577  WRITE(83,*)'coordinates'
1578  DO i = 1,nodespermeshcomponent(1)
1579  WRITE(83,'(e12.5)')nodeuvalue(i)
1580  ENDDO
1581  CLOSE(83)
1582 
1583  !==================
1584  filename="./output/"//name//".scl4"
1585  OPEN(unit=84, file=char(filename),status='unknown')
1586  WRITE(84,*)'Y Velocity'
1587  WRITE(84,*)'part'
1588  WRITE(84,*)' 1'
1589  WRITE(84,*)'coordinates'
1590  DO i = 1,nodespermeshcomponent(1)
1591  WRITE(84,'(e12.5)')nodevvalue(i)
1592  ENDDO
1593  CLOSE(84)
1594 
1595  !==================
1596  filename="./output/"//name//".scl5"
1597  OPEN(unit=85, file=char(filename),status='unknown')
1598  WRITE(85,*)'Z Velocity'
1599  WRITE(85,*)'part'
1600  WRITE(85,*)' 1'
1601  WRITE(85,*)'coordinates'
1602  DO i = 1,nodespermeshcomponent(1)
1603  WRITE(85,'(e12.5)')nodewvalue(i)
1604  ENDDO
1605  CLOSE(85)
1606  !==================
1607  filename="./output/"//name//".vel"
1608  OPEN(unit=86, file=char(filename),status='unknown')
1609  WRITE(86,*)'Velocity'
1610  WRITE(86,*)'part'
1611  WRITE(86,*)' 1'
1612  WRITE(86,*)'coordinates'
1613  DO i = 1,nodespermeshcomponent(1)
1614  WRITE(86,'(e12.5)')nodeuvalue(i)
1615  ENDDO
1616  DO i = 1,nodespermeshcomponent(1)
1617  WRITE(86,'(e12.5)')nodevvalue(i)
1618  ENDDO
1619  DO i = 1,nodespermeshcomponent(1)
1620  WRITE(86,'(e12.5)')nodewvalue(i)
1621  ENDDO
1622  CLOSE(86)
1623  !==================
1624  CALL write_string(general_output_type,"Writing Encas data...",err,error,*999)
1625  RETURN
1626 999 errorsexits("FLUID_MECHANICS_IO_WRITE_DATA_ENCAS",err,error)
1627  RETURN
1628 
1630 
1631  ! OK
1632  !================================================================================================================================
1633  !
1634 
1635 
1637  SUBROUTINE fluid_mechanics_io_write_master_encas_ppe(NAME,start_time_step,number_of_timesteps,time_increment)
1639  IMPLICIT NONE
1640 
1641  INTEGER:: I,J,start_time_step,number_of_timesteps
1642  DOUBLE PRECISION:: time_increment,time
1643 
1644  CHARACTER(9), INTENT(IN) :: NAME
1645  TYPE(varying_string) :: FILENAME
1646 ! CHARACTER :: FILENAME !<the prefix name of file.
1647  INTEGER(INTG) :: ERR
1648  TYPE(varying_string):: ERROR
1649 
1650 
1651  filename="./output/"//name//".case"
1652  OPEN(unit=87, file=char(filename),status='unknown')
1653 
1654  WRITE(87,*)'FORMAT'
1655  WRITE(87,*)'type: ensight gold'
1656  WRITE(87,*)''
1657  WRITE(87,*)'GEOMETRY'
1658  WRITE(87,*)'model: ./OpenCMISS.geo'
1659  WRITE(87,*)''
1660  WRITE(87,*)'VARIABLE'
1661  WRITE(87,*)''
1662  WRITE(87,*)'scalar per node: 1 Magnitude OpenCMISS.Magnitude_**'
1663  WRITE(87,*)'scalar per node: 1 Speed_SumSquares OpenCMISS.Speed_SumSquares_**'
1664  WRITE(87,*)'scalar per node: 1 Pressure OpenCMISS.Pressure_**'
1665  WRITE(87,*)'vector per node: 1 Velocity OpenCMISS.Velocity_**'
1666  WRITE(87,*)''
1667  WRITE(87,*)'TIME'
1668  WRITE(87,*)'time set: 1 Model'
1669  WRITE(87,*)'number of steps: ', number_of_timesteps
1670  WRITE(87,*)'filename start number: ', start_time_step
1671  WRITE(87,*)'filename increment: 1'
1672  WRITE(87,'(" time values:")',advance="NO")
1673  j=1
1674  DO i=1,number_of_timesteps
1675  j=j+1
1676  time=i*time_increment
1677  WRITE(87,'(e13.5)',advance="NO")time
1678  IF(j==8) THEN
1679  WRITE(87,*) ' '
1680  j=0
1681  ENDIF
1682  ENDDO
1683 
1684  CLOSE(87)
1685 
1686 
1687  CALL write_string(general_output_type,"Writing Encas PPE master...",err,error,*999)
1688  RETURN
1689 999 errorsexits("FLUID_MECHANICS_IO_WRITE_MASTER_ENCAS_PPE",err,error)
1690  RETURN
1691 
1693 
1694 
1695  ! OK
1696  !================================================================================================================================
1697  !
1698 
1699 
1701  SUBROUTINE fluid_mechanics_io_write_master_encas(NAME,start_time_step,number_of_timesteps,time_increment)
1703  IMPLICIT NONE
1704 
1705  INTEGER:: I,J,start_time_step,number_of_timesteps
1706  DOUBLE PRECISION:: time_increment,time
1707 
1708  CHARACTER(14), INTENT(IN) :: NAME
1709  TYPE(varying_string) :: FILENAME
1710 ! CHARACTER :: FILENAME !<the prefix name of file.
1711  INTEGER(INTG) :: ERR
1712  TYPE(varying_string):: ERROR
1713 
1714 
1715  filename="./output/"//name//".encas"
1716  OPEN(unit=87, file=char(filename),status='unknown')
1717 
1718  WRITE(87,*)'FORMAT'
1719  WRITE(87,*)'type: ensight gold'
1720  WRITE(87,*)''
1721  WRITE(87,*)'GEOMETRY'
1722  WRITE(87,*)'model: 1 ./TIME_STEP_****.geo'
1723  WRITE(87,*)''
1724  WRITE(87,*)'VARIABLE'
1725  WRITE(87,*)'scalar per node: 1 pressure ./TIME_STEP_****.scl1'
1726  WRITE(87,*)'scalar per node: 1 velocity-magnitude ./TIME_STEP_****.scl2'
1727  WRITE(87,*)'scalar per node: 1 velocity-u ./TIME_STEP_****.scl3'
1728  WRITE(87,*)'scalar per node: 1 velocity-v ./TIME_STEP_****.scl4'
1729  WRITE(87,*)'scalar per node: 1 velocity-w ./TIME_STEP_****.scl5'
1730  WRITE(87,*)'vector per node: 1 velocity ./TIME_STEP_****.vel'
1731  WRITE(87,*)''
1732  WRITE(87,*)'TIME'
1733  WRITE(87,*)'time set: 1 Model'
1734  WRITE(87,*)'number of steps: ', number_of_timesteps
1735  WRITE(87,*)'filename start number: ', start_time_step
1736  WRITE(87,*)'filename increment: 1'
1737  WRITE(87,'(" time values:")',advance="NO")
1738  j=1
1739  DO i=1,number_of_timesteps
1740  j=j+1
1741  time=i*time_increment
1742  WRITE(87,'(e13.5)',advance="NO")time
1743  IF(j==8) THEN
1744  WRITE(87,*) ' '
1745  j=0
1746  ENDIF
1747  ENDDO
1748 
1749  CLOSE(87)
1750 
1751 
1752  CALL write_string(general_output_type,"Writing Encas master...",err,error,*999)
1753  RETURN
1754 999 errorsexits("FLUID_MECHANICS_IO_WRITE_MASTER_ENCAS",err,error)
1755  RETURN
1756 
1758 
1759  ! OK
1760  !================================================================================================================================
1761  !
1762 
1764  SUBROUTINE fluid_mechanics_io_write_encas_block(REGION, EQUATIONS_SET_GLOBAL_NUMBER, NAME, ERR, ERROR,*)
1766  !Argument variables
1767  TYPE(region_type), POINTER :: REGION
1768 ! TYPE(VARYING_STRING), INTENT(IN) :: NAME !<the prefix name of file.
1769  CHARACTER(14) :: NAME
1770  INTEGER(INTG) :: EQUATIONS_SET_GLOBAL_NUMBER
1771  INTEGER(INTG) :: ERR
1772  TYPE(varying_string):: ERROR
1773  !Local Variables
1774  INTEGER(INTG):: I,J,K !,icompartment
1775 ! INTEGER(INTG):: MATERIAL_INTERPOLATION_TYPE
1776  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1777 ! TYPE(FIELD_TYPE), POINTER :: EQUATIONS_SET_FIELD_FIELD !<A pointer to the equations set field
1778 ! INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
1779 
1780  enters("FLUID_MECHANICS_IO_WRITE_ENCAS_BLOCK",err,error,*999)
1781 
1782  IF (ALLOCATED(nodesperelement)) DEALLOCATE(nodesperelement)
1783  IF (ALLOCATED(nodespermeshcomponent)) DEALLOCATE(nodespermeshcomponent)
1784  IF (ALLOCATED(xi_coordinates)) DEALLOCATE(xi_coordinates)
1785  IF (ALLOCATED(coordinates)) DEALLOCATE(coordinates)
1786  IF (ALLOCATED(nodexvalue)) DEALLOCATE(nodexvalue)
1787  IF (ALLOCATED(nodeyvalue)) DEALLOCATE(nodeyvalue)
1788  IF (ALLOCATED(nodezvalue)) DEALLOCATE(nodezvalue)
1789  IF (ALLOCATED(nodeuvalue)) DEALLOCATE(nodeuvalue)
1790  IF (ALLOCATED(nodevvalue)) DEALLOCATE(nodevvalue)
1791  IF (ALLOCATED(nodewvalue)) DEALLOCATE(nodewvalue)
1792  IF (ALLOCATED(nodeuvalueorg)) DEALLOCATE(nodeuvalueorg)
1793  IF (ALLOCATED(nodewvalueorg)) DEALLOCATE(nodevvalueorg)
1794  IF (ALLOCATED(nodevvalueorg)) DEALLOCATE(nodewvalueorg)
1795  IF (ALLOCATED(nodepvalue)) DEALLOCATE(nodepvalue)
1796  IF (ALLOCATED(nodemuvalue)) DEALLOCATE(nodemuvalue)
1797  IF (ALLOCATED(nodelabelvalue)) DEALLOCATE(nodelabelvalue)
1798  IF (ALLOCATED(noderhovalue)) DEALLOCATE(noderhovalue)
1799  IF (ALLOCATED(nodekappavalue)) DEALLOCATE(nodekappavalue)
1800  IF (ALLOCATED(elementnodesscales)) DEALLOCATE(elementnodesscales)
1801  IF (ALLOCATED(elementnodes)) DEALLOCATE(elementnodes)
1802 
1803  knot = '0'
1804  nms(1) = '1'
1805  nms(2) = '2'
1806  nms(3) = '3'
1807  nms(4) = '4'
1808  nms(5) = '5'
1809  nms(6) = '6'
1810  nms(7) = '7'
1811  nms(8) = '8'
1812  nms(9) = '9'
1813 
1814  k = 9
1815  DO i = 1,9
1816  k = k + 1
1817  nms(k) = trim(nms(i))//trim(knot)
1818  DO j = 1,9
1819  k = k + 1
1820  nms(k) = trim(nms(i))//trim(nms(j))
1821  END DO
1822  END DO
1823 
1824  equations_set => region%equations_sets%equations_sets(equations_set_global_number)%ptr
1825 
1826 !---tob
1827 ! FIELD_VAR_TYPE=EQUATIONS_SET%EQUATIONS%EQUATIONS_MAPPING%LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE%VARIABLE_TYPE
1828 ! ! '1' associated with linear matrix
1829 
1830  var_idx = 1
1831  field_var_type = field_u_variable_type
1832  SELECT CASE(equations_set%SPECIFICATION(1))
1834  SELECT CASE(equations_set%SPECIFICATION(2))
1836  SELECT CASE(equations_set%SPECIFICATION(3))
1839  var_idx = 3
1840  field_var_type = field_v_variable_type
1841  END SELECT
1842  END SELECT
1843  END SELECT
1844 
1845  parameter_set_idx = 1
1846  SELECT CASE(equations_set%SPECIFICATION(1))
1848  SELECT CASE(equations_set%SPECIFICATION(2))
1850  SELECT CASE(equations_set%SPECIFICATION(3))
1853 ! parameter_set_idx = 3 ! of shared dependent variable field
1854  END SELECT
1855  END SELECT
1856  END SELECT
1857 !---toe
1858 
1859 ! NumberOfFields=REGION%fields%number_of_fields
1860 ! Hack for ALE... to be removed later
1861  numberoffields=3
1862  numberofdimensions=region%coordinate_system%number_of_dimensions
1863 ! NumberOfVariableComponents=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%dependent%dependent_field% &
1864 ! & variables(1)%number_of_components
1865  numberofvariablecomponents=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1866  & variables(var_idx)%number_of_components
1867 
1868  numberofmaterialcomponents=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
1869  & variables(1)%number_of_components
1870  numberofelements=region%meshes%meshes(1)%ptr%number_of_elements
1871  numberofmeshcomponents=region%meshes%meshes(1)%ptr%number_of_components
1872 ! IF(.NOT.ALLOCATED(NodesPerElement)) ALLOCATE(NodesPerElement(NumberOfMeshComponents))
1873 
1874  IF(.NOT.ALLOCATED(nodesperelement)) ALLOCATE(nodesperelement(max(numberofmeshcomponents,numberofelements)))
1875 
1878 
1879  DO i=1,numberofmeshcomponents
1880  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
1881  & %ptr%topology%elements%elements(1)%basis%number_of_element_parameters
1882  nodespermeshcomponent(i)=region%meshes%meshes(1)%ptr%topology(i)%ptr%nodes%numberOfNodes
1883  END DO
1884 
1885 
1886 ! MaxNodesPerElement=NodesPerElement(1)
1888 
1889 
1890  DO i=1,numberofelements
1891  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
1892  & %ptr%topology%elements%elements(i)%basis%number_of_element_parameters
1894  END DO
1895 
1896 
1897  IF(.NOT.ALLOCATED(xi_coordinates)) ALLOCATE(xi_coordinates(numberofdimensions))
1898  IF(.NOT.ALLOCATED(coordinates)) ALLOCATE(coordinates(numberofdimensions))
1899  IF(.NOT.ALLOCATED(nodexvalue)) ALLOCATE(nodexvalue(nodespermeshcomponent(1)))
1900  IF(.NOT.ALLOCATED(nodeyvalue)) ALLOCATE(nodeyvalue(nodespermeshcomponent(1)))
1901  IF(.NOT.ALLOCATED(nodezvalue)) ALLOCATE(nodezvalue(nodespermeshcomponent(1)))
1902  IF(.NOT.ALLOCATED(nodeuvalue)) ALLOCATE(nodeuvalue(nodespermeshcomponent(1)))
1903  IF(.NOT.ALLOCATED(nodevvalue)) ALLOCATE(nodevvalue(nodespermeshcomponent(1)))
1904  IF(.NOT.ALLOCATED(nodewvalue)) ALLOCATE(nodewvalue(nodespermeshcomponent(1)))
1905  IF(.NOT.ALLOCATED(nodeuvalueorg)) ALLOCATE(nodeuvalueorg(nodespermeshcomponent(1)))
1906  IF(.NOT.ALLOCATED(nodevvalueorg)) ALLOCATE(nodevvalueorg(nodespermeshcomponent(1)))
1907  IF(.NOT.ALLOCATED(nodewvalueorg)) ALLOCATE(nodewvalueorg(nodespermeshcomponent(1)))
1908  IF(.NOT.ALLOCATED(nodepvalue)) ALLOCATE(nodepvalue(nodespermeshcomponent(1)))
1909  IF(.NOT.ALLOCATED(nodemuvalue)) ALLOCATE(nodemuvalue(nodespermeshcomponent(1)))
1910  IF(.NOT.ALLOCATED(nodelabelvalue)) ALLOCATE(nodelabelvalue(nodespermeshcomponent(1)))
1911  IF(.NOT.ALLOCATED(noderhovalue)) ALLOCATE(noderhovalue(nodespermeshcomponent(1)))
1912  IF(.NOT.ALLOCATED(nodekappavalue)) ALLOCATE(nodekappavalue(nodespermeshcomponent(1)))
1913  IF(.NOT.ALLOCATED(elementnodesscales)) ALLOCATE(elementnodesscales(numberofelements,nodesperelement(1)))
1914  IF(.NOT.ALLOCATED(elementnodes)) ALLOCATE(elementnodes(numberofelements,nodesperelement(1)))
1915 
1916  enters("CMGUI OUTPUT",err,error,*999)
1917 
1918  field=>region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field
1919 
1920  NULLIFY(interpolation_parameters)
1921  NULLIFY(interpolated_point)
1922 
1923  CALL field_interpolation_parameters_initialise(field,interpolation_parameters,err,error,*999)
1924  CALL field_interpolated_points_initialise(interpolation_parameters,interpolated_point,err,error,*999)
1925 
1926  DO i=1,numberofelements
1927 
1928  nodesperelement(i)=region%fields%fields(1)%ptr%geometric_field%decomposition%domain(1) &
1929  & %ptr%topology%elements%elements(i)%basis%number_of_element_parameters
1930 
1931 ! DO J=1,NodesPerElement(1)
1932  DO j=1,nodesperelement(i)
1933 
1934  element_number=i
1935 ! ! ! XI_COORDINATES(1)=(REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
1936 ! ! ! & geometric_interp_parameters(FIELD_U_VARIABLE_TYPE)%ptr%bases(1)%ptr%node_position_index(J,1)-1.0)/(REGION% &
1937 ! ! ! & equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
1938 ! ! ! & geometric_interp_parameters(FIELD_U_VARIABLE_TYPE)%ptr%bases(1)%ptr%number_of_nodes_xi(1)-1.0)
1939 ! ! ! XI_COORDINATES(2)=(REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
1940 ! ! ! & geometric_interp_parameters(FIELD_U_VARIABLE_TYPE)%ptr%bases(1)%ptr%node_position_index(J,2)-1.0)/(REGION% &
1941 ! ! ! & equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
1942 ! ! ! & geometric_interp_parameters(FIELD_U_VARIABLE_TYPE)%ptr%bases(1)%ptr%number_of_nodes_xi(2)-1.0)
1943 ! ! ! XI_COORDINATES(3)=(REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
1944 ! ! ! & geometric_interp_parameters(FIELD_U_VARIABLE_TYPE)%ptr%bases(1)%ptr%node_position_index(J,3)-1.0)/(REGION% &
1945 ! ! ! & equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%equations%interpolation% &
1946 ! ! ! & geometric_interp_parameters(FIELD_U_VARIABLE_TYPE)%ptr%bases(1)%ptr%number_of_nodes_xi(3)-1.0)
1947 ! ! ! IF(NumberOfDimensions==2)THEN
1948 ! ! ! STOP 'Encas format only available for 3D hex and tets'
1949 ! ! ! END IF
1950 ! ! !
1951  !K is global node number
1952  k=region%meshes%meshes(1)%ptr%topology(1)%ptr%elements%elements(i)%global_element_nodes(j)
1953 
1954  coordinates(1:3)=[1,1,1]
1955 ! ! !
1956 ! ! ! CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER, &
1957 ! ! ! & INTERPOLATION_PARAMETERS(FIELD_VAR_TYPE)%ptr,ERR,ERROR,*999)
1958 ! ! ! CALL FIELD_INTERPOLATE_XI(NO_PART_DERIV,XI_COORDINATES,INTERPOLATED_POINT(FIELD_VAR_TYPE)%ptr,ERR,ERROR,*999)
1959  nodexvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field%variables(1) &
1960  & %parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k)
1961  nodeyvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field%variables(1) &
1962  & %parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+nodespermeshcomponent(1))
1963 
1964  IF(numberofdimensions==3)THEN
1965  nodezvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%geometry%geometric_field% &
1966  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(k+2*nodespermeshcomponent(1))
1967  END IF
1968 
1969  IF((equations_set%SPECIFICATION(1)==equations_set_classical_field_class) &
1970  & .AND.(equations_set%SPECIFICATION(2)==equations_set_poisson_equation_type) &
1971  & .AND.(equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype))THEN
1972  nodeuvalue(k)=0.0_dp
1973  nodevvalue(k)=0.0_dp
1974  ELSE
1975 ! NodeUValue(K)=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%source%source_field%variables(1) &
1976  nodeuvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%source%source_field% &
1977  & variables(var_idx)%parameter_sets%parameter_sets(2)%ptr%parameters%cmiss%data_dp(k)
1978 ! ! NodeVValue(K)=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%source%source_field%variables(1) &
1979  nodevvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%source%source_field% &
1980  & variables(var_idx)%parameter_sets%parameter_sets(2)%ptr%parameters%cmiss%data_dp(k+nodespermeshcomponent(1))
1981  ENDIF
1982 
1983  parameter_set_idx = 2
1984  nodeuvalueorg(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1985  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
1986  parameter_set_idx = 3
1987  nodevvalueorg(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1988  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
1989 
1990  IF(numberofdimensions==3)THEN
1991  parameter_set_idx = 4
1992  nodewvalueorg(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
1993  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
1994 
1995  IF((equations_set%SPECIFICATION(1)==equations_set_classical_field_class) &
1996  & .AND.(equations_set%SPECIFICATION(2)==equations_set_poisson_equation_type) &
1997  & .AND.(equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype))THEN
1998  nodewvalue(k)=0.0_dp
1999  else!
2000 ! ! NodeWValue(K)=REGION%equations_sets%equations_sets(EQUATIONS_SET_GLOBAL_NUMBER)%ptr%source%source_field%variables(1) &
2001  nodewvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%source%source_field% &
2002  & variables(var_idx)%parameter_sets%parameter_sets(2)%ptr%parameters%cmiss%data_dp(k+2*nodespermeshcomponent(1))
2003  END IF
2004  END IF
2005 
2006 
2007 ! ! ! NodeUValue(K)=INTERPOLATED_POINT%VALUES(1,1)
2008 ! ! ! NodeVValue(K)=INTERPOLATED_POINT%VALUES(2,1)
2009 ! ! ! NodeWValue(K)=INTERPOLATED_POINT%VALUES(3,1)
2010 
2011  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
2012  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
2013  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
2014  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
2015  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype) &
2016  & .OR. (equations_set%SPECIFICATION(1)==equations_set_classical_field_class) &
2017  & .AND.(equations_set%SPECIFICATION(2)==equations_set_poisson_equation_type) &
2018  & .AND.((equations_set%SPECIFICATION(3)==equations_set_linear_pressure_poisson_subtype) &
2019  & .OR.(equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype) &
2020  & .OR.(equations_set%SPECIFICATION(3)==equations_set_ale_pressure_poisson_subtype) &
2021  & .OR.(equations_set%SPECIFICATION(3)==equations_set_nonlinear_pressure_poisson_subtype))) THEN
2022  parameter_set_idx = 1
2023  nodepvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2024 ! & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(K)
2025  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
2026  END IF
2027 
2028  nodemuvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
2029  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(1)
2030  noderhovalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%materials%materials_field% &
2031  & variables(1)%parameter_sets%parameter_sets(1)%ptr%parameters%cmiss%data_dp(2)
2032 
2033  parameter_set_idx = 5
2034  nodelabelvalue(k)=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2035  & variables(var_idx)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(k)
2036 
2037  END DO
2038  END DO
2039 
2040 ! NodeMUValue=0.0_DP
2041 ! NodeRHOValue=0.0_DP
2042 
2043  IF( numberofdimensions==3 )THEN
2044  !For 3D, the following call works ...
2045  lagrange_simplex=region%equations_sets%equations_sets(equations_set_global_number)%ptr%equations% &
2046  & interpolation%geometric_interp_parameters(field_u_variable_type)%ptr%bases(1)%ptr%type
2047 ! lagrange_simplex=2
2048  ELSE
2049  !chrm, 20.08.09:
2050  ! ... but the above call does not work for 2D.
2051  !Thus, for 2D, we hard-wire it to 'quad':
2054  ELSE IF(maxnodesperelement==3.OR.maxnodesperelement==6.OR.maxnodesperelement==10) THEN
2056  ENDIF
2057  END IF
2058 
2059  !This is for Poisson-Flow problems only
2063  IF(equations_set%SPECIFICATION(1)==equations_set_classical_field_class &
2064  & .AND.equations_set%SPECIFICATION(2)==equations_set_poisson_equation_type &
2065  & .AND.equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype) THEN
2067  ENDIF
2068 
2069 
2073 
2074 
2075  DO i=1,numberofelements
2076 ! DO J=1,NodesPerElement(1)
2077  DO j=1,nodesperelement(i)
2078  elementnodes(i,j)=region%meshes%meshes(1)%ptr%topology(1)% &
2079  & ptr%elements%elements(i)%global_element_nodes(j)
2080  elementnodesscales(i,j)=1.0000000000000000e+00
2081  END DO
2082  END DO
2083 
2084  CALL fluid_mechanics_io_write_nodes_cmgui(name,equations_set)
2087 
2088 
2089  exits("FLUID_MECHANICS_IO_WRITE_ENCAS_BLOCK")
2090  RETURN
2091 999 errorsexits("FLUID_MECHANICS_IO_WRITE_ENCAS_BLOCK",err,error)
2092  RETURN
2093 
2095 
2096 
2097  ! OK
2098  !================================================================================================================================
2099  !
2100 
2101  SUBROUTINE fluid_mechanics_io_write_fitted_field(REGION, EQUATIONS_SET_GLOBAL_NUMBER, NAME, ERR, ERROR,*)
2103  !Argument variables
2104  TYPE(region_type), POINTER :: REGION
2105 ! TYPE(VARYING_STRING), INTENT(IN) :: NAME !<the prefix name of file.
2106  CHARACTER(7) :: NAME
2107  INTEGER(INTG) :: EQUATIONS_SET_GLOBAL_NUMBER
2108  INTEGER(INTG) :: ERR
2109  TYPE(varying_string):: ERROR
2110  !Local Variables
2111  INTEGER(INTG):: NODES_PER_COMPONENT,I!,J,K !,icompartment
2112 ! INTEGER(INTG):: MATERIAL_INTERPOLATION_TYPE
2113 ! TYPE(EQUATIONS_SET_TYPE), POINTER :: EQUATIONS_SET !<A pointer to the equations set
2114 ! TYPE(FIELD_TYPE), POINTER :: EQUATIONS_SET_FIELD_FIELD !<A pointer to the equations set field
2115 ! INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
2116  TYPE(varying_string) :: FILENAME
2117  DOUBLE PRECISION:: fit
2118 
2119  enters("FLUID_MECHANICS_IO_WRITE_FITTED_FIELD",err,error,*999)
2120 
2121  nodes_per_component=region%equations_sets%equations_sets(1)%ptr%dependent% &
2122  & dependent_field%variables(1)%components(1)%param_to_dof_map%node_param2dof_map%number_of_node_parameters
2123 
2124  filename="./output/VEL_"//name//".fit"
2125  OPEN(unit=81, file=char(filename),status='unknown')
2126  WRITE(81,*)nodes_per_component*3
2127  parameter_set_idx = 1
2128  DO i = 1,nodes_per_component*3
2129  fit=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2130  & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(i)
2131  WRITE(81,'(e12.5)') fit
2132  ENDDO
2133  CLOSE(81)
2134 
2135  filename="./output/LAM_"//name//".err"
2136  OPEN(unit=81, file=char(filename),status='unknown')
2137  WRITE(81,*)nodes_per_component
2138  parameter_set_idx = 1
2139  DO i = nodes_per_component*3+1,nodes_per_component*4
2140  fit=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2141  & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(i)
2142  WRITE(81,'(e12.5)') fit
2143  ENDDO
2144  CLOSE(81)
2145 
2146  filename="./output/U_"//name//".fit"
2147  OPEN(unit=81, file=char(filename),status='unknown')
2148 ! WRITE(81,*)NodesPerMeshComponent(1)
2149  parameter_set_idx = 1
2150  DO i = 1,nodes_per_component
2151  fit=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2152  & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(i)
2153  WRITE(81,'(e12.5)') fit
2154  ENDDO
2155  CLOSE(81)
2156 
2157  filename="./output/V_"//name//".fit"
2158  OPEN(unit=81, file=char(filename),status='unknown')
2159 ! WRITE(81,*)NodesPerMeshComponent(1)
2160  parameter_set_idx = 1
2161  DO i = nodes_per_component+1,nodes_per_component*2
2162  fit=region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2163  & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(i)
2164  WRITE(81,'(e12.5)') fit
2165  ENDDO
2166  CLOSE(81)
2167 
2168  filename="./output/W_"//name//".fit"
2169  OPEN(unit=81, file=char(filename),status='unknown')
2170 ! WRITE(81,*)NodesPerMeshComponent(1)
2171  parameter_set_idx = 1
2172  DO i = nodes_per_component*2+1,nodes_per_component*3
2173  WRITE(81,'(e12.5)') region%equations_sets%equations_sets(equations_set_global_number)%ptr%dependent%dependent_field% &
2174  & variables(1)%parameter_sets%parameter_sets(parameter_set_idx)%ptr%parameters%cmiss%data_dp(i)
2175  ENDDO
2176  CLOSE(81)
2177 
2178  exits("FLUID_MECHANICS_IO_WRITE_FITTED_FIELD")
2179  RETURN
2180 999 errorsexits("FLUID_MECHANICS_IO_WRITE_FITTED_FIELD",err,error)
2181  RETURN
2182 
2184 
2185 
2186  ! OK
2187  !================================================================================================================================
2188  !
2189 
2193  IMPLICIT NONE
2194  CHARACTER(14), INTENT(IN) :: NAME
2195  TYPE(varying_string) :: FILENAME
2196 ! CHARACTER :: FILENAME !<the prefix name of file.
2197  INTEGER:: I
2198  INTEGER(INTG) :: ERR
2199  TYPE(varying_string):: ERROR
2200 ! DOUBLE PRECISION:: velocity_magnitude
2201 ! CHARACTER(80):: OUTPUT_FILE
2202 
2203  !==================
2204 
2205  filename="./output/"//name//".scl1"
2206  OPEN(unit=81, file=char(filename),status='unknown')
2207  WRITE(81,*)'Absolute Pressure'
2208  WRITE(81,*)'part'
2209  WRITE(81,*)' 1'
2210  WRITE(81,*)'block'
2211  DO i = 1,nodespermeshcomponent(1)
2212  WRITE(81,'(e12.5)')nodepvalue(i)
2213  ENDDO
2214  CLOSE(81)
2215 
2216  CALL write_string(general_output_type,"Writing Encas data...",err,error,*999)
2217  RETURN
2218 999 errorsexits("FLUID_MECHANICS_IO_WRITE_DATA_ENCAS_BLOCK",err,error)
2219  RETURN
2220 
2222 
2223 
2224 
2225  ! OK
2226  !================================================================================================================================
2227  !
2228 
2230  SUBROUTINE fluid_mechanics_io_write_nodes_cmgui(NAME,EQUATIONS_SET)
2232  IMPLICIT NONE
2233 
2234  CHARACTER(len=14), INTENT(IN) :: NAME
2235  TYPE(varying_string) :: FILENAME
2236  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2237  INTEGER(INTG):: I
2238  INTEGER(INTG) :: ERR, ierror
2239  TYPE(varying_string):: ERROR
2240  LOGICAL:: ANALYTIC
2241 
2242  analytic=.false.
2243  IF( darcy%ANALYTIC ) analytic=.true.
2244 
2245  IF( darcy%ANALYTIC ) THEN
2248  END IF
2249 
2250  IF(ASSOCIATED(equations_set%ANALYTIC)) analytic=.true.
2251 
2252  filename="./output/"//name//".exnode"
2253  OPEN(unit=14, file=char(filename), status='unknown', iostat=ierror)
2254 
2255 ! WRITING HEADER INFORMATION
2256 
2257  WRITE(14,*) 'Group name: OpenCMISS'
2258 
2259  IF( analytic ) THEN
2260  WRITE(14,*) '#Fields=',trim(nms(numberoffields + 2))
2261  ELSE
2262  WRITE(14,*) '#Fields=',trim(nms(numberoffields))
2263  END IF
2264 
2265  valueindex=1
2266  WRITE(14,*) ' 1) coordinates, coordinate, rectangular cartesian, #Components=',trim(nms(numberofdimensions))
2267 
2268  DO i=1,numberofdimensions
2269  IF(i==1) THEN
2270  WRITE(14,*) ' x. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2271  ELSE IF(i==2) THEN
2272  WRITE(14,*) ' y. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2273  ELSE
2274  WRITE(14,*) ' z. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2275  END IF
2277  END DO
2278 
2279  WRITE(14,*) ' 2) general, field, rectangular cartesian, #Components=',trim(nms(numberofvariablecomponents))
2280 
2282  WRITE(14,*) ' ',trim(nms(i)),'. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2284  END DO
2285 
2286  IF( output_material) THEN
2287  WRITE(14,*) ' 3) material, field, rectangular cartesian, #Components=',trim(nms(numberofmaterialcomponents))
2288 
2290  WRITE(14,*) ' ',trim(nms(i)),'. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2292  END DO
2293  ENDIF
2294 
2295  IF( output_source ) THEN !Watch out that no numbering conflict occurs with Analytic: 4.)
2296  WRITE(14,*) ' 4) source, field, rectangular cartesian, #Components=',trim(nms(numberofsourcecomponents))
2297 
2299  WRITE(14,*) ' ',trim(nms(i)),'. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2301  END DO
2302  END IF
2303 
2304 
2305  IF( analytic ) THEN
2306  WRITE(14,*) ' 4) exact, field, rectangular cartesian, #Components=',trim(nms(numberofvariablecomponents))
2308  WRITE(14,*) ' ',trim(nms(i)),'. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2310  END DO
2311 
2312  WRITE(14,*) ' 5) error, field, rectangular cartesian, #Components=',trim(nms(numberofvariablecomponents))
2314  WRITE(14,*) ' ',trim(nms(i)),'. Value index= ',trim(nms(valueindex)),', #Derivatives= 0'
2316  END DO
2317  END IF
2318 
2319 
2320 ! NOW WRITE NODE INFORMATION
2321 
2322  DO i = 1,nodespermeshcomponent(1)
2323  !DO I = 1,DofsPerMeshComponent(1)
2324 
2325  WRITE(14,*) ' Node: ',i
2326  WRITE(14,'(" ", es25.16 )')nodexvalue(i)
2327 
2328  IF(numberofdimensions==2 .OR. numberofdimensions==3) THEN
2329  WRITE(14,'(" ", es25.16 )')nodeyvalue(i)
2330  END IF
2331 
2332  IF(numberofdimensions==3) THEN
2333  WRITE(14,'(" ", es25.16 )')nodezvalue(i)
2334  END IF
2335 
2336  WRITE(14,'(" ", es25.16 )')nodeuvalue(i)
2337 
2338  IF(numberofdimensions==2 .OR. numberofdimensions==3) THEN
2339  WRITE(14,'(" ", es25.16 )')nodevvalue(i)
2340  ENDIF
2341 
2342  IF(numberofdimensions==3) THEN
2343  WRITE(14,'(" ", es25.16 )')nodewvalue(i)
2344  END IF
2345 
2346  IF(numberofdimensions/=1) THEN
2347  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
2348  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
2349  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
2350  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
2351  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype))THEN
2352  WRITE(14,'(" ", es25.16 )')nodepvalue(i)
2353  END IF
2354  IF((equations_set%SPECIFICATION(1)==equations_set_classical_field_class) &
2355  & .AND.(equations_set%SPECIFICATION(2)==equations_set_poisson_equation_type) &
2356  & .AND.(equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype)) THEN
2357  WRITE(14,'(" ", es25.16 )')nodeuvalueorg(i)
2358  WRITE(14,'(" ", es25.16 )')nodevvalueorg(i)
2359  WRITE(14,'(" ", es25.16 )')nodewvalueorg(i)
2360  END IF
2361  IF((equations_set%SPECIFICATION(1)==equations_set_classical_field_class) &
2362  & .AND.(equations_set%SPECIFICATION(2)==equations_set_poisson_equation_type) &
2363  & .AND.((equations_set%SPECIFICATION(3)==equations_set_linear_pressure_poisson_subtype) &
2364  & .OR.(equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype) &
2365  & .OR.(equations_set%SPECIFICATION(3)==equations_set_ale_pressure_poisson_subtype) &
2366  & .OR.(equations_set%SPECIFICATION(3)==equations_set_nonlinear_pressure_poisson_subtype))) THEN
2367  WRITE(14,'(" ", es25.16 )')nodepvalue(i)
2368  WRITE(14,'(" ", es25.16 )')nodelabelvalue(i)
2369  END IF
2370  END IF
2371 
2372  !---tob: Mass increase for coupled elasticity Darcy INRIA model
2373  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
2374  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) &
2375  & .AND.(equations_set%SPECIFICATION(3)==equations_set_elasticity_darcy_inria_model_subtype) )THEN
2376  IF(numberofdimensions==3)THEN
2377  WRITE(14,'(" ", es25.16 )')nodemivalue(i)
2378  END IF
2379  END IF
2380  !---toe
2381 
2382  IF(output_material) THEN
2383  WRITE(14,'(" ", es25.16 )')nodemuvalue(i)
2384  WRITE(14,'(" ", es25.16 )')noderhovalue(i)
2385 
2386  IF( (equations_set%specification(1)==equations_set_fluid_mechanics_class) &
2387  & .AND.(equations_set%specification(2)==equations_set_transient1d_navier_stokes_subtype)) THEN
2388  WRITE(14,'(" ", es25.16 )')nodeevalue(i)
2389  WRITE(14,'(" ", es25.16 )')nodeh0value(i)
2390  WRITE(14,'(" ", es25.16 )')nodea0value(i)
2391  WRITE(14,'(" ", es25.16 )')nodea1value(i)
2392  WRITE(14,'(" ", es25.16 )')nodea2value(i)
2393  END IF
2394 
2395  IF(equations_set%specification(1)==equations_set_elasticity_class)THEN
2396  IF(equations_set%specification(2)==equations_set_finite_elasticity_type)THEN
2397  IF( (equations_set%specification(3)==equations_set_compressible_finite_elasticity_subtype) &
2398  & .OR.(equations_set%specification(3)==equations_set_elasticity_darcy_inria_model_subtype) &
2399  & .OR.(equations_set%specification(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) &
2400  & .OR. (equations_set%specification(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) )THEN
2401  WRITE(14,'(" ", es25.16 )')nodekappavalue(i)
2402  END IF
2403  END IF
2404  END IF
2405 
2406  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
2407  & .AND.(equations_set%SPECIFICATION(2)==equations_set_darcy_equation_type) &
2408  & .AND.(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) )THEN
2409  !--- Remaining tensor material data of permeability tensor
2410  WRITE(14,'(" ", es25.16 )')nodeperm2value(i)
2411  WRITE(14,'(" ", es25.16 )')nodeperm3value(i)
2412  WRITE(14,'(" ", es25.16 )')nodeperm4value(i)
2413  WRITE(14,'(" ", es25.16 )')nodeperm5value(i)
2414  WRITE(14,'(" ", es25.16 )')nodeperm6value(i)
2415 
2416  !source field
2417  IF( output_source ) THEN
2418  WRITE(14,'(" ", es25.16 )')nodesourcevalue1(i)
2419  WRITE(14,'(" ", es25.16 )')nodesourcevalue2(i)
2420  WRITE(14,'(" ", es25.16 )')nodesourcevalue3(i)
2421  WRITE(14,'(" ", es25.16 )')nodesourcevalue4(i)
2422  END IF
2423  ENDIF
2424  END IF
2425 
2426 
2427  IF( analytic ) THEN
2428  WRITE(14,'(" ", es25.16 )')nodeuvalue_analytic(i)
2429  IF(numberofdimensions==2 .OR. numberofdimensions==3) THEN
2430  WRITE(14,'(" ", es25.16 )')nodevvalue_analytic(i)
2431  END IF
2432  IF(numberofdimensions==3) THEN
2433  WRITE(14,'(" ", es25.16 )')nodewvalue_analytic(i)
2434  END IF
2435  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
2436  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
2437  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
2438  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
2439  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype) )THEN
2440  WRITE(14,'(" ", es25.16 )')nodepvalue_analytic(i)
2441  END IF
2442 
2443  !Calculate Analytic/Numeric error values
2448 
2449  WRITE(14,'(" ", es25.16 )')nodeuvalue_error(i)
2450  IF(numberofdimensions==2 .OR. numberofdimensions==3) THEN
2451  WRITE(14,'(" ", es25.16 )')nodevvalue_error(i)
2452  END IF
2453  IF(numberofdimensions==3) THEN
2454  WRITE(14,'(" ", es25.16 )')nodewvalue_error(i)
2455  END IF
2456  IF( (equations_set%SPECIFICATION(1)==equations_set_fluid_mechanics_class) &
2457  & .OR.(equations_set%SPECIFICATION(1)==equations_set_elasticity_class) &
2458  & .AND.(equations_set%SPECIFICATION(2)==equations_set_finite_elasticity_type) &
2459  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_compressible_finite_elasticity_subtype) &
2460  & .AND.(equations_set%SPECIFICATION(3).NE.equations_set_elasticity_darcy_inria_model_subtype) )THEN
2461  WRITE(14,'(" ", es25.16 )')nodepvalue_error(i)
2462  END IF
2463  END IF
2464 
2465  END DO
2466 
2467  WRITE(14,*) ' '
2468  CLOSE(unit=14)
2469 
2470  IF( darcy%ANALYTIC ) THEN
2472  END IF
2473 
2474 ! !test output for DN only
2475 ! IF(NumberOfDimensions==2 .OR. NumberOfDimensions==3) THEN
2476 ! IF(DN) THEN
2477 ! FILENAME="./output/"//NAME//".davidn"
2478 ! OPEN(UNIT=14, FILE=CHAR(FILENAME),STATUS='unknown')
2479 ! WRITE(14,*) NodesPerMeshComponent(1),NodesPerMeshComponent(1),NodesPerMeshComponent(2)
2480 ! DO I=1,NodesPerMeshComponent(1)
2481 ! WRITE(14,'(3(" ", es25.16 ))')NodeXValue(I),NodeYValue(I),NodeZValue(I)
2482 ! ENDDO
2483 ! DO I=1,NodesPerMeshComponent(1)
2484 ! WRITE(14,'(6(" ", es25.16 ))')NodeXValue(I),NodeYValue(I),NodeZValue(I),NodeUValue(I),NodeVValue(I),NodeWValue(I)
2485 ! ENDDO
2486 ! DO I=1,NodesPerMeshComponent(2)
2487 ! WRITE(14,'(6(" ", es25.16 ))')NodeXValue(I),NodeYValue(I),NodeZValue(I),NodePValue2(I)
2488 ! ENDDO
2489 ! CLOSE(14)
2490 ! ENDIF
2491 ! END IF
2492 
2493 ! IF(NumberOfDimensions==1) THEN
2494 ! IF(DN) THEN
2495 ! FILENAME="./output/"//NAME//".davidn"
2496 ! OPEN(UNIT=14, FILE=CHAR(FILENAME),STATUS='unknown')
2497 ! WRITE(14,*) NodesPerMeshComponent(1),NodesPerMeshComponent(1)
2498 ! DO I=1,NodesPerMeshComponent(1)
2499 ! WRITE(14,'(3(" ", es25.16 ))')NodeXValue(I),NodeYValue(I),NodeZValue(I)
2500 ! ENDDO
2501 ! DO I=1,NodesPerMeshComponent(1)
2502 ! WRITE(14,'(6(" ", es25.16 ))')NodeXValue(I),NodeYValue(I),NodeZValue(I),NodeUValue(I),NodeVValue(I),NodeWValue(I)
2503 ! ENDDO
2504 ! CLOSE(14)
2505 ! ENDIF
2506 ! END IF
2507 
2508  CALL write_string(general_output_type,"Writing Nodes...",err,error,*999)
2509  RETURN
2510 999 errorsexits("FLUID_MECHANICS_IO_WRITE_NODES_CMGUI",err,error)
2511  RETURN
2512 
2514 
2515  ! OK
2516  !================================================================================================================================
2517  !
2518 
2522 ! TYPE(VARYING_STRING), INTENT(IN) :: NAME !<the prefix name of file.
2523  CHARACTER(len=14), INTENT(IN) :: NAME
2524  TYPE(varying_string) :: FILENAME
2525 ! CHARACTER :: FILENAME !<the prefix name of file.
2526  ! CHARACTER*60 ELEM_TYPE
2527  INTEGER(INTG):: I,J,K,KK
2528  INTEGER(INTG) :: ERR,ierror
2529  TYPE(varying_string):: ERROR
2530  LOGICAL:: OUTPUT_FLAG
2531 
2532  filename="./output/"//name//".exelem"
2533  OPEN(unit=5, file=char(filename),status='unknown', iostat=ierror)
2534  WRITE(5,*) 'Group name: OpenCMISS'
2535 
2536 
2537  DO kk = 1,numberofelements
2538 
2539  !write out the running header only once #NodesPerElement is changing
2540  IF( kk == 1 ) THEN
2541  output_flag = .true.
2542  ELSE IF( kk > 1 .AND. nodesperelement(kk).NE.nodesperelement(kk-1) ) THEN
2543  output_flag = .true.
2544  ELSE
2545  output_flag = .false.
2546  END IF
2547 
2548  IF( output_flag ) THEN
2549 
2550  IF(lagrange_simplex==2) THEN
2551  IF(numberofdimensions==2) THEN
2552  WRITE(5,*) 'Shape. Dimension=',trim(nms(numberofdimensions)),', simplex(2)*simplex'
2553  IF(maxnodesperelement==3) THEN
2554  WRITE(5,*) '#Scale factor sets= 1'
2555 ! WRITE(5,*) ' l.simplex(2)*l.simplex, #Scale factors= ', NodesPerElement(1)
2556  WRITE(5,*) ' l.simplex(2)*l.simplex, #Scale factors= ', nodesperelement(kk)
2557  ELSE IF(maxnodesperelement==6) THEN
2558  WRITE(5,*) '#Scale factor sets= 1'
2559 ! WRITE(5,*) ' l.simplex(2)*l.simplex, #Scale factors= ', NodesPerElement(1)
2560  WRITE(5,*) ' l.simplex(2)*l.simplex, #Scale factors= ', nodesperelement(kk)
2561  ELSE IF (maxnodesperelement== 10 ) THEN
2562  WRITE(5,*) '#Scale factor sets= 1'
2563 ! WRITE(5,*) ' q.simplex(2)*q.simplex, #Scale factors= ', NodesPerElement(1)
2564  WRITE(5,*) ' q.simplex(2)*q.simplex, #Scale factors= ', nodesperelement(kk)
2565  ENDIF
2566  ELSE IF(numberofdimensions==3) THEN
2567  WRITE(5,*) 'Shape. Dimension=',trim(nms(numberofdimensions)),', simplex(2;3)*simplex*simplex'
2568  IF(maxnodesperelement==4) THEN
2569  WRITE(5,*) '#Scale factor sets= 1'
2570 ! WRITE(5,*) ' l.simplex(2;3)*l.simplex*l.simplex, #Scale factors= ', NodesPerElement(1)
2571  WRITE(5,*) ' l.simplex(2;3)*l.simplex*l.simplex, #Scale factors= ', nodesperelement(kk)
2572  ELSE IF (maxnodesperelement== 10 ) THEN
2573  WRITE(5,*) '#Scale factor sets= 1'
2574 ! WRITE(5,*) ' q.simplex(2;3)*q.simplex*q.simplex, #Scale factors= ', NodesPerElement(1)
2575  WRITE(5,*) ' q.simplex(2;3)*q.simplex*q.simplex, #Scale factors= ', nodesperelement(kk)
2576  ELSE IF(maxnodesperelement==20) THEN
2577  WRITE(5,*) '#Scale factor sets= 1'
2578 ! WRITE(5,*) ' q.simplex(2;3)*q.simplex*q.simplex, #Scale factors= ', NodesPerElement(1)
2579  WRITE(5,*) ' q.simplex(2;3)*q.simplex*q.simplex, #Scale factors= ', nodesperelement(kk)
2580  ENDIF
2581  ELSE
2582  WRITE(5,*) '#Scale factor sets= 0'
2583  END IF
2584  ELSE IF (lagrange_simplex==1) THEN
2585  WRITE(5,*) 'Shape. Dimension= ',trim(nms(numberofdimensions))
2586  WRITE(5,*) '#Scale factor sets= 1'
2587  IF(numberofdimensions==1) THEN
2588 ! WRITE(5,*) 'q.Lagrange, #Scale factors=',NodesPerElement(1)
2589  WRITE(5,*) 'q.Lagrange, #Scale factors=',nodesperelement(kk)
2590  ELSE IF (numberofdimensions==2) THEN
2591  IF(maxnodesperelement==4) THEN
2592 ! WRITE(5,*) 'l.Lagrange*l.Lagrange, #Scale factors=',NodesPerElement(1)
2593  WRITE(5,*) 'l.Lagrange*l.Lagrange, #Scale factors=',nodesperelement(kk)
2594  ELSE IF(maxnodesperelement==9) THEN
2595 ! WRITE(5,*) 'q.Lagrange*q.Lagrange, #Scale factors=',NodesPerElement(1)
2596  WRITE(5,*) 'q.Lagrange*q.Lagrange, #Scale factors=',nodesperelement(kk)
2597  ELSE IF(maxnodesperelement==16) THEN
2598 ! WRITE(5,*) 'c.Lagrange*c.Lagrange, #Scale factors=',NodesPerElement(1)
2599  WRITE(5,*) 'c.Lagrange*c.Lagrange, #Scale factors=',nodesperelement(kk)
2600  END IF
2601  ELSE
2602  IF(maxnodesperelement==8) THEN
2603 ! WRITE(5,*) 'l.Lagrange*l.Lagrange*l.Lagrange, #Scale factors=',NodesPerElement(1)
2604  WRITE(5,*) 'l.Lagrange*l.Lagrange*l.Lagrange, #Scale factors=',nodesperelement(kk)
2605  ELSE IF(maxnodesperelement==27) THEN
2606 ! WRITE(5,*) 'q.Lagrange*q.Lagrange*q.Lagrange, #Scale factors=',NodesPerElement(1)
2607  WRITE(5,*) 'q.Lagrange*q.Lagrange*q.Lagrange, #Scale factors=',nodesperelement(kk)
2608  ELSE IF(maxnodesperelement==64) THEN
2609 ! WRITE(5,*) 'c.Lagrange*c.Lagrange*c.Lagrange, #Scale factors=',NodesPerElement(1)
2610  WRITE(5,*) 'c.Lagrange*c.Lagrange*c.Lagrange, #Scale factors=',nodesperelement(kk)
2611  END IF
2612  END IF
2613  END IF
2614 
2615 ! WRITE(5,*) '#Nodes= ',TRIM(NMs(NodesPerElement(1)))
2616  WRITE(5,*) '#Nodes= ',trim(nms(nodesperelement(kk)))
2617  WRITE(5,*) '#Fields= ',trim(nms(numberoffields))
2618 
2619  DO i=1,numberoffields
2620  IF(i==1)THEN
2621  WRITE(5,*)' 1) coordinates, coordinate, rectangular cartesian, #Components= ',trim(nms(numberofdimensions))
2622  ELSE IF(i==2) THEN
2623  WRITE(5,*)' 2) general, field, rectangular cartesian, #Components= ',trim(nms(numberofvariablecomponents))
2624  ELSE IF(i==3) THEN
2625  WRITE(5,*)' 3) material, field, rectangular cartesian, #Components= ',trim(nms(numberofmaterialcomponents))
2626  ELSE IF(i==4) THEN
2627  WRITE(5,*)' 4) source, field, rectangular cartesian, #Components= ',trim(nms(numberofsourcecomponents))
2628  END IF
2629 
2630  DO j=1,numberoffieldcomponent(i)
2631  IF(numberofdimensions==1) THEN
2632  IF(i==1)THEN
2633  IF(j==1) THEN
2634  WRITE(5,*)' x. q.Lagrange, no modify, standard node based.'
2635  ELSE IF(j==2) THEN
2636  WRITE(5,*)' y. q.Lagrange, no modify, standard node based.'
2637  ELSE IF(j==3) THEN
2638  WRITE(5,*)' z. q.Lagrange, no modify, standard node based.'
2639  END IF
2640  ELSE
2641  WRITE(5,*)' ',trim(nms(j)),'. q.Lagrange, no modify, standard node based.'
2642  END IF
2643  ELSE IF(numberofdimensions==2) THEN
2644  IF(i==1)THEN
2645  IF(j==1) THEN
2646  IF(maxnodesperelement==4)THEN
2647  WRITE(5,*)' x. l.Lagrange*l.Lagrange, no modify, standard node based.'
2648  ELSE IF(maxnodesperelement==9) THEN
2649  WRITE(5,*)' x. q.Lagrange*q.Lagrange, no modify, standard node based.'
2650  ELSE IF(maxnodesperelement==16) THEN
2651  WRITE(5,*)' x. c.Lagrange*c.Lagrange, no modify, standard node based.'
2652 
2653  ELSE IF(maxnodesperelement==3) THEN
2654  WRITE(5,*)' x. l.simplex(2)*l.simplex, no modify, standard node based.'
2655  ELSE IF(maxnodesperelement==6) THEN
2656  WRITE(5,*)' x. q.simplex(2)*q.simplex, no modify, standard node based.'
2657  ELSE IF(maxnodesperelement==10) THEN
2658  WRITE(5,*)' x. c.simplex(2)*c.simplex, no modify, standard node based.'
2659  END IF
2660  ELSE IF(j==2) THEN
2661  IF(maxnodesperelement==4) THEN
2662  WRITE(5,*)' y. l.Lagrange*l.Lagrange, no modify, standard node based.'
2663  ELSE IF(maxnodesperelement==9) THEN
2664  WRITE(5,*)' y. q.Lagrange*q.Lagrange, no modify, standard node based.'
2665  ELSE IF(maxnodesperelement==16) THEN
2666  WRITE(5,*)' y. c.Lagrange*c.Lagrange, no modify, standard node based.'
2667  ELSE IF(maxnodesperelement==3) THEN
2668  WRITE(5,*)' y. l.simplex(2)*l.simplex, no modify, standard node based.'
2669  ELSE IF(maxnodesperelement==6) THEN
2670  WRITE(5,*)' y. q.simplex(2)*q.simplex, no modify, standard node based.'
2671  ELSE IF(maxnodesperelement==10) THEN
2672  WRITE(5,*)' y. c.simplex(2)*c.simplex, no modify, standard node based.'
2673  END IF
2674  ELSE IF(j==3) THEN
2675  IF(maxnodesperelement==4) THEN
2676  WRITE(5,*)' z. l.Lagrange*l.Lagrange, no modify, standard node based.'
2677  ELSE IF(maxnodesperelement==9) THEN
2678  WRITE(5,*)' z. q.Lagrange*q.Lagrange, no modify, standard node based.'
2679  ELSE IF(maxnodesperelement==16) THEN
2680  WRITE(5,*)' z. c.Lagrange*c.Lagrange, no modify, standard node based.'
2681  ELSE IF(maxnodesperelement==3) THEN
2682  WRITE(5,*)' z. l.simplex(2)*l.simplex, no modify, standard node based.'
2683  ELSE IF(maxnodesperelement==6) THEN
2684  WRITE(5,*)' z. q.simplex(2)*q.simplex, no modify, standard node based.'
2685  ELSE IF(maxnodesperelement==10) THEN
2686  WRITE(5,*)' z. c.simplex(2)*c.simplex, no modify, standard node based.'
2687  END IF
2688  END IF
2689  ELSE
2690  IF(maxnodesperelement==4) THEN
2691  WRITE(5,*)' ',trim(nms(j)),'. l.Lagrange*l.Lagrange, no modify, standard node based.'
2692  ELSE IF(maxnodesperelement==9) THEN
2693  WRITE(5,*)' ',trim(nms(j)),'. q.Lagrange*q.Lagrange, no modify, standard node based.'
2694  ELSE IF(maxnodesperelement==16) THEN
2695  WRITE(5,*)' ',trim(nms(j)),'. c.Lagrange*c.Lagrange, no modify, standard node based.'
2696  ELSE IF(maxnodesperelement==3) THEN
2697  WRITE(5,*)' ',trim(nms(j)),'. l.simplex(2)*l.simplex, no modify, standard node based.'
2698  ELSE IF(maxnodesperelement==6) THEN
2699  WRITE(5,*)' ',trim(nms(j)),'. q.simplex(2)*q.simplex, no modify, standard node based.'
2700  ELSE IF(maxnodesperelement==10) THEN
2701  WRITE(5,*)' ',trim(nms(j)),'. c.simplex(2)*c.simplex, no modify, standard node based.'
2702  END IF
2703  END IF
2704  ELSE IF(numberofdimensions==3) THEN
2705  IF(i==1)THEN
2706  IF(j==1) THEN
2707  IF(maxnodesperelement==8) THEN
2708  WRITE(5,*)' x. l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.'
2709  ELSE IF(maxnodesperelement==27) THEN
2710  WRITE(5,*)' x. q.Lagrange*q.Lagrange*q.Lagrange, no modify, standard node based.'
2711  ELSE IF(maxnodesperelement==64) THEN
2712  WRITE(5,*)' x. c.Lagrange*c.Lagrange*c.Lagrange, no modify, standard node based.'
2713  ELSE IF(maxnodesperelement==4) THEN
2714  WRITE(5,*)' x. l.simplex(2;3)*l.simplex*l.simplex, no modify, standard node based.'
2715  ELSE IF(maxnodesperelement==10) THEN
2716  WRITE(5,*)' x. q.simplex(2;3)*q.simplex*q.simplex, no modify, standard node based.'
2717  ELSE IF(maxnodesperelement==20) THEN
2718  WRITE(5,*)' x. c.simplex(2;3)*c.simplex*c.simplex, no modify, standard node based.'
2719  END IF
2720  ELSE IF(j==2) THEN
2721  IF(maxnodesperelement==8) THEN
2722  WRITE(5,*)' y. l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.'
2723  ELSE IF(maxnodesperelement==27) THEN
2724  WRITE(5,*)' y. q.Lagrange*q.Lagrange*q.Lagrange, no modify, standard node based.'
2725  ELSE IF(maxnodesperelement==64) THEN
2726  WRITE(5,*)' y. c.Lagrange*c.Lagrange*c.Lagrange, no modify, standard node based.'
2727  ELSE IF(maxnodesperelement==4) THEN
2728  WRITE(5,*)' y. l.simplex(2;3)*l.simplex*l.simplex, no modify, standard node based.'
2729  ELSE IF(maxnodesperelement==10) THEN
2730  WRITE(5,*)' y. q.simplex(2;3)*q.simplex*q.simplex, no modify, standard node based.'
2731  ELSE IF(maxnodesperelement==20) THEN
2732  WRITE(5,*)' y. c.simplex(2;3)*c.simplex*c.simplex, no modify, standard node based.'
2733  END IF
2734  ELSE IF(j==3) THEN
2735  IF(maxnodesperelement==8) THEN
2736  WRITE(5,*)' z. l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.'
2737  ELSE IF(maxnodesperelement==27) THEN
2738  WRITE(5,*)' z. q.Lagrange*q.Lagrange*q.Lagrange, no modify, standard node based.'
2739  ELSE IF(maxnodesperelement==64) THEN
2740  WRITE(5,*)' z. c.Lagrange*c.Lagrange*c.Lagrange, no modify, standard node based.'
2741  ELSE IF(maxnodesperelement==4) THEN
2742  WRITE(5,*)' z. l.simplex(2;3)*l.simplex*l.simplex, no modify, standard node based.'
2743  ELSE IF(maxnodesperelement==10) THEN
2744  WRITE(5,*)' z. q.simplex(2;3)*q.simplex*q.simplex, no modify, standard node based.'
2745  ELSE IF(maxnodesperelement==20) THEN
2746  WRITE(5,*)' z. c.simplex(2;3)*c.simplex*c.simplex, no modify, standard node based.'
2747  END IF
2748  END IF
2749  ELSE
2750  IF(maxnodesperelement==8) THEN
2751  WRITE(5,*)' ',trim(nms(j)),'. l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based.'
2752  ELSE IF(maxnodesperelement==27) THEN
2753  WRITE(5,*)' ',trim(nms(j)),'. q.Lagrange*q.Lagrange*q.Lagrange, no modify, standard node based.'
2754  ELSE IF(maxnodesperelement==64) THEN
2755  WRITE(5,*)' ',trim(nms(j)),'. c.Lagrange*c.Lagrange*c.Lagrange, no modify, standard node based.'
2756  ELSE IF(maxnodesperelement==4) THEN
2757  WRITE(5,*)' ',trim(nms(j)),'. l.simplex(2;3)*l.simplex*l.simplex, no modify, standard node based.'
2758  ELSE IF(maxnodesperelement==10) THEN
2759  WRITE(5,*)' ',trim(nms(j)),'. q.simplex(2;3)*q.simplex*q.simplex, no modify, standard node based.'
2760  ELSE IF(maxnodesperelement==20) THEN
2761  WRITE(5,*)' ',trim(nms(j)),'. c.simplex(2;3)*c.simplex*c.simplex, no modify, standard node based.'
2762  END IF
2763  END IF
2764  END IF
2765 
2766  WRITE(5,*) ' #Nodes= ',trim(nms(maxnodesperelement))
2767 
2768  DO k = 1,maxnodesperelement
2769  WRITE(5,*) ' ',trim(nms(k)),'. #Values=1'
2770  WRITE(5,*) ' Value indices: 1'
2771  WRITE(5,*) ' Scale factor indices: ',trim(nms(k))
2772  END DO
2773  END DO
2774  END DO
2775 
2776  END IF !write out the running header only once #NodesPerElement is changing
2777 
2778 
2779  IF(lagrange_simplex==2) THEN
2780 ! IF(.NOT.ALLOCATED(SimplexOutputHelp)) ALLOCATE(SimplexOutputHelp(NodesPerElement(1)))
2781  IF(.NOT.ALLOCATED(simplexoutputhelp)) ALLOCATE(simplexoutputhelp(nodesperelement(kk)))
2782 
2783 ! DO K = 1,NumberOfElements
2784  k = kk
2785 
2786  IF(numberofdimensions==2)THEN
2793  ELSE IF(numberofdimensions==3) THEN
2804  END IF
2805 
2806  WRITE(5,*) 'Element: ', k,' 0 0'
2807  WRITE(5,*) ' Nodes:'
2808  WRITE(5,*) ' ', simplexoutputhelp
2809  WRITE(5,*) ' Scale factors:'
2810 ! WRITE(5,*) ' ',ElementNodesScales(K,1:NodesPerElement(1))
2811  WRITE(5,*) ' ',elementnodesscales(k,1:nodesperelement(k))
2812 ! END DO
2813 
2814  ELSE IF (lagrange_simplex==1) THEN
2815 
2816 ! DO K = 1,NumberOfElements
2817  k = kk
2818 
2819  WRITE(5,*) 'Element: ', k,' 0 0'
2820  WRITE(5,*) ' Nodes:'
2821 ! WRITE(5,*) ' ', ElementNodes(K,1:NodesPerElement(1))
2822  WRITE(5,*) ' ', elementnodes(k,1:nodesperelement(k))
2823  WRITE(5,*) ' Scale factors:'
2824 ! WRITE(5,*) ' ',ElementNodesScales(K,1:NodesPerElement(1))
2825  WRITE(5,*) ' ',elementnodesscales(k,1:nodesperelement(k))
2826 ! END DO
2827  END IF
2828 
2829 
2830  ENDDO
2831 
2832 
2833 
2834  WRITE(5,*) ' '
2835  CLOSE(unit=5)
2836  CALL write_string(general_output_type,"Writing Elements...",err,error,*999)
2837  RETURN
2838 999 errorsexits("FLUID_MECHANICS_IO_WRITE_ELEMENTS_CMGUI",err,error)
2839  RETURN
2840 
2842 
2843  ! OK
2844  !================================================================================================================================
2845  !
2846 
2848  SUBROUTINE fluid_mechanics_io_read_cmheart1(EXPORT,ERR)
2849  !DLLEXPORT(FLUID_MECHANICS_IO_READ_CMHEART1)
2850 
2851  !Argument variables
2852  type(export_container):: export
2853  INTEGER(INTG) :: ERR
2854  TYPE(varying_string):: ERROR
2855  !Local Variables
2856 ! TYPE (EXPORT_CONTAINER):: TMP
2857  ! INTEGER(INTG):: test
2858 
2859  enters("FLUID_MECHANICS_IO_READ_CMHEART1",err,error,*999)
2860  WRITE(*,*)' '
2861  WRITE(*,*)'Importing CMHEART information...'
2862  WRITE(*,*)' '
2863 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2864 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2865 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Run from bin directory if input needed.",ERR,ERROR,*999)
2866 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2867 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"REMEMBER: M,V,P need to be defined as required by cmHeart!",ERR,ERROR,*999)
2868 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2869 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Press ENTER to start.",ERR,ERROR,*999)
2870 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2871 ! READ(*,*)
2872  OPEN(unit=42, file='./input/CMHEART.inp',status='old')
2873 
2877 
2878  CLOSE(42)
2879 
2881  & stat=alloc_error)
2883  & stat=alloc_error)
2885  & stat=alloc_error)
2886 
2889  & numberofnodesperelement(1),1)
2891  & numberofnodesperelement(2),2)
2893  & numberofnodesperelement(3),3)
2894  WRITE(*,*)' '
2895  WRITE(*,*)'Import finished successfully...'
2896  WRITE(*,*)' '
2897 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Export finished successfully...",ERR,ERROR,*999)
2898 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2899 
2900  IF(alloc_error.NE.0) THEN
2901  CALL flagerror("Error during allocation.",err,error,*999)
2902  END IF
2903 
2907  ALLOCATE(tmp%N(totalnumberofnodes,3),stat=alloc_error)
2908 
2913  tmp%D=dimen
2914  tmp%F=base_info%n_B
2915  tmp%ID_M=1
2916  tmp%ID_V=2
2917  tmp%ID_P=3
2921 
2922  IF (base_info%HEXA==1) THEN
2923  !LAGRANGIAN BASIS
2924  tmp%IT_T=1
2925  ELSE
2926  ! SIMPLEX BASIS
2927  tmp%IT_T=2
2928  END IF
2929 
2937  tmp%EN_T=tmp%EN_M+tmp%EN_V+tmp%EN_P
2938  tmp%N_M=arrayofnodesdefined(1)
2939  tmp%N_V=arrayofnodesdefined(2)
2940  tmp%N_P=arrayofnodesdefined(3)
2941  tmp%N_T=tmp%N_M+tmp%N_V+tmp%N_P
2942 
2943  export=tmp
2944 
2945  IF(alloc_error.NE.0) THEN
2946  CALL flagerror("Error during allocation.",err,error,*999)
2947  END IF
2948 
2949  exits("FLUID_MECHANICS_IO_READ_CMHEART1")
2950  RETURN
2951 999 errorsexits("FLUID_MECHANICS_IO_READ_CMHEART1",err,error)
2952  RETURN
2953 
2954  END SUBROUTINE fluid_mechanics_io_read_cmheart1
2955 
2956 
2957  ! OK
2958  !================================================================================================================================
2959  !
2960 
2961 
2963  SUBROUTINE fluid_mechanics_io_read_cmheart2(EXPORT,BC,ERR)
2964  !DLLEXPORT(FLUID_MECHANICS_IO_READ_CMHEART2)
2965 
2966  !Argument variables
2967  type(export_container):: export
2968  type(boundary_parameters):: bc
2969  INTEGER(INTG) :: ERR
2970  TYPE(varying_string):: ERROR
2971  !Local Variables
2972 ! TYPE (EXPORT_CONTAINER):: TMP
2973  ! INTEGER(INTG):: test
2974 
2975  enters("FLUID_MECHANICS_IO_READ_CMHEART2",err,error,*999)
2976  WRITE(*,*)' '
2977  WRITE(*,*)'Importing CMHEART information...'
2978  WRITE(*,*)' '
2979 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2980 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2981 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Run from bin directory if input needed.",ERR,ERROR,*999)
2982 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2983 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"REMEMBER: M,V,P need to be defined as required by cmHeart!",ERR,ERROR,*999)
2984 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2985 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Press ENTER to start.",ERR,ERROR,*999)
2986 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
2987 ! READ(*,*)
2988  OPEN(unit=42, file='./input/CMHEART_MESH1.inp',status='old')
2989 
2993 
2994  CLOSE(42)
2995 
2997  & stat=alloc_error)
2999  & stat=alloc_error)
3001  & stat=alloc_error)
3002 
3005  & numberofnodesperelement(1),1)
3007  & numberofnodesperelement(2),2)
3009  & numberofnodesperelement(3),3)
3010  WRITE(*,*)' '
3011  WRITE(*,*)'Import finished successfully...'
3012  WRITE(*,*)' '
3013 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Export finished successfully...",ERR,ERROR,*999)
3014 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
3015 
3016  IF(alloc_error.NE.0) THEN
3017  CALL flagerror("Error during allocation.",err,error,*999)
3018  END IF
3019 
3023  ALLOCATE(tmp%N(totalnumberofnodes,3),stat=alloc_error)
3024 
3029  tmp%D=dimen
3030  tmp%F=base_info%n_B
3031  tmp%ID_M=1
3032  tmp%ID_V=2
3033  tmp%ID_P=3
3037 
3038  IF (base_info%HEXA==1) THEN
3039  !LAGRANGIAN BASIS
3040  tmp%IT_T=1
3041  ELSE
3042  ! SIMPLEX BASIS
3043  tmp%IT_T=2
3044  END IF
3045 
3053  tmp%EN_T=tmp%EN_M+tmp%EN_V+tmp%EN_P
3054  tmp%N_M=arrayofnodesdefined(1)
3055  tmp%N_V=arrayofnodesdefined(2)
3056  tmp%N_P=arrayofnodesdefined(3)
3057  tmp%N_T=tmp%N_M+tmp%N_V+tmp%N_P
3058 
3059  export=tmp
3060 
3061  IF(alloc_error.NE.0) THEN
3062  CALL flagerror("Error during allocation.",err,error,*999)
3063  END IF
3064 
3065  !Now read boundary information from CMHEART
3066  OPEN(unit=42, file='./input/CMHEART_BC.inp',status='old')
3068  bc=bc_tmp
3069  CLOSE(42)
3070 
3071  exits("FLUID_MECHANICS_IO_READ_CMHEART2")
3072  RETURN
3073 999 errorsexits("FLUID_MECHANICS_IO_READ_CMHEART2",err,error)
3074  RETURN
3075 
3076  END SUBROUTINE fluid_mechanics_io_read_cmheart2
3077 
3078 
3079  ! OK
3080  !================================================================================================================================
3081  !
3082 
3083 
3085  SUBROUTINE fluid_mechanics_io_read_cmheart3(EXPORT1,EXPORT2,EXPORT3,CONNECT,BC,ERR)
3086  !DLLEXPORT(FLUID_MECHANICS_IO_READ_CMHEART3)
3087 
3088  !Argument variables
3089  type(export_container):: export1,export2,export3
3090  type(coupling_parameters):: connect
3091  type(boundary_parameters):: bc
3092  INTEGER(INTG) :: ERR
3093  TYPE(varying_string):: ERROR
3094  !Local Variables
3095 ! TYPE (EXPORT_CONTAINER):: TMP1
3096  INTEGER(INTG):: I
3097 
3098  enters("FLUID_MECHANICS_IO_READ_CMHEART3",err,error,*999)
3099  WRITE(*,*)' '
3100  WRITE(*,*)'Importing CMHEART information...'
3101  WRITE(*,*)' '
3102 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
3103 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
3104 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Run from bin directory if input needed.",ERR,ERROR,*999)
3105 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
3106 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"REMEMBER: M,V,P need to be defined as required by cmHeart!",ERR,ERROR,*999)
3107 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
3108 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Press ENTER to start.",ERR,ERROR,*999)
3109 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE," ",ERR,ERROR,*999)
3110 ! READ(*,*)
3111 
3112  DO i=1,3
3113  IF(i==1) THEN
3114  OPEN(unit=42, file='./input/CMHEART_MESH1.inp',status='old')
3115  ELSE IF(i==2) THEN
3116  OPEN(unit=42, file='./input/CMHEART_MESH2.inp',status='old')
3117  ELSE IF(i==3) THEN
3118  OPEN(unit=42, file='./input/CMHEART_INTERFACE.inp',status='old')
3119  ENDIF
3120 
3124 
3125  CLOSE(42)
3126 
3127  IF(ALLOCATED(opencmiss_elem_m)) DEALLOCATE(opencmiss_elem_m)
3128  IF(ALLOCATED(opencmiss_elem_v)) DEALLOCATE(opencmiss_elem_v)
3129  IF(ALLOCATED(opencmiss_elem_p)) DEALLOCATE(opencmiss_elem_p)
3130 
3132  & stat=alloc_error)
3134  & stat=alloc_error)
3136  & stat=alloc_error)
3137 
3140  & numberofnodesperelement(1),1)
3142  & numberofnodesperelement(2),2)
3144  & numberofnodesperelement(3),3)
3145  WRITE(*,*)' '
3146  WRITE(*,*)'Import finished successfully...', i
3147  WRITE(*,*)' '
3148  CALL write_string(general_output_type,"Export finished successfully...",err,error,*999)
3149  CALL write_string(general_output_type," ",err,error,*999)
3150 
3151  IF(alloc_error.NE.0) THEN
3152  CALL flagerror("Error during allocation.",err,error,*999)
3153  END IF
3154 
3158  ALLOCATE(tmp%N(totalnumberofnodes,3),stat=alloc_error)
3159 
3164  tmp%D=dimen
3165  tmp%F=base_info%n_B
3166  tmp%ID_M=1
3167  tmp%ID_V=2
3168  tmp%ID_P=3
3172 
3173  IF (base_info%HEXA==1) THEN
3174  !LAGRANGIAN BASIS
3175  tmp%IT_T=1
3176  ELSE
3177  ! SIMPLEX BASIS
3178  tmp%IT_T=2
3179  END IF
3180 
3188  tmp%EN_T=tmp%EN_M+tmp%EN_V+tmp%EN_P
3189  tmp%N_M=arrayofnodesdefined(1)
3190  tmp%N_V=arrayofnodesdefined(2)
3191  tmp%N_P=arrayofnodesdefined(3)
3192  tmp%N_T=tmp%N_M+tmp%N_V+tmp%N_P
3193 
3194  IF(i==1) THEN
3195  export1=tmp
3196  ELSE IF (i==2) THEN
3197  export2=tmp
3198  ELSE IF (i==3) THEN
3199  export3=tmp
3200  ENDIF
3201 
3202  IF(alloc_error.NE.0) THEN
3203  CALL flagerror("Error during allocation.",err,error,*999)
3204  END IF
3205 
3206  ENDDO
3207 
3208  !Now read connectivity information from CMHEART
3209  OPEN(unit=42, file='./input/CMHEART_CONNECT.inp',status='old')
3211  connect=connect_tmp
3212  CLOSE(42)
3213 
3214 
3215 
3216  !Now read boundary information from CMHEART
3217  OPEN(unit=42, file='./input/CMHEART_BC.inp',status='old')
3219  bc=bc_tmp
3220  CLOSE(42)
3221 
3222 
3223  exits("FLUID_MECHANICS_IO_READ_CMHEART3")
3224  RETURN
3225 999 errorsexits("FLUID_MECHANICS_IO_READ_CMHEART3",err,error)
3226  RETURN
3227 
3228  END SUBROUTINE fluid_mechanics_io_read_cmheart3
3229 
3230  ! OK
3231  !================================================================================================================================
3232  !
3233 
3237  !Argument variables
3238 
3239  !Local Variables
3240 
3241  ! INTEGER(INTG):: test
3242 
3243  WRITE(*,*)' '
3244  WRITE(*,*)'Free all CMHEART information memory...'
3245  WRITE(*,*)' '
3246 
3247  IF(ALLOCATED(opencmiss_elem_m)) DEALLOCATE(opencmiss_elem_m)
3248  IF(ALLOCATED(opencmiss_elem_v)) DEALLOCATE(opencmiss_elem_v)
3249  IF(ALLOCATED(opencmiss_elem_p)) DEALLOCATE(opencmiss_elem_p)
3250  DEALLOCATE(tmp%M)
3251  DEALLOCATE(tmp%V)
3252  DEALLOCATE(tmp%P)
3253  DEALLOCATE(tmp%N)
3254  DEALLOCATE(base_info%B)
3255 
3256  DEALLOCATE(mesh_info(1)%X)
3257  DEALLOCATE(mesh_info(2)%X)
3258  DEALLOCATE(mesh_info(3)%X)
3259  DEALLOCATE(opencmiss_node_coord)
3260 
3261  DEALLOCATE(mesh_info(1)%T)
3262  DEALLOCATE(mesh_info(2)%T)
3263  DEALLOCATE(mesh_info(3)%T)
3264 
3265  exits("FLUID_MECHANICS_IO_READ_CMHEART_FINISH")
3266 
3268 
3269  ! OK
3270  !================================================================================================================================
3271  !
3272 
3274  SUBROUTINE fluid_mechanics_io_ppe_masking(BC_MASKING)
3276  !Argument variables
3277  INTEGER, DIMENSION(TMP%E_T):: ELEMENT_MASKING
3278  REAL, DIMENSION(TMP%N_V):: NODE_MASKING
3279  INTEGER, POINTER:: BC_MASKING(:)
3280  INTEGER:: I,J
3281  INTEGER:: CONTROL_INSIDE, CONTROL_BETWEEN, CONTROL_OUTSIDE
3282  REAL :: DUMMY
3283 
3284  !Local Variables
3285 
3286  ! INTEGER(INTG):: test
3287 
3288  WRITE(*,*)' '
3289  WRITE(*,*)'Perform enhanced PPE masking...'
3290  WRITE(*,*)' '
3291 
3292  bc_masking=1
3293  element_masking=1 !0=outside, 1=inside, 2=between
3294 
3295  !IDENTIFY FLUID DOMAIN / BOUNDARY_CONDITIONS_POISSON
3296  OPEN(unit=42, file="./input/data/ORI_DATA_01.dat",status='unknown')
3297  READ(42,*)dummy
3298  DO i=1,tmp%N_V !NUMBER_OF_NODES_PRESSURE
3299  READ(42,*)node_masking(i)
3300  ENDDO
3301  CLOSE(42)
3302  DO i=1,tmp%E_T !TOTAL_NUMBER_OF_ELEMENTS
3303  DO j=1, tmp%EN_V !NUMBER_OF_ELEMENT_NODES_PRESSURE
3304  IF(node_masking(tmp%M(i,j))<0.5_dp) THEN
3305  element_masking(i)=0
3306  ENDIF
3307  ENDDO
3308  ENDDO
3309  DO i=1,tmp%E_T !TOTAL_NUMBER_OF_ELEMENTS
3310  IF(element_masking(i)==0) THEN
3311  DO j=1, tmp%EN_V !NUMBER_OF_ELEMENT_NODES_PRESSURE
3312  IF(node_masking(tmp%M(i,j))>0.5_dp) THEN
3313  element_masking(i)=2
3314  ENDIF
3315  ENDDO
3316  ENDIF
3317  ENDDO
3318  !CONTROL
3319  control_inside=0
3320  control_between=0
3321  control_outside=0
3322  DO i=1,tmp%E_T !TOTAL_NUMBER_OF_ELEMENTS
3323  IF(element_masking(i)==0) THEN
3324  control_outside=control_outside+1
3325  ELSE IF(element_masking(i)==1) THEN
3326  control_inside=control_inside+1
3327  ELSE IF(element_masking(i)==2) THEN
3328  control_between=control_between+1
3329  ENDIF
3330  ENDDO
3331 ! ! ! WRITE(*,*)'CONTROL_INSIDE ',CONTROL_INSIDE
3332 ! ! ! WRITE(*,*)'CONTROL_OUTSIDE ',CONTROL_OUTSIDE
3333 ! ! ! WRITE(*,*)'CONTROL_BETWEEN ',CONTROL_BETWEEN
3334  !LOOP OVER ALL OUTSIDE ELEMENTS FIRST
3335  DO i=1,tmp%E_T !TOTAL_NUMBER_OF_ELEMENTS
3336  IF(element_masking(i)==0) THEN
3337  DO j=1, tmp%EN_V
3338  bc_masking(tmp%M(i,j))=0
3339  ENDDO
3340  ENDIF
3341  ENDDO
3342  !THEN LOOP OVER ALL "BETWEEN" ELEMENTS
3343  DO i=1,tmp%E_T !TOTAL_NUMBER_OF_ELEMENTS
3344  IF(element_masking(i)==2) THEN
3345  DO j=1, tmp%EN_V
3346  bc_masking(tmp%M(i,j))=2
3347  ENDDO
3348  ENDIF
3349  ENDDO
3350  !CORRECT BC LABEL BY LOOPING OVER ALL INSIDE ELEMENTS NOW
3351  DO i=1,tmp%E_T !TOTAL_NUMBER_OF_ELEMENTS
3352  IF(element_masking(i)==1) THEN
3353  DO j=1, tmp%EN_V
3354  bc_masking(tmp%M(i,j))=1
3355  ENDDO
3356  ENDIF
3357  ENDDO
3358 
3359  exits("FLUID_MECHANICS_IO_PPE_MASKING")
3360 
3361  END SUBROUTINE fluid_mechanics_io_ppe_masking
3362 
3363 
3364  ! OK
3365  !================================================================================================================================
3366  !
3367 
3369  SUBROUTINE fluid_mechanics_io_read_aux
3371  INTEGER(INTG):: I
3372  INTEGER(INTG) :: ERR
3373  TYPE(varying_string):: ERROR
3374 
3375  READ(42,*) nimz
3376  nimz = trim(nimz); base_info%n_B = 0; base_info%HEXA = 0; base_info%DM = 3
3377  OPEN(unit=1,file=nimz,status='old',action='read') ! Read base file for initial parameters
3378 
3379  base_info%n_B=0
3380  base_info%n_pts=0
3381  base_info%VL=0
3382  base_info%n_pts_f=0
3383  base_info%FACES=0
3384  base_info%FNODES=0
3385  base_info%HEXA=0
3386  base_info%DM=0
3387  base_info%TRI_BASIS=0
3388  base_info%TET_BASIS=0
3389  base_info%QUAD_BASIS=0
3390  base_info%HEX_BASIS=0
3391 
3392 
3393  DO WHILE (0 < 1)
3394  READ(1,*,end=50) in_char
3395  IF (index(in_char,'no_fields!') == 1) READ(1,*) base_info%n_B
3396  IF (index(in_char,'no_gauss!') == 1) READ(1,*) base_info%n_pts
3397  IF (index(in_char,'volume!') == 1) READ(1,*) base_info%VL
3398  IF (index(in_char,'no_gauss_f!') == 1) READ(1,*) base_info%n_pts_f
3399  IF (index(in_char,'no_ele_faces!') == 1) READ(1,*) base_info%FACES
3400  IF (index(in_char,'no_ele_nodes_f!') == 1) READ(1,*) base_info%FNODES
3401  IF (index(in_char,'hexa_basis!') == 1) base_info%HEXA = 1
3402  IF (index(in_char,'domain_dimension!') == 1) READ(1,*) base_info%DM
3403  IF (index(in_char,'TRI_BASIS!') == 1) base_info%TRI_BASIS = 1
3404  IF (index(in_char,'TET_BASIS!') == 1) base_info%TET_BASIS = 1
3405  IF (index(in_char,'QUAD_BASIS!') == 1) base_info%QUAD_BASIS = 1
3406  IF (index(in_char,'HEX_BASIS!') == 1) base_info%HEX_BASIS = 1
3407  END DO
3408 
3409 50 CLOSE(1)
3410 
3411  dimen=base_info%DM
3412 
3413  ALLOCATE(base_info%B(base_info%n_B),stat=alloc_error)
3414  OPEN(unit=1,file=nimz,status='old',action='read') ! Read base file for initial parameters
3415 
3416  DO WHILE (0 < 1)
3417  READ(1,*,end=52) in_char
3418  IF (index(in_char,'no_basis_M!') == 1) READ(1,*) base_info%B(1)%n
3419  IF (index(in_char,'no_basis_V!') == 1) READ(1,*) base_info%B(2)%n
3420  IF (index(in_char,'no_basis_P!') == 1) READ(1,*) base_info%B(3)%n
3421  IF (index(in_char,'dim_field_M!') == 1) READ(1,*) base_info%B(1)%DM
3422  IF (index(in_char,'dim_field_V!') == 1) READ(1,*) base_info%B(2)%DM
3423  IF (index(in_char,'dim_field_P!') == 1) READ(1,*) base_info%B(3)%DM
3424  END DO
3425 
3426 52 CLOSE(1)
3427 
3428  IF (base_info%QUAD_BASIS/= 1.AND.base_info%HEX_BASIS /= 1.AND.base_info%TRI_BASIS/= 1.AND.base_info%TET_BASIS /= 1)THEN
3429  CALL write_string(general_output_type,"Cubic Hermite not implemented yet...",err,error,*999)
3430  ELSE
3431 
3432  DO i=1,3
3433  IF (base_info%TRI_BASIS== 1.OR.base_info%TET_BASIS == 1)THEN
3434  IF(dimen==2) THEN
3435  IF(base_info%B(i)%n==3) THEN
3438  ELSE IF(base_info%B(i)%n==6) THEN
3441  ELSE IF(base_info%B(i)%n==10) THEN
3444  ELSE
3445  stop
3446  END IF
3447  ELSE IF(dimen==3) THEN
3448  IF(base_info%B(i)%n==4) THEN
3451  ELSE IF(base_info%B(i)%n==10) THEN
3454  ELSE IF(base_info%B(i)%n==20) THEN
3457  ELSE
3458  stop
3459  END IF
3460  ELSE
3461  stop
3462  END IF
3463  END IF
3464 
3465  IF (base_info%QUAD_BASIS== 1.OR.base_info%HEX_BASIS == 1)THEN
3466  IF(base_info%B(i)%n==2) THEN
3467  !2D/3D LINEAR LAGRANGE
3469  IF(dimen==2) THEN
3471  ELSE IF(dimen==3) THEN
3473  ELSE
3474  stop
3475  END IF
3476  ELSE IF(base_info%B(i)%n==3) THEN
3477  !2D/3D QUADRATIC LAGRANGE
3479  IF(dimen==2) THEN
3481  ELSE IF(dimen==3) THEN
3483  ELSE
3484  stop
3485  END IF
3486  ELSE IF(base_info%B(i)%n==4) THEN
3487  !2D/3D CUBIC LAGRANGE
3489  IF(dimen==2) THEN
3491  ELSE IF(dimen==3) THEN
3493  ELSE
3494  stop
3495  END IF
3496  ELSE
3497  stop
3498  END IF
3499  END IF
3500  END DO
3501  END IF
3502 
3503  IF(alloc_error.NE.0) THEN
3504  stop 'Error during allocation'
3505  END IF
3506  RETURN
3507 999 errorsexits("FLUID_MECHANICS_IO_READ_AUX",err,error)
3508  RETURN
3509 
3510  END SUBROUTINE fluid_mechanics_io_read_aux
3511 
3512  ! OK
3513  !================================================================================================================================
3514  !
3515 
3517  SUBROUTINE fluid_mechanics_io_order_numbering(NEW,OLD,n,m,I)
3519  INTEGER(INTG):: I,J,M,N
3520  INTEGER(INTG)::NEW(n,m),OLD(n,m)
3521 ! INTEGER(INTG) :: ERR
3522 ! TYPE(VARYING_STRING):: ERROR
3523 
3524 ! NEW=OLD
3525  DO j=1,n
3526  IF (base_info%QUAD_BASIS == 1) THEN
3527  IF(base_info%B(i)%n==2) THEN
3528  !2D HEX LINEAR
3529  new(j,1)=old(j,1)
3530  new(j,2)=old(j,2)
3531  new(j,3)=old(j,3)
3532  new(j,4)=old(j,4)
3533  ELSE IF(base_info%B(i)%n==3) THEN
3534  !2D HEX QUADR
3535  new(j,1)=old(j,1)
3536  new(j,2)=old(j,5)
3537  new(j,3)=old(j,2)
3538  new(j,4)=old(j,6)
3539  new(j,5)=old(j,7)
3540  new(j,6)=old(j,8)
3541  new(j,7)=old(j,3)
3542  new(j,8)=old(j,9)
3543  new(j,9)=old(j,4)
3544  ELSE IF(base_info%B(i)%n==4) THEN
3545  !2D HEX CUB
3546  new(j,1)=old(j,1)
3547  new(j,2)=old(j,5)
3548  new(j,3)=old(j,6)
3549  new(j,4)=old(j,2)
3550  new(j,5)=old(j,7)
3551  new(j,6)=old(j,8)
3552  new(j,7)=old(j,9)
3553  new(j,8)=old(j,10)
3554  new(j,9)=old(j,11)
3555  new(j,10)=old(j,12)
3556  new(j,11)=old(j,13)
3557  new(j,12)=old(j,14)
3558  new(j,13)=old(j,3)
3559  new(j,14)=old(j,15)
3560  new(j,15)=old(j,16)
3561  new(j,16)=old(j,4)
3562  ELSE
3563  stop
3564  END IF
3565  ELSE IF (base_info%HEX_BASIS == 1) THEN
3566  IF(base_info%B(i)%n==2) THEN
3567  !3D HEX LINEAR
3568  new(j,1)=old(j,1)
3569  new(j,2)=old(j,2)
3570  new(j,3)=old(j,3)
3571  new(j,4)=old(j,4)
3572  new(j,5)=old(j,5)
3573  new(j,6)=old(j,6)
3574  new(j,7)=old(j,7)
3575  new(j,8)=old(j,8)
3576  ELSE IF(base_info%B(i)%n==3) THEN
3577  !3D HEX QUADR
3578  new(j,1)=old(j,1)
3579  new(j,2)=old(j,9)
3580  new(j,3)=old(j,2)
3581  new(j,4)=old(j,10)
3582  new(j,5)=old(j,11)
3583  new(j,6)=old(j,12)
3584  new(j,7)=old(j,3)
3585  new(j,8)=old(j,13)
3586  new(j,9)=old(j,4)
3587  new(j,10)=old(j,14)
3588  new(j,11)=old(j,15)
3589  new(j,12)=old(j,16)
3590  new(j,13)=old(j,17)
3591  new(j,14)=old(j,18)
3592  new(j,15)=old(j,19)
3593  new(j,16)=old(j,20)
3594  new(j,17)=old(j,21)
3595  new(j,18)=old(j,22)
3596  new(j,19)=old(j,5)
3597  new(j,20)=old(j,23)
3598  new(j,21)=old(j,6)
3599  new(j,22)=old(j,24)
3600  new(j,23)=old(j,25)
3601  new(j,24)=old(j,26)
3602  new(j,25)=old(j,7)
3603  new(j,26)=old(j,27)
3604  new(j,27)=old(j,8)
3605  ELSE IF(base_info%B(i)%n==4) THEN
3606  !3D HEX CUB
3607  new(j,1)=old(j,1)
3608  new(j,2)=old(j,9)
3609  new(j,3)=old(j,10)
3610  new(j,4)=old(j,2)
3611  new(j,5)=old(j,11)
3612  new(j,6)=old(j,12)
3613  new(j,7)=old(j,13)
3614  new(j,8)=old(j,14)
3615  new(j,9)=old(j,15)
3616  new(j,10)=old(j,16)
3617  new(j,11)=old(j,17)
3618  new(j,12)=old(j,18)
3619  new(j,13)=old(j,3)
3620  new(j,14)=old(j,19)
3621  new(j,15)=old(j,20)
3622  new(j,16)=old(j,4)
3623  new(j,17)=old(j,21)
3624  new(j,18)=old(j,22)
3625  new(j,19)=old(j,23)
3626  new(j,20)=old(j,24)
3627  new(j,21)=old(j,25)
3628  new(j,22)=old(j,26)
3629  new(j,23)=old(j,27)
3630  new(j,24)=old(j,28)
3631  new(j,25)=old(j,29)
3632  new(j,26)=old(j,30)
3633  new(j,27)=old(j,31)
3634  new(j,28)=old(j,32)
3635  new(j,29)=old(j,33)
3636  new(j,30)=old(j,34)
3637  new(j,31)=old(j,35)
3638  new(j,32)=old(j,36)
3639  new(j,33)=old(j,37)
3640  new(j,34)=old(j,38)
3641  new(j,35)=old(j,39)
3642  new(j,36)=old(j,40)
3643  new(j,37)=old(j,41)
3644  new(j,38)=old(j,42)
3645  new(j,39)=old(j,43)
3646  new(j,40)=old(j,44)
3647  new(j,41)=old(j,45)
3648  new(j,42)=old(j,46)
3649  new(j,43)=old(j,47)
3650  new(j,44)=old(j,48)
3651  new(j,45)=old(j,49)
3652  new(j,46)=old(j,50)
3653  new(j,47)=old(j,51)
3654  new(j,48)=old(j,52)
3655  new(j,49)=old(j,5)
3656  new(j,50)=old(j,53)
3657  new(j,51)=old(j,54)
3658  new(j,52)=old(j,6)
3659  new(j,53)=old(j,55)
3660  new(j,54)=old(j,56)
3661  new(j,55)=old(j,57)
3662  new(j,56)=old(j,58)
3663  new(j,57)=old(j,59)
3664  new(j,58)=old(j,60)
3665  new(j,59)=old(j,61)
3666  new(j,60)=old(j,62)
3667  new(j,61)=old(j,7)
3668  new(j,62)=old(j,63)
3669  new(j,63)=old(j,64)
3670  new(j,64)=old(j,8)
3671  ELSE
3672  stop
3673  END IF
3674  ELSE IF (base_info%TRI_BASIS == 1) THEN
3675  IF(base_info%B(i)%n==3) THEN
3676  !2D TET LINEAR
3677  new(j,1)=old(j,1)
3678  new(j,2)=old(j,2)
3679  new(j,3)=old(j,3)
3680  ELSE IF(base_info%B(i)%n==6) THEN
3681  !2D TET QUAD
3682  new(j,1)=old(j,1)
3683  new(j,2)=old(j,2)
3684  new(j,3)=old(j,3)
3685  new(j,4)=old(j,4)
3686  new(j,5)=old(j,6)
3687  new(j,6)=old(j,5)
3688  ELSE IF(base_info%B(i)%n==10) THEN
3689  !2D TET CUB
3690  new(j,1)=old(j,1)
3691  new(j,2)=old(j,2)
3692  new(j,3)=old(j,3)
3693  new(j,4)=old(j,4)
3694  new(j,5)=old(j,5)
3695  new(j,6)=old(j,8)
3696  new(j,7)=old(j,9)
3697  new(j,8)=old(j,7)
3698  new(j,9)=old(j,6)
3699  new(j,10)=old(j,10)
3700  ELSE
3701  stop
3702  END IF
3703  ELSE IF (base_info%TET_BASIS == 1) THEN
3704  IF(base_info%B(i)%n==4) THEN
3705  !3D TET LINEAR
3706  new(j,1)=old(j,1)
3707  new(j,2)=old(j,2)
3708  new(j,3)=old(j,3)
3709  new(j,4)=old(j,4)
3710  ELSE IF(base_info%B(i)%n==10) THEN
3711  !3D TET QUAD
3712  new(j,1)=old(j,1)
3713  new(j,2)=old(j,2)
3714  new(j,3)=old(j,3)
3715  new(j,4)=old(j,4)
3716  new(j,5)=old(j,5)
3717  new(j,6)=old(j,6)
3718  new(j,7)=old(j,7)
3719  new(j,8)=old(j,8)
3720  new(j,9)=old(j,10)
3721  new(j,10)=old(j,9)
3722  ELSE IF(base_info%B(i)%n==20) THEN
3723  !3D TET CUB
3724  new(j,1)=old(j,1)
3725  new(j,2)=old(j,2)
3726  new(j,3)=old(j,3)
3727  new(j,4)=old(j,4)
3728  new(j,5)=old(j,5)
3729  new(j,6)=old(j,6)
3730  new(j,7)=old(j,7)
3731  new(j,8)=old(j,8)
3732  new(j,9)=old(j,9)
3733  new(j,10)=old(j,10)
3734  new(j,11)=old(j,11)
3735  new(j,12)=old(j,12)
3736  new(j,13)=old(j,15)
3737  new(j,14)=old(j,16)
3738  new(j,15)=old(j,13)
3739  new(j,16)=old(j,14)
3740  new(j,17)=old(j,17)
3741  new(j,18)=old(j,18)
3742  new(j,19)=old(j,19)
3743  new(j,20)=old(j,20)
3744  ELSE
3745  stop
3746  END IF
3747  ELSE
3748  stop
3749  END IF
3750  END DO
3751 ! RETURN
3752 ! 999 ERRORSEXITS("FLUID_MECHANICS_IO_ORDER_NUMBERING",ERR,ERROR)
3753 ! RETURN
3754  END SUBROUTINE fluid_mechanics_io_order_numbering
3755 
3756  ! OK
3757  !================================================================================================================================
3758  !
3759 
3763 ! INTEGER(INTG) :: ERR
3764 ! TYPE(VARYING_STRING):: ERROR
3765 
3766  mesh_info(1)%T=mesh_info(1)%T
3769 
3771  ! copy all node numbers from 2 -> 1
3772  mesh_info(2)%T(:,:)= mesh_info(1)%T(:,:)
3773  ELSE
3774  IF (base_info%TRI_BASIS== 1) THEN
3775  mesh_info(2)%T(:,1:3)=mesh_info(1)%T(:,1:3)
3776  ELSE IF (base_info%TET_BASIS == 1)THEN
3777  mesh_info(2)%T(:,1:4)=mesh_info(1)%T(:,1:4)
3778  ELSE IF (base_info%QUAD_BASIS == 1)THEN
3779  mesh_info(2)%T(:,1:4)=mesh_info(1)%T(:,1:4)
3780  ELSE IF (base_info%HEX_BASIS == 1)THEN
3781  mesh_info(2)%T(:,1:8)=mesh_info(1)%T(:,1:8)
3782  ELSE
3783  stop
3784  END IF
3785  END IF
3786 
3788  ! copy all node numbers from 2 -> 1
3789  mesh_info(3)%T(:,:)=mesh_info(1)%T(:,:)
3790  ELSE IF(arrayofnodesdefined(2)==arrayofnodesdefined(3)) THEN
3791  mesh_info(3)%T(:,:)=mesh_info(2)%T(:,:)
3792  ELSE
3793  IF (base_info%TRI_BASIS== 1) THEN
3794  mesh_info(3)%T(:,1:3)=mesh_info(1)%T(:,1:3)
3795  ELSE IF (base_info%TET_BASIS == 1)THEN
3796  mesh_info(3)%T(:,1:4)=mesh_info(1)%T(:,1:4)
3797  ELSE IF (base_info%QUAD_BASIS == 1)THEN
3798  mesh_info(3)%T(:,1:4)=mesh_info(1)%T(:,1:4)
3799  ELSE IF (base_info%HEX_BASIS == 1)THEN
3800  mesh_info(3)%T(:,1:8)=mesh_info(1)%T(:,1:8)
3801  ELSE
3802  stop
3803  END IF
3804  END IF
3805 
3806 ! RETURN
3807 ! 999 ERRORSEXITS("FLUID_MECHANICS_IO_MAKE_UNIQUE",ERR,ERROR)
3808 ! RETURN
3809 
3810  END SUBROUTINE fluid_mechanics_io_make_unique
3811 
3812  ! OK
3813  !================================================================================================================================
3814  !
3815 
3819  INTEGER(INTG):: I,J
3820  INTEGER(INTG):: a,b
3821 ! INTEGER(INTG) :: ERR
3822 ! TYPE(VARYING_STRING):: ERROR
3823  REAL(DP) :: TEMP(3)
3824 ! REAL(DP),DIMENSION(832,3):: sebo_test_array
3825 ! sebo_test_array=0.0
3826 
3827  READ(42,*) namz
3828  OPEN(unit = 1, file=namz,status='old')
3829  READ(1,*) arrayofnodesdefined(1:3)
3830 
3832 ! ALLOCATE AND READ MESH NODE INFORMATION
3833 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Reading Nodes...",ERR,ERROR,*999)
3834  WRITE(*,*)'Reading Nodes...'
3835  DO i=1,3
3837  ALLOCATE(mesh_info(i)%X(mesh_info(i)%Lx,3),stat=alloc_error)
3838  DO j = 1,mesh_info(i)%Lx
3839  READ(1,*,end=35) temp(1:3)
3840  mesh_info(i)%X(j,1:3)=temp(1:3)
3841 ! WRITE(*,*) MESH_INFO(I)%X(J,1:3)
3842 ! READ(1,*,END=35) sebo_test_array(J,1:3)
3843 ! sebo_test_array(J,1:3)=[1,2,3]
3844  END DO
3845  END DO
3846  CLOSE(1)
3847 
3848  IF(ALLOCATED(opencmiss_node_coord)) DEALLOCATE(opencmiss_node_coord)
3849  IF(.NOT.ALLOCATED(opencmiss_node_coord)) ALLOCATE(opencmiss_node_coord(totalnumberofnodes,3),stat=alloc_error)
3850  a=1
3851  b=0
3852  DO i=1,3
3853  a=b+1
3854  b=b+arrayofnodesdefined(i)
3855  opencmiss_node_coord(a:b,1:3)=mesh_info(i)%X(1:arrayofnodesdefined(i),1:3)
3856  END DO
3857  RETURN
3858 
3859 35 print *, 'FAILS'
3860  stop
3861  IF(alloc_error.NE.0) THEN
3862  stop 'Error during allocation'
3863  END IF
3864 
3865  RETURN
3866 ! 999 ERRORSEXITS("FLUID_MECHANICS_IO_READ_NODES",ERR,ERROR)
3867  RETURN
3868 
3869  END SUBROUTINE fluid_mechanics_io_read_nodes
3870 
3871  ! OK
3872  !================================================================================================================================
3873  !
3874 
3878  INTEGER(INTG):: I
3879 
3880  !Read in the interface connectivity mapping
3881 
3882  READ(42,*) nimz
3883  nimz = trim(nimz)
3884  OPEN(unit=79,file=nimz,status='old',action='read') ! Read base file for initial parameters
3885 
3886  READ(79,*) connect_tmp%NUMBER_OF_COUPLINGS
3887  ALLOCATE(connect_tmp%INTERFACE_ELEMENT_NUMBER(connect_tmp%NUMBER_OF_COUPLINGS))
3888  ALLOCATE(connect_tmp%INTERFACE_ELEMENT_LOCAL_NODE(connect_tmp%NUMBER_OF_COUPLINGS))
3889  ALLOCATE(connect_tmp%MESH1_ELEMENT_NUMBER(connect_tmp%NUMBER_OF_COUPLINGS))
3890  ALLOCATE(connect_tmp%MESH1_ELEMENT_XI(connect_tmp%NUMBER_OF_COUPLINGS,3))
3891  ALLOCATE(connect_tmp%MESH2_ELEMENT_NUMBER(connect_tmp%NUMBER_OF_COUPLINGS))
3892  ALLOCATE(connect_tmp%MESH2_ELEMENT_XI(connect_tmp%NUMBER_OF_COUPLINGS,3))
3893 
3894  DO i=1,connect_tmp%NUMBER_OF_COUPLINGS
3895  READ(79,*) connect_tmp%INTERFACE_ELEMENT_NUMBER(i), &
3896  & connect_tmp%INTERFACE_ELEMENT_LOCAL_NODE(i),connect_tmp%MESH1_ID,connect_tmp%MESH1_ELEMENT_NUMBER(i), &
3897  & connect_tmp%MESH1_ELEMENT_XI(i,1:3), connect_tmp%MESH2_ID,connect_tmp%MESH2_ELEMENT_NUMBER(i), &
3898  & connect_tmp%MESH2_ELEMENT_XI(i,1:3)
3899  ENDDO
3900  WRITE(*,*)' '
3901  WRITE(*,*)'Import finished successfully... IMC'
3902  WRITE(*,*)' '
3903  CLOSE(79)
3904 
3905  END SUBROUTINE fluid_mechanics_io_read_connect
3906 
3907 
3908 
3909  ! OK
3910  !================================================================================================================================
3911  !
3912 
3913 
3915  SUBROUTINE fluid_mechanics_io_read_bc
3917  INTEGER(INTG):: NumberOfFaces, NumberOfNodes, HighestNodeNumber
3918  INTEGER(INTG):: FLAG,I1,I2,I3,I4,I5,I,J
3919  LOGICAL :: TEST
3920 
3921  INTEGER(INTG), DIMENSION(:), ALLOCATABLE:: DUMMY,TMP1,TMP2,TMP3,TMP4,TMP5
3922 
3923  READ(42,*) namz
3924  OPEN(unit = 79, file=namz,status='old')
3925  READ(79,*) numberoffaces, numberofnodes
3926 
3927  ALLOCATE(dummy(numberofnodes+1),stat=alloc_error)
3928 
3929  bc_tmp%NUMBER_OF_BC=0
3930  bc_tmp%NUMBER_OF_BC1=0
3931  bc_tmp%NUMBER_OF_BC2=0
3932  bc_tmp%NUMBER_OF_BC3=0
3933  bc_tmp%NUMBER_OF_BC4=0
3934  bc_tmp%NUMBER_OF_BC5=0
3935  highestnodenumber=1
3936 
3937  DO i=1,numberoffaces
3938  READ(79,*) dummy(1:numberofnodes+1),flag
3939  DO j=2,numberofnodes+1
3940  IF(dummy(j)>highestnodenumber) highestnodenumber=dummy(j)
3941  ENDDO
3942  ENDDO
3943 
3944  CLOSE(79)
3945 
3946  ALLOCATE(tmp1(highestnodenumber))
3947  ALLOCATE(tmp2(highestnodenumber))
3948  ALLOCATE(tmp3(highestnodenumber))
3949  ALLOCATE(tmp4(highestnodenumber))
3950  ALLOCATE(tmp5(highestnodenumber))
3951 
3952  tmp1=0
3953  tmp2=0
3954  tmp3=0
3955  tmp4=0
3956  tmp5=0
3957 
3958  OPEN(unit = 79, file=namz,status='old')
3959  READ(79,*) numberoffaces, numberofnodes
3960 
3961  DO i=1,numberoffaces
3962  READ(79,*) dummy(1:numberofnodes+1),flag
3963  IF(flag==1) THEN
3964  DO j=2,numberofnodes+1
3965  IF(tmp1(dummy(j))==1) THEN
3966  !do nothing
3967  ELSE
3968  tmp1(dummy(j))=1
3969  bc_tmp%NUMBER_OF_BC1=bc_tmp%NUMBER_OF_BC1+1
3970  ENDIF
3971  ENDDO
3972  ELSE IF(flag==2) THEN
3973  DO j=2,numberofnodes+1
3974  IF(tmp2(dummy(j))==1) THEN
3975  !do nothing
3976  ELSE
3977  tmp2(dummy(j))=1
3978  bc_tmp%NUMBER_OF_BC2=bc_tmp%NUMBER_OF_BC2+1
3979  ENDIF
3980  ENDDO
3981  ELSE IF(flag==3) THEN
3982  DO j=2,numberofnodes+1
3983  IF(tmp3(dummy(j))==1) THEN
3984  !do nothing
3985  ELSE
3986  tmp3(dummy(j))=1
3987  bc_tmp%NUMBER_OF_BC3=bc_tmp%NUMBER_OF_BC3+1
3988  ENDIF
3989  ENDDO
3990  ELSE IF(flag==4) THEN
3991  DO j=2,numberofnodes+1
3992  IF(tmp4(dummy(j))==1) THEN
3993  !do nothing
3994  ELSE
3995  tmp4(dummy(j))=1
3996  bc_tmp%NUMBER_OF_BC4=bc_tmp%NUMBER_OF_BC4+1
3997  ENDIF
3998  ENDDO
3999  ELSEIF(flag==5) THEN
4000  DO j=2,numberofnodes+1
4001  IF(tmp5(dummy(j))==1) THEN
4002  !do nothing
4003  ELSE
4004  tmp5(dummy(j))=1
4005  bc_tmp%NUMBER_OF_BC5=bc_tmp%NUMBER_OF_BC5+1
4006  ENDIF
4007  ENDDO
4008  ENDIF
4009  ENDDO
4010 
4011  IF(bc_tmp%NUMBER_OF_BC1>0) ALLOCATE(bc_tmp%BC_NODE_NUMBER1(bc_tmp%NUMBER_OF_BC1))
4012  IF(bc_tmp%NUMBER_OF_BC2>0) ALLOCATE(bc_tmp%BC_NODE_NUMBER2(bc_tmp%NUMBER_OF_BC2))
4013  IF(bc_tmp%NUMBER_OF_BC3>0) ALLOCATE(bc_tmp%BC_NODE_NUMBER3(bc_tmp%NUMBER_OF_BC3))
4014  IF(bc_tmp%NUMBER_OF_BC4>0) ALLOCATE(bc_tmp%BC_NODE_NUMBER4(bc_tmp%NUMBER_OF_BC4))
4015  IF(bc_tmp%NUMBER_OF_BC5>0) ALLOCATE(bc_tmp%BC_NODE_NUMBER5(bc_tmp%NUMBER_OF_BC5))
4016 
4017  i1=0
4018  i2=0
4019  i3=0
4020  i4=0
4021  i5=0
4022 
4023  test=.true.
4024 
4025  DO i=1,highestnodenumber
4026  IF(tmp1(i)==1) THEN
4027  i1=i1+1
4028  bc_tmp%BC_NODE_NUMBER1(i1)=i
4029  bc_tmp%NUMBER_OF_BC=bc_tmp%NUMBER_OF_BC+1
4030  test=.false.
4031  ENDIF
4032  IF(tmp2(i)==1) THEN
4033  i2=i2+1
4034  bc_tmp%BC_NODE_NUMBER2(i2)=i
4035  bc_tmp%NUMBER_OF_BC=bc_tmp%NUMBER_OF_BC+1
4036  ENDIF
4037  IF(tmp3(i)==1) THEN
4038  i3=i3+1
4039  bc_tmp%BC_NODE_NUMBER3(i3)=i
4040  bc_tmp%NUMBER_OF_BC=bc_tmp%NUMBER_OF_BC+1
4041  ENDIF
4042  IF(tmp4(i)==1) THEN
4043  i4=i4+1
4044  bc_tmp%BC_NODE_NUMBER4(i4)=i
4045  bc_tmp%NUMBER_OF_BC=bc_tmp%NUMBER_OF_BC+1
4046  ENDIF
4047  IF(tmp5(i)==1) THEN
4048  i5=i5+1
4049  bc_tmp%BC_NODE_NUMBER5(i5)=i
4050  bc_tmp%NUMBER_OF_BC=bc_tmp%NUMBER_OF_BC+1
4051  ENDIF
4052  ENDDO
4053 
4054  CLOSE(79)
4055 
4056  WRITE(*,*)' '
4057  WRITE(*,*)'Import bc successfully... BC'
4058  WRITE(*,*)' '
4059 
4060  DEALLOCATE(dummy)
4061  DEALLOCATE(tmp1)
4062  DEALLOCATE(tmp2)
4063  DEALLOCATE(tmp3)
4064  DEALLOCATE(tmp4)
4065  DEALLOCATE(tmp5)
4066 
4067 
4068  RETURN
4069 
4070  END SUBROUTINE fluid_mechanics_io_read_bc
4071 
4072 
4073  ! OK
4074  !================================================================================================================================
4075  !
4077  SUBROUTINE fluid_mechanics_io_read_boundary_conditions_iteration(SOLVER_TYPE,BOUNDARY_VALUES,BOUNDARY_NODES, &
4078  & number_of_dimensions,boundary_condition,option,iteration)
4079  INTEGER(INTG):: SOLVER_TYPE,I,NUMBER_OF_TIME_STEPS,OPTION
4080  REAL(DP), POINTER :: BOUNDARY_VALUES(:)
4081  INTEGER(INTG), POINTER :: BOUNDARY_NODES(:)
4082  INTEGER(INTG):: NUMBER_OF_DIMENSIONS,BOUNDARY_CONDITION,NUM_BC_NODES,ITERATION
4083 
4084  CHARACTER(34) :: INPUT_FILE
4085 
4086  INTEGER(INTG):: ENDI
4087 
4088  option=1
4089  ENDI=1
4090  number_of_time_steps=1
4091  number_of_dimensions=42
4092 
4093  IF(solver_type==1) THEN !LINEAR
4094  IF(boundary_condition==1)THEN !SCALAR BOUNDARY CONDITION - SET THE NUMBER OF DIMENSIONS TO BE NUMBER OF DEPENDENT FIELD COMPONENTS (this will usually be one)!
4095  !ENDI=SIZE(BOUNDARY_VALUES)/NUMBER_OF_DIMENSIONS
4096 ! IF(J<10) THEN
4097  WRITE(input_file,'("./input/BC/BC_VALUES_",I0,".dat")') iteration
4098  !WRITE (*,*) INPUT_FILE
4099 ! ELSE IF(J<100) THEN
4100 ! WRITE(INPUT_FILE,'("./input/BC/BC_VALUES_",F10,".dat")') J
4101 ! ENDIF
4102  OPEN(unit=1, file=input_file,status='unknown')
4103  READ(1,*) num_bc_nodes
4104  ALLOCATE(boundary_nodes(num_bc_nodes))
4105  ALLOCATE(boundary_values(num_bc_nodes))
4106  DO i=1,num_bc_nodes
4107  READ(1,*) boundary_nodes(i)
4108  ENDDO
4109  DO i=1,num_bc_nodes
4110  READ(1,*) boundary_values(i)
4111  ENDDO
4112  boundary_values=boundary_values
4113  boundary_nodes=boundary_nodes
4114  !WRITE(*,*)'1! BOUNDARY_VALUES=BOUNDARY_VALUES'
4115  CLOSE(1)
4116  END IF
4117  ENDIF
4118 
4119  RETURN
4120 
4122  !
4123  !================================================================================================================================
4124  !
4126  SUBROUTINE fluid_mechanics_io_read_boundary_conditions(SOLVER_TYPE,BOUNDARY_VALUES,NUMBER_OF_DIMENSIONS,BOUNDARY_CONDITION, &
4127  & option,time_step,time,length_scale)
4129  INTEGER(INTG):: SOLVER_TYPE,I,OPTION, TIME_STEP
4130  REAL(DP), POINTER :: BOUNDARY_VALUES(:)
4131  INTEGER(INTG):: NUMBER_OF_DIMENSIONS,BOUNDARY_CONDITION,NUM_BC_NODES
4132  REAL(DP):: LENGTH_SCALE,TIME
4133 
4134  CHARACTER(34) :: INPUT_FILE
4135 
4136  !We assume a cardiac cycle to be 1s resolved with 20 time-steps which
4137  !gives us a time-step size of 0.05s
4138 
4139  INTEGER(INTG):: ENDI
4140 
4141  ENDI=42
4142 
4143 
4144  IF(solver_type==1) THEN !LINEAR
4145  IF(boundary_condition==5)THEN !MOVED WALL
4146  IF(option==0) THEN
4147  !do nothing (default)
4148  ELSE IF(option==1) THEN
4149  ENDI=SIZE(BOUNDARY_VALUES)
4150  IF(time_step<10) THEN
4151  WRITE(input_file,'("./input/motion/DISPLACEMENT_0",I0,".dat")') time_step
4152  ELSE IF(time_step<100) THEN
4153  WRITE(input_file,'("./input/motion/DISPLACEMENT_",I0,".dat")') time_step
4154  ENDIF
4155  OPEN(unit=time_step,file=input_file,status='unknown')
4156  DO i=1,endi
4157  READ(time_step,*) boundary_values(i)
4158  ENDDO
4159  boundary_values=boundary_values*length_scale
4160  WRITE(*,*)'1! BOUNDARY_VALUES=BOUNDARY_VALUES'
4161  CLOSE(time_step)
4162  ELSEIF(option==2) THEN
4163  ENDI=SIZE(BOUNDARY_VALUES)/number_of_dimensions
4164  !U,X COMPONENT
4165  boundary_values(1:endi)=0.0_dp
4166  !V,Y COMPONENT
4167  boundary_values(endi+1:endi+endi)=-0.05_dp*cos(2.0_dp*pi*1.0_dp/40.0_dp*time)
4168  !W,Z COMPONENT
4169  boundary_values(endi+endi+1:endi+endi+endi)=0.0_dp
4170  ELSEIF(option==3) THEN
4171  ENDI=SIZE(BOUNDARY_VALUES)/number_of_dimensions
4172  !U,X COMPONENT
4173  boundary_values(1:endi)=-1.0_dp
4174  !V,Y COMPONENT
4175  boundary_values(endi+1:endi+endi)=-1.0_dp
4176  !W,Z COMPONENT
4177  boundary_values(endi+endi+1:endi+endi+endi)=-1.0_dp
4178  ELSEIF(option==4) THEN
4179  ENDI=SIZE(BOUNDARY_VALUES)/number_of_dimensions
4180  !U,X COMPONENT
4181  boundary_values(1:endi)=0.0_dp
4182  !V,Y COMPONENT
4183  boundary_values(endi+1:endi+endi)=0.0_dp
4184  !W,Z COMPONENT
4185  boundary_values(endi+endi+1:endi+endi+endi)=-0.1_dp*cos(2.0_dp*pi*1.0_dp/40.0_dp*time)
4186  ELSEIF(option==69) THEN
4187  ENDI=SIZE(BOUNDARY_VALUES)/number_of_dimensions
4188  !U,X COMPONENT
4189  boundary_values(1:endi)=0.0_dp
4190  !V,Y COMPONENT
4191  boundary_values(endi+1:endi+endi)=0.0_dp
4192  !W,Z COMPONENT
4193  boundary_values(endi+endi+1:endi+endi+endi)=-100.0_dp*time
4194  ELSE
4195  stop 'Error during boundary input'
4196  ENDIF
4197  ELSEIF(boundary_condition==1)THEN !SCALAR BOUNDARY CONDITION - SET THE NUMBER OF DIMENSIONS TO BE NUMBER OF DEPENDENT FIELD COMPONENTS (this will usually be one)!
4198  !ENDI=SIZE(BOUNDARY_VALUES)/NUMBER_OF_DIMENSIONS
4199 ! IF(J<10) THEN
4200  WRITE(input_file,'("./input/BC/BC_VALUES_",F10.7,".dat")') time
4201  !WRITE (*,*) INPUT_FILE
4202 ! ELSE IF(J<100) THEN
4203 ! WRITE(INPUT_FILE,'("./input/BC/BC_VALUES_",F10,".dat")') J
4204 ! ENDIF
4205  OPEN(unit=1, file=input_file,status='unknown')
4206  READ(1,*) num_bc_nodes
4207  ALLOCATE(boundary_values(num_bc_nodes))
4208  DO i=1,num_bc_nodes
4209  READ(1,*) boundary_values(i)
4210  ENDDO
4211  boundary_values=boundary_values
4212  WRITE(*,*)'1! BOUNDARY_VALUES=BOUNDARY_VALUES'
4213  CLOSE(1)
4214  END IF
4215  ELSE IF(solver_type==2) THEN !NONLINEAR
4216  IF(boundary_condition==2)THEN !FIXED INLET
4217  IF(option==69) THEN
4218  ENDI=SIZE(BOUNDARY_VALUES)/number_of_dimensions
4219  !U,X COMPONENT
4220  boundary_values(1:endi)=0.0_dp
4221  !V,Y COMPONENT
4222  boundary_values(endi+1:endi+endi)=0.0_dp
4223  !W,Z COMPONENT
4224  boundary_values(endi+endi+1:endi+endi+endi)=1000.0_dp*time
4225  ELSE
4226  stop 'Error during boundary input'
4227  ENDIF
4228  END IF
4229  ELSE IF(solver_type==3) THEN
4230  IF(boundary_condition==2)THEN !FIXED INLET
4231  !do nothing
4232  END IF
4233  ENDIF
4234 ! WRITE(*,*)'TEST_OUTPUT',BOUNDARY_VALUES(ENDI+1),PI,TIME
4235  RETURN
4236 
4238  !
4239  !================================================================================================================================
4240  !
4242 
4243  SUBROUTINE fluidmechanics_io_updateboundaryconditionupdatenodes(GeometricField,SolverType,InletNodes, &
4244  & boundaryvalues,boundarycondition,option,time,stoptime)
4246  INTEGER(INTG) :: Err
4247  TYPE(varying_string) :: Error
4248  TYPE(field_type), POINTER :: GeometricField
4249  INTEGER(INTG) :: SolverType,Option
4250  REAL(DP) :: Value
4251  REAL(DP), ALLOCATABLE, INTENT(OUT) :: BoundaryValues(:)
4252  INTEGER(INTG) :: BoundaryCondition,MaterialSpecification,arraySize,NodeNumber,I,ComponentNumber
4253  INTEGER(INTG), PARAMETER :: Blood=1
4254  INTEGER(INTG), PARAMETER :: Plate2D=2
4255  INTEGER(INTG),PARAMETER :: Plate3D=3
4256  INTEGER(Intg), ALLOCATABLE, INTENT(OUT) :: InletNodes(:)
4257  REAL(DP) :: Time,StopTime,Y
4258 
4259  materialspecification=option
4260 
4261  IF(solvertype==3) THEN
4262  SELECT CASE(materialspecification)
4263  CASE(blood)
4264  ALLOCATE(inletnodes(5))
4265  inletnodes=(/29,7,43,13,58/)
4266  ALLOCATE(boundaryvalues(5))
4267  y=0.0_dp
4268  DO i=1,SIZE(inletnodes)
4269  y=y+0.5_dp
4270  ! BoundaryValues(I)=ABS(SIN(Y*PI/3.0_DP)*SIN(PI*Time/4.0_DP)*0.75_DP) ! 0.75 cm/s in artery with 2 cm diameter
4271  ! BoundaryValues(I)=ABS(SIN(Y*PI/3.0_DP)*SIN(PI*Time/4.0_DP)*0.075_DP) ! 0.75 cm/s in artery with 2 cm diameter
4272  boundaryvalues(i)=0.1_dp/(1.0_dp+exp(-0.5_dp*time+2.0_dp)) ! m / s
4273  ENDDO
4274  CASE(plate2d,plate3d)
4275  IF(materialspecification==plate2d) THEN
4276  OPEN (unit=1, file='./Input/2D/BC/Plate2DinletBC.txt', &
4277  & status='old', action='read')
4278  componentnumber=2
4279  ELSE
4280  OPEN (unit=1, file='./Input/3D/BC/Plate3DinletBC.txt', &
4281  & status='old', action='read')
4282  componentnumber=3
4283  ENDIF
4284  READ(1,*) arraysize
4285  IF(arraysize==0) stop "Number of boundary conditions for fluid domain is zero. Invalid."
4286  ALLOCATE(inletnodes(arraysize))
4287  ALLOCATE(boundaryvalues(arraysize))
4288  READ(1,*) inletnodes(:)
4289  arraysize=0
4290  CLOSE(1)
4291  boundaryvalues=0.0_dp
4292  DO i=1,SIZE(inletnodes)
4293  nodenumber=inletnodes(i)
4294  CALL field_parameter_set_get_node(geometricfield,field_u_variable_type, &
4295  & field_values_set_type,1,1,nodenumber,componentnumber,Value,err,error,*999)
4296  ! BoundaryValues(I)=0.1_DP*Time/StopTime! m / s
4297  boundaryvalues(i)=0.01_dp/(1.0_dp+exp(-0.5_dp*time*1.0_dp+2.0_dp)) ! m / s
4298  ! BoundaryValues(I)=0.1_DP/(1.0_DP+exp(-0.5_DP*Time*0.01_DP+5.0_DP)) ! m / s
4299  ! BoundaryValues(I)=0.005_DP/(1.0_DP+exp(-0.5_DP*Time*0.5+2.0_DP)) ! m / s
4300  ENDDO
4301  CASE DEFAULT
4302  stop "Invalid input option for position dependent boundary conditions in FSI problem."
4303  END SELECT
4304  ELSE
4305  stop "Not implemented."
4306  ENDIF
4307 
4308  RETURN
4309 999 RETURN
4311 
4312  !
4313  !================================================================================================================================
4314  !
4315 
4317  SUBROUTINE fluid_mechanics_io_read_data(SOLVER_TYPE,INPUT_VALUES,NUMBER_OF_DIMENSIONS,INPUT_TYPE, &
4318  & input_option,time_step,length_scale)
4320  INTEGER(INTG):: SOLVER_TYPE,I,INPUT_OPTION,CHECK
4321  INTEGER(INTG) :: ERR
4322  INTEGER(INTG) :: NUMBER_OF_DIMENSIONS
4323  TYPE(varying_string):: ERROR
4324  REAL(DP), POINTER :: INPUT_VALUES(:)
4325  INTEGER(INTG):: INPUT_TYPE
4326  REAL(DP) :: LENGTH_SCALE
4327 
4328  INTEGER(INTG):: ENDI,TIME_STEP
4329  CHARACTER(35) :: INPUT_FILE
4330  CHARACTER(29) :: UVEL_FILE
4331 
4332  i=solver_type
4333  number_of_dimensions=42
4334  ENDI=42
4335 
4336 
4337 ! IF(SOLVER_TYPE==1) THEN !LINEAR
4338  IF(input_type==1)THEN !POISSON VECTOR SOURCE TEMPORARY
4339  ENDI=SIZE(INPUT_VALUES)
4340 
4341 ! WRITE(*,*) "TIME_STEP", TIME_STEP
4342 
4343  IF(input_option==1) THEN
4344  IF(time_step<10) THEN
4345  WRITE(uvel_file,'("./input/data/VEL_DATA_0",I0,".dat")') time_step
4346  ELSE IF(time_step<100) THEN
4347  WRITE(uvel_file,'("./input/data/VEL_DATA_",I0,".dat")') time_step
4348  ENDIF
4349  ELSE IF(input_option==2) THEN
4350  IF(time_step<=10) THEN
4351  WRITE(uvel_file,'("./input/data/VEL_DATA_0",I0,".dat")') time_step-1
4352  ELSE IF(time_step<100) THEN
4353  WRITE(uvel_file,'("./input/data/VEL_DATA_",I0,".dat")') time_step-1
4354  ENDIF
4355  ELSE IF(input_option==3) THEN
4356  IF(time_step<10) THEN
4357  WRITE(uvel_file,'("./input/data/ORI_DATA_0",I0,".dat")') time_step
4358  ELSE IF(time_step<100) THEN
4359  WRITE(uvel_file,'("./input/data/ORI_DATA_",I0,".dat")') time_step
4360  ENDIF
4361  ELSE IF(input_option==4) THEN
4362  IF(time_step<10) THEN
4363  WRITE(uvel_file,'("./input/data/U_DATA_0",I0,".dat")') time_step
4364  ELSE IF(time_step<100) THEN
4365  WRITE(uvel_file,'("./input/data/U_DATA_",I0,".dat")') time_step
4366  ENDIF
4367  ELSE IF(input_option==5) THEN
4368  IF(time_step<10) THEN
4369  WRITE(uvel_file,'("./input/data/V_DATA_0",I0,".dat")') time_step
4370  ELSE IF(time_step<100) THEN
4371  WRITE(uvel_file,'("./input/data/V_DATA_",I0,".dat")') time_step
4372  ENDIF
4373  ELSE IF(input_option==6) THEN
4374  IF(time_step<10) THEN
4375  WRITE(uvel_file,'("./input/data/W_DATA_0",I0,".dat")') time_step
4376  ELSE IF(time_step<100) THEN
4377  WRITE(uvel_file,'("./input/data/W_DATA_",I0,".dat")') time_step
4378  ENDIF
4379  ENDIF
4380  CALL write_string(general_output_type,uvel_file,err,error,*999)
4381  OPEN(unit=42, file=uvel_file,status='unknown')
4382  READ(42,*) check
4383  IF(check/=endi) THEN
4384  stop 'Error during data input - probably wrong Lagrangian/Hermite input file!'
4385  ENDIF
4386  DO i=1,endi
4387  READ(42,*) input_values(i)
4388  ENDDO
4389  CLOSE(42)
4390 
4391  ELSE IF(input_type==42) THEN
4392 ! do nothing for now
4393  IF(input_option==0) THEN
4394  !do nothing (default)
4395  ELSE IF(input_option==1) THEN
4396  ENDI=SIZE(INPUT_VALUES)
4397  IF(time_step<=10) THEN
4398  WRITE(input_file,'("./input/motion/DISPLACEMENT_00",I0,".dat")') time_step
4399  ELSE IF(time_step<100) THEN
4400  WRITE(input_file,'("./input/motion/DISPLACEMENT_0",I0,".dat")') time_step
4401  ELSE IF(time_step<1000) THEN
4402  WRITE(input_file,'("./input/motion/DISPLACEMENT_",I0,".dat")') time_step
4403  ENDIF
4404  OPEN(unit=time_step, file=input_file,status='unknown')
4405  READ(time_step,*) check
4406  IF(check/=endi) THEN
4407  stop 'Error during data input - probably wrong Lagrangian/Hermite input file!'
4408  ENDIF
4409  DO i=1,endi
4410  READ(time_step,*) input_values(i)
4411  ENDDO
4412 ! ! TESTETSTEST
4413  input_values=input_values/length_scale
4414  WRITE(*,*)'1! INPUT_VALUES=INPUT_VALUES/LENGTH_SCALE'
4415  CLOSE(time_step)
4416  ELSE IF(input_option==2) THEN ! For Darcy, invoke the length scale (consistent with reading in the geometry data)
4417  ENDI=SIZE(INPUT_VALUES)
4418  IF(time_step<=10) THEN
4419  WRITE(input_file,'("./input/motion/DISPLACEMENT_00",I0,".dat")') time_step
4420  ELSE IF(time_step<100) THEN
4421  WRITE(input_file,'("./input/motion/DISPLACEMENT_0",I0,".dat")') time_step
4422  ELSE IF(time_step<1000) THEN
4423  WRITE(input_file,'("./input/motion/DISPLACEMENT_",I0,".dat")') time_step
4424  ENDIF
4425  OPEN(unit=time_step, file=input_file,status='unknown')
4426  DO i=1,endi
4427  READ(time_step,*) input_values(i)
4428  input_values(i) = length_scale * input_values(i)
4429  ENDDO
4430  CLOSE(time_step)
4431  ELSE IF(input_option==3) THEN
4432  ENDI=SIZE(INPUT_VALUES)
4433  IF(time_step<=10) THEN
4434  WRITE(input_file,'("./input/motion/DISPLACEMENT_00",I0,".dat")') time_step
4435  ELSE IF(time_step<100) THEN
4436  WRITE(input_file,'("./input/motion/DISPLACEMENT_0",I0,".dat")') time_step
4437  ELSE IF(time_step<500) THEN
4438  WRITE(input_file,'("./input/motion/DISPLACEMENT_",I0,".dat")') time_step
4439  ELSE IF(time_step<1000) THEN
4440  WRITE(input_file,'("./input/motion/DISPLACEMENT_",I0,".dat")') 1000-time_step
4441  ENDIF
4442  OPEN(unit=time_step, file=input_file,status='unknown')
4443  DO i=1,endi
4444  READ(time_step,*) input_values(i)
4445  input_values(i) = length_scale * input_values(i)
4446  ENDDO
4447  CLOSE(time_step)
4448  ELSE
4449  stop 'Error during data input'
4450  ENDIF
4451  ENDIF
4452  RETURN
4453 999 errorsexits("FLUID_MECHANICS_IO_READ_DATA",err,error)
4454  RETURN
4455  END SUBROUTINE fluid_mechanics_io_read_data
4456 
4457  ! OK
4458  !================================================================================================================================
4459  !
4460 
4464  INTEGER(INTG):: I,J
4465 ! INTEGER(INTG) :: ERR
4466 ! TYPE(VARYING_STRING):: ERROR
4467  INTEGER(INTG) :: TEMP(200)
4468 
4469  READ(42,*) namz
4470 
4471  OPEN(unit = 1, file=namz,status='old')
4472  READ(1,*) numberofelementsdefined(1:3)
4473 
4474 ! ALLOCATE AND READ MESH ELEMENT INFORMATION
4476 ! ! ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Reading Elements...",ERR,ERROR,*999)
4477  WRITE(*,*)'Reading Elements...'
4478  DO i=1,3
4480  ALLOCATE(mesh_info(i)%T(mesh_info(i)%Lt,numberofnodesperelement(i)),stat=alloc_error)
4481  DO j = 1,mesh_info(i)%Lt
4482  READ(1,*,end=30) temp(1:numberofnodesperelement(i))
4484  END DO
4485  END DO
4486  CLOSE(1)
4487  RETURN
4488 
4489  30 print *, 'FAILS'
4490  stop
4491 
4492  IF(alloc_error.NE.0) THEN
4493  stop 'Error during allocation'
4494  END IF
4495 
4496  RETURN
4497 ! 999 ERRORSEXITS("FLUID_MECHANICS_IO_READ_ELEMENTS",ERR,ERROR)
4498  RETURN
4499 
4500  END SUBROUTINE fluid_mechanics_io_read_elements
4501 
4502  ! OK
4503  !================================================================================================================================
4504  !
4505 
4509  INTEGER(INTG):: I
4510 ! INTEGER(INTG) :: ERR
4511 ! TYPE(VARYING_STRING):: ERROR
4512 
4513  DO i = 1,totalnumberofnodes
4514  WRITE(*,'("Node ",(I0,4x),1000( F6.3,2x ))')i,opencmiss_node_coord(i,1:3)
4515  END DO
4516 
4517  ! where are the element nodes stored -> 3 MATRICES
4518 
4519  WRITE(*,*)
4520  WRITE(*,*)
4521  DO i = 1,numberofelementsdefined(1)
4522  WRITE(*,'("M-Elements: ", (I0,3x), (1000(I0, 1x)) )')i, &
4524  END DO
4525  WRITE(*,*)
4526  DO i = 1,numberofelementsdefined(2)
4527  WRITE(*,'("V-Elements: ", (I0,3x), (1000(I0, 1x)) )')i, &
4529  END DO
4530  WRITE(*,*)
4531  DO i = 1,numberofelementsdefined(3)
4532  WRITE(*,'("P-Elements: ", (I0,3x), (1000(I0, 1x)) )')i, &
4534  END DO
4535  WRITE(*,*)
4536  WRITE(*,*)
4537  WRITE(*,*)
4538 
4539 ! RETURN
4540 ! 999 ERRORSEXITS("FLUID_MECHANICS_IO_PRINT_ON_SCREEN",ERR,ERROR)
4541 
4542  END SUBROUTINE fluid_mechanics_io_print_on_screen
4543 
4544  ! OK
4545  !================================================================================================================================
4546  !
4547 
4548  !================================================================================================================================
4549  !
4550  ! Below are the Darcy routines
4551  ! chrm, 20.08.09
4552  !
4553  !================================================================================================================================
4554 
4555  ! OK
4556  !================================================================================================================================
4557  !
4558 
4562  IMPLICIT NONE
4563 
4564  CHARACTER*90 DARCY_PARAM_FILE
4565 
4566  darcy_param_file='./input/Darcy_parameters.inp'
4567 
4568  WRITE(*,*)'Reading Darcy parameters.'
4569 
4570  OPEN(unit=37,file=darcy_param_file,status='old',action='read') ! Read base file for initial parameters
4571 
4572 
4573  DO WHILE (0 < 1)
4574  READ(37,*,end=50) in_char
4575  IF (index(in_char,'TESTCASE:') == 1) READ(37,*) darcy%TESTCASE
4576 
4577  IF (index(in_char,'STAB:') == 1) READ(37,*) darcy%STAB
4578  IF (index(in_char,'ANALYTIC:') == 1) READ(37,*) darcy%ANALYTIC
4579  IF (index(in_char,'DEBUG:') == 1) READ(37,*) darcy%DEBUG_MODE
4580 
4581  IF (index(in_char,'LENGTH:') == 1) READ(37,*) darcy%LENGTH
4582 
4583  IF (index(in_char,'GEOM_TOL:') == 1) READ(37,*) darcy%GEOM_TOL
4584  IF (index(in_char,'X1:') == 1) READ(37,*) darcy%X1
4585  IF (index(in_char,'X2:') == 1) READ(37,*) darcy%X2
4586  IF (index(in_char,'Y1:') == 1) READ(37,*) darcy%Y1
4587  IF (index(in_char,'Y2:') == 1) READ(37,*) darcy%Y2
4588  IF (index(in_char,'Z1:') == 1) READ(37,*) darcy%Z1
4589  IF (index(in_char,'Z2:') == 1) READ(37,*) darcy%Z2
4590 
4591  IF (index(in_char,'PERM:') == 1) READ(37,*) darcy%PERM
4592  IF (index(in_char,'VIS:') == 1) READ(37,*) darcy%VIS
4593  IF (index(in_char,'P_SINK:') == 1) READ(37,*) darcy%P_SINK
4594 
4595  IF (index(in_char,'BC_NUMBER_OF_WALL_NODES:') == 1) READ(37,*) darcy%BC_NUMBER_OF_WALL_NODES
4596  IF (index(in_char,'NUMBER_OF_BCS:') == 1) READ(37,*) darcy%NUMBER_OF_BCS
4597  END DO
4598 
4599 50 CLOSE(37)
4600 
4601  IF( darcy%DEBUG_MODE ) THEN
4602  write(*,*)'Read Darcy parameters from the following file = ',darcy_param_file
4603  write(*,*)'Press ENTER to continue.'
4604  read(*,*)
4605  END IF
4606 
4607  IF( abs(darcy%VIS) > 1.0e-14 ) THEN
4608  darcy%PERM_OVER_VIS =