Choral
diffusion.F90
Go to the documentation of this file.
1 
23 
24 module diffusion
25 
28  use real_type
29  use basic_tools
30  use abstract_interfaces, only: r3tor, r3xr3xr3tor
31  use algebra_lin
32  use cell_mod
33  use graph_mod
34  use mesh_mod
35  use mesh_interfaces
36  use fe_mod
37  use fespace_mod
38  use quad_mod
39  use quadmesh_mod
40  use geotsf_mod
41  use csr_mod
42  use csr_fromgraph
43 
44  implicit none
45  private
46 
47  public :: diffusion_matrix_pattern
48  public :: diffusion_massmat ! tested
49  public :: diffusion_stiffmat ! tested
50  public :: diffusion_massmat_vect ! tested
51  public :: diffusion_mixed_divmat ! tested
52  public :: l2_product ! tested
53  public :: diffusion_neumann_rhs ! tested
54  public :: diffusion_dirichlet
55 
56  ! %----------------------------------------%
57  ! | |
58  ! | GENERIc SUBROUTINES |
59  ! | |
60  ! %----------------------------------------%
61 
64  interface l2_product
65  module procedure fespace_l2_product
66  end interface l2_product
67 
68 
69 contains
70 
71 
95  subroutine diffusion_matrix_pattern(g, X_h, qdm)
96  type(graph) , intent(inout):: g
97  type(fespace) , intent(in) :: X_h
98  type(quadmesh), intent(in) :: qdm
99 
100  logical, dimension(X_h%mesh%nbCl) :: b
101  type(graph) :: g0
102  integer :: ii
103 
104  if (choral_verb>2) write(*,*) &
105  & 'diffusion : diffusion_matrix_pattern'
106 
107  if(.NOT.valid(x_h)) call quit( &
108  &"diffusion: diffusion_matrix_pattern:&
109  & fe space 'X_h' not valid" )
110 
111  if(.NOT.valid(qdm)) call quit( &
112  &"diffusion: diffusion_matrix_pattern:&
113  & quad mesh 'qdmh' not valid" )
114 
115  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
116  &"diffusion: diffusion_matrix_pattern:&
117  & 'X_h' and 'qdm' associated with different meshes" )
118 
119  do ii=1, x_h%mesh%nbCl
120  if ( (qdm%quadType(ii)==quad_none) .OR.&
121  & (x_h%feType(ii)==fe_none) ) then
122  b(ii) = .false.
123  else
124  b(ii) = .true.
125  end if
126  end do
127 
128  call graph_extract(g0, x_h%ClToDof, b)
129  call prod_tab(g, g0, g0)
130 
131  end subroutine diffusion_matrix_pattern
132 
133 
134  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135  !!
136  !!
137  !! SCALAR FINITE ELEMENTS
138  !!
139  !!
140  !!
141  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
142 
167  subroutine diffusion_massmat(mass, a, X_h, qdm, dofToDof)
168  type(csr) , intent(inout):: mass
169  procedure(r3tor) :: a
170  type(fespace) , intent(in) :: X_h
171  type(quadmesh), intent(in) :: qdm
172  type(graph) , intent(in), optional :: dofToDof
173 
174  type(mesh), pointer :: m
175  type(graph) :: gph
176  integer :: dim, nbDof, nn, ft, qt
177 
178  if (choral_verb>1) write(*,*) &
179  & "diffusion : massMat"
180  if (choral_verb>2) write(*,*)&
181  & " pattern provided = ",&
182  & present(doftodof)
183 
184  if(.NOT.valid(x_h)) call quit( &
185  &"diffusion: diffusion_massMat: &
186  & fe space 'X_h' not valid" )
187 
188  if(.NOT.valid(qdm)) call quit( &
189  &"diffusion: diffusion_massMat:&
190  & quad mesh 'qdmh' not valid" )
191 
192  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
193  &"diffusion: diffusion_massMat:&
194  & 'X_h' and 'qdm' associated with different meshes" )
195 
196 
197  !! connectivity graph
198  !! DOF --> DOF
199  if (present(doftodof)) then
200 
201  if ( valid(doftodof) < 0 ) call quit( &
202  & "diffusion: diffusion_massMat:&
203  & matrix pattern 'dofToDof' not valid " )
204 
205  mass = csr(doftodof)
206  else
207  call diffusion_matrix_pattern(gph, x_h, qdm)
208  mass = csr(gph)
209  call clear(gph)
210  end if
211  m => x_h%mesh
212 
213  do ft=1, fe_tot_nb
214  if (x_h%fe_count(ft)==0) cycle
215 
216  do qt=1, quad_tot_nb
217  if (qdm%quad_count(qt)==0) cycle
218  if (quad_geo(qt) /= fe_geo(ft) ) cycle
219 
220  dim = quad_dim(qt) ! dimension of K_ref
221  nbdof = fe_nbdof(ft) ! fe number of DOF
222  nn = quad_nbnodes(qt) ! quad number of nodes
223 
224  call cell_loop()
225 
226  end do
227  end do
228 
229  contains
230 
231  subroutine cell_loop()
233  real(RP), dimension(nbDof, nbDof) :: mat
234  real(RP), dimension(nbDof, nbDof, nn):: mat0
235  real(RP), dimension(nbDof) :: u
236  integer , dimension(nbDof) :: p
237  integer , dimension(CELL_MAX_NBNODES) :: p_cl
238  real(RP), dimension(3, CELL_MAX_NBNODES) :: X
239  real(RP), dimension( nn) :: wgt
240  real(RP), dimension(dim, nn) :: y
241 
242  type(geotsf):: g
243  integer :: cl, ii, jj, ll, cl_nd
244  real(RP) :: s
245 
246  !! quad nodes and weights
247  !!
248  wgt = quad_wgt(qt)%y
249  y = quad_coord(qt)%y
250 
251  !! u = basis functions evaluated at quad node ll
252  !!
253  do ll=1, nn
254  call scal_fe(u, nbdof, y(:,ll), dim, ft)
255  do ii=1, nbdof
256  do jj=1, nbdof
257  mat0(ii,jj,ll) = u(ii)*u(jj)*wgt(ll)
258  end do
259  end do
260  end do
261 
262  !! geometricat transformation K_ref --> K
263  !!
264  g = geotsf(dim, nn, y)
265 
266  !! CELL LOOP
267  !!
268  do cl=1, m%nbCl
269 
270  jj = qdm%quadType(cl)
271  if (jj/=qt) cycle
272 
273  ii = x_h%feType(cl)
274  if (ii/=ft) cycle
275 
276  ! cell node coordinates
277  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
278  do ll=1, cl_nd
279  x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
280  end do
281 
282  ! transformation T : K_ref --> K = T(K_ref)
283  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
284  call compute_j(g)
285 
286  ! elementary mass matrix on cell cl
287  mat = 0._rp
288 
289  ! loop on quad nodes
290  do ll=1, nn
291  s = a(g%Ty(:,ll))*g%Jy(ll)
292 
293  do ii=1, nbdof
294  do jj=1, ii
295  mat(ii,jj) = mat(ii,jj) + &
296  & s* mat0(ii,jj,ll)
297  end do
298  end do
299  end do
300 
301  ! completion by symmetry
302  do ii=2, nbdof
303  do jj=1, ii-1
304  mat(jj,ii) = mat(ii,jj)
305  end do
306  end do
307 
308  ! dof associated with cell cl indexes
309  call getrow(p, nbdof, x_h%clToDof, cl)
310 
311  ! add elementary mass matrix to global mass matrix
312  do jj=1, nbdof
313  do ll=1, nbdof
314  call addentry_2(mass, p(jj), p(ll), mat(jj,ll))
315  end do
316  end do
317 
318  end do
319 
320  end subroutine cell_loop
321 
322  end subroutine diffusion_massmat
323 
324 
365  subroutine diffusion_stiffmat(stiff, b, X_h, qdm, dofToDof)
366  type(csr) , intent(inout) :: stiff
367  procedure(r3xr3xr3tor) :: b
368  type(fespace) , intent(in) :: X_h
369  type(quadmesh), intent(in) :: qdm
370  type(graph) , intent(in), optional :: dofToDof
371 
372  type(mesh), pointer :: m
373  type(graph) :: gph
374  integer :: dim, nbDof, nn, ft, qt
375 
376  if (choral_verb>1) write(*,*) &
377  & "diffusion : stiffMat"
378  if (choral_verb>2) write(*,*) &
379  & " pattern provided = ",&
380  & present(doftodof)
381 
382  if(.NOT.valid(x_h)) call quit( &
383  &"diffusion: diffusion_stiffMat: &
384  & fe space 'X_h' not valid" )
385 
386  if(.NOT.valid(qdm)) call quit( &
387  &"diffusion: diffusion_stiffMat:&
388  & quad mesh 'qdmh' not valid" )
389 
390  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
391  &"diffusion: diffusion_stiffMat:&
392  & 'X_h' and 'qdm' associated with different meshes" )
393 
394  !! connectivity graph
395  !! DOF --> DOF
396  if (present(doftodof)) then
397 
398  if ( valid(doftodof) < 0 ) call quit( &
399  & "diffusion: diffusion_stiffMat:&
400  & matrix pattern 'dofToDof' not valid " )
401 
402  stiff = csr(doftodof)
403  else
404  call diffusion_matrix_pattern(gph, x_h, qdm)
405  stiff = csr(gph)
406  call clear(gph)
407  end if
408 
409  m => x_h%mesh
410 
411  do ft=1, fe_tot_nb
412  if (x_h%fe_count(ft)==0) cycle
413 
414  do qt=1, quad_tot_nb
415  if (qdm%quad_count(qt)==0) cycle
416  if (quad_geo(qt) /= fe_geo(ft) ) cycle
417 
418  dim = quad_dim(qt) ! dimension of K_ref
419  nbdof = fe_nbdof(ft) ! fe number of DOF
420  nn = quad_nbnodes(qt) ! quad number of nodes
421 
422  call cell_loop()
423 
424  end do
425  end do
426 
427  contains
428 
429  subroutine cell_loop()
430 
431  real(RP), dimension(dim, nbDof, nn) :: grad_ref
432  real(RP), dimension(3 , nbDof) :: grad
433  real(RP), dimension(nbDof, nbDof) :: mat
434  integer , dimension(nbDof) :: p
435  integer , dimension(CELL_MAX_NBNODES) :: p_cl
436  real(RP), dimension(3, CELL_MAX_NBNODES):: X
437  real(RP), dimension(nn) :: sqrt_wgt
438  real(RP), dimension(dim, nn) :: y
439 
440  type(geotsf):: g
441  integer :: cl, ii, jj, ll, cl_nd
442  real(RP) :: s, inv_J
443  real(RP), dimension(3) :: Ty
444 
445  !! square roots of the quad weights
446  !!
447  sqrt_wgt = quad_wgt(qt)%y
448  !!
449  !! negative weights not allowed
450  if ( minval(sqrt_wgt) < 0.0_rp ) call quit( &
451  & "diffusion: diffusion_stiffMat:&
452  & quadrature rule with negative weights not allowed" )
453  !!
454  sqrt_wgt = sqrt(sqrt_wgt)
455  y = quad_coord(qt)%y
456 
457 
458  !! gradient of the basis functions
459  !! at the quadrature nodes
460  !!
461  do ll=1, nn
462  call scalgrad_fe(grad_ref(:,:,ll), nbdof, &
463  & y(:,ll), dim, ft)
464  grad_ref(:,:,ll) = grad_ref(:,:,ll) * sqrt_wgt(ll)
465  end do
466 
467  !! geometricat transformation K_ref --> K
468  !!
469  g = geotsf(dim, nn, y)
470 
471  !! CELL LOOP
472  !!
473  do cl=1, m%nbCl
474 
475  jj = qdm%quadType(cl)
476  if (jj/=qt) cycle
477 
478  ii = x_h%feType(cl)
479  if (ii/=ft) cycle
480 
481  ! cell node coordinates
482  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
483  do ll=1, cl_nd
484  x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
485  end do
486 
487  ! transformation T : K_ref --> K = T(K_ref)
488  call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
489  call compute_jc(g)
490 
491  !! elementary stiffness matrix on cell cl
492  !!
493  mat = 0._rp
494 
495  !! loop on quad nodes
496  do ll=1, nn
497 
498  !! node coordinates
499  ty = g%Ty(:,ll)
500 
501  !! Jacobian inverse
502  inv_j = 1.0_rp / g%Jy(ll)
503 
504  !! grad(:,ii) = Cy(:,:,ll) * grad_ref(:,ii,ll)
505  !!
506  select case(dim)
507  case(3)
508  do ii=1, nbdof
509  call matvecprod_3x3( grad(:,ii), &
510  & g%Cy(:,:,ll), grad_ref(:,ii,ll))
511  end do
512 
513  case(2)
514  do ii=1, nbdof
515  call matvecprod_3x2( grad(:,ii), &
516  & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
517  end do
518 
519  case(1)
520  do ii=1, nbdof
521  call matvecprod_3x1( grad(:,ii), &
522  & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
523  end do
524 
525  end select
526 
527  do ii=1, nbdof
528  do jj=1, ii
529  s = b( ty, grad(:,ii), grad(:,jj) )
530  mat(ii,jj) = mat(ii,jj) + s * inv_j
531  end do
532  end do
533 
534  end do
535  !!
536  !! completion by symmetry
537  do ii=2, nbdof
538  do jj=1, ii-1
539  mat(jj,ii) = mat(ii,jj)
540  end do
541  end do
542 
543  ! dof associated with cell cl indexes
544  call getrow(p, nbdof, x_h%clToDof, cl)
545 
546  ! add elementary stiffness matrix to global stiffness matrix
547  do jj=1, nbdof
548  do ll=1, nbdof
549  call addentry_2(stiff, p(jj), p(ll), mat(jj,ll))
550  end do
551  end do
552  end do
553 
554  end subroutine cell_loop
555 
556  end subroutine diffusion_stiffmat
557 
558 
559 
589  subroutine diffusion_neumann_rhs(rhs, g, X_h, quad_type, f)
590  real(RP), dimension(:), allocatable :: rhs
591  procedure(r3tor) :: g
592  type(fespace) , intent(in) :: X_h
593  integer , intent(in) :: quad_type
594  procedure(r3tor) , optional :: f
595 
596  type(quadmesh) :: qdm
597 
598  if (choral_verb>1) write(*,*) &
599  & "diffusion : diffusion_Neumann_rhs"
600 
601  if(.NOT.valid(x_h)) call quit( &
602  &"diffusion: diffusion_Neumann_rhs: &
603  & fe space 'X_h' not valid" )
604  if ( quad_dim(quad_type) /= x_h%mesh%dim-1 ) call quit(&
605  &"diffusion: diffusion_Neumann_rhs: &
606  & unproper quadrature method" )
607 
608  qdm = quadmesh(x_h%mesh)
609  call set(qdm, quad_type, f)
610  call assemble(qdm)
611  call l2_product(rhs, g, x_h, qdm)
612 
613  end subroutine diffusion_neumann_rhs
614 
615 
616 
655  subroutine diffusion_dirichlet(K, rhs, g, X_h, rho, f)
656  type(csr) , intent(inout) :: K
657  real(RP), dimension(:), intent(inout) :: rhs
658  procedure(r3tor) :: g
659  type(fespace) , intent(in) :: X_h
660  real(RP) , intent(in) :: rho
661  procedure(r3tor) , optional :: f
662 
663  logical , dimension(:), allocatable :: flag
664  real(RP), dimension(:), allocatable :: g_h
665 
666  integer :: ii, jj
667 
668  if (choral_verb>1) write(*,*) &
669  & "diffusion : diffusion_Dirichlet"
670 
671  if(.NOT.valid(x_h)) call quit( &
672  &"diffusion: diffusion_Dirichlet: &
673  & fe space 'X_h' not valid" )
674 
675  if(.NOT.valid(k)) call quit( &
676  &"diffusion: diffusion_Dirichlet: &
677  & csr matrix 'K' not valid" )
678 
679  if(k%nl/=x_h%nbDof) call quit( &
680  &"diffusion: diffusion_Dirichlet: &
681  & csr matrix 'K' has a wrong size" )
682 
683  if ( size(rhs,1) /= x_h%nbDof ) call quit( &
684  &"diffusion: diffusion_Dirichlet: &
685  & 'rhs' has a wrong size" )
686 
687  call flag_x_h_dof(flag, x_h, x_h%mesh%dim-1, f)
688 
689  !! interpolate the boundary source 'g'
690  !! on all finite element nodes
691  call interp_scal_func(g_h, g, x_h)
692 
693  !! penalise rhs(ii) and K(ii,ii)
694  !! on index ii such that
695  !! such that flag(ii) = .TRUE.
696  call getdiag(k)
697  do ii=1, x_h%nbDof
698 
699  if ( .NOT. flag(ii) ) cycle
700 
701  !! penalise rhs(ii)
702  rhs(ii) = rhs(ii) + g_h(ii) * rho
703 
704  !! penalise K(ii,ii)
705  jj = k%diag(ii)
706  k%a(jj) = k%a(jj) + rho
707 
708  end do
709 
710  end subroutine diffusion_dirichlet
711 
712 
713 
740  subroutine fespace_l2_product(FV, f, X_h, qdm)
741  real(RP), dimension(:), allocatable :: FV
742  type(fespace) , intent(in) :: X_h
743  type(quadmesh) , intent(in) :: qdm
744  procedure(r3tor) :: f
745 
746  type(mesh), pointer :: m
747  integer :: dim, nbDof, nn, ft, qt
748 
749  if (choral_verb>1) write(*,*) &
750  & "diffusion : feSpace_L2_product"
751 
752  if(.NOT.valid(x_h)) call quit( &
753  &"diffusion: feSpace_L2_product: &
754  & fe space 'X_h' not valid" )
755 
756  if(.NOT.valid(qdm)) call quit( &
757  &"diffusion: feSpace_L2_product:&
758  & quad mesh 'qdmh' not valid" )
759 
760  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
761  &"diffusion: feSpace_L2_product:&
762  & 'X_h' and 'qdm' associated with different meshes" )
763 
764  call allocmem(fv, x_h%nbDof)
765  fv = 0._rp
766 
767  m => x_h%mesh
768 
769  do ft=1, fe_tot_nb
770  if (x_h%fe_count(ft)==0) cycle
771 
772  do qt=1, quad_tot_nb
773  if (qdm%quad_count(qt)==0) cycle
774  if (quad_geo(qt) /= fe_geo(ft) ) cycle
775 
776  dim = quad_dim(qt) ! dimension of K_ref
777  nbdof = fe_nbdof(ft) ! fe number of DOF
778  nn = quad_nbnodes(qt) ! quad number of nodes
779 
780  call cell_loop()
781 
782  end do
783  end do
784 
785  contains
786 
787  subroutine cell_loop()
788 
789  real(RP), dimension(nbDof, nn) :: val
790  real(RP), dimension(nbDof) :: xx, u
791  integer , dimension(nbDof) :: p
792  integer , dimension(CELL_MAX_NBNODES) :: p_cl
793  real(RP), dimension(3, CELL_MAX_NBNODES) :: X
794  real(RP), dimension( nn) :: wgt
795  real(RP), dimension(dim, nn) :: y
796 
797  type(geotsf):: g
798  integer :: cl, ii, jj, ll, cl_nd
799  real(RP) :: fx
800 
801  !! quad nodes and weights
802  !!
803  wgt = quad_wgt(qt)%y
804  y = quad_coord(qt)%y
805 
806  !! u = basis functions evalated at quad node ll
807  !!
808  do ll=1, nn
809  call scal_fe(u, nbdof, y(:,ll), dim, ft)
810  val(:,ll) = u * wgt(ll)
811  end do
812 
813  !! geometricat transformation K_ref --> K
814  !!
815  g = geotsf(dim, nn, y)
816 
817  !! CELL LOOP
818  !!
819  do cl=1, m%nbCl
820 
821  jj = qdm%quadType(cl)
822  if (jj/=qt) cycle
823 
824  ii = x_h%feType(cl)
825  if (ii/=ft) cycle
826 
827  ! cell node coordinates
828  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
829  do ll=1, cl_nd
830  x(1:3,ll) = m%nd( :, p_cl(ll) )
831  end do
832 
833  ! transformation T : K_ref --> K = T(K_ref)
834  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
835  call compute_j(g)
836 
837  ! dof associated with cell cl indexes
838  call getrow(p, nbdof, x_h%clToDof, cl)
839 
840  xx = 0._rp
841  ! loop on quadrature nodes
842  do ll=1, nn
843  fx = f(g%Ty(:,ll))
844  fx = fx * g%Jy(ll)
845  xx = xx + val(:,ll)*fx
846  end do
847 
848  fv(p) = fv(p) + xx
849 
850  end do
851 
852  end subroutine cell_loop
853 
854  end subroutine fespace_l2_product
855 
856  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
857  !!
858  !!
859  !! VECTOR FINITE ELEMENTS
860  !!
861  !!
862  !!
863  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
864 
865 
898  subroutine diffusion_massmat_vect(mass, b, X_h, qdm, dofToDof)
899  type(csr) , intent(inout) :: mass
900  procedure(r3xr3xr3tor) :: b
901  type(fespace) , intent(in) :: X_h
902  type(quadmesh), intent(in) :: qdm
903  type(graph) , intent(in), optional :: dofToDof
904 
905  type(mesh), pointer :: m
906  type(graph) :: gph
907  integer :: dim, nbDof, nn, ft, qt
908 
909  if (choral_verb>1) write(*,*) &
910  & "diffusion : massMat_vect"
911  if (choral_verb>2) write(*,*) &
912  & " pattern provided = ",&
913  & present(doftodof)
914 
915  m => x_h%mesh
916 
917  if(.NOT.valid(x_h)) call quit( &
918  &"diffusion: diffusion_massMat_vect: &
919  & fe space 'X_h' not valid" )
920 
921  if(.NOT.valid(qdm)) call quit( &
922  &"diffusion: diffusion_massMat_vect:&
923  & quad mesh 'qdmh' not valid" )
924 
925  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
926  &"diffusion: diffusion_massMat_vect:&
927  & 'X_h' and 'qdm' associated with different meshes" )
928 
929  if (m%nbItf<=0) call define_interfaces(m)
930 
931  !! connectivity graph
932  !! DOF --> DOF
933  if (present(doftodof)) then
934 
935  if ( valid(doftodof) < 0 ) call quit( &
936  & "diffusion: diffusion_massMat_vect:&
937  & matrix pattern 'dofToDof' not valid " )
938 
939  mass = csr(doftodof)
940  else
941  call diffusion_matrix_pattern(gph, x_h, qdm)
942  mass = csr(gph)
943  call clear(gph)
944  end if
945 
946  do ft=1, fe_tot_nb
947  if (x_h%fe_count(ft)==0) cycle
948 
949  do qt=1, quad_tot_nb
950  if (qdm%quad_count(qt)==0) cycle
951  if (quad_geo(qt) /= fe_geo(ft) ) cycle
952 
953  dim = quad_dim(qt) ! dimension of K_ref
954  nbdof = fe_nbdof(ft) ! fe number of DOF
955  nn = quad_nbnodes(qt) ! quad number of nodes
956 
957  call cell_loop()
958 
959  end do
960  end do
961 
962  contains
963 
964  subroutine cell_loop()
965 
966  real(RP), dimension(dim, nbDof, nn) :: phi
967  real(RP), dimension(3 , nbDof, nn) :: DT_phi
968  real(RP), dimension(nbDof, nbDof) :: mat
969  integer , dimension(nbDof) :: p
970  integer , dimension(CELL_MAX_NBNODES) :: p_cl
971  real(RP), dimension(3, CELL_MAX_NBNODES):: X
972  real(RP), dimension( nn) :: wgt
973  real(RP), dimension(dim, nn) :: y
974  integer , dimension(CELL_MAX_NBITF) :: eps
975 
976  type(geotsf):: g
977  integer :: cl, ii, jj, ll, cl_nd, nb_itf
978  real(RP) :: s
979 
980  !! quad nodes and weights
981  !!
982  wgt = quad_wgt(qt)%y
983  y = quad_coord(qt)%y
984 
985  !! negative weights not allowed
986  if ( minval(wgt) < 0.0_rp ) call quit( &
987  & "diffusion: diffusion_massMat_vect:&
988  & quadrature rule with negative weights not allowed" )
989 
990  do ll=1, nn
991  call vect_fe(phi(:,:,ll), nbdof, y(:,ll), dim, ft)
992  phi(:,:,ll) = phi(:,:,ll) * sqrt(wgt(ll))
993  end do
994 
995  !! geometricat transformation K_ref --> K
996  !!
997  g = geotsf(dim, nn, y)
998 
999  !! CELL LOOP
1000  !!
1001  do cl=1, m%nbCl
1002 
1003  jj = qdm%quadType(cl)
1004  if (jj/=qt) cycle
1005 
1006  ii = x_h%feType(cl)
1007  if (ii/=ft) cycle
1008 
1009  ! cell node coordinates
1010  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
1011  do ll=1, cl_nd
1012  x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
1013  end do
1014 
1015  ! transformation T : K_ref --> K = T(K_ref)
1016  call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
1017  call compute_dtj(g)
1018 
1019  !! DT_phi(:,ii) \in \R^3 =
1020  !! DT(:,:,ll).phi(:,ii,ll)
1021  !! with DT(:,:,ll) = DT at node y(:,ll)
1022  !!
1023  select case(dim)
1024 
1025  case(3)
1026  do ll=1, nn
1027  do ii=1, nbdof
1028  call matvecprod_3x3( dt_phi(:,ii,ll), &
1029  & g%DTy(:,:,ll), phi(:,ii,ll) )
1030  end do
1031  end do
1032 
1033  case(2)
1034  do ll=1, nn
1035  do ii=1, nbdof
1036  call matvecprod_3x2( dt_phi(:,ii,ll), &
1037  & g%DTy(:,:,ll), phi(:,ii,ll) )
1038  end do
1039  end do
1040 
1041  case(1)
1042  do ll=1, nn
1043  do ii=1, nbdof
1044  call matvecprod_3x1( dt_phi(:,ii,ll), &
1045  & g%DTy(:,:,ll), phi(:,ii,ll) )
1046  end do
1047  end do
1048 
1049  end select
1050 
1051  !! elementary mass matrix on cell cl
1052  !!
1053  mat = 0._rp
1054  !!
1055  !! Jacobian inverse
1056  g%Jy = 1._rp / g%Jy
1057  !!
1058  !! loop on quad nodes
1059  do ll=1, nn
1060 
1061  do ii=1, nbdof
1062  do jj=1, ii
1063  s = b( g%Ty(:,ll), dt_phi(:,ii,ll), &
1064  & dt_phi(:,jj,ll) )
1065  mat(ii,jj) = mat(ii,jj) + s * g%Jy(ll)
1066  end do
1067  end do
1068 
1069  end do
1070  !!
1071  !! completion by symmetry
1072  do ii=2, nbdof
1073  do jj=1, ii-1
1074  mat(jj,ii) = mat(ii,jj)
1075  end do
1076  end do
1077 
1078  !! Modification of the elementary mass matrix
1079  !! to take into account the interface orientation
1080 #if DBG>0
1081  !!
1082  !! checks the cell orientation
1083  if (.NOT.cell_orientation(m, cl)) call quit( &
1084  & "diffusion: diffusion_massMat_vect:&
1085  & cell orientation error" )
1086 #endif
1087  !!
1088  !! orientation of the interfaces of cell cl
1089  call interface_orientation(nb_itf, &
1090  & eps, cell_max_nbitf, m, cl)
1091  !!
1092  select case(ft)
1093  case(fe_rt0_1d, fe_rt0_2d)
1094  do ii = 1, nbdof
1095  if (eps(ii)==1) cycle
1096  mat(ii, :) = -mat(ii, :)
1097  mat( :,ii) = -mat( :,ii)
1098  end do
1099 
1100  case(fe_rt1_1d)
1101  do ii = 1, 2
1102  if (eps(ii)==1) cycle
1103  mat(ii, :) = -mat(ii, :)
1104  mat( :,ii) = -mat( :,ii)
1105  end do
1106 
1107  case(fe_rt1_2d_2)
1108  do ii = 1, 3
1109  if (eps(ii)==1) cycle
1110 
1111  mat( 2*ii-1, :) = -mat( 2*ii-1, :)
1112  mat( :, 2*ii-1) = -mat( :, 2*ii-1)
1113  mat( 2*ii , :) = -mat( 2*ii , :)
1114  mat( :, 2*ii ) = -mat( :, 2*ii )
1115  end do
1116 
1117  case default
1118  call quit("diffusion: diffusion_massMat_vect:&
1119  & fe type not supported" )
1120  end select
1121 
1122  ! dof associated with cell cl indexes
1123  call getrow(p, nbdof, x_h%clToDof, cl)
1124 
1125  ! add elementary mass matrix to global mass matrix
1126  do jj=1, nbdof
1127  do ll=1, nbdof
1128  call addentry_2(mass, p(jj), p(ll), mat(jj,ll))
1129  end do
1130  end do
1131 
1132  end do
1133 
1134  end subroutine cell_loop
1135 
1136  end subroutine diffusion_massmat_vect
1137 
1138 
1163  subroutine diffusion_mixed_divmat(divMat, X_s, X_v, qdm)
1164  type(csr) , intent(inout) :: divMat
1165  type(fespace) , intent(in) :: X_s, X_v
1166  type(quadmesh), intent(in) :: qdm
1167 
1168  type(graph) :: dof_sToDof_v
1169 
1170  type(mesh), pointer :: m
1171  integer :: dim, nbDof_s, nbDof_v, nn, ft_s, ft_v, qt
1172 
1173  if (choral_verb>1) write(*,*) &
1174  & "diffusion : mixed_divMat"
1175 
1176  m => x_s%mesh
1177 
1178  if(.NOT.valid(x_s)) call quit( &
1179  &"diffusion: diffusion_mixed_divMat:&
1180  & scalar fe space 'X_s' not valid" )
1181 
1182  if(.NOT.valid(x_v)) call quit( &
1183  &"diffusion: diffusion_mixed_divMat:&
1184  & vector fe space 'X_v' not valid" )
1185 
1186  if(.NOT.valid(qdm)) call quit( &
1187  &"diffusion: diffusion_mixed_divMat:&
1188  & quad mesh 'qdmh' not valid" )
1189 
1190  if(.NOT.associated( qdm%mesh, x_s%mesh)) call quit( &
1191  &"diffusion: diffusion_mixed_divMat:&
1192  & 'X_s' and 'qdm' associated with different meshes" )
1193 
1194  if(.NOT.associated( qdm%mesh, x_v%mesh)) call quit( &
1195  &"diffusion: diffusion_mixed_divMat:&
1196  & 'X_v' and 'qdm' associated with different meshes" )
1197 
1198  !! connectivity graph
1199  !! scalar DOF --> vector DOF
1200  call prod_tab(dof_stodof_v, x_s%clToDof, x_v%clToDof)
1201  divmat = csr(dof_stodof_v)
1202  call clear(dof_stodof_v)
1203 
1204  do ft_s=1, fe_tot_nb
1205  if (x_s%fe_count(ft_s)==0) cycle
1206 
1207  do ft_v=1, fe_tot_nb
1208  if (x_v%fe_count(ft_v)==0) cycle
1209 
1210  do qt=1, quad_tot_nb
1211  if (qdm%quad_count(qt)==0) cycle
1212  if (quad_geo(qt) /= fe_geo(ft_s) ) cycle
1213  if (quad_geo(qt) /= fe_geo(ft_v) ) cycle
1214 
1215  dim = quad_dim(qt) ! dimension of K_ref
1216  nbdof_s = fe_nbdof(ft_s) ! scalar fe number of DOF
1217  nbdof_v = fe_nbdof(ft_v) ! vector fe number of DOF
1218  nn = quad_nbnodes(qt) ! quad number of nodes
1219 
1220  call cell_loop()
1221 
1222  end do
1223  end do
1224  end do
1225 
1226  contains
1227 
1228  subroutine cell_loop()
1229 
1230  real(RP), dimension(nbDof_v, nn) :: div
1231  real(RP), dimension(nbDof_s, nn) :: u
1232  real(RP), dimension(nbDof_s, nbDof_v) :: mat_0, mat
1233  integer , dimension(nbDof_s) :: p_s
1234  integer , dimension(nbDof_v) :: p_v
1235  real(RP), dimension( nn) :: wgt
1236  real(RP), dimension(dim, nn) :: y
1237  integer , dimension(CELL_MAX_NBITF) :: eps
1238 
1239  integer :: cl, ii, jj, ll, nb_itf
1240  real(RP) :: s
1241 
1242  !! quad nodes and weights
1243  !!
1244  wgt = quad_wgt(qt)%y
1245  y = quad_coord(qt)%y
1246 
1247  !! Scalar basis functions and
1248  !! divergence of the vector basis functions
1249  !! evaluated at the quad nodes
1250  !!
1251  do ll=1, nn
1252  call scal_fe(u(:,ll), nbdof_s, y(:,ll), dim, ft_s)
1253  call vectdiv_fe(div(:,ll), nbdof_v, y(:,ll), dim, ft_v)
1254  end do
1255 
1256  !! elementary matrix
1257  !!
1258  mat_0 = 0._rp
1259  !!
1260  !! loop on quad nodes
1261  do ll=1, nn
1262 
1263  do ii=1, nbdof_s
1264  do jj=1, nbdof_v
1265  s = u(ii,ll) * div(jj,ll) * wgt(ll)
1266  mat_0(ii,jj) = mat_0(ii,jj) + s
1267  end do
1268  end do
1269 
1270  end do
1271 
1272 
1273  !! CELL LOOP
1274  !!
1275  do cl=1, m%nbCl
1276 
1277  jj = qdm%quadType(cl)
1278  if (jj/=qt) cycle
1279 
1280  ii = x_s%feType(cl)
1281  if (ii/=ft_s) cycle
1282 
1283  ii = x_v%feType(cl)
1284  if (ii/=ft_v) cycle
1285 
1286  !! elementary matrix on cell cl
1287  !!
1288  mat = mat_0
1289  !!
1290  !! Modification of the elementary matrix
1291  !! to take into account the interface orientation
1292 #if DBG>0
1293  !!
1294  !! checks the cell orientation
1295  if (.NOT.cell_orientation(m, cl)) call quit( &
1296  &"diffusion: diffusion_mixed_divMat:&
1297  & cell orientation error" )
1298 #endif
1299  !!
1300  !! orientation of the interfaces of cell cl
1301  call interface_orientation(nb_itf, &
1302  & eps, cell_max_nbitf, m, cl)
1303  !!
1304  select case(ft_v)
1306  do ii = 1, nbdof_v
1307  if (eps(ii)==1) cycle
1308 
1309  mat(:, ii) = -mat(:, ii)
1310  end do
1311 
1312  case(fe_rt1_1d, fe_rt1_1d_dual, &
1314  do ii = 1, 2
1315  if (eps(ii)==1) cycle
1316 
1317  mat(:, ii) = -mat(:, ii)
1318  end do
1319 
1320  case(fe_rt1_2d_2)
1321  do ii = 1, 3
1322  if (eps(ii)==1) cycle
1323 
1324  mat( :, 2*ii-1) = -mat( :, 2*ii-1)
1325  mat( :, 2*ii ) = -mat( :, 2*ii )
1326  end do
1327 
1328  case default
1329  call quit("diffusion: diffusion_mixed_divMat:&
1330  & fe type not supported" )
1331  end select
1332 
1333  ! dof associated with cell cl indexes
1334  call getrow(p_s, nbdof_s, x_s%clToDof, cl)
1335  call getrow(p_v, nbdof_v, x_v%clToDof, cl)
1336 
1337  ! add elementary matrix to global divMat matrix
1338  do ii=1, nbdof_s
1339  do jj=1, nbdof_v
1340  call addentry_2(divmat, p_s(ii), p_v(jj), mat(ii,jj))
1341  end do
1342  end do
1343 
1344  end do
1345 
1346  end subroutine cell_loop
1347 
1348  end subroutine diffusion_mixed_divmat
1349 
1350 
1351 end module diffusion
DERIVED TYPE geoTsf: geometrical transformation of reference cells
Definition: geoTsf_mod.F90:92
integer, parameter fe_rt0_2d
Derived type mesh.
Definition: mesh_mod.F90:92
subroutine, public diffusion_matrix_pattern(g, X_h, qdm)
Define the sparsity pattern for diffusion matrices
Definition: diffusion.F90:96
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
DERIVED TYPE graph: sparse matrices of boolean in CSR format
Definition: graph_mod.F90:33
LINEAR ALGEBRA TOOLS
Definition: algebra_lin.F90:9
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine, public matvecprod_3x2(y, A, X)
matrix vector product R2 –> R3
L2 product of a function with the basis functions.
Definition: diffusion.F90:64
subroutine, public compute_dtj(g)
Compute DTy(:,:,i) and Jy(i) at the nodes y(:,i)
Definition: geoTsf_mod.F90:539
type(r_1d), dimension(quad_tot_nb), public quad_wgt
quad weights
Definition: quad_mod.f90:88
QUADRATURE RULES ON REFERENCE CELLS
Definition: quad_mod.f90:31
subroutine, public flag_x_h_dof(dof_flag, X_h, dim, f)
OUTPUT: flag(:) array of logical with size X_hnbDof.
subroutine, public matvecprod_3x3(y, A, X)
matrix vector product R3 –> R3
DEFINITION OF FINITE ELEMENT METHODS
Definition: fe_mod.F90:19
integer, parameter fe_rt1_1d_dual
subroutine, public diffusion_mixed_divmat(divMat, X_s, X_v, qdm)
Assemble the matrix of the bilinear product:
Definition: diffusion.F90:1164
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public graph_extract(g2, g1, line)
line i of g2 = line i of g1 if line(i) = .TRUE. = void otherwise
Definition: graph_mod.F90:598
subroutine, public define_interfaces(m)
Defines the mesh interfaces.
integer, parameter quad_none
integer, parameter fe_rt0_1d_dual
subroutine, public diffusion_massmat(mass, a, X_h, qdm, dofToDof)
Assemble the mass matrix of the bilinear product:
Definition: diffusion.F90:168
integer, parameter fe_rt1_2d_2
integer, parameter fe_tot_nb
Number of FE methods.
destructor
Definition: real_type.F90:127
subroutine, public interface_orientation(nb_itf, eps, n, m, cl)
eps(ii) = +1 if the interface ii of the cell cl in mesh m is oriented outwards the cell cl = -1 other...
Definition: mesh_mod.F90:824
integer, dimension(quad_tot_nb), public quad_dim
Reference cell dimension for each quad method.
Definition: quad_mod.f90:79
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
subroutine, public scalgrad_fe(val, nbDof, x, dim, ft)
Evaluate the gradient of the FE basis functions at point x.
Definition: fe_mod.F90:166
Geometrical transformation .
Definition: geoTsf_mod.F90:130
DERIVED TYPE quadMesh: integration methods on meshes
subroutine, public compute_j(g)
Compute Jy(i) at the nodes y(:;i)
Definition: geoTsf_mod.F90:383
DERIVED TYPE feSpace: finite element spaces
Definition: feSpace_mod.F90:58
type(r_2d), dimension(quad_tot_nb), public quad_coord
quad node coordinates
Definition: quad_mod.f90:85
integer, parameter quad_tot_nb
Total number of quadrature rules.
MODULE FOR DIFFUSION PROBLEMS
Definition: diffusion.F90:24
integer, dimension(fe_tot_nb), public fe_nbdof
Number of DOF for each fe method.
Definition: fe_mod.F90:68
integer, dimension(quad_tot_nb), public quad_nbnodes
Number of nodes for each quad method.
Definition: quad_mod.f90:73
CHORAL CONSTANTS
Derived type for sparse matrices in CSR format.
Definition: csr_mod.F90:62
subroutine, public vect_fe(val, nbDof, x, dim, ft)
Evaluate vector FE basis functions at point x.
Definition: fe_mod.F90:195
subroutine, public matvecprod_3x1(y, A, X)
matrix vector product R1 –> R3
subroutine, public diffusion_neumann_rhs(rhs, g, X_h, quad_type, f)
L2 scalar product of a scalar function with the basis functions fo the finite element space on a...
Definition: diffusion.F90:590
subroutine fespace_l2_product(FV, f, X_h, qdm)
L2 Product of a scalar function with the basis functions of a scalar finite element space...
Definition: diffusion.F90:741
subroutine, public diffusion_dirichlet(K, rhs, g, X_h, rho, f)
DIRICHLET BOUNDARY CONDITION FOR A DIFFUSION PROBLEM
Definition: diffusion.F90:656
integer, parameter, public cell_max_nbitf
Cell maximal number of interfaces.
DEFINE THE INTERFACES AND THE BOUNDARY CELLS OF A mesh
DERIVED TYPE mesh
Definition: mesh_mod.F90:57
allocate memory for real(RP) arrays
Definition: real_type.F90:158
subroutine, public getdiag(A)
Computes the array Adiag of the indices in Acol for the diagonal coefficients.
Definition: csr_mod.F90:835
Derived type feSpace: finite element space on a mesh.
integer choral_verb
Verbosity level.
integer, parameter fe_rt1_1d_2_dual
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
Definition: geoTsf_mod.F90:428
subroutine, public diffusion_stiffmat(stiff, b, X_h, qdm, dofToDof)
Assemble the stiffness matrix of the bilinear product:
Definition: diffusion.F90:366
subroutine, public interp_scal_func(uh, u, X_h)
Interpolation of a scalar function to a scalar finite element function.
integer, dimension(fe_tot_nb), public fe_geo
Reference cell geometry for each fe method.
Definition: fe_mod.F90:71
integer, parameter fe_rt1_1d
USE graph TO DEFINE csr MATRICES
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
integer, parameter fe_none
integer, parameter fe_rt1_1d_2
subroutine cell_loop()
Definition: diffusion.F90:232
The type quadMesh defines integration methods on meshes.
subroutine, public vectdiv_fe(val, nbDof, x, dim, ft)
Evaluate the divergence of the vector FE basis functions at point x.
Definition: fe_mod.F90:225
subroutine, public diffusion_massmat_vect(mass, b, X_h, qdm, dofToDof)
Assemble the mass matrix of the bilinear product:
Definition: diffusion.F90:899
DERIVED TYPE csr for sparse matrices
Definition: csr_mod.F90:13
subroutine, public addentry_2(m, rw, cl, val)
Filling matrix m entries.
Definition: csr_mod.F90:418
logical function, public cell_orientation(m, cl)
has cell &#39;cl&#39; the direct orientation ?
Definition: mesh_mod.F90:777
The type graph stores sparse matrices of boolean in CSR format.
Definition: graph_mod.F90:78
subroutine, public prod_tab(g, A, B)
graph product
Definition: graph_mod.F90:774
DEFINITION OF GEOMETRICAL CELLS (for meshes)
Definition: cell_mod.f90:82
integer, parameter fe_rt0_1d
integer, dimension(quad_tot_nb), public quad_geo
Reference cell geometry for each quad method.
Definition: quad_mod.f90:76
subroutine, public scal_fe(val, nbDof, x, dim, ft)
Evaluate FE basis functions at point x.
Definition: fe_mod.F90:136
extract a graph row
Definition: graph_mod.F90:149
integer, parameter, public cell_max_nbnodes
Cell maximal number of nodes.