Feap installation with visual studio

If you ever require a FEM system for quick fem for educational or research purpose. Feap is one of the easiest to get started with.

Here’s a rundown with screenshots on how to build it with visual studio with Intel Fortran.

Hope this helps.

Step 0:

Download from:

http://projects.ce.berkeley.edu/feap/feappv/feappv41.zip

Project Page:

http://projects.ce.berkeley.edu/feap/feappv/

Two steps

  1. Build a library
  2. Build the program
  1. Select New Project

  1. Select New Project
    1. Select Library: Select Static Library
    2. Name library e.g. lib22

  1. Under the Projects tab select Add Existing Item
    1. Add all subroutines in directories: Elements, Plot, Program, User, and Windows (do not include Unix, Include, or Main).

Under the Projects tab select Properties

Select Fortran then General

b.) Set additional include path to point to the feappv include

directory (e.g. c:\users\xxx\feappv\ver22\include) and the appropriate

    directory for 32-bit or 64-bit pointers (e.g.

    c:\users\xxx\feappv\ver22\include\integer4 or

    c:\users\xxx\feappv\ver22\include\integer8)

Build library

Building the Main program

Open Visual Studio or if open new project

a.) Select QuickWin Application, select QuickWin option

(not standard graphics QuickWin)

b.) Name main program e.g. feappv

At top select Release build (as opposed to Debug)

3. Under the Projects tab select Add Existing Item

a.) Set show all files in window and add library (e.g. lib22)

Visual Studio normally places this in

c:\users\xxx\documents\visual studio\projects\lib22\lib22\release

b.) Add feappv.f from the subdirectory Main.

4. Under the Projects tab select Properties

a.) Select Fortran then General

b.) Set additional include path to point to the feappv include

directory (e.g. c:\users\xxx\feappv\ver22\include) and the appropriate

    directory for 32-bit or 64-bit pointers (e.g.

    c:\users\xxx\feappv\ver22\include\integer4 or

    c:\users\xxx\feappv\ver22\include\integer8)

Add the libraries and the library path

Build… If you get error like this, you have not included all forttan file sin the liberay include and come pback..



Updated Binary STL File Reader

bent_plate

A recent comment on the post titled binary stl file reader in numpy prompted me to revisit the code and while I was looking at it , I updated the same for code to view the loaded stl file.

The above picture shows the result.

The code as usual available on github from here

from binarySTLreader import BinarySTL,ShowSTLFile
h,p,n,v1,v2,v3=BinarySTL('bent_plate.stl')
ShowSTLFile(v1,v2,v3)

OD Command in Windows with Python

I have seen OD command for last two years but had never bothered to look it up. It always remained under the surface.

Life went on well without it, until now.

During a recent debugging crisis, had to resort to this command for an efficient way to check binary results. All the while had to go to Linux and copy the files. This is when I knew I have to see what this OD is?

OD is octal dump. Simple functionality, don’t know if there is a windows equivalent. So before searching for a windows equivalent, I knew I can replicate this in python and here’s my first attempt.

Works for me and suits my requirement. Not an entire mimic of actual OD command but if you are looking for long and double precision values, this does the job.

I have kept the command line arguments similar to the actual OD command. Here’s a simple demo.

As always, the code is available from my github repo here. Hope to add more functionality as and when I get more time. Please free to fork and play around.

Special thanks to Andy for introducing the command and helping me make sense of the important arguments.

Quick Gantt Chart with Matplotlib

ProjectPlan_GANTT_CHart_MatplotlibThe problem: You need a quick Gantt chart for a quick proposal report and you dint have any project planner software installed.

Solution:

While there are many different ways, you only have access to python, well Here’s a simple Gantt chart plotter with just matplotlib.

Not for rigorous use but a good substitute to make quick Gantt plot for quick report.

Continue reading

Quick Intro to 2D and 3D Visualization in Python

contour_visualization_with_python
Posted a quick tutorial introducing visualization in python.

This tutorial offers a quick look at some common 2d and 3d Visualization available in python. This tutorial is meant to give anyone a quick start in matplotlib and mayavi.mlab.

The tutorial was written as part of CSRP program of AeSIAA which I am volunteering as a mentor to few aeronautical students.

View the tutorial by clicking this Quick Intro to 2D and 3D Visualization in Python

Analysing Trusses – a python program

truss_in_python

Was digging into my laptop and found this Truss program written in python. Computes, displacement, stresses and reactions.


"""
Created on Thu May 08 07:07:24 2014

@author: Sukhbinder Singh

Truss FEM

"""
from __future__ import division
import numpy as np

modE=30.0e6
Area=2.0


#numElem=3
#numNodes=4

elemNodes=np.array([[0,1],[0,2],[0,3]])
nodeCords=np.array([[0.0,0.0],[0.0,120.0],[120.0,120.0],[120.0,0]])

#'''
elemNodes=np.array([[0,1],[0,2],[1,2],[1,3],
                    [0,3],[2,3],[2,5],[3,4],[3,5],[2,4],[4,5]])

nodeCords=np.array([
[0.0,0.0],[0.0,3000.0],
[3000.0,0.0],[3000.0,3000.0],
[6000.0,0.0],[6000.0,3000.0]
])


modE=70000
Area=300
#'''

numElem=elemNodes.shape[0]
numNodes=nodeCords.shape[0]

xx=nodeCords[:,0]
yy=nodeCords[:,1]

EA=modE*Area
tdof = 2*numNodes #total number of degrees of freedom
disps=np.zeros((tdof,1))
force=np.zeros((tdof,1))
sigma=np.zeros((numElem,1))
stiffness=np.zeros((tdof,tdof))
np.set_printoptions(precision=3)

'''
force[1]=-10000.0
presDof=np.arange(2,9)

'''
force[3]=-50000.0
force[7]=-100000.0
force[11]=-50000.0

presDof=np.array([0,1,9])
#'''

for e in xrange(numElem):
    indice= elemNodes[e,:]
    elemDof=np.array([indice[0]*2, indice[0]*2+1, indice[1]*2, indice[1]*2+1 ])
    xa=xx[indice[1]]-xx[indice[0]]
    ya=yy[indice[1]]-yy[indice[0]]
    len_elem=np.sqrt(xa*xa+ya*ya)
    c=xa/len_elem
    s=ya/len_elem
    k1=(EA/len_elem)* np.array([[c*c,c*s,-c*c, -c*s],
                                [c*s,s*s,-c*s ,-s*s],
                                [-c*c,-c*s,c*c,c*s],
                                [-c*s,-s*s,c*s,s*s]])
    stiffness[np.ix_(elemDof,elemDof)] +=k1


actDof=np.setdiff1d(np.arange(tdof),presDof)

disp1=np.linalg.solve(stiffness[np.ix_(actDof,actDof)],force[np.ix_(actDof)]);
disps[np.ix_(actDof)]=disp1


# stresses at elements

for e in xrange(numElem):
    indice= elemNodes[e,:]
    elemDof=np.array([indice[0]*2, indice[0]*2+1, indice[1]*2, indice[1]*2+1 ])
    xa=xx[indice[1]]-xx[indice[0]]
    ya=yy[indice[1]]-yy[indice[0]]
    len_elem=np.sqrt(xa*xa+ya*ya)
    c=xa/len_elem
    s=ya/len_elem
    sigma[e] = (modE/len_elem) * np.dot(np.array([-c,-s,c,s]),disps[np.ix_(elemDof)])

        
print disps
print sigma

react = np.dot(stiffness,disps)
print react.reshape((numNodes,2))



Map Your Google Location History in 3D

Map Your Google Location History in 3D

Inspired by this post “Map Your Google Location History”, I re-mashed the excellent tutorial and added a new function to view google location history on a 3D globe.

Basemap is cool but VTK visualisation toolkit is cooler. Check the results.

And here’s the python notebook to generate the above result.

Cool isn’t it?

Voronoi with Python

Voronoi with Python

A few weeks back, a colleague was searching for matlab for a task to get the Voronoi diagram from some points he had. I suggested why not use #python? He was not convinced and he got hold of matlab license and did his job.

But I had this itch for trying Voronoi in python, so here it is. Finally got some spare time. It was so easy in python, that I built the example around a tkinter gui.

But before we go to code, for everyone’s benefit What is Voronoi Diagram?

In mathematics, a Voronoi diagram is a way of dividing space into a number of regions. A set of points (called seeds, sites, or generators) is specified beforehand and for each seed there will be a corresponding region consisting of all points closer to that seed than to any other. The regions are called Voronoi cells. It is dual to the Delaunay triangulation.

via wikipedia

And here’s the code.

Continue reading

Fun Experiments with Markov Chain

independancedayA stochastic process is a mathematical model that evolves over time in a probabilistic manner. Markov chain is a special kind of stochastic process, where the outcome of the experiment depends only on the outcome of the previous experiment. That is the next state of the system depends only on the present state, not on preceding states. Markov chains are named after the Russian mathematician A. A. Markov who stated the theory of stochastic processes.

Now what if we mix markov chain and PM Modi’s independence day speech, find the hilarious results in this ipython notebook.

Markov Chain and Modi’s Speech

After this experiment, I know how some of the comment spam is being generated and is fooling the comment filters on wordpress.

Exploring AirWire with Pandas – A Dashboard from Session Histroy

airwire internet usage analyticsLast week received my internet bill and my wife wondered how much data we are using and if someone was snooping at our wifi?

Assured her, this isn’t the case but her doubt lingered.

Went to the ISP website, downloaded the session data and showed her. Still not convinced, so imported the data into pandas.

And that exploration.. produced the above PDF. A dashboard of information mined from the csv file. My wife loved it.

Here’s the script that will produce this pdf. Will be directly useful to anyone using airwire services and for other’s it’s a good starting point.

Continue reading

Little Experiments with BeautifulSoup

cricketwithball
Look at any conference related to python, web scraping is one of the topic that is always there.

This prompted me to BeautifulSoup. I am toying with it for few months now, but yesterday got a real chance to use it.

Here’s a small function that I used yesterday to scrape out all the test match and one day match results from espn website. Someone needed the data.


import urllib2
from bs4 import BeautifulSoup
import csv
def cricketScrap(url,f):
    page = urllib2.urlopen(url).read()
    soup = BeautifulSoup(page)
    table = soup.find('table', attrs={'class': 'engineTable'})
    x = (len(table.findAll('tr')) )
    if x==1:
        return
    
    for row in table.findAll('tr')[1:x]:
        col = row.findAll('td')
        team1 = col[0].getText()
        team2 = col[1].getText()
        winner = col[2].getText()
        margin = col[3].getText()
        ground = col[4].getText()
        match_date = col[5].getText()
        match = (team1, team2, winner, margin, ground, match_date )
        f.writerow(match)

Here’s how to use it.

'''
For one day match

'''

lis=["Team 1","Team 2","Winner","Margin","Ground","Match Date"]
f = csv.writer(open("oneindiancricketscrape.csv", "w"))
f.writerow(lis)

onedayurl=r"http://stats.espncricinfo.com/indian-premier-league-2014/engine/records/team/match_results.html?class=2;team=6;type=year;id="
for year in range(1974,2015):
    cricketScrap(onedayurl+str(year),f)


And here’s the interesting one, scape all results of the matches played between 1772 to 2014


allmatchurl=r"http://stats.espncricinfo.com/indian-premier-league-2014/engine/records/team/match_results.html?class=0;team=0;type=year;id="

lis=["Team 1","Team 2","Winner","Margin","Ground","Match Date"]
f = csv.writer(open("allCricketMatchResults1900.csv", "w"))
f.writerow(lis)


for year in range(1772,2015):
    cricketScrap(allmatchurl+str(year),f)


And if you are just interested in one day and test matches results played by india till now, then click the following links to download the csv files.

Foundation of FEM

Spring_FEM_in_Python
“Spring is the foundation of FEM. Understand this well and you are through in this course.” said the IISC professor as we settled on the hard benches on the first day of our 6 month FEM course at IISC, bangalore.

The course proved very useful. Learnt a lot. Fortunately , being in NAL had access to matlab and so converted the learning into matlab scripts. Today those scripts must be lying in some system in NAL.

Two months back, working on a problem in office, revisited all those notes and had this impulse to revisit the spring FEM problem in python…

So here’s the result. Two problem defined and solved in python.


"""
Created on Tue April 06 18:46:01 2014

@author: Sukhbinder Singh

Finite Element Analysis of spring system

problem 1

               #---^^^^^----#
#              #
#-----^^^^^----# ===> P
#              #
               #---^^^^^----#

1              2            3,4     

Spring Problem

P=100
ki=10
u1=u3=u4 =0


problem 2

#                             #
#------^^^^^^^^^^^^^^^^-------#
#                             #
#         |---^^^^--|         #
#---^^^^--|==>P     |--^^^^---#
#         |---^^^^--|         #
          |                   #
          |------^^^^^^-------#

1         2         3         4

k1=10,k2=15,k3=20,k4=25,k5=30,k6=35
P=100
u1=u4=0


Solution for problem 1
Displacement
0 0.0
1 3.33333333333
2 0.0
3 0.0

Reactions
0 -33.3333333333
2 -33.3333333333
3 -33.3333333333


Solution for problem 2
Displacement
0 0.0
1 1.54589371981
2 0.869565217391
3 0.0

Reactions
0 -23.1884057971
3 -76.8115942029


#usage of np.ix_
import numpy as np
A = np.array([ 1, 2, 3, 4, 4,
      2, 3, 4, 5, 3,
      4, 5, 6, 7, 2,
      5, 6, 7, 8, 9,
      6, 7, 8, 9, 0 ]).reshape(5, 5)

B = np.array([60, 70, 80, 90]).reshape(2, 2)
A[np.ix_([2, 4], [2, 4])] = B  


# usage of np.setdiff1d(ar1,ar2)
# Return the sorted, unique values in `ar1` that are not in `ar2`.
a = np.array([1, 2, 3, 2, 4, 1])
b = np.array([3, 4, 5, 6])
np.setdiff1d(a, b)

"""
import numpy as np

p1 = True # True for problem 1

if p1: # problem 1
    elemNodes= np.array([[0,1],[1,2],[1,3]])
    numElem = elemNodes.shape[0]
    numNodes=4
    presDof=np.array([0,2,3])
    k=np.array([10,10,10]) # Spring Stiffness
else: # problem 2
    elemNodes= np.array([[0,3],[0,1],[1,2],[1,2],[1,3],[2,3]])
    numElem = elemNodes.shape[0]
    numNodes=4
    presDof=np.array([0,3])
    k=np.array([10,15,20,25,30,35]) # Spring Stiffness

disp=np.zeros((numNodes,1))
force=np.zeros((numNodes,1))
stiffness=np.zeros((numNodes,numNodes))

force[1]=100.0

for e in xrange(numElem):
    stiffness[np.ix_(elemNodes[e,:],elemNodes[e,:])]+=np.array([[k[e],-k[e]],[-k[e],k[e]]])

actDof=np.setdiff1d(np.arange(numNodes),presDof)

disp1=np.linalg.solve(stiffness[np.ix_(actDof,actDof)],force[np.ix_(actDof)]);
disp[np.ix_(actDof)]=disp1

print "\nDisplacement"
for i,d in enumerate(disp):
    print i,d[0]

# reactions
react = np.dot(stiffness,disp)
reaction = react[presDof]

print "\nReactions"
for (p,r) in zip(presDof,reaction):
    print p,r[0]


How to Read HDF5 with Python – a simple introduction

HDF keeps coming into my world from one place or the other. A colleague is working on it, then it showed up in a side project i am working and now I stumbled on a video praising hdf5. So finally decided to dig into it and find more about it.

  1. HDF stands for Hierarchical Data Format. This (HDF, HDF4, or HDF5) is a set of file formats and libraries designed to store and organize large amounts of numerical data. Originally developed at the National Center for Supercomputing Applications, it is supported by the non-profit HDF Group.
  2. HDF5 format supports any kind of data for digital storage regardless of their origin and some.
  3. HDF5 stores data in a highly organised and hierarchical format
  4. HDF5 is highly efficient chunked input output operations.
  5. HDF5 allows inclusion of metadata and attribute.
  6. The format is platform independent and is widely used in scientific codes.
  7. HDF5 data model consists of two major objects, Datasets and Groups.

Enough theory? Here’s a super simple introduction to reading hdf5 data via python.

pytables_intro_how_to_read_hdf5_with_python

enjoy!

Read more about HDF5 here [pdf] and here

Simple functions to collect data from excel files using xlrd and python

excel_as_a_form_use_xlrd_to_read_data_from_excelIt’s no secret. Excel is abused and misused in so many different ways and the most rough abuse is disregarding Excel’s inherent 2D table structure and using it to collect data as an excel based form.

Recently had a task where I faced two hundred or so such forms, from where needed to pull the data. Turned to my favourite hammer, Python.

And with xlrd, this was as easy as cutting a cake. Here’s three functions that did the job.

GetFilenames collects all the filenames including path given a top level folder. Getdata collects the required data from the given row, col and writeOut dumps all the collected data into a csv file as a 2D table.

from os import walk
from os.path import join
def GetFilenames(mypath):
    return (join(dirpath,f) for (dirpath, dirnames, filenames) in walk(mypath) for f in filenames if f.endswith(('.xls','.xlsx')))

import xlrd
def getdata(fname,listx,listy):
"""
Reads data from excel sheet

fname excel file name
listx=[2,4,6,8,10,12,14,17,19,21,27,29,2,6,8,10,4] # Rows index
listy=[1,1,1,1,1,1,1,1,1,1,1,1,4,4,4,4,4]          # Column index

from which data has to be extracted

returns dataList containing the extracted info from excel.

"""
    dataList=[]   # initial an empty list
    book=xlrd.open_workbook(fname) # open workbook
    sheet=book.sheet_by_index(0) # Get first sheet
    dataList.append(fname) 
    for x,y in zip(listx,listy): # walk to all the row,col combination and collect data
        if sheet.cell_type(x,y) <> 3: # if its not a data data
           dataList.append(sheet.cell_value(x,y))
        else:                         # if its a date data
          datetup =xlrd.xldate_as_tuple(sheet.cell_value(x,y),book.datemode)
          dataList.append(datetup[0])
          dataList.append(datetup[1])
          dataList.append(datetup[2])
    return dataList


import csv
def writeOut(fname,files):
    f = open(fname,"wb")
    writer=csv.writer(f)
    for f in files:
        writer.writerow(getdata(f))
    print fname,"written !!"    
    return None

My Netas with PANDAS

Election_2014_Candidates_pandasWas surfing internet and found this affidavit data filled by candidates of election 2014 in a convenient comma separated file format.

Couldn’t resist myself. The tinkerer in me wanted to play with this data, so tinkering began.

#
# Data Source https://github.com/datameet/india-election-data
#
import numpy as np
import pandas as pd
import matplotlib

%matplotlib inline

fname=r"/Users/Sukhbinder/myneta2014.csv"

data =pd.read_csv(fname)
data.describe()

Using Pandas library, loaded the data and first two results brought some surprises. So spent the next couple hours digging the data.

Ok answer this

  • Which party has the highest count of registered criminal cases ?
  • Which party do the top two most criminal cases candidate belong to?

Answer before you proceed. I was surprised by the answer this data revealed!

I also found answers to these questions.

  • How many candidates were running in the election?
  • How many parties were contesting elections?
  • What is the Education Status of candidates in different party?
  • Who is the richest contestant by total assets?

and many more.

Find out. And while you are at it, you might learn some pandas. I am still a beginner in pandas and might not have done full justice to the awesome pandas library!!

Here’s the notebook to check the code and the results. And if you prefer PDF, then click here.

Visit my github to play with the data and the code.

Simple and Minimal File Selection GUI with Standard Python

I have a simple script which parses an excel file and does some data analysis on it’s contents and reports it as a html output. I don’t need a UI, but right now I’m prompting the user for the file to parse using raw_input which is very inconvenient. So needed to find a quick and easy way to present a file selection dialog to the user with which they can select the excel file.

My first choice was wx. But then realised i needed something thats part of standard python and that landed me to this page File Tkinter dialogs (Python recipe)

Here’s the super convenient and simple way to get a minimal file selection GUI from standard python modules.

import Tkinter,tkFileDialog
root = Tkinter.Tk()
root.withdraw()
filename = tkFileDialog.askopenfile(parent=root,mode='r',filetypes=[(“Excel file”,”*.xls”)],title='Choose an excel file')
if filename != None:
    print “This excel file has been selected”, file

Thanks internet!!

And here’s a good presentation on Tkinter titled “Tkinter Does Not Suck”. I agree…

Python App to perform ICP alignment finally out

Raptor: Python ICP with VTK and WX PythonBefore reading further, if you haven’t seen the demo yet, please visit this page. Go now! there’s a Raptor waiting for you. Done that.. good..

As everybody and everything is in github now a days, so have posted the code in github at this url. PythonICPvtk Check it out.

Work in progress and will try to update it soon. Have lots of todo’s for this code.

Go have a look, just 470+ odd lines of python code including comments and you will find the Raptors too!!