Choral
mesh_interfaces.F90
Go to the documentation of this file.
1 
29 
31 
32 
35  use real_type
36  use algebra_set
37  use basic_tools
38  use graph_mod
39  use cell_mod
40  use mesh_mod
41 
42  implicit none
43  private
44 
45 
46  ! %----------------------------------------%
47  ! | |
48  ! | PUBLIC DATA |
49  ! | |
50  ! %----------------------------------------%
51  public :: define_interfaces
52  public :: define_domain_boundary !! TESTED
53 
54 contains
55 
56 
59  subroutine define_interfaces(m)
60  type(mesh), intent(inout) :: m
61 
62  if (choral_verb>0) write(*,*) &
63  & "mesh_interfaces : define_interfaces"
64 
65  select case(m%dim)
66  case(1)
67  call define_1dmesh_vertexes(m)
68 
69  case(2)
70  call define_2dmesh_edges(m)
71 
72  case(3)
73  call define_3dmesh_faces(m)
74 
75  end select
76 
77  if (choral_verb>1) write(*,*) &
78  &" Number of interfaces =", m%nbItf
79 
80  end subroutine define_interfaces
81 
82 
83 
90  subroutine define_domain_boundary(m)
91  type(mesh), intent(inout) :: m
92 
93  integer, dimension(:,:), allocatable :: tab, clTab
94  integer, dimension(:) , allocatable :: clType
95  integer, dimension(CELL_MAX_NBNODES) :: cl_nodes
96 
97  integer :: cpt1, cpt2, cpt3
98  integer :: ii, cl, idx, j1, j2, nNd
99 
100  if (choral_verb>1) write(*,*)&
101  &"mesh_interfaces : define_domain_boundary"
102 
103  if (.NOT.valid(m)) call quit(&
104  & "mesh_interfaces: define_domain_boundary: mesh 'm' not valid")
105 
106  if (m%nbItf<=0) call define_interfaces(m)
107 
108  call boundary_interface_scan(tab, m)
109 
110  !! Flag the boundary cells already present in the mesh
111  !! and checks wether boundary cells are missing
112  !!
113  cpt1 = 0
114  cpt2 = 0
115  do ii=1, m%nbItf
116 
117  !! internal interface
118  if (tab(1,ii)==0) cycle
119 
120  !! boundary interface
121  cl = tab(3,ii)
122  if (cl==0) then
123  !! interface ii is a boundary cell
124  !! that does not presently belong to the mesh
125  cpt2 = cpt2 + 1
126 
127  else
128  !! interface ii is a boundary cell
129  !! that belongs to the mesh
130  cpt1 = cpt1 +1
131 
132  end if
133 
134  end do
135 
136  !! All the boundary cells already belong to the mesh
137  !!
138  if (cpt2==0) then
139  call freemem(tab)
140 
141  !! Boundary cells need to be added
142  !! The mesh is rebuilt including these cells.
143  !!
144  !! For the new cells :
145  !! cell types will be in clType
146  !! cell nodes will be in clTab
147  !!
148  else
149 
150  call allocmem(cltab , cell_max_nbnodes, cpt2 )
151  call allocmem(cltype, cpt2 )
152  !!
153  cpt3 = 0
154  do ii=1, m%nbItf
155 
156  !! internal interface
157  if (tab(1,ii)==0) cycle
158 
159  !! interface ii is a boundary cell
160  !! that belongs to the mesh
161  if (tab(3,ii)/=0) cycle
162 
163  !! interface ii as a new cell
164  !! whose index will be 'cpt3'
165  cpt3 = cpt3 + 1
166 
167  !! the co-boundary of the interface 'ii'
168  !! is the cell 'cl' with 'hNd' nodes
169  !! and with nodes 'cl_nodes(1:nNd)'
170  cl = tab(1,ii)
171  j1 = m%clToNd%row(cl)
172  j2 = m%clToNd%row(cl+1)
173  nnd = j2-j1
174  cl_nodes(1:nnd) = m%clToNd%col(j1: j2-1)
175  !!
176  !! The interface 'ii' is the interface 'idx'
177  !! of the cell 'cl' for the local indexing
178  idx = tab(2,ii)
179  !!
180  select case(m%clType(cl))
181  case(cell_edg, cell_edg_2)
182  cltype(cpt3) = cell_vtx
183 
184  cltab(1,cpt3) = cl_nodes(idx)
185 
186  case(cell_trg)
187  cltype(cpt3) = cell_edg
188 
189  cltab(1:2,cpt3) = cl_nodes( cell_ed_vtx(1:2, idx, cell_trg) )
190 
191  case(cell_tet)
192  cltype(cpt3) = cell_trg
193 
194  cltab(1:3,cpt3) = cl_nodes( cell_fc_vtx(1:3, idx, cell_tet) )
195 
196  case default
197  call quit("mesh_interfaces: define_domain_boundary:&
198  & cell type not implemented: "&
199  & // trim(cell_name(tab(2,ii))))
200 
201  end select
202 
203  end do
204  call freemem(tab)
205 
206  !! Add the new cells
207  !!
208  call add_cells(m, cltype, cltab)
209 
210  end if
211 
212  if (choral_verb>1) write(*,*)&
213  &" Mesh domain boundary =", cpt1+cpt2, " cells"
214 
215  end subroutine define_domain_boundary
216 
243  subroutine boundary_interface_scan(tab, m)
244  type(mesh) , intent(inout) :: m
245  integer, dimension(:,:), allocatable :: tab
246 
247  integer :: itf, itf2, itf_idx, itf_nb_vtx
248  integer :: j1, j2, jj, ll
249  integer :: dim, cl, cl_type, wdt
250 
251  integer, dimension(CELL_MAX_FC_NBVTX) :: itf_vtx
252  integer, dimension(:), allocatable :: t0, t1, t2
253 
254  if (choral_verb>2) write(*,*)&
255  &"mesh_interfaces : boundary_interface_scan "
256 
257  if (m%nbItf<=0) call define_interfaces(m)
258 
259  if (.NOT.valid(m)) call quit(&
260  & "mesh_interfaces: boundary_interface_scan: mesh 'm' not valid")
261 
262  dim = m%dim
263  call allocmem(tab, 3, m%nbItf)
264  tab = 0
265 
266  !! max number of cell-neighbour for a node
267  !!
268  wdt = m%ndToCl%width
269  allocate(t0(wdt), t1(wdt), t2(wdt))
270 
271 
272  !! First determines the boundary interfaces
273  !! and tab(1,itf) = co-boundary of these interfaces
274  !!
275  do itf=1, m%itfToCl%nl
276  j1 = m%itfToCl%row(itf)
277  j2 = m%itfToCl%row(itf+1)
278  j2 = j2-j1
279 
280  !! boundary interface
281  if (j2==1) then
282  cl = m%itfToCl%col(j1)
283  tab(1, itf) = cl
284  end if
285 
286  end do
287 
288 
289  !! Second determines if there is a second cell
290  !! with vertexes the vertexes of a boundary interface
291  !!
292  do itf=1, m%itfToCl%nl
293  if (tab(1, itf) == 0) cycle
294 
295  !! interface itf co-boundary is cell cl
296  !!
297  cl = tab(1, itf)
298  cl_type = m%clType(cl)
299 
300  !! interface itf has index itf_idx
301  !! in the local numbering associated to cell cl
302  !!
303  j1 = m%clToItf%row(cl)
304  j2 = m%clToItf%row(cl+1) - 1
305  loop_1: do jj= j1, j2
306  itf2 = m%clToItf%col(jj)
307 
308  if (itf == itf2) then
309  itf_idx = jj-j1+1
310  exit loop_1
311  end if
312  end do loop_1
313  tab(2, itf) = itf_idx
314 
315  !! interface itf has vertexes itf_vtx(1:itf_nb_vtx)
316  !!
317  select case(dim)
318  case(1)
319  j1 = m%clToNd%row(cl) - 1
320  j2 = j1 + itf_idx
321  itf_vtx(1) = m%clToNd%col(j2)
322 
323  case(2)
324  itf_vtx(1:2) = cell_ed_vtx(:,itf_idx, cl_type)
325 
326  j1 = m%clToNd%row(cl) - 1
327 
328  j2 = j1 + itf_vtx(1)
329  itf_vtx(1) = m%clToNd%col(j2)
330 
331  j2 = j1 + itf_vtx(2)
332  itf_vtx(2) = m%clToNd%col(j2)
333 
334  case(3)
335  itf_nb_vtx = cell_fc_nbvtx(itf_idx, cl_type)
336 
337  itf_vtx(1:itf_nb_vtx) = &
338  & cell_fc_vtx(1:itf_nb_vtx,itf_idx, cl_type)
339 
340  j1 = m%clToNd%row(cl) - 1
341 
342  do jj=1, itf_nb_vtx
343  j2 = j1 + itf_vtx(jj)
344  itf_vtx(jj) = m%clToNd%col(j2)
345  end do
346 
347  end select
348 
349  !! determine what cell in the mesh has
350  !! the nodes itf_vtx(1:itf_nb_vtx) among its nodes
351  !!
352  select case(dim)
353  case(1)
354  call getrow(jj, t0, wdt, m%ndToCl, itf_vtx(1))
355 
356  case(2)
357  call getrow(j1, t1, wdt, m%ndToCl, itf_vtx(1))
358  call getrow(j2, t2, wdt, m%ndToCl, itf_vtx(2))
359 
360  call cap_sorted_set(jj, t0, wdt, t1(1:j1), j1, t2(1:j2), j2)
361 
362  case(3)
363  call getrow(j1, t1, wdt, m%ndToCl, itf_vtx(1))
364 
365  loop_2: do ll=2, itf_nb_vtx
366  call getrow(j2, t2, wdt, m%ndToCl, itf_vtx(ll))
367 
368  call cap_sorted_set(jj, t0, wdt, t1(1:j1), j1, t2(1:j2), j2)
369 
370 
371  if (ll==itf_nb_vtx) exit loop_2
372  j1 = jj
373  t1(1:j1) = t0(1:j1)
374 
375  end do loop_2
376 
377  end select
378 
379 #if DBG>0
380  if ((jj==0) .OR. (jj>2)) call quit(&
381  & "mesh_interfaces: boundary_interface_scan: 3")
382 #endif
383 
384  if (jj==2) then
385 
386  if (t0(1)==cl) then
387  tab(3, itf) = t0(2)
388  else
389  tab(3, itf) = t0(1)
390  end if
391 
392 #if DBG>0
393  j1 = m%clType(tab(3, itf) )
394  j2 = cell_dim( j1 )
395  if (j2 /= dim-1) call quit(&
396  & "mesh_interfaces: boundary_interface_scan: 4")
397 #endif
398  end if
399 
400  end do
401 
402  end subroutine boundary_interface_scan
403 
404  subroutine define_2dmesh_edges(m)
405  type(mesh), intent(inout) :: m
406 
407  integer, dimension(:,:), allocatable :: edTab
408  integer, dimension(:) , allocatable :: nnz
409  integer :: cl, j1, j2
410 
411  integer :: cl_t, nb_vtx, nb_ed, width, nb_itf
412  real(DP) :: t1, t2
413 
414  if (choral_verb>1) write(*,*) &
415  & "mesh_tools : define_2DMesh_edges"
416 
417  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
418  !!
419  !! Computation of the connectivity graph
420  !! cells --> interfaces
421  !!
422  t1 = clock()
423 
424  !! STEP 1: compute the graph pattern
425  !!
426  call allocmem(nnz, m%nbCl)
427  do cl=1, m%nbCl
428  cl_t = m%clType(cl)
429  if ( cell_dim(cl_t) ==2 ) then
430  nnz(cl) = cell_nbed(cl_t)
431  else
432  nnz(cl) = 0
433  end if
434  end do
435  m%clToItf = graph(nnz, m%nbCl, 0)
436  call freemem(nnz)
437 
438  !! STEP 2 : create edTab
439  !!
440  call allocmem(edtab, 3, m%clToItf%nnz)
441  width = m%ndToCl%width
442  do cl_t=1, cell_tot_nb
443  if ( cell_dim(cl_t) /= 2 ) cycle
444  if ( m%cell_count(cl_t) == 0 ) cycle
445 
446  nb_vtx = cell_nbvtx(cl_t)
447  nb_ed = cell_nbed(cl_t)
448  call cell_loop()
449 
450  end do
451 
452  !! STEP 3 : number interfaces
453  !!
454  call number_edges(nb_itf)
455  m%nbItf = nb_itf
456  m%clToItf%nc = nb_itf
457 
458  !! fill in the conectivity graph 2D-cells --> edges
459  !!
460  do cl=1, m%nbCl
461 
462  j1 = m%clToItf%row(cl)
463  j2 = m%clToItf%row(cl+1)-1
464 
465  m%clToItf%col(j1:j2) = edtab(3,j1:j2)
466 
467  end do
468 
469  call freemem(edtab)
470  t2 = clock() - t1
471  if (choral_verb>1) write(*,*) &
472  & " 2D-cells --> edges graph, CPU = ", &
473  & real(t2, sp)
474 
475  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
476  !!
477  !! Computation of the connectivity graph
478  !! interfaces --> cells
479  !!
480  t1 = clock()
481  call transpose(m%itfToCl, m%clToItf)
482  t2 = clock() - t1
483  if (choral_verb>1) write(*,*) &
484  & " edges --> 2D-cells graph, CPU = ", &
485  & real(t2, sp)
486 
487  contains
488 
489  !! Array of all the 2D-cell edges,
490  !! internal edges are counted twice
491  !!
492  !! Step 1 = get the co-boundary of all 2D-cell edgees
493  !!
494  subroutine cell_loop()
496  integer, dimension(nb_vtx) :: vtx_nbCl
497  integer, dimension(width, nb_vtx) :: vtx_cl
498  integer, dimension(nb_vtx) :: cl_vtx
499  integer, dimension(3) :: coBord
500  integer :: ed, v1, v2
501  integer :: l1, l2
502 
503  !$OMP PARALLEL PRIVATE(ed, v1, v2, j1, j2, l1, l2, vtx_nbCl, vtx_cl, cl_vtx, coBord)
504  !$OMP DO
505  do cl=1, m%nbCl
506  l1 = m%clType(cl)
507  if (l1 /= cl_t ) cycle
508 
509  !! cell cl vertexes
510  !!
511  j1 = m%clToNd%row(cl)
512  j2 = j1 + nb_vtx - 1
513  cl_vtx = m%clToNd%col( j1:j2 )
514 
515  !! neighbour cells for each vertex
516  !! of the cell cl
517  !!
518  do v1=1, nb_vtx
519  l1 = cl_vtx(v1)
520  j1 = m%ndToCl%row(l1)
521  j2 = m%ndToCl%row(l1+1)
522 
523  l1 = j2-j1
524  vtx_nbcl(v1) = l1
525 
526  !! vtx_cl is sorted !
527  vtx_cl(1:l1, v1) = m%ndToCl%col(j1:j2-1)
528 
529  end do
530 
531  !! build the array of edges
532  !!
533  j1 = m%clToItf%row(cl)
534  do ed=1, nb_ed
535 
536  edtab(1, j1) = cl
537 
538  !! list the cells sharing the vertexes v1 and v2
539  !! coBord is sorted !
540  v1 = cell_ed_vtx(1, ed, cl_t)
541  v2 = cell_ed_vtx(2, ed, cl_t)
542  l1 = vtx_nbcl(v1)
543  l2 = vtx_nbcl(v2)
544  call cap_sorted_set(j2, cobord, 3, &
545  & vtx_cl(1:l1, v1), l1, vtx_cl(1:l2, v2), l2)
546 
547  !! j2 = number of cells sharing this edge
548  select case(j2)
549 
550  case(1) !! boundary edge
551  edtab(2, j1) = cl
552 
553  case(2)
554  !! In this case the edge can be
555  !! either internal or not,
556 
557  if (cobord(1)==cl) then
558  !! it is not necessary to know whether
559  !! the edge is internal or not here
560  !!
561  edtab(2, j1) = cobord(2)
562 
563  else
564  !! it is necessary to know whether
565  !! the edge is internal or not here
566  !!
567  !! For this we check the dimension
568  !! of the second cell
569  l1 = cell_dim( m%clType(cobord(1)) )
570 
571  if ( l1==2 ) then !! internal edge
572  edtab(2, j1) = cobord(1)
573 
574  else !! boundary edge
575  edtab(2, j1) = cl
576 
577  end if
578  end if
579 
580  case(3)
581  !! In this case the edge is internal,
582  !! but one of the three cells in the
583  !! co-boundary is of dimension 1.
584  !!
585  l1 = cell_dim( m%clType(cobord(1)) )
586  if (l1==1) then
587 
588  if (cobord(2)==cl) then
589  edtab(2, j1) = cobord(3)
590  else
591  edtab(2, j1) = cobord(2)
592  end if
593 
594  else
595 
596  l1 = cell_dim( m%clType(cobord(2)) )
597  if (l1==1) then
598 
599  if (cobord(1)==cl) then
600  edtab(2, j1) = cobord(3)
601  else
602  edtab(2, j1) = cobord(1)
603  end if
604 
605  else
606 
607  if (cobord(1)==cl) then
608  edtab(2, j1) = cobord(2)
609  else
610  edtab(2, j1) = cobord(1)
611  end if
612 
613  end if
614  end if
615 
616  end select
617 
618  !! increment the edge
619  j1 = j1 + 1
620  end do
621  end do
622  !$OMP END DO
623  !$OMP END PARALLEL
624 
625  end subroutine cell_loop
626 
627 
628  !! Array of all the 2D-cell edges,
629  !! internal edges are counted twice
630  !!
631  !! Step 2 = global numberingg of all the 2D-cell edgees
632  !!
633  subroutine number_edges(cpt)
634  integer, intent(out) :: cpt
635 
636  integer :: jj, cl, cl2, ll
637 
638  cpt = 0
639  do jj=1, m%clToItf%nnz
640 
641  cl = edtab(1,jj)
642  cl2 = edtab(2,jj)
643 
644  if ( cl <= cl2 ) then
645 
646  cpt = cpt + 1
647  edtab(3,jj) = cpt
648 
649  else
650  ll = m%clToItf%row(cl2)
651 
652  do while( edtab(2,ll) /= cl )
653  ll = ll + 1
654  end do
655  edtab(3,jj) = edtab(3,ll)
656 
657  end if
658 
659  end do
660 
661  end subroutine number_edges
662 
663  end subroutine define_2dmesh_edges
664 
665 
666  subroutine define_3dmesh_faces(m)
667  type(mesh), intent(inout) :: m
668 
669  integer, dimension(:,:), allocatable :: fcTab
670  integer, dimension(:) , allocatable :: nnz
671  integer :: cl, j1, j2
672 
673  integer :: cl_t, nb_vtx, nb_fc, width, nb_itf
674  real(DP) :: t1, t2
675 
676  if (choral_verb>1) write(*,*) &
677  & "mesh_tools : define_3DMesh_faces"
678 
679  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
680  !!
681  !! Computation of the connectivity graph
682  !! cells --> interfaces
683  !!
684  t1 = clock()
685 
686  !! STEP 1: compute the graph pattern
687  !!
688  call allocmem(nnz, m%nbCl)
689  do cl=1, m%nbCl
690  cl_t = m%clType(cl)
691  if ( cell_dim(cl_t) ==3 ) then
692  nnz(cl) = cell_nbfc(cl_t)
693  else
694  nnz(cl) = 0
695  end if
696  end do
697  m%clToItf = graph(nnz, m%nbCl, 0)
698  call freemem(nnz)
699 
700  !! STEP 2 : create fcTab
701  !!
702  call allocmem(fctab, 3, m%clToItf%nnz)
703  width = m%ndToCl%width
704  do cl_t=1, cell_tot_nb
705  if ( cell_dim(cl_t) /= 3 ) cycle
706  if ( m%cell_count(cl_t) == 0 ) cycle
707 
708  nb_vtx = cell_nbvtx(cl_t)
709  nb_fc = cell_nbfc(cl_t)
710  call cell_loop()
711 
712  end do
713 
714  !! STEP 3 : number faces
715  !!
716  call number_faces(nb_itf)
717  m%nbItf = nb_itf
718  m%clToItf%nc = nb_itf
719 
720  !! fill in the conectivity graph 2D-cells --> edges
721  !!
722  do cl=1, m%nbCl
723 
724  j1 = m%clToItf%row(cl)
725  j2 = m%clToItf%row(cl+1)-1
726 
727  m%clToItf%col(j1:j2) = fctab(3,j1:j2)
728 
729  end do
730 
731  call freemem(fctab)
732  t2 = clock() - t1
733  if (choral_verb>1) write(*,*) &
734  & " cells --> faces graph, CPU = ", &
735  & real(t2, sp)
736 
737  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
738  !!
739  !! Computation of the connectivity graph
740  !! interfaces --> cells
741  !!
742  t1 = clock()
743  call transpose(m%itfToCl, m%clToItf)
744  t2 = clock() - t1
745  if (choral_verb>1) write(*,*) &
746  & " faces --> cells graph, CPU = ", &
747  & real(t2, sp)
748 
749  contains
750 
751  !! Array of all the 2D-cell edges,
752  !! internal edges are counted twice
753  !!
754  !! Step 1 = get the co-boundary of all 2D-cell edgees
755  !!
756  subroutine cell_loop()
757 
758  integer, dimension(nb_vtx) :: vtx_nbCl
759  integer, dimension(width, nb_vtx) :: vtx_cl
760  integer, dimension(nb_vtx) :: cl_vtx
761  integer, dimension(width) :: coBord, aux
762  integer :: fc, v1, v2
763  integer :: l1, l2, ii
764 
765  !$OMP PARALLEL PRIVATE(fc, v1, v2, j1, j2, l1, l2, vtx_nbCl, vtx_cl, cl_vtx, coBord, aux, ii)
766  !$OMP DO
767  do cl=1, m%nbCl
768  l1 = m%clType(cl)
769  if (l1 /= cl_t ) cycle
770 
771  !! cell cl vertexes
772  !!
773  j1 = m%clToNd%row(cl)
774  j2 = j1 + nb_vtx - 1
775  cl_vtx = m%clToNd%col( j1:j2 )
776 
777  !! neighbour cells for each vertex
778  !! of the cell cl
779  !!
780  do v1=1, nb_vtx
781  l1 = cl_vtx(v1)
782  j1 = m%ndToCl%row(l1)
783  j2 = m%ndToCl%row(l1+1)
784 
785  l1 = j2-j1
786  vtx_nbcl(v1) = l1
787 
788  !! vtx_cl is sorted !
789  vtx_cl(1:l1, v1) = m%ndToCl%col(j1:j2-1)
790 
791  end do
792 
793  !! build the array of faces
794  !!
795  j1 = m%clToItf%row(cl)
796  do fc=1, nb_fc
797 
798  fctab(1, j1) = cl
799 
800  !! list the cells sharing the vertexes of the
801  !! face fc
802  v1 = cell_fc_vtx(1, fc, cl_t)
803  v2 = cell_fc_vtx(2, fc, cl_t)
804  l1 = vtx_nbcl(v1)
805  l2 = vtx_nbcl(v2)
806  call cap_sorted_set(j2, cobord, width, &
807  & vtx_cl(1:l1, v1), l1, vtx_cl(1:l2, v2), l2)
808 
809  do ii=3, cell_fc_nbvtx(fc, cl_t)
810  v2 = cell_fc_vtx(ii, fc, cl_t)
811  l2 = vtx_nbcl(v2)
812 
813  l1 = j2
814  aux(1:l1) = cobord(1:l1)
815 
816  call cap_sorted_set(j2, cobord(1:l1), l1, &
817  & aux(1:l1), l1, vtx_cl(1:l2, v2), l2)
818 
819  end do
820 
821  !! j2 = number of cells sharing this face
822  select case(j2)
823  case(1) !! boundary face
824  fctab(2, j1) = cl
825 
826  case(2)
827  !! In this case the face can be
828  !! either internal or not,
829 
830  if (cobord(1)==cl) then
831  !! it is not necessary to know whether
832  !! the face is internal or not here
833  !!
834  fctab(2, j1) = cobord(2)
835 
836  else
837  !! it is necessary to know whether
838  !! the face is internal or not here
839  !!
840  !! For this we check the dimension
841  !! of the second cell
842  l1 = cell_dim( m%clType(cobord(1)) )
843 
844  if ( l1==3 ) then !! internal edge
845  fctab(2, j1) = cobord(1)
846 
847  else !! boundary edge
848  fctab(2, j1) = cl
849 
850  end if
851  end if
852 
853  case(3)
854 
855  !! In this case the edge is internal,
856  !! but one of the three cells in the
857  !! co-boundary is of dimension 2.
858  !!
859  l1 = cell_dim( m%clType(cobord(1)) )
860  if (l1==2) then
861 
862  if (cobord(2)==cl) then
863  fctab(2, j1) = cobord(3)
864  else
865  fctab(2, j1) = cobord(2)
866  end if
867 
868  else
869 
870  l1 = cell_dim( m%clType(cobord(2)) )
871  if (l1==2) then
872 
873  if (cobord(1)==cl) then
874  fctab(2, j1) = cobord(3)
875  else
876  fctab(2, j1) = cobord(1)
877  end if
878 
879  else
880 
881  if (cobord(1)==cl) then
882  fctab(2, j1) = cobord(2)
883  else
884  fctab(2, j1) = cobord(1)
885  end if
886 
887  end if
888  end if
889 
890  end select
891 
892  !! increment the edge
893  j1 = j1 + 1
894  end do
895  end do
896  !$OMP END DO
897  !$OMP END PARALLEL
898 
899  end subroutine cell_loop
900 
901 
902  !! Array of all the 3D-cell faces,
903  !! internal faces are counted twice
904  !!
905  !! Step 2 = global numberingg of all the 3D-cell faces
906  !!
907  subroutine number_faces(cpt)
908  integer, intent(out) :: cpt
909 
910  integer :: jj, cl, cl2, ll
911 
912  cpt = 0
913  do jj=1, m%clToItf%nnz
914 
915  cl = fctab(1,jj)
916  cl2 = fctab(2,jj)
917 
918  if ( cl <= cl2 ) then
919 
920  cpt = cpt + 1
921  fctab(3,jj) = cpt
922 
923  else
924  ll = m%clToItf%row(cl2)
925 
926  do while( fctab(2,ll) /= cl )
927  ll = ll + 1
928  end do
929  fctab(3,jj) = fctab(3,ll)
930 
931  end if
932 
933  end do
934 
935  end subroutine number_faces
936 
937  end subroutine define_3dmesh_faces
938 
939 
940 
941  subroutine define_1dmesh_vertexes(m)
942  type(mesh), intent(inout) :: m
943 
944  integer, dimension(:,:), allocatable :: vtxTab
945  integer, dimension(:) , allocatable :: nnz
946  integer :: cl, j1, j2, cl_t
947 
948  integer :: width, nb_itf
949  real(DP) :: t1, t2
950 
951  if (choral_verb>1) write(*,*) &
952  & "mesh_tools : define_1DMesh_vertexes"
953 
954  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
955  !!
956  !! Computation of the connectivity graph
957  !! cells --> interfaces
958  !!
959  t1 = clock()
960 
961  !! STEP 1: compute the graph pattern
962  !!
963  call allocmem(nnz, m%nbCl)
964  do cl=1, m%nbCl
965  cl_t = m%clType(cl)
966  if ( cell_dim(cl_t) == 1 ) then
967  nnz(cl) = 2
968  else
969  nnz(cl) = 0
970  end if
971  end do
972  m%clToItf= graph(nnz, m%nbCl, 0)
973  call freemem(nnz)
974 
975  !! STEP 2 : create vtxTab
976  !!
977  call allocmem(vtxtab, 3, m%clToItf%nnz)
978  width = m%ndToCl%width
979  call cell_loop()
980 
981 
982  !! STEP 3 : number interfaces
983  !!
984  call number_vertexes(nb_itf)
985  m%nbItf = nb_itf
986  m%clToItf%nc = nb_itf
987 
988  !! fill in the conectivity graph 1D-cells --> vertexes
989  !!
990  do cl=1, m%nbCl
991 
992  j1 = m%clToItf%row(cl)
993  j2 = m%clToItf%row(cl+1)-1
994 
995  m%clToItf%col(j1:j2) = vtxtab(3,j1:j2)
996 
997  end do
998 
999  call freemem(vtxtab)
1000  t2 = clock() - t1
1001  if (choral_verb>1) write(*,*) &
1002  & " 1D-cells --> vertex graph, CPU = ", &
1003  & real(t2, sp)
1004 
1005  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006  !!
1007  !! Computation of the connectivity graph
1008  !! interfaces --> cells
1009  !!
1010  t1 = clock()
1011  call transpose(m%itfToCl, m%clToItf)
1012  t2 = clock() - t1
1013  if (choral_verb>1) write(*,*) &
1014  & " vertex --> 1D-cells graph, CPU = ", &
1015  & real(t2, sp)
1016 
1017  contains
1018 
1019  !! Array of all the 1D-cell vertexes,
1020  !! internal vertexes are counted twice
1021  !!
1022  !! Step 1 = get the co-boundary of all 1D-cell vertexes
1023  !!
1024  subroutine cell_loop()
1025 
1026  integer, dimension(2) :: cl_vtx
1027  integer, dimension(width) :: coBord
1028  integer :: vtx, v, l1, l2
1029 
1030  do cl=1, m%nbCl
1031  l1 = m%clType(cl)
1032  if ( cell_dim(l1) /= 1 ) cycle
1033 
1034  !! cell cl vertexes
1035  !!
1036  j1 = m%clToNd%row(cl)
1037  j2 = j1 + 1
1038  cl_vtx = m%clToNd%col( j1:j2 )
1039 
1040  !! build the array of vertexes
1041  !!
1042  j1 = m%clToItf%row(cl)
1043  do vtx=1, 2
1044 
1045  vtxtab(1, j1) = cl
1046 
1047  v = cl_vtx(vtx)
1048  l1 = m%ndToCl%row(v)
1049  l2 = m%ndToCl%row(v+1)
1050 
1051  j2 = l2-l1
1052  cobord(1:j2) = m%ndToCl%col(l1:l2-1)
1053 
1054  !! j2 = number of cells sharing this vertex
1055  select case(j2)
1056 
1057  case(1) !! boundary vertex
1058  vtxtab(2, j1) = cl
1059 
1060  case(2)
1061  !! In this case the vertex can be
1062  !! either internal or not,
1063 
1064  if (cobord(1)==cl) then
1065  !! it is not necessary to know whether
1066  !! the vertex is internal or not here
1067  !!
1068  vtxtab(2, j1) = cobord(2)
1069 
1070  else
1071  !! it is necessary to know whether
1072  !! the vertex is internal or not here
1073  !!
1074  !! For this we check the dimension
1075  !! of the second cell
1076  l1 = cell_dim( m%clType(cobord(1)) )
1077 
1078  if ( l1==1 ) then !! internal vertex
1079  vtxtab(2, j1) = cobord(1)
1080 
1081  else !! boundary vertex
1082  vtxtab(2, j1) = cl
1083 
1084  end if
1085  end if
1086 
1087  case(3)
1088  !! In this case the vertex is internal,
1089  !! but one of the three cells in the
1090  !! co-boundary is of dimension 1.
1091  !!
1092  l1 = cell_dim( m%clType(cobord(1)) )
1093  if (l1==0) then
1094 
1095  if (cobord(2)==cl) then
1096  vtxtab(2, j1) = cobord(3)
1097  else
1098  vtxtab(2, j1) = cobord(2)
1099  end if
1100 
1101  else
1102 
1103  l1 = cell_dim( m%clType(cobord(2)) )
1104  if (l1==0) then
1105 
1106  if (cobord(1)==cl) then
1107  vtxtab(2, j1) = cobord(3)
1108  else
1109  vtxtab(2, j1) = cobord(1)
1110  end if
1111 
1112  else
1113 
1114  if (cobord(1)==cl) then
1115  vtxtab(2, j1) = cobord(2)
1116  else
1117  vtxtab(2, j1) = cobord(1)
1118  end if
1119 
1120  end if
1121  end if
1122 
1123  case default
1124  call quit("mesh_tools: define_1DMesh_vertexes: non conformal 1D mesh")
1125 
1126  end select
1127 
1128  !! increment the vertex
1129  j1 = j1 + 1
1130  end do
1131  end do
1132 
1133  end subroutine cell_loop
1134 
1135 
1136  !! Array of all the 2D-cell vertexes,
1137  !! internal vertexes are counted twice
1138  !!
1139  !! Step 2 = global numberingg of all the 2D-cell vertexes
1140  !!
1141  subroutine number_vertexes(cpt)
1142  integer, intent(out) :: cpt
1143 
1144  integer :: jj, cl, cl2, ll
1145 
1146  cpt = 0
1147  do jj=1, m%clToItf%nnz
1148 
1149  cl = vtxtab(1,jj)
1150  cl2 = vtxtab(2,jj)
1151 
1152  if ( cl <= cl2 ) then
1153 
1154  cpt = cpt + 1
1155  vtxtab(3,jj) = cpt
1156 
1157  else
1158  ll = m%clToItf%row(cl2)
1159 
1160  do while( vtxtab(2,ll) /= cl )
1161  ll = ll + 1
1162  end do
1163  vtxtab(3,jj) = vtxtab(3,ll)
1164 
1165  end if
1166 
1167  end do
1168 
1169  end subroutine number_vertexes
1170 
1171  end subroutine define_1dmesh_vertexes
1172 
1184  subroutine add_cells(m, type, nodes)
1185  type(mesh), intent(inout) :: m
1186  integer, dimension(:) , intent(in) :: type
1187  integer, dimension(:,:), intent(in) :: nodes
1188 
1189  integer, dimension(:) , allocatable :: clType
1190  integer, dimension(:,:), allocatable :: clTab
1191 
1192  integer :: nn, j1, j2, nNd, cl
1193 
1194  nn = size(type,1)
1195  if (choral_verb>1) write(*,*) &
1196  &"mesh_tools : add_cells =", nn, " new cells"
1197 
1198  !! For the new mesh:
1199  !!
1200  !! cell types will be in clType
1201  !! cell nodes will be in clTab
1202  !!
1203  call allocmem(cltab , cell_max_nbnodes, nn + m%nbCl)
1204  call allocmem(cltype, nn + m%nbCl)
1205 
1206  !! initialise clType and clTab
1207  !! with the cells already present in the mesh
1208  !!
1209  cltype(nn+1:) = m%clType
1210  do cl = 1, m%nbCl
1211 
1212  j1 = m%clToNd%row(cl)
1213  j2 = m%clToNd%row(cl+1)
1214  nnd = j2-j1
1215 
1216  cltab(1:nnd, cl + nn) = m%clToNd%col(j1: j2-1)
1217 
1218  end do
1219 
1220 
1221  !! Add the new boundary cells to
1222  !! clType(1:cpt2) and to clTab(:, 1:cpt2)
1223  !!
1224  cltype(1:nn) = type
1225  cltab(:,1:nn) = nodes
1226 
1227  !! Update m%clType
1228  !!
1229  call allocmem(m%clType, nn + m%nbCl)
1230  m%clType = cltype
1231  call freemem(cltype)
1232 
1233  !! build the new mesh including the boundary cells
1234  !!
1235  call mesh_clear_2(m)
1236  call mesh_create_cltab(m, cltab)
1237  call freemem(cltab)
1238  call mesh_create_end(m)
1239 
1240  if (.NOT.valid(m)) call quit('mesh_interfaces: &
1241  & add_cells: construction not valid')
1242 
1243  end subroutine add_cells
1244 
1245 end module mesh_interfaces
subroutine boundary_interface_scan(tab, m)
BUILDS Tab(1:3,1:mnbItf)
Derived type mesh.
Definition: mesh_mod.F90:92
is the csr matrix well defined ?
Definition: csr_mod.F90:101
DERIVED TYPE graph: sparse matrices of boolean in CSR format
Definition: graph_mod.F90:33
integer, parameter cell_edg_2
Quadratic edge.
integer, parameter cell_tot_nb
Number of CELL types.
integer, parameter cell_tet
Tetrahedron.
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
integer, dimension(cell_tot_nb), public cell_nbvtx
Number of vertexes for each cell type.
Definition: cell_mod.f90:118
double precision function, public clock()
Returns the internal clock time
Definition: basic_tools.F90:94
subroutine, public mesh_clear_2(m)
Destructor for mesh type clears the data tn the mesh built after &#39;call create(...)&#39;.
Definition: mesh_mod.F90:184
integer, dimension(cell_tot_nb), public cell_nbfc
Number of faces for each cell type.
Definition: cell_mod.f90:124
subroutine add_cells(m, type, nodes)
Adds N cells to the mesh These new cells are added in first position.
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public define_interfaces(m)
Defines the mesh interfaces.
integer, parameter, public sp
simple precision for real numbers
Definition: real_type.F90:60
subroutine define_3dmesh_faces(m)
integer, parameter cell_edg
Edge (line segment)
integer, dimension(cell_tot_nb), public cell_dim
Dimension for each cell type.
Definition: cell_mod.f90:130
integer, dimension(cell_max_nbvtx, cell_max_nbfc, cell_tot_nb), public cell_fc_vtx
CELL_FC_VTX(1:n, fc, cl) = vertexes for the face fc of the cell cl, with n = CELL_FC_NBVTX(fc, cl)
Definition: cell_mod.f90:146
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
integer, dimension(cell_max_nbfc, cell_tot_nb), public cell_fc_nbvtx
CELL_FC_NBVTX(fc, cl) = number of vertexes for the face fc of the cell cl.
Definition: cell_mod.f90:141
subroutine define_1dmesh_vertexes(m)
integer, parameter cell_trg
Triangle.
subroutine number_vertexes(cpt)
CHORAL CONSTANTS
subroutine number_edges(cpt)
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
integer, dimension(cell_tot_nb), public cell_nbed
Number of edges for each cell type.
Definition: cell_mod.f90:121
DEFINE THE INTERFACES AND THE BOUNDARY CELLS OF A mesh
DERIVED TYPE mesh
Definition: mesh_mod.F90:57
subroutine, public define_domain_boundary(m)
For a mesh of dimension d The boundary of the mesh domain is a collection of cells of dimension d-1...
allocate memory for real(RP) arrays
Definition: real_type.F90:158
character(len=4), dimension(cell_tot_nb), public cell_name
CELL_ARRAY Arrays describing cells.
Definition: cell_mod.f90:112
integer choral_verb
Verbosity level.
subroutine, public mesh_create_end(m)
mesh constructor : final step
Definition: mesh_mod.F90:567
subroutine, public mesh_create_cltab(m, clTab)
mesh constructor, second step, cell types and cell nodes are known cell types are in mclType cell nod...
Definition: mesh_mod.F90:380
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
subroutine define_2dmesh_edges(m)
subroutine cell_loop()
Definition: diffusion.F90:232
integer, parameter cell_vtx
Vertex.
subroutine, public cap_sorted_set(cpt, t, n, t1, n1, t2, n2)
t(1:cpt) = t1(1:n1) \cap t2(1:n2)
integer, dimension(2, cell_max_nbed, cell_tot_nb), public cell_ed_vtx
CELL_ED_VTX(1:2, ed, cl) = vertexes for the edge ed of the cell cl.
Definition: cell_mod.f90:137
The type graph stores sparse matrices of boolean in CSR format.
Definition: graph_mod.F90:78
DEFINITION OF GEOMETRICAL CELLS (for meshes)
Definition: cell_mod.f90:82
subroutine number_faces(cpt)
integer, parameter, public cell_max_nbnodes
Cell maximal number of nodes.