A PSE for Automatic Matlab 3D Finite Element Code Generation and Simplified Grid Computing

1
A PSE for Automatic Matlab 3DFinite Element Code
Generation and Simplified Grid Computing

• Zhou Jun, and Yukio Umetani
• Shizuoka University of Japan

2
Outline

• Motivation
• Previous work
• Automatic Matlab 3D FEM code generation
• Execution of Matlab script via UNICORE
• Conclusions and future work

3
Motivation

• Interface to Matlab
• popular for matrix computation
• visualization facility
• verifying simulation results of PSILAB
• Reduce Matlab FEM programming effort
• coding time
• learning time
• programming errors
• Run Matlab script on the Grid via UNICORE

4
GridPSi overview
PSILAB
PSILAB MATLAB
Local Machine
Internet
Unicore Client
GridPSi Client
GridPSi Server
Unicore Gateway
Remote Server
Usite
NJS
GridPSi Server
TSI
GridPSi Server
GridPSi client connects GridPSi server in three
ways.
5
PDE wizard-based problem modeling
6
GridPSi client arcitecture
MATLAB Script Generator Self-designed Matlab
Functions
PSILAB Script Generator
7
Automatic generation of Matlab 3D FEM code

• MATLAB Script Generator
• Self-coded Matlab functions

PSILAB Script Generator
Matlab Script Generator
Generic Data structures
Self-coded Matlab functions
8
Example problem
To investigate the vertical movements at the
bridge area of a 3D piano soundboard of
simplified shape that is subjected to vertical
pressures caused by the vibrations of the
strings. (The soundboard is in anisotropic
material)
9
Mathematical model
Governing equations MÜ CÙ KU R (1)
External forces
(2)
Boundary conditions

• Fixed Dirichlet conditions on the four sides
• Generalized Neumann conditions on the top and
  bottom

Initial conditions
No displacements at the initial state
10
Finite element solution
Solid linear tetrahedron elements
Central difference method
11
Matlab FEM data structures
Nodal coordinate matrix Element connectivity
matrices Problem domain Boundary regions
stiffness, mass, and damping matrices Nodal
force vector Surface tractions (Neumann boundary
condition) Unknown field on boundaries (Dirichlet
condition)
12
Build nodal coordinate and element connectivity
matrices
function node, setS, setV gmsh2matlab(meshFil
e)

• node a matrix of nodal coordinates (n-by-3)
• setS a set of element connectivity matrices of
  boundary surfaces
• setV a set of element connectivity matrices of
  the volume(s)
• meshFile The input mesh file

13
Build stiffness matrix
function K buildStiffnessMatrix3D(node,
elemV, E)
(3)
E the consistent stress-strain compliance
matrix function E buildComplianceMatriceElasti
c3D(material properties)
nG, Wk integration point and
weight function W,Q quadrature(order,type,dime
nsion) Jk the determinant of the Jacobian
matrix Bk the bar matrix of the element
function N,dNdxilagrange_basis('H4',pt)
return the Lagrange interpolant basis

and its gradients
14
function K buildStiffnessMatrix3D(node,
elemV, E) n size(node,1) get the total
node number n. m size (ElemV,1) get the
total element number m. K sparse (3n,3n)
create the sparse global stiffness matrix
K. for e (1m) loops over the
elements. str ElemV(e,) get the four
nodes comprising the current element. the
following is to build the Bar array strBK,
which is used to locate the positions of
the four nodes in the global stiffness matrix.
strBk str(1) str(1)n str(1)2n
str(2) str(2)n str(2)2n ...
str(3) str(3)n str(3)2n str(4)
str(4)n str(4)2n W,Q
quadrature(1,'TRIANGULAR',3) A function call
returns the quadratrue weight W and Q

