107 logical,
dimension(:),
allocatable :: b
113 &
'elasticity : matrix_pattern' 116 &
"elasticity: elasticity_matrix_pattern:& 117 & feSpacexk 'Y' not valid" )
119 allocate( b(x_h%mesh%nbCl) )
122 &
"elasticity: elasticity_matrix_pattern:& 123 & quad mesh 'qdmh' not valid" )
125 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
126 &
"elasticity: elasticity_matrix_pattern:& 127 & 'Y' and 'qdmh' associated to different meshes" )
129 do ii=1, x_h%mesh%nbCl
131 & (x_h%feType(ii)==
fe_none) )
then 175 type(
csr) ,
intent(inout):: mass
176 procedure(r3tor) :: a
179 type(
graph) ,
intent(in),
optional :: dofToDof
182 type(
mesh) ,
pointer :: m
184 integer :: dim, nbDof, nn, ft, qt
187 &
"elasticity : massMat" 189 &
" pattern provided = ",&
193 &
"elasticity: elasticity_massMat:& 194 & feSpacexk 'Y' not valid" )
202 &
"elasticity: elasticity_massMat:& 203 & quad mesh 'qdmh' not valid" )
205 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
206 &
"elasticity: elasticity_massMat:& 207 & 'Y' and 'qdmh' associated to different meshes" )
211 if(y%k /= m%dim)
call quit( &
212 &
"elasticity: elasticity_massMat:& 213 & exponent 'Y%k' must be equal to mesh%dim" )
217 if (
present(doftodof))
then 220 &
"elasticity: elasticity_massMat:& 221 & matrix pattern 'dofToDof' not valid " )
231 if (x_h%fe_count(ft)==0) cycle
234 if (qdm%quad_count(qt)==0) cycle
241 if(dim /= m%dim)
call quit( &
242 &
'elasticity: elasticity_massMat:& 243 & integration restricted to cell& 244 & with dimension mesh%dim' )
255 real(RP),
dimension(nbDof, nbDof) :: mat
256 real(RP),
dimension(nbDof, nbDof, nn) :: mat0
257 real(RP),
dimension(nbDof) :: u
258 integer ,
dimension(nbDof) :: p
259 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
260 real(RP),
dimension(3, CELL_MAX_NBNODES) :: X
261 real(RP),
dimension( nn) :: wgt
262 real(RP),
dimension(dim, nn) :: y
265 integer :: cl, ii, jj, ll, cl_nd, j2, l2, component
276 call scal_fe(u, nbdof, y(:,ll), dim, ft)
279 mat0(ii,jj,ll) = u(ii)*u(jj)*wgt(ll)
292 jj = qdm%quadType(cl)
301 x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
305 call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
313 s = a(g%Ty(:,ll))*g%Jy(ll)
317 mat(ii,jj) = mat(ii,jj) + &
326 mat(jj,ii) = mat(ii,jj)
331 call getrow(p, nbdof, x_h%clToDof, cl)
338 j2 = ( p(jj) - 1) * dim + component
339 l2 = ( p(ll) - 1) * dim + component
402 type(
csr) ,
intent(inout):: stiff
403 procedure(r3tor) :: lambda, mu
406 type(
graph) ,
intent(in),
optional :: dofToDof
408 type(
mesh) ,
pointer :: m
411 integer :: dim, nbDof, nn, ft, qt, nbDof2
414 &
"elasticity : stiffMat" 416 &
" pattern provided = ",&
420 &
"elasticity: elasticity_stiffMat:& 421 & feSpacexk 'Y' not valid" )
429 &
"elasticity: elasticity_stiffMat:& 430 & quad mesh 'qdmh' not valid" )
432 if(.NOT.
associated( qdm%mesh, x_h%mesh))
call quit( &
433 &
"elasticity: elasticity_stiffMat:& 434 & 'Y' and 'qdmh' associated to different meshes" )
438 if(y%k /= m%dim)
call quit( &
439 &
"elasticity: elasticity_stiffMat:& 440 & exponent 'Y%k' must be equal to mesh%dim" )
444 if (
present(doftodof))
then 447 &
"elasticity: elasticity_stiffMat:& 448 & matrix pattern 'dofToDof' not valid " )
450 stiff =
csr(doftodof)
458 if (x_h%fe_count(ft)==0) cycle
461 if (qdm%quad_count(qt)==0) cycle
468 if(dim /= m%dim)
call quit( &
469 &
'elasticity: elasticity_stiffMat:& 470 & integration restricted to cell& 471 & with dimension mesh%dim' )
484 real(RP),
dimension(dim, dim, dim, nbDof):: eu, Aeu
485 real(RP),
dimension(dim, nbDof, nn) :: grad_ref
486 real(RP),
dimension(3 , nbDof) :: grad
487 real(RP),
dimension(nbDof2, nbDof2) :: mat
488 real(RP),
dimension(dim, dim) :: prod
489 integer ,
dimension(nbDof2) :: p
490 integer ,
dimension(CELL_MAX_NBNODES) :: p_cl
491 real(RP),
dimension(3, CELL_MAX_NBNODES) :: X
492 real(RP),
dimension(nn) :: sqrt_wgt
495 integer :: cl, ii, jj, ll, cc, i2, c2, cl_nd, i, j
496 real(RP) :: inv_J, lbda_Ty, mu_Tyx2, val
503 if ( minval(sqrt_wgt) < 0.0_rp )
call quit( &
504 &
"elasticity: elasticity_stiffMat:& 505 & quadrature rule with negative weights not allowed" )
507 sqrt_wgt = sqrt(sqrt_wgt)
515 grad_ref(:,:,ll) = grad_ref(:,:,ll) * sqrt_wgt(ll)
526 jj = qdm%quadType(cl)
535 x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
539 call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
552 inv_j = 1.0_rp / g%Jy(ll)
555 lbda_ty = lambda(g%Ty(:,ll))
556 mu_tyx2 = 2._rp * mu(g%Ty(:,ll))
564 & g%Cy(:,:,ll), grad_ref(:,ii,ll))
570 & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
576 & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
589 eu(:,:,1,ii) =
e3(grad(1:3,ii), 1)
590 eu(:,:,2,ii) =
e3(grad(1:3,ii), 2)
591 eu(:,:,3,ii) =
e3(grad(1:3,ii), 3)
596 eu(:,:,1,ii) =
e2(grad(1:2,ii), 1)
597 eu(:,:,2,ii) =
e2(grad(1:2,ii), 2)
602 eu(1,1,1,ii) = grad(1,ii)
611 val = lbda_ty * grad(cc,ii)
614 aeu(jj,jj,cc,ii) = aeu(jj,jj,cc,ii) + val
630 prod = aeu(:,:,cc,ii) * eu(:,:,c2,i2)
631 val = sum( prod ) * inv_j
633 mat(i,j) = mat(i,j) + val
642 call getrow(p, nbdof2, y%clToDof2, cl)
660 function e2(v, component)
661 real(RP),
dimension(2,2) :: e2
662 real(RP),
dimension(2) ,
intent(in) :: v
663 integer ,
intent(in) :: component
667 select case(component)
686 call quit(
"elasticity: e2: 1 <= component <=2" )
694 function e3(v, component)
695 real(RP),
dimension(3,3) :: e3
696 real(RP),
dimension(3) ,
intent(in) :: v
697 integer ,
intent(in) :: component
701 select case(component)
742 call quit(
"elasticity: e3: 1 <= component <=3" )
786 real(RP),
dimension(:),
allocatable :: FV
789 procedure(r3tor) :: f_1, f_2
790 procedure(r3tor) ,
optional :: f_3
793 integer :: i, dim, jj, ll
794 real(RP),
dimension(:),
allocatable :: FV_i
797 &
"elasticity : feSpacexk_L2_product" 800 &
"elasticity: feSpacexk_L2_product: & 801 & fe space 'Y' not valid" )
805 if (y%k /= dim)
call quit(&
806 &
"elasticity: feSpacexk_L2_product: & 807 & fespacsexk exponent 'Y%k' incorrect" )
808 if ( dim==1 )
call quit(&
809 &
"elasticity: feSpacexk_L2_product: & 810 & dim 1 case, use diffusion::L2_product" )
811 if ( (dim==3).AND.(.NOT.
present(f_3)) )
call quit(&
812 &
"elasticity: feSpacexk_L2_product: & 813 & argument 'f_3' missing" )
822 ll = (jj-1) * y%k + i
830 ll = (jj-1) * y%k + i
839 ll = (jj-1) * y%k + i
886 real(RP),
dimension(:),
allocatable :: rhs
888 integer ,
intent(in) :: quad_type
889 procedure(r3tor) :: g_1, g_2
890 procedure(r3tor) ,
optional :: g_3, f
894 integer :: i, dim, jj, ll
895 real(RP),
dimension(:),
allocatable :: rhs_i
898 &
"elasticity : elasticity_Neumann_rhs" 901 &
"elasticity: elasticity_Neumann_rhs: & 902 & fe space 'Y' not valid" )
907 if (y%k /= dim)
call quit(&
908 &
"elasticity: elasticity_Neumann_rhs: & 909 & fespacsexk exponent 'Y%k' incorrect" )
910 if ( dim==1 )
call quit(&
911 &
"elasticity: elasticity_Neumann_rhs: & 912 & dim 1 case, use diffusion::diffusion_Neumann_rhs" )
913 if ( (dim==3).AND.(.NOT.
present(g_3)) )
call quit(&
914 &
"elasticity: elasticity_Neumann_rhs: & 915 & argument 'g_3' missing" )
917 &
"elasticity: elasticity_Neumann_rhs: & 918 & unproper quadrature method" )
921 call set(qdm, quad_type, f)
932 ll = (jj-1) * y%k + i
940 ll = (jj-1) * y%k + i
949 ll = (jj-1) * y%k + i
992 type(
csr) ,
intent(inout) :: K
993 real(RP),
dimension(:),
intent(inout) :: rhs
995 procedure(r3tor) :: g1, g2
996 procedure(r3tor) ,
optional :: g3, f
998 logical ,
dimension(:),
allocatable :: flag
999 real(RP),
dimension(:),
allocatable :: g1_h, g2_h, g3_h
1002 integer :: dim, ii, ll, j1, j2
1005 &
"elasticity : elasticity_Dirichlet" 1008 &
"elasticity: elasticity_Dirichlet_rhs: & 1009 & fe space 'Y' not valid" )
1014 if (y%k /= dim)
call quit(&
1015 &
"elasticity: elasticity_Dirichlet_rhs: & 1016 & fespacsexk exponent 'Y%k' incorrect" )
1017 if ( dim==1 )
call quit(&
1018 &
"elasticity: elasticity_Dirichlet_rhs: & 1019 & dim 1 case, use diffusion::diffusion_Dirichlet_rhs" )
1020 if ( (dim==3).AND.(.NOT.
present(g3)) )
call quit(&
1021 &
"elasticity: elasticity_Dirichlet_rhs: & 1022 & argument 'g3' missing" )
1025 &
"elasticity: elasticity_Dirichlet: & 1026 & csr matrix 'K' not valid" )
1028 if(k%nl/=y%nbDof2)
call quit( &
1029 &
"elasticity: elasticity_Dirichlet: & 1030 & csr matrix 'K' has a wrong size" )
1032 if (
size(rhs,1) /= y%nbDof2 )
call quit( &
1033 &
"elasticity: elasticity_Dirichlet: & 1034 & 'rhs' has a wrong size" )
1051 if ( .NOT. flag(ii) ) cycle
1055 rhs(ll+1) = g1_h(ii)
1056 rhs(ll+2) = g2_h(ii)
1057 if (dim==3) rhs(ll+3) = g3_h(ii)
1061 if (dim==2) j2 = k%row(ll+3)-1
1062 if (dim==3) j2 = k%row(ll+4)-1
DERIVED TYPE geoTsf: geometrical transformation of reference cells
subroutine, public elasticity_massmat(mass, a, Y, qdm, dofToDof)
Assemble the mass matrix of the bilinear product:
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
type(r_1d), dimension(quad_tot_nb), public quad_wgt
quad weights
subroutine fespacexk_l2_product(FV, Y, qdm, f_1, f_2, f_3)
L2 Product of a vector function with the basis functions of a finite element space ...
QUADRATURE RULES ON REFERENCE CELLS
real(rp) function, dimension(2, 2) e2(v, component)
To compute the symmetrised gradient in dim 2.
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
subroutine, public graph_extract(g2, g1, line)
line i of g2 = line i of g1 if line(i) = .TRUE. = void otherwise
integer, parameter quad_none
The type feSpacexk defines for a finite element space.
integer, parameter fe_tot_nb
Number of FE methods.
integer, dimension(quad_tot_nb), public quad_dim
Reference cell dimension for each quad method.
DERIVED TYPE feSpacexk: define for a finite element space.
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)
subroutine, public elasticity_dirichlet(K, rhs, Y, g1, g2, g3, f)
DIRICHLET BOUNDARY CONDITION FOR AN ELASTICITY PROBLEM
DERIVED TYPE feSpace: finite element spaces
type(r_2d), dimension(quad_tot_nb), public quad_coord
quad node coordinates
subroutine, public elasticity_stiffmat(stiff, lambda, mu, Y, qdm, dofToDof)
Assemble the stiffness matrix of the bilinear product:
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 matvecprod_3x1(y, A, X)
matrix vector product R1 –> R3
MODULE FOR LINEAR ELASTICITY PROBLEMS
allocate memory for real(RP) arrays
subroutine, public getdiag(A)
Computes the array Adiag of the indices in Acol for the diagonal coefficients.
L2 product of a function with the basis functions.
Derived type feSpace: finite element space on a mesh.
real(rp) function, dimension(3, 3) e3(v, component)
To compute the symmetrised gradient in dim 3.
integer choral_verb
Verbosity level.
subroutine, public elasticity_neumann_rhs(rhs, Y, quad_type, g_1, g_2, g_3, f)
L2 scalar product of a vector function with the basis functions fo the finite element space on a...
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
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.
USE graph TO DEFINE csr MATRICES
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
integer, parameter fe_none
The type quadMesh defines integration methods on meshes.
DERIVED TYPE csr for sparse matrices
subroutine, public addentry_2(m, rw, cl, val)
Filling matrix m entries.
subroutine, public elasticity_matrix_pattern(g, Y, qdm)
Define the sparsity pattern for elasticity matrices
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, 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.