1!$Id:$ 2 subroutine interp1(nr,xs,side,is,ns,ndm,shp, x) 3 4! * * F E A P * * A Finite Element Analysis Program 5 6!.... Copyright (c) 1984-2017: Regents of the University of California 7! All rights reserved 8 9!-----[--.----+----.----+----.-----------------------------------------] 10! Purpose: Construct one dimensional Lagrange interpolation for coords. 11 12! Inputs: 13! nr - Number of increments on side 14! xs(3,*) - Nodal values of interpolation function 15! side - side number (check sign) 16! is(*) - List of side nodes 17! ns - Order of Lagrange polynomial for side 18! ndm - Spatial dimension of mesh 19! shp(*) - Shape functions for nodal values 20 21! Outputs: 22! x(ndm,ip) - Coordinates of points 23!-----[--.----+----.----+----.-----------------------------------------] 24 25 implicit none 26 27 integer i,i1, j, m,m1,m2, n1,n2,n3,nr,ns,ndm, side, is(*) 28 real*8 den, rn, xi, xid,xii,xij, xs(3,*), shp(*), x(ndm,*) 29 30 save 31 32! Set order of using nodes on side 33 34 if(side.gt.0) then 35 m1 = 1 36 m2 = 2 37 n1 = 3 38 n2 = ns 39 n3 = 1 40 elseif(side.lt.0) then 41 m1 = 2 42 m2 = 1 43 n1 = ns 44 n2 = 3 45 n3 = -1 46 endif 47 48! Left and right ends 49 50 do j = 1,ndm 51 x(j, 1) = xs(j,is(m1)) 52 x(j,nr+1) = xs(j,is(m2)) 53 end do ! j 54 55 xid = 1.0d0/dble(nr) 56 rn = 1.0d0/dble(ns-1) 57 58! Loop through interior points 59 60 do m = 2,nr 61 62 xi = dble(m-1)*xid 63 64! Compute Lagrange interpolation function at point 'xi' 65 66! End nodes 67 68 shp(1) = 1.0d0 69 den = 1.0d0 70 xii = 0.0d0 71 do j = 2,ns 72 xij = dble(j-1)*rn 73 shp(1) = shp(1)*(xi - xij) 74 den = den * (xii - xij) 75 end do ! j 76 shp(1) = shp(1)/den 77 78 shp(2) = 1.0d0 79 den = 1.0d0 80 xii = 1.0d0 81 do j = 1,ns-1 82 xij = dble(j-1)*rn 83 shp(2) = shp(2)*(xi - xij) 84 den = den * (xii - xij) 85 end do ! j 86 shp(2) = shp(2)/den 87 88 do i = 3,ns 89 90 shp(i) = 1.0d0 91 den = 1.0d0 92 xii = dble(i-2)*rn 93 do j = 1,ns 94 if(i-1.ne.j) then 95 xij = dble(j-1)*rn 96 shp(i) = shp(i)*(xi - xij) 97 den = den * (xii - xij) 98 end if 99 end do ! j 100 shp(i) = shp(i)/den 101 102 end do ! i 103 104! Perform interpolation for current shape function 105 106 do j = 1,ndm 107 x(j,m) = shp(1)*xs(j,is(m1)) + shp(2)*xs(j,is(m2)) 108 end do ! j 109 110 i1 = 2 111 do i = n1,n2,n3 112 i1 = i1 + 1 113 do j = 1,ndm 114 x(j,m) = x(j,m) + shp(i1)*xs(j,is(i)) 115 end do ! j 116 end do ! i 117 end do ! m 118 119 end 120