Skip to content

Commit

Permalink
add process sim class
Browse files Browse the repository at this point in the history
  • Loading branch information
aradermacher committed Sep 12, 2023
1 parent c24222a commit cdc02e7
Showing 1 changed file with 64 additions and 60 deletions.
124 changes: 64 additions & 60 deletions usecases/TrussArc/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,13 +119,13 @@ def pd_fkt(path_time: np.ndarray | list) -> np.ndarray | list:
return l_active

# copied from fenicsconcrete test_am_layer maybe put into lib fenicsxconcrete extended to 3D
def define_path(prob, t_diff, t_0=0):
def define_path(prob, params, t_diff, t_0=0):
"""create path as layer wise at quadrature space
one layer by time
prob: problem
param: parameter dictionary
param: parameter dictionary in original measure (conform with mesh nodes not required to be base units!!)
t_diff: time difference between each layer
t_0: start time for all (0 if static computation)
(-end_time last layer if dynamic computation)
Expand All @@ -142,22 +142,26 @@ def define_path(prob, t_diff, t_0=0):
x = positions.evaluate()
dof_map = np.reshape(x.flatten(), [len(q_path), 3])

# get parameters in conform unit to mesh!!
num_layers = params["num_layers"].magnitude
layer_height = params["layer_height"].magnitude

# select layers only by layer height - z coordinate
h_CO = np.array(dof_map)[:, 2]
h_min = np.arange(0, prob.p["num_layers"] * prob.p["layer_height"], prob.p["layer_height"])
h_min = np.arange(0, num_layers * layer_height, layer_height)
h_max = np.arange(
prob.p["layer_height"],
(prob.p["num_layers"] + 1) * prob.p["layer_height"],
prob.p["layer_height"],
layer_height,
(num_layers + 1) * layer_height,
layer_height,
)
# print("h_CO", h_CO)
# print("h_min", h_min)
# print("h_max", h_max)
print("h_CO", h_CO)
print("h_min", h_min)
print("h_max", h_max)
new_path = np.zeros_like(q_path)
EPS = 1e-8
for i in range(0, len(h_min)):
layer_index = np.where((h_CO > h_min[i] - EPS) & (h_CO <= h_max[i] + EPS))
new_path[layer_index] = t_0 + (prob.p["num_layers"] - 1 - i) * t_diff
new_path[layer_index] = t_0 + (num_layers - 1 - i) * t_diff

q_path = new_path

Expand Down Expand Up @@ -382,61 +386,61 @@ def run_amprocess(param):
param["dt"] = time_layer / 2
param["load_time"] = 1 * param["dt"] # interval where load is applied linear over time

# problem = ConcreteAM(
# experiment,
# param,
# nonlinear_problem=ConcreteThixElasticModel,
# pv_name=f"sim_01_{mesh_path.stem}",
# pv_path=mesh_path.parent,
# )
problem = ConcreteAM(
experiment,
param,
nonlinear_problem=ConcreteThixElasticModel,
pv_name=f"sim_01_{mesh_path.stem}",
pv_path=mesh_path.parent,
)

# # initial path function describing layer activation
# path_activation = define_path(
# problem, time_layer.magnitude, t_0=-(param["num_layers"].magnitude - 1) * time_layer.magnitude
# )
# problem.set_initial_path(path_activation)
# initial path function describing layer activation
path_activation = define_path(
problem, param, time_layer.magnitude, t_0=-(param["num_layers"].magnitude - 1) * time_layer.magnitude
)
problem.set_initial_path(path_activation)

total_time = param["num_layers"] * time_layer
# while problem.time <= total_time.to_base_units().magnitude:
# problem.solve()
# problem.pv_plot()
# print(f"computed disp t={problem.time}, u_max={problem.fields.displacement.x.array[:].max()}")
#


# test layer activation
V = df.fem.FunctionSpace(experiment.mesh, ("DG", 0)) # one value per element
path = df.fem.Function(V, name="path time")
dens = df.fem.Function(V, name="pseudo density")
path_out = define_path_time(path, param, time_layer)
print('path out', path_out.x.array[:])
while problem.time <= total_time.to_base_units().magnitude:
problem.solve()
problem.pv_plot()
print(f"computed disp t={problem.time}, u_max={problem.fields.displacement.x.array[:].max()}")

# output file
with df.io.XDMFFile(experiment.mesh.comm, mesh_path.parent / f"sim_01_{mesh_path.stem}.xdmf", "w") as f:
f.write_mesh(experiment.mesh)

t = 0
dt = param["dt"].to_base_units().magnitude
while t <= total_time.to_base_units().magnitude:
print(f"time {t} compute path and density")
# compute density based on path function
print('path out', path_out.x.array[:].min(), path_out.x.array[:].max())
dens_list = pd_fkt(path_out.x.array[:])
dens.x.array[:] = dens_list[:]
print(dens_list.min(), dens_list.max())
print('dens', dens.x.array[:].min(), dens.x.array[:].max())


# export path & density into xdmf for t
dens.name = "pseudo density"
path_out.name = "path time"
with df.io.XDMFFile(experiment.mesh.comm, mesh_path.parent / f"sim_01_{mesh_path.stem}.xdmf", "a") as f:
f.write_function(dens, t)
f.write_function(path_out, t)

# update path and time
path_out.x.array[:] += dt
t += dt
# # test layer activation
# V = df.fem.FunctionSpace(experiment.mesh, ("DG", 0)) # one value per element
# path = df.fem.Function(V, name="path time")
# dens = df.fem.Function(V, name="pseudo density")
# path_out = define_path_time(path, param, time_layer)
# print('path out', path_out.x.array[:])
#
# # output file
# with df.io.XDMFFile(experiment.mesh.comm, mesh_path.parent / f"sim_01_{mesh_path.stem}.xdmf", "w") as f:
# f.write_mesh(experiment.mesh)
#
# t = 0
# dt = param["dt"].to_base_units().magnitude
# while t <= total_time.to_base_units().magnitude:
# print(f"time {t} compute path and density")
# # compute density based on path function
# print('path out', path_out.x.array[:].min(), path_out.x.array[:].max())
# dens_list = pd_fkt(path_out.x.array[:])
# dens.x.array[:] = dens_list[:]
# print(dens_list.min(), dens_list.max())
# print('dens', dens.x.array[:].min(), dens.x.array[:].max())
#
#
# # export path & density into xdmf for t
# dens.name = "pseudo density"
# path_out.name = "path time"
# with df.io.XDMFFile(experiment.mesh.comm, mesh_path.parent / f"sim_01_{mesh_path.stem}.xdmf", "a") as f:
# f.write_function(dens, t)
# f.write_function(path_out, t)
#
# # update path and time
# path_out.x.array[:] += dt
# t += dt



Expand Down Expand Up @@ -466,7 +470,7 @@ def run_amprocess(param):
'degree': 2 * ureg(""),
'top_displacement': -50.0*ureg("mm"),
'layer_height': 10*ureg("mm"), # for process simulation to activate layer by layer
'num_layers': 10*ureg(""), # 5 in truss example 10 in geo test
'num_layers': 5*ureg(""), # 5 in truss example 10 in geo test
}

# run_amstructure(params)
Expand Down

0 comments on commit cdc02e7

Please sign in to comment.