Patran Neutral Mesh File Writer updated

Sphere

A recent comment from a blog reader made me relook at the patran neutral mesh file writer routine that I had on this blog.

Since the code was pasted on this blog, it lacked proper update and ability to have a test case or example file.

So spent this weekend, creating an example file and have uploaded the code and the example on github.

The entire writer is now availabe at https://github.com/sukhbinder/PatranMeshWriter

Hope to add and update the code more often. 🙂

 

Advertisements

Few Python Regex for FORTRAN Code

 

Most of my day job involves Fortran code and finding quick details is a regular need. Here are four regex in python that I have been using for some time now. Sharing it here.

To Find Subroutines and Functions

 import re procedule_name=re.compile(r'\s*(RECURSIVE)?\s+(SUBROUTINE|FUNCTION)\s+\S+\(.*',re.IGNORECASE) 

 

To Find All Call Statements
			

 callname=re.compile(r'^\s*call\s*\S*[(]?\S*[)]?',re.IGNORECASE) 

    

 

To Find Variables (Integer, Double Precision, Logical, Real)
			

 variregex=re.compile(r'^\s+(logical|integer|double.+precision|real(\*\s*\d)?)\s+\S*',re.IGNORECASE) 

 

To Find Character Variables
			

 characteregex=re.compile(r'^\s+(character+([\*\s*\d+]|\(\S*)*)\s+\S*',re.IGNORECASE) 

 

Simple Python code to Extract Fortran Routines and its Argument

fortranRoutinesWithPython
Suppose you got a long fortran file full of code, and want to parse it and get the subroutine/functions defined in the file along with its argument.

Here’s a simple python code to extract all this info in one go

import re
from collections import defaultdict

def processsubroutine(line):
    if ')' in line:
        return line.split("(")[1].split(')')[0].split(',')
    

def collectsubroutines(fname):
    procnam=defaultdict()
    fort_proc_def = re.compile(r'\s*(RECURSIVE)?\s*(SUBROUTINE|FUNCTION)\s+\S+\(*', re.IGNORECASE)
    fp = open(fname,"r")
    datalines=fp.readlines()
    fp.close()
    for i,line in enumerate(datalines):
        if fort_proc_def.match(line):
            procnam[(line.split()[1].split('(')[0])]=processsubroutine(line)
    return procnam

def main():
    for key,item  in procnam.items():
        print key,item

Lights Out Game In Fortran

Here’s an implementation of the Lights out game in fortran. Wrote this in 2011 and found it yesterday. Works as it is and can be improved. Checkout the screenshot.

LightsOut_game_in_fortran

For those who don’t know, here’s a brief description from wikipedia

The game consists of a 5 by 5 grid of lights. When the game starts, a random number or a stored pattern of these lights is switched on. Pressing any of the lights will toggle it and the four adjacent lights. The goal of the puzzle is to switch all the lights off, preferably in as few button presses as possible.



module lightsout
implicit none
integer :: elements(0:25)
integer :: iswon,moves
logical,parameter :: debug =.false.
 character*30 :: sdata(4)

 contains

subroutine printtable1()
 integer :: i,j,ii,jj

 sdata(1)='        _\|/_        '
 sdata(2)='        (o o)        '
 sdata(3)='+----oOO-{_}-OOo----+'
 
 do i=1,3
  write(*,'(a30)') sdata(i)
 end do
    write(*,'(6(2x,a))')  ' ','a','b','c','d','e'
    ii=0
    jj=4
    do j=1,5
!      write(*,'(6i3)') j,(elements(i),i=ii,jj)
      write(*,'(i3,5a3)') j,(achar(elements(i)+46),i=ii,jj)
      ii=jj+1
      jj=ii+4
    end do

!    if(debug) then
!    print *,"Degub mode:"
!    write(*,'(5i8)')(elements(i),i=0,24)
!    end if
 
end subroutine



subroutine printtable()
 integer :: i,j,ii,jj

 sdata(1)='        _\|/_        '
 sdata(2)='        (o o)        '
 sdata(3)='+----oOO-{_}-OOo----+'
 
 do i=1,3
  write(*,'(a30)') sdata(i)
 end do
    write(*,'(6(1x,a,2x))')  ' ','a','b','c','d','e'
    ii=1
    
    write(*,'(7(5a))')'    ','___','.','___','.','___','.','___','.','___','.'
    do i=1,5
       do j=1,5
        
        if(elements(ii) .eq. 0) then
          write(*,'(4a)',advance='no') achar(32),achar(32),achar(32),'|'
        else
          write(*,'(4a)',advance='no')achar(219),achar(219),achar(219),'|'
        end if
         ii=ii+1
       end do
       write(*,*)
       write(*,'(i3,5(5a))')i,'.','___','|','___','|','___','|','___','|','___','|'
    end do   

!    if(debug) then
!    print *,"Degub mode:"
!    write(*,'(5i8)')(elements(i),i=0,24)
!    end if
 
end subroutine

subroutine fillall()
integer :: i
 do i =0 , 24
   elements(i)=1
 end do
end subroutine

subroutine clearall()
integer :: i
 do i = 0, 24
   elements(i)=0
 end do
end subroutine

subroutine level1()
 iswon=0
 moves=0
 call clearall()
 elements(10)=1
 elements(12)=1
 elements(14)=1
end subroutine

subroutine level2()
 iswon=0
 moves=0
 call clearall();
 elements (12)=1;
 elements (16)=1;
 elements (17)=1;
 elements (18)=1;
 elements (20)=1;
 elements (21)=1;
 elements (22)=1;
 elements (23)=1;
 elements (24)=1;
end subroutine

subroutine level3()
 iswon=0
 moves=0
 call fillall();
 elements (4)=0;
 elements (6)=0;
 elements (7)=0;
 elements (8)=0;
 elements (11)=0;
 elements (12)=0;
 elements (13)=0;
 elements (16)=0;
 elements (17)=0;
 elements (18)=0;
 elements (24)=0;
end subroutine

subroutine level4()
 iswon=0
 moves=0
 call clearall();
 elements (2)=1;
 elements (6)=1;
 elements (8)=1;
 elements (10)=1;
 elements (12)=1;
 elements (14)=1;
 elements (16)=1;
 elements (18)=1;
 elements (22)=1;
end subroutine

subroutine level5()
 iswon=0
 moves=0
 call clearall();
 elements (11)=1;
 elements (16)=1;
 elements (21)=1;
end subroutine

subroutine checkall() 
 integer :: win,i
 win=1;
 do i = 0 ,24
   if (elements (i)==0) then
      win=0;
   end if
 end do

 if (win==1) then
    call printtable()
    write(*,*) "You Won!!"
    write(*,*) "You did it in ",moves, "moves"
    iswon=1
 end if
if(debug) print *,'7',iswon
end subroutine

subroutine check(v) 
integer :: q,w,row,intv,a,b,c,d,kc,kd,kv,is,v
 q=5; 
 w=1;
 moves=moves+1
 row=Int(v/q)+w;
 intv =Int(v);
 a=intv+q;
 b=intv-q;
 c=intv+w;
 d=intv-w;
 if(debug) print *, 'check',iswon
 if(debug) print *, 'row a b c d', row,a,b,c,d
 if (a<0 .or. a>24) then
    a=25;
 end if
 if (b<0 .or. b>24) then
   b=25;
 end if
 if (c<0 .or. c>24) then
   c=25;
 end if
 if (d<0 .or. d>24) then
   d=25;
 end if
 kc = (Int(c/q)+w)
 kd = (Int(d/q)+w)
 kv = row;
 if (kc /= kv) then
   c=25;
 end if
 if (kd /= kv) then
   d=25;
 end if
 if (v==5) d=25;
 if(debug) print *,'a b c d second', a,b,c,d 
 if ( elements (a) == 1) then 
   elements (a)=0;
 else
   elements (a)=1;
 end if

 if ( elements (b) == 1) then 
  elements (b)=0;
 else
  elements (b)=1;
 end if

 if ( elements (c) == 1) then
   elements (c)=0;
 else 
   elements (c)=1;
