Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SourceIMP/utilities/RocfracPrep/2Dmesh.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 SUBROUTINE mesh2d(nprocs,iProcs,ichr4)
54 
55  use meshdata
56  USE linked_list2
57 
58  IMPLICIT NONE
59 
60  include 'roccomf90.h'
61 
62  CHARACTER*4 :: ichr4
63  INTEGER :: new
64  INTEGER :: i,k
65  INTEGER :: icounter
66  INTEGER :: iprocs
67  INTEGER :: nprocs
68  INTEGER :: numnpnew
69  INTEGER :: glbndnum
70  LOGICAL :: surfaceelonproc
71  INTEGER :: testnumsurfel
72  INTEGER, POINTER, DIMENSION(:,:) :: nodeflag2d
73  INTEGER :: itestcnt
74 
75 ! User-defined list element
76 ! The Link_Type field MUST be the FIRST in the user-defined list element
77 ! Note pointer to data so as to easily create sublists
78 
79  TYPE user_type_surfndlist
80  TYPE(link_type) :: link
81  TYPE(user_data_type_procnodelist), POINTER :: data
82  END TYPE user_type_surfndlist
83 
84  TYPE user_data_type_surfndlist
85  INTEGER :: glbndnum
86  END TYPE user_data_type_surfndlist
87 
88 ! Auxilliary data type required for the transfer function
89  TYPE user_ptr_type_surfndlist
90  TYPE(user_type_surfndlist), POINTER :: p
91  END TYPE user_ptr_type_surfndlist
92 
93  TYPE(user_ptr_type_surfndlist) :: user_surfndlist
94  TYPE(list_type) :: surfndlist
95 
96  TYPE(link_ptr_type) :: link
97 
98  TYPE(surfmesh_tri3_ptr),POINTER :: ptrtri3
99  TYPE(surfmesh_tri6_ptr), POINTER :: ptrtri6
100  TYPE(surfmesh_hex8_ptr), POINTER :: ptrhex8
101 
102  INTEGER, POINTER, DIMENSION(:) :: bcsurfflagzero,bcsurfflagone,bcsurfflagtwo
103 
104 
105  CHARACTER(*), PARAMETER :: surwin = "sfrac"
106 
107  REAL*8, POINTER, DIMENSION(:,:) :: meshcoor
108 
109  INTEGER, POINTER, DIMENSION(:) :: elflag_list
110 
111  INTEGER :: write_attr, set_option, ierrflg, sur_all
112  INTEGER :: inb, ibu, ini
113 
114 ! obtain function handle ------------------------------------------------------
115 
116  write_attr = com_get_function_handle( 'OUT.write_attribute')
117  set_option = com_get_function_handle( 'OUT.set_option')
118 
119 ! do not append process rank -----------------
120 
121  CALL com_call_function( set_option, 2, 'rankwidth', '0')
122 ! CALL COM_call_function( set_option, 2, 'mode', 'w')
123 
124  CALL com_new_window( surwin )
125 
126 ! --------------------------------------------------
127 ! ------------------------------------------------
128 ! (1) Ignitable interface (BURNING)
129 ! ---------------------------------------
130 ! -------------------------------------
131 
132  inb = iprocs*100 + 1
133  ibu = iprocs*100 + 2
134  ini = iprocs*100 + 3
135 
136  numnpnew = 0
137 
138  IF(meshtype2d.EQ.6)THEN
139 
140 ! renumber elements nodes
141 
142  ptrtri6 => surfmesh_tri6_sf_head
143  DO WHILE(ASSOCIATED(ptrtri6))
144 ! mark that changed surface mesh by changing it to a negative
145  DO i = 1, 6
146  glbndnum = ptrtri6%ElemData(i)
147  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
148  numnpnew = numnpnew + 1
149  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
150  ENDIF
151  ENDDO
152  ptrtri6 => ptrtri6%next
153  ENDDO
154 
155 ! check for special situation where there is just an edge on the surface
156 
157  ptrtri6 => surfmesh_tri6_sf_head
158  DO WHILE(ASSOCIATED(ptrtri6))
159  DO i = 1, 6
160  glbndnum = ptrtri6%ElemData(i)
161  DO k = 1,maxnumberofprocstosharenode
162  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
163  print*,'Detected surface sliver for pressure'
164  numnpnew = numnpnew + 1
165  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
166  ENDIF
167  ENDDO
168  ENDDO
169  ptrtri6 => ptrtri6%next
170  ENDDO
171 
172 ! write the surface mesh files
173 
174  IF(numeltet2d(1,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
175  print*,'ERROR, found in surface mesh for Solid/Fluid Ignitable interface'
176  print*,'Found ',numeltet2d(1,iprocs),'triangles but no nodes'
177  print*,'Stopping'
178  stop
179  ENDIF
180 
181  IF(numeltet2d(1,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
182  print*,'Warning: found in sufrace mesh for Solid/Fluid Ignitable interface'
183  print*,'Found ',numnpnew,' nodes but no triangles'
184  ENDIF
185 
186 
187 ! IF(NumNPNew.NE.0.AND.NumEltet2D(1,iProcs).NE.0) &
188 ! write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumEltet2D(1,iProcs)*4,'F=FEPOINT, ET=TRIANGLE'
189 
190  nodeflag(:)= abs( nodeflag(:)) ! reset
191 
192  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
193 
194  numnpnew = 0
195 ! renumber elements nodes
196 
197  ptrtri6 => surfmesh_tri6_sf_head
198  DO WHILE(ASSOCIATED(ptrtri6))
199 ! mark that changed surface mesh by changing it to a negative
200  DO i = 1, 6
201  glbndnum = ptrtri6%ElemData(i)
202  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
203  numnpnew = numnpnew + 1
204 !!$ WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
205 !!$ WRITE(4003,*) coor(1:3,glbNdNum)
206  nodeflag2d(1,numnpnew) = glbndnum
207  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
208  nodeflag(glbndnum) = - numnpnew
209  ENDIF
210  ENDDO
211  ptrtri6 => ptrtri6%next
212  ENDDO
213  ptrtri6 => surfmesh_tri6_sf_head
214  DO WHILE(ASSOCIATED(ptrtri6))
215  DO i = 1, 6
216  glbndnum = ptrtri6%ElemData(i)
217  DO k = 1,maxnumberofprocstosharenode
218  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
219  numnpnew = numnpnew + 1
220  print*,'Dectected surface sliver for pressure, writing'
221 ! WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
222  nodeflag2d(1,numnpnew) = glbndnum
223  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
224  nodeflag(glbndnum) = - numnpnew
225  ENDIF
226  ENDDO
227  ENDDO
228  ptrtri6 => ptrtri6%next
229  ENDDO
230 
231  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
232  CALL com_set_size( surwin//'.nc', ibu, numnpnew )
233  CALL com_allocate_array(surwin//'.nc', ibu, meshcoor, 3)
234 
235  DO i = 1, numnpnew
236  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
237  END DO
238 
239 ! nodeflag to VolumeNode
240 
241  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, 'm')
242  CALL com_set_array( surwin//'.bv', ibu, nodeflag2d(2,1), 2)
243 
244 
245 ! connectiveity
246  testnumsurfel = numeltet2d(1,iprocs)
247 
248  CALL com_set_size( surwin//'.:t6', ibu, testnumsurfel)
249  CALL com_allocate_array( surwin//'.:t6', ibu, elconntable, 6)
250 
251 
252 ! elemflag to VolumeElem
253 
254  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
255  CALL com_allocate_array(surwin//'.bf2c', ibu, elflag_list, 1)
256 
257  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
258  CALL com_set_size( surwin//'.bcflag', ibu, 1)
259  CALL com_allocate_array(surwin//'.bcflag', ibu, bcsurfflagone, 1)
260 
261  bcsurfflagone(1) = 1
262 
263  testnumsurfel = 0
264  ptrtri6 => surfmesh_tri6_sf_head
265  DO WHILE(ASSOCIATED(ptrtri6))
266  surfaceelonproc = .false.
267 !!$ do i = 1, 6
268 !!$ IF(NodeFlag(ptrtri6%ElemData(i)).eq.0) THEN
269 !!$ SurfaceElOnProc = .false.
270 !!$ exit
271 !!$ endif
272 !!$ enddo
273 
274  IF(epart(ptrtri6%ElemData(7)).EQ.iprocs) surfaceelonproc = .true.
275  IF(surfaceelonproc) THEN
276 !!$ WRITE(4001,'(10i10)') ABS(NodeFlag(ptrtri6%ElemData(1:6))),&
277 !!$ ElFlag(ptrtri6%ElemData(7)),1,1 ! last two are not really used
278 
279  testnumsurfel = testnumsurfel + 1
280  elconntable(1:6,testnumsurfel) = abs(nodeflag(ptrtri6%ElemData(1:6)))
281  elflag_list(testnumsurfel) = elflag(ptrtri6%ElemData(7))
282 
283 !!$ write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(1))),ABS(NodeFlag(ptrtri6%ElemData(4))),ABS(NodeFlag(ptrtri6%ElemData(6)))
284 !!$ write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(4))),ABS(NodeFlag(ptrtri6%ElemData(2))),ABS(NodeFlag(ptrtri6%ElemData(5)))
285 !!$ write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(6))),ABS(NodeFlag(ptrtri6%ElemData(4))),ABS(NodeFlag(ptrtri6%ElemData(5)))
286 !!$ write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(3))),ABS(NodeFlag(ptrtri6%ElemData(6))),ABS(NodeFlag(ptrtri6%ElemData(5)))
287  ENDIF
288  ptrtri6 => ptrtri6%next
289  ENDDO
290 
291  IF(testnumsurfel.NE.numeltet2d(1,iprocs))THEN
292  print*,'ERROR, number of triangles from link list in for ignitable surfaces'
293  print*,' different then that in read_patran'
294  print*,testnumsurfel,numeltet2d(1,iprocs)
295  stop
296  ENDIF
297 
298  nullify(ptrtri6)
299 
300 
301  ! obtain function handle ------------------------------------------------------
302 
303 ! write_attr = COM_get_function_handle( 'OUT.write_attribute')
304 ! set_option = COM_get_function_handle( 'OUT.set_option')
305 
306 ! IF(iProcs.EQ.1)THEN
307 
308  CALL com_call_function( set_option, 2, 'mode', 'w')
309 
310 ! ELSE
311 !
312 ! CALL COM_call_function( set_option, 2, 'mode', 'a')
313 !
314 ! ENDIF
315 
316 ! do not append process rank -----------------
317 
318 ! CALL COM_call_function( set_option, 2, 'rankwidth', '0')
319 
320 ! write surface window ------------------------
321 
322 ! sur_all = Com_get_attribute_handle( surWin//'.all')
323 ! CALL COM_call_function( write_attr, 4, 'Rocin/SurfMesh_FS_Ignt.hdf', sur_all,&
324 ! "solid_surf","00.000000")
325 
326 ! CALL COM_delete_window( surWin)
327 
328 
329  ELSE IF(meshtype2d.EQ.3)THEN
330 
331 ! renumber elements nodes
332 
333  ptrtri3 => surfmesh_tri3_sf_head
334  do while(associated(ptrtri3))
335 ! mark that changed surface mesh by changing it to a negative
336  DO i = 1, 3
337  glbndnum = ptrtri3%ElemData(i)
338  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
339  numnpnew = numnpnew + 1
340  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
341  endif
342  enddo
343  ptrtri3 => ptrtri3%next
344  ENDDO
345 
346 ! check for special situation where there is just an edge on the surface
347 
348  ptrtri3 => surfmesh_tri3_sf_head
349  do while(associated(ptrtri3))
350  DO i = 1, 3
351  glbndnum = ptrtri3%ElemData(i)
352  DO k = 1,maxnumberofprocstosharenode
353  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
354  print*,'Dectected surface sliver for pressure'
355  numnpnew = numnpnew + 1
356  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
357  ENDIF
358  ENDDO
359  ENDDO
360  ptrtri3 => ptrtri3%next
361  ENDDO
362 
363 ! write the surface mesh files
364 
365  IF(numeltet2d(1,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
366  print*,'ERROR, found in surface mesh for Solid/Fluid interface'
367  print*,'Found ',numeltet2d(1,iprocs),'triangles but no nodes'
368  print*,'Stopping'
369  stop
370  ENDIF
371 
372  IF(numeltet2d(1,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
373  print*,'Warning: found in sufrace mesh for Solid/Fluid interface'
374  print*,'Found ',numnpnew,' nodes but no triangles'
375  ENDIF
376 
377  nodeflag(:)= abs( nodeflag(:)) ! reset
378 
379  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
380 
381  numnpnew = 0
382 ! renumber elements nodes
383 
384  ptrtri3 => surfmesh_tri3_sf_head
385  do while(associated(ptrtri3))
386 ! mark that changed surface mesh by changing it to a negative
387  DO i = 1, 3
388  glbndnum = ptrtri3%ElemData(i)
389  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
390  numnpnew = numnpnew + 1
391 !!$ WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
392 !!$ write(4003,*) coor(1:3,glbNdNum)
393  nodeflag2d(1,numnpnew) = glbndnum
394  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
395  nodeflag(glbndnum) = - numnpnew
396  endif
397  enddo
398  ptrtri3 => ptrtri3%next
399  ENDDO
400  ptrtri3 => surfmesh_tri3_sf_head
401  do while(associated(ptrtri3))
402  DO i = 1, 3
403  glbndnum = ptrtri3%ElemData(i)
404  DO k = 1,maxnumberofprocstosharenode
405  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
406  numnpnew = numnpnew + 1
407  print*,'Dectected surface sliver for pressure, writing'
408 ! WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
409  nodeflag2d(1,numnpnew) = glbndnum
410  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
411  nodeflag(glbndnum) = - numnpnew
412  ENDIF
413  ENDDO
414  ENDDO
415  ptrtri3 => ptrtri3%next
416  ENDDO
417 
418  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
419  CALL com_set_size( surwin//'.nc', ibu, numnpnew )
420  CALL com_allocate_array(surwin//'.nc', ibu, meshcoor, 3)
421 
422  DO i = 1, numnpnew
423  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
424  END DO
425 
426 ! nodeflag to VolumeNode
427 
428  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, 'm')
429  CALL com_set_array( surwin//'.bv', ibu, nodeflag2d(2,1), 2)
430 
431 
432 ! connectiveity
433  testnumsurfel = numeltet2d(1,iprocs) ! 0
434 
435  CALL com_set_size( surwin//'.:t3', ibu, testnumsurfel)
436  CALL com_allocate_array( surwin//'.:t3', ibu, elconntable, 3)
437 
438 
439 ! elemflag to VolumeElem
440 
441  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
442  CALL com_allocate_array(surwin//'.bf2c', ibu, elflag_list, 1)
443 
444  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
445  CALL com_set_size( surwin//'.bcflag', ibu, 1)
446  CALL com_allocate_array(surwin//'.bcflag', ibu, bcsurfflagone, 1)
447 
448  bcsurfflagone(1) = 1
449 
450  testnumsurfel = 0
451  ptrtri3 => surfmesh_tri3_sf_head
452  do while(associated(ptrtri3))
453  surfaceelonproc = .false.
454 !!$ do i = 1, 6
455 !!$ IF(NodeFlag(ptrtri6%ElemData(i)).eq.0) THEN
456 !!$ SurfaceElOnProc = .false.
457 !!$ exit
458 !!$ endif
459 !!$ enddo
460 
461  IF(epart(ptrtri3%ElemData(4)).EQ.iprocs) surfaceelonproc = .true.
462  IF(surfaceelonproc) THEN
463 !!$ WRITE(4001,'(10i10)') ABS(NodeFlag(ptrtri3%ElemData(1:3))),&
464 !!$ ElFlag(ptrtri3%ElemData(4)),1,1 ! last two are not really used
465 
466  testnumsurfel = testnumsurfel + 1
467  elconntable(1:3,testnumsurfel) = abs(nodeflag(ptrtri3%ElemData(1:3)))
468  elflag_list(testnumsurfel) = elflag(ptrtri3%ElemData(4))
469 
470 ! write(4003,*) ABS(NodeFlag(ptrtri3%ElemData(1))),ABS(NodeFlag(ptrtri3%ElemData(2))),ABS(NodeFlag(ptrtri3%ElemData(3)))
471  ENDIF
472  ptrtri3 => ptrtri3%next
473  ENDDO
474 
475  IF(testnumsurfel.NE.numeltet2d(1,iprocs))THEN
476  print*,'ERROR, number of triangles from link list in Ignitable surface mesh'
477  print*,' different then that in read_patran'
478  print*,' For Fluid Solid Inteface mesh'
479  print*,testnumsurfel,numeltet2d(1,iprocs)
480  stop
481  ENDIF
482 
483  nullify(ptrtri3)
484  ! obtain function handle ------------------------------------------------------
485 
486 ! write_attr = COM_get_function_handle( 'OUT.write_attribute')
487 ! set_option = COM_get_function_handle( 'OUT.set_option')
488 
489 ! IF(iProcs.EQ.1)THEN
490 
491  CALL com_call_function( set_option, 2, 'mode', 'w')
492 
493 ! ELSE
494 !
495 ! CALL COM_call_function( set_option, 2, 'mode', 'a')
496 !
497 ! ENDIF
498 
499 ! do not append process rank -----------------
500 
501 ! CALL COM_call_function( set_option, 2, 'rankwidth', '0')
502 
503 ! write surface window ------------------------
504 
505 ! sur_all = Com_get_attribute_handle( surWin//'.all')
506 ! CALL COM_call_function( write_attr, 4, 'Rocin/SurfMesh_FS_Ignt.hdf', sur_all,&
507 ! "solid_surf","00.000000")
508 
509 ! CALL COM_delete_window( surWin)
510 
511 
512 
513  ELSE IF(meshtype2d.EQ.4)THEN ! quad surface mesh ! 8 node brick
514 
515 
516 ! renumber elements nodes
517 
518  ptrhex8 => surfmesh_hex8_sf_head
519  DO WHILE(ASSOCIATED(ptrhex8))
520 ! mark that changed surface mesh by changing it to a negative
521  DO i = 1, 4
522  glbndnum = ptrhex8%ElemData(i)
523  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
524  numnpnew = numnpnew + 1
525  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
526  ENDIF
527  ENDDO
528  ptrhex8 => ptrhex8%next
529  ENDDO
530 
531 ! write the surface mesh files
532 
533  IF(numelhex2d(1,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
534  print*,'ERROR, found in surface mesh for Solid/Fluid Ignitable interface'
535  print*,'Found ',numelhex2d(1,iprocs),'quads but no nodes'
536  print*,'Stopping'
537  stop
538  ENDIF
539 
540  IF(numelhex2d(1,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
541  print*,'Warning: found in sufrace mesh for Solid/Fluid Ignitable interface'
542  print*,'Found ',numnpnew,' nodes but no quads'
543  ENDIF
544 
545 
546 ! write the surface mesh files
547 
548  ! WRITE(4001,'(2i10)') NumNpNew, NumElhex2D(1,iProcs)
549 
550 !!$ IF(NumNPNew.NE.0.AND.NumElhex2D(1,iProcs).NE.0) &
551 !!$ write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumElhex2D(1,iProcs),'F=FEPOINT, ET=QUADRILATERAL'
552 
553  nodeflag(:)= abs( nodeflag(:)) ! reset
554 
555  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
556 
557  numnpnew = 0
558 ! renumber elements nodes
559 
560  testnumsurfel = 0
561 
562  ptrhex8 => surfmesh_hex8_sf_head
563  do while(associated(ptrhex8))
564 ! mark that changed surface mesh by changing it to a negative
565  DO i = 1, 4
566  glbndnum = ptrhex8%ElemData(i)
567  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
568  numnpnew = numnpnew + 1
569 ! WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
570 ! write(4003,*) coor(1:3,glbNdNum)
571  nodeflag2d(1,numnpnew) = glbndnum
572  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
573  nodeflag(glbndnum) = - numnpnew
574  endif
575  enddo
576  ptrhex8 => ptrhex8%next
577  ENDDO
578 
579  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
580  CALL com_set_size( surwin//'.nc', ibu, numnpnew )
581  CALL com_allocate_array(surwin//'.nc', ibu, meshcoor, 3)
582 
583  DO i = 1, numnpnew
584  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
585  END DO
586 
587 ! nodeflag to VolumeNode
588 
589  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, 'm')
590  CALL com_set_array( surwin//'.bv', ibu, nodeflag2d(2,1), 2)
591 
592 ! connectiveity
593 
594  testnumsurfel = numelhex2d(1,iprocs)
595 
596  CALL com_set_size( surwin//'.:q4', ibu, testnumsurfel)
597  CALL com_allocate_array( surwin//'.:q4', ibu, elconntable, 4)
598 
599 
600 ! elemflag to VolumeElem
601 
602  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
603  CALL com_allocate_array(surwin//'.bf2c', ibu, elflag_list, 1)
604 
605  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
606  CALL com_set_size( surwin//'.bcflag', ibu, 1)
607  CALL com_allocate_array(surwin//'.bcflag', ibu, bcsurfflagone, 1)
608 
609  bcsurfflagone(1) = 1
610 
611  testnumsurfel = 0
612  ptrhex8 => surfmesh_hex8_sf_head
613  DO WHILE(ASSOCIATED(ptrhex8))
614  surfaceelonproc = .false.
615 !!$ do i = 1, 4
616 !!$ IF(NodeFlag(ptrhex8%ElemData(i)).eq.0) THEN
617 !!$ SurfaceElOnProc = .false.
618 !!$ exit
619 !!$ endif
620 !!$ enddo
621  IF(epart(ptrhex8%ElemData(5)).EQ.iprocs) surfaceelonproc = .true.
622  IF(surfaceelonproc)THEN
623 ! WRITE(4001,'(10i10)') ABS(NodeFlag(ptrhex8%ElemData(1:4))),&
624 ! ElFlag(ptrhex8%ElemData(5)),1,1 ! last two are not really used
625 
626  testnumsurfel = testnumsurfel + 1
627  elconntable(1:4,testnumsurfel) = abs(nodeflag(ptrhex8%ElemData(1:4)))
628  elflag_list(testnumsurfel) = elflag(ptrhex8%ElemData(5))
629 ! write(4003,*) ABS(NodeFlag(ptrhex8%ElemData(1:4)))
630  ENDIF
631 
632  ptrhex8 => ptrhex8%next
633  ENDDO
634 
635  IF(testnumsurfel.NE.numelhex2d(1,iprocs))THEN
636  print*,'ERROR, number of quads from link list in Ignitable interface'
637  print*,' different then that in read_patran'
638  print*,testnumsurfel,numelhex2d(1,iprocs)
639  stop
640  ENDIF
641 
642  nullify(ptrhex8)
643 
644  CALL com_call_function( set_option, 2, 'mode', 'w')
645 
646  ENDIF
647  CALL com_call_function( set_option, 2, 'mode', 'w') ! why would this be called twice ? fix
648 
649  sur_all = com_get_attribute_handle( surwin//'.all')
650 
651  CALL com_call_function( write_attr, 4, 'Rocin/SurfMesh.'//ichr4, sur_all,&
652  "isolid","00.000000")
653 
654 
655  CALL com_delete_window( surwin)
656 
657 ! --------------------------------------------------
658 ! ------------------------------------------------
659 ! (2) Non Ignitable interface (INTERACTING)
660 ! ---------------------------------------
661 ! -------------------------------------
662 
663  CALL com_new_window( surwin )
664 
665 
666 ! Reset NodeFlag for non solid/fluid interface
667 
668  IF(numnpnew.NE.0)THEN
669 
670  DO i = 1, numnpnew
671  nodeflag(nodeflag2d(1,i)) = nodeflag2d(2,i)
672  enddo
673  ENDIF
674  IF(associated(nodeflag2d)) deallocate(nodeflag2d)
675 
676 !!$ do i= 1, numnp_prmry
677 !!$ write(1001,*) i,NodeFlag(i)
678 !!$ enddo
679 
680  IF(meshtype2d.EQ.6)THEN
681 
682 ! renumber elements nodes
683 
684  ptrtri6 => surfmesh_tri6_sf_nonignt_head
685  DO WHILE(ASSOCIATED(ptrtri6))
686 ! mark that changed surface mesh by changing it to a negative
687  DO i = 1, 6
688  glbndnum = ptrtri6%ElemData(i)
689  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
690  numnpnew = numnpnew + 1
691  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
692  ENDIF
693  ENDDO
694  ptrtri6 => ptrtri6%next
695  ENDDO
696 
697 ! check for special situation where there is just an edge on the surface
698 
699  ptrtri6 => surfmesh_tri6_sf_nonignt_head
700  DO WHILE(ASSOCIATED(ptrtri6))
701  DO i = 1, 6
702  glbndnum = ptrtri6%ElemData(i)
703  DO k = 1,maxnumberofprocstosharenode
704  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
705  print*,'Detected surface sliver for pressure'
706  numnpnew = numnpnew + 1
707  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
708  ENDIF
709  ENDDO
710  ENDDO
711  ptrtri6 => ptrtri6%next
712  ENDDO
713 
714 ! write the surface mesh files
715 
716  IF(numeltet2d(3,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
717  print*,'ERROR, found in surface mesh for Solid/Fluid Ignitable interface'
718  print*,'Found ',numeltet2d(3,iprocs),'triangles but no nodes'
719  print*,'Stopping'
720  stop
721  ENDIF
722 
723  IF(numeltet2d(3,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
724  print*,'Warning: found in sufrace mesh for Solid/Fluid Ignitable interface'
725  print*,'Found ',numnpnew,' nodes but no triangles'
726  ENDIF
727 
728 
729 
730 !!$ IF(NumNPNew.NE.0.AND.NumEltet2D(3,iProcs).NE.0) &
731 !!$ write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumEltet2D(3,iProcs)*4,'F=FEPOINT, ET=TRIANGLE'
732 
733  nodeflag(:)= abs( nodeflag(:)) ! reset
734 
735  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
736 
737  numnpnew = 0
738 ! renumber elements nodes
739 
740  ptrtri6 => surfmesh_tri6_sf_nonignt_head
741  DO WHILE(ASSOCIATED(ptrtri6))
742 ! mark that changed surface mesh by changing it to a negative
743  DO i = 1, 6
744  glbndnum = ptrtri6%ElemData(i)
745  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
746  numnpnew = numnpnew + 1
747  !WRITE(4001,'(3(1X,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
748  !WRITE(4003,*) coor(1:3,glbNdNum)
749  nodeflag2d(1,numnpnew) = glbndnum
750  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
751  nodeflag(glbndnum) = - numnpnew
752  ENDIF
753  ENDDO
754  ptrtri6 => ptrtri6%next
755  ENDDO
756  ptrtri6 => surfmesh_tri6_sf_nonignt_head
757  DO WHILE(ASSOCIATED(ptrtri6))
758  DO i = 1, 6
759  glbndnum = ptrtri6%ElemData(i)
760  DO k = 1,maxnumberofprocstosharenode
761  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
762  numnpnew = numnpnew + 1
763  print*,'Dectected surface sliver for pressure, writing'
764  !WRITE(4001,'(3(1X,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
765  nodeflag2d(1,numnpnew) = glbndnum
766  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
767  nodeflag(glbndnum) = - numnpnew
768  ENDIF
769  ENDDO
770  ENDDO
771  ptrtri6 => ptrtri6%next
772  ENDDO
773 
774  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
775  CALL com_set_size( surwin//'.nc', inb, numnpnew )
776  CALL com_allocate_array(surwin//'.nc', inb, meshcoor, 3)
777 
778  DO i = 1, numnpnew
779  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
780  END DO
781 
782 
783  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, '')
784 ! CALL COM_set_size( surWin//'.bv', iNB, NumNpNew )
785  CALL com_set_array(surwin//'.bv', inb, nodeflag2d(2,1), 2)
786 
787 ! connectiveity
788 
789 
790  testnumsurfel = numeltet2d(3,iprocs)
791 
792  CALL com_set_size( surwin//'.:t6', inb, testnumsurfel)
793  CALL com_allocate_array( surwin//'.:t6', inb, elconntable, 6)
794 
795 
796  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
797  CALL com_allocate_array(surwin//'.bf2c', inb, elflag_list, 1)
798 
799  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
800  CALL com_set_size( surwin//'.bcflag', inb, 1)
801  CALL com_allocate_array(surwin//'.bcflag', inb, bcsurfflagzero, 1)
802 
803  bcsurfflagzero(1) = 0
804 
805 
806 
807  testnumsurfel = 0
808  ptrtri6 => surfmesh_tri6_sf_nonignt_head
809  DO WHILE(ASSOCIATED(ptrtri6))
810  surfaceelonproc = .false.
811 
812  IF(epart(ptrtri6%ElemData(7)).EQ.iprocs) surfaceelonproc = .true.
813  IF(surfaceelonproc) THEN
814 !!$ WRITE(4001,'(10i10)') ABS(NodeFlag(ptrtri6%ElemData(1:6))),&
815 !!$ ElFlag(ptrtri6%ElemData(7)),1,1 ! last two are not really used
816 
817  testnumsurfel = testnumsurfel + 1
818  elconntable(1:6,testnumsurfel) = abs(nodeflag(ptrtri6%ElemData(1:6)))
819  elflag_list(testnumsurfel) = elflag(ptrtri6%ElemData(7))
820  !write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(1))),ABS(NodeFlag(ptrtri6%ElemData(4))),ABS(NodeFlag(ptrtri6%ElemData(6)))
821  !write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(4))),ABS(NodeFlag(ptrtri6%ElemData(2))),ABS(NodeFlag(ptrtri6%ElemData(5)))
822  !write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(6))),ABS(NodeFlag(ptrtri6%ElemData(4))),ABS(NodeFlag(ptrtri6%ElemData(5)))
823  !write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(3))),ABS(NodeFlag(ptrtri6%ElemData(6))),ABS(NodeFlag(ptrtri6%ElemData(5)))
824  ENDIF
825  ptrtri6 => ptrtri6%next
826  ENDDO
827 
828  IF(testnumsurfel.NE.numeltet2d(3,iprocs))THEN
829  print*,'ERROR, number of triangles from link list in Mesh2d'
830  print*,' different then that in read_patran'
831  print*,testnumsurfel,numeltet2d(3,iprocs)
832  stop
833  ENDIF
834 
835  nullify(ptrtri6)
836 
837 
838  ! obtain function handle ------------------------------------------------------
839 
840 ! write_attr = COM_get_function_handle( 'OUT.write_attribute')
841 ! set_option = COM_get_function_handle( 'OUT.set_option')
842 
843 ! IF(iProcs.EQ.1)THEN
844 
845 ! CALL COM_call_function( set_option, 2, 'mode', 'w')
846 
847 ! ELSE
848 
849 ! CALL COM_call_function( set_option, 2, 'mode', 'a')
850 
851 ! ENDIF
852 
853 ! do not append process rank -----------------
854 
855 ! CALL COM_call_function( set_option, 2, 'rankwidth', '0')
856 
857 ! write surface window ------------------------
858 
859 ! sur_all = Com_get_attribute_handle( surWin//'.all')
860 
861 ! CALL COM_call_function( write_attr, 4, 'Rocin/SurfMesh_FS_NonIgnt.hdf', sur_all,&
862 ! "solid_surf","00.000000")
863 
864 ! CALL COM_delete_window( surWin)
865 
866 
867  ELSE IF(meshtype2d.EQ.3)THEN
868 
869 ! renumber elements nodes
870 
871  ptrtri3 => surfmesh_tri3_sf_nonignt_head
872  do while(associated(ptrtri3))
873 ! mark that changed surface mesh by changing it to a negative
874  DO i = 1, 3
875  glbndnum = ptrtri3%ElemData(i)
876  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
877  numnpnew = numnpnew + 1
878  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
879  endif
880  enddo
881  ptrtri3 => ptrtri3%next
882  ENDDO
883 
884 ! check for special situation where there is just an edge on the surface
885 
886  ptrtri3 => surfmesh_tri3_sf_nonignt_head
887  do while(associated(ptrtri3))
888  DO i = 1, 3
889  glbndnum = ptrtri3%ElemData(i)
890  DO k = 1,maxnumberofprocstosharenode
891  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
892  print*,'Dectected surface sliver for pressure'
893  numnpnew = numnpnew + 1
894  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
895 
896  ENDIF
897  ENDDO
898  ENDDO
899  ptrtri3 => ptrtri3%next
900  ENDDO
901 
902 ! write the surface mesh files
903 
904 
905  !WRITE(4001,'(2i10)') NumNpNew, NumEltet2D(3,iProcs)
906 
907  IF(numeltet2d(3,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
908  print*,'ERROR, found in surface mesh for Solid/Fluid interface'
909  print*,'Found ',numeltet2d(3,iprocs),'triangles but no nodes'
910  print*,'Stopping'
911  stop
912  ENDIF
913 
914  IF(numeltet2d(3,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
915  print*,'Warning: found in sufrace mesh for Solid/Fluid interface'
916  print*,'Found ',numnpnew,' nodes but no triangles'
917  ENDIF
918 
919 ! IF(NumNPNew.NE.0.AND.NumEltet2D(3,iProcs).NE.0) &
920 ! write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumEltet2D(3,iProcs),'F=FEPOINT, ET=TRIANGLE'
921 
922  nodeflag(:)= abs( nodeflag(:)) ! reset
923 
924  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
925 
926  numnpnew = 0
927 ! renumber elements nodes
928 
929  ptrtri3 => surfmesh_tri3_sf_nonignt_head
930  do while(associated(ptrtri3))
931 ! mark that changed surface mesh by changing it to a negative
932  DO i = 1, 3
933  glbndnum = ptrtri3%ElemData(i)
934  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
935  numnpnew = numnpnew + 1
936 ! WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
937 ! write(4003,*) coor(1:3,glbNdNum)
938  nodeflag2d(1,numnpnew) = glbndnum
939  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
940  nodeflag(glbndnum) = - numnpnew
941  endif
942  enddo
943  ptrtri3 => ptrtri3%next
944  ENDDO
945  ptrtri3 => surfmesh_tri3_sf_nonignt_head
946  do while(associated(ptrtri3))
947  DO i = 1, 3
948  glbndnum = ptrtri3%ElemData(i)
949  DO k = 1,maxnumberofprocstosharenode
950  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
951  numnpnew = numnpnew + 1
952  print*,'Dectected surface sliver for pressure, writing'
953 ! WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
954  nodeflag2d(1,numnpnew) = glbndnum
955  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
956  nodeflag(glbndnum) = - numnpnew
957  ENDIF
958  ENDDO
959  ENDDO
960  ptrtri3 => ptrtri3%next
961  ENDDO
962 
963  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
964  CALL com_set_size( surwin//'.nc', inb, numnpnew )
965  CALL com_allocate_array(surwin//'.nc', inb, meshcoor, 3)
966 
967  DO i = 1, numnpnew
968  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
969  END DO
970 
971 
972  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, '')
973 ! CALL COM_set_size( surWin//'.bv', iNB, NumNpNew )
974  CALL com_set_array(surwin//'.bv', inb, nodeflag2d(2,1), 2)
975 ! connectiveity
976 
977  testnumsurfel = numeltet2d(3,iprocs)
978 
979  CALL com_set_size( surwin//'.:t3', inb, testnumsurfel)
980  CALL com_allocate_array( surwin//'.:t3', inb, elconntable, 3)
981 
982 
983  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
984  CALL com_allocate_array(surwin//'.bf2c', inb, elflag_list, 1)
985 
986  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
987  CALL com_set_size( surwin//'.bcflag', inb, 1)
988  CALL com_allocate_array(surwin//'.bcflag', inb, bcsurfflagzero, 1)
989 
990  bcsurfflagzero(1) = 0
991 
992  testnumsurfel = 0
993  ptrtri3 => surfmesh_tri3_sf_nonignt_head
994  do while(associated(ptrtri3))
995  surfaceelonproc = .false.
996 !!$ do i = 1, 6
997 !!$ IF(NodeFlag(ptrtri6%ElemData(i)).eq.0) THEN
998 !!$ SurfaceElOnProc = .false.
999 !!$ exit
1000 !!$ endif
1001 !!$ enddo
1002 
1003  IF(epart(ptrtri3%ElemData(4)).EQ.iprocs) surfaceelonproc = .true.
1004  IF(surfaceelonproc) THEN
1005 ! WRITE(4001,'(10i10)') ABS(NodeFlag(ptrtri3%ElemData(1:3))),&
1006 ! ElFlag(ptrtri3%ElemData(4)),1,1 ! last two are not really used
1007  testnumsurfel = testnumsurfel + 1
1008  elconntable(1:3,testnumsurfel) = abs(nodeflag(ptrtri3%ElemData(1:3)))
1009  elflag_list(testnumsurfel) = elflag(ptrtri3%ElemData(4))
1010  ENDIF
1011  ptrtri3 => ptrtri3%next
1012  ENDDO
1013 
1014  IF(testnumsurfel.NE.numeltet2d(3,iprocs))THEN
1015  print*,'ERROR, number of triangles from link list in Mesh2d'
1016  print*,' different then that in read_patran'
1017  print*,' For Fluid Solid Inteface mesh'
1018  print*,testnumsurfel,numeltet2d(3,iprocs)
1019  stop
1020  ENDIF
1021 
1022  nullify(ptrtri3)
1023 
1024  ! obtain function handle ------------------------------------------------------
1025 
1026 ! write_attr = COM_get_function_handle( 'OUT.write_attribute')
1027 ! set_option = COM_get_function_handle( 'OUT.set_option')
1028 
1029 ! IF(iProcs.EQ.1)THEN
1030 
1031 ! CALL COM_call_function( set_option, 2, 'mode', 'w')
1032 
1033 ! ELSE
1034 
1035 ! CALL COM_call_function( set_option, 2, 'mode', 'a')
1036 
1037 ! ENDIF
1038 
1039 ! do not append process rank -----------------
1040 
1041 ! CALL COM_call_function( set_option, 2, 'rankwidth', '0')
1042 
1043 ! write surface window ------------------------
1044 
1045 ! sur_all = Com_get_attribute_handle( surWin//'.all')
1046 
1047 ! CALL COM_call_function( write_attr, 4, 'Rocin/SurfMesh_FS_NonIgnt.hdf', sur_all,&
1048 ! "solid_surf","00.000000")
1049 
1050 ! CALL COM_delete_window( surWin)
1051 
1052  ELSE IF(meshtype2d.EQ.4)THEN ! quad surface mesh ! 8 node brick
1053 
1054 
1055 ! renumber elements nodes
1056 
1057  ptrhex8 => surfmesh_hex8_sf_nonignt_head
1058  DO WHILE(ASSOCIATED(ptrhex8))
1059 ! mark that changed surface mesh by changing it to a negative
1060  DO i = 1, 4
1061  glbndnum = ptrhex8%ElemData(i)
1062  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
1063  numnpnew = numnpnew + 1
1064  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
1065  ENDIF
1066  ENDDO
1067  ptrhex8 => ptrhex8%next
1068  ENDDO
1069 
1070 ! write the surface mesh files
1071 
1072  IF(numelhex2d(3,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
1073  print*,'ERROR, found in surface mesh for Solid/Fluid Ignitable interface'
1074  print*,'Found ',numelhex2d(3,iprocs),'quads but no nodes'
1075  print*,'Stopping'
1076  stop
1077  ENDIF
1078 
1079  IF(numelhex2d(3,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
1080  print*,'Warning: found in sufrace mesh for Solid/Fluid Ignitable interface'
1081  print*,'Found ',numnpnew,' nodes but no quads'
1082  ENDIF
1083 
1084 ! WRITE(4001,'(2i10)') NumNpNew, NumElhex2D(3,iProcs)
1085 
1086 ! IF(NumNPNew.NE.0.AND.NumElhex2D(3,iProcs).NE.0) &
1087 ! write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumElhex2D(3,iProcs),'F=FEPOINT, ET=QUADRILATERAL'
1088 
1089  nodeflag(:)= abs( nodeflag(:)) ! reset
1090 
1091  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
1092 
1093  numnpnew = 0
1094 ! renumber elements nodes
1095  ptrhex8 => surfmesh_hex8_sf_nonignt_head
1096  do while(associated(ptrhex8))
1097 ! mark that changed surface mesh by changing it to a negative
1098  DO i = 1, 4
1099  glbndnum = ptrhex8%ElemData(i)
1100  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
1101  numnpnew = numnpnew + 1
1102 ! WRITE(4001,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
1103 ! write(4003,*) coor(1:3,glbNdNum)
1104  nodeflag2d(1,numnpnew) = glbndnum
1105  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
1106  nodeflag(glbndnum) = - numnpnew
1107  endif
1108  enddo
1109  ptrhex8 => ptrhex8%next
1110  ENDDO
1111 
1112  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
1113  CALL com_set_size( surwin//'.nc', inb, numnpnew )
1114  CALL com_allocate_array(surwin//'.nc', inb, meshcoor, 3)
1115 
1116  DO i = 1, numnpnew
1117  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
1118  END DO
1119 
1120  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, '')
1121 ! CALL COM_set_size( surWin//'.bv', iNB, NumNpNew )
1122  CALL com_set_array(surwin//'.bv', inb, nodeflag2d(2,1), 2)
1123 
1124 ! connectiveity
1125 
1126 
1127  testnumsurfel = numelhex2d(3,iprocs)
1128 
1129  CALL com_set_size( surwin//'.:q4', inb, testnumsurfel)
1130  CALL com_allocate_array( surwin//'.:q4', inb, elconntable, 4)
1131 
1132 
1133  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
1134  CALL com_allocate_array(surwin//'.bf2c', inb, elflag_list, 1)
1135 
1136  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
1137  CALL com_set_size( surwin//'.bcflag', inb, 1)
1138  CALL com_allocate_array(surwin//'.bcflag', inb, bcsurfflagzero, 1)
1139 
1140  bcsurfflagzero(1) = 0
1141 
1142 
1143  testnumsurfel = 0
1144  ptrhex8 => surfmesh_hex8_sf_nonignt_head
1145  DO WHILE(ASSOCIATED(ptrhex8))
1146  surfaceelonproc = .false.
1147 
1148  IF(epart(ptrhex8%ElemData(5)).EQ.iprocs) surfaceelonproc = .true.
1149  IF(surfaceelonproc)THEN
1150 ! WRITE(4001,'(10i10)') ABS(NodeFlag(ptrhex8%ElemData(1:4))),&
1151 ! ElFlag(ptrhex8%ElemData(5)),1,1 ! last two are not really used
1152 
1153  testnumsurfel = testnumsurfel + 1
1154  elconntable(1:4,testnumsurfel) = abs(nodeflag(ptrhex8%ElemData(1:4)))
1155  elflag_list(testnumsurfel) = elflag(ptrhex8%ElemData(5))
1156 
1157 ! write(4003,*) ABS(NodeFlag(ptrhex8%ElemData(1:4)))
1158  ENDIF
1159  ptrhex8 => ptrhex8%next
1160  ENDDO
1161 
1162  IF(testnumsurfel.NE.numelhex2d(3,iprocs))THEN
1163  print*,'ERROR, number of quads from link list in Non-Ignitable surface mesh'
1164  print*,' different then that in read_patran'
1165  print*,testnumsurfel,numelhex2d(3,iprocs)
1166  stop
1167  ENDIF
1168 
1169  nullify(ptrhex8)
1170 
1171  ENDIF
1172 
1173 ! Reset NodeFlag for non solid/fluid interface
1174 
1175  IF(numnpnew.NE.0)THEN
1176 
1177  DO i = 1, numnpnew
1178  nodeflag(nodeflag2d(1,i)) = nodeflag2d(2,i)
1179  enddo
1180  ENDIF
1181 
1182  ! write surface window ------------------------
1183 
1184  CALL com_call_function( set_option, 2, 'mode', 'a')
1185 
1186  sur_all = com_get_attribute_handle( surwin//'.all')
1187 
1188  CALL com_call_function( write_attr, 4, 'Rocin/SurfMesh.'//ichr4, sur_all,&
1189  "isolid","00.000000")
1190 
1191 
1192  CALL com_delete_window( surwin)
1193 
1194  IF(associated(nodeflag2d)) deallocate(nodeflag2d)
1195 
1196 ! --------------------------------------------------
1197 ! ------------------------------------------------
1198 ! (3) Non Fluid/Structure interface
1199 ! ---------------------------------------
1200 ! -------------------------------------
1201 
1202 ! renumber elements nodes
1203 
1204  CALL com_new_window( surwin)
1205 
1206  IF(meshtype2d.EQ.6)THEN
1207 
1208  numnpnew = 0
1209  nullify(ptrtri6)
1210 ! renumber elements nodes
1211  ptrtri6 => surfmesh_tri6_s_head
1212 
1213  DO WHILE(ASSOCIATED(ptrtri6))
1214 ! mark that changed surface mesh by changing it to a negative
1215  DO i = 1, 6
1216  glbndnum = ptrtri6%ElemData(i)
1217  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
1218  numnpnew = numnpnew + 1
1219  nodeflag(glbndnum) = -nodeflag(glbndnum)
1220  ENDIF
1221  ENDDO
1222  ptrtri6 => ptrtri6%next
1223  ENDDO
1224 
1225 
1226 ! check for special situation where there is just an edge on the surface
1227 
1228  ptrtri6 => surfmesh_tri6_s_head
1229  do while(associated(ptrtri6))
1230  DO i = 1, 6
1231  glbndnum = ptrtri6%ElemData(i)
1232  DO k = 1,maxnumberofprocstosharenode
1233  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
1234  print*,'Dectected surface sliver for pressure'
1235  numnpnew = numnpnew + 1
1236  nodeflag(glbndnum) = -nodeflag(glbndnum)
1237  ENDIF
1238  ENDDO
1239  ENDDO
1240  ptrtri6 => ptrtri6%next
1241  ENDDO
1242 
1243 ! WRITE(4002,'(2i10)') NumNpNew, NumEltet2D(2,iProcs)
1244 ! write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumEltet2D(2,iProcs),'F=FEPOINT, ET=TRIANGLE'
1245 
1246  nodeflag(:)= abs( nodeflag(:)) ! reset
1247 
1248  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
1249 
1250  numnpnew = 0
1251 ! renumber elements nodes
1252 
1253  ptrtri6 => surfmesh_tri6_s_head
1254 
1255  do while(associated(ptrtri6))
1256 ! mark that changed surface mesh by changing it to a negative
1257  DO i = 1, 6
1258  glbndnum = ptrtri6%ElemData(i)
1259  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
1260  numnpnew = numnpnew + 1
1261 ! WRITE(4002,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
1262 ! write(4003,*) coor(1:3,glbNdNum)
1263  nodeflag2d(1,numnpnew) = glbndnum
1264  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
1265  nodeflag(ptrtri6%ElemData(i)) = - numnpnew
1266  endif
1267  enddo
1268  ptrtri6 => ptrtri6%next
1269  ENDDO
1270  ptrtri6 => surfmesh_tri6_s_head
1271  do while(associated(ptrtri6))
1272  DO i = 1, 6
1273  glbndnum = ptrtri6%ElemData(i)
1274  DO k = 1,maxnumberofprocstosharenode
1275  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
1276  print*,'Dectected surface sliver for pressure'
1277  numnpnew = numnpnew + 1
1278 ! WRITE(4002,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
1279  nodeflag2d(1,numnpnew) = glbndnum
1280  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
1281  nodeflag(glbndnum) = - numnpnew
1282  ENDIF
1283  ENDDO
1284  ENDDO
1285  ptrtri6 => ptrtri6%next
1286  ENDDO
1287 
1288 
1289  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
1290  CALL com_set_size( surwin//'.nc', ini, numnpnew )
1291  CALL com_allocate_array(surwin//'.nc', ini, meshcoor, 3)
1292 
1293  DO i = 1, numnpnew
1294  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
1295  END DO
1296 
1297 
1298  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, '')
1299 ! CALL COM_set_size( surWin//'.bv', iNI, NumNpNew )
1300  CALL com_set_array( surwin//'.bv', ini, nodeflag2d(2,1), 2)
1301 
1302 
1303 ! connectiveity
1304  testnumsurfel = numeltet2d(2,iprocs)
1305  print*,'Number of Non-Interacting nodes and elements',numnpnew,testnumsurfel
1306 
1307  CALL com_set_size( surwin//'.:t6', ini, testnumsurfel)
1308  CALL com_allocate_array( surwin//'.:t6', ini, elconntable, 6)
1309 
1310 
1311  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
1312  CALL com_allocate_array(surwin//'.bf2c', ini, elflag_list, 1)
1313 
1314  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
1315  CALL com_set_size( surwin//'.bcflag', ini, 1)
1316  CALL com_allocate_array(surwin//'.bcflag', ini, bcsurfflagtwo, 1)
1317 
1318  bcsurfflagtwo(1) = 2
1319 
1320  testnumsurfel = 0
1321  ptrtri6 => surfmesh_tri6_s_head
1322  do while(associated(ptrtri6))
1323  surfaceelonproc = .false.
1324 !!$ do i = 1, 6
1325 !!$ IF(NodeFlag(ptrtri6%ElemData(i)).eq.0)THEN
1326 !!$ SurfaceElOnProc = .false.
1327 !!$ exit
1328 !!$ endif
1329 !!$ enddo
1330  IF(epart(ptrtri6%ElemData(7)).EQ.iprocs) surfaceelonproc = .true.
1331  IF(surfaceelonproc) THEN
1332  ! WRITE(4002,'(10i10)') ABS(NodeFlag(ptrtri6%ElemData(1:6))),&
1333 ! ElFlag(ptrtri6%ElemData(7)),2,1 ! last two are not really used
1334 !!$ IF(SurfaceElOnProc) THEN
1335 !!$ write(4003,*) ABS(NodeFlag(ptrtri6%ElemData(1:6)))
1336 !!$ ENDIF
1337  testnumsurfel = testnumsurfel + 1
1338  elconntable(1:6,testnumsurfel) = abs(nodeflag(ptrtri6%ElemData(1:6)))
1339  ENDIF
1340  ptrtri6 => ptrtri6%next
1341  ENDDO
1342 
1343  IF(testnumsurfel.NE.numeltet2d(2,iprocs))THEN
1344  print*,'ERROR, number of triangles from link list in Mesh2d'
1345  print*,' different then that in read_patran'
1346  print*,testnumsurfel,numeltet2d(2,iprocs)
1347  stop
1348  ENDIF
1349 
1350  nullify(ptrtri6) ! obtain function handle ------------------------------------------------------
1351 
1352 ! write_attr = COM_get_function_handle( 'OUT.write_attribute')
1353 ! set_option = COM_get_function_handle( 'OUT.set_option')
1354 
1355 ! IF(iProcs.EQ.1)THEN
1356 
1357 ! CALL COM_call_function( set_option, 2, 'mode', 'w')
1358 
1359 ! ELSE
1360 
1361 ! CALL COM_call_function( set_option, 2, 'mode', 'a')
1362 
1363 ! ENDIF
1364 
1365 ! do not append process rank -----------------
1366 
1367 ! CALL COM_call_function( set_option, 2, 'rankwidth', '0')
1368 
1369 ! write surface window ------------------------
1370 
1371 ! sur_all = Com_get_attribute_handle( surWin//'.all')
1372 ! CALL COM_call_function( write_attr, 4, 'Rocin/SurfMesh_S.hdf', sur_all,&
1373 ! "solid_surf","00.000000")
1374 
1375 
1376 
1377  ELSE IF(meshtype2d.EQ.3)THEN
1378 
1379  numnpnew = 0
1380  nullify(ptrtri3)
1381 ! renumber elements nodes
1382 
1383  ptrtri3 => surfmesh_tri3_s_head
1384 
1385  do while(associated(ptrtri3))
1386 ! mark that changed surface mesh by changing it to a negative
1387  DO i = 1, 3
1388  glbndnum = ptrtri3%ElemData(i)
1389  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
1390  numnpnew = numnpnew + 1
1391  nodeflag(glbndnum) = -nodeflag(glbndnum)
1392  endif
1393  enddo
1394  ptrtri3 => ptrtri3%next
1395  ENDDO
1396 
1397 
1398 ! check for special situation where there is just an edge on the surface
1399 
1400  ptrtri3 => surfmesh_tri3_s_head
1401  do while(associated(ptrtri3))
1402  DO i = 1, 3
1403  glbndnum = ptrtri3%ElemData(i)
1404  DO k = 1,maxnumberofprocstosharenode
1405  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
1406  print*,'Dectected surface sliver for pressure'
1407  numnpnew = numnpnew + 1
1408  nodeflag(glbndnum) = -nodeflag(glbndnum)
1409  ENDIF
1410  ENDDO
1411  ENDDO
1412  ptrtri3 => ptrtri3%next
1413  ENDDO
1414 ! WRITE(4002,'(2i10)') NumNpNew, NumEltet2D(2,iProcs)
1415 ! write(4003,*) 'ZONE N=', NumNPNew, 'E=',NumEltet2D(2,iProcs),'F=FEPOINT, ET=TRIANGLE'
1416 
1417  nodeflag(:)= abs( nodeflag(:)) ! reset
1418 
1419  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
1420 
1421  numnpnew = 0
1422 ! renumber elements nodes
1423 
1424  ptrtri3 => surfmesh_tri3_s_head
1425 
1426  do while(associated(ptrtri3))
1427 ! mark that changed surface mesh by changing it to a negative
1428  DO i = 1, 3
1429  glbndnum = ptrtri3%ElemData(i)
1430  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
1431  numnpnew = numnpnew + 1
1432 ! WRITE(4002,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
1433 ! write(4003,*) coor(1:3,glbNdNum)
1434  nodeflag2d(1,numnpnew) = glbndnum
1435  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
1436  nodeflag(ptrtri3%ElemData(i)) = - numnpnew
1437  endif
1438  enddo
1439  ptrtri3 => ptrtri3%next
1440  ENDDO
1441  ptrtri3 => surfmesh_tri3_s_head
1442  do while(associated(ptrtri3))
1443  DO i = 1, 3
1444  glbndnum = ptrtri3%ElemData(i)
1445  DO k = 1,maxnumberofprocstosharenode
1446  IF(procndlist(glbndnum,k).EQ.iprocs.AND.nodeflag(glbndnum).GT.0)THEN
1447  print*,'Dectected surface sliver for pressure'
1448  numnpnew = numnpnew + 1
1449 ! WRITE(4002,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
1450  nodeflag2d(1,numnpnew) = glbndnum
1451  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
1452  nodeflag(ptrtri3%ElemData(i)) = - numnpnew
1453  ENDIF
1454  ENDDO
1455  ENDDO
1456  ptrtri3 => ptrtri3%next
1457  ENDDO
1458 
1459 
1460  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
1461  CALL com_set_size( surwin//'.nc', ini, numnpnew )
1462  CALL com_allocate_array(surwin//'.nc', ini, meshcoor, 3)
1463 
1464  DO i = 1, numnpnew
1465  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
1466  END DO
1467 
1468  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, '')
1469 ! CALL COM_set_size( surWin//'.bv', iNI, NumNpNew )
1470  CALL com_set_array( surwin//'.bv', ini, nodeflag2d(2,1), 2)
1471 
1472 
1473 ! connectiveity
1474  testnumsurfel = numeltet2d(2,iprocs) !0
1475 
1476  CALL com_set_size( surwin//'.:t3', ini, testnumsurfel)
1477  CALL com_allocate_array( surwin//'.:t3', ini, elconntable, 3)
1478 
1479 
1480  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
1481  CALL com_allocate_array(surwin//'.bf2c', ini, elflag_list, 1)
1482 
1483 
1484  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
1485  CALL com_set_size( surwin//'.bcflag', ini,1)
1486  CALL com_allocate_array(surwin//'.bcflag', ini, bcsurfflagtwo, 1)
1487 
1488  bcsurfflagtwo(1) = 2
1489 
1490  testnumsurfel = 0
1491  ptrtri3 => surfmesh_tri3_s_head
1492  do while(associated(ptrtri3))
1493  surfaceelonproc = .false.
1494 !!$ do i = 1, 3
1495 !!$ IF(NodeFlag(ptrtri3%ElemData(i)).eq.0)THEN
1496 !!$ SurfaceElOnProc = .false.
1497 !!$ exit
1498 !!$ endif
1499 !!$ enddo
1500  IF(epart(ptrtri3%ElemData(4)).EQ.iprocs) surfaceelonproc = .true.
1501  IF(surfaceelonproc) THEN
1502  ! WRITE(4002,'(10i10)') ABS(NodeFlag(ptrtri3%ElemData(1:3))),&
1503  ! ElFlag(ptrtri3%ElemData(4)),2,1 ! last two are not really used
1504 !!$ IF(SurfaceElOnProc) THEN
1505 !!$ write(4003,*) ABS(NodeFlag(ptrtri3%ElemData(1:6)))
1506 !!$ ENDIF
1507  testnumsurfel = testnumsurfel + 1
1508  elconntable(1:3,testnumsurfel) = abs(nodeflag(ptrtri3%ElemData(1:3)))
1509  elflag_list(testnumsurfel) = elflag(ptrtri3%ElemData(4))
1510  ENDIF
1511  ptrtri3 => ptrtri3%next
1512  ENDDO
1513 
1514  IF(testnumsurfel.NE.numeltet2d(2,iprocs))THEN
1515  print*,'ERROR, number of triangles from link list in Mesh2d'
1516  print*,' different then that in read_patran'
1517  print*,testnumsurfel,numeltet2d(2,iprocs)
1518  stop
1519  ENDIF
1520 
1521 ! write_attr = COM_get_function_handle( 'OUT.write_attribute')
1522 ! set_option = COM_get_function_handle( 'OUT.set_option')
1523 
1524 ! IF(iProcs.EQ.1)THEN
1525 
1526 ! CALL COM_call_function( set_option, 2, 'mode', 'w')
1527 
1528 ! ELSE
1529 
1530 ! CALL COM_call_function( set_option, 2, 'mode', 'a')
1531 
1532 ! ENDIF
1533 
1534 ! do not append process rank -----------------
1535 
1536 ! CALL COM_call_function( set_option, 2, 'rankwidth', '0')
1537 
1538 ! write surface window ------------------------
1539 
1540 ! sur_all = Com_get_attribute_handle( surWin//'.all')
1541 ! CALL COM_call_function( write_attr, 4, 'Rocin/SurfMesh_S.hdf', sur_all,&
1542 ! "solid_surf","00.000000")
1543 
1544 
1545 
1546  ELSE IF(meshtype2d.EQ.4)THEN
1547 
1548 
1549  numnpnew = 0
1550  nullify(ptrhex8)
1551 ! renumber elements nodes
1552 
1553  ptrhex8 => surfmesh_hex8_s_head
1554 
1555  DO WHILE(ASSOCIATED(ptrhex8))
1556 ! mark that changed surface mesh by changing it to a negative
1557  DO i = 1, 4
1558  glbndnum = ptrhex8%ElemData(i)
1559  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor, neg means already renumber node
1560  numnpnew = numnpnew + 1
1561  nodeflag(glbndnum) = -nodeflag(glbndnum) ! overwrite nodeflag here but i need it for non S/F interface
1562  ENDIF
1563  ENDDO
1564  ptrhex8 => ptrhex8%next
1565  ENDDO
1566 
1567 
1568  IF(numelhex2d(2,iprocs).NE.0.AND.numnpnew.EQ.0)THEN
1569  print*,'ERROR, found in surface mesh for Solid/Fluid interface'
1570  print*,'Found ',numelhex2d(2,iprocs),'quads but no nodes'
1571  print*,'Stopping'
1572  stop
1573  ENDIF
1574 
1575  IF(numelhex2d(2,iprocs).EQ.0.AND.numnpnew.NE.0)THEN
1576  print*,'Warning: found in sufrace mesh for Solid/Fluid interface'
1577  print*,'Found ',numnpnew,' nodes but no quads'
1578  ENDIF
1579 
1580 ! write the surface mesh files
1581 
1582 ! WRITE(4002,'(2i10)') NumNpNew, NumElhex2D(2,iProcs)
1583 
1584  nodeflag(:)= abs( nodeflag(:)) ! reset
1585 
1586  ALLOCATE(nodeflag2d(1:2,1:numnpnew)) ! for storing what I'm overwriting
1587  numnpnew = 0
1588 ! renumber elements nodes
1589 
1590  ptrhex8 => surfmesh_hex8_s_head
1591  do while(associated(ptrhex8))
1592 ! mark that changed surface mesh by changing it to a negative
1593  DO i = 1, 4
1594  glbndnum = ptrhex8%ElemData(i)
1595  IF( nodeflag(glbndnum).GT.0) THEN ! if zero means not on this processor
1596  numnpnew = numnpnew + 1
1597 ! WRITE(4002,'(3(x,e16.9),1i10)') coor(1:3,glbNdNum),NodeFlag(glbNdNum)
1598  nodeflag2d(1,numnpnew) = glbndnum
1599  nodeflag2d(2,numnpnew) = nodeflag(glbndnum)
1600  nodeflag(glbndnum) = - numnpnew
1601  endif
1602  enddo
1603  ptrhex8 => ptrhex8%next
1604  ENDDO
1605 
1606  CALL com_new_attribute( surwin//'.nc', 'n', com_double, 3, 'm')
1607  CALL com_set_size( surwin//'.nc', ini, numnpnew )
1608  CALL com_allocate_array(surwin//'.nc', ini, meshcoor, 3)
1609 
1610  DO i = 1, numnpnew
1611  meshcoor(1:3,i) = coor(1:3,nodeflag2d(1,i))
1612  END DO
1613 
1614 
1615  CALL com_new_attribute( surwin//'.bv', 'n', com_integer, 1, '')
1616 ! CALL COM_set_size( surWin//'.bv', iNI, NumNpNew )
1617  CALL com_set_array( surwin//'.bv', ini, nodeflag2d(2,1), 2)
1618 
1619 
1620 ! connectivity
1621  testnumsurfel = numelhex2d(2,iprocs)
1622  print*,'Number of Non-Interacting nodes and elements',numnpnew,testnumsurfel
1623 
1624  CALL com_set_size( surwin//'.:q4', ini, testnumsurfel)
1625  CALL com_allocate_array( surwin//'.:q4', ini, elconntable, 4)
1626 
1627 
1628  CALL com_new_attribute( surwin//'.bf2c', 'e', com_integer, 1, '')
1629  CALL com_allocate_array(surwin//'.bf2c', ini, elflag_list, 1)
1630 
1631  CALL com_new_attribute( surwin//'.bcflag', 'p', com_integer, 1, '')
1632  CALL com_set_size( surwin//'.bcflag', ini, 1)
1633  CALL com_allocate_array(surwin//'.bcflag', ini, bcsurfflagtwo, 1)
1634 
1635 
1636  bcsurfflagtwo(1) = 2
1637 ! connectivity
1638  nullify(ptrhex8)
1639  testnumsurfel = 0
1640 
1641  ptrhex8 => surfmesh_hex8_s_head
1642  do while(associated(ptrhex8))
1643  surfaceelonproc = .false.
1644  IF(epart(ptrhex8%ElemData(5)).EQ.iprocs) surfaceelonproc = .true.
1645  IF(surfaceelonproc)THEN
1646 ! WRITE(4002,'(10i10)') ABS(NodeFlag(ptrhex8%ElemData(1:4))),&
1647 ! ElFlag(ptrhex8%ElemData(5)),1,1 ! last two are not really used
1648  testnumsurfel = testnumsurfel + 1
1649 
1650  elconntable(1:4,testnumsurfel) = abs(nodeflag(ptrhex8%ElemData(1:4)))
1651  ENDIF
1652  ptrhex8 => ptrhex8%next
1653  ENDDO
1654 
1655  IF(testnumsurfel.NE.numelhex2d(2,iprocs))THEN
1656  print*,'Error: Number of Solid Surface hex elements in linked list'
1657  print*,' does not match that given by NumElhex2D(2,iProcs)'
1658  print*,' In linked list =',testnumsurfel
1659  print*,' NumElhex2D(2,iProcs) =',numelhex2d(2,iprocs)
1660  print*,'Stopping'
1661  ENDIF
1662 
1663  nullify(ptrhex8)
1664 
1665  IF(associated(nodeflag2d)) deallocate(nodeflag2d)
1666 
1667  ENDIF
1668 
1669  ! write surface window ------------------------
1670  CALL com_call_function( set_option, 2, 'mode', 'a')
1671 
1672  sur_all = com_get_attribute_handle( surwin//'.all')
1673 
1674  CALL com_call_function( write_attr, 4, 'Rocin/SurfMesh.'//ichr4, sur_all,&
1675  "isolid","00.000000")
1676 
1677 
1678  CALL com_delete_window( surwin)
1679 
1680 
1681  IF(associated(nodeflag2d)) deallocate(nodeflag2d)
1682 
1683 
1684 
1685 end subroutine mesh2d
1686 
void set_option(const char *option_name, const char *option_val)
Set an option for Rocout, such as controlling the output format.
Definition: Rocout.C:552
j indices k indices k
Definition: Indexing.h:6
static void write_attr(std::ostream &os, const COM::Attribute *attr, int i)
blockLoc i
Definition: read.cpp:79
subroutine mesh2d(nprocs, iProcs, ichr4)
virtual std::ostream & print(std::ostream &os) const