Télécharger volshb.eso

Retour à la liste

volshb
C VOLSHB    SOURCE    CB215821  20/11/25    13:42:36     10792          
      subroutine volshb(ipt1,epai,mchpoi,ip1,meleme)
      implicit real*8(a-h,o-z)
      implicit integer(i-n)
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMELEME
-INC SMCOORD
      segment icpr(nbpts)
      segment ipos(nbp)
      segment vec(3,nbp)
      segment icprp(nbpts)
      segact ipt1
      if( ipt1.lisous(/1).ne.0.or.ipt1.itypel.ne.8) then
       call erreur (1002)
       return
      endif
      icle=1
      if(mchpoi.ne.0) then
       icle=2
       segact mchpoi
       if(ipchp(/1).ne.1) then
          call erreur(180)
          return
       endif
       msoupo=ipchp(1)
       segact msoupo
       if(nocomp(/2).ne.1) then
          call erreur(180)
          return
       endif
       mpoval=ipoval
       ipt2=igeoc
       segact ipt2,mpoval
       segini icprp
       do i=1,ipt2.num(/2)
          ia=ipt2.num(1,i)
          icprp(ia)=i
       enddo
       segdes ipt2
      endif
*
* pour chaque point on regarde quels sont les éléments qui le touchent
*
       xp1=xcoor((ip1-1)*(idim+1)+1)
       yp1=xcoor((ip1-1)*(idim+1)+2)
       zp1=xcoor((ip1-1)*(idim+1)+3)
       nbp=0
       segini icpr
       do i=1,ipt1.num(/2)
          do j=1,4
           ia= ipt1.num(j,i)
           if(icpr(ia).eq.0) then
             nbp=nbp+1
             icpr(ia)=nbp
           endif
          enddo
       enddo
       segini ipos,vec
*
* on calcule pour chaque point le vecteur normale
*
       isens=0
       do i=1,ipt1.num(/2)
           ia=(ipt1.num(1,i)-1)*4
           ib=(ipt1.num(2,i)-1)*4
           ic=(ipt1.num(3,i)-1)*4
           xa=xcoor(ia+1)
           ya=xcoor(ia+2)
           za=xcoor(ia+3)
           xb=xcoor(ib+1)
           yb=xcoor(ib+2)
           zb=xcoor(ib+3)
           xc=xcoor(ic+1)
           yc=xcoor(ic+2)
           zc=xcoor(ic+3)
           xn=(yb-ya)*(zc-za)-(yc-ya)*(zb-za)
           yn=(zb-za)*(xc-xa)-(zc-za)*(xb-xa)
           zn=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya)
           xyz=sqrt(xn*xn+yn*yn+zn*zn)
           xn=xn/xyz
           yn=yn/xyz
           zn=zn/xyz
           do j=1,4
             ia=ipt1.num(j,i)
             ib=icpr(ia)
             ipos(ib)=ipos(ib)+1
             vec(1,ib)=vec(1,ib)+xn
             vec(2,ib)=vec(2,ib)+yn
             vec(3,ib)=vec(3,ib)+zn
           enddo
        enddo
        do i=1,nbp
          na= ipos(i)
          ipos(i)=0
          vec(1,i)=vec(1,i)/na
          vec(2,i)=vec(2,i)/na
          vec(3,i)=vec(3,i)/na
          xv= vec(1,i)*vec(1,i)+vec(2,i)*vec(2,i)+vec(3,i)*vec(3,i)
          xv=sqrt(xv)
          vec(1,i)=vec(1,i)/xv
          vec(2,i)=vec(2,i)/xv
          vec(3,i)=vec(3,i)/xv
          xps= vec(1,i)*xp1+vec(2,i)*yp1+vec(3,i)*zp1
          if(xps.ge.0.d0) then
            if( isens.lt.0) then
              call erreur (19)
              return
            endif
            isens=1
          elseif(xps.le.0.d0) then
            if( isens.gt.0) then
              call erreur (19)
              return
            endif
            isens=-1
          endif
        enddo
*
* il faut maintenant créer le maillage
*
        nbnn=8
        nbelem=ipt1.num(/2)
        nbref=2
        nbsous=0
        segini meleme
        itypel=14
        ep=epai
        segact mcoord*mod
        nbpini=nbpts
        nbpts=nbpts+2*nbp
        segadj mcoord
        nbnn=4
        nbref=0
        segini ipt5,ipt6
        lisref(1)=ipt5
        lisref(2)=ipt6
        ipt5.itypel=8
        ipt6.itypel=8
        do i=1,ipt1.num(/2)
          do j=1,4
           ia=ipt1.num(j,i)
           ib=icpr(ia)
* faut-il créer les points ?
           if( ipos(ib).eq.0) then
             iadec=(ia-1)*(idim+1)
             if(icle.eq.2) then
               ic=icprp(ia)
               if(ic.eq.0) then
                 meleme=0
                 call erreur(19)
                 return
               endif
               ep=vpocha(ic,1)
             endif
             nbpini=nbpini+1
             nbpdec=(nbpini-1)*(idim+1)
             ie=nbpini
             eps=ep*isens
             do k=1,3
               xcoor(nbpdec+k)= xcoor(iadec+k)-vec(k,ib)*eps
             enddo
             nbpini=nbpini+1
             nbpdec=nbpdec+4
             do k=1,3
               xcoor(nbpdec+k)= xcoor(iadec+k)+vec(k,ib)*eps
             enddo
             ipos(ib)=ie
           else
             ie=ipos(ib)
           endif
           num(J,i)=ie
           num(J+4,i)=ie+1
           ipt5.num(j,i)=ie
           ipt6.num(j,i)=ie+1
          enddo
        enddo
        if( icle.eq.2)then
          segdes mpoval,ipt2,msoupo,mchpoi
          segsup icprp
        endif
        segsup icpr,ipos,vec
        segdes ipt1
        segdes meleme,ipt5,ipt6
        return
        end
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales