Sources Esope | Cast3M

ravc
C RAVC      SOURCE    CB215821  20/11/25    13:38:29     10792          
      subroutine ravc
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMELEME
-INC SMCOORD
      segment icpr(nbpts)
      call lirobj ( 'CHPOINT',MCHPO1, 1,iretou)
      if(ierr . ne.0) return
      call lirobj ( 'CHPOINT',MCHPO2, 1,iretou)
      if(ierr . ne.0) return
      call lirobj ( 'CHPOINT',MCHPO3, 1,iretou)
      if(ierr . ne.0) return
      segact mchpo1, mchpo2,mchpo3
      if(mchpo1.ipchp(/1). ne . 1) then
*       write(6,*) ' erreur 1 '
       call erreur(19)
       return
      endif
      if(mchpo2.ipchp(/1). ne . 1) then
*       write(6,*) ' erreur 2 '
       call erreur(19)
       return
      endif
      if(mchpo3.ipchp(/1). ne . 1) then
*       write(6,*) ' erreur 3 '
       call erreur(19)
       return
      endif
      msoup1= mchpo1.ipchp(1)
      msoup2= mchpo2.ipchp(1)
      msoup3= mchpo3.ipchp(1)
      segact msoup1,msoup2,msoup3
      ipt1=msoup1.igeoc
      ipt2=msoup2.igeoc
      ipt3=msoup3.igeoc
      segact ipt1,ipt2,ipt3
      segini icpr
      n1 = ipt1.num(/2)
      n2 = ipt2.num(/2)
      n3 = ipt3.num(/2)
      ipo=0
      do 100 i=1,n1
       ia=ipt1.num(1,i)
       if(icpr(ia).eq.0) then
         ipo=ipo+1
         icpr(ia)=ipo
       endif
  100 continue
      do 101 i=1,n2
       ia=ipt2.num(1,i)
       if(icpr(ia).eq.0) then
         ipo=ipo+1
         icpr(ia)=ipo
       endif
  101 continue
      do 102 i=1,n3
       ia=ipt3.num(1,i)
       if(icpr(ia).eq.0) then
         ipo=ipo+1
         icpr(ia)=ipo
       endif
  102 continue
      nbelem=ipo
      nbnn=1
      nbsous=0
      nbref=0
      segini meleme
      itypel=1
      do 103 i=1,icpr(/1)
       if(icpr(I).ne.0) then
        num(1,icpr(i))=i
       endif
  103 continue
 
*      write(6,*) ' ipt1',(ipt1.num(1,iou),iou=1,ipt1.num(/2))
*      write(6,*) ' ipt2',(ipt2.num(1,iou),iou=1,ipt2.num(/2))
*      write(6,*) ' ipt3',(ipt3.num(1,iou),iou=1,ipt3.num(/2))
*      if( (n1.ne.n2). or .(n1.ne.n3))then
*          write(6,*) ' erreur 4 '
*          call erreur(19)
*          return
*      endif
      if( msoup1.noharm(/1) . ne . 1) then
*          write(6,*) ' erreur 5 '
       call erreur(19)
       return
      endif
      if( msoup2.noharm(/1) . ne . 1) then
*          write(6,*) ' erreur 6 '
       call erreur(19)
       return
      endif
      if( msoup3.noharm(/1) . ne . 1) then
*          write(6,*) ' erreur 7 '
       call erreur(19)
       return
      endif
      segini,mchpoi=mchpo1
      segini,msoupo=msoup1
      ipchp(1)=msoupo
      mpova1=msoup1.ipoval
      mpova2=msoup2.ipoval
      mpova3=msoup3.ipoval
      segact mpova1,mpova2,mpova3
      n=nbelem
      nc=1
      segini mpoval
      ipoval=mpoval
      igeoc=meleme
*      write(6,*) ' meleme' , ( num(1,iou),iou=1,nbelem)
*
* on cherche les max de chaque champs
*
*      amma=0.
*      bmma=0.
*      cmma=0.
*      do 10 i=1,nbelem
*      if(amma.le.(abs(mpova1.vpocha(i,1))))amma=abs(mpova1.vpocha(i,1))
*      if(bmma.le.(abs(mpova2.vpocha(i,1))))bmma=abs(mpova2.vpocha(i,1))
*      if(cmma.le.(abs(mpova3.vpocha(i,1))))cmma=abs(mpova3.vpocha(i,1))
*   10 continue
*
*    calcul des racines
*
      do 1 j=1,nbelem
      ia = num(1,j)
*      write(6,*) ' ia ' , ia
       a = 0.d0
       do 8 i=1,n1
         if( ipt1.num(1,i).eq.ia)then
         a = mpova1.vpocha(i,1)
          go to 9
         endif
   8   continue
   9   continue
       b=0.d0
       do 2 k=1,n2
        if(ipt2.num(1,k).eq.ia) then
          b = mpova2.vpocha(k,1)
          go to 3
        endif
   2  continue
   3  continue
      c=0.d0
      do 4 l=1,n3
        if(ipt3.num(1,l).eq.ia)then
        c = mpova3.vpocha(l,1)
        go to 5
        endif
   4  continue
   5  continue
 
* test sur les grandeurs
*     write (6,*) ' a b c ',a,b,c
      if( a . eq. 0.d0. or . abs(a).lt. (abs(b)*1.e-13))then
* on resoud  une seule equation
        if( b . ne . 0.d0)  then
          xx = -c / b
          if( xx .ge. 0.d0 ) xx = 0.d0
        else
          xx = 0.d0
        endif
      else
        xx=0.d0
        delta2 =  ( b*b - 4*a*c)
       if (delta2.lt. b*b*1d-13) then
         xx = -b/(2.d0*a)
         xx = min(xx,0.d0)
         xx = 0.d0
**       write (6,*) ' ravc discriminant negatif ',delta2,b
       else
        delta2=max(0.d0,delta2)
        delta = sqrt (delta2)
        s1= (-b + delta )/ (2.*a)
        s2= (-b - delta) / (2.*a)
          s1=min(s1,0.d0)
          s2=min(s2,0.d0)
*        write (6,*) ' s1 s2 ' , s1,s2
        if( s1. gt . s2) then
           if( s1 . lt. 0.d0) then
             xx = s1
           else
           if(s2.lt.0.d0) xx = s2
           endif
        else
           if( s2 . lt. 0.d0) then
             xx = s2
           else
           if(s1.lt.0.d0) xx = s1
           endif
        endif
       endif
      endif
**    write (6,*) ' xx ',xx
      vpocha(j,1)=xx
   1  continue
      segdes mpova1,mpova2,mpova3,mpoval
      segdes ipt1,ipt2,ipt3,meleme
      segdes msoupo,msoup1,msoup2,msoup3
      segdes mchpoi,mchpo1,mchpo2,mchpo3
      segsup icpr
      call ecrobj('CHPOINT',mchpoi)
      return
      end