end if

 if ( elements (d) == 1) then 
  elements (d)=0;
 else 
  elements (d)=1;
 end if

 if ( elements (v) == 1) then 
  elements (v)=0;
 else 
  elements (v)=1;
 end if

 if(debug) print *,'6',iswon
 call checkall()
 if(debug) print *, 'check',iswon
end subroutine

end module


program test
use lightsout
 character*2 ::ss
 integer :: ival,ilev
 if(debug) call fillall()
 if(debug) call printtable()
 if(debug) print *, "*******"
 if(debug) call clearall()
 if(debug)  call printtable()

 if(debug)  print *, "*******"
 if(debug)  call level1()
 if(debug)  call printtable()
 if(debug)  print *, "*******"
 if(debug)  call level2()
 if(debug)  call printtable()
 if(debug)  print *, "*******"
 if(debug)  call level3()
 if(debug)  call printtable() 
 if(debug)  print *, "*******"
 if(debug)  call level4()
 if(debug)  call printtable() 
 if(debug)  print *, "*******"
 if(debug)  call level5()
 if(debug)  call printtable() 

 ival=0
 ilev=1
 call level1()
 
 do while(.true.)
  if(debug)  print *,iswon,ilev
  if(iswon .ne. 0) then
    ilev=ilev+1
    print *, "Entering level :" ,ilev
    if (ilev == 2) call level2
    if (ilev == 3) call level3
    if (ilev == 4) call level4
    if (ilev == 5) call level5
  end if
  
  call printtable()
   if(debug) print *,'2',iswon,ilev
  write(*,'(a1)',advance="no") ":"
  read(*,'(2a)') ss
  call givevalue(ss,ival)
  if(debug)  print *, ival
  if (ival .eq. 999) exit
   if(debug) print *,'3',iswon,ilev
   if(ival .ne. -1) call check(ival)
   if(debug) print *,'4',iswon,ilev

 end do
 contains
 
 subroutine givevalue(sstext,value)
  character*2 :: sstext
  integer :: value,i,k
  
  if (sstext(1:1) .eq. 'q') then
    value=999
  end if
  
  select case(sstext(1:1))
  
   case('a')
   i=1
   case ('b')
   i=2
   case ('c')
   i=3
   case ('d')
   i=4
   case ('e')
   i=5
   case default 
   i=-1
   value=-1
  end select
  
  
  read(sstext(2:2),'(i1)',err=100) k
  
    
  
  if (k <7 .and. k .ne. 0) then
      value=(i-1)+(k-1)*5
  else
    value=-1
    k=-1
  end if
  if(debug) print *, k,i,value
  return  
100  value =999

 end subroutine

      
end program


Fortran from python – Python Binding for Fortran Code

fortran from python
As an aeronautical engineer, I received fortran as a legacy. From NAL to working for Rolls-Royce, work life revolves around fortran. This  is the language that I have been programming for 12+ years now. That explains this many fortran related posts here and here.

But from couple of years, I am falling in love with Python. Here in this post, I will walk through the steps to use fortran and python together.

There are many times, when we don’t want to use the methods included in Python and need to use the trusted, optimised, time tested fortran library. This post will show how to create a python binding of the fortran code and use it in python.

I will present the complete code first and then explain all the parts of it.

I will use this simple NACA code library.

naca4_sym(c,t,n) evaluates y(x) for a NACA symmetric 4-digit airfoil. C the chord length. T the maximum relative thickness. N the number of sample points.

file name: naca.f

This is the actual fortran code to be used from python.

subroutine naca4_symmetric ( t, c, n, x, y )

