137 real(RP),
dimension(3, CELL_MAX_NBNODES) :: t
140 real(RP),
dimension(:,:) ,
allocatable :: y
143 real(RP),
dimension(:,:) ,
allocatable :: ty
146 real(RP),
dimension(:,:,:),
allocatable :: dty
149 real(RP),
dimension(:) ,
allocatable :: jy
152 real(RP),
dimension(:,:,:),
allocatable :: cy
217 integer ,
intent(in) :: dim, nn
218 real(RP),
dimension(dim,nn),
intent(in) :: y
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')
236 call allocmem(g%DTy, 3, g%dim, g%nn)
237 call allocmem(g%Cy , 3, g%dim, g%nn)
251 integer ,
intent(in) :: nNd, type
252 real(RP),
dimension(3,nNd),
intent(in) :: X
258 &
call quit(
'geoTsf_mod: geoTsf_assemble: wrong cell type')
260 &
call quit(
'geoTsf_mod: geoTsf_assemble: wrong number of nodes')
269 g%T(1:3,1) = x(1:3,1)
270 g%T(1:3,2) = x(1:3,2) - x(1:3,1)
273 g%Ty(1:3,ii) = g%T(1:3,1) + g%T(1:3,2)*g%y(1,ii)
277 g%T(1:3,1) = x(1:3,1)
279 g%T(1:3,2) = -3._rp*x(1:3,1) - 1._rp*x(1:3,2) &
282 g%T(1:3,3) = 2._rp*x(1:3,1) + 2._rp*x(1:3,2) &
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
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)
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)
301 g%T(1:3,1) = x(1:3,1)
304 g%T(1:3,2) = -3._rp*x(1:3,1) - x(1:3,2) &
308 g%T(1:3,3) = -3._rp*x(1:3,1) - x(1:3,3) &
312 g%T(1:3,4) = 2._rp*x(1:3,1) + 2._rp*x(1:3,2) &
316 g%T(1:3,5) = 2._rp*x(1:3,1) + 2._rp*x(1:3,3) &
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)
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)
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)
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)
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)
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)
375 call quit(
"geoTsf_mod: geoTsf_assemble: cell type not supported" )
399 s = abs(
det_3x3( g%T(1:3,2:4) ) )
405 g%Jy(ii) =
nrm2_3d(g%DTy(:,1,ii))
411 g%Jy(ii) =
det_3x2( g%DTy(:,:,ii) )
417 g%Jy(ii) = abs(
det_3x3( g%DTy(:,:,ii) ) )
421 call quit(
'geoTsf_mod: compute_J: cell type not supported' )
436 g%Cy(1:3,1,1) = g%T(1:3,2) / g%Jy(1)
440 g%Cy(1:3,1,ii) = g%Cy(1:3,1,1)
444 g%Jy(1) =
det_3x2( g%T(1:3,2:3) )
448 g%Cy(1:3,1:2,ii) = g%Cy(1:3,1:2,1)
452 g%Jy(1) = abs(
det_3x3( g%T(1:3,2:4) ) )
456 g%Cy(1:3,1:3,ii) = g%Cy(1:3,1:3,1)
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)
469 g%Jy(ii) =
det_3x2( g%DTy(1:3,1:2,ii) )
476 g%Jy(ii) = abs(
det_3x3( g%DTy(1:3,1:3,ii) ))
481 call quit(
'geoTsf_mod: compute_JC: cell type not supported' )
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)
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)
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) )
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)
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) ) )
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)
530 call quit(
'geoTsf_mod: compute_DTJC: cell type not supported' )
546 g%DTy(1:3,1,1) = g%T(1:3,2)
547 g%Jy(1) =
nrm2_3d( g%T(1:3,2) )
551 g%DTy(1:3,1,ii) = g%DTy(1:3,1,1)
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) )
560 g%DTy(1:3,1:2,ii) = g%DTy(1:3,1:2,1)
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) ) )
569 g%DTy(1:3,1:3,ii) = g%DTy(1:3,1:3,1)
575 g%Jy(ii) =
nrm2_3d( g%DTy(1:3,1,ii) )
581 g%Jy(ii) =
det_3x2( g%DTy(1:3,1:2,ii) )
587 g%Jy(ii) = abs(
det_3x3( g%DTy(1:3,1:3,ii) ) )
591 call quit(
'geoTsf_mod: compute_DTJ: cell type not supported' )
604 real(RP),
dimension(3) ,
intent(in) :: X
605 real(RP),
dimension(3,n),
intent(inout) :: Y
606 integer ,
intent(in) :: n, ct
608 real(RP),
dimension(3) :: z
613 &
call quit(
'geoTsf_mod: belongsToCell: wrong cell type')
615 &
call quit(
'geoTsf_mod: belongsToCell: wrong number of nodes')
626 if ( z(1) < -tol )
return 629 if ( z(3) > tol )
return 631 z = y(:,1) - z(1)*y(:,2)
632 b = ( maxval( abs( z ) ) < tol )
636 if ( (z(1) < -tol) )
return 638 if ( z(2) < -tol )
return 640 z(3) = z(1) + z(2) - 1.0_rp
641 if ( z(3)>tol )
return 644 z = y(:,1) - z(1)*y(:,2) - z(2)*y(:,3)
645 b = ( maxval( abs( z ) ) < tol )
649 if ( z(1) < -tol )
return 651 if ( z(2) < -tol )
return 653 if ( z(3) < -tol )
return 655 z(1) = sum(z) - 1.0_rp
659 call quit(
"geoTsf_mod: belongsToCell: cell type not supported" )
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
692 &
call quit(
'geoTsf_mod: geoTsf_T_Inv: wrong cell type')
694 &
call quit(
'geoTsf_mod: geoTsf_T_Inv: wrong number of nodes')
700 y(:,ii) = y(:,ii) - y(:,1)
707 z(1) = sum(y(:,1) * y(:,2)) / sum(y(:,2)**2)
716 call quit(
"geoTsf_mod : geoTsf_T_Inv: invalid cell type" )
724 real(RP),
dimension(3) :: v1, v
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
756 real(RP),
dimension(3) :: v1, v2
762 &
"geoTsf_mod: cell_ms_itf_n: wrong number of interfaces" )
786 v2 = g%T(1:3,3) - g%T(1:3,2)
794 n(:,ii) = n(:,ii) / itf_ms(ii)
798 call quit(
"geoTsf_mod: cell_ms_itf_n: cell type not supported" )
811 g%DTy(1:3,1, ii) = g%T(1:3,2) + 2._rp*g%T(1:3,3)*g%y(1,ii)
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)
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)
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)
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)
DERIVED TYPE geoTsf: geometrical transformation of reference cells
subroutine compute_dt_cell_edg_2(g)
integer, parameter cell_edg_2
Quadratic edge.
integer, parameter cell_tot_nb
Number of CELL types.
integer, parameter cell_tet
Tetrahedron.
real(rp) function, dimension(3, 2), public cplmtmat_3x2(M)
N = 'complementary' matrix of M.
subroutine geotsf_assemble(g, type, X, nNd)
Define the transformation T: K_ref –> K.
deallocate memory for real(RP) arrays
subroutine, public compute_dtj(g)
Compute DTy(:,:,i) and Jy(i) at the nodes y(:,i)
subroutine compute_dt_cell_tet_2(g)
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.
subroutine geotsf_clear(g)
Destructor for geoTsf type.
1- Define the transformation T: K_ref –> K i.e. compute the ytansformation 'T' coefficients ...
subroutine, public solve_3x3(x, A, b)
Solve A x = b for 3x3 matrix.
integer, parameter cell_edg
Edge (line segment)
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Geometrical transformation .
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...
subroutine compute_dt_cell_trg_2(g)
subroutine, public compute_j(g)
Compute Jy(i) at the nodes y(:;i)
integer, dimension(cell_tot_nb), public cell_nbnodes
Number of nodes for each cell type.
integer, parameter cell_trg
Triangle.
integer, dimension(cell_tot_nb), public cell_nbitf
Number of interfaces for each cell type.
integer, parameter cell_tet_2
Quadratic tetrahedron.
integer, parameter cell_trg_2
Quadratic triangle.
allocate memory for real(RP) arrays
subroutine, public compute_dtjc(g)
Compute DTy(:,:,i), Jy(i) and Cy(:,:,i) at the nodes y(:,i)
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)
subroutine, public compute_jc(g)
Computes Jy(i) and Cy(:,:,i) at the nodes y(:;i)
real(rp) function, dimension(3, 3), public cplmtmat_3x3(M)
N = complementary matrix of M = transpose of the cofactor matrix.
real(rp) function, public det_3x3(M)
determinant of 3x3 matrix
real(rp), parameter, public real_tol
type(geotsf) function geotsf_create(dim, nn, y)
Constructor for the type geoTsf
real(rp) function, dimension(3), public crossprod_3d(u, v)
Cross product in R**3.
real(rp) function, public nrm2_3d(v)
Euclidian norm.
DEFINITION OF GEOMETRICAL CELLS (for meshes)
real(rp) function, public det_3x2(M)
'determinant' of 3x2 matrix