Hi
I am very new to Fortran, I am trying to find if there is a simple way to do partial derivative, like in MATLAB or python for instance-. Is there a package that I can use ?
Thank you !
Welcome to the forum. Some threads on automatic differentiation are
Welcome @dka!
Could you be more specific what type of partial derivatives do you need? Is it related to finite differences, finite elements, or a perhaps just a Jacobian calculation?
Hi @ivanpribec
Yes it is related to finite elements. I need to calculate the Jacobian and B-matrix for a specific element.
For analytical jacobian or partial derivatives, I think you would require symbolic computing. But for the numerical ones - simple forward, backward or central difference would do.
For which element are you doing ? Are you working on solid/structural elements - beam, plate, shell …?
Hi @Ashok
it is for a 4 nodes isoparametric element
For FEM, the Jacobian and other derivatives depend on the shape functions you are using. The Jacobian to map derivatives between iso-parametric and physical coordinates can be computed analytically (they are just polynomials or rational functions such as NURBS). The resulting equations are the same in any programming language so I’m not sure why you think Fortran would be any different. Any FEM textbook written in the last 50 years will discuss how to compute the derivatives and there are a multitude of open source FEM codes in various programming languages you can use as an example.
There are as many FEM textbooks as there are grains of sand on a beach but here are three that I find useful.
Thomas Hughes,“The Finite Element Method - Linear, Static, and Dynamic Finite Element Analysis”, Dover Publications
Covers the basics well and is a lot cheaper than other FEM textbooks which can be over $200 for the most recent editions.
Smith and Griffiths," Programming the Finite Element Method", Wiley
There is companion Fortran code for this book
Gennady Nikishkov, “Programming Finite Elements in Java”, Springer
A useful book if you want to take an object oriented approach because the author defines the various classes you need to create a FEM code and Java is a relatively easy language to translate into Fortran.
For the open source code, a web search or searching github will show you a wealth of working codes.
One can also search “finite elements” at my Fortran-related books list, which has links to codes when available.
For a 4 node isoparametric element, the shape functions are:
N1 = (1-s)(1-t)/4
N2 = (1-s)(1+t)/4
N3 = (1+s)(1+t)/4
N4 = (1+s)(1-t)/4
The B matrix would be as shown in the attached image. Replace x and y in the figure with s and t (isoparametric coordinates)
The partial derivatives can be calculated simply by hand and for each element you have to pass the nodal coordinates - to evaluate the B matrix for that element.
@Ashok @Beliavsky @rwmsu Thank you very much! I am actually aware of the mathematical part, I was wondering if there is a tool that does this automatically like diff() in matlab.
MATLAB has two diff
routines:
-
diff(X)
to calculate differences and approximate derivatives -
diff(f,var)
to differentiate a symbolic expression or function.
I guess the second one is what you had in mind. Here’s a MATLAB script, which generates the Jacobian for the shape functions given by @Ashok:
syms s t real
N1 = (1-s)*(1-t)/4;
N2 = (1-s)*(1+t)/4;
N3 = (1+s)*(1+t)/4;
N4 = (1+s)*(1-t)/4;
J = jacobian([N1,N2,N3,N4],[s, t]);
fortran(J,'file','four_node_jac.inc');
The resulting Fortran code:
t2 = s/4.0D0
t3 = t/4.0D0
t4 = -t2
t5 = -t3
A0(1,1) = t3-1.0D0/4.0D0
A0(1,2) = t2-1.0D0/4.0D0
A0(2,1) = t5-1.0D0/4.0D0
A0(2,2) = t4+1.0D0/4.0D0
A0(3,1) = t3+1.0D0/4.0D0
A0(3,2) = t2+1.0D0/4.0D0
A0(4,1) = t5+1.0D0/4.0D0
A0(4,2) = t4-1.0D0/4.0D0
You’d probably want to edit this to match expected variables names, or change the precision.
SymPy has a more customizable Fortran code printer compared to MATLAB; see Fortran printing in the SymPy documentation. A previous post of mine gives a complete example: Code generation using SymPy.