Skip to content
This repository has been archived by the owner on Feb 13, 2024. It is now read-only.

Commit

Permalink
v0.91
Browse files Browse the repository at this point in the history
Added readme.txt
Example files renamed
Added overview pdf
  • Loading branch information
spacether committed Dec 23, 2014
1 parent deb74d6 commit bd1f004
Show file tree
Hide file tree
Showing 10 changed files with 235 additions and 46 deletions.
Binary file removed Pycalculix_2014-12-11_Boston.pdf
Binary file not shown.
Binary file removed Pycalculix_2014-12-11_Boston.pptx
Binary file not shown.
Binary file added Pycalculix_Overview_2014-12-22.pdf
Binary file not shown.
108 changes: 108 additions & 0 deletions README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@
Description:
pycalculix is a Python 3 library to automate and build finite element analysis (FEA) models in Calculix.
Meshing uses Calculix or GMSH.
Website: http://justinablack.com/pycalculix/
Source Code: https://github.com/spacether/pycalculix

Usefull applications of Pycalculix:
-Trade studies for plane stress, splane strain, or axisymmetric parts
-Quick Kt analysis of 2D geometry
-Leaning finite element analyis (FEA) and Python

Notes:
I build a chunker in python which tries to cut big areas (> 5 sides) which
cgx can't mesh into smaller areas (<= 5 sides) which are meshable in cgx.
The chunker may not always work. I'd like to integrate gmsh in the future
so mesh quality can be more controllable, and less error prone.

License:
See LICENSE.txt (GPL v2)

Creator:
Justin Black, justin.a.black[at-sine]gmail[dot]com
Initial Release: December 2014

Elements Supported:
Axisymmetric, plane stress, and plane strain elements are supported.
First and second order triangles and quadrilaterals are supported.
First order elments only have corner nodes
Second order elements have midside nodes
Second order elements produce more accurate results
Setting element divisions on lines is supported

Geometry Building:
One can build separate parts made of points, lines, arcs, and areas.
Straight lines and arcs are currently supported.
One can draw a part made of straight lines, then smooth out corners by adding
blends/fillets with the part method: part.fillet_lines(L1, L2, arc_radius)
The new filleet will be tangent to both adjacent lines.

Loading:
Force loading, constant pressure, linearly varying pressure, gravity, and rotational speed forces are implemented.
Displacement constraints are also supported.
Loads are stored on geometry primitives (points lines etc) and can be applied before or after meshing.

Getting Started:
Please read through the documented example files below:
example1_dam.py
example2_hole_in_plate.py
example3_compr_rotor.py
example4_hole_kt.py
example5_times_dam.py

Files Produced:
Meshing and solving are done in the background using cgx or gmsh for meshing, and Calculix ccx for solving.
Files Used:
*.fbd (Calculix cgx gemetry file)
*.inp (Calculix solver input file, or mesh definition)
*.geo (Gmsh geometry file)
*.msh (Gmsh native mesh file)
*.frd (Calculix ccx nodal results file, values are at nodes and were created by interpolating element integration point results back to the nodes)
*.dat (Calculix ccx element results file, includes integration point results)

Installation:
Install the below required software.
Then move pycalculix.py to the directory you want to work in.
Any new models should be created in separate .py files. See the above example files in the 'Getting Started' section.

Required Software:
Calculix:
Free and very full-featured preprocessor (cgx) and finite element analysis (FEA) solver (ccx)
http://www.calculix.de/
Linux Version: http://www.dhondt.de/
Windows Version: http://www.bconverged.com/download.php#calculix
Gmsh:
Free and full featured mesher
http://geuz.org/gmsh/#Download

Anaconda:
An installation package that includes Python and many often used python libraries
Anaconda includes the below Python3+, Numpy, and Matplotlib.
If you are a Python beginner, I suggest downloading and installing it rather than the below separate installers.
Url:
http://continuum.io/downloads#py34

Python 3+:
https://www.python.org/downloads/release/python-342/
Numpy:
http://sourceforge.net/projects/numpy/files/NumPy/1.9.1/
Matplotlib:
http://matplotlib.org/downloads.html

Optional Software:
Suggested IDE (program to edit and run python programs):
Wing IDE:
http://wingware.com/downloads/wingide-101

