Source

PetIGA / src / petigageo.f90.in

Full commit
Lisandro Dalcin efb6552 

Lisandro Dalcin 8192bd7 
Lisandro Dalcin efb6552 



Lisandro Dalcin 8192bd7 
Lisandro Dalcin efb6552 

Lisandro Dalcin ca3d144 














Lisandro Dalcin efb6552 
Lisandro Dalcin ca3d144 








Lisandro Dalcin efb6552 

Lisandro Dalcin 8192bd7 
Lisandro Dalcin efb6552 
Lisandro Dalcin 8192bd7 
Lisandro Dalcin efb6552 




Lisandro Dalcin 8192bd7 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 









Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 








Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 












Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 














Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 










Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 



















Lisandro Dalcin 791c6f2 
Lisandro Dalcin efb6552 



















Lisandro Dalcin 8192bd7 
Lisandro Dalcin efb6552 
Lisandro Dalcin 5e19e94 
! -*- 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_INTEGER_KIND),parameter   :: dim = 1,2,3
  integer(kind=IGA_INTEGER_KIND), intent(in)  :: order
  integer(kind=IGA_INTEGER_KIND), intent(in)  :: nen
  real   (kind=IGA_REAL_KIND   ), intent(in)  :: X(dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N0(            nen)
  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N1(        dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N2(    dim,dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N3(dim,dim,dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: R0(            nen)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: R1(        dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: R2(    dim,dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: R3(dim,dim,dim,nen)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: DetG
  real   (kind=IGA_REAL_KIND   ), intent(out) :: Grad(dim,dim)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: InvG(dim,dim)

  integer(kind=IGA_INTEGER_KIND)  :: idx
  integer(kind=IGA_INTEGER_KIND)  :: i, j, k, l
  integer(kind=IGA_INTEGER_KIND)  :: a, b, c, d
  real   (kind=IGA_REAL_KIND   )  :: X1(dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: X2(dim,dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: X3(dim,dim,dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: E1(dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: E2(dim,dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: 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'