97 type(
fespace) ,
intent(in) :: X_h
100 logical,
dimension(X_h%mesh%nbCl) :: b
105 &
'diffusion : diffusion_matrix_pattern' 108 &
"diffusion: diffusion_matrix_pattern:& 109 & fe space 'X_h' not valid" )
112 &
"diffusion: diffusion_matrix_pattern:& 113 & quad mesh 'qdmh' not valid" )
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" )
119 do ii=1, x_h%mesh%nbCl
121 & (x_h%feType(ii)==
fe_none) )
then 168 type(
csr) ,
intent(inout):: mass
169 procedure(r3tor) :: a
170 type(
fespace) ,
intent(in) :: X_h
172 type(
graph) ,
intent(in),
optional :: dofToDof
174 type(
mesh),
pointer :: m
176 integer :: dim, nbDof, nn, ft, qt
179 &
"diffusion : massMat" 181 &
" pattern provided = ",&
185 &
"diffusion: diffusion_massMat: & 186 & fe space 'X_h' not valid" )
189 &
"diffusion: diffusion_massMat:& 190 & quad mesh 'qdmh' not valid" )
192 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
193 &
"diffusion: diffusion_massMat:& 194 & 'X_h' and 'qdm' associated with different meshes" )
199 if (
present(doftodof))
then 202 &
"diffusion: diffusion_massMat:& 203 & matrix pattern 'dofToDof' not valid " )
214 if (x_h%fe_count(ft)==0) cycle
217 if (qdm%quad_count(qt)==0) cycle
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
243 integer :: cl, ii, jj, ll, cl_nd
254 call scal_fe(u, nbdof, y(:,ll), dim, ft)
257 mat0(ii,jj,ll) = u(ii)*u(jj)*wgt(ll)
270 jj = qdm%quadType(cl)
279 x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
283 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
291 s = a(g%Ty(:,ll))*g%Jy(ll)
295 mat(ii,jj) = mat(ii,jj) + &
304 mat(jj,ii) = mat(ii,jj)
309 call getrow(p, nbdof, x_h%clToDof, cl)
314 call addentry_2(mass, p(jj), p(ll), mat(jj,ll))
366 type(
csr) ,
intent(inout) :: stiff
367 procedure(r3xr3xr3tor) :: b
368 type(
fespace) ,
intent(in) :: X_h
370 type(
graph) ,
intent(in),
optional :: dofToDof
372 type(
mesh),
pointer :: m
374 integer :: dim, nbDof, nn, ft, qt
377 &
"diffusion : stiffMat" 379 &
" pattern provided = ",&
383 &
"diffusion: diffusion_stiffMat: & 384 & fe space 'X_h' not valid" )
387 &
"diffusion: diffusion_stiffMat:& 388 & quad mesh 'qdmh' not valid" )
390 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
391 &
"diffusion: diffusion_stiffMat:& 392 & 'X_h' and 'qdm' associated with different meshes" )
396 if (
present(doftodof))
then 399 &
"diffusion: diffusion_stiffMat:& 400 & matrix pattern 'dofToDof' not valid " )
402 stiff =
csr(doftodof)
412 if (x_h%fe_count(ft)==0) cycle
415 if (qdm%quad_count(qt)==0) cycle
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
441 integer :: cl, ii, jj, ll, cl_nd
443 real(RP),
dimension(3) :: Ty
450 if ( minval(sqrt_wgt) < 0.0_rp )
call quit( &
451 &
"diffusion: diffusion_stiffMat:& 452 & quadrature rule with negative weights not allowed" )
454 sqrt_wgt = sqrt(sqrt_wgt)
464 grad_ref(:,:,ll) = grad_ref(:,:,ll) * sqrt_wgt(ll)
475 jj = qdm%quadType(cl)
484 x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
488 call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
502 inv_j = 1.0_rp / g%Jy(ll)
510 & g%Cy(:,:,ll), grad_ref(:,ii,ll))
516 & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
522 & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
529 s = b( ty, grad(:,ii), grad(:,jj) )
530 mat(ii,jj) = mat(ii,jj) + s * inv_j
539 mat(jj,ii) = mat(ii,jj)
544 call getrow(p, nbdof, x_h%clToDof, cl)
549 call addentry_2(stiff, p(jj), p(ll), mat(jj,ll))
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
599 &
"diffusion : diffusion_Neumann_rhs" 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" )
609 call set(qdm, quad_type, 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
663 logical ,
dimension(:),
allocatable :: flag
664 real(RP),
dimension(:),
allocatable :: g_h
669 &
"diffusion : diffusion_Dirichlet" 672 &
"diffusion: diffusion_Dirichlet: & 673 & fe space 'X_h' not valid" )
676 &
"diffusion: diffusion_Dirichlet: & 677 & csr matrix 'K' not valid" )
679 if(k%nl/=x_h%nbDof)
call quit( &
680 &
"diffusion: diffusion_Dirichlet: & 681 & csr matrix 'K' has a wrong size" )
683 if (
size(rhs,1) /= x_h%nbDof )
call quit( &
684 &
"diffusion: diffusion_Dirichlet: & 685 & 'rhs' has a wrong size" )
699 if ( .NOT. flag(ii) ) cycle
702 rhs(ii) = rhs(ii) + g_h(ii) * rho
706 k%a(jj) = k%a(jj) + rho
741 real(RP),
dimension(:),
allocatable :: FV
742 type(
fespace) ,
intent(in) :: X_h
744 procedure(r3tor) :: f
746 type(
mesh),
pointer :: m
747 integer :: dim, nbDof, nn, ft, qt
750 &
"diffusion : feSpace_L2_product" 753 &
"diffusion: feSpace_L2_product: & 754 & fe space 'X_h' not valid" )
757 &
"diffusion: feSpace_L2_product:& 758 & quad mesh 'qdmh' not valid" )
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" )
770 if (x_h%fe_count(ft)==0) cycle
773 if (qdm%quad_count(qt)==0) cycle
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
798 integer :: cl, ii, jj, ll, cl_nd
809 call scal_fe(u, nbdof, y(:,ll), dim, ft)
810 val(:,ll) = u * wgt(ll)
821 jj = qdm%quadType(cl)
830 x(1:3,ll) = m%nd( :, p_cl(ll) )
834 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
838 call getrow(p, nbdof, x_h%clToDof, cl)
845 xx = xx + val(:,ll)*fx
899 type(
csr) ,
intent(inout) :: mass
900 procedure(r3xr3xr3tor) :: b
901 type(
fespace) ,
intent(in) :: X_h
903 type(
graph) ,
intent(in),
optional :: dofToDof
905 type(
mesh),
pointer :: m
907 integer :: dim, nbDof, nn, ft, qt
910 &
"diffusion : massMat_vect" 912 &
" pattern provided = ",&
918 &
"diffusion: diffusion_massMat_vect: & 919 & fe space 'X_h' not valid" )
922 &
"diffusion: diffusion_massMat_vect:& 923 & quad mesh 'qdmh' not valid" )
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" )
933 if (
present(doftodof))
then 936 &
"diffusion: diffusion_massMat_vect:& 937 & matrix pattern 'dofToDof' not valid " )
947 if (x_h%fe_count(ft)==0) cycle
950 if (qdm%quad_count(qt)==0) cycle
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
977 integer :: cl, ii, jj, ll, cl_nd, nb_itf
986 if ( minval(wgt) < 0.0_rp )
call quit( &
987 &
"diffusion: diffusion_massMat_vect:& 988 & quadrature rule with negative weights not allowed" )
991 call vect_fe(phi(:,:,ll), nbdof, y(:,ll), dim, ft)
992 phi(:,:,ll) = phi(:,:,ll) * sqrt(wgt(ll))
1003 jj = qdm%quadType(cl)
1012 x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
1016 call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
1029 & g%DTy(:,:,ll), phi(:,ii,ll) )
1037 & g%DTy(:,:,ll), phi(:,ii,ll) )
1045 & g%DTy(:,:,ll), phi(:,ii,ll) )
1063 s = b( g%Ty(:,ll), dt_phi(:,ii,ll), &
1065 mat(ii,jj) = mat(ii,jj) + s * g%Jy(ll)
1074 mat(jj,ii) = mat(ii,jj)
1084 &
"diffusion: diffusion_massMat_vect:& 1085 & cell orientation error" )
1095 if (eps(ii)==1) cycle
1096 mat(ii, :) = -mat(ii, :)
1097 mat( :,ii) = -mat( :,ii)
1102 if (eps(ii)==1) cycle
1103 mat(ii, :) = -mat(ii, :)
1104 mat( :,ii) = -mat( :,ii)
1109 if (eps(ii)==1) cycle
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 )
1118 call quit(
"diffusion: diffusion_massMat_vect:& 1119 & fe type not supported" )
1123 call getrow(p, nbdof, x_h%clToDof, cl)
1128 call addentry_2(mass, p(jj), p(ll), mat(jj,ll))
1164 type(
csr) ,
intent(inout) :: divMat
1165 type(
fespace) ,
intent(in) :: X_s, X_v
1168 type(
graph) :: dof_sToDof_v
1170 type(
mesh),
pointer :: m
1171 integer :: dim, nbDof_s, nbDof_v, nn, ft_s, ft_v, qt
1174 &
"diffusion : mixed_divMat" 1179 &
"diffusion: diffusion_mixed_divMat:& 1180 & scalar fe space 'X_s' not valid" )
1183 &
"diffusion: diffusion_mixed_divMat:& 1184 & vector fe space 'X_v' not valid" )
1187 &
"diffusion: diffusion_mixed_divMat:& 1188 & quad mesh 'qdmh' not valid" )
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" )
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" )
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)
1205 if (x_s%fe_count(ft_s)==0) cycle
1208 if (x_v%fe_count(ft_v)==0) cycle
1211 if (qdm%quad_count(qt)==0) cycle
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
1239 integer :: cl, ii, jj, ll, nb_itf
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)
1265 s = u(ii,ll) * div(jj,ll) * wgt(ll)
1266 mat_0(ii,jj) = mat_0(ii,jj) + s
1277 jj = qdm%quadType(cl)
1296 &
"diffusion: diffusion_mixed_divMat:& 1297 & cell orientation error" )
1307 if (eps(ii)==1) cycle
1309 mat(:, ii) = -mat(:, ii)
1315 if (eps(ii)==1) cycle
1317 mat(:, ii) = -mat(:, ii)
1322 if (eps(ii)==1) cycle
1324 mat( :, 2*ii-1) = -mat( :, 2*ii-1)
1325 mat( :, 2*ii ) = -mat( :, 2*ii )
1329 call quit(
"diffusion: diffusion_mixed_divMat:& 1330 & fe type not supported" )
1334 call getrow(p_s, nbdof_s, x_s%clToDof, cl)
1335 call getrow(p_v, nbdof_v, x_v%clToDof, cl)
1340 call addentry_2(divmat, p_s(ii), p_v(jj), mat(ii,jj))
DERIVED TYPE geoTsf: geometrical transformation of reference cells
integer, parameter fe_rt0_2d
subroutine, public diffusion_matrix_pattern(g, X_h, qdm)
Define the sparsity pattern for diffusion matrices
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
DERIVED TYPE graph: sparse matrices of boolean in CSR format
subroutine, public matvecprod_3x2(y, A, X)
matrix vector product R2 –> R3
L2 product of a function with the basis functions.
subroutine, public compute_dtj(g)
Compute DTy(:,:,i) and Jy(i) at the nodes y(:,i)
type(r_1d), dimension(quad_tot_nb), public quad_wgt
quad weights
QUADRATURE RULES ON REFERENCE CELLS
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
integer, parameter fe_rt1_1d_dual
subroutine, public diffusion_mixed_divmat(divMat, X_s, X_v, qdm)
Assemble the matrix of the bilinear product:
subroutine, public graph_extract(g2, g1, line)
line i of g2 = line i of g1 if line(i) = .TRUE. = void otherwise
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:
integer, parameter fe_rt1_2d_2
integer, parameter fe_tot_nb
Number of FE methods.
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...
integer, dimension(quad_tot_nb), public quad_dim
Reference cell dimension for each quad method.
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
subroutine, public scalgrad_fe(val, nbDof, x, dim, ft)
Evaluate the gradient of the FE basis functions at point x.
Geometrical transformation .
DERIVED TYPE quadMesh: integration methods on meshes
subroutine, public compute_j(g)
Compute Jy(i) at the nodes y(:;i)
DERIVED TYPE feSpace: finite element spaces
type(r_2d), dimension(quad_tot_nb), public quad_coord
quad node coordinates
integer, parameter quad_tot_nb
Total number of quadrature rules.
MODULE FOR DIFFUSION PROBLEMS
integer, dimension(fe_tot_nb), public fe_nbdof
Number of DOF for each fe method.
integer, dimension(quad_tot_nb), public quad_nbnodes
Number of nodes for each quad method.
Derived type for sparse matrices in CSR format.
subroutine, public vect_fe(val, nbDof, x, dim, ft)
Evaluate vector FE basis functions at point x.
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...
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...
subroutine, public diffusion_dirichlet(K, rhs, g, X_h, rho, f)
DIRICHLET BOUNDARY CONDITION FOR A DIFFUSION PROBLEM
integer, parameter, public cell_max_nbitf
Cell maximal number of interfaces.
DEFINE THE INTERFACES AND THE BOUNDARY CELLS OF A mesh
allocate memory for real(RP) arrays
subroutine, public getdiag(A)
Computes the array Adiag of the indices in Acol for the diagonal coefficients.
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)
subroutine, public diffusion_stiffmat(stiff, b, X_h, qdm, dofToDof)
Assemble the stiffness matrix of the bilinear product:
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.
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
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.
subroutine, public diffusion_massmat_vect(mass, b, X_h, qdm, dofToDof)
Assemble the mass matrix of the bilinear product:
DERIVED TYPE csr for sparse matrices
subroutine, public addentry_2(m, rw, cl, val)
Filling matrix m entries.
logical function, public cell_orientation(m, cl)
has cell 'cl' the direct orientation ?
The type graph stores sparse matrices of boolean in CSR format.
subroutine, public prod_tab(g, A, B)
graph product
DEFINITION OF GEOMETRICAL CELLS (for meshes)
integer, parameter fe_rt0_1d
integer, dimension(quad_tot_nb), public quad_geo
Reference cell geometry for each quad method.
subroutine, public scal_fe(val, nbDof, x, dim, ft)
Evaluate FE basis functions at point x.
integer, parameter, public cell_max_nbnodes
Cell maximal number of nodes.