PetIGA / src / petigageo.f90.in

! -*- f90 -*-

pure subroutine GeometryMap(&
     order,&
     nen,X,&
     N0,N1,N2,N3,&
     R0,R1,R2,R3,&
     DetG,Grad,InvG)
  use PetIGA
  implicit none
  !integer(kind=IGA_INT), parameter   :: dim = 1,2,3
  integer(kind=IGA_INT ), intent(in)  :: order
  integer(kind=IGA_INT ), intent(in)  :: nen
  real   (kind=IGA_REAL), intent(in)  :: X(dim,nen)
  real   (kind=IGA_REAL), intent(in)  :: N0(            nen)
  real   (kind=IGA_REAL), intent(in)  :: N1(        dim,nen)
  real   (kind=IGA_REAL), intent(in)  :: N2(    dim,dim,nen)
  real   (kind=IGA_REAL), intent(in)  :: N3(dim,dim,dim,nen)
  real   (kind=IGA_REAL), intent(out) :: R0(            nen)
  real   (kind=IGA_REAL), intent(out) :: R1(        dim,nen)
  real   (kind=IGA_REAL), intent(out) :: R2(    dim,dim,nen)
  real   (kind=IGA_REAL), intent(out) :: R3(dim,dim,dim,nen)
  real   (kind=IGA_REAL), intent(out) :: DetG
  real   (kind=IGA_REAL), intent(out) :: Grad(dim,dim)
  real   (kind=IGA_REAL), intent(out) :: InvG(dim,dim)

  integer(kind=IGA_INT ) :: idx
  integer(kind=IGA_INT ) :: i, j, k, l
  integer(kind=IGA_INT ) :: a, b, c, d
  real   (kind=IGA_REAL) :: X1(dim,dim)
  real   (kind=IGA_REAL) :: X2(dim,dim,dim)
  real   (kind=IGA_REAL) :: X3(dim,dim,dim,dim)
  real   (kind=IGA_REAL) :: E1(dim,dim)
  real   (kind=IGA_REAL) :: E2(dim,dim,dim)
  real   (kind=IGA_REAL) :: E3(dim,dim,dim,dim)

  ! gradient of the mapping
  Grad = matmul(N1,transpose(X))
  DetG = Determinant(dim,Grad)
  InvG = Inverse(dim,DetG,Grad)

  ! 0th derivatives
  R0 = N0

  ! 1st derivatives
  X1 = transpose(Grad)
  E1 = transpose(InvG)
  R1 = 0
  do idx = 1,nen
     do i = 1,dim
        do a = 1,dim
           R1(i,idx) = N1(a,idx)*E1(a,i) +  R1(i,idx)
        end do
     end do
  end do

  ! 2nd derivatives
  if (order < 2) return
  X2 = 0
  do b = 1,dim
     do a = 1,dim
        do i = 1,dim
           do idx = 1,nen
              X2(i,a,b) = X(i,idx)*N2(a,b,idx) + X2(i,a,b)
           end do
        end do
     end do
  end do
  E2 = 0
  do j = 1,dim
     do i = 1,dim
        do c = 1,dim
           do b = 1,dim
              do a = 1,dim
                 do k = 1,dim
                    E2(c,i,j) = - X2(k,a,b)*E1(a,i)*E1(b,j)*E1(c,k) + E2(c,i,j)
                 end do
              end do
           end do
        end do
     end do
  end do
  R2 = 0
  do idx = 1,nen
     do j = 1,dim
        do i = 1,dim
           do b = 1,dim
              do a = 1,dim
                 R2(i,j,idx) = N2(a,b,idx)*E1(a,i)*E1(b,j) + R2(i,j,idx)
              end do
              R2(i,j,idx) = N1(b,idx)*E2(b,i,j) + R2(i,j,idx)
           end do
        end do
     end do
  end do

  ! 3rd derivatives
  if (order < 3) return
  X3 = 0
  do c = 1,dim
     do b = 1,dim
        do a = 1,dim
           do i = 1,dim
              do idx = 1,nen
                 X3(i,a,b,c) = X(i,idx)*N3(a,b,c,idx) + X3(i,a,b,c)
              end do
           end do
        end do
     end do
  end do
  E3 = 0
  do k = 1,dim
     do j = 1,dim
        do i = 1,dim
           do d = 1,dim
              do a = 1,dim
                 do b = 1,dim
                    do l = 1,dim
                       do c = 1,dim
                          E3(d,i,j,k) = - X3(l,c,b,a)*E1(c,i)*E1(b,j)*E1(a,k)*E1(d,l) + E3(d,i,j,k)
                       end do
                       E3(d,i,j,k) = - X2(l,b,a)*( E1(a,j)*E2(b,i,k) &
                            + E1(a,k)*E2(b,i,j) &
                            + E1(b,i)*E2(a,j,k) )*E1(d,l) + E3(d,i,j,k)
                    end do
                 end do
              end do
           end do
        end do
     end do
  end do
  R3 = 0
  do idx = 1,nen
     do k = 1,dim
        do j = 1,dim
           do i = 1,dim
              do a = 1,dim
                 do b = 1,dim
                    do c = 1,dim
                       R3(i,j,k,idx) = N3(c,b,a,idx)*E1(c,i)*E1(b,j)*E1(a,k) + R3(i,j,k,idx)
                    end do
                    R3(i,j,k,idx) = N2(b,a,idx)*( E1(a,j)*E2(b,i,k) &
                         + E1(a,k)*E2(b,i,j) &
                         + E1(b,i)*E2(a,j,k) ) + R3(i,j,k,idx)
                 end do
                 R3(i,j,k,idx) = N1(a,idx)*E3(a,i,j,k) + R3(i,j,k,idx)
              end do
           end do
        end do
     end do
  end do

end subroutine GeometryMap

include 'petigainv.f90.in'
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.