for q 1size(W,1) pt
Q(q,) wk W(q)
N,dNdxilagrange_basis('H4',pt) A function
call returns the lagrange interpolant basis.
Jacobian node(str,)'dNdxi get the
Jacobian matrix. dNdX dNdxi
inv(Jacobian) a 4-by-3 matrix of the shape
functions. construct the
stress-displacement matrix BK (a 6-by-12 matrix)
Bk dNdX(1,1) 0 0
dNdX(2,1) 0 0 dNdX(3,1) 0 0 dNdX(4,1) 0
0 0 dNdX(1,2) 0 0
dNdX(2,2) 0 0 dNdX(3,2) 0 0 dNdX(4,2) 0
0 0 dNdX(1,3) 0 0 dNdX(2,3)
0 0 dNdX(3,3) 0 0 dNdX(4,3)
dNdX(1,2) dNdX(1,1) 0 dNdX(2,2) dNdX(2,1) 0
dNdX(3,2) dNdX(3,1) 0 dNdX(4,2) dNdX(4,1) 0
0 dNdX(1,3) dNdX(1,2) 0 dNdX(2,3) dNdX(2,2) 0
dNdX(3,3) dNdX(3,2) 0 dNdX(4,3) dNdX(4,2)
dNdX(1,3) 0 dNdX(1,1) dNdX(2,3) 0 dNdX(2,1)
dNdX(3,3) 0 dNdX(3,1) dNdX(4,3) 0 dNdX(4,1)
K(strBk,strBk) K(strBk,strBk) Bk E Bk
wk det(Jacobian) assemble the global
matrix. end end
15
Build mass and damping matrices
function M buildMassMatrix3D(node, elemV
,rho)
(4)
(the integral over the element)
16
Enforce boundary conditions

• Dirichlet boundary conditions

function v applyEssentialBCs(v,node,elemS,i
d,dx,dy,dz)
u t1 u t 2?tv update displacement vector
(5)

• Neumann boundary conditions

function fcalculateSurfaceTraction(node,elem
S,id,px,py,pz)
R f body f surface traction1 f surface
traction2
17
Initial Conditions
function Uold,Unow initalizeDisVectors(Uold
,Unow,initU,initV)

• Uold the displacement vector at time t-1
• Unow the displacement vector at time t
• initU the initial value at time -1
• initV the initial value at time 0

18
Solution Scheme
while (CT lt500) t tdelt CT
CT1 Apply the Neumann conditions
on given surfaces. f1
caculateSurfaceTraction2(node,elemS30,'30',0,0,-sz
) f2 caculateSurfaceTraction2(node,el
emS8,'8',0,0,0) R f f1 f2
Unew (a0Ma1C)\(Feffective -
(K-a2M)Unow - (a0M-a1C)Uold)
Vnow a1 (Unew - Uold)
Apply the Dirichlet conditions on given surfaces.
Vnow applyEssentialBCs(Vn
ow,node,elemS17,'17',0,0,0) Vnow
applyEssentialBCs(Vnow,node,elemS21,'21',0,0,0)
Vnow applyEssentialBCs(Vnow,node,elemS2
5,'25',0,0,0) Vnow
applyEssentialBCs(Vnow,node,elemS29,'29',0,0,0)
Unew a3 Vnow Uold
Uold Unow Unow Unew Uzz
Unew(2n13n,1) get the displacements in Z
direction. end fid fopen('matlabResult
s.txt','w') fwrite(fid,Uzz) fclose(fid)
19
(No Transcript)
20
Executing Matlab script via UNICORE

• Use of Script Task plug-in interface of UNIOCRE
  client
• Preparation of a Job containing two tasks with
  dependency
• 1. script task
• command matlab -nosplash -nodesktop
  -nodisplay -r soundboard
• import files the generated script,
  predefined Matlab functions, the mesh file
• 2. Export task retrieve the simulation
  results back

21
Job transferring and monitoring
22
Visual Effects of the Results
Results solved by MATLAB
Results solved by PSILAB
23
Correctness
In the same iteration step, at the same
integration point, the values of the small
displacements on the top of the soundboard solved
by Matlab and PSILAB are in the same order
(10-10) the first three digitals of the
values are almost the same
24
Related Works
PDE Toolbox (MathWorks) Only supports 1D and
2D PDE-based problems Geodise (University of
Southampton) Uses MATLAB as a scripting
language engine Users need to program
MATLAB. FEMLAB (COMSOL Inc.) Does not
support Grid access GridPSi Supports
real-world 2D and 3D PDE-based problems
Automatic generation MATLAB script without
programming Grid accessible.
25
Conclusions

• Layered architecture and component-based design
  makes our system easily interface to Matlab.
• Automatic Matlab FEM code generation greatly
  reduces users work and programming errors.
• Simple Grid access pattern lowers barriers to
  enter for users who wish to run their physical
  simulations on the Grid.

26
Future Work

• The real-time visualization through Unicore/VISIT
  and AVS/Express
• The design of nonlinear PDE Wizards

27
Thank You !