#!/usr/bin/env python3
import math
import sys

import pycalculix as pyc

# We'll be modeling a masonry gravity dam, the Beetaloo dam in Australia
# This time, we'll include multiple time steps
# make model
model_name = 'dam-times'
model = pyc.FeaModel(model_name)
# this sets dist units to meters, labels our consistent units
model.set_units('m')

# set whether or not to show gui plots
show_gui = True
if '-nogui' in sys.argv:
    show_gui = False
# set element shape
eshape = 'quad'
if '-tri' in sys.argv:
    eshape = 'tri'

# Problem constants
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 part, coordinates are x, y = radial, axial
part = pyc.Part(model)
part.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] = part.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] = part.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] = part.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] = part.draw_line_to(x, y)
    air_lines.append(L1)
part.draw_line_to(0, 0)
# view the points, lines, and areas
model.plot_geometry(model_name+'_geom', display=show_gui)

# set part material
mat = pyc.Material('concrete')
mat.set_mech_props(2300, 30000*(10**6), 0.2)
model.set_matl(mat, part)

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

# set loads and constraints, all loads at first are at Time = 0s
model.set_load('press', air_lines, press_atm)
model.set_fluid_press(water_lines, dens_water, grav, water_ht_m,
                      press_atm)
model.set_gravity(grav, part)
model.set_constr('fix', part.bottom, 'x')
model.set_constr('fix', part.bottom, 'y')
model.plot_pressures(model_name+'_press_1', display=show_gui)
model.plot_constraints(model_name+'_constr', display=show_gui)

# Time = 2s, ambient pressure and gravity
model.set_time(2.0)
model.set_load('press', water_lines+air_lines, press_atm)
model.plot_pressures(model_name+'_press_2', display=show_gui)

# Time = 3s, gravity only
model.set_time(3.0)
model.set_load('press', water_lines+air_lines, 0.0)
model.plot_pressures(model_name+'_press_3', display=show_gui)

# make model and solve it
prob = pyc.Problem(model, 'struct')
prob.solve()

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

# store max and min values so plots can have the same max and min
# over all times
max_val, min_val = {}, {}
for field in fields:
    max_vals = [prob.rfile.get_emax(field, time) for time in prob.rfile.steps]
    min_vals = [prob.rfile.get_emin(field, time) for time in prob.rfile.steps]
    max_val[field] = max(max_vals)
    min_val[field] = min(min_vals)

#store titles for each time point
titles= {1.0: 'gravity + air pressure + water pressure',
         2.0: 'gravity + air pressure',
         3.0: 'gravity only'}
# plot results
for time in prob.rfile.steps:
    prob.rfile.set_time(time)
    for field in fields:
        fname = '%s_%s_%i' % (model_name, field, int(time))
        vmax, vmin = max_val[field], min_val[field]
        prob.rfile.eplot(field, fname, display=disp,
                         max_val=vmax, min_val=vmin, title=titles[time])

b5fd8a5a56660cd8e1df2a8f4d28ab7e5535f8e3
<mask token>

<mask token>
model.set_units('m')
<mask token>
if '-nogui' in sys.argv:
    show_gui = False
<mask token>
if '-tri' in sys.argv:
    eshape = 'tri'
<mask token>
part.goto(0.0, 0.0)
<mask token>
for [x, y] in pts_ft_water:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    water_lines.append(L1)
for [x, y] in pts_ft_air:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    air_lines.append(L1)
<mask token>
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] = part.draw_arc(x, y, xc, yc)
    air_lines.append(L1)
<mask token>
for [x, y] in pts_ft_other:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    air_lines.append(L1)
part.draw_line_to(0, 0)
model.plot_geometry(model_name + '_geom', display=show_gui)
<mask token>
mat.set_mech_props(2300, 30000 * 10 ** 6, 0.2)
model.set_matl(mat, part)
model.set_eshape('quad', 2)
model.set_etype('plstrain', part, thickness)
model.set_ediv('L8', 2)
model.mesh(0.5, 'gmsh')
model.plot_elements(model_name + '_elem', display=show_gui)
model.set_load('press', air_lines, press_atm)
model.set_fluid_press(water_lines, dens_water, grav, water_ht_m, press_atm)
model.set_gravity(grav, part)
model.set_constr('fix', part.bottom, 'x')
model.set_constr('fix', part.bottom, 'y')
model.plot_pressures(model_name + '_press_1', display=show_gui)
model.plot_constraints(model_name + '_constr', display=show_gui)
model.set_time(2.0)
model.set_load('press', water_lines + air_lines, press_atm)
model.plot_pressures(model_name + '_press_2', display=show_gui)
model.set_time(3.0)
model.set_load('press', water_lines + air_lines, 0.0)
model.plot_pressures(model_name + '_press_3', display=show_gui)
<mask token>
prob.solve()
<mask token>
for field in fields:
    max_vals = [prob.rfile.get_emax(field, time) for time in prob.rfile.steps]
    min_vals = [prob.rfile.get_emin(field, time) for time in prob.rfile.steps]
    max_val[field] = max(max_vals)
    min_val[field] = min(min_vals)
