Choral
feSpacexk_mod.f90
Go to the documentation of this file.
1 
29 !
32 
34  use real_type
35  use basic_tools
36  use abstract_interfaces, only: r3tor, r3tor3
37  use graph_mod
38  use mesh_mod
39  use fespace_mod
40  use quadmesh_mod
41  use integral
42 
43  implicit none
44  private
45 
46 
47  ! %----------------------------------------%
48  ! | |
49  ! | PUBLIC DATA |
50  ! | |
51  ! %----------------------------------------%
52  public :: fespacexk
53 
54  public :: clear
55  public :: valid
56  public :: print
57  public :: interp_vect_func !! tested
58  public :: extract_component !! tested
59  public :: l2_dist
60  public :: l2_dist_grad
61 
62  ! %----------------------------------------%
63  ! | |
64  ! | DERIVED TYPE |
65  ! | |
66  ! %----------------------------------------%
94  type fespacexk
95 
97  type(fespace), pointer :: x_h=>null()
98 
100  integer :: k
101 
103  integer :: nbdof2=0
104 
106  integer , dimension(:), allocatable :: dof2todof
107 
108 
109  !! cell dof local indexing --> dof2 global indexing
110  !! cell to dof connectivity
111  type(graph) :: cltodof2
112 
113  contains
114 
117 
118  end type fespacexk
119 
120 
121  ! %----------------------------------------%
122  ! | |
123  ! | GENERIc SUBROUTINES |
124  ! | |
125  ! %----------------------------------------%
127  interface clear
128  module procedure fespacexk_clear
129  end interface clear
130 
131 
132  !! constructor
133  interface fespacexk
134  module procedure fespacexk_create
135  end interface fespacexk
136 
137  interface valid
138  module procedure fespacexk_valid
139  end interface valid
140 
142  module procedure fespacexk_interp_vect_func
143  end interface interp_vect_func
144 
146  interface print
147  module procedure fespacexk_print
148  end interface print
149 
152  interface l2_dist
153  module procedure l2_dist_2
154  end interface l2_dist
155  interface l2_dist_grad
156  module procedure l2_dist_grad_2
157  end interface l2_dist_grad
158 
159 contains
160 
162  subroutine fespacexk_clear(Y)
163  type(fespacexk), intent(inout) :: Y
164 
165  y%X_h => null()
166  y%nbDof2 = 0
167  call freemem(y%dof2ToDof)
168  call clear( y%clToDof2)
169 
170  end subroutine fespacexk_clear
171 
183  function fespacexk_create(X_h, k) result(Y)
184  type(fespacexk) :: Y
185  type(fespace) , target :: X_h
186  integer , intent(in) :: k
187 
188  integer :: cl, ii, jj, ll, wdt, cpt, nn
189  integer, dimension(:), allocatable :: nnz, row, row2
190 
191  if (choral_verb>0) write(*,*) &
192  & 'feSpacexk_mod : create k =', int(k,1)
193  call clear(y)
194  if (k<1) call quit( &
195  & "feSpacexk_mod: feSpacexk_create:&
196  & wrong argument 'k'" )
197  if (.NOT.valid(x_h)) call quit( &
198  &"feSpacexk_mod: feSpacexk_create:&
199  & fe space 'X_h' not valid" )
200 
201  !! parameters
202  y%X_h => x_h
203  y%k = k
204  y%nbDof2 = x_h%nbDof * k
205 
206  !! correspondance dof of Y --> dof of X_h
207  call allocmem( y%dof2ToDof, y%nbDof2)
208  do ii=1, y%nbDof2
209  y%dof2ToDof(ii) = (ii-1)/k + 1
210  end do
211 
212  !! connectivity cell --> dof of Y
213  !!
214  call allocmem( nnz, y%X_h%clToDof%nl )
215  do cl=1, x_h%clToDof%nl
216 
217  jj = x_h%clToDof%row(cl)
218  ll = x_h%clToDof%row(cl+1)
219  nnz(cl) = ll - jj
220 
221  end do
222  nnz = nnz * k
223  y%clToDof2 = graph(nnz, size(nnz,1), y%nbDof2)
224  call freemem(nnz)
225  !!
226  !! fill in the graph
227  !!
228  wdt = x_h%clToDof%width
229  call allocmem(row , wdt)
230  call allocmem(row2, wdt*k)
231  do cl=1, x_h%clToDof%nl
232 
233  call getrow(ll, row, wdt, x_h%clToDof, cl)
234  cpt = 0
235  do ii=1, ll
236  nn = (row(ii)-1) * k
237  do jj=1, k
238  cpt = cpt + 1
239  row2(cpt) = nn + jj
240  end do
241  end do
242  nn = ll * k
243  call setrow(y%clToDof2, row2(1:nn), nn, cl)
244 
245  end do
246 
247  if (.NOT.valid(y)) call quit(&
248  & 'feSpacexk_mod : feSpacexk_create:&
249  & assembling not valid' )
250 
251  end function fespacexk_create
252 
253 
256  function fespacexk_valid(Y) result(b)
257  type(fespacexk), intent(in) :: y
258  logical :: b
259 
260  b = .false.
261  if(.NOT. valid(y%X_h) ) return
262  if( valid(y%clToDof2) < 0 ) return
263  if(.NOT. allocated(y%dof2ToDof) ) return
264 
265  b = y%nbdof2>0
266  b = b .AND. ( y%k >0 )
267  b = b .AND. ( (y%k * y%X_h%nbDof ) == y%nbdof2 )
268  b = b .AND. ( size( y%dof2ToDof, 1) == y%nbDof2 )
269 
270  end function fespacexk_valid
271 
272 
273 
281  subroutine fespacexk_interp_vect_func(u_h, Y, u_1, u_2, u_3)
282  real(RP), dimension(:), allocatable :: u_h
283  type(fespacexk) , intent(in) :: Y
284  procedure(r3tor) :: u_1, u_2
285  procedure(r3tor) , optional :: u_3
286 
287  type(fespace), pointer :: X_h
288  integer :: i, dim, jj, ll
289  real(RP), dimension(:), allocatable :: u_ih
290 
291  if (choral_verb>1) write(*,*) &
292  & "feSpacexk_mod : interp_vetc_func"
293 
294  if (.NOT.valid(y)) call quit(&
295  &"feSpacexk_mod: feSpacexk_interp_vect_func:&
296  & fespacsexk 'Y' not valid" )
297 
298  x_h => y%X_h
299  dim = x_h%mesh%dim
300  if (y%k /= dim) call quit(&
301  &"feSpacexk_mod: feSpacexk_interp_vect_func:&
302  & fespacsexk exponent 'Y%k' incorrect" )
303  if ( dim==1 ) call quit(&
304  &"feSpacexk_mod: feSpacexk_interp_vect_func:&
305  & dim 1 case, use feSpace_mod::interp_scal_func" )
306  if ( (dim==3).AND.(.NOT.present(u_3)) ) call quit(&
307  &"feSpacexk_mod: feSpacexk_interp_vect_func:&
308  & argument 'u3' missing" )
309 
310  call allocmem(u_h , y%nbDof2)
311 
312  !! first component
313  i = 1
314  call interp_scal_func(u_ih, u_1, x_h)
315  do jj=1, x_h%nbDof
316  ll = (jj-1) * y%k + i
317  u_h(ll) = u_ih(jj)
318  end do
319 
320  !! second component
321  i = 2
322  call interp_scal_func(u_ih, u_2, x_h)
323  do jj=1, x_h%nbDof
324  ll = (jj-1) * y%k + i
325  u_h(ll) = u_ih(jj)
326  end do
327  if (dim==2) return
328 
329  !! third component
330  i = 3
331  call interp_scal_func(u_ih, u_3, x_h)
332  do jj=1, x_h%nbDof
333  ll = (jj-1) * y%k + i
334  u_h(ll) = u_ih(jj)
335  end do
336 
337  end subroutine fespacexk_interp_vect_func
338 
341  subroutine fespacexk_print(Y)
342  type(fespacexk), intent(in) :: Y
343 
344  write(*,*)"feSpacexk_mod : print"
345  write(*,*)" Number of DOF =", y%nbDof2
346  write(*,*)" Exponent =", y%k
347  write(*,*)" Valid =", valid(y)
348 
349  end subroutine fespacexk_print
350 
351 
374  subroutine extract_component(u_c, u, Y, c)
375  real(RP), dimension(:), allocatable :: u_c
376  real(RP), dimension(:), intent(in) :: u
377  type(fespacexk) , intent(in) :: Y
378  integer , intent(in) :: c
379 
380  integer :: k, jj, ll
381 
382  if (choral_verb>1) write(*,*) &
383  & "feSpacexk_mod : extract_component"
384 
385  if (.NOT.valid(y)) call quit(&
386  &"feSpacexk_mod: extract_component:&
387  & fespacsexk 'Y' not valid" )
388 
389  if ( (c <1) .OR. (c> y%k) ) call quit(&
390  &"feSpacexk_mod: extract_component:&
391  & component 'c' incorrect" )
392 
393  if ( size(u,1) /= y%nbDof2 ) call quit(&
394  &"feSpacexk_mod: extract_component:&
395  & wrong size for the finite element function 'u'" )
396 
397  call allocmem(u_c, y%X_h%nbDof)
398 
399  k = y%k
400  ll = c
401  do jj=1, y%X_h%nbDof
402  u_c(jj) = u(ll)
403  ll = ll + k
404  end do
405 
406  end subroutine extract_component
407 
408 
430  function l2_dist_2(uh, Y, qdm, u_1, u_2, u_3) result(dist)
431  real(RP) :: dist
432  real(RP), dimension(:), intent(in) :: uh
433  type(fespacexk) , intent(in) :: Y
434  type(quadmesh) , intent(in) :: qdm
435  procedure(r3tor) :: u_1, u_2
436  procedure(r3tor) , optional :: u_3
437 
438  real(RP), dimension(:), allocatable :: uh_i
439  real(RP) :: itg
440 
441  if (.NOT.valid(y)) call quit(&
442  &"feSpacexk_mod: L2_dist_2:&
443  & fespacsexk 'Y' not valid" )
444  if (.NOT.valid(qdm)) call quit(&
445  &"feSpacexk_mod: L2_dist_2:&
446  & quadMesh 'qdm' not valid" )
447  if(.NOT.associated( qdm%mesh, y%X_h%mesh)) call quit( &
448  &"feSpacexk_mod: L2_dist_2:&
449  & 'Y' and 'qdmh' associated to different meshes" )
450  if (size(uh,1) /= y%nbDof2) call quit( &
451  &"feSpacexk_mod: L2_dist_2:&
452  & wrong size for the finite element function 'uh'")
453  if ( y%k==1 ) call quit(&
454  &"feSpacexk_mod: L2_dist_2:&
455  & dim 1 case, use integral::L2_dist" )
456  if ( y%k>3 ) call quit(&
457  &"feSpacexk_mod: L2_dist_2: dim > 3 " )
458  if (( y%k==3 ) .AND. (.NOT.present(u_3))) call quit(&
459  &"feSpacexk_mod: L2_dist_2:&
460  & argument 'u_3' not provided " )
461 
462  if (choral_verb>1) write(*,*) &
463  & 'feSpacexk_mod : L2_dist_2'
464 
465  !! first component
466  call extract_component(uh_i, uh, y, 1)
467  itg = l2_dist(u_1, uh_i, y%X_h, qdm)
468  dist = itg**2
469 
470  !! second component
471  call extract_component(uh_i, uh, y, 2)
472  itg = l2_dist(u_2, uh_i, y%X_h, qdm)
473  dist = dist + itg**2
474 
475  !! third component
476  if (y%k==3) then
477  call extract_component(uh_i, uh, y, 3)
478  itg = l2_dist(u_3, uh_i, y%X_h, qdm)
479  dist = dist + itg**2
480  end if
481 
482  dist = sqrt(dist)
483 
484  end function l2_dist_2
485 
486 
512  function l2_dist_grad_2(uh, Y, qdm, grad_u1,&
513  & grad_u2, grad_u3) result(dist)
514  real(RP) :: dist
515  real(RP), dimension(:), intent(in) :: uh
516  type(fespacexk) , intent(in) :: Y
517  type(quadmesh) , intent(in) :: qdm
518  procedure(r3tor3) :: grad_u1, grad_u2
519  procedure(r3tor3) , optional :: grad_u3
520 
521  real(RP), dimension(:), allocatable :: uh_i
522  real(RP) :: itg
523 
524  if (.NOT.valid(y)) call quit(&
525  &"feSpacexk_mod: L2_dist_grad_2:&
526  & fespacsexk 'Y' not valid" )
527  if (.NOT.valid(qdm)) call quit(&
528  &"feSpacexk_mod: L2_dist_grad_2:&
529  & quadMesh 'qdm' not valid" )
530  if(.NOT.associated( qdm%mesh, y%X_h%mesh)) call quit( &
531  &"feSpacexk_mod: L2_dist_grad_2:&
532  & 'Y' and 'qdmh' associated to different meshes" )
533  if (size(uh,1) /= y%nbDof2) call quit( &
534  &"feSpacexk_mod: L2_dist_grad_2:&
535  & wrong size for the finite element function 'uh'")
536  if ( y%k==1 ) call quit(&
537  &"feSpacexk_mod: L2_dist_grad_2:&
538  & dim 1 case, use integral::L2_dist_grad" )
539  if ( y%k>3 ) call quit(&
540  &"feSpacexk_mod: L2_dist_grad_2: dim > 3 " )
541  if (( y%k==3 ) .AND. (.NOT.present(grad_u3))) call quit(&
542  &"feSpacexk_mod: L2_dist_grad_2:&
543  & argument 'grad_u3' not provided " )
544 
545  if (choral_verb>1) write(*,*) &
546  & 'feSpacexk_mod : L2_dist_grad_2'
547 
548  !! first component
549  call extract_component(uh_i, uh, y, 1)
550  itg = l2_dist_grad(grad_u1, uh_i, y%X_h, qdm)
551  dist = itg**2
552 
553  !! second component
554  call extract_component(uh_i, uh, y, 2)
555  itg = l2_dist_grad(grad_u2, uh_i, y%X_h, qdm)
556  dist = dist + itg**2
557 
558  !! third component
559  if (y%k==3) then
560  call extract_component(uh_i, uh, y, 3)
561  itg = l2_dist_grad(grad_u3, uh_i, y%X_h, qdm)
562  dist = dist + itg**2
563  end if
564 
565  dist = sqrt(dist)
566 
567  end function l2_dist_grad_2
568 
569 
570 
571 end module fespacexk_mod
set a graph row
Definition: graph_mod.F90:144
subroutine fespacexk_interp_vect_func(u_h, Y, u_1, u_2, u_3)
Interpolate a function u : R^3 –> R^d given by its components u_1, u_2 and u_3 on Y = [X_h]^d...
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
DERIVED TYPE graph: sparse matrices of boolean in CSR format
Definition: graph_mod.F90:33
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
subroutine, public quit(message)
Stop program execution, display an error messahe.
print a short description
The type feSpacexk defines for a finite element space.
DERIVED TYPE feSpacexk: define for a finite element space.
subroutine fespacexk_print(Y)
Print a short description.
integral L2 distance
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
DERIVED TYPE quadMesh: integration methods on meshes
subroutine, public extract_component(u_c, u, Y, c)
Extract the component of a finite element function .
DERIVED TYPE feSpace: finite element spaces
Definition: feSpace_mod.F90:58
type(fespacexk) function fespacexk_create(X_h, k)
Constructor for the type feSpacexk
logical function fespacexk_valid(Y)
Check if the structure content is correct.
real(rp) function l2_dist_2(uh, Y, qdm, u_1, u_2, u_3)
Returns .
DERIVED TYPE mesh
Definition: mesh_mod.F90:57
allocate memory for real(RP) arrays
Definition: real_type.F90:158
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 fespacexk_clear(Y)
Destructor for feSpacexk type.
subroutine, public interp_scal_func(uh, u, X_h)
Interpolation of a scalar function to a scalar finite element function.
real(rp) function l2_dist_grad_2(uh, Y, qdm, grad_u1, grad_u2, grad_u3)
Returns .
subroutine k(y, x)
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
The type quadMesh defines integration methods on meshes.
The type graph stores sparse matrices of boolean in CSR format.
Definition: graph_mod.F90:78
extract a graph row
Definition: graph_mod.F90:149