sortho31.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    sortho31   ../engine/source/elements/solid/solide/sortho31.F
!||--- called by ------------------------------------------------------
!||    srota6     ../engine/source/output/anim/generate/srota6.F
!||    srotorth   ../engine/source/elements/solid/srotorth.F
!||    thsol      ../engine/source/output/th/thsol.F
!||====================================================================
      SUBROUTINE sortho31(
     .   X1, X2, X3, X4, X5, X6, X7, X8,
     .   Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,
     .   Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8,
     .   E1X, E2X, E3X, E1Y, E2Y, E3Y, E1Z, E2Z, E3Z )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      my_real
     .   x1, x2, x3, x4, x5, x6, x7, x8,
     .   y1, y2, y3, y4, y5, y6, y7, y8,  
     .   z1, z2, z3, z4, z5, z6, z7, z8,  
     .   e1x, e1y, e1z,
     .   e2x, e2y, e2z,
     .   e3x, e3y, e3z
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER N,NITER
      my_real
     .   X17 , X28 , X35 , X46,
     .   Y17 , Y28 , Y35 , Y46,
     .   z17 , z28 , z35 , z46,
     .   a17 , a28 ,
     .   b17 , b28 ,
     .   c17 , c28 ,
     .   rx , ry , rz ,
     .   sx , sy , sz ,
     .   tx , ty , tz 
      my_real
     .   aa,bb
      DATA niter/3/
C-----------------------------------------------
      x17=x7-x1
      x28=x8-x2
      x35=x5-x3
      x46=x6-x4
      y17=y7-y1
      y28=y8-y2
      y35=y5-y3
      y46=y6-y4
      z17=z7-z1
      z28=z8-z2
      z35=z5-z3
      z46=z6-z4
      rx=x17+x28-x35-x46
      ry=y17+y28-y35-y46
      rz=z17+z28-z35-z46
      a17=x17+x46
      a28=x28+x35
      b17=y17+y46
      b28=y28+y35
      c17=z17+z46
      c28=z28+z35
      sx=a17+a28
      sy=b17+b28
      sz=c17+c28
      tx=a17-a28
      ty=b17-b28
      tz=c17-c28
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c     norme r s t
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      aa = sqrt(rx*rx + ry*ry + rz*rz)
      if ( aa/=zero) aa = one / aa
      rx = rx * aa
      ry = ry * aa
      rz = rz * aa
      aa = sqrt(sx*sx + sy*sy + sz*sz)
      if ( aa/=zero) aa = one / aa
      sx = sx * aa
      sy = sy * aa
      sz = sz * aa
      aa = sqrt(tx*tx + ty*ty + tz*tz)
      if ( aa/=zero) aa = one / aa
      tx = tx * aa
      ty = ty * aa
      tz = tz * aa
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c     iterations
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      n=0
111   CONTINUE
      n=n+1
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       e1x = sy * tz - sz * ty + rx
       e1y = sz * tx - sx * tz + ry
       e1z = sx * ty - sy * tx + rz
c
       e2x = ty * rz - tz * ry + sx
       e2y = tz * rx - tx * rz + sy
       e2z = tx * ry - ty * rx + sz
c
       e3x = ry * sz - rz * sy + tx
       e3y = rz * sx - rx * sz + ty
       e3z = rx * sy - ry * sx + tz
c
       bb = sqrt(e1x*e1x + e1y*e1y + e1z*e1z)
       if ( bb/=zero) bb = one / bb
       rx = e1x * bb
       ry = e1y * bb
       rz = e1z * bb
c
       bb = sqrt(e2x*e2x + e2y*e2y + e2z*e2z)
       if ( bb/=zero) bb = one / bb
       sx = e2x * bb
       sy = e2y * bb
       sz = e2z * bb
c
       bb = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
       if ( bb/=zero) bb = one / bb
       tx = e3x * bb
       ty = e3y * bb
       tz = e3z * bb
c
      IF (n<niter) GOTO 111
c     norme et orthogonalisation
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      e1x = rx
      e1y = ry
      e1z = rz
c
      e3x = e1y * sz - e1z * sy
      e3y = e1z * sx - e1x * sz
      e3z = e1x * sy - e1y * sx
c
      aa = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
      if ( aa/=zero) aa = one / aa
      e3x = e3x * aa
      e3y = e3y * aa
      e3z = e3z * aa
c
      e2x = e3y * e1z - e3z * e1y
      e2y = e3z * e1x - e3x * e1z
      e2z = e3x * e1y - e3y * e1x
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      RETURN
      END