Isoparametric subdivision of one element's geometry node coordinates into its four children. The parent geometry is the degree-nGeo Lagrange interpolant through parentCoords (sampled at geomInterp's control points, i.e. the mesh geometry nodes); each child node coordinate is that interpolant evaluated at the corresponding point of the parent reference square. Exact for any polynomial geometry of degree <= nGeo and for any control-node type, so straight-sided and curved (isoparametric) elements are both handled without approximation.
childCoords(d,i,j,c) = sum_{ii,jj} H(ii,i,ax) H(jj,j,ay) parentCoords(d,ii,jj)
where H(ii,i,a) = l_ii( halfmap(s_i,a) ), halfmap(s,a) = (s + (2a-1))/2 maps the child node reference coordinate s_i into the parent half [a-1,a], and l_ii is the parent's Lagrange basis. d = 1,2 are the physical coordinate components.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(Lagrange), | intent(in) | :: | geomInterp | |||
| integer, | intent(in) | :: | nGeo | |||
| real(kind=prec), | intent(in) | :: | parentCoords(1:2,1:nGeo+1,1:nGeo+1) | |||
| real(kind=prec), | intent(out) | :: | childCoords(1:2,1:nGeo+1,1:nGeo+1,1:4) |
subroutine SubdivideNodeCoords(geomInterp,nGeo,parentCoords,childCoords)
!! Isoparametric subdivision of one element's geometry node coordinates into its four
!! children. The parent geometry is the degree-nGeo Lagrange interpolant through
!! parentCoords (sampled at geomInterp's control points, i.e. the mesh geometry nodes); each
!! child node coordinate is that interpolant evaluated at the corresponding point of the
!! parent reference square. Exact for any polynomial geometry of degree <= nGeo and for any
!! control-node type, so straight-sided and curved (isoparametric) elements are both handled
!! without approximation.
!!
!! childCoords(d,i,j,c) = sum_{ii,jj} H(ii,i,ax) H(jj,j,ay) parentCoords(d,ii,jj)
!!
!! where H(ii,i,a) = l_ii( halfmap(s_i,a) ), halfmap(s,a) = (s + (2a-1))/2 maps the child
!! node reference coordinate s_i into the parent half [a-1,a], and l_ii is the parent's
!! Lagrange basis. d = 1,2 are the physical coordinate components.
implicit none
type(Lagrange),intent(in) :: geomInterp
integer,intent(in) :: nGeo
real(prec),intent(in) :: parentCoords(1:2,1:nGeo+1,1:nGeo+1)
real(prec),intent(out) :: childCoords(1:2,1:nGeo+1,1:nGeo+1,1:4)
! Local
integer :: i,ii,j,jj,c,a,ax,ay
integer,parameter :: axc(1:4) = [0,1,1,0]
integer,parameter :: ayc(1:4) = [0,0,1,1]
real(prec) :: s
real(prec) :: H(1:nGeo+1,1:nGeo+1,0:1)
real(prec) :: rowsum(1:2),acc(1:2)
! Precompute the half-interval interpolation weights H(ii,i,a) for the lower/left (a=0) and
! upper/right (a=1) halves. H(:,i,a) are the parent basis values at the mapped position of
! child node i.
do a = 0,1
do i = 1,nGeo+1
s = geomInterp%controlPoints(i)
H(1:nGeo+1,i,a) = geomInterp%CalculateLagrangePolynomials(0.5_prec*(s+real(2*a-1,prec)))
enddo
enddo
do c = 1,4
ax = axc(c)
ay = ayc(c)
do j = 1,nGeo+1
do i = 1,nGeo+1
acc(1:2) = 0.0_prec
do jj = 1,nGeo+1
rowsum(1:2) = 0.0_prec
do ii = 1,nGeo+1
rowsum(1:2) = rowsum(1:2)+H(ii,i,ax)*parentCoords(1:2,ii,jj)
enddo
acc(1:2) = acc(1:2)+H(jj,j,ay)*rowsum(1:2)
enddo
childCoords(1:2,i,j,c) = acc(1:2)
enddo
enddo
enddo
endsubroutine SubdivideNodeCoords