diff --git a/usecases/TrussArc/simulation.py b/usecases/TrussArc/simulation.py index 1ddd55e..73d034f 100644 --- a/usecases/TrussArc/simulation.py +++ b/usecases/TrussArc/simulation.py @@ -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) @@ -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 @@ -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 @@ -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)