<mask token>
for time in prob.rfile.steps:
    prob.rfile.set_time(time)
    for field in fields:
        fname = '%s_%s_%i' % (model_name, field, int(time))
        vmax, vmin = max_val[field], min_val[field]
        prob.rfile.eplot(field, fname, display=disp, max_val=vmax, min_val=
            vmin, title=titles[time])

<mask token>
model_name = 'dam-times'
model = pyc.FeaModel(model_name)
model.set_units('m')
show_gui = True
if '-nogui' in sys.argv:
    show_gui = False
eshape = 'quad'
if '-tri' in sys.argv:
    eshape = 'tri'
grav = 9.81
dens_water = 1000
press_atm = 101325
dam_ht_ft = 115
water_ht_ft = 109
length_dam_ft = 580
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
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
part = pyc.Part(model)
part.goto(0.0, 0.0)
water_lines = []
air_lines = []
for [x, y] in pts_ft_water:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    water_lines.append(L1)
for [x, y] in pts_ft_air:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    air_lines.append(L1)
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] = part.draw_arc(x, y, xc, yc)
    air_lines.append(L1)
pts_ft_other = [[2, 110], [0, 110]]
for [x, y] in pts_ft_other:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    air_lines.append(L1)
part.draw_line_to(0, 0)
model.plot_geometry(model_name + '_geom', display=show_gui)
mat = pyc.Material('concrete')
mat.set_mech_props(2300, 30000 * 10 ** 6, 0.2)
model.set_matl(mat, part)
model.set_eshape('quad', 2)
model.set_etype('plstrain', part, thickness)
model.set_ediv('L8', 2)
model.mesh(0.5, 'gmsh')
model.plot_elements(model_name + '_elem', display=show_gui)
model.set_load('press', air_lines, press_atm)
model.set_fluid_press(water_lines, dens_water, grav, water_ht_m, press_atm)
model.set_gravity(grav, part)
model.set_constr('fix', part.bottom, 'x')
model.set_constr('fix', part.bottom, 'y')
model.plot_pressures(model_name + '_press_1', display=show_gui)
model.plot_constraints(model_name + '_constr', display=show_gui)
model.set_time(2.0)
model.set_load('press', water_lines + air_lines, press_atm)
model.plot_pressures(model_name + '_press_2', display=show_gui)
model.set_time(3.0)
model.set_load('press', water_lines + air_lines, 0.0)
model.plot_pressures(model_name + '_press_3', display=show_gui)
prob = pyc.Problem(model, 'struct')
prob.solve()
disp = False
fields = 'S1,S2,S3,Seqv,Sx,Sy'
fields = fields.split(',')
max_val, min_val = {}, {}
for field in fields:
    max_vals = [prob.rfile.get_emax(field, time) for time in prob.rfile.steps]
    min_vals = [prob.rfile.get_emin(field, time) for time in prob.rfile.steps]
    max_val[field] = max(max_vals)
    min_val[field] = min(min_vals)
titles = {(1.0): 'gravity + air pressure + water pressure', (2.0):
    'gravity + air pressure', (3.0): 'gravity only'}
for time in prob.rfile.steps:
    prob.rfile.set_time(time)
    for field in fields:
        fname = '%s_%s_%i' % (model_name, field, int(time))
        vmax, vmin = max_val[field], min_val[field]
        prob.rfile.eplot(field, fname, display=disp, max_val=vmax, min_val=
            vmin, title=titles[time])

import math
import sys
import pycalculix as pyc
model_name = 'dam-times'
model = pyc.FeaModel(model_name)
model.set_units('m')
show_gui = True
if '-nogui' in sys.argv:
    show_gui = False
eshape = 'quad'
if '-tri' in sys.argv:
    eshape = 'tri'
grav = 9.81
dens_water = 1000
press_atm = 101325
dam_ht_ft = 115
water_ht_ft = 109
length_dam_ft = 580
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
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
part = pyc.Part(model)
part.goto(0.0, 0.0)
water_lines = []
air_lines = []
for [x, y] in pts_ft_water:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    water_lines.append(L1)
for [x, y] in pts_ft_air:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    air_lines.append(L1)
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] = part.draw_arc(x, y, xc, yc)
    air_lines.append(L1)
pts_ft_other = [[2, 110], [0, 110]]
for [x, y] in pts_ft_other:
    [x, y] = [x * 0.3048, y * 0.3048]
    [L1, p1, p2] = part.draw_line_to(x, y)
    air_lines.append(L1)
part.draw_line_to(0, 0)
model.plot_geometry(model_name + '_geom', display=show_gui)
mat = pyc.Material('concrete')
mat.set_mech_props(2300, 30000 * 10 ** 6, 0.2)
model.set_matl(mat, part)
model.set_eshape('quad', 2)
model.set_etype('plstrain', part, thickness)
model.set_ediv('L8', 2)
model.mesh(0.5, 'gmsh')
model.plot_elements(model_name + '_elem', display=show_gui)
model.set_load('press', air_lines, press_atm)
model.set_fluid_press(water_lines, dens_water, grav, water_ht_m, press_atm)
model.set_gravity(grav, part)
model.set_constr('fix', part.bottom, 'x')
model.set_constr('fix', part.bottom, 'y')
model.plot_pressures(model_name + '_press_1', display=show_gui)
model.plot_constraints(model_name + '_constr', display=show_gui)
model.set_time(2.0)
model.set_load('press', water_lines + air_lines, press_atm)
model.plot_pressures(model_name + '_press_2', display=show_gui)
model.set_time(3.0)
model.set_load('press', water_lines + air_lines, 0.0)
model.plot_pressures(model_name + '_press_3', display=show_gui)
prob = pyc.Problem(model, 'struct')
prob.solve()
disp = False
fields = 'S1,S2,S3,Seqv,Sx,Sy'
fields = fields.split(',')
max_val, min_val = {}, {}
for field in fields:
    max_vals = [prob.rfile.get_emax(field, time) for time in prob.rfile.steps]
    min_vals = [prob.rfile.get_emin(field, time) for time in prob.rfile.steps]
    max_val[field] = max(max_vals)
    min_val[field] = min(min_vals)
titles = {(1.0): 'gravity + air pressure + water pressure', (2.0):
    'gravity + air pressure', (3.0): 'gravity only'}
for time in prob.rfile.steps:
    prob.rfile.set_time(time)
    for field in fields:
        fname = '%s_%s_%i' % (model_name, field, int(time))
        vmax, vmin = max_val[field], min_val[field]
        prob.rfile.eplot(field, fname, display=disp, max_val=vmax, min_val=
            vmin, title=titles[time])

#!/usr/bin/env python3
import math
import sys

import pycalculix as pyc

# We'll be modeling a masonry gravity dam, the Beetaloo dam in Australia
# This time, we'll include multiple time steps
# make model
model_name = 'dam-times'
model = pyc.FeaModel(model_name)
# this sets dist units to meters, labels our consistent units
model.set_units('m')

# set whether or not to show gui plots
show_gui = True
if '-nogui' in sys.argv:
    show_gui = False
# set element shape
eshape = 'quad'
if '-tri' in sys.argv:
    eshape = 'tri'

# Problem constants
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 part, coordinates are x, y = radial, axial
part = pyc.Part(model)
part.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] = part.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] = part.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] = part.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] = part.draw_line_to(x, y)
    air_lines.append(L1)
part.draw_line_to(0, 0)
# view the points, lines, and areas
model.plot_geometry(model_name+'_geom', display=show_gui)

# set part material
mat = pyc.Material('concrete')
mat.set_mech_props(2300, 30000*(10**6), 0.2)
model.set_matl(mat, part)

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

# set loads and constraints, all loads at first are at Time = 0s
model.set_load('press', air_lines, press_atm)
model.set_fluid_press(water_lines, dens_water, grav, water_ht_m,
                      press_atm)
model.set_gravity(grav, part)
model.set_constr('fix', part.bottom, 'x')
model.set_constr('fix', part.bottom, 'y')
model.plot_pressures(model_name+'_press_1', display=show_gui)
model.plot_constraints(model_name+'_constr', display=show_gui)

# Time = 2s, ambient pressure and gravity
model.set_time(2.0)
model.set_load('press', water_lines+air_lines, press_atm)
model.plot_pressures(model_name+'_press_2', display=show_gui)

# Time = 3s, gravity only
model.set_time(3.0)
model.set_load('press', water_lines+air_lines, 0.0)
model.plot_pressures(model_name+'_press_3', display=show_gui)

# make model and solve it
prob = pyc.Problem(model, 'struct')
prob.solve()

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

# store max and min values so plots can have the same max and min
# over all times
max_val, min_val = {}, {}
for field in fields:
    max_vals = [prob.rfile.get_emax(field, time) for time in prob.rfile.steps]
    min_vals = [prob.rfile.get_emin(field, time) for time in prob.rfile.steps]
    max_val[field] = max(max_vals)
    min_val[field] = min(min_vals)

#store titles for each time point
titles= {1.0: 'gravity + air pressure + water pressure',
         2.0: 'gravity + air pressure',
         3.0: 'gravity only'}
# plot results
for time in prob.rfile.steps:
    prob.rfile.set_time(time)
    for field in fields:
        fname = '%s_%s_%i' % (model_name, field, int(time))
        vmax, vmin = max_val[field], min_val[field]
        prob.rfile.eplot(field, fname, display=disp,
                         max_val=vmax, min_val=vmin, title=titles[time])
