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