!*****************************************************************************80
!
!! NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 21 May 2014
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Eastman Jacobs, Kenneth Ward, Robert Pinkerton,
! "The characteristics of 78 related airfoil sections from tests in
! the variable-density wind tunnel",
! NACA Report 460, 1933.
!
! Parameters:
!
! Input, real ( kind = 8 ) T, the maximum relative thickness.
!
! Input, real ( kind = 8 ) C, the chord length.
!
! Input, integer ( kind = 4 ) N, the number of sample points.
!
! Input, real ( kind = 8 ) X(N), points along the chord length.
! 0.0 <= X(*) <= C.
!
! Output, real ( kind = 8 ) Y(N), for each value of X, the corresponding
! value of Y so that (X,Y) is on the upper wing surface, and (X,-Y) is on the
! lower wing surface.
!
implicit none

integer ( kind = 4 ) n

real ( kind = 8 ) c
real ( kind = 8 ) t
real ( kind = 8 ) x(n)
real ( kind = 8 ) y(n)

y(1:n) = 5.0D+00 * t * c * ( &
0.2969D+00 * sqrt ( x(1:n) / c ) &
+ (((( &
– 0.1015D+00 ) * ( x(1:n) / c ) &
+ 0.2843D+00 ) * ( x(1:n) / c ) &
– 0.3516D+00 ) * ( x(1:n) / c ) &
– 0.1260D+00 ) * ( x(1:n) / c ) )

return
end

file name: aefoilnaca.py
This is the python wrapper using the nacalib created from fortran code


def naca4_sym(c,t,n):
	"""
	NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
	C the chord length.
	T the maximum relative thickness.
	N the number of sample points.

	python wrapper for the code by  John Burkardt

	by Sukhbinder Singh

	"""
	import numpy as np
	import nacalib

	x = np.linspace(0.0,c,n)
	y=np.zeros(n)
	y=nacalib.naca4_symmetric(t,c,x,n)
	xy=np.zeros((2,2*n))
	xy[0,0:n]=x
	xy[0,n:]=x[n::-1]
	xy[1,0:n]=-y
	xy[1,n:]=y[n::-1]
	return xy

file name: test.py
This file tests the python wrapper module.


"""
Testing aefoilnaca

"""

import aefoilnaca
import matplotlib.pyplot as plt

c = 10.0 # Chord length
t=0.15   # the maximum relative thickness
n=80     # number of point
m=0.02   # the maximum camber
p=0.4    # the location of maximum camber

y = aefoilnaca.naca4_sym(c,t,n)

plt.plot(y[0,:],y[1,:])
plt.show()

The steps:
First, we create a signature file from naca.f by running

f2py naca.f -m nacalib -h naca.pyf

The signature file is saved to naca.pyf

Second step edit naca.pyf file.

Next we’ll open the naca.pyf and tell F2PY that the argument t,c,x,n are input arguments (use intent(in) attribute) and that the result, i.e. the contents of y after calling Fortran function naca4_symmetric should be returned to Python (use intent(out) attribute). In addition, an array a should be created dynamically using the size given by the input argument n (use depend(n) attribute to indicate dependence relation).

This is the content of the modified naca.pyf

file name: naca.pyf


!    -*- f90 -*-
! Note: the context of this file is case sensitive.

python module nacalib
    interface
        subroutine naca4_symmetric(t,c,n,x,y)
            real(kind=8),intent(in) :: t
            real(kind=8),intent(in) :: c
            integer(kind=4), optional,check(len(x)>=n),depend(x) :: n=len(x)
            real(kind=8) dimension(n),intent(in) :: x
            real(kind=8) dimension(n),depend(n),intent(out) :: y
        end subroutine naca4_symmetric
    end interface
end python module nacalib

! This file was auto-generated with f2py (version:2).
! See http://cens.ioc.ee/projects/f2py2e/

And finally, we build the extension module by running

f2py -c naca.pyf naca.f

This creates the naca.so file and it can be used as shown in the aefoilnaca.py

	import nacalib

	x = np.linspace(0.0,c,n)
	y=np.zeros(n)
	y=nacalib.naca4_symmetric(t,c,x,n)

Hope this is clear and simple, if you have any questions, feel free to comment.

All this code is available at my github account

For further reading refer to this two excellent reference
1. Python as a glue
2. Using f2py – The smart way

