85 type(
mesh),
pointer :: m
86 integer :: qt, dim, nn
89 &
'integral : integ_scal_func' 92 &
"integral: integ_scal_func:& 93 & quad mesh 'qdmh' not valid" )
99 if (qdm%quad_count(qt)==0) cycle
112 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
113 real(RP),
dimension(3, CELL_MAX_NBNODES):: X
114 real(RP),
dimension(nn) :: wgt, u_x
115 real(RP),
dimension(dim, nn) :: y
118 integer :: cl, jj, ll, cl_nd
133 jj = qdm%quadType(cl)
139 x(:,ll) = m%nd( :, p_cl(ll) )
143 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
148 u_x(ll) = u(g%Ty(:,ll))
152 int2 = int2 + sum(u_x)
194 procedure(r3xrxrtor) :: E
195 procedure(r3tor) :: u
196 real(RP),
dimension(:),
intent(in) :: uh
197 type(
fespace) ,
intent(in) :: X_h
200 type(
mesh),
pointer :: m
201 integer :: dim, nbDof, nn, ft, qt
204 &
'integral : integ_scal_fe' 207 &
"integral: integ_scal_fe:& 208 & fe space 'X_h' not valid" )
211 &
"integral: integ_scal_fe:& 212 & quad mesh 'qdmh' not valid" )
214 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
215 &
"integral: integ_scal_fe:& 216 & 'X_h' and 'qdmh' associated to different meshes" )
218 if(
size(uh,1) /= x_h%nbDof)
call quit( &
219 &
"integral: integ_scal_fe:& 220 & fe function 'uh' has not the expected size" )
227 if (x_h%fe_count(ft)==0) cycle
230 if (qdm%quad_count(qt)==0) cycle
246 real(RP),
dimension(nbDof, nn) :: val
247 integer ,
dimension(nbDof) :: p
248 real(RP),
dimension(nbDof) :: v
249 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
250 real(RP),
dimension(3, CELL_MAX_NBNODES):: X
251 real(RP),
dimension( nn) :: wgt, E_x
252 real(RP),
dimension(dim, nn) :: y
255 integer :: cl, ii, jj, ll, cl_nd
256 real(RP) :: int2, uh_x, u_x
268 call scal_fe(val(:,ll), nbdof, y(:,ll), dim, ft)
279 jj = qdm%quadType(cl)
288 x(:,ll) = m%nd( :, p_cl(ll) )
292 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
296 call getrow(p, nbdof, x_h%clToDof, cl)
305 uh_x = dot_product(val(:,ll), v)
307 e_x(ll)= e(g%Ty(:,ll), u_x, uh_x)
314 int2 = int2 + sum(e_x)
346 procedure(r3xr3xr3tor) :: E
347 real(RP) ,
dimension(:),
intent(in) :: uh
348 type(
fespace) ,
intent(in) :: X_h
350 procedure(r3tor3) :: phi
352 type(
mesh),
pointer :: m
353 integer :: dim, nbDof, nn, ft, qt
356 &
'integral : Integ_scal_fe_grad' 359 &
"integral: integ_scal_fe_grad:& 360 & fe space 'X_h' not valid" )
363 &
"integral: integ_scal_fe_grad:& 364 & quad mesh 'qdmh' not valid" )
366 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
367 &
"integral: integ_scal_fe_grad:& 368 & 'X_h' and 'qdmh' associated to different meshes" )
370 if(
size(uh,1) /= x_h%nbDof)
call quit( &
371 &
"integral: integ_scal_fe_grad:& 372 & fe function 'uh' has not the expected size" )
379 if (x_h%fe_count(ft)==0) cycle
382 if (qdm%quad_count(qt)==0) cycle
397 real(RP),
dimension(dim, nbDof, nn) :: val
398 integer ,
dimension(nbDof) :: p
399 real(RP),
dimension(nbDof) :: v
400 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
401 real(RP),
dimension(3, CELL_MAX_NBNODES):: X
402 real(RP),
dimension(3):: grad_uh_x, w, phi_x
403 real(RP),
dimension( nn) :: wgt
404 real(RP),
dimension(dim, nn) :: y
407 integer :: cl, ii, jj, ll, cl_nd
408 real(RP) :: int2, E_x
421 call scalgrad_fe(val(:,:,ll), nbdof, y(:,ll), dim, ft)
432 jj = qdm%quadType(cl)
441 x(:,ll) = m%nd( :, p_cl(ll) )
445 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
449 call getrow(p, nbdof, x_h%clToDof, cl)
464 & g%Cy(:,:,ll), val(:,jj,ll))
465 grad_uh_x = grad_uh_x + v(jj) * w
471 & g%Cy(:,:,ll), val(:,jj,ll) )
472 grad_uh_x = grad_uh_x + v(jj) * w
478 & g%Cy(:,:,ll), val(:,jj,ll) )
479 grad_uh_x = grad_uh_x + v(jj) * w
483 grad_uh_x = grad_uh_x / g%Jy(ll)
486 phi_x = phi(g%Ty(:,ll))
489 e_x = e(g%Ty(:,ll), phi_x, grad_uh_x)
490 int2 = int2 + e_x * wgt(ll) * g%Jy(ll)
512 procedure(r3xr3xr3tor) :: E
513 real(RP) ,
dimension(:),
intent(in) :: uh
514 type(
fespace) ,
intent(in) :: X_h
516 procedure(r3tor3) :: phi
518 type(
mesh),
pointer :: m
519 integer :: dim, nbDof, nn, ft, qt
522 &
'integral : Integ_scal_fe_grad_proj' 525 &
"integral: integ_scal_fe_grad_proj:& 526 & fe space 'X_h' not valid" )
529 &
"integral: integ_scal_fe_grad_proj:& 530 & quad mesh 'qdmh' not valid" )
532 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
533 &
"integral: integ_scal_fe_grad_proj:& 534 & 'X_h' and 'qdmh' associated to different meshes" )
536 if(
size(uh,1) /= x_h%nbDof)
call quit( &
537 &
"integral: integ_scal_fe_grad_proj:& 538 & fe function 'uh' has not the expected size" )
545 if (x_h%fe_count(ft)==0) cycle
548 if (qdm%quad_count(qt)==0) cycle
564 real(RP),
dimension(dim, nbDof, nn) :: val
565 integer ,
dimension(nbDof) :: p
566 real(RP),
dimension(nbDof) :: v
567 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
568 real(RP),
dimension(3, CELL_MAX_NBNODES):: X
569 real(RP),
dimension(3):: grad_uh_x, w, phi_x
570 real(RP),
dimension( nn) :: wgt
571 real(RP),
dimension(dim, nn) :: y
574 integer :: cl, ii, jj, ll, cl_nd
575 real(RP) :: int2, E_x
589 call scalgrad_fe(val(:,:,ll), nbdof, y(:,ll), dim, ft)
600 jj = qdm%quadType(cl)
609 x(:,ll) = m%nd( :, p_cl(ll) )
613 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
617 call getrow(p, nbdof, x_h%clToDof, cl)
632 & g%Cy(:,:,ll), val(:,jj,ll))
633 grad_uh_x = grad_uh_x + v(jj) * w
639 & g%Cy(:,:,ll), val(:,jj,ll) )
640 grad_uh_x = grad_uh_x + v(jj) * w
646 & g%Cy(:,:,ll), val(:,jj,ll) )
647 grad_uh_x = grad_uh_x + v(jj) * w
651 grad_uh_x = grad_uh_x / g%Jy(ll)
654 phi_x = phi(g%Ty(:,ll))
660 else if (dim==2)
then 663 phi_x = phi_x - sum(phi_x * w)*w
665 else if (dim==1)
then 666 phi_x = sum( g%Cy(:,1,ll)*phi_x ) * g%Cy(:,1,ll)
671 e_x = e(g%Ty(:,ll), phi_x, grad_uh_x)
672 int2 = int2 + e_x * wgt(ll) * g%Jy(ll)
703 procedure(r3xr3xr3tor) :: E
704 real(RP) ,
dimension(:),
intent(in) :: uh1, uh2
705 type(
fespace) ,
intent(in) :: X_h
708 type(
mesh),
pointer :: m
709 integer :: dim, nbDof, nn, ft, qt
712 &
'integral : Integ_scal_fe_grad_2' 715 &
"integral: integ_scal_fe_grad_2:& 716 & fe space 'X_h' not valid" )
719 &
"integral: integ_scal_fe_grad_2:& 720 & quad mesh 'qdmh' not valid" )
722 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
723 &
"integral: integ_scal_fe_grad_2:& 724 & 'X_h' and 'qdmh' associated to different meshes" )
726 if(
size(uh1,1) /= x_h%nbDof)
call quit( &
727 &
"integral: integ_scal_fe_grad_2:& 728 & fe function 'uh1' has not the expected size" )
730 if(
size(uh2,1) /= x_h%nbDof)
call quit( &
731 &
"integral: integ_scal_fe_grad_2:& 732 & fe function 'uh2' has not the expected size" )
740 if (x_h%fe_count(ft)==0) cycle
743 if (qdm%quad_count(qt)==0) cycle
759 real(RP),
dimension(dim, nbDof, nn) :: val
760 integer ,
dimension(nbDof) :: p
761 real(RP),
dimension(nbDof) :: v1, v2
762 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
763 real(RP),
dimension(3, CELL_MAX_NBNODES):: X
764 real(RP),
dimension(3):: grad_uh1_x, grad_uh2_x, w
765 real(RP),
dimension( nn) :: wgt
766 real(RP),
dimension(dim, nn) :: y
769 integer :: cl, ii, jj, ll, cl_nd
770 real(RP) :: int2, E_x
784 call scalgrad_fe(val(:,:,ll), nbdof, y(:,ll), dim, ft)
795 jj = qdm%quadType(cl)
804 x(:,ll) = m%nd( :, p_cl(ll) )
808 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
812 call getrow(p, nbdof, x_h%clToDof, cl)
829 & g%Cy(:,:,ll), val(:,jj,ll))
830 grad_uh1_x = grad_uh1_x + v1(jj) * w
831 grad_uh2_x = grad_uh2_x + v2(jj) * w
837 & g%Cy(:,:,ll), val(:,jj,ll) )
838 grad_uh1_x = grad_uh1_x + v1(jj) * w
839 grad_uh2_x = grad_uh2_x + v2(jj) * w
846 & g%Cy(:,:,ll), val(:,jj,ll) )
847 grad_uh1_x = grad_uh1_x + v1(jj) * w
848 grad_uh2_x = grad_uh2_x + v2(jj) * w
853 grad_uh1_x = grad_uh1_x / g%Jy(ll)
854 grad_uh2_x = grad_uh2_x / g%Jy(ll)
857 e_x = e(g%Ty(:,ll), grad_uh1_x, grad_uh2_x)
858 int2 = int2 + e_x * wgt(ll) * g%Jy(ll)
874 function e_l2(x, u1, u2)
result(r)
876 real(RP),
dimension(3),
intent(in) :: x
877 real(RP),
intent(in) :: u1, u2
885 real(RP),
dimension(3),
intent(in) :: x, p1, p2
886 r = sum( (p1-p2)**2 )
906 real(RP),
dimension(:),
intent(in) :: uh
907 type(
fespace) ,
intent(in) :: X_h
909 procedure(r3tor) :: u
912 &
'integral : L2_dist_scalFE' 935 real(RP),
dimension(:),
intent(in) :: uh
936 type(
fespace) ,
intent(in) :: X_h
938 procedure(r3tor3) :: phi
941 &
'integral : L2_dist_grad_scalFE' 957 real(RP),
dimension(:),
intent(in) :: uh
958 type(
fespace) ,
intent(in) :: X_h
960 procedure(r3tor3) :: phi
963 &
'integral : L2_dist_grad_proj' 985 function l2_dist_vect(phi, phi_h, X_h, qdm)
result(dist)
987 procedure(r3tor3) :: phi
988 real(RP),
dimension(:),
intent(in) :: phi_h
989 type(
fespace) ,
intent(in) :: X_h
992 type(
mesh),
pointer :: m
993 integer :: dim, nbDof, nn, ft, qt
996 &
"integral : L2_dist_vect" 999 &
"integral: L2_dist_vect:& 1000 & fe space 'X_h' not valid" )
1003 &
"integral: L2_dist_vect:& 1004 & quad mesh 'qdmh' not valid" )
1006 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
1007 &
"integral: L2_dist_vect:& 1008 & 'X_h' and 'qdmh' associated to different meshes" )
1010 if(
size(phi_h,1) /= x_h%nbDof)
call quit( &
1011 &
"integral: L2_dist_vect:& 1012 & fe function 'phi_h' has not the expected size" )
1015 if (m%nbItf<=0)
call quit(
"integral: L2_dist_vect:& 1016 & mesh interfaces not defined")
1021 if (x_h%fe_count(ft)==0) cycle
1024 if (qdm%quad_count(qt)==0) cycle
1042 real(RP),
dimension(dim, nbDof, nn) :: base
1043 real(RP),
dimension(3 , nbDof, nn) :: DT_base
1044 real(RP),
dimension(nbDof) :: v
1045 integer ,
dimension(nbDof) :: p
1046 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
1047 real(RP),
dimension(3, CELL_MAX_NBNODES):: X
1048 real(RP),
dimension( nn) :: wgt
1049 real(RP),
dimension(dim, nn) :: y
1050 integer ,
dimension(CELL_MAX_NBITF) :: eps
1051 real(RP),
dimension(3) :: phi_x, phi_h_x
1054 integer :: cl, ii, jj, ll, cl_nd, nb_itf
1063 call vect_fe(base(:,:,ll), nbdof, y(:,ll), dim, ft)
1074 jj = qdm%quadType(cl)
1081 call getrow(p, nbdof, x_h%clToDof, cl)
1092 &
"integral : L2_dist_vect:& 1093 & cell orientation error" )
1103 if (eps(jj)==1) cycle
1109 if (eps(jj)==1) cycle
1115 if (eps(jj)==1) cycle
1117 v(ll:ll+1) = -v(ll:ll+1)
1121 call quit(
"integral: L2_dist_vect:& 1122 & fe type not supported" )
1130 x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
1134 call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
1147 & g%DTy(:,:,ll), base(:,ii,ll) )
1155 & g%DTy(:,:,ll), base(:,ii,ll) )
1163 & g%DTy(:,:,ll), base(:,ii,ll) )
1175 phi_h_x = dt_base(:,1,ll)*v(1)
1177 phi_h_x = phi_h_x + dt_base(:,jj,ll)*v(jj)
1179 s = 1.0_rp / g%Jy(ll)
1180 phi_h_x = phi_h_x * s
1184 phi_x = phi(g%Ty(:,ll))
1185 phi_x = (phi_h_x - phi_x)**2
1186 dist = dist + sum(phi_x) * wgt(ll) * g%Jy(ll)
DERIVED TYPE geoTsf: geometrical transformation of reference cells
integer, parameter fe_rt0_2d
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
real(rp) function integ_scal_func(u, qdm)
Returns .
DERIVED TYPE graph: sparse matrices of boolean in CSR format
real(rp) function integ_scal_fe_grad_2(E, uh1, uh2, X_h, qdm)
Returns .
subroutine, public matvecprod_3x2(y, A, X)
matrix vector product R2 –> R3
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 matvecprod_3x3(y, A, X)
matrix vector product R3 –> R3
DEFINITION OF FINITE ELEMENT METHODS
integer, parameter, public rp
real(kind=RP) = real precision in the code REAL_TOL = epsilon to test real equality ...
real(rp) function, public l2_dist_vect(phi, phi_h, X_h, qdm)
Returns .
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(rp) function l2_dist_scalfe(u, uh, X_h, qdm)
Returns .
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
real(rp) function e_l2(x, u1, u2)
integer, parameter quad_tot_nb
Total number of quadrature rules.
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.
subroutine, public vect_fe(val, nbDof, x, dim, ft)
Evaluate vector FE basis functions at point x.
Integration of scalar functions.
real(rp) function, public l2_dist_grad_proj(phi, uh, X_h, qdm)
Returns .
subroutine, public matvecprod_3x1(y, A, X)
matrix vector product R1 –> R3
integer, parameter, public cell_max_nbitf
Cell maximal number of interfaces.
real(rp) function integ_scal_fe_grad(E, phi, uh, X_h, qdm)
Returns .
real(rp) function integ_scal_fe(E, u, uh, X_h, qdm)
Returns .
real(rp) function l2_dist_grad_scalfe(phi, uh, X_h, qdm)
Returns .
subroutine, public compute_dtjc(g)
Compute DTy(:,:,i), Jy(i) and Cy(:,:,i) at the nodes y(:,i)
Derived type feSpace: finite element space on a mesh.
COMPUTATION OF INTEGRALS using a quadrature methods on a mesh.
integer choral_verb
Verbosity level.
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
real(rp) function e_l2_vect(x, p1, p2)
integer, dimension(fe_tot_nb), public fe_geo
Reference cell geometry for each fe method.
integer, parameter fe_rt1_1d
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
integer, parameter fe_rt1_1d_2
real(rp) function, public integ_scal_fe_grad_proj(E, phi, uh, X_h, qdm)
Returns .
The type quadMesh defines integration methods on meshes.
real(rp) function, dimension(3), public crossprod_3d(u, v)
Cross product in R**3.
logical function, public cell_orientation(m, cl)
has cell 'cl' the direct orientation ?
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.