Choral
geoTsf_mod.F90
Go to the documentation of this file.
1 
91 
92 module geotsf_mod
93 
94  use choral_constants
95  use real_type
96  use basic_tools
97  use algebra_lin
98  use cell_mod
99 
100  implicit none
101  private
102 
103  ! %----------------------------------------%
104  ! | |
105  ! | PUBLIC DATA |
106  ! | |
107  ! %----------------------------------------%
108  public :: geotsf ! TYPE
109 
110  public :: clear ! TESTED
111  public :: assemble ! TESTED
112  public :: compute_j ! TESTED
113  public :: compute_jc ! TESTED
114  public :: compute_dtj ! TESTED
115  public :: compute_dtjc ! TESTED
116  public :: belongstocell ! TESTED
117  public :: geotsf_t_inv ! TESTED
118 
119  public :: cell_ms_itf_n
120 
121  ! %----------------------------------------%
122  ! | |
123  ! | DERIVED TYPE |
124  ! | |
125  ! %----------------------------------------%
130  type geotsf
131 
132  integer :: type
133  integer :: dim
134  integer :: nn
135 
137  real(RP), dimension(3, CELL_MAX_NBNODES) :: t
138 
140  real(RP), dimension(:,:) , allocatable :: y
141 
143  real(RP), dimension(:,:) , allocatable :: ty
144 
146  real(RP), dimension(:,:,:), allocatable :: dty
147 
149  real(RP), dimension(:) , allocatable :: jy
150 
152  real(RP), dimension(:,:,:), allocatable :: cy
153 
154  contains
155 
157  final :: geotsf_clear
158 
159  end type geotsf
160 
161  ! %----------------------------------------%
162  ! | |
163  ! | GENERIc SUBROUTINES |
164  ! | |
165  ! %----------------------------------------%
167  interface clear
168  module procedure geotsf_clear
169  end interface clear
170 
171  !! Constructor
172  interface geotsf
173  module procedure geotsf_create
174  end interface geotsf
175 
180  interface assemble
181  module procedure geotsf_assemble
182  end interface assemble
183 
184 contains
185 
187  subroutine geotsf_clear(g)
188  type(geotsf) , intent(inout) :: g
189 
190  call freemem(g%y)
191  call freemem(g%Ty)
192  call freemem(g%DTy)
193  call freemem(g%Jy)
194  call freemem(g%Cy)
195 
196  end subroutine geotsf_clear
197 
198 
215  function geotsf_create(dim, nn, y) result(g)
216  type(geotsf) :: g
217  integer , intent(in) :: dim, nn
218  real(RP), dimension(dim,nn), intent(in) :: y
219 
220 #if DBG>0
221  if ( (dim > 3).OR. (dim < 0) ) call quit( &
222  & 'geoTsf_mod: geoTsf_create: wrong dimension')
223  if ( nn < 0 ) call quit( &
224  & 'geoTsf_mod: geoTsf_create: wrnong number of points')
225 #endif
226 
227  call clear(g)
228 
229  g%dim = dim
230  g%nn = nn
231  call allocmem(g%y,g%dim, g%nn)
232  g%y = y
233 
234  call allocmem(g%Ty , 3, g%nn)
235  call allocmem(g%Jy , g%nn)
236  call allocmem(g%DTy, 3, g%dim, g%nn)
237  call allocmem(g%Cy , 3, g%dim, g%nn)
238 
239  end function geotsf_create
240 
249  subroutine geotsf_assemble(g, type, X, nNd)
250  type(geotsf) , intent(inout):: g
251  integer , intent(in) :: nNd, type
252  real(RP), dimension(3,nNd), intent(in) :: X
253 
254  integer :: ii
255 
256 #if DBG>0
257  if ( (type < 0).OR.(type>cell_tot_nb) ) &
258  & call quit('geoTsf_mod: geoTsf_assemble: wrong cell type')
259  if ( nnd /= cell_nbnodes(type) ) &
260  & call quit('geoTsf_mod: geoTsf_assemble: wrong number of nodes')
261 #endif
262 
263  g%type = type
264 
265  ! define the transformation
266  select case(g%type)
267 
268  case(cell_edg)
269  g%T(1:3,1) = x(1:3,1)
270  g%T(1:3,2) = x(1:3,2) - x(1:3,1)
271 
272  do ii=1, g%nn
273  g%Ty(1:3,ii) = g%T(1:3,1) + g%T(1:3,2)*g%y(1,ii)
274  end do
275 
276  case(cell_edg_2)
277  g%T(1:3,1) = x(1:3,1)
278 
279  g%T(1:3,2) = -3._rp*x(1:3,1) - 1._rp*x(1:3,2) &
280  & + 4._rp*x(1:3,3)
281 
282  g%T(1:3,3) = 2._rp*x(1:3,1) + 2._rp*x(1:3,2) &
283  & - 4._rp*x(1:3,3)
284 
285  do ii=1, g%nn
286  g%Ty(1:3,ii) = g%T(1:3,1) + g%T(1:3,2)*g%y(1,ii) &
287  & + g%T(1:3,3)*g%y(1,ii)**2
288  end do
289 
290  case(cell_trg)
291  g%T(1:3,1) = x(1:3,1)
292  g%T(1:3,2) = x(1:3,2) - x(1:3,1)
293  g%T(1:3,3) = x(1:3,3) - x(1:3,1)
294 
295  do ii=1, g%nn
296  g%Ty(1:3,ii) = g%T(1:3,1) + g%T(1:3,2)*g%y(1,ii) &
297  & + g%T(1:3,3)*g%y(2,ii)
298  end do
299 
300  case(cell_trg_2)
301  g%T(1:3,1) = x(1:3,1)
302 
303  ! coef x
304  g%T(1:3,2) = -3._rp*x(1:3,1) - x(1:3,2) &
305  & + 4._rp*x(1:3,4)
306 
307  ! coef y
308  g%T(1:3,3) = -3._rp*x(1:3,1) - x(1:3,3) &
309  & + 4._rp*x(1:3,6)
310 
311  ! coef x**2
312  g%T(1:3,4) = 2._rp*x(1:3,1) + 2._rp*x(1:3,2) &
313  & - 4._rp*x(1:3,4)
314 
315  ! coef y**2
316  g%T(1:3,5) = 2._rp*x(1:3,1) + 2._rp*x(1:3,3) &
317  & - 4._rp*x(1:3,6)
318 
319  ! coef x*y
320  g%T(1:3,6) = 4._rp*x(1:3,5) - 4._rp*g%T(1:3,1) &
321  & - 2._rp*g%T(1:3,2) - 2._rp*g%T(1:3,3) &
322  & - g%T(1:3,4) - g%T(1:3,5)
323 
324  do ii=1, g%nn
325  g%Ty(1:3,ii) = g%T(1:3,1) &
326  & + g%T(1:3,2)*g%y(1,ii) &
327  & + g%T(1:3,3)*g%y(2,ii) &
328  & + g%T(1:3,4)*g%y(1,ii)**2 &
329  & + g%T(1:3,5)*g%y(2,ii)**2 &
330  & + g%T(1:3,6)*g%y(1,ii)*g%y(2,ii)
331  end do
332 
333  case(cell_tet)
334  g%T(1:3,1) = x(1:3,1)
335  g%T(1:3,2) = x(1:3,2) - x(1:3,1)
336  g%T(1:3,3) = x(1:3,3) - x(1:3,1)
337  g%T(1:3,4) = x(1:3,4) - x(1:3,1)
338 
339  do ii=1, g%nn
340  g%Ty(1:3,ii) = g%T(1:3,1) &
341  & + g%T(1:3,2)*g%y(1,ii) &
342  & + g%T(1:3,3)*g%y(2,ii) &
343  & + g%T(1:3,4)*g%y(3,ii)
344  end do
345 
346  case(cell_tet_2)
347 
348  g%T(1:3,1 ) = x(1:3,1)
349  g%T(1:3,2 ) = -3._rp*x(1:3,1) - x(1:3,2) + 4._rp*x(1:3,5)
350  g%T(1:3,3 ) = -3._rp*x(1:3,1) - x(1:3,3) + 4._rp*x(1:3,7)
351  g%T(1:3,4 ) = -3._rp*x(1:3,1) - x(1:3,4) + 4._rp*x(1:3,8)
352  g%T(1:3,5 ) = 2._rp*x(1:3,1) + 2._rp*x(1:3,2) - 4._rp*x(1:3,5)
353  g%T(1:3,6 ) = 2._rp*x(1:3,1) + 2._rp*x(1:3,3) - 4._rp*x(1:3,7)
354  g%T(1:3,7 ) = 2._rp*x(1:3,1) + 2._rp*x(1:3,4) - 4._rp*x(1:3,8)
355  g%T(1:3,8 ) = 4._rp*x(1:3,1) - 4._rp*x(1:3,7) - 4._rp*x(1:3,8) + 4._rp*x(1:3,9)
356  g%T(1:3,9 ) = 4._rp*x(1:3,1) - 4._rp*x(1:3,5) - 4._rp*x(1:3,8) + 4._rp*x(1:3,10)
357  g%T(1:3,10) = 4._rp*x(1:3,1) - 4._rp*x(1:3,5) - 4._rp*x(1:3,7) + 4._rp*x(1:3,6)
358 
359  do ii=1, g%nn
360  g%Ty(1:3,ii) = g%T(1:3,1) &
361  & + g%T(1:3,2) * g%y(1,ii) &
362  & + g%T(1:3,3) * g%y(2,ii) &
363  & + g%T(1:3,4) * g%y(3,ii) &
364  & + g%T(1:3,5) * g%y(1,ii)**2 &
365  & + g%T(1:3,6) * g%y(2,ii)**2 &
366  & + g%T(1:3,7) * g%y(3,ii)**2 &
367  & + g%T(1:3,8) * g%y(2,ii)*g%y(3,ii) &
368  & + g%T(1:3,9) * g%y(3,ii)*g%y(1,ii) &
369  & + g%T(1:3,10)* g%y(1,ii)*g%y(2,ii)
370 
371  end do
372 
373 
374  case default
375  call quit( "geoTsf_mod: geoTsf_assemble: cell type not supported" )
376 
377  end select
378 
379  end subroutine geotsf_assemble
380 
382  subroutine compute_j(g)
383  type(geotsf), intent(inout) :: g
384 
385  integer :: ii
386  real(RP) :: s
387 
388  select case(g%Type)
389 
390  case(cell_edg)
391  s = nrm2_3d(g%T(1:3,2))
392  g%Jy = s
393 
394  case(cell_trg)
395  s = det_3x2( g%T(1:3,2:3) )
396  g%Jy = s
397 
398  case(cell_tet)
399  s = abs( det_3x3( g%T(1:3,2:4) ) )
400  g%Jy = s
401 
402  case(cell_edg_2)
403  call compute_dt_cell_edg_2(g)
404  do ii=1, g%nn
405  g%Jy(ii) = nrm2_3d(g%DTy(:,1,ii))
406  end do
407 
408  case(cell_trg_2)
409  call compute_dt_cell_trg_2(g)
410  do ii=1, g%nn
411  g%Jy(ii) = det_3x2( g%DTy(:,:,ii) )
412  end do
413 
414  case(cell_tet_2)
415  call compute_dt_cell_tet_2(g)
416  do ii=1, g%nn
417  g%Jy(ii) = abs( det_3x3( g%DTy(:,:,ii) ) )
418  end do
419 
420  case default
421  call quit( 'geoTsf_mod: compute_J: cell type not supported' )
422  end select
423 
424  end subroutine compute_j
425 
427  subroutine compute_jc(g)
428  type(geotsf), intent(inout) :: g
429 
430  integer :: ii
431 
432  select case(g%Type)
433 
434  case(cell_edg)
435  g%Jy(1) = nrm2_3d(g%T(1:3,2))
436  g%Cy(1:3,1,1) = g%T(1:3,2) / g%Jy(1)
437 
438  do ii=2, g%nn
439  g%Jy(ii) = g%Jy(1)
440  g%Cy(1:3,1,ii) = g%Cy(1:3,1,1)
441  end do
442 
443  case(cell_trg)
444  g%Jy(1) = det_3x2( g%T(1:3,2:3) )
445  g%Cy(1:3,1:2,1) = cplmtmat_3x2( g%T(1:3,2:3) )
446  do ii=2, g%nn
447  g%Jy(ii) = g%Jy(1)
448  g%Cy(1:3,1:2,ii) = g%Cy(1:3,1:2,1)
449  end do
450 
451  case(cell_tet)
452  g%Jy(1) = abs( det_3x3( g%T(1:3,2:4) ) )
453  g%Cy(1:3,1:3,1) = cplmtmat_3x3( g%T(1:3,2:4) )
454  do ii=2, g%nn
455  g%Jy(ii) = g%Jy(1)
456  g%Cy(1:3,1:3,ii) = g%Cy(1:3,1:3,1)
457  end do
458 
459  case(cell_edg_2)
460  call compute_dt_cell_edg_2(g)
461  do ii=1, g%nn
462  g%Jy(ii) = nrm2_3d(g%DTy(1:3,1,ii))
463  g%Cy(1:3,1,ii) = g%DTy(1:3,1,ii) / g%Jy(ii)
464  end do
465 
466  case(cell_trg_2)
467  call compute_dt_cell_trg_2(g)
468  do ii=1, g%nn
469  g%Jy(ii) = det_3x2( g%DTy(1:3,1:2,ii) )
470  g%Cy(1:3,1:2,ii) = cplmtmat_3x2( g%DTy(1:3,1:2,ii) )
471  end do
472 
473  case(cell_tet_2)
474  call compute_dt_cell_tet_2(g)
475  do ii=1, g%nn
476  g%Jy(ii) = abs(det_3x3( g%DTy(1:3,1:3,ii) ))
477  g%Cy(1:3,1:3,ii) = cplmtmat_3x3( g%DTy(1:3,1:3,ii) )
478  end do
479 
480  case default
481  call quit( 'geoTsf_mod: compute_JC: cell type not supported' )
482  end select
483 
484  end subroutine compute_jc
485 
486 
489  subroutine compute_dtjc(g)
490  type(geotsf), intent(inout) :: g
491 
492  integer :: ii
493 
494  select case(g%Type)
495 
496  case(cell_edg)
497  g%DTy(1:3,1,1) = g%T(1:3,2)
498  g%Jy(1) = nrm2_3d( g%T(1:3,2) )
499  g%Cy(1:3,1,1) = g%T(1:3,2) / g%Jy(1)
500  do ii=2, g%nn
501  g%Jy(ii) = g%Jy(1)
502  g%Cy( 1:3,1,ii) = g%Cy( 1:3,1,1)
503  g%DTy(1:3,1,ii) = g%DTy(1:3,1,1)
504  end do
505 
506  case(cell_trg)
507  g%DTy(1:3,1:2,1) = g%T(1:3,2:3)
508  g%Jy(1) = det_3x2( g%T(1:3,2:3) )
509  g%Cy( 1:3,1:2,1) = cplmtmat_3x2( g%T(1:3,2:3) )
510  do ii=2, g%nn
511  g%Jy(ii) = g%Jy(1)
512  g%Cy( 1:3,1:2,ii) = g%Cy( 1:3,1:2,1)
513  g%DTy(1:3,1:2,ii) = g%DTy(1:3,1:2,1)
514  end do
515 
516  case(cell_tet)
517  g%DTy(1:3,1:3,1) = g%T(1:3,2:4)
518  g%Jy(1) = abs( det_3x3( g%T(1:3,2:4) ) )
519  g%Cy( 1:3,1:3,1) = cplmtmat_3x3( g%T(1:3,2:4) )
520  do ii=2, g%nn
521  g%Jy(ii) = g%Jy(1)
522  g%Cy( 1:3,1:3,ii) = g%Cy( 1:3,1:3,1)
523  g%DTy(1:3,1:3,ii) = g%DTy(1:3,1:3,1)
524  end do
525 
526  case(cell_edg_2, cell_trg_2, cell_tet_2)
527  call compute_jc(g)
528 
529  case default
530  call quit( 'geoTsf_mod: compute_DTJC: cell type not supported' )
531  end select
532 
533  end subroutine compute_dtjc
534 
535 
536 
538  subroutine compute_dtj(g)
539  type(geotsf), intent(inout) :: g
540 
541  integer :: ii
542 
543  select case(g%Type)
544 
545  case(cell_edg)
546  g%DTy(1:3,1,1) = g%T(1:3,2)
547  g%Jy(1) = nrm2_3d( g%T(1:3,2) )
548 
549  do ii=2, g%nn
550  g%Jy(ii) = g%Jy(1)
551  g%DTy(1:3,1,ii) = g%DTy(1:3,1,1)
552  end do
553 
554  case(cell_trg)
555  g%DTy(1:3,1:2,1) = g%T(1:3,2:3)
556  g%Jy(1) = det_3x2( g%T(1:3,2:3) )
557 
558  do ii=2, g%nn
559  g%Jy(ii) = g%Jy(1)
560  g%DTy(1:3,1:2,ii) = g%DTy(1:3,1:2,1)
561  end do
562 
563  case(cell_tet)
564  g%DTy(1:3,1:3,1) = g%T(1:3,2:4)
565  g%Jy(1) = abs( det_3x3( g%T(1:3,2:4) ) )
566 
567  do ii=2, g%nn
568  g%Jy(ii) = g%Jy(1)
569  g%DTy(1:3,1:3,ii) = g%DTy(1:3,1:3,1)
570  end do
571 
572  case(cell_edg_2)
573  call compute_dt_cell_edg_2(g)
574  do ii=1, g%nn
575  g%Jy(ii) = nrm2_3d( g%DTy(1:3,1,ii) )
576  end do
577 
578  case(cell_trg_2)
579  call compute_dt_cell_trg_2(g)
580  do ii=1, g%nn
581  g%Jy(ii) = det_3x2( g%DTy(1:3,1:2,ii) )
582  end do
583 
584  case(cell_tet_2)
585  call compute_dt_cell_tet_2(g)
586  do ii=1, g%nn
587  g%Jy(ii) = abs( det_3x3( g%DTy(1:3,1:3,ii) ) )
588  end do
589 
590  case default
591  call quit( 'geoTsf_mod: compute_DTJ: cell type not supported' )
592  end select
593 
594  end subroutine compute_dtj
595 
596 
602  function belongstocell(x, Y, n, ct) result(b)
603  logical :: b
604  real(RP), dimension(3) , intent(in) :: X
605  real(RP), dimension(3,n), intent(inout) :: Y
606  integer , intent(in) :: n, ct
607 
608  real(RP), dimension(3) :: z
609  real(RP) :: tol
610 
611 #if DBG>0
612  if ( (ct < 0).OR.(ct>cell_tot_nb) ) &
613  & call quit('geoTsf_mod: belongsToCell: wrong cell type')
614  if ( n /= cell_nbnodes(ct) ) &
615  & call quit('geoTsf_mod: belongsToCell: wrong number of nodes')
616 #endif
617 
618  call geotsf_t_inv(z, x, y, n, ct)
619  tol = real_tol * 1e2_rp
620 
621  b = .false.
622  select case(ct)
623 
624  case(cell_edg)
625 
626  if ( z(1) < -tol ) return
627 
628  z(3) = z(1) - 1.0_rp
629  if ( z(3) > tol ) return
630 
631  z = y(:,1) - z(1)*y(:,2) ! AX - z(1)*AB
632  b = ( maxval( abs( z ) ) < tol )
633 
634  case(cell_trg)
635 
636  if ( (z(1) < -tol) ) return
637 
638  if ( z(2) < -tol ) return
639 
640  z(3) = z(1) + z(2) - 1.0_rp
641  if ( z(3)>tol ) return
642 
643  ! AX - z(1)*AB - z(2)*AC
644  z = y(:,1) - z(1)*y(:,2) - z(2)*y(:,3)
645  b = ( maxval( abs( z ) ) < tol )
646 
647  case(cell_tet)
648 
649  if ( z(1) < -tol ) return
650 
651  if ( z(2) < -tol ) return
652 
653  if ( z(3) < -tol ) return
654 
655  z(1) = sum(z) - 1.0_rp
656  b = ( z(1) < tol )
657 
658  case default
659  call quit( "geoTsf_mod: belongsToCell: cell type not supported" )
660 
661  end select
662 
663  end function belongstocell
664 
665 
682  subroutine geotsf_t_inv(z, X, Y, n, ct)
683  real(RP), dimension(3) , intent(out) :: z
684  real(RP), dimension(3) , intent(in) :: X
685  real(RP), dimension(3,n), intent(inout) :: Y
686  integer , intent(in) :: n, ct
687 
688  integer :: ii
689 
690 #if DBG>0
691  if ( (ct < 0).OR.(ct>cell_tot_nb) ) &
692  & call quit('geoTsf_mod: geoTsf_T_Inv: wrong cell type')
693  if ( n /= cell_nbnodes(ct) ) &
694  & call quit('geoTsf_mod: geoTsf_T_Inv: wrong number of nodes')
695 #endif
696 
697  z = 0.0_rp
698 
699  do ii=2, n
700  y(:,ii) = y(:,ii) - y(:,1)
701  end do
702  y(:,1) = x - y(:,1)
703 
704  select case(ct)
705 
706  case(cell_edg)
707  z(1) = sum(y(:,1) * y(:,2)) / sum(y(:,2)**2)
708 
709  case(cell_trg)
710  call t_inv_trg()
711 
712  case(cell_tet)
713  call solve_3x3(z, y(:,2:4), y(:,1))
714 
715  case default
716  call quit( "geoTsf_mod : geoTsf_T_Inv: invalid cell type" )
717 
718  end select
719 
720  contains
721 
722  subroutine t_inv_trg()
723 
724  real(RP), dimension(3) :: v1, v
725  real(RP) :: a
726 
727  v1 = crossprod_3d(y(:,2), y(:,3) ) ! AB^AC
728  a = sum(v1**2) ! AB^AC . AB^AC
729 
730  v = crossprod_3d(y(:,1), y(:,3) ) ! AX^AC
731  z(1) = sum(v*v1) / a
732 
733  v = crossprod_3d(y(:,1), y(:,2) ) ! AX^AB
734  z(2) =-sum(v*v1) / a
735 
736  end subroutine t_inv_trg
737 
738  end subroutine geotsf_t_inv
739 
740 
748  subroutine cell_ms_itf_n(ms, n, itf_ms, nbItf, g)
749  type(geotsf) , intent(in) :: g
750  integer , intent(in) :: nbItf
751  real(RP) , intent(out) :: ms
752  real(RP), dimension(3,nbItf), intent(out) :: n
753  real(RP), dimension(nbItf) , intent(out) :: itf_ms
754 
755 
756  real(RP), dimension(3) :: v1, v2
757  real(RP) :: J
758  integer :: ii
759 
760 #if DBG>0
761  if ( nbitf /= cell_nbitf(g%type)) call quit( &
762  & "geoTsf_mod: cell_ms_itf_n: wrong number of interfaces" )
763 #endif
764 
765  select case(g%type)
766 
767  case(cell_edg)
768 
769  v1 = g%T(1:3,2)
770  ms = nrm2_3d(v1)
771 
772  v1 = v1 / ms
773 
774  n(:,1) = -v1
775  n(:,2) = v1
776 
777  itf_ms(1:2) = 1.0_rp
778 
779  case(cell_trg)
780  v1 = crossprod_3d( g%T(1:3,2), g%T(1:3,3) )
781 
782  j = nrm2_3d(v1)
783  ms = j * 0.5_rp
784  v1 = v1 / j
785 
786  v2 = g%T(1:3,3) - g%T(1:3,2)
787 
788  n(:,1) = crossprod_3d( g%T(1:3,2), v1 )
789  n(:,2) = crossprod_3d( v2 , v1 )
790  n(:,3) = crossprod_3d(-g%T(1:3,3), v1 )
791 
792  do ii=1, 3
793  itf_ms(ii) = nrm2_3d(n(:,ii))
794  n(:,ii) = n(:,ii) / itf_ms(ii)
795  end do
796 
797  case default
798  call quit( "geoTsf_mod: cell_ms_itf_n: cell type not supported" )
799 
800  end select
801 
802  end subroutine cell_ms_itf_n
803 
804 
805  subroutine compute_dt_cell_edg_2(g)
806  type(geotsf) , intent(inout) :: g
807 
808  integer :: ii
809 
810  do ii=1, g%nn
811  g%DTy(1:3,1, ii) = g%T(1:3,2) + 2._rp*g%T(1:3,3)*g%y(1,ii)
812  end do
813 
814  end subroutine compute_dt_cell_edg_2
815 
816 
817  subroutine compute_dt_cell_trg_2(g)
818  type(geotsf) , intent(inout) :: g
819 
820  integer :: ii
821 
822  do ii=1, g%nn
823  g%DTy(1:3,1,ii) = g%T(1:3,2) + 2._rp*g%T(1:3,4)*g%y(1,ii) &
824  & + g%T(1:3,6)*g%y(2,ii)
825  g%DTy(1:3,2,ii) = g%T(1:3,3) + 2._rp*g%T(1:3,5)*g%y(2,ii) &
826  & + g%T(1:3,6)*g%y(1,ii)
827  end do
828 
829  end subroutine compute_dt_cell_trg_2
830 
831 
832  subroutine compute_dt_cell_tet_2(g)
833  type(geotsf) , intent(inout) :: g
834 
835  integer :: ii
836 
837  do ii=1, g%nn
838  g%DTy(1:3,1,ii) = g%T(1:3,2) + 2._rp*g%T(1:3,5)*g%y(1,ii) &
839  & + g%T(1:3,9) * g%y(3,ii) &
840  & + g%T(1:3,10)* g%y(2,ii)
841 
842  g%DTy(1:3,2,ii) = g%T(1:3,3) + 2._rp*g%T(1:3,6)*g%y(2,ii) &
843  & + g%T(1:3,8) * g%y(3,ii) &
844  & + g%T(1:3,10)* g%y(1,ii)
845 
846  g%DTy(1:3,3,ii) = g%T(1:3,4) + 2._rp*g%T(1:3,7)*g%y(3,ii) &
847  & + g%T(1:3,8) * g%y(2,ii) &
848  & + g%T(1:3,9) * g%y(1,ii)
849  end do
850 
851  end subroutine compute_dt_cell_tet_2
852 
853 
854 end module geotsf_mod
855 
geotsf_mod
DERIVED TYPE geoTsf: geometrical transformation of reference cells
Definition: geoTsf_mod.F90:92
geotsf_mod::compute_dt_cell_edg_2
subroutine compute_dt_cell_edg_2(g)
Definition: geoTsf_mod.F90:806
choral_constants::cell_edg_2
integer, parameter cell_edg_2
Quadratic edge.
Definition: choral_constants.f90:33
algebra_lin
LINEAR ALGEBRA TOOLS
Definition: algebra_lin.F90:9
choral_constants::cell_tot_nb
integer, parameter cell_tot_nb
Number of CELL types.
Definition: choral_constants.f90:40
choral_constants::cell_tet
integer, parameter cell_tet
Tetrahedron.
Definition: choral_constants.f90:36
algebra_lin::cplmtmat_3x2
real(rp) function, dimension(3, 2), public cplmtmat_3x2(M)
N = &#39;complementary&#39; matrix of M.
Definition: algebra_lin.F90:247
basic_tools
BASIC TOOLS
Definition: basic_tools.F90:8
geotsf_mod::geotsf_assemble
subroutine geotsf_assemble(g, type, X, nNd)
Define the transformation T: K_ref –> K.
Definition: geoTsf_mod.F90:250
real_type::freemem
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
geotsf_mod::compute_dtj
subroutine, public compute_dtj(g)
Compute DTy(:,:,i) and Jy(i) at the nodes y(:,i)
Definition: geoTsf_mod.F90:539
geotsf_mod::compute_dt_cell_tet_2
subroutine compute_dt_cell_tet_2(g)
Definition: geoTsf_mod.F90:833
geotsf_mod::belongstocell
logical function, public belongstocell(x, Y, n, ct)
Tests if the node x belongs to the cell with node coordinates Y and of type ct.
Definition: geoTsf_mod.F90:603
geotsf_mod::clear
Destructor.
Definition: geoTsf_mod.F90:167
geotsf_mod::geotsf_clear
subroutine geotsf_clear(g)
Destructor for geoTsf type.
Definition: geoTsf_mod.F90:188
basic_tools::quit
subroutine, public quit(message)
Stop program execution, display an error messahe.
Definition: basic_tools.F90:106
geotsf_mod::assemble
1- Define the transformation T: K_ref –> K i.e. compute the ytansformation &#39;T&#39; coefficients ...
Definition: geoTsf_mod.F90:180
algebra_lin::solve_3x3
subroutine, public solve_3x3(x, A, b)
Solve A x = b for 3x3 matrix.
Definition: algebra_lin.F90:322
choral_constants::cell_edg
integer, parameter cell_edg
Edge (line segment)
Definition: choral_constants.f90:32
real_type
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
geotsf_mod::geotsf
Geometrical transformation .
Definition: geoTsf_mod.F90:130
geotsf_mod::cell_ms_itf_n
subroutine, public cell_ms_itf_n(ms, n, itf_ms, nbItf, g)
INPUT :: g = geoTsf, must have been set nbItf = number of interfaces OUTPUT : ms = measure of T(K_ref...
Definition: geoTsf_mod.F90:749
geotsf_mod::compute_dt_cell_trg_2
subroutine compute_dt_cell_trg_2(g)
Definition: geoTsf_mod.F90:818
geotsf_mod::compute_j
subroutine, public compute_j(g)
Compute Jy(i) at the nodes y(:;i)
Definition: geoTsf_mod.F90:383
t_inv_trg
subroutine t_inv_trg()
Definition: geoTsf_mod.F90:723
cell_mod::cell_nbnodes
integer, dimension(cell_tot_nb), public cell_nbnodes
Number of nodes for each cell type.
Definition: cell_mod.f90:115
choral_constants::cell_trg
integer, parameter cell_trg
Triangle.
Definition: choral_constants.f90:34
cell_mod::cell_nbitf
integer, dimension(cell_tot_nb), public cell_nbitf
Number of interfaces for each cell type.
Definition: cell_mod.f90:127
choral_constants::cell_tet_2
integer, parameter cell_tet_2
Quadratic tetrahedron.
Definition: choral_constants.f90:37
choral_constants
CHORAL CONSTANTS
Definition: choral_constants.f90:13
choral_constants::cell_trg_2
integer, parameter cell_trg_2
Quadratic triangle.
Definition: choral_constants.f90:35
real_type::allocmem
allocate memory for real(RP) arrays
Definition: real_type.F90:158
geotsf_mod::compute_dtjc
subroutine, public compute_dtjc(g)
Compute DTy(:,:,i), Jy(i) and Cy(:,:,i) at the nodes y(:,i)
Definition: geoTsf_mod.F90:490
geotsf_mod::geotsf_t_inv
subroutine, public geotsf_t_inv(z, X, Y, n, ct)
x = 3D node 3D cell of type ct with nodes = A, B, C, ... given by Y(:,1)=A, Y(:,2)=B, ..., Y(:,n)
Definition: geoTsf_mod.F90:683
geotsf_mod::compute_jc
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
Definition: geoTsf_mod.F90:428
algebra_lin::cplmtmat_3x3
real(rp) function, dimension(3, 3), public cplmtmat_3x3(M)
N = complementary matrix of M = transpose of the cofactor matrix.
Definition: algebra_lin.F90:233
algebra_lin::det_3x3
real(rp) function, public det_3x3(M)
determinant of 3x3 matrix
Definition: algebra_lin.F90:191
real_type::real_tol
real(rp), parameter, public real_tol
Definition: real_type.F90:91
geotsf_mod::geotsf_create
type(geotsf) function geotsf_create(dim, nn, y)
Constructor for the type geoTsf
Definition: geoTsf_mod.F90:216
algebra_lin::crossprod_3d
real(rp) function, dimension(3), public crossprod_3d(u, v)
Cross product in R**3.
Definition: algebra_lin.F90:178
algebra_lin::nrm2_3d
real(rp) function, public nrm2_3d(v)
Euclidian norm.
Definition: algebra_lin.F90:271
cell_mod
DEFINITION OF GEOMETRICAL CELLS (for meshes)
Definition: cell_mod.f90:82
algebra_lin::det_3x2
real(rp) function, public det_3x2(M)
&#39;determinant&#39; of 3x2 matrix
Definition: algebra_lin.F90:210
Generated on Mon May 17 2021 15:14:13 for Choral by   doxygen 1.8.13