1 | include(dom.inc)
2 |
3 | SUBROUTINE writeinfiles(path, ndirs, s)
4 |
5 | ! ===============================================================!
6 | ! !
7 | ! writeinfiles.F : Writes out the data base of the geometry. !
8 | ! Format: V2.04 !
9 | ! !
10 | ! in : The complete Data Base !
11 | ! !
12 | ! author : J. AMAYA (june 2007) !
13 | ! !
14 | ! ===============================================================!
15 |
16 | USE datas
17 |
18 | IMPLICIT NONE
19 |
20 | INCLUDE 'dom_constants.h'
21 |
22 | ! IN
23 | CHARACTER*80 :: path
24 | DOM_INT :: ndirs
25 | DOM_REAL :: s(4, ndirs)
26 |
27 | ! OUT
28 |
29 | ! LOCAL
30 | DOM_INT, allocatable, dimension(:) :: fatnodelist
31 | DOM_INT, allocatable, dimension(:) :: neighbour
32 | DOM_INT, allocatable, dimension(:) :: pathway
33 | DOM_INT, allocatable, dimension(:) :: added
34 | DOM_INT, allocatable, dimension(:) :: knownf
35 | DOM_INT, allocatable, dimension(:) :: ninfaces
36 | DOM_INT, allocatable, dimension(:) :: bcfaces
37 | DOM_INT, allocatable, dimension(:) :: cellnodes
38 | DOM_INT, allocatable, dimension(:) :: cellfaces
39 | DOM_INT, allocatable, dimension(:) :: patch
40 | DOM_INT :: i, j, i_fatnodes, parent2,parent
41 | DOM_INT :: nb_faces, thiscell_id, nbdyfaces
42 | DOM_INT :: ndim, nbnodes, thisnode_nb
43 | DOM_INT :: thisface_id, k, neparent
44 | DOM_INT :: cellcount, myneigh
45 | DOM_INT :: pathbeg, pathend, ipath
46 | DOM_REAL :: xc, yc, zc
47 | DOM_REAL :: kabs_gray, k_scat, Dns
48 | DOM_REAL :: nmin, Tw
49 | DOM_REAL, allocatable, dimension(:) :: n1, n2, n3
50 | DOM_REAL, allocatable, dimension(:) :: xfc, yfc, zfc
51 | DOM_REAL, allocatable, dimension(:) :: area
52 | character*80 :: c2cfile, ccellsfile, cfacesfile
53 | character*80 :: normfile, volafile, gdatafile
54 | character*80 :: quadfile, clfacefile, clpropfile
55 | character*80 :: progfile, nodesfile, clnodefile
56 | character*80 :: c2facfile, facfile
57 | type(ptr_cell), allocatable, dimension(:) :: pt_pathway
58 | type(ptr_cell), allocatable, dimension(:) :: pt_cellist
59 | type(cell), pointer :: pt_cell, auxcell_ptr
60 | type(face), pointer :: pt_face
61 | type(fatnode), pointer :: pt_fatnode
62 |
63 | ! --------------------------------------!
64 | ! Set file names and open them to write !
65 | ! --------------------------------------!
66 |
67 | nodesfile = path(1:len_trim(path))//'/Nodelist.in'
68 | c2cfile = path(1:len_trim(path))//'/Cell2cells.in'
69 | c2facfile = path(1:len_trim(path))//'/Cell2faces.in'
70 | facfile = path(1:len_trim(path))//'/Facelist.in'
71 | ccellsfile = path(1:len_trim(path))//'/Centercells.in'
72 | cfacesfile = path(1:len_trim(path))//'/Centerfaces.in'
73 | clnodefile = path(1:len_trim(path))//'/Cellnodes.in'
74 | normfile = path(1:len_trim(path))//'/Normals.in'
75 | volafile = path(1:len_trim(path))//'/Volumesareas.in'
76 | clpropfile = path(1:len_trim(path))//'/CLProperties.in'
77 | clfacefile = path(1:len_trim(path))//'/CLFaces.in'
78 | quadfile = path(1:len_trim(path))//'/Quadrature.in'
79 | progfile = path(1:len_trim(path))//'/Progress.in'
80 | gdatafile = path(1:len_trim(path))//'/Global.in'
81 |
82 | OPEN(FILE_NODES, FILE=nodesfile , FORM='UNFORMATTED')
83 | OPEN(FILE_C2C , FILE=c2cfile , FORM='UNFORMATTED')
84 | OPEN(FILE_C2FAC, FILE=c2facfile , FORM='UNFORMATTED')
85 | OPEN(FILE_FAC , FILE=facfile , FORM='UNFORMATTED')
86 | OPEN(FILE_CCELL, FILE=ccellsfile, FORM='UNFORMATTED')
87 | OPEN(FILE_CFACE, FILE=cfacesfile, FORM='UNFORMATTED')
88 | OPEN(FILE_CLNOD, FILE=clnodefile, FORM='UNFORMATTED')
89 | OPEN(FILE_NORM , FILE=normfile , FORM='UNFORMATTED')
90 | OPEN(FILE_VOLA , FILE=volafile , FORM='UNFORMATTED')
91 | OPEN(FILE_CLPRO, FILE=clpropfile, FORM='UNFORMATTED')
92 | OPEN(FILE_CLFAC, FILE=clfacefile, FORM='UNFORMATTED')
93 | OPEN(FILE_QUADR, FILE=quadfile , FORM='UNFORMATTED')
94 | OPEN(FILE_PROG , FILE=progfile , FORM='UNFORMATTED')
95 | OPEN(FILE_GDATA, FILE=gdatafile , FORM='UNFORMATTED')
96 |
97 | nbdyfaces = 0
98 |
99 | ! ---------------------------!
100 | ! Write 'Quadrature.in' file !
101 | ! ---------------------------!
102 |
103 | WRITE(*,*) " Writing Quadrature..."
104 | DO i=1, ndirs
105 | WRITE(FILE_QUADR) (s(j,i),j=1,4)
106 | ENDDO
107 |
108 | CLOSE(FILE_QUADR)
109 |
110 | ! -------------------------!
111 | ! Write 'Progress.in' file !
112 | ! -------------------------!
113 |
114 | WRITE(*,*) " Writing Pathway over all directions:"
115 |
116 | IF (ALLOCATED(pathway)) DEALLOCATE(pathway)
117 | IF (ALLOCATED(added)) DEALLOCATE(added)
118 | IF (ALLOCATED(ninfaces)) DEALLOCATE(ninfaces)
119 | IF (ALLOCATED(bcfaces)) DEALLOCATE(bcfaces)
120 | IF (ALLOCATED(knownf)) DEALLOCATE(knownf)
121 |
122 | ALLOCATE( pathway(i_ncells))
123 | ALLOCATE( added(i_ncells))
124 | ALLOCATE( knownf(i_ncells))
125 | ALLOCATE( ninfaces(i_ncells))
126 | ALLOCATE( bcfaces(i_ncells))
127 | ALLOCATE(pt_pathway(i_ncells))
128 | ALLOCATE(pt_cellist(i_ncells))
129 |
130 | pathway = 0
131 | pt_cell => cell_list
132 | DO i = 1, i_ncells
133 | NULLIFY(pt_pathway(1)%cell_ptr)
134 | pt_cellist(i)%cell_ptr => pt_cell
135 | pt_cell => pt_cell%next_cell
136 | ENDDO
137 |
138 | ! ------------------------------!
139 | ! Do one register per direction !
140 | ! ------------------------------!
141 | DO i=1, ndirs
142 |
143 | added = 0
144 | cellcount = 0
145 | knownf = 0
146 | ninfaces = 0
147 | bcfaces = 0
148 | pathbeg = 1
149 |
150 | WRITE(*,'(1x,A,i3,A,1x,f6.4,1x,f6.4,1x,f6.4)') &
151 | & " direction ",i,":", s(1:3,i)
152 |
153 | ! -------------------------------------------------------!
154 | ! Detect 'in' and 'in'+boundary faces for this direction !
155 | ! -------------------------------------------------------!
156 | pt_cell => cell_list
157 |
158 | DO WHILE (ASSOCIATED(pt_cell))
159 |
160 | j = pt_cell%cell_id
161 |
162 | DO k=1, pt_cell%i_nbfaces
163 |
164 | pt_face => pt_cell%cell_face(k)%face_ptr
165 | parent = 1
166 | IF (pt_face%parent_cell(1).ne.j) parent = 2
167 |
168 | ! ---------------------------------------!
169 | ! If face is in, normal*direction is < 0 !
170 | ! ---------------------------------------!
171 | Dns = pt_face%face_normal(1,parent)*s(1,i) + &
172 | & pt_face%face_normal(2,parent)*s(2,i) + &
173 | & pt_face%face_normal(3,parent)*s(3,i)
174 |
175 | IF (Dns.lt.0.) THEN
176 |
177 | ninfaces(j) = ninfaces(j) + 1
178 |
179 | IF (pt_face%i_patchnb.ne.0) THEN
180 | bcfaces(j) = bcfaces(j) + 1
181 | ENDIF
182 |
183 | ENDIF
184 |
185 | ENDDO
186 |
187 | ! --------------------------------------------!
188 | ! Add the first 'row' of cells to the pathway !
189 | ! --------------------------------------------!
190 |
191 | IF ( (added(j).eq.0).and.(bcfaces(j).eq.ninfaces(j)) ) THEN
192 |
193 | added(j) = 1
194 | cellcount = cellcount + 1
195 | pathway(cellcount) = j
196 | pt_pathway(cellcount)%cell_ptr => pt_cell
197 |
198 | ENDIF
199 |
200 | pt_cell => pt_cell%next_cell
201 |
202 | ENDDO
203 |
204 | pathend = cellcount
205 |
206 | ! ---------------------------------------!
207 | ! Fill pathway vector in multiple cycles !
208 | ! ---------------------------------------!
209 | DO WHILE(cellcount.lt.i_ncells)
210 |
211 | DO ipath = pathbeg, pathend
212 |
213 | pt_cell => pt_pathway(ipath)%cell_ptr
214 |
215 | j = pt_cell%cell_id
216 |
217 | DO k=1, pt_cell%i_nbfaces
218 |
219 | pt_face => pt_cell%cell_face(k)%face_ptr
220 | parent = 1
221 | IF (pt_face%parent_cell(1).ne.j) parent = 2
222 |
223 | ! ---------------------------------------!
224 | ! If face is 'out', normal*direction > 0 !
225 | ! ---------------------------------------!
226 | Dns = pt_face%face_normal(1,parent)*s(1,i) + &
227 | & pt_face%face_normal(2,parent)*s(2,i) + &
228 | & pt_face%face_normal(3,parent)*s(3,i)
229 |
230 | IF ((Dns.gt.nmin).and.(pt_face%i_patchnb.eq.0)) &
231 | & THEN
232 | IF (pt_face%parent_cell(1).eq.j) THEN
233 | myneigh = pt_face%parent2_ptr%cell_id
234 | ELSE
235 | myneigh = pt_face%parent1_ptr%cell_id
236 | ENDIF
237 | knownf(myneigh) = knownf(myneigh) + 1
238 | ENDIF
239 |
240 | IF ( (added(myneigh).eq.0).and. &
241 | & (knownf(myneigh)+bcfaces(myneigh).eq. &
242 | & ninfaces(myneigh)) ) THEN
243 |
244 | added(myneigh) = 1
245 | cellcount = cellcount + 1
246 | pathway(cellcount) = myneigh
247 | pt_pathway(cellcount)%cell_ptr => &
248 | & pt_cellist(myneigh)%cell_ptr
249 |
250 | ENDIF
251 |
252 | ENDDO
253 |
254 | ENDDO
255 |
256 | pathbeg = pathend + 1
257 | pathend = cellcount
258 |
259 | ENDDO
260 |
261 | WRITE(FILE_PROG) (pathway(j),j=1,i_ncells)
262 |
263 | ENDDO
264 |
265 | DEALLOCATE( pathway)
266 | DEALLOCATE( added)
267 | DEALLOCATE( knownf)
268 | DEALLOCATE( ninfaces)
269 | DEALLOCATE( bcfaces)
270 | DEALLOCATE(pt_pathway)
271 | DEALLOCATE(pt_cellist)
272 |
273 | CLOSE(FILE_PROG)
274 |
275 | ! -------------------------!
276 | ! Write 'Facelist.in' file !
277 | ! -------------------------!
278 |
279 | pt_face => face_list
280 | DO i=1, i_nfaces
281 |
282 | IF (pt_face%i_nbparents.eq.2) THEN
283 | IF (pt_face%parent2_ptr%cell_id.eq.0) THEN
284 | WRITE(*,*) " Fatal error: face ", pt_face%face_id, &
285 | & " has lost one parent"
286 | ENDIF
287 | ENDIF
288 |
289 | WRITE(FILE_FAC) pt_face%face_id, pt_face%i_nbnodes, &
290 | & (pt_face%face_point(j,1), &
291 | & j=1,pt_face%i_nbnodes)
292 | pt_face => pt_face%next_face
293 | ENDDO
294 |
295 | ! -------------------------!
296 | ! Write 'Nodelist.in' file !
297 | ! -------------------------!
298 |
299 | DO i=1, i_nnodes
300 | WRITE(FILE_NODES) i, (node_list(j,i),j=1,3)
301 | ENDDO
302 |
303 | ! -------------------------------------------!
304 | ! Save global mesh and initial solution data !
305 | ! -------------------------------------------!
306 |
307 | WRITE(FILE_GDATA) i_ncells,i_nnodes,i_nfaces,i_nfacesmax,ndirs
308 |
309 | ! --------------------!
310 | ! Loop over all cells !
311 | ! --------------------!
312 |
313 | ALLOCATE(neighbour(1:2*i_nfacesmax))
314 | ALLOCATE(cellnodes(MAX_NNODES_CELL))
315 | ALLOCATE(cellfaces(MAX_NFACES_CELL))
316 | ALLOCATE(xfc (MAX_NFACES_CELL))
317 | ALLOCATE(yfc (MAX_NFACES_CELL))
318 | ALLOCATE(zfc (MAX_NFACES_CELL))
319 | ALLOCATE(n1 (MAX_NFACES_CELL))
320 | ALLOCATE(n2 (MAX_NFACES_CELL))
321 | ALLOCATE(n3 (MAX_NFACES_CELL))
322 | ALLOCATE(area (MAX_NFACES_CELL))
323 | ALLOCATE(patch (MAX_NFACES_CELL))
324 |
325 | WRITE(*,*) " Looping over all the cells", i_ncells
326 |
327 | pt_cell => cell_list
328 |
329 | DO i=1, i_ncells
330 |
331 | thiscell_id = pt_cell%cell_id
332 | nb_faces = pt_cell%i_nbfaces
333 | neighbour = 0
334 | ! print*, "C: ", thiscell_id, nb_faces
335 |
336 | ! -------------------------!
337 | ! Detect cell's properties !
338 | ! -------------------------!
339 |
340 | SELECTCASE(pt_cell%celltype)
341 | CASE(EL_TRI)
342 | ndim = 2
343 | nbnodes = 3
344 | CASE(EL_QUAD)
345 | ndim = 2
346 | nbnodes = 4
347 | CASE(EL_TETRA)
348 | ndim = 3
349 | nbnodes = 4
350 | CASE(EL_PYRAM)
351 | ndim = 3
352 | nbnodes = 5
353 | CASE(EL_PRISM)
354 | ndim = 3
355 | nbnodes = 6
356 | CASE(EL_HEXA)
357 | ndim = 3
358 | nbnodes = 8
359 | CASE DEFAULT
360 | WRITE(*,*) "Error: element type not recognised ", &
361 | & pt_cell%celltype
362 | STOP
363 | ENDSELECT
364 |
365 | ! --------------------------------------------------------!
366 | ! Detect the neighbour cells (0 if no neigh. at the face) !
367 | ! --------------------------------------------------------!
368 |
369 | DO j=1, nb_faces
370 |
371 | pt_face => pt_cell%cell_face(j)%face_ptr
372 | thisface_id = pt_face%face_id
373 |
374 | IF (pt_face%i_nbparents.eq.2) THEN
375 | IF (pt_face%parent2_ptr%cell_id.eq.0) THEN
376 | WRITE(*,*) " Face id: ", thisface_id
377 | WRITE(*,*) " Fatal error: face ",pt_face%face_id, &
378 | & " has lost one parent"
379 | ENDIF
380 | ENDIF
381 |
382 |
383 | ! -----------------------------------!
384 | ! Detect the adjacent cell to face j !
385 | ! -----------------------------------!
386 | IF (pt_face%i_nbparents.eq.2) THEN
387 |
388 | neparent = 1
389 | auxcell_ptr => pt_face%parent1_ptr
390 |
391 | IF (pt_face%parent_cell(1).eq.thiscell_id) THEN
392 | neparent = 2
393 | auxcell_ptr => pt_face%parent2_ptr
394 | ENDIF
395 |
396 | ! print*, " nparent=",pt_face%i_nbparents
397 |
398 | IF (auxcell_ptr%cell_id.eq.0) THEN
399 | print*, " F: ",thisface_id,"neparent: ", neparent
400 | print*, " neighid=",auxcell_ptr%cell_id
401 | print*, " parents:", pt_face%i_nbparents
402 | print*, " parent1:", pt_face%parent1_ptr%cell_id
403 | print*, " parent2:", pt_face%parent2_ptr%cell_id
404 | WRITE(*,*) " Fatal error, face ",pt_face%face_id, &
405 | & " points to no parent cell"
406 | STOP
407 | ENDIF
408 |
409 | ! ------------------------------------------------!
410 | ! Set the face rank in the neighbour cell's array !
411 | ! ------------------------------------------------!
412 | neighbour(2*j-1) = pt_face%parent_cell(neparent)
413 | DO k=1, nb_faces
414 | ! print*, " k=", k
415 | ! print*, " neighid=",auxcell_ptr%cell_id
416 | ! print*, "id:",auxcell_ptr%cell_face(k)%face_ptr%face_id
417 | IF (auxcell_ptr%cell_face(k)%face_ptr% &
418 | & face_id.eq.thisface_id) EXIT
419 | ENDDO
420 | neighbour(2*j) = k
421 | ELSE
422 | ! print*, " boundary face!"
423 | nbdyfaces = nbdyfaces + 1
424 | neighbour(2*j-1) = 0
425 | neighbour (2*j) = 0
426 | ENDIF
427 |
428 | ENDDO
429 |
430 | ! ---------------------------!
431 | ! Write 'Cell2cells.in' file !
432 | ! ---------------------------!
433 |
434 | ! print*, " Cell2cells.in"
435 | WRITE(FILE_C2C) thiscell_id, &
436 | & nb_faces, &
437 | & neighbour(1:2*nb_faces)
438 |
439 | ! --------------------------------------------------!
440 | ! Fill cellnodes vector and calculate cell's center !
441 | ! --------------------------------------------------!
442 |
443 | cellnodes = 0
444 |
445 | xc = 0
446 | yc = 0
447 | zc = 0
448 |
449 | DO j=1, nbnodes
450 | thisnode_nb = pt_cell%cellnodes(j)
451 | cellnodes(j) = thisnode_nb
452 | xc = xc + node_list(1, thisnode_nb)
453 | yc = yc + node_list(2, thisnode_nb)
454 | zc = zc + node_list(3, thisnode_nb)
455 | ENDDO
456 |
457 | xc = xc / nbnodes
458 | yc = yc / nbnodes
459 | zc = zc / nbnodes
460 |
461 | ! ------------------------------------------------!
462 | ! Write 'Centercells.in' and 'Cellnodes.in' files !
463 | ! ------------------------------------------------!
464 |
465 | ! print*, " Centercells.in"
466 | WRITE(FILE_CCELL) thiscell_id, xc, yc, zc
467 | WRITE(FILE_CLNOD) thiscell_id, nbnodes, &
468 | & (cellnodes(j),j=1,nbnodes)
469 |
470 | ! --------------------!
471 | ! Compute faces' data !
472 | ! --------------------!
473 |
474 | cellfaces = 0
475 | xfc = 0.
476 | yfc = 0.
477 | zfc = 0.
478 | n1 = 0.
479 | n2 = 0.
480 | n3 = 0.
481 | area = 0.
482 | patch = 0
483 |
484 | DO j=1, nb_faces
485 |
486 | ! ----------------!
487 | ! Get face center !
488 | ! ----------------!
489 | pt_face => pt_cell%cell_face(j)%face_ptr
490 | xfc(j) = pt_face%x
491 | yfc(j) = pt_face%y
492 | zfc(j) = pt_face%z
493 |
494 | ! ------------------------------------------!
495 | ! Get face normal, area and B.C. properties !
496 | ! ------------------------------------------!
497 | parent = 1
498 | IF (pt_face%parent_cell(1).ne.thiscell_id) parent = 2
499 | n1(j) = pt_face%face_normal(1,parent)
500 | n2(j) = pt_face%face_normal(2,parent)
501 | n3(j) = pt_face%face_normal(3,parent)
502 |
503 | area (j) = pt_face%area
504 | cellfaces(j) = pt_face%face_id
505 | patch (j) = pt_face%i_patchnb
506 |
507 | ENDDO
508 |
509 | ! ------------------------------------------------!
510 | ! Write 'Centerfaces.in' and 'Cell2faces.in' file !
511 | ! ------------------------------------------------!
512 |
513 | WRITE(FILE_CFACE) thiscell_id, nb_faces, &
514 | & (xfc(j), yfc(j), zfc(j), &
515 | & j=1, nb_faces)
516 | WRITE(FILE_C2FAC) thiscell_id, nb_faces, &
517 | & (cellfaces(j),j=1, nb_faces)
518 |
519 | ! -------------------!
520 | ! Write 'Normals.in' !
521 | ! -------------------!
522 |
523 | WRITE(FILE_NORM) thiscell_id, nb_faces, &
524 | & (n1(j), n2(j), n3(j), &
525 | & j=1, nb_faces)
526 |
527 | ! ------------------------!
528 | ! Write 'Volumesareas.in' !
529 | ! ------------------------!
530 |
531 | WRITE(FILE_VOLA) thiscell_id,nb_faces,pt_cell%volume, &
532 | & (area(j),j=1, nb_faces)
533 |
534 | ! ------------------------!
535 | ! Write 'CLProperties.in' !
536 | ! ------------------------!
537 |
538 | Tw = 0.
539 | WRITE(FILE_CLPRO) thiscell_id, nb_faces, &
540 | & (Tw,patch(j),j=1, nb_faces)
541 | ! print*, thiscell_id,":",nb_faces,":",Tw,patch
542 |
543 | ! ----------------!
544 | ! Go to next cell !
545 | ! ----------------!
546 |
547 | pt_cell => pt_cell%next_cell
548 |
549 | ENDDO
550 |
551 | ! -------------------!
552 | ! Write 'CLFaces.in' !
553 | ! -------------------!
554 |
555 | WRITE(FILE_CLFAC) nbdyfaces
556 |
557 | ! ----------------!
558 | ! Close all files !
559 | ! ----------------!
560 |
561 | WRITE(*,*) " Closing files"
562 |
563 | CLOSE(FILE_NODES)
564 | CLOSE(FILE_C2C )
565 | CLOSE(FILE_C2FAC)
566 | CLOSE(FILE_FAC )
567 | CLOSE(FILE_CCELL)
568 | CLOSE(FILE_CFACE)
569 | CLOSE(FILE_CLNOD)
570 | CLOSE(FILE_NORM )
571 | CLOSE(FILE_VOLA )
572 | CLOSE(FILE_CLPRO)
573 | CLOSE(FILE_CLFAC)
574 | CLOSE(FILE_CFACE)
575 | CLOSE(FILE_GDATA)
576 |
577 | END SUBROUTINE writeinfiles
writeinfiles.F could be called by: