1 subroutine f3lhs (shpb, shglb, xlb, flhsl, 2 & fnrml, sgn ) 3c 4c---------------------------------------------------------------------- 5c 6c This subroutine computes the element LHS matrix and the normal 7c to the boundary for computation of (output) boundary fluxes. 8c 9c input: 10c shpb (nen,nintg) : boundary element shape-functions 11c shglb (nsd,nen,nintg) : boundary element grad-shape-functions 12c wghtb (nintg) : boundary element weight 13c xlb (npro,nenl,nsd) : nodal coordinates 14c sgn (npro,nshl) : mode signs for hierarchic basis 15c 16c output: 17c flhsl (npro,nenl,1) : element lumped lhs on flux boundary 18c fnrml (npro,nenl,nsd) : RHS of LS projection of normal to 19c flux boundary 20c 21c 22c Note: Special lumping technique is used to compute the LHS. 23c See T.J.R. Hughes, "The Finite Element Method: Linear 24c Static and Dynamic Finite Element Analysis", page 445. 25c 26c Note: Least-squares projection is used to compute the normal to 27c the boundary at the nodes. This routine provides the element 28c contribution to the RHS of the projection linear system. 29c 30c 31c Zdenek Johan, Summer 1991. 32c---------------------------------------------------------------------- 33c 34 include "common.h" 35c 36 dimension shpb(nshl,ngaussb), shglb(nsd,nshl,ngaussb), 37 & xlb(npro,nenl,nsd), 38 & flhsl(npro,nshl,1), fnrml(npro,nshl,nsd) 39c 40 dimension WdetJb(npro), 41 & bnorm(npro,nsd), fmstot(npro), 42 & temp(npro), temp1(npro), 43 & temp2(npro), temp3(npro) 44 45 dimension sgn(npro,nshl), shape(npro,nshl), 46 & shdrv(npro,nsd,nshl), 47 & v1(npro,nsd), v2(npro,nsd) 48c 49c.... integrate the lumped LHS matrix and normal 50c 51 fmstot = zero 52c 53c 54c.... compute the normal to the boundary 55c 56 57 v1 = xlb(:,2,:) - xlb(:,1,:) 58 v2 = xlb(:,3,:) - xlb(:,1,:) 59 60 if (lcsyst .eq. 1) then 61 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 62 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 63 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 64 else 65 temp1 = - v1(:,2) * v2(:,3) + v2(:,2) * v1(:,3) 66 temp2 = - v2(:,1) * v1(:,3) + v1(:,1) * v2(:,3) 67 temp3 = - v1(:,1) * v2(:,2) + v2(:,1) * v1(:,2) 68 endif 69c 70 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 71 bnorm(:,1) = temp1 * temp 72 bnorm(:,2) = temp2 * temp 73 bnorm(:,3) = temp3 * temp 74 75 do intp = 1, ngaussb 76c 77c.... get the hierarchic shape functions at this int point 78c 79 call getshp(shpb, shglb, sgn, 80 & shape, shdrv) 81c 82 83c$$$ WdetJb = Qwtb(lcsyst,intp) / (four*temp) 84 if (lcsyst .eq. 1) then 85 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 86 elseif (lcsyst .eq. 2) then 87 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 88 elseif (lcsyst .eq. 3) then 89 WdetJb = Qwtb(lcsyst,intp) / (two*temp) 90 elseif (lcsyst .eq. 4) then 91 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 92 elseif (lcsyst .eq. 5) then 93 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 94 elseif (lcsyst .eq. 6) then 95 WdetJb = Qwtb(lcsyst,intp) / (two*temp) 96 endif 97 98c 99c.... compute the lumped LHS and normal 100c 101 do n = 1, nenl ! when changed to nshl solution degraded ipord 10 102 flhsl(:,n,1) = flhsl(:,n,1) + WdetJb * shape(:,n) 103 104c for curved geometries the below construct for the normals has to be used 105 fnrml(:,n,1) = fnrml(:,n,1) + WdetJb * bnorm(:,1) 106 & * shape(:,n) 107 fnrml(:,n,2) = fnrml(:,n,2) + WdetJb * bnorm(:,2) 108 & * shape(:,n) 109 fnrml(:,n,3) = fnrml(:,n,3) + WdetJb * bnorm(:,3) 110 & * shape(:,n) 111 enddo 112c 113c To best represent this case it should be assigned to the vertex 114c modes and higher entities should get zero as is done below 115c 116 fmstot = fmstot + WdetJb 117c 118 enddo 119 120c$$$ do i=1,nenl 121c$$$ fnrml(:,i,:)=bnorm(:,:) 122c$$$ enddo 123 if(ipord.gt.1) fnrml(:,nenl:nshl,:)=zero 124c 125c.... scale the LHS matrix contribution 126c 127 temp = zero 128 do n = 1, nshl 129 temp = temp + flhsl(:,n,1) 130 enddo 131c 132 do n = 1, nshl 133 flhsl(:,n,1) = flhsl(:,n,1) * fmstot / temp 134 enddo 135c 136c.... return 137c 138 return 139 end 140