Future Goals:
-Plotting
Add element results plotting
-Add compression supports
-Make lists for lines and signed lines (need to write a signed lines class)
-Add struct-thermal and thermal support
-Auto-detect contact regions between parts
-CAD import of brep and igs via gmsh
-CAD export via gmsh (step, brep)
-Add mp4 maker
-Bolted joint example perhaps, nodal thickness on bolt and nut areas
-Interactive selection?
94 changes: 94 additions & 0 deletions example1_dam.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
from pycalculix import FeaModel
import math

# We'll be modeling a masonry gravity dam, the Beetaloo dam in Australia
# Problem constants
proj_name = 'example1_dam'
grav = 9.81 # m/s^2
dens_water = 1000 #kg/m^3
press_atm = 101325 # Pascals = N/m^2 = kg/(m-s^2)
dam_ht_ft = 115
water_ht_ft = 109
length_dam_ft = 580

# derived dims
water_ax = (water_ht_ft-22)/math.tan(math.radians(89)) + 12
top_ax = (dam_ht_ft-22)/math.tan(math.radians(89)) + 12
thickness = length_dam_ft*0.3048 # conversion to metric

# derived dimensions, tan = o/a --> o/tan = a
pts_ft_water = [[8,0],[22,12],[water_ht_ft,water_ax]]
pts_ft_air = [[dam_ht_ft,top_ax], [dam_ht_ft, top_ax + 14]]
water_ht_m = water_ht_ft*0.3048

# make model
a = FeaModel(proj_name)
a.set_units('m') # this sets dist units to meters, labels our consistent units

# make part, coordinates are x, y = radial, axial
b = a.PartMaker()
b.goto(0.0,0.0)
water_lines = []
air_lines = []
for [x,y] in pts_ft_water:
[x,y] = [x*0.3048,y*0.3048] # conversion to metric
[L1,p1,p2] = b.draw_line_to(x, y)
water_lines.append(L1)
for [x,y] in pts_ft_air:
[x,y] = [x*0.3048,y*0.3048] # conversion to metric
[L1,p1,p2] = b.draw_line_to(x, y)
air_lines.append(L1)
# make the two arcs
pts_ft_arcs = [ [[22,73],[146,208]], [[14,98],[41,93]] ]
for [[x,y],[xc,yc]] in pts_ft_arcs:
[x,y] = [x*0.3048,y*0.3048]
[xc,yc] = [xc*0.3048,yc*0.3048]
[L1,p1,p2] = b.draw_arc(x,y,xc,yc)
air_lines.append(L1)
# make the last 3 lines after the arcs
pts_ft_other = [ [2,110], [0, 110] ]
for [x,y] in pts_ft_other:
[x,y] = [x*0.3048,y*0.3048] # conversion to metric
[L1,p1,p2] = b.draw_line_to(x, y)
air_lines.append(L1)
b.draw_line_to(0, 0)
a.plot_geometry(proj_name+'_geom', b) # view the points, lines, and areas

# set part material
mat = a.MatlMaker('concrete')
mat.set_mech_props(2300, 30000*(10**6), 0.2)
a.set_matl(mat, b)

# set the element type, line division, and mesh the database
a.set_eshape('quad', 2)
a.set_etype(b, 'plstrain', thickness)
b.get_item('L8').set_ediv(2)
a.mesh(0.5, 'gmsh') # mesh with 1.0 or less fineness, smaller is finer
a.plot_elements(proj_name+'_elem') # plot the part elements

# set loads and constraints
a.set_load('press', air_lines, press_atm)
a.set_fluid_press(water_lines, dens_water, grav, water_ht_m, press_atm)
a.set_gravity(grav, b)
a.set_constr('fix', b.bottom, 'x')
a.set_constr('fix', b.bottom, 'y')
a.plot_pressures(proj_name+'_press_1')
a.plot_constraints(proj_name+'_constr')

# make model and solve it
mod = a.ModelMaker(b, 'struct')
mod.solve()

# query results and store them
disp = False # turn off display plotting
fields = 'Seqv,Sx,Sy,Sz,S1,S2,S3,ux,uy,utot' # store the fields to write
fields = fields.split(',')

for field in fields:
fname = proj_name+'_'+field
mod.rfile.nplot(field, fname, display=disp)
smax = mod.rfile.get_nmax('Seqv')
[fx, fy, fz] = mod.rfile.get_fsum(a.get_item('L9'))
print('Seqv_max= %3.2f' % (smax))
print('Reaction forces (fx,fy,fz) = (%12.10f, %12.10f, %12.10f)' % (fx, fy, fz))

