Choral
integral.F90
Go to the documentation of this file.
1 
20 
21 module integral
22 
25  use real_type, only: rp
26  use basic_tools
27  use abstract_interfaces, only: r3xrxrtor, r3tor, r3xr3xr3tor, r3tor3
28  use algebra_lin
29  use cell_mod
30  use fe_mod
31  use graph_mod
32  use mesh_mod
33  use fespace_mod
34  use quad_mod
35  use quadmesh_mod
36  use geotsf_mod
37 
38  implicit none
39  private
40 
41  public :: integ ! partially tested
42  public :: integ_scal_fe_grad_proj
44  public :: l2_dist_vect
45 
46  ! %----------------------------------------%
47  ! | |
48  ! | GENERIc SUBROUTINES |
49  ! | |
50  ! %----------------------------------------%
51 
54  interface integ
55  module procedure integ_scal_func ! tested
56  module procedure integ_scal_fe
57  module procedure integ_scal_fe_grad ! tested
58  module procedure integ_scal_fe_grad_2
59  end interface integ
60 
63  interface l2_dist
64  module procedure l2_dist_scalfe
65  end interface l2_dist
66  interface l2_dist_grad
67  module procedure l2_dist_grad_scalfe
68  end interface l2_dist_grad
69 
70 contains
71 
72 
80  function integ_scal_func(u, qdm) result(int)
81  real(RP) :: int
82  procedure(r3tor) :: u
83  type(quadmesh), intent(in) :: qdm
84 
85  type(mesh), pointer :: m
86  integer :: qt, dim, nn
87 
88  if (choral_verb>1) write(*,*) &
89  & 'integral : integ_scal_func'
90 
91  if(.NOT.valid(qdm)) call quit( &
92  &"integral: integ_scal_func:&
93  & quad mesh 'qdmh' not valid" )
94 
95  int = 0.0_rp
96  m => qdm%mesh
97 
98  do qt=1, quad_tot_nb
99  if (qdm%quad_count(qt)==0) cycle
100 
101  dim = quad_dim(qt) ! dimension of K_ref
102  nn = quad_nbnodes(qt) ! quad number of nodes
103 
104  call cell_loop()
105 
106  end do
107 
108  contains
109 
110  subroutine cell_loop()
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
116 
117  type(geotsf):: g
118  integer :: cl, jj, ll, cl_nd
119  real(RP) :: int2
120 
121  int2 = 0._rp
122 
123  !! quad nodes and weights
124  !!
125  wgt = quad_wgt(qt)%y
126  y = quad_coord(qt)%y
127 
128  g = geotsf(dim, nn, y)
129 
130  ! Loop on mesh cells
131  do cl=1, m%nbCl
132 
133  jj = qdm%quadType(cl)
134  if (jj/=qt) cycle
135 
136  ! cell node coordinates
137  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
138  do ll=1, cl_nd
139  x(:,ll) = m%nd( :, p_cl(ll) )
140  end do
141 
142  ! transformation T : K_ref --> K = T(K_ref)
143  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
144  call compute_j(g)
145 
146  ! loop on quadrature nodes
147  do ll=1, nn
148  u_x(ll) = u(g%Ty(:,ll))
149  end do
150  u_x = u_x * g%Jy
151  u_x = u_x * wgt
152  int2 = int2 + sum(u_x)
153 
154  end do
155 
156  int = int + int2
157 
158  end subroutine cell_loop
159 
160  end function integ_scal_func
161 
162 
163 
164 
165 
166  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167  !!
168  !!
169  !! SCALAR FINITE ELEMENTS
170  !!
171  !!
172  !!
173  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174 
175 
192  function integ_scal_fe(E, u, uh, X_h, qdm) result(int)
193  real(RP) :: int
194  procedure(r3xrxrtor) :: E
195  procedure(r3tor) :: u
196  real(RP), dimension(:), intent(in) :: uh
197  type(fespace) , intent(in) :: X_h
198  type(quadmesh) , intent(in) :: qdm
199 
200  type(mesh), pointer :: m
201  integer :: dim, nbDof, nn, ft, qt
202 
203  if (choral_verb>1) write(*,*) &
204  & 'integral : integ_scal_fe'
205 
206  if(.NOT.valid(x_h)) call quit( &
207  &"integral: integ_scal_fe:&
208  & fe space 'X_h' not valid" )
209 
210  if(.NOT.valid(qdm)) call quit( &
211  &"integral: integ_scal_fe:&
212  & quad mesh 'qdmh' not valid" )
213 
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" )
217 
218  if(size(uh,1) /= x_h%nbDof) call quit( &
219  &"integral: integ_scal_fe:&
220  & fe function 'uh' has not the expected size" )
221 
222  int = 0._rp
223 
224  m => x_h%mesh
225 
226  do ft=1, fe_tot_nb
227  if (x_h%fe_count(ft)==0) cycle
228 
229  do qt=1, quad_tot_nb
230  if (qdm%quad_count(qt)==0) cycle
231  if (quad_geo(qt) /= fe_geo(ft) ) cycle
232 
233  dim = quad_dim(qt) ! dimension of K_ref
234  nbdof = fe_nbdof(ft) ! fe number of DOF
235  nn = quad_nbnodes(qt) ! quad number of nodes
236 
237  call cell_loop()
238 
239  end do
240  end do
241 
242  contains
243 
244  subroutine cell_loop()
245 
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
253 
254  type(geotsf):: g
255  integer :: cl, ii, jj, ll, cl_nd
256  real(RP) :: int2, uh_x, u_x
257 
258  int2 = 0._rp
259 
260  !! quad nodes and weights
261  !!
262  wgt = quad_wgt(qt)%y
263  y = quad_coord(qt)%y
264 
265  !! fe basis functions at quad nodes
266  !!
267  do ll=1, nn
268  call scal_fe(val(:,ll), nbdof, y(:,ll), dim, ft)
269  end do
270 
271  !! geometricat transformation K_ref --> K
272  !!
273  g = geotsf(dim, nn, y)
274 
275  !! CELL LOOP
276  !!
277  do cl=1, m%nbCl
278 
279  jj = qdm%quadType(cl)
280  if (jj/=qt) cycle
281 
282  ii = x_h%feType(cl)
283  if (ii/=ft) cycle
284 
285  ! cell node coordinates
286  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
287  do ll=1, cl_nd
288  x(:,ll) = m%nd( :, p_cl(ll) )
289  end do
290 
291  ! transformation T : K_ref --> K = T(K_ref)
292  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
293  call compute_j(g)
294 
295  ! dof associated with cell cl indexes
296  call getrow(p, nbdof, x_h%clToDof, cl)
297 
298  ! local dof of the fe function uh
299  v = uh(p)
300 
301  ! loop on quad nodes
302  do ll=1, nn
303  ! x = T(node ll)
304 
305  uh_x = dot_product(val(:,ll), v) ! u(x)
306  u_x = u(g%Ty(:,ll)) ! uh(x)
307  e_x(ll)= e(g%Ty(:,ll), u_x, uh_x) ! E(x, u(x), uh(x))
308 
309  end do
310 
311  ! update integral
312  e_x = e_x * g%Jy
313  e_x = e_x * wgt
314  int2 = int2 + sum(e_x)
315 
316  end do
317 
318  int = int + int2
319 
320  end subroutine cell_loop
321 
322  end function integ_scal_fe
323 
324 
325 
326 
327 
344  function integ_scal_fe_grad(E, phi, uh, X_h, qdm) result(int)
345  real(RP) :: int
346  procedure(r3xr3xr3tor) :: E
347  real(RP) , dimension(:), intent(in) :: uh
348  type(fespace) , intent(in) :: X_h
349  type(quadmesh) , intent(in) :: qdm
350  procedure(r3tor3) :: phi
351 
352  type(mesh), pointer :: m
353  integer :: dim, nbDof, nn, ft, qt
354 
355  if (choral_verb>1) write(*,*) &
356  & 'integral : Integ_scal_fe_grad'
357 
358  if(.NOT.valid(x_h)) call quit( &
359  &"integral: integ_scal_fe_grad:&
360  & fe space 'X_h' not valid" )
361 
362  if(.NOT.valid(qdm)) call quit( &
363  &"integral: integ_scal_fe_grad:&
364  & quad mesh 'qdmh' not valid" )
365 
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" )
369 
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" )
373 
374  int = 0._rp
375 
376  m => x_h%mesh
377 
378  do ft=1, fe_tot_nb
379  if (x_h%fe_count(ft)==0) cycle
380 
381  do qt=1, quad_tot_nb
382  if (qdm%quad_count(qt)==0) cycle
383  if (quad_geo(qt) /= fe_geo(ft) ) cycle
384 
385  dim = quad_dim(qt) ! dimension of K_ref
386  nbdof = fe_nbdof(ft) ! fe number of DOF
387  nn = quad_nbnodes(qt) ! quad number of nodes
388 
389  call cell_loop()
390 
391  end do
392  end do
393 
394  contains
395 
396  subroutine cell_loop()
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
405 
406  type(geotsf):: g
407  integer :: cl, ii, jj, ll, cl_nd
408  real(RP) :: int2, E_x
409 
410  int2 = 0._rp
411 
412  !! quad nodes and weights
413  !!
414  wgt = quad_wgt(qt)%y
415  y = quad_coord(qt)%y
416 
417  !! gradient of the fe basis functions
418  !! evaluated at quad nodes
419  !!
420  do ll=1, nn
421  call scalgrad_fe(val(:,:,ll), nbdof, y(:,ll), dim, ft)
422  end do
423 
424  !! geometricat transformation K_ref --> K
425  !!
426  g = geotsf(dim, nn, y)
427 
428  !! CELL LOOP
429  !!
430  do cl=1, m%nbCl
431 
432  jj = qdm%quadType(cl)
433  if (jj/=qt) cycle
434 
435  ii = x_h%feType(cl)
436  if (ii/=ft) cycle
437 
438  ! cell node coordinates
439  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
440  do ll=1, cl_nd
441  x(:,ll) = m%nd( :, p_cl(ll) )
442  end do
443 
444  ! transformation T : K_ref --> K = T(K_ref)
445  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
446  call compute_jc(g)
447 
448  ! dof associated with cell cl indexes
449  call getrow(p, nbdof, x_h%clToDof, cl)
450 
451  ! local dof of the fe function uh
452  v = uh(p)
453 
454  ! loop on quad nodes
455  do ll=1, nn
456  ! x = T(node ll)
457 
458  ! compute the gradient
459  grad_uh_x(:) = 0._rp
460  select case(dim)
461  case(3)
462  do jj=1, nbdof
463  call matvecprod_3x3( w, &
464  & g%Cy(:,:,ll), val(:,jj,ll))
465  grad_uh_x = grad_uh_x + v(jj) * w
466  end do
467 
468  case(2)
469  do jj=1, nbdof
470  call matvecprod_3x2( w, &
471  & g%Cy(:,:,ll), val(:,jj,ll) )
472  grad_uh_x = grad_uh_x + v(jj) * w
473  end do
474 
475  case(1)
476  do jj=1, nbdof
477  call matvecprod_3x1( w, &
478  & g%Cy(:,:,ll), val(:,jj,ll) )
479  grad_uh_x = grad_uh_x + v(jj) * w
480  end do
481 
482  end select
483  grad_uh_x = grad_uh_x / g%Jy(ll)
484 
485  ! phi(x)
486  phi_x = phi(g%Ty(:,ll))
487 
488  ! update integral
489  e_x = e(g%Ty(:,ll), phi_x, grad_uh_x)
490  int2 = int2 + e_x * wgt(ll) * g%Jy(ll)
491 
492  end do
493  end do
494 
495  int = int + int2
496 
497  end subroutine cell_loop
498 
499  end function integ_scal_fe_grad
500 
501 
510  function integ_scal_fe_grad_proj(E, phi, uh, X_h, qdm) result(int)
511  real(RP) :: int
512  procedure(r3xr3xr3tor) :: E
513  real(RP) , dimension(:), intent(in) :: uh
514  type(fespace) , intent(in) :: X_h
515  type(quadmesh) , intent(in) :: qdm
516  procedure(r3tor3) :: phi
517 
518  type(mesh), pointer :: m
519  integer :: dim, nbDof, nn, ft, qt
520 
521  if (choral_verb>1) write(*,*) &
522  & 'integral : Integ_scal_fe_grad_proj'
523 
524  if(.NOT.valid(x_h)) call quit( &
525  &"integral: integ_scal_fe_grad_proj:&
526  & fe space 'X_h' not valid" )
527 
528  if(.NOT.valid(qdm)) call quit( &
529  &"integral: integ_scal_fe_grad_proj:&
530  & quad mesh 'qdmh' not valid" )
531 
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" )
535 
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" )
539 
540  int = 0._rp
541 
542  m => x_h%mesh
543 
544  do ft=1, fe_tot_nb
545  if (x_h%fe_count(ft)==0) cycle
546 
547  do qt=1, quad_tot_nb
548  if (qdm%quad_count(qt)==0) cycle
549  if (quad_geo(qt) /= fe_geo(ft) ) cycle
550 
551  dim = quad_dim(qt) ! dimension of K_ref
552  nbdof = fe_nbdof(ft) ! fe number of DOF
553  nn = quad_nbnodes(qt) ! quad number of nodes
554 
555  call cell_loop()
556 
557  end do
558  end do
559 
560  contains
561 
562  subroutine cell_loop()
563 
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
572 
573  type(geotsf):: g
574  integer :: cl, ii, jj, ll, cl_nd
575  real(RP) :: int2, E_x
576 
577  int2 = 0._rp
578 
579  !! quad nodes and weights
580  !!
581  wgt = quad_wgt(qt)%y
582  y = quad_coord(qt)%y
583 
584 
585  !! gradient of the fe basis functions
586  !! evaluated at quad nodes
587  !!
588  do ll=1, nn
589  call scalgrad_fe(val(:,:,ll), nbdof, y(:,ll), dim, ft)
590  end do
591 
592  !! geometricat transformation K_ref --> K
593  !!
594  g = geotsf(dim, nn, y)
595 
596  !! CELL LOOP
597  !!
598  do cl=1, m%nbCl
599 
600  jj = qdm%quadType(cl)
601  if (jj/=qt) cycle
602 
603  ii = x_h%feType(cl)
604  if (ii/=ft) cycle
605 
606  ! cell node coordinates
607  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
608  do ll=1, cl_nd
609  x(:,ll) = m%nd( :, p_cl(ll) )
610  end do
611 
612  ! transformation T : K_ref --> K = T(K_ref)
613  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
614  call compute_dtjc(g)
615 
616  ! dof associated with cell cl indexes
617  call getrow(p, nbdof, x_h%clToDof, cl)
618 
619  ! local dof of the fe function uh
620  v = uh(p)
621 
622  ! loop on quad nodes
623  do ll=1, nn
624  ! x = T(node ll)
625 
626  ! compute the gradient
627  grad_uh_x(:) = 0._rp
628  select case(dim)
629  case(3)
630  do jj=1, nbdof
631  call matvecprod_3x3( w, &
632  & g%Cy(:,:,ll), val(:,jj,ll))
633  grad_uh_x = grad_uh_x + v(jj) * w
634  end do
635 
636  case(2)
637  do jj=1, nbdof
638  call matvecprod_3x2( w, &
639  & g%Cy(:,:,ll), val(:,jj,ll) )
640  grad_uh_x = grad_uh_x + v(jj) * w
641  end do
642 
643  case(1)
644  do jj=1, nbdof
645  call matvecprod_3x1( w, &
646  & g%Cy(:,:,ll), val(:,jj,ll) )
647  grad_uh_x = grad_uh_x + v(jj) * w
648  end do
649 
650  end select
651  grad_uh_x = grad_uh_x / g%Jy(ll)
652 
653  ! phi(x)
654  phi_x = phi(g%Ty(:,ll))
655 
656  ! tangential component phi_T(x)
657  if (dim==3) then
658  ! nothing
659 
660  else if (dim==2) then
661  w = crossprod_3d( g%DTy(:,1,ll), g%DTy(:,2,ll))
662  w = w / norm2(w)
663  phi_x = phi_x - sum(phi_x * w)*w
664 
665  else if (dim==1) then
666  phi_x = sum( g%Cy(:,1,ll)*phi_x ) * g%Cy(:,1,ll)
667 
668  end if
669 
670  ! update integral
671  e_x = e(g%Ty(:,ll), phi_x, grad_uh_x)
672  int2 = int2 + e_x * wgt(ll) * g%Jy(ll)
673 
674  end do
675  end do
676 
677  int = int + int2
678 
679  end subroutine cell_loop
680 
681  end function integ_scal_fe_grad_proj
682 
683 
684 
685 
686 
701  function integ_scal_fe_grad_2(E, uh1, uh2, X_h, qdm) result(int)
702  real(RP) :: int
703  procedure(r3xr3xr3tor) :: E
704  real(RP) , dimension(:), intent(in) :: uh1, uh2
705  type(fespace) , intent(in) :: X_h
706  type(quadmesh) , intent(in) :: qdm
707 
708  type(mesh), pointer :: m
709  integer :: dim, nbDof, nn, ft, qt
710 
711  if (choral_verb>1) write(*,*) &
712  & 'integral : Integ_scal_fe_grad_2'
713 
714  if(.NOT.valid(x_h)) call quit( &
715  &"integral: integ_scal_fe_grad_2:&
716  & fe space 'X_h' not valid" )
717 
718  if(.NOT.valid(qdm)) call quit( &
719  &"integral: integ_scal_fe_grad_2:&
720  & quad mesh 'qdmh' not valid" )
721 
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" )
725 
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" )
729 
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" )
733 
734 
735  int = 0._rp
736 
737  m => x_h%mesh
738 
739  do ft=1, fe_tot_nb
740  if (x_h%fe_count(ft)==0) cycle
741 
742  do qt=1, quad_tot_nb
743  if (qdm%quad_count(qt)==0) cycle
744  if (quad_geo(qt) /= fe_geo(ft) ) cycle
745 
746  dim = quad_dim(qt) ! dimension of K_ref
747  nbdof = fe_nbdof(ft) ! fe number of DOF
748  nn = quad_nbnodes(qt) ! quad number of nodes
749 
750  call cell_loop()
751 
752  end do
753  end do
754 
755  contains
756 
757  subroutine cell_loop()
758 
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
767 
768  type(geotsf):: g
769  integer :: cl, ii, jj, ll, cl_nd
770  real(RP) :: int2, E_x
771 
772  int2 = 0._rp
773 
774  !! quad nodes and weights
775  !!
776  wgt = quad_wgt(qt)%y
777  y = quad_coord(qt)%y
778 
779 
780  !! gradient of the fe basis functions
781  !! evaluated at quad nodes
782  !!
783  do ll=1, nn
784  call scalgrad_fe(val(:,:,ll), nbdof, y(:,ll), dim, ft)
785  end do
786 
787  !! geometricat transformation K_ref --> K
788  !!
789  g = geotsf(dim, nn, y)
790 
791  !! CELL LOOP
792  !!
793  do cl=1, m%nbCl
794 
795  jj = qdm%quadType(cl)
796  if (jj/=qt) cycle
797 
798  ii = x_h%feType(cl)
799  if (ii/=ft) cycle
800 
801  ! cell node coordinates
802  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
803  do ll=1, cl_nd
804  x(:,ll) = m%nd( :, p_cl(ll) )
805  end do
806 
807  ! transformation T : K_ref --> K = T(K_ref)
808  call assemble(g, m%clType(cl), x(1:3,1:cl_nd), cl_nd)
809  call compute_jc(g)
810 
811  ! dof associated with cell cl indexes
812  call getrow(p, nbdof, x_h%clToDof, cl)
813 
814  ! local dof of the fe function uh
815  v1 = uh1(p)
816  v2 = uh2(p)
817 
818  ! loop on quad nodes
819  do ll=1, nn
820  ! x = T(node ll)
821 
822  ! compute the gradient
823  grad_uh1_x = 0._rp
824  grad_uh2_x = 0._rp
825  select case(dim)
826  case(3)
827  do jj=1, nbdof
828  call matvecprod_3x3( w, &
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
832  end do
833 
834  case(2)
835  do jj=1, nbdof
836  call matvecprod_3x2( 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
840 
841  end do
842 
843  case(1)
844  do jj=1, nbdof
845  call matvecprod_3x1( 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
849 
850  end do
851 
852  end select
853  grad_uh1_x = grad_uh1_x / g%Jy(ll)
854  grad_uh2_x = grad_uh2_x / g%Jy(ll)
855 
856  ! update integral
857  e_x = e(g%Ty(:,ll), grad_uh1_x, grad_uh2_x)
858  int2 = int2 + e_x * wgt(ll) * g%Jy(ll)
859 
860  end do
861  end do
862 
863  int = int + int2
864 
865  end subroutine cell_loop
866 
867  end function integ_scal_fe_grad_2
868 
869 
870 
871 
874  function e_l2(x, u1, u2) result(r)
875  real(RP) :: r
876  real(RP), dimension(3), intent(in) :: x
877  real(RP), intent(in) :: u1, u2
878  r = (u1 - u2)**2
879  end function e_l2
880 
883  function e_l2_vect(x, p1, p2) result(r)
884  real(RP) :: r
885  real(RP), dimension(3), intent(in) :: x, p1, p2
886  r = sum( (p1-p2)**2 )
887  end function e_l2_vect
888 
889 
904  function l2_dist_scalfe(u, uh, X_h, qdm) result(dist)
905  real(RP) :: dist
906  real(RP), dimension(:), intent(in) :: uh
907  type(fespace) , intent(in) :: X_h
908  type(quadmesh) , intent(in) :: qdm
909  procedure(r3tor) :: u
910 
911  if (choral_verb>1) write(*,*) &
912  & 'integral : L2_dist_scalFE'
913  dist = integ(e_l2, u, uh, x_h, qdm)
914  dist = sqrt(dist)
915 
916  end function l2_dist_scalfe
917 
918 
933  function l2_dist_grad_scalfe(phi, uh, X_h, qdm) result(dist)
934  real(RP) :: dist
935  real(RP), dimension(:), intent(in) :: uh
936  type(fespace) , intent(in) :: X_h
937  type(quadmesh) , intent(in) :: qdm
938  procedure(r3tor3) :: phi
939 
940  if (choral_verb>1) write(*,*) &
941  & 'integral : L2_dist_grad_scalFE'
942  dist = integ(e_l2_vect, phi, uh, x_h, qdm)
943  dist = sqrt(dist)
944 
945  end function l2_dist_grad_scalfe
946 
955  function l2_dist_grad_proj(phi, uh, X_h, qdm) result(dist)
956  real(RP) :: dist
957  real(RP), dimension(:), intent(in) :: uh
958  type(fespace) , intent(in) :: X_h
959  type(quadmesh) , intent(in) :: qdm
960  procedure(r3tor3) :: phi
961 
962  if (choral_verb>1) write(*,*) &
963  & 'integral : L2_dist_grad_proj'
964  dist = integ_scal_fe_grad_proj(e_l2_vect, phi, uh, x_h, qdm)
965  dist = sqrt(dist)
966 
967  end function l2_dist_grad_proj
968 
969 
970 
985  function l2_dist_vect(phi, phi_h, X_h, qdm) result(dist)
986  real(RP) :: dist
987  procedure(r3tor3) :: phi
988  real(RP), dimension(:), intent(in) :: phi_h
989  type(fespace) , intent(in) :: X_h
990  type(quadmesh) , intent(in) :: qdm
991 
992  type(mesh), pointer :: m
993  integer :: dim, nbDof, nn, ft, qt
994 
995  if (choral_verb>1) write(*,*) &
996  & "integral : L2_dist_vect"
997 
998  if(.NOT.valid(x_h)) call quit( &
999  &"integral: L2_dist_vect:&
1000  & fe space 'X_h' not valid" )
1001 
1002  if(.NOT.valid(qdm)) call quit( &
1003  &"integral: L2_dist_vect:&
1004  & quad mesh 'qdmh' not valid" )
1005 
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" )
1009 
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" )
1013 
1014  m => x_h%mesh
1015  if (m%nbItf<=0) call quit( "integral: L2_dist_vect:&
1016  & mesh interfaces not defined")
1017 
1018  dist = 0.0_rp
1019 
1020  do ft=1, fe_tot_nb
1021  if (x_h%fe_count(ft)==0) cycle
1022 
1023  do qt=1, quad_tot_nb
1024  if (qdm%quad_count(qt)==0) cycle
1025  if (quad_geo(qt) /= fe_geo(ft) ) cycle
1026 
1027  dim = quad_dim(qt) ! dimension of K_ref
1028  nbdof = fe_nbdof(ft) ! fe number of DOF
1029  nn = quad_nbnodes(qt) ! quad number of nodes
1030 
1031  call cell_loop()
1032 
1033  end do
1034  end do
1035 
1036  dist = sqrt( dist )
1037 
1038  contains
1039 
1040  subroutine cell_loop()
1041 
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
1052 
1053  type(geotsf):: g
1054  integer :: cl, ii, jj, ll, cl_nd, nb_itf
1055  real(RP) :: s
1056 
1057  !! quad nodes and weights
1058  !!
1059  wgt = quad_wgt(qt)%y
1060  y = quad_coord(qt)%y
1061 
1062  do ll=1, nn
1063  call vect_fe(base(:,:,ll), nbdof, y(:,ll), dim, ft)
1064  end do
1065 
1066  !! geometricat transformation K_ref --> K
1067  !!
1068  g = geotsf(dim, nn, y)
1069 
1070  !! CELL LOOP
1071  !!
1072  do cl=1, m%nbCl
1073 
1074  jj = qdm%quadType(cl)
1075  if (jj/=qt) cycle
1076 
1077  ii = x_h%feType(cl)
1078  if (ii/=ft) cycle
1079 
1080  !! p = dof indexes associated with cell cl
1081  call getrow(p, nbdof, x_h%clToDof, cl)
1082 
1083  !! v = local dof
1084  v = phi_h( p )
1085 
1086  !! Modification of the local dof
1087  !! to take into account the interface orientation
1088 #if DBG>0
1089  !!
1090  !! checks the cell orientation
1091  if (.NOT.cell_orientation(m, cl)) call quit( &
1092  & "integral : L2_dist_vect:&
1093  & cell orientation error" )
1094 #endif
1095  !!
1096  !! orientation of the interfaces of cell cl
1097  call interface_orientation(nb_itf, &
1098  & eps, cell_max_nbitf, m, cl)
1099  !!
1100  select case(ft)
1101  case(fe_rt0_1d, fe_rt0_2d)
1102  do jj=1, nbdof
1103  if (eps(jj)==1) cycle
1104  v(jj) = -v(jj)
1105  end do
1106 
1107  case(fe_rt1_1d, fe_rt1_1d_2)
1108  do jj=1, 2
1109  if (eps(jj)==1) cycle
1110  v(jj) = -v(jj)
1111  end do
1112 
1113  case(fe_rt1_2d_2)
1114  do jj=1, 3
1115  if (eps(jj)==1) cycle
1116  ll = 2*jj-1
1117  v(ll:ll+1) = -v(ll:ll+1)
1118  end do
1119 
1120  case default
1121  call quit( "integral: L2_dist_vect:&
1122  & fe type not supported" )
1123 
1124  end select
1125 
1126 
1127  ! cell node coordinates
1128  call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
1129  do ll=1, cl_nd
1130  x(1:3,ll) = m%nd( 1:3, p_cl(ll) )
1131  end do
1132 
1133  ! transformation T : K_ref --> K = T(K_ref)
1134  call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)
1135  call compute_dtj(g)
1136 
1137  !! DT_base(:,ii) \in R^3 =
1138  !! DT(:,:,ll).base(:,ii,ll)
1139  !! with DT(:,:,ll) = DT at node y(:,ll)
1140  !!
1141  select case(dim)
1142 
1143  case(3)
1144  do ll=1, nn
1145  do ii=1, nbdof
1146  call matvecprod_3x3( dt_base(:,ii,ll), &
1147  & g%DTy(:,:,ll), base(:,ii,ll) )
1148  end do
1149  end do
1150 
1151  case(2)
1152  do ll=1, nn
1153  do ii=1, nbdof
1154  call matvecprod_3x2( dt_base(:,ii,ll), &
1155  & g%DTy(:,:,ll), base(:,ii,ll) )
1156  end do
1157  end do
1158 
1159  case(1)
1160  do ll=1, nn
1161  do ii=1, nbdof
1162  call matvecprod_3x1( dt_base(:,ii,ll), &
1163  & g%DTy(:,:,ll), base(:,ii,ll) )
1164  end do
1165  end do
1166 
1167  end select
1168 
1169 
1170  !! Update integral
1171  !!
1172  do ll=1, nn
1173  ! x = Ty(:,ll)
1174  ! phi_h( x )
1175  phi_h_x = dt_base(:,1,ll)*v(1)
1176  do jj=2, nbdof
1177  phi_h_x = phi_h_x + dt_base(:,jj,ll)*v(jj)
1178  end do
1179  s = 1.0_rp / g%Jy(ll)
1180  phi_h_x = phi_h_x * s
1181 
1182 
1183  !! phi_x = phi(x)
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)
1187 
1188  end do
1189 
1190  end do
1191 
1192  end subroutine cell_loop
1193 
1194  end function l2_dist_vect
1195 
1196 
1197 end module integral
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
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
real(rp) function integ_scal_func(u, qdm)
Returns .
Definition: integral.F90:81
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
real(rp) function integ_scal_fe_grad_2(E, uh1, uh2, X_h, qdm)
Returns .
Definition: integral.F90:702
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)
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 matvecprod_3x3(y, A, X)
matrix vector product R3 –> R3
DEFINITION OF FINITE ELEMENT METHODS
Definition: fe_mod.F90:19
integer, parameter, public rp
real(kind=RP) = real precision in the code REAL_TOL = epsilon to test real equality ...
Definition: real_type.F90:90
real(rp) function, public l2_dist_vect(phi, phi_h, X_h, qdm)
Returns .
Definition: integral.F90:986
subroutine, public quit(message)
Stop program execution, display an error messahe.
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...
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(rp) function l2_dist_scalfe(u, uh, X_h, qdm)
Returns .
Definition: integral.F90:905
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
real(rp) function e_l2(x, u1, u2)
Definition: integral.F90:875
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.
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
integral L2 distance
Definition: integral.F90:63
CHORAL CONSTANTS
subroutine, public vect_fe(val, nbDof, x, dim, ft)
Evaluate vector FE basis functions at point x.
Definition: fe_mod.F90:195
Integration of scalar functions.
Definition: integral.F90:54
real(rp) function, public l2_dist_grad_proj(phi, uh, X_h, qdm)
Returns .
Definition: integral.F90:956
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 .
Definition: integral.F90:345
DERIVED TYPE mesh
Definition: mesh_mod.F90:57
real(rp) function integ_scal_fe(E, u, uh, X_h, qdm)
Returns .
Definition: integral.F90:193
real(rp) function l2_dist_grad_scalfe(phi, uh, X_h, qdm)
Returns .
Definition: integral.F90:934
subroutine, public compute_dtjc(g)
Compute DTy(:,:,i), Jy(i) and Cy(:,:,i) at the nodes y(:,i)
Definition: geoTsf_mod.F90:490
Derived type feSpace: finite element space on a mesh.
COMPUTATION OF INTEGRALS using a quadrature methods on a mesh.
Definition: integral.F90:21
integer choral_verb
Verbosity level.
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
Definition: geoTsf_mod.F90:428
real(rp) function e_l2_vect(x, p1, p2)
Definition: integral.F90:884
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
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 .
Definition: integral.F90:511
subroutine cell_loop()
Definition: diffusion.F90:232
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 &#39;cl&#39; the direct orientation ?
Definition: mesh_mod.F90:777
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.