A Simple Fortran code to Simulate Stochastic Stock Price Movement

Stochastic_Stock_fortran_ProgramA few months ago I posted this pgplot powered fortran animation program on calculating the approximate value of PI by Monte Carlo method using random numbers.

And the other found myself reading about random or stochastic phenomenons. After reading about stochastic processes, I was invariably reminded of the stock market. The most stochastic place one can imagine.

After a bit of googling and a bit of tinkering, created the following program.

Given the annual drift, volatility, initial price and number of trials, the program simulates the movement of a stock price.


     program MonteCarloStocks
     double precision a,b,c
     character(len=100) buffer
     integer i,n,noarg,IARGC
     
     noarg=IARGC()

     if(noarg == 4) then     
       CALL GETARG (1, buffer)
       read(buffer,'(f6.2)') a
       CALL GETARG (2, buffer)
       read(buffer,'(f6.2)') b
       CALL GETARG (3, buffer)
       read(buffer,'(f6.2)') c
       CALL GETARG (4, buffer)
       read(buffer,'(i7)') i
       print *, a,b,c,i
     else
       write(*,'(a)',advance='no') "Please enter annual drift :"
       read(*,'(f8.2)')a
       write(*,'(a)',advance='no') "Please enter annual Volatility :"
       read(*,'(f8.2)')b
       write(*,'(a)',advance='no') "Please enter Initial Price :"
       read(*,'(f8.2)')c
       write(*,'(a)',advance='no') "Please enter no of iteration :"
       read(*,'(i7)')i
     end if
     call monte(a,b,c,i)
     call system("stock.xls")
     
     end program
     
     subroutine monte(driftyearly,volatility,initialprice,not)
     double precision driftyearly,volatility,initialprice
     double precision driftdaily,voldaily,p,pri(252),PPND16
     integer,parameter :: nod=252
     integer not
     
     driftyearly = driftyearly/100
     volatility  = volatility/100
     
     driftdaily=driftyearly/nod
     voldaily=volatility/(nod**0.5)
     open(200,file="fort.200")   
     call random_seed
     if(not <= 0) not=1
     if(not >=25000) not =25000
     do j=1,not
       pri(1)=initialprice
        do i=2,nod
         call random_number(p)
         p=PPND16 (p, IFAULT)
         pri(i)=pri(i-1)*EXP(driftdaily+voldaily*p)
         end do
         write(200,'(252f8.2)')pri
     end do    
     close(200)
     end subroutine
     
 

	DOUBLE PRECISION FUNCTION PPND16 (P, IFAULT)
!
!	ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
!
!	Produces the normal deviate Z corresponding to a given lower
!	tail area of P; Z is accurate to about 1 part in 10**16.
!
!	The hash sums below are the sums of the mantissas of the
!	coefficients.   They are included for use in checking
!	transcription.
!
	DOUBLE PRECISION ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1,   &
     		CONST2, A0, A1,	A2, A3, A4, A5, A6, A7, B1, B2, B3, &
                B4, B5, B6, B7,                                 &
     		C0, C1, C2, C3, C4, C5, C6, C7,	D1, D2, D3, D4, D5, &
     		D6, D7, E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, &
     		F4, F5, F6, F7, P, Q, R
	PARAMETER (ZERO = 0.D0, ONE = 1.D0, HALF = 0.5D0,           &
     		   SPLIT1 = 0.425D0, SPLIT2 = 5.D0,                 &
     		   CONST1 = 0.180625D0, CONST2 = 1.6D0)
!
!	Coefficients for P close to 0.5
!
	PARAMETER (A0 = 3.3871328727963666080D0, &
     		   A1 = 1.3314166789178437745D+2,&
		A2=1.9715909503065514427D+3,&
		A3=1.3731693765509461125D+4,&
		A4=4.5921953931549871457D+4,&
		A5=6.7265770927008700853D+4,&
		A6=3.3430575583588128105D+4,&
		A7=2.5090809287301226727D+3,&
		B1=4.2313330701600911252D+1,&
		B2=6.8718700749205790830D+2,&
		B3=5.3941960214247511077D+3,&
		B4=2.1213794301586595867D+4,&
		B5=3.9307895800092710610D+4,&
		B6=2.8729085735721942674D+4,&
		B7=5.2264952788528545610D+3)
!	HASH SUM AB    55.88319 28806 14901 4439
!
!	Coefficients for P not close to 0, 0.5 or 1.
!
	PARAMETER (C0 = 1.42343711074968357734D0,&
		C1=4.63033784615654529590D0,&
		C2=5.76949722146069140550D0,&
		C3=3.64784832476320460504D0,&
		C4=1.27045825245236838258D0,&
		C5=2.41780725177450611770D-1,&
        C6=2.27238449892691845833D-2,&
		C7=7.74545014278341407640D-4,&
		D1=2.05319162663775882187D0,&
		D2=1.67638483018380384940D0,&
		D3=6.89767334985100004550D-1,&
		D4=1.48103976427480074590D-1,&
		D5=1.51986665636164571966D-2,&
		D6=5.47593808499534494600D-4,&
		D7=1.05075007164441684324D-9)
!	HASHSUMCD49.33206503301610289036
!
!	CoefficientsforPnear0or1.
!
	PARAMETER(E0=6.65790464350110377720D0,&
		E1=5.46378491116411436990D0,&
		E2=1.78482653991729133580D0,&
		E3=2.96560571828504891230D-1,&
		E4=2.65321895265761230930D-2,&
		E5=1.24266094738807843860D-3,&
		E6=2.71155556874348757815D-5,&
		E7=2.01033439929228813265D-7,&
		F1=5.99832206555887937690D-1,&
		F2=1.36929880922735805310D-1,&
		F3=1.48753612908506148525D-2,&
		F4=7.86869131145613259100D-4,&
		F5=1.84631831751005468180D-5,&
		F6=1.42151175831644588870D-7,&
		F7=2.04426310338993978564D-15)
!	HASH SUM EF    47.52583 31754 92896 71629
!
	IFAULT = 0
	Q = P - HALF
	IF (ABS(Q) .LE. SPLIT1) THEN
	  R = CONST1 - Q * Q
	  PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3)  &
     			* R + A2) * R + A1) * R + A0) /  &
     		      (((((((B7 * R + B6) * R + B5) * R + B4) * R + B3)  &
     			* R + B2) * R + B1) * R + ONE)
	  RETURN
	ELSE
	  IF (Q .LT. ZERO) THEN
	    R = P
	  ELSE
	    R = ONE - P
	  END IF
	  IF (R .LE. ZERO) THEN
	    IFAULT = 1
	    PPND16 = ZERO
	    RETURN
	  END IF
	  R = SQRT(-LOG(R))
	  IF (R .LE. SPLIT2) THEN
	    R = R - CONST2
	    PPND16 = (((((((C7 * R + C6) * R + C5) * R + C4) * R + C3)  &
     			* R + C2) * R + C1) * R + C0) /  &
     		     (((((((D7 * R + D6) * R + D5) * R + D4) * R + D3) &
     			* R + D2) * R + D1) * R + ONE)
	  ELSE
	    R = R - SPLIT2
	    PPND16 = (((((((E7 * R + E6) * R + E5) * R + E4) * R + E3)  &
     			* R + E2) * R + E1) * R + E0) /  &
     		     (((((((F7 * R + F6) * R + F5) * R + F4) * R + F3)  &
     			* R + F2) * R + F1) * R + ONE)
	  END IF
	  IF (Q .LT. ZERO) PPND16 = - PPND16
	  RETURN
	END IF
	END

Cool I would say.

An example of Monte Carlo method in action.

Knight in Fortran

Knight Move Game in Fortran 90
This is a brain game that you might enjoy. In the chess board you can move the knight as if you’re playing a regular game of chess, but you need to cover the whole board with the least possible moves, a good score would be 65 moves or less. This is a very good training for strategy games.

Written in fortran. Just for fun!! As always work in progress..

Continue reading