Choral
elasticity.F90
Go to the documentation of this file.
1 
30 
31 module elasticity
32 
35  use real_type
36  use basic_tools
37  use abstract_interfaces, only: r3tor
38  use algebra_lin
39  use cell_mod
40  use graph_mod
41  use mesh_mod, only: mesh, valid
42  use fe_mod
43  use fespace_mod
44  use quad_mod
45  use quadmesh_mod
46  use geotsf_mod
47  use csr_mod
48  use csr_fromgraph
49  use fespacexk_mod
50  use diffusion
51 
52  implicit none
53  private
54 
55  public :: elasticity_matrix_pattern
56  public :: elasticity_massmat !! tested
57  public :: elasticity_stiffmat !! tested
58  public :: l2_product !! tested
59  public :: elasticity_neumann_rhs
60  public :: elasticity_dirichlet
61 
62  ! %----------------------------------------%
63  ! | |
64  ! | GENERIc SUBROUTINES |
65  ! | |
66  ! %----------------------------------------%
67 
70  interface l2_product
71  module procedure fespacexk_l2_product
72  end interface l2_product
73 
74 contains
75 
76 
102  subroutine elasticity_matrix_pattern(g, Y, qdm)
103  type(graph) , intent(inout):: g
104  type(fespacexk), intent(in) :: Y
105  type(quadmesh) , intent(in) :: qdm
106 
107  logical, dimension(:), allocatable :: b
108  type(fespace), pointer :: X_h
109  type(graph) :: g0
110  integer :: ii
111 
112  if (choral_verb>2) write(*,*) &
113  & 'elasticity : matrix_pattern'
114 
115  if(.NOT.valid(y)) call quit( &
116  &"elasticity: elasticity_matrix_pattern:&
117  & feSpacexk 'Y' not valid" )
118  x_h => y%X_h
119  allocate( b(x_h%mesh%nbCl) )
120 
121  if(.NOT.valid(qdm)) call quit( &
122  &"elasticity: elasticity_matrix_pattern:&
123  & quad mesh 'qdmh' not valid" )
124 
125  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
126  &"elasticity: elasticity_matrix_pattern:&
127  & 'Y' and 'qdmh' associated to different meshes" )
128 
129  do ii=1, x_h%mesh%nbCl
130  if ( (qdm%quadType(ii)==quad_none) .OR.&
131  & (x_h%feType(ii)==fe_none) ) then
132  b(ii) = .false.
133  else
134  b(ii) = .true.
135  end if
136  end do
137 
138  call graph_extract(g0, y%ClToDof2, b)
139  deallocate(b)
140 
141  call prod_tab(g, g0, g0)
142 
143  end subroutine elasticity_matrix_pattern
144 
145 
146 
174  subroutine elasticity_massmat(mass, a, Y, qdm, dofToDof)
175  type(csr) , intent(inout):: mass
176  procedure(r3tor) :: a
177  type(fespacexk), intent(in) :: Y
178  type(quadmesh) , intent(in) :: qdm
179  type(graph) , intent(in), optional :: dofToDof
180 
181  type(fespace), pointer :: X_h
182  type(mesh) , pointer :: m
183  type(graph) :: gph
184  integer :: dim, nbDof, nn, ft, qt
185 
186  if (choral_verb>1) write(*,*) &
187  & "elasticity : massMat"
188  if (choral_verb>2) write(*,*) &
189  & " pattern provided = ",&
190  & present(doftodof)
191 
192  if(.NOT.valid(y)) call quit( &
193  &"elasticity: elasticity_massMat:&
194  & feSpacexk 'Y' not valid" )
195 
196  !! shortcuts
197  !!
198  x_h => y%X_h
199  m => x_h%mesh
200 
201  if(.NOT.valid(qdm)) call quit( &
202  &"elasticity: elasticity_massMat:&
203  & quad mesh 'qdmh' not valid" )
204 
205  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
206  &"elasticity: elasticity_massMat:&
207  & 'Y' and 'qdmh' associated to different meshes" )
208 
209  !! Y = [X_h]^k, with k = dimension in space
210  !!
211  if(y%k /= m%dim) call quit( &
212  &"elasticity: elasticity_massMat:&
213  & exponent 'Y%k' must be equal to mesh%dim" )
214 
215  !! connectivity graph
216  !! DOF2 --> DOF2
217  if (present(doftodof)) then
218 
219  if (valid(doftodof)<0) call quit( &
220  & "elasticity: elasticity_massMat:&
221  & matrix pattern 'dofToDof' not valid " )
222 
223  mass = csr(doftodof)
224  else
225  call elasticity_matrix_pattern(gph, y, qdm)
226  mass = csr(gph)
227  call clear(gph)
228  end if
229 
230  do ft=1, fe_tot_nb
231  if (x_h%fe_count(ft)==0) cycle
232 
233  do qt=1, quad_tot_nb
234  if (qdm%quad_count(qt)==0) cycle
235  if (quad_geo(qt) /= fe_geo(ft) ) cycle
236 
237  dim = quad_dim(qt) ! dimension of K_ref
238  nbdof = fe_nbdof(ft) ! fe number of DOF
239  nn = quad_nbnodes(qt) ! quad number of nodes
240 
241  if(dim /= m%dim) call quit( &
242  &'elasticity: elasticity_massMat:&
243  & integration restricted to cell&
244  & with dimension mesh%dim' )
245 
246  call cell_loop()
247 
248  end do
249  end do
250 
251  contains
252 
253  subroutine cell_loop()
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
263 
264  type(geotsf):: g
265  integer :: cl, ii, jj, ll, cl_nd, j2, l2, component
266  real(RP) :: s
267 
268  !! quad nodes and weights
269  !!
270  wgt = quad_wgt(qt)%y
271  y = quad_coord(qt)%y
272 
273  !! u = basis functions evaluated at quad node ll
274  !!
275  do ll=1, nn
276  call scal_fe(u, nbdof, y(:,ll), dim, ft)
277  do ii=1, nbdof
278  do jj=1, nbdof
279  mat0(ii,jj,ll) = u(ii)*u(jj)*wgt(ll)
280  end do
281  end do
282  end do
283 
284  !! geometricat transformation K_ref --> K
285  !!
286  g = geotsf(dim, nn, y)
287 
288  !! CELL LOOP
289  !!
290  do cl=1, m%nbCl
291 
292  jj = qdm%quadType(cl)
293  if (jj/=qt) cycle
294 
295  ii = x_h%feType(cl)
296  if (ii/=ft) cycle
297 
298  ! cell node coordinates
299  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
300  do ll=1, cl_nd
301  x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
302  end do
303 
304  ! transformation T : K_ref --> K = T(K_ref)
305  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
306  call compute_j(g)
307 
308  ! elementary mass matrix on cell cl
309  mat = 0._rp
310 
311  ! loop on quad nodes
312  do ll=1, nn
313  s = a(g%Ty(:,ll))*g%Jy(ll)
314 
315  do ii=1, nbdof
316  do jj=1, ii
317  mat(ii,jj) = mat(ii,jj) + &
318  & s* mat0(ii,jj,ll)
319  end do
320  end do
321  end do
322 
323  ! completion by symmetry
324  do ii=2, nbdof
325  do jj=1, ii-1
326  mat(jj,ii) = mat(ii,jj)
327  end do
328  end do
329 
330  ! dof associated with cell cl indexes
331  call getrow(p, nbdof, x_h%clToDof, cl)
332 
333  ! add elementary mass matrix to global mass matrix
334  do jj=1, nbdof
335  do ll=1, nbdof
336  do component=1, dim
337 
338  j2 = ( p(jj) - 1) * dim + component
339  l2 = ( p(ll) - 1) * dim + component
340  call addentry_2(mass, j2, l2, mat(jj,ll))
341 
342  end do
343  end do
344  end do
345 
346  end do
347 
348  end subroutine cell_loop
349 
350  end subroutine elasticity_massmat
351 
352 
353 
400  subroutine elasticity_stiffmat(stiff, lambda, mu, &
401  & Y, qdm, dofToDof)
402  type(csr) , intent(inout):: stiff
403  procedure(r3tor) :: lambda, mu
404  type(fespacexk), intent(in) :: Y
405  type(quadmesh) , intent(in) :: qdm
406  type(graph) , intent(in), optional :: dofToDof
407 
408  type(mesh) , pointer :: m
409  type(fespace), pointer :: X_h
410  type(graph) :: gph
411  integer :: dim, nbDof, nn, ft, qt, nbDof2
412 
413  if (choral_verb>1) write(*,*) &
414  & "elasticity : stiffMat"
415  if (choral_verb>2) write(*,*) &
416  & " pattern provided = ",&
417  & present(doftodof)
418 
419  if(.NOT.valid(y)) call quit( &
420  &"elasticity: elasticity_stiffMat:&
421  & feSpacexk 'Y' not valid" )
422 
423  !! shortcuts
424  !!
425  x_h => y%X_h
426  m => x_h%mesh
427 
428  if(.NOT.valid(qdm)) call quit( &
429  &"elasticity: elasticity_stiffMat:&
430  & quad mesh 'qdmh' not valid" )
431 
432  if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit( &
433  &"elasticity: elasticity_stiffMat:&
434  & 'Y' and 'qdmh' associated to different meshes" )
435 
436  !! Y = [X_h]^k, with k = dimension in space
437  !!
438  if(y%k /= m%dim) call quit( &
439  &"elasticity: elasticity_stiffMat:&
440  & exponent 'Y%k' must be equal to mesh%dim" )
441 
442  !! connectivity graph
443  !! DOF2 --> DOF2
444  if (present(doftodof)) then
445 
446  if (valid(doftodof)<0) call quit( &
447  & "elasticity: elasticity_stiffMat:&
448  & matrix pattern 'dofToDof' not valid " )
449 
450  stiff = csr(doftodof)
451  else
452  call elasticity_matrix_pattern(gph, y, qdm)
453  stiff = csr(gph)
454  call clear(gph)
455  end if
456 
457  do ft=1, fe_tot_nb
458  if (x_h%fe_count(ft)==0) cycle
459 
460  do qt=1, quad_tot_nb
461  if (qdm%quad_count(qt)==0) cycle
462  if (quad_geo(qt) /= fe_geo(ft) ) cycle
463 
464  dim = quad_dim(qt) ! dimension of K_ref
465  nbdof = fe_nbdof(ft) ! fe number of DOF
466  nn = quad_nbnodes(qt) ! quad number of nodes
467 
468  if(dim /= m%dim) call quit( &
469  &'elasticity: elasticity_stiffMat:&
470  & integration restricted to cell&
471  & with dimension mesh%dim' )
472 
473  nbdof2 = nbdof*dim
474 
475  call cell_loop()
476 
477  end do
478  end do
479 
480  contains
481 
482  subroutine cell_loop()
483 
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
493 
494  type(geotsf):: g
495  integer :: cl, ii, jj, ll, cc, i2, c2, cl_nd, i, j
496  real(RP) :: inv_J, lbda_Ty, mu_Tyx2, val
497 
498  !! square roots of the quad weights
499  !!
500  sqrt_wgt = quad_wgt(qt)%y
501  !!
502  !! negative weights not allowed
503  if ( minval(sqrt_wgt) < 0.0_rp ) call quit( &
504  & "elasticity: elasticity_stiffMat:&
505  & quadrature rule with negative weights not allowed" )
506  !!
507  sqrt_wgt = sqrt(sqrt_wgt)
508 
509  !! gradient of the basis functions
510  !! at the quadrature nodes
511  !!
512  do ll=1, nn
513  call scalgrad_fe(grad_ref(:,:,ll), nbdof, &
514  & quad_coord(qt)%y(:,ll), dim, ft)
515  grad_ref(:,:,ll) = grad_ref(:,:,ll) * sqrt_wgt(ll)
516  end do
517 
518  !! geometricat transformation K_ref --> K
519  !!
520  g = geotsf( dim, nn, quad_coord(qt)%y)
521 
522  !! CELL LOOP
523  !!
524  do cl=1, m%nbCl
525 
526  jj = qdm%quadType(cl)
527  if (jj/=qt) cycle
528 
529  ii = x_h%feType(cl)
530  if (ii/=ft) cycle
531 
532  ! cell node coordinates
533  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
534  do ll=1, cl_nd
535  x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
536  end do
537 
538  ! transformation T : K_ref --> K = T(K_ref)
539  call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
540  call compute_jc(g)
541 
542  !! elementary stiffness matrix on cell cl
543  !!
544  mat = 0._rp
545 
546  !! loop on quad nodes
547  !!
548  !!
549  do ll=1, nn
550 
551  !! Jacobian inverse
552  inv_j = 1.0_rp / g%Jy(ll)
553 
554  !! Lamé coefficients
555  lbda_ty = lambda(g%Ty(:,ll))
556  mu_tyx2 = 2._rp * mu(g%Ty(:,ll))
557 
558  !! grad(:,ii) = Cy(:,:,ll) * grad_ref(:,ii,ll)
559  !!
560  select case(dim)
561  case(3)
562  do ii=1, nbdof
563  call matvecprod_3x3( grad(:,ii), &
564  & g%Cy(:,:,ll), grad_ref(:,ii,ll))
565  end do
566 
567  case(2)
568  do ii=1, nbdof
569  call matvecprod_3x2( grad(:,ii), &
570  & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
571  end do
572 
573  case(1)
574  do ii=1, nbdof
575  call matvecprod_3x1( grad(:,ii), &
576  & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
577  end do
578 
579  end select
580 
581  !! With the assumption Omega \subset \R^d
582  !! we have grda(:,ii) \in \R^d
583  !!
584  !! We now compute the symmetrised gradient
585  !!
586  select case(dim)
587  case(3)
588  do ii=1, nbdof
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)
592  end do
593 
594  case(2)
595  do ii=1, nbdof
596  eu(:,:,1,ii) = e2(grad(1:2,ii), 1)
597  eu(:,:,2,ii) = e2(grad(1:2,ii), 2)
598  end do
599 
600  case(1)
601  do ii=1, nbdof
602  eu(1,1,1,ii) = grad(1,ii)
603  end do
604  end select
605 
606  !! Hooke tensor
607  !!
608  aeu = mu_tyx2 * eu
609  do ii=1, nbdof
610  do cc=1, dim
611  val = lbda_ty * grad(cc,ii)
612 
613  do jj=1, dim
614  aeu(jj,jj,cc,ii) = aeu(jj,jj,cc,ii) + val
615  end do
616 
617  end do
618  end do
619 
620  !! elementary stiffness matrix
621  !!
622  do ii=1, nbdof
623  do cc=1, dim
624  i = (ii-1)*dim + cc
625 
626  do i2=1, nbdof
627  do c2=1, dim
628  j = (i2-1)*dim + c2
629 
630  prod = aeu(:,:,cc,ii) * eu(:,:,c2,i2)
631  val = sum( prod ) * inv_j
632 
633  mat(i,j) = mat(i,j) + val
634  end do
635  end do
636  end do
637  end do
638 
639  end do
640 
641  !! dof associated with cell cl
642  call getrow(p, nbdof2, y%clToDof2, cl)
643 
644  !! add elementary stiffness matrix to global stiffness matrix
645  do i=1, nbdof2
646  do j=1, nbdof2
647  call addentry_2(stiff, p(i), p(j), mat(i,j))
648  end do
649  end do
650 
651  end do
652 
653  end subroutine cell_loop
654 
655  end subroutine elasticity_stiffmat
656 
657 
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
664 
665  real(RP) :: a
666 
667  select case(component)
668  case(1)
669  a = v(2) * 0.5_rp
670 
671  e2(1,1) = v(1)
672  e2(1,2) = a
673  e2(2,1) = a
674  e2(2,2) = 0.0_rp
675 
676  case(2)
677  a = v(1) * 0.5_rp
678 
679  e2(1,1) = 0.0_rp
680  e2(1,2) = a
681  e2(2,1) = a
682  e2(2,2) = v(2)
683 
684 #if DBG>0
685  case default
686  call quit("elasticity: e2: 1 <= component <=2" )
687 #endif
688  end select
689 
690  end function e2
691 
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
698 
699  real(RP) :: a, b
700 
701  select case(component)
702  case(1)
703  a = v(2) * 0.5_rp
704  b = v(3) * 0.5_rp
705 
706  e3(1 ,1) = v(1)
707  e3(2 ,1) = a
708  e3(3 ,1) = b
709  e3(1 ,2) = a
710  e3(2:3,2) = 0.0_rp
711  e3(1 ,3) = b
712  e3(2:3,3) = 0.0_rp
713 
714  case(2)
715  a = v(1) * 0.5_rp
716  b = v(3) * 0.5_rp
717 
718  e3(1,1) = 0.0_rp
719  e3(2,1) = a
720  e3(3,1) = 0.0_rp
721  e3(1,2) = a
722  e3(2,2) = v(2)
723  e3(3,2) = b
724  e3(1,3) = 0.0_rp
725  e3(2,3) = b
726  e3(3,3) = 0.0_rp
727 
728  case(3)
729  a = v(1) * 0.5_rp
730  b = v(2) * 0.5_rp
731 
732  e3(1:2,1) = 0.0_rp
733  e3(3 ,1) = a
734  e3(1:2,2) = 0.0_rp
735  e3(3 ,2) = b
736  e3(1 ,3) = a
737  e3(2 ,3) = b
738  e3(3 ,3) = v(3)
739 
740 #if DBG>0
741  case default
742  call quit("elasticity: e3: 1 <= component <=3" )
743 #endif
744  end select
745 
746  end function e3
747 
748 
749 
785  subroutine fespacexk_l2_product(FV, Y, qdm, f_1, f_2, f_3)
786  real(RP), dimension(:), allocatable :: FV
787  type(fespacexk) , intent(in) :: Y
788  type(quadmesh) , intent(in) :: qdm
789  procedure(r3tor) :: f_1, f_2
790  procedure(r3tor) , optional :: f_3
791 
792  type(fespace), pointer :: X_h
793  integer :: i, dim, jj, ll
794  real(RP), dimension(:), allocatable :: FV_i
795 
796  if (choral_verb>1) write(*,*) &
797  & "elasticity : feSpacexk_L2_product"
798 
799  if(.NOT.valid(y)) call quit( &
800  &"elasticity: feSpacexk_L2_product: &
801  & fe space 'Y' not valid" )
802 
803  x_h => y%X_h
804  dim = x_h%mesh%dim
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" )
814 
815  call allocmem(fv, y%nbDof2)
816  fv = 0._rp
817 
818  !! first component
819  i = 1
820  call l2_product(fv_i, f_1, x_h, qdm)
821  do jj=1, x_h%nbDof
822  ll = (jj-1) * y%k + i
823  fv(ll) = fv_i(jj)
824  end do
825 
826  !! second component
827  i = 2
828  call l2_product(fv_i, f_2, x_h, qdm)
829  do jj=1, x_h%nbDof
830  ll = (jj-1) * y%k + i
831  fv(ll) = fv_i(jj)
832  end do
833  if (dim==2) return
834 
835  !! third component
836  i = 3
837  call l2_product(fv_i, f_3, x_h, qdm)
838  do jj=1, x_h%nbDof
839  ll = (jj-1) * y%k + i
840  fv(ll) = fv_i(jj)
841  end do
842 
843  end subroutine fespacexk_l2_product
844 
845 
885  subroutine elasticity_neumann_rhs(rhs, Y, quad_type, g_1, g_2, g_3, f)
886  real(RP), dimension(:), allocatable :: rhs
887  type(fespacexk) , intent(in) :: Y
888  integer , intent(in) :: quad_type
889  procedure(r3tor) :: g_1, g_2
890  procedure(r3tor) , optional :: g_3, f
891 
892  type(quadmesh) :: qdm
893  type(fespace), pointer :: X_h
894  integer :: i, dim, jj, ll
895  real(RP), dimension(:), allocatable :: rhs_i
896 
897  if (choral_verb>1) write(*,*) &
898  & "elasticity : elasticity_Neumann_rhs"
899 
900  if(.NOT.valid(y)) call quit( &
901  &"elasticity: elasticity_Neumann_rhs: &
902  & fe space 'Y' not valid" )
903 
904  x_h => y%X_h
905  dim = x_h%mesh%dim
906 
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" )
916  if ( quad_dim(quad_type) /= dim-1 ) call quit(&
917  &"elasticity: elasticity_Neumann_rhs: &
918  & unproper quadrature method" )
919 
920  qdm = quadmesh(x_h%mesh)
921  call set(qdm, quad_type, f)
922  call assemble(qdm)
923 
924  call allocmem(rhs, y%nbDof2)
925  rhs = 0._rp
926 
927 
928  !! first component
929  i = 1
930  call l2_product(rhs_i, g_1, x_h, qdm)
931  do jj=1, x_h%nbDof
932  ll = (jj-1) * y%k + i
933  rhs(ll) = rhs_i(jj)
934  end do
935 
936  !! second component
937  i = 2
938  call l2_product(rhs_i, g_2, x_h, qdm)
939  do jj=1, x_h%nbDof
940  ll = (jj-1) * y%k + i
941  rhs(ll) = rhs_i(jj)
942  end do
943  if (dim==2) return
944 
945  !! third component
946  i = 3
947  call l2_product(rhs_i, g_3, x_h, qdm)
948  do jj=1, x_h%nbDof
949  ll = (jj-1) * y%k + i
950  rhs(ll) = rhs_i(jj)
951  end do
952 
953  end subroutine elasticity_neumann_rhs
954 
955 
956 
991  subroutine elasticity_dirichlet(K, rhs, Y, g1, g2, g3, f)
992  type(csr) , intent(inout) :: K
993  real(RP), dimension(:), intent(inout) :: rhs
994  type(fespacexk) , intent(in) :: Y
995  procedure(r3tor) :: g1, g2
996  procedure(r3tor) , optional :: g3, f
997 
998  logical , dimension(:), allocatable :: flag
999  real(RP), dimension(:), allocatable :: g1_h, g2_h, g3_h
1000 
1001  type(fespace), pointer :: X_h
1002  integer :: dim, ii, ll, j1, j2
1003 
1004  if (choral_verb>1) write(*,*) &
1005  & "elasticity : elasticity_Dirichlet"
1006 
1007  if(.NOT.valid(y)) call quit( &
1008  &"elasticity: elasticity_Dirichlet_rhs: &
1009  & fe space 'Y' not valid" )
1010 
1011  x_h => y%X_h
1012  dim = x_h%mesh%dim
1013 
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" )
1023 
1024  if(.NOT.valid(k)) call quit( &
1025  &"elasticity: elasticity_Dirichlet: &
1026  & csr matrix 'K' not valid" )
1027 
1028  if(k%nl/=y%nbDof2) call quit( &
1029  &"elasticity: elasticity_Dirichlet: &
1030  & csr matrix 'K' has a wrong size" )
1031 
1032  if ( size(rhs,1) /= y%nbDof2 ) call quit( &
1033  &"elasticity: elasticity_Dirichlet: &
1034  & 'rhs' has a wrong size" )
1035 
1036  !! Flag the finite element nodes in {\Gamma + f >=0}
1037  !!
1038  call flag_x_h_dof(flag, x_h, x_h%mesh%dim-1, f)
1039 
1040  !! interpolate the boundary source 'g'
1041  !! on all finite element nodes
1042  call interp_scal_func(g1_h, g1, x_h)
1043  call interp_scal_func(g2_h, g2, x_h)
1044  if (dim==3) call interp_scal_func(g3_h, g3, x_h)
1045 
1046  !! modifies rhs(ll) and K(ll,ll)
1047  !! on index ii such that flag(ii) = .TRUE.
1048  call getdiag(k)
1049  do ii=1, x_h%nbDof
1050 
1051  if ( .NOT. flag(ii) ) cycle
1052 
1053  !! modifies rhs(ll)
1054  ll = (ii-1)*dim
1055  rhs(ll+1) = g1_h(ii)
1056  rhs(ll+2) = g2_h(ii)
1057  if (dim==3) rhs(ll+3) = g3_h(ii)
1058 
1059  !! modifies K(ll,:)
1060  j1 = k%row(ll+1)
1061  if (dim==2) j2 = k%row(ll+3)-1
1062  if (dim==3) j2 = k%row(ll+4)-1
1063  k%a(j1:j2) = 0.0_rp
1064 
1065  j1 = k%diag(ll+1)
1066  k%a(j1) = 1.0_rp
1067  j1 = k%diag(ll+2)
1068  k%a(j1) = 1.0_rp
1069  if (dim==3) then
1070  j1 = k%diag(ll+3)
1071  k%a(j1) = 1.0_rp
1072  end if
1073 
1074  end do
1075 
1076  end subroutine elasticity_dirichlet
1077 
1078 
1079 end module elasticity
DERIVED TYPE geoTsf: geometrical transformation of reference cells
Definition: geoTsf_mod.F90:92
Derived type mesh.
Definition: mesh_mod.F90:92
subroutine, public elasticity_massmat(mass, a, Y, qdm, dofToDof)
Assemble the mass matrix of the bilinear product:
Definition: elasticity.F90:175
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
type(r_1d), dimension(quad_tot_nb), public quad_wgt
quad weights
Definition: quad_mod.f90:88
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 ...
Definition: elasticity.F90:786
QUADRATURE RULES ON REFERENCE CELLS
Definition: quad_mod.f90:31
real(rp) function, dimension(2, 2) e2(v, component)
To compute the symmetrised gradient in dim 2.
Definition: elasticity.F90:661
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
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
integer, parameter quad_none
The type feSpacexk defines for a finite element space.
integer, parameter fe_tot_nb
Number of FE methods.
destructor
Definition: real_type.F90:127
integer, dimension(quad_tot_nb), public quad_dim
Reference cell dimension for each quad method.
Definition: quad_mod.f90:79
DERIVED TYPE feSpacexk: define for a finite element space.
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
subroutine, public elasticity_dirichlet(K, rhs, Y, g1, g2, g3, f)
DIRICHLET BOUNDARY CONDITION FOR AN ELASTICITY PROBLEM
Definition: elasticity.F90:992
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
subroutine, public elasticity_stiffmat(stiff, lambda, mu, Y, qdm, dofToDof)
Assemble the stiffness matrix of the bilinear product:
Definition: elasticity.F90:402
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 matvecprod_3x1(y, A, X)
matrix vector product R1 –> R3
MODULE FOR LINEAR ELASTICITY PROBLEMS
Definition: elasticity.F90:31
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
L2 product of a function with the basis functions.
Definition: elasticity.F90:70
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.
Definition: elasticity.F90:695
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...
Definition: elasticity.F90:886
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
Definition: geoTsf_mod.F90:428
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
USE graph TO DEFINE csr MATRICES
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
integer, parameter fe_none
subroutine cell_loop()
Definition: diffusion.F90:232
The type quadMesh defines integration methods on meshes.
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
subroutine, public elasticity_matrix_pattern(g, Y, qdm)
Define the sparsity pattern for elasticity matrices
Definition: elasticity.F90:103
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, 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.