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 | myneigh = 1
210 |
211 | DO WHILE(cellcount.lt.i_ncells)
212 |
213 | DO ipath = pathbeg, pathend
214 |
215 | pt_cell => pt_pathway(ipath)%cell_ptr
216 |
217 | j = pt_cell%cell_id
218 |
219 | DO k=1, pt_cell%i_nbfaces
220 |
221 | pt_face => pt_cell%cell_face(k)%face_ptr
222 | parent = 1
223 | IF (pt_face%parent_cell(1).ne.j) parent = 2
224 |
225 | ! ---------------------------------------!
226 | ! If face is 'out', normal*direction > 0 !
227 | ! ---------------------------------------!
228 | Dns = pt_face%face_normal(1,parent)*s(1,i) + &
229 | & pt_face%face_normal(2,parent)*s(2,i) + &
230 | & pt_face%face_normal(3,parent)*s(3,i)
231 |
232 | IF ((Dns.gt.nmin).and.(pt_face%i_patchnb.eq.0)) &
233 | & THEN
234 | IF (pt_face%parent_cell(1).eq.j) THEN
235 | myneigh = pt_face%parent2_ptr%cell_id
236 | ELSE
237 | myneigh = pt_face%parent1_ptr%cell_id
238 | ENDIF
239 | knownf(myneigh) = knownf(myneigh) + 1
240 | ENDIF
241 |
242 | IF ( (added(myneigh).eq.0).and. &
243 | & (knownf(myneigh)+bcfaces(myneigh).eq. &
244 | & ninfaces(myneigh)) ) THEN
245 |
246 | added(myneigh) = 1
247 | cellcount = cellcount + 1
248 | pathway(cellcount) = myneigh
249 | pt_pathway(cellcount)%cell_ptr => &
250 | & pt_cellist(myneigh)%cell_ptr
251 |
252 | ENDIF
253 |
254 | ENDDO
255 |
256 | ENDDO
257 |
258 | pathbeg = pathend + 1
259 | pathend = cellcount
260 |
261 | ENDDO
262 |
263 | WRITE(FILE_PROG) (pathway(j),j=1,i_ncells)
264 |
265 | ENDDO
266 |
267 | DEALLOCATE( pathway)
268 | DEALLOCATE( added)
269 | DEALLOCATE( knownf)
270 | DEALLOCATE( ninfaces)
271 | DEALLOCATE( bcfaces)
272 | DEALLOCATE(pt_pathway)
273 | DEALLOCATE(pt_cellist)
274 |
275 | CLOSE(FILE_PROG)
276 |
277 | ! -------------------------!
278 | ! Write 'Facelist.in' file !
279 | ! -------------------------!
280 |
281 | pt_face => face_list
282 | DO i=1, i_nfaces
283 |
284 | IF (pt_face%i_nbparents.eq.2) THEN
285 | IF (pt_face%parent2_ptr%cell_id.eq.0) THEN
286 | WRITE(*,*) " Fatal error: face ", pt_face%face_id, &
287 | & " has lost one parent"
288 | ENDIF
289 | ENDIF
290 |
291 | WRITE(FILE_FAC) pt_face%face_id, pt_face%i_nbnodes, &
292 | & (pt_face%face_point(j,1), &
293 | & j=1,pt_face%i_nbnodes)
294 | pt_face => pt_face%next_face
295 | ENDDO
296 |
297 | ! -------------------------!
298 | ! Write 'Nodelist.in' file !
299 | ! -------------------------!
300 |
301 | DO i=1, i_nnodes
302 | WRITE(FILE_NODES) i, (node_list(j,i),j=1,3)
303 | ENDDO
304 |
305 | ! -------------------------------------------!
306 | ! Save global mesh and initial solution data !
307 | ! -------------------------------------------!
308 |
309 | WRITE(FILE_GDATA) i_ncells,i_nnodes,i_nfaces,i_nfacesmax,ndirs
310 |
311 | ! --------------------!
312 | ! Loop over all cells !
313 | ! --------------------!
314 |
315 | ALLOCATE(neighbour(1:2*i_nfacesmax))
316 | ALLOCATE(cellnodes(MAX_NNODES_CELL))
317 | ALLOCATE(cellfaces(MAX_NFACES_CELL))
318 | ALLOCATE(xfc (MAX_NFACES_CELL))
319 | ALLOCATE(yfc (MAX_NFACES_CELL))
320 | ALLOCATE(zfc (MAX_NFACES_CELL))
321 | ALLOCATE(n1 (MAX_NFACES_CELL))
322 | ALLOCATE(n2 (MAX_NFACES_CELL))
323 | ALLOCATE(n3 (MAX_NFACES_CELL))
324 | ALLOCATE(area (MAX_NFACES_CELL))
325 | ALLOCATE(patch (MAX_NFACES_CELL))
326 |
327 | WRITE(*,*) " Looping over all the cells", i_ncells
328 |
329 | pt_cell => cell_list
330 |
331 | DO i=1, i_ncells
332 |
333 | thiscell_id = pt_cell%cell_id
334 | nb_faces = pt_cell%i_nbfaces
335 | neighbour = 0
336 | ! print*, "C: ", thiscell_id, nb_faces
337 |
338 | ! -------------------------!
339 | ! Detect cell's properties !
340 | ! -------------------------!
341 |
342 | SELECTCASE(pt_cell%celltype)
343 | CASE(EL_TRI)
344 | ndim = 2
345 | nbnodes = 3
346 | CASE(EL_QUAD)
347 | ndim = 2
348 | nbnodes = 4
349 | CASE(EL_TETRA)
350 | ndim = 3
351 | nbnodes = 4
352 | CASE(EL_PYRAM)
353 | ndim = 3
354 | nbnodes = 5
355 | CASE(EL_PRISM)
356 | ndim = 3
357 | nbnodes = 6
358 | CASE(EL_HEXA)
359 | ndim = 3
360 | nbnodes = 8
361 | CASE DEFAULT
362 | WRITE(*,*) "Error: element type not recognised ", &
363 | & pt_cell%celltype
364 | STOP
365 | ENDSELECT
366 |
367 | ! --------------------------------------------------------!
368 | ! Detect the neighbour cells (0 if no neigh. at the face) !
369 | ! --------------------------------------------------------!
370 |
371 | DO j=1, nb_faces
372 |
373 | pt_face => pt_cell%cell_face(j)%face_ptr
374 | thisface_id = pt_face%face_id
375 |
376 | IF (pt_face%i_nbparents.eq.2) THEN
377 | IF (pt_face%parent2_ptr%cell_id.eq.0) THEN
378 | WRITE(*,*) " Face id: ", thisface_id
379 | WRITE(*,*) " Fatal error: face ",pt_face%face_id, &
380 | & " has lost one parent"
381 | ENDIF
382 | ENDIF
383 |
384 |
385 | ! -----------------------------------!
386 | ! Detect the adjacent cell to face j !
387 | ! -----------------------------------!
388 | IF (pt_face%i_nbparents.eq.2) THEN
389 |
390 | neparent = 1
391 | auxcell_ptr => pt_face%parent1_ptr
392 |
393 | IF (pt_face%parent_cell(1).eq.thiscell_id) THEN
394 | neparent = 2
395 | auxcell_ptr => pt_face%parent2_ptr
396 | ENDIF
397 |
398 | ! print*, " nparent=",pt_face%i_nbparents
399 |
400 | IF (auxcell_ptr%cell_id.eq.0) THEN
401 | print*, " F: ",thisface_id,"neparent: ", neparent
402 | print*, " neighid=",auxcell_ptr%cell_id
403 | print*, " parents:", pt_face%i_nbparents
404 | print*, " parent1:", pt_face%parent1_ptr%cell_id
405 | print*, " parent2:", pt_face%parent2_ptr%cell_id
406 | WRITE(*,*) " Fatal error, face ",pt_face%face_id, &
407 | & " points to no parent cell"
408 | STOP
409 | ENDIF
410 |
411 | ! ------------------------------------------------!
412 | ! Set the face rank in the neighbour cell's array !
413 | ! ------------------------------------------------!
414 | neighbour(2*j-1) = pt_face%parent_cell(neparent)
415 | DO k=1, nb_faces
416 | ! print*, " k=", k
417 | ! print*, " neighid=",auxcell_ptr%cell_id
418 | ! print*, "id:",auxcell_ptr%cell_face(k)%face_ptr%face_id
419 | IF (auxcell_ptr%cell_face(k)%face_ptr% &
420 | & face_id.eq.thisface_id) EXIT
421 | ENDDO
422 | neighbour(2*j) = k
423 | ELSE
424 | ! print*, " boundary face!"
425 | nbdyfaces = nbdyfaces + 1
426 | neighbour(2*j-1) = 0
427 | neighbour (2*j) = 0
428 | ENDIF
429 |
430 | ENDDO
431 |
432 | ! ---------------------------!
433 | ! Write 'Cell2cells.in' file !
434 | ! ---------------------------!
435 |
436 | ! print*, " Cell2cells.in"
437 | WRITE(FILE_C2C) thiscell_id, &
438 | & nb_faces, &
439 | & neighbour(1:2*nb_faces)
440 |
441 | ! --------------------------------------------------!
442 | ! Fill cellnodes vector and calculate cell's center !
443 | ! --------------------------------------------------!
444 |
445 | cellnodes = 0
446 |
447 | xc = 0
448 | yc = 0
449 | zc = 0
450 |
451 | DO j=1, nbnodes
452 | thisnode_nb = pt_cell%cellnodes(j)
453 | cellnodes(j) = thisnode_nb
454 | xc = xc + node_list(1, thisnode_nb)
455 | yc = yc + node_list(2, thisnode_nb)
456 | zc = zc + node_list(3, thisnode_nb)
457 | ENDDO
458 |
459 | xc = xc / nbnodes
460 | yc = yc / nbnodes
461 | zc = zc / nbnodes
462 |
463 | ! ------------------------------------------------!
464 | ! Write 'Centercells.in' and 'Cellnodes.in' files !
465 | ! ------------------------------------------------!
466 |
467 | ! print*, " Centercells.in"
468 | WRITE(FILE_CCELL) thiscell_id, xc, yc, zc
469 | WRITE(FILE_CLNOD) thiscell_id, nbnodes, &
470 | & (cellnodes(j),j=1,nbnodes)
471 |
472 | ! --------------------!
473 | ! Compute faces' data !
474 | ! --------------------!
475 |
476 | cellfaces = 0
477 | xfc = 0.
478 | yfc = 0.
479 | zfc = 0.
480 | n1 = 0.
481 | n2 = 0.
482 | n3 = 0.
483 | area = 0.
484 | patch = 0
485 |
486 | DO j=1, nb_faces
487 |
488 | ! ----------------!
489 | ! Get face center !
490 | ! ----------------!
491 | pt_face => pt_cell%cell_face(j)%face_ptr
492 | xfc(j) = pt_face%x
493 | yfc(j) = pt_face%y
494 | zfc(j) = pt_face%z
495 |
496 | ! ------------------------------------------!
497 | ! Get face normal, area and B.C. properties !
498 | ! ------------------------------------------!
499 | parent = 1
500 | IF (pt_face%parent_cell(1).ne.thiscell_id) parent = 2
501 | n1(j) = pt_face%face_normal(1,parent)
502 | n2(j) = pt_face%face_normal(2,parent)
503 | n3(j) = pt_face%face_normal(3,parent)
504 |
505 | area (j) = pt_face%area
506 | cellfaces(j) = pt_face%face_id
507 | patch (j) = pt_face%i_patchnb
508 |
509 | ENDDO
510 |
511 | ! ------------------------------------------------!
512 | ! Write 'Centerfaces.in' and 'Cell2faces.in' file !
513 | ! ------------------------------------------------!
514 |
515 | WRITE(FILE_CFACE) thiscell_id, nb_faces, &
516 | & (xfc(j), yfc(j), zfc(j), &
517 | & j=1, nb_faces)
518 | WRITE(FILE_C2FAC) thiscell_id, nb_faces, &
519 | & (cellfaces(j),j=1, nb_faces)
520 |
521 | ! -------------------!
522 | ! Write 'Normals.in' !
523 | ! -------------------!
524 |
525 | WRITE(FILE_NORM) thiscell_id, nb_faces, &
526 | & (n1(j), n2(j), n3(j), &
527 | & j=1, nb_faces)
528 |
529 | ! ------------------------!
530 | ! Write 'Volumesareas.in' !
531 | ! ------------------------!
532 |
533 | WRITE(FILE_VOLA) thiscell_id,nb_faces,pt_cell%volume, &
534 | & (area(j),j=1, nb_faces)
535 |
536 | ! ------------------------!
537 | ! Write 'CLProperties.in' !
538 | ! ------------------------!
539 |
540 | Tw = 0.
541 | WRITE(FILE_CLPRO) thiscell_id, nb_faces, &
542 | & (Tw,patch(j),j=1, nb_faces)
543 | ! print*, thiscell_id,":",nb_faces,":",Tw,patch
544 |
545 | ! ----------------!
546 | ! Go to next cell !
547 | ! ----------------!
548 |
549 | pt_cell => pt_cell%next_cell
550 |
551 | ENDDO
552 |
553 | ! -------------------!
554 | ! Write 'CLFaces.in' !
555 | ! -------------------!
556 |
557 | WRITE(FILE_CLFAC) nbdyfaces
558 |
559 | ! ----------------!
560 | ! Close all files !
561 | ! ----------------!
562 |
563 | WRITE(*,*) " Closing files"
564 |
565 | CLOSE(FILE_NODES)
566 | CLOSE(FILE_C2C )
567 | CLOSE(FILE_C2FAC)
568 | CLOSE(FILE_FAC )
569 | CLOSE(FILE_CCELL)
570 | CLOSE(FILE_CFACE)
571 | CLOSE(FILE_CLNOD)
572 | CLOSE(FILE_NORM )
573 | CLOSE(FILE_VOLA )
574 | CLOSE(FILE_CLPRO)
575 | CLOSE(FILE_CLFAC)
576 | CLOSE(FILE_CFACE)
577 | CLOSE(FILE_GDATA)
578 |
579 | END SUBROUTINE writeinfiles
writeinfiles.F could be called by: