Exploring MFEM

With little time I have after office and playing with kids , I am devoting some to this MFEM library.

Really good piece of software. Exploration continues.

So What is MFEM?

No, it’s not Ministry of Finance and Economics Management, it’s …….

MFEM is a lightweight, general, scalable C++ library for finite element methods and available at http://mfem.googlecode.com

The goal of MFEM is to enable research and development of scalable
finite element discretization and solver algorithms through general
finite element abstractions, accurate and flexible visualization, and
tight integration with the hypre linear solvers library. Its features
include:
– 2D and 3D, arbitrary high-order H1, H(curl), H(div), L2 and NURBS
elements.
– Parallel version scalable to hundreds of thousands of MPI cores.
– Conforming or nonconforming adaptive mesh refinement (AMR),
including anisotropic refinement.
– Galerkin, mixed, isogeometric, DG and DPG discretizations.
– Support for triangular, quadrilateral, tetrahedral and hexahedral
elements with curved boundaries.
– Lightweight interactive OpenGL visualization with GLVis,
http://glvis.googlecode.com.

An interactive documentation of MFEM’s serial and parallel example
codes can be found here
MFEM is freely available under LGPL 2.1.

 

Advertisement

FEM in Python A Simple Start Guide

Wrote this a couple of months back. Yet another tutorial in python, if you are interested in finite element analysis. Nicely goes with this excellent tutorial on FEM.

Click the below link to view the tutorial. FEM with Python

Displaying Truss

The truss program that was posted two weeks back, looked incomplete without a proper visualisations, so here’s a small VTK function written to display the truss with any quantity like stress/reactions or displacement.

Foundation of FEM

“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]

Brief Intro to FEA or Finite Element Method (FEM)

Recently someone asked me for a simple and short primer or refresher material on Finite element Method.

Here’s a super simple, brief, to the point , no nonsense PDF on FEA theory from the university of victoria.

Brief and simple. Just enough material to spark interest for the curious!!!

Brief Intro to FEA or Finite Element Method (FEM) [PDF]

Do you know any such good, brief and useful resource? Please share.

Unofficial History of ANSYS

A few days back I was searching for a paper and stumbled upon this interesting page. The page listed an interesting account of unofficial history of ANSYS. I don’t know why, but i like reading accounts like this, so here’s part of it that i liked…

In 1963, Dr. John Swanson worked at Westinghouse Astronuclear Labs in Pittsburgh, responsible for stress analysis of the components in NERVA nuclear reactor rockets. He used computer codes to model and predict transient stresses and displacements of the reactor system due to thermal and pressure loads.

Swanson continued to develop 3-D analysis, plate bending, nonlinear analysis for plasticity and creep, and transient dynamic analysis, in the next several years, using a finite element heat conduction program that was developed by Wilson at Aerojet. The old Westinghouse codes included a 2d/axisymmetric one also, possibly called FEATS (according to Kohnke).

John wanted to combine these codes to remove the duplication, like equation solvers and some postprocessing Swanson believed an integrated, general-purpose FEA code could be used to do complex calculations that engineers typically did manually, such as heat transfer analysis. It would save money and time for Westinghouse and other companies.

Westinghouse didn’t support the idea, and Swanson left the company in 1969. Before he left, he made sure that all code work had been sent to COSMIC, so that he could pick it up again from the outside.

Swanson Analysis Systems, Inc was incorporated in the middle of 1970 at Swanson’s home. The offices were part of Swanson’s home (There was no garage) in Pittsburgh. at the same time, Westinghouse realized that they needed John, so they hired him as a consultant. John said sure, but with the proviso that whatever he put into STASYS, the Westinghouse code, he could also put into ANSYS. Westinghouse had no trouble with this, as they just wanted to solve their problems. So this consulting kept bread on the table for the Swansons, and at the same time brought forth further improvements to ANSYS.

He developed his program using a keypunch and a time-shared mainframe at U.S. Steel. The first version of ANSYS was coded by the end of 1970, and the ANSYS program was first leased soon after that. Westinghouse was the first customer, running as a data center. The data center was at the Telecommunications Center on the Parkway East, on the east side of  Pittsburgh.

According to Dr. Swanson, the name ANSYS was because the copyright lawyers assured Swanson that ANSYS was just a name, and did not stand for anything. Understandable, during that period all programs were "written" on punch card. When installing the program on the customer’s computer, it meant carried a relatively big case of punch cards to the customer’s place, and fed them into the machine.

Dr. Peter Kohnke met John Swanson first about early 1971. Swanson offered Peter a job in the fall of 71, but Peter did not accept. As that time Peter was a brand new father at the time, and Westinghouse looked a lot more secure than SAS. Peter told John he was interested and would accept in the spring of 72, but then he had hired Gabe DeSalvo, and did not have the resources to hire Peter also. But that fall John did, and then hired Peter. Dr. Peter Kohnke’s start date was 1/1/73. (Dr. Kohnke still worked at ANSYS as of December of2005).

When Peter started work, he asked what John wanted him to do. John told Peter that he developed code, did technical support, wrote manuals, gave seminars and did systems work so that the program would run on a variety of systems. John said he needed relief, so Peter should pick one or two of them. Ultimately, Peter did all of them, except systems work.

In around 1970, users can ran ANSYS 2.x on a CDC 6600 machine over the Cybernet timesharing network. That time only fixed format input was available. The users would work up the input listing off-line, key it onto a tape cassette, log on, submit the run about quitting time for the best computer rates and stop by the CDC data center next morning to find out what went wrong. In 1973, ANSYS ran on three kinds of hardware: CDC, Univac, and IBM. And around 1973 the USS mainframe that they used to develop code was the US Steel CDC 6500.

The first minicomputer that ANSYS ran on was a MODCOMP 4 (or IV?) VAX came later. Being a small company, everyone did everything. When a ‘mini’ computer was delivered, everyone helped wrestle it off of the truck. When printout paper was delivered, everyone helped unload the
boxes.

Lot more to read at Unofficial history of ANSYS

Plate Mesh File

A few months back, i needed a plate model. Searched the net and tried creating it with a CAD software, but then said to myself. What the heck!! All I need is a mesh with quad elements and so wrote this program.

Not a polished code, but if you need a patran neurtal mesh file  for a simple plate with custom amount of 4 noded elements, then this is the program for you. Run this and import the generated file any CAD/CAE program.

 

        IMPLICIT NONE
        INTEGER N_elts_x,D_Node_Number,i,j,N_elts_y
        REAL*8 coord(121,2),dx,dy,PlateLength,PlateHeight
        INTEGER conec(100,4),nnodes,nele,ishap,D_Node_Numbery,ielem

        CHARACTER(len=80) title,datetime
        DATA title/"Created by fortran code"/
        DATA datetime/" Developed by sukhbinder"/



        PlateLength = 0.25
        PlateHeight = 0.25
! number of elements in x dir
        N_elts_x = 3
! number of elements in y dir
        N_elts_y = 4


        dx = PlateLength / DBLE(N_elts_x)
        dy = PlateHeight / DBLE(N_elts_y)

        D_Node_Number = N_elts_x + 1
        D_Node_Numbery = N_elts_y + 1
        DO i = 1 ,(N_elts_y + 1)
          DO j =1 , (N_elts_x + 1)
            COORD((i-1)*D_Node_Number+j , 1 ) = (j-1) *  dx
            COORD((i-1)*D_Node_Number+j , 2 ) = (i-1) * dy
          END DO
        END DO
!
          ielem=0
          DO j = 1, N_elts_y
            DO i = 1, N_elts_x

              ielem = ielem + 1

              conec(ielem,1) = ( j - 1 ) * ( N_elts_x + 1 ) + i
              conec(ielem,2) = ( j - 1 ) * ( N_elts_x + 1 ) + i + 1
              conec(ielem,3) =   j       * ( N_elts_x + 1 ) + i + 1
              conec(ielem,4) =   j       * ( N_elts_x + 1 ) + i


            END DO
          END DO

        DO i=1,121
          WRITE(*,'(i,2f12.8)') i,coord(i,:)
        END DO

        DO i=1,100
        WRITE(*,'(5i)') i,conec(i,:)
        END DO

!  Writing the file

        nnodes=(N_elts_x+1)*(N_elts_y+1)
        nele=N_elts_x*N_elts_y
        ishap =4

        OPEN(13,file="c:\temp\temp.mesh")

        WRITE(13,9000) 25,0,0,1,0,0,0,0,0
        WRITE(13,'(A)') title
        WRITE(13,9000) 26,0,0,1,nnodes,nele,0,0,0
        WRITE(13,'(A)') datetime

        DO i=1, nnodes
           WRITE(13,9005) i,0,2,0,0,0,0,0,coord(i,1),coord(i,2),0.0d0
        END DO

        DO i=1,nele
           WRITE(13,9010) i, ishap, ishap+1,0,0,0,0,0, ishap, 1,1,0,0.d0,0.d0,0.d0,conec(i,1:ishap)
        END DO
        WRITE(13,9000) 99, 0, 0, 1, 0, 0, 0, 0, 0



9000    format(i2, 8i8)
9005    format( ' 1', 8i8,  /&
           1p,3e16.9 / &
             '1G       6       0       0  000000')
9010  format( ' 2', 8i8, / &
                   4I8,3e16.9 /&
                 (10i8 ))
        CLOSE(13)
        END

Meshing Fortran and Matlab together with Mex

A few days back colleague of mine needed a help in calling Fortran from matlab. At that time I knew about the call system command.

But as we tried solving his problem, I learnt about Mex. I knew Mex was matlab executable but had never explored them or used or created them.

My colleague’s problem was solved once we got matlab and the c compiler. That day and few days after, I continued exploring the dark alleys of Mex.

Learnt a lot.  Here’s the list of links that were great at explaining the meshing of fortran and c with matlab with Mex.

• Creating, compiling and linking mex files with fortran and c
• Creating Fortran MEX files from Mathworks
• Writing C Functions in MATLAB (MEX-Files)
• Writing Matlab MEX files using Fortran 90 / 95 (and the Intel Fortran Compiler)
• Interfacing Fortran routines from matlab in a simple and efficient way [pdf]

Perhaps someday I will use them in some project, but till then let them be kept here.

Meet fable and fortwarp

I met fable and fortwarp a few weeks back. I am beginning to know them and before I get too comfy with them, I thought of introducing them to the readers of this blog.

If you work with c,c++ and fortran, you will definitely like their company!

Fable
fable converts fixed-format Fortran sources to C++. The generated C++ code is designed to be human-readable and suitable for further development. In many cases it can be compiled and run without manual intervention.

fable comes with the C++ fem Fortran EMulation library.

fable is written in Python (version 2.x)

The name “fable” is short for “Fortran ABLEitung”. “Ableitung” is a german word which can mean both “derivative” and “branching off”.

FortWarp
FortWrap is a python script that parses Fortran 90/95/200X source files and generates wrapper code for interfacing with the original Fortran code from C++.

FortWrap is compatible with both g95
and gfortran and is intended to be used with Fortran code that takes an object oriented approach and makes use of Fortran derived types.

The resulting wrapper code provides a C++ interface that wraps the Fortran derived types with C++ proxy classes.

Behavior of plates

I seriously studied FEM basics when I had attended the finite element method course in IISc, Bangalore. After that everything went back stage as work assignments took over.

Now as I need to revisit this topic again, I am again starting from the basics.

In this line, a colleague of mine passed this excellent PowerPoint on behavior of plates as used in FEM. It’s really good if you want to brush off the dust on your fea knowledge.

Do check it out.