18 changes: 7 additions & 11 deletions hole_model.py → example2_hole_in_plate.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
from pycalculix import FeaModel

# Vertical hole in plate model, make model
proj_name = 'hole_model'
proj_name = 'example2_hole_in_plate'
a = FeaModel(proj_name)
a.set_units('m') # this sets dist units to meters, labels our consistent units

Expand Down Expand Up @@ -38,7 +38,7 @@

# set part material
mat = a.MatlMaker('steel')
mat.set_mech_props(7800, 210000, 0.3)
mat.set_mech_props(7800, 210*(10**9), 0.3)
a.set_matl(mat, b)

# set the element type and mesh database
Expand All @@ -62,13 +62,9 @@

# Plot results
disp = False
mod.rfile.nplot('Sx', proj_name+'_Sx', display=disp)
mod.rfile.nplot('Sy', proj_name+'_Sy', display=disp)
mod.rfile.nplot('S1', proj_name+'_S1', display=disp)
mod.rfile.nplot('S2', proj_name+'_S2', display=disp)
mod.rfile.nplot('S3', proj_name+'_S3', display=disp)
mod.rfile.nplot('Seqv', proj_name+'_Seqv', display=disp)
mod.rfile.nplot('ux', proj_name+'_ux', display=disp)
mod.rfile.nplot('uy', proj_name+'_uy', display=disp)
mod.rfile.nplot('utot', proj_name+'_utot', display=disp)
fields = 'Sx,Sy,S1,S2,S3,Seqv,ux,uy,utot' # store the fields to plot
fields = fields.split(',')
for field in fields:
fname = proj_name+'_'+field
mod.rfile.nplot(field, fname, display=disp)

19 changes: 6 additions & 13 deletions compr_rotor_stage.py → example3_compr_rotor.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,8 @@
import math

# We'll be modeling a rotating jet engine part

# Problem constants
proj_name = 'compr_rotor_stage'
proj_name = 'example3_compr_rotor'

# problem + geometry constants
rpm = 1000 # rotor speed in rpm
Expand Down Expand Up @@ -94,14 +93,8 @@
# view the results
disp = False # turn off display plotting
g = 0.0 # turn off adding results displacement to the model
mod.rfile.nplot('Sx', proj_name+'_Sx', display=disp, gmult=g)
mod.rfile.nplot('Sy', proj_name+'_Sy', display=disp, gmult=g)
mod.rfile.nplot('Sz', proj_name+'_Sz', display=disp, gmult=g)
mod.rfile.nplot('S1', proj_name+'_S1', display=disp, gmult=g)
mod.rfile.nplot('S2', proj_name+'_S2', display=disp, gmult=g)
mod.rfile.nplot('S3', proj_name+'_S3', display=disp, gmult=g)
mod.rfile.nplot('Seqv', proj_name+'_Seqv', display=disp, gmult=g)
mod.rfile.nplot('ux', proj_name+'_ux', display=disp, gmult=g)
mod.rfile.nplot('uy', proj_name+'_uy', display=disp, gmult=g)
mod.rfile.nplot('uz', proj_name+'_uz', display=disp, gmult=g)
mod.rfile.nplot('utot', proj_name+'_utot', display=disp, gmult=g)
fields = 'Sx,Sy,Sz,S1,S2,S3,Seqv,ux,uy,uz,utot' # store the fields to plot
fields = fields.split(',')
for field in fields:
fname = proj_name+'_'+field
mod.rfile.nplot(field, fname, display=disp, gmult=g)
13 changes: 7 additions & 6 deletions hole_kt.py → example4_hole_kt.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from pycalculix import FeaModel, frange
from pycalculix import FeaModel
import matplotlib.pyplot as plt
import math
import numpy as np
Expand Down Expand Up @@ -34,7 +34,8 @@ def kt_peterson(ratio):
left = right - rad

# vertical hole in plate model, make model
a = FeaModel('hole_kt')
proj_name = 'example4_hole_kt'
a = FeaModel(proj_name)
a.set_units('m') # this sets dist units to meters, labels our consistent units

# make part, coordinates are x, y = radial, axial
Expand All @@ -48,7 +49,7 @@ def kt_peterson(ratio):
b.draw_line_ax(-bot)
# b.plot_geometry('hole_kt_prechunk', display=disp) # view the geometry, points, lines, and areas
b.chunk()
a.plot_geometry('hole_kt_chunked', display=disp) # view the geometry, points, lines, and areas
a.plot_geometry(proj_name+'_chunked', display=disp) # view the geometry, points, lines, and areas

# set loads and constraints
a.set_load('press',b.top,-1*stress_val)
Expand All @@ -66,9 +67,9 @@ def kt_peterson(ratio):

a.set_eshape('tri', 2)
a.set_etype(b, 'plstress', thickness)
a.mesh(1, 'gmsh') # mesh with 1.0 fineness, smaller is finer
a.plot_elements('hole_kt_elem_%.3f' % (ratio), display=disp) # plot the part elements
a.plot_pressures('hole_kt_press', display=disp)
a.mesh(1.0, 'gmsh') # mesh with 1.0 fineness, smaller is finer
a.plot_elements('%s_elem_%.3f' % (proj_name, ratio), display=disp)
a.plot_pressures('%s_press' % (proj_name), display=disp)

# make model and solve it
mod = a.ModelMaker(b, 'struct')
Expand Down
27 changes: 12 additions & 15 deletions dam.py → example5_times_dam.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@
from pycalculix import FeaModel, frange
import matplotlib.pyplot as plt
from pycalculix import FeaModel
import math

# We'll be modeling a masonry gravity dam, the Beetaloo dam in Australia

# This time, we'll include multiple time steps
# Problem constants
proj_name = 'dam'
proj_name = 'example5_dam'
grav = 9.81 # m/s^2
dens_water = 1000 #kg/m^3
press_atm = 101325 # Pascals = N/m^2 = kg/(m-s^2)
Expand Down Expand Up @@ -54,7 +53,7 @@
[L1,p1,p2] = b.draw_line_to(x, y)
air_lines.append(L1)
b.draw_line_to(0, 0)
a.plot_geometry(proj_name+'_geom', b) # view the geometry, points, lines, and areas
a.plot_geometry(proj_name+'_geom', b) # view the points, lines, and areas

# set part material
mat = a.MatlMaker('concrete')
Expand All @@ -65,7 +64,7 @@
a.set_eshape('quad', 2)
a.set_etype(b, 'plstrain', thickness)
b.get_item('L8').set_ediv(2)
a.mesh(0.5, 'gmsh') # mesh with 1.0 fineness, smaller is finer
a.mesh(0.5, 'gmsh') # mesh with 1.0 or less fineness, smaller is finer
a.plot_elements(proj_name+'_elem') # plot the part elements

# set loads and constraints
Expand All @@ -77,13 +76,13 @@
a.plot_pressures(proj_name+'_press_1')
a.plot_constraints(proj_name+'_constr')

# Time = 2, ambient pressure and gravity
a.set_time(2.0)
# ambient pressure and gravity
a.set_load('press', water_lines+air_lines, press_atm)
a.plot_pressures(proj_name+'_press_2')

# Time = 3, gravity only
a.set_time(3.0)
# gravity only
a.set_load('press', water_lines+air_lines, 0.0)
a.plot_pressures(proj_name+'_press_3')

Expand All @@ -101,10 +100,8 @@
for field in fields:
fname = '%s_%i_%s' % (proj_name, int(time), field)
mod.rfile.nplot(field, fname, display=disp)

'''
smax = mod.rfile.get_nmax('Seqv')
[fx, fy, fz] = mod.rfile.get_fsum(a.get_item('L9'))
print('Seqv_max= %3.2f' % (smax))
print('Reaction forces (fx,fy,fz) = (%12.10f, %12.10f, %12.10f)' % (fx, fy, fz))
'''
smax = mod.rfile.get_nmax('Seqv')
[fx, fy, fz] = mod.rfile.get_fsum(a.get_item('L9'))
print('Seqv_max= %3.2f' % (smax))
print('Reaction forces (fx,fy,fz) = (%12.10f, %12.10f, %12.10f)' % (fx, fy, fz))

2 changes: 1 addition & 1 deletion pycalculix.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/python
"""
calculix.py
pycalculix.py
Description:
This module provides a python interface to the calculix FEA preprocessor and
solver. Simplicity is emphasized.
Expand Down

0 comments on commit bd1f004

Please sign in to comment.