Skip to content

Commit

Permalink
chaneg simulation to compression boundary conditions and add that in …
Browse files Browse the repository at this point in the history 
…wall example
  • Loading branch information
aradermacher committed Feb 13, 2024
1 parent d3e1841 commit c01553a
Show file tree
Hide file tree
Showing 6 changed files with 194 additions and 112 deletions.
104 changes: 82 additions & 22 deletions amworkflow/simulation/experiment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,59 +276,93 @@ def create_displacement_boundary(
# define boundary conditions generator
bc_generator = BoundaryConditions(self.mesh, V)

if self.p["bc_setting"] == "fixed_y_bottom":
# Attention:
# normally geometries are defined as x/y and z in height due to AM production z is the direction perpendicular to the layers with is not usually the loading direction

if self.p["bc_setting"] == "compr_disp_y":
# loading displacement controlled in y direction at max_y surface, whereas y-z surface at min_y is fixed
if self.p["dim"] == 3:
# normally geometries are defined as x/y and z in height due to AM production here y is the height!!
# get mesh_points to define boundaries
mesh_points = self.mesh.geometry.x
min_y = mesh_points[:, 1].min()
max_y = mesh_points[:, 1].max()

# min_x = mesh_points[:, 0].min()
# max_x = mesh_points[:, 0].max()
# min_z = mesh_points[:, 2].min()
# max_z = mesh_points[:, 2].max()

print("fix at y=", min_y)
print("check points at y_min", mesh_points[mesh_points[:, 1] == min_y])
self.logger.info("fix at y= %s", min_y)
self.logger.debug("check points at y_min %s", mesh_points[mesh_points[:, 1] == min_y])
# dofs at bottom y_min fixed in x, y and z direction
bc_generator.add_dirichlet_bc(
np.array([0.0, 0.0, 0.0], dtype=ScalarType),
boundary=plane_at(min_y, 1),
method="geometrical",
entity_dim=self.mesh.topology.dim - 1, # surface
)
# dy = 0.001
# bc_generator.add_dirichlet_bc(
# np.array([0.0, 0.0, 0.0], dtype=ScalarType),
# boundary=within_range(
# [min_x, min_y - dy, min_z], [max_x, min_y + dy, max_z]
# ),
# method="geometrical",
# entity_dim=self.mesh.topology.dim - 2, # points
# )

# top displacement at top (==max_y) in y direction
# look for max y values in mesh
print("apply disp where", max_y)
print("check points at max_y", mesh_points[mesh_points[:, 1] == max_y])
self.logger.info("apply disp where y= %s", max_y)
self.logger.debug("check points at max_y: %s", mesh_points[mesh_points[:, 1] == max_y])
bc_generator.add_dirichlet_bc(
self.top_displacement,
boundary=plane_at(max_y, 1),
sub=1,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1,
)

# # maybe add some threshold to find bc points:
# min_x = mesh_points[:, 0].min()
# max_x = mesh_points[:, 0].max()
# min_z = mesh_points[:, 2].min()
# max_z = mesh_points[:, 2].max()
# dy = 0.001
# bc_generator.add_dirichlet_bc(
# np.array([0.0, 0.0, 0.0], dtype=ScalarType),
# boundary=within_range(
# [min_x, max_y - dy, min_z], [max_x, max_y + dy, max_z]
# [min_x, min_y - dy, min_z], [max_x, min_y + dy, max_z]
# ),
# method="geometrical",
# entity_dim=self.mesh.topology.dim - 2, # points
# )
# define displacement sensor position for this set-up
min_x = mesh_points[:, 0].min()
min_z = mesh_points[:, 2].min()
self.sensor_location_corner_top = [min_x, max_y, min_z] # corner point
self.sensor_location_middle_endge = [min_x, (max_y-min_y)/2., min_z] # point in the middle of the one edge


elif self.p["bc_setting"] == "compr_disp_x":
# loading displacement controlled in x direction at max_x surface, whereas x-z surface at min_x is fixed
if self.p["dim"] == 3:
# get mesh_points to define boundaries
mesh_points = self.mesh.geometry.x
min_x = mesh_points[:, 0].min()
max_x = mesh_points[:, 0].max()

self.logger.info("fix at x=%s", min_x)
self.logger.debug("check points at x_min %s", mesh_points[mesh_points[:, 0] == min_x])
# dofs at bottom x_min fixed in x, y and z direction
bc_generator.add_dirichlet_bc(
np.array([0.0, 0.0, 0.0], dtype=ScalarType),
boundary=plane_at(min_x, 0),
method="geometrical",
entity_dim=self.mesh.topology.dim - 1, # surface
)

# top displacement at top (==max_x) in x direction
self.logger.info("apply disp where x=%s", max_x)
self.logger.debug("check points at max_x: %s", mesh_points[mesh_points[:, 0] == max_x])
bc_generator.add_dirichlet_bc(
self.top_displacement,
boundary=plane_at(max_x, 0),
sub=0,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1,
)
# define displacement sensor position for this set-up
min_y = mesh_points[:, 1].min()
min_z = mesh_points[:, 2].min()
self.sensor_location_corner_top = [max_x, min_y, min_z] # corner point
self.sensor_location_middle_endge = [(max_x-min_x)/2, min_y, min_z] # point in the middle of the one edge
else:
raise ValueError(f"Wrong boundary setting: {self.p['bc_setting']}")

Expand All @@ -344,6 +378,32 @@ def apply_displ_load(self, top_displacement: pint.Quantity | float) -> None:
top_displacement.ito_base_units()
self.top_displacement.value = top_displacement.magnitude

def boundary_bottom(self):
"""specifies boundary: plane at bottom because different for different bc_setting
Returns: fct defining if dof is at boundary
"""
if self.p["bc_setting"] == "compr_disp_y":
if self.p["dim"] == 3:
# get mesh_points to define boundaries
mesh_points = self.mesh.geometry.x
min_y = mesh_points[:, 1].min()

return plane_at(min_y, "y")

elif self.p["bc_setting"] == "compr_disp_x":
if self.p["dim"] == 3:
# get mesh_points to define boundaries
mesh_points = self.mesh.geometry.x
min_x = mesh_points[:, 0].min()

return plane_at(min_x, "x")

else:
raise ValueError(f"Wrong boundary setting: {self.p['bc_setting']}")


# def create_body_force(self, v: ufl.argument.Argument) -> ufl.form.Form:
# """defines body force
#
Expand All @@ -363,4 +423,4 @@ def apply_displ_load(self, top_displacement: pint.Quantity | float) -> None:
# f = df.fem.Constant(self.mesh, ScalarType(force_vector))
# L = ufl.dot(f, v) * ufl.dx
#
# return L
# return L
59 changes: 50 additions & 9 deletions amworkflow/simulation/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,13 @@

import dolfinx as df
import numpy as np
import pandas as pd
from fenicsxconcrete.finite_element_problem import ConcreteAM, ConcreteThixElasticModel
from fenicsxconcrete.finite_element_problem.linear_elasticity import LinearElasticity
from fenicsxconcrete.sensor_definition.reaction_force_sensor import ReactionForceSensor
from fenicsxconcrete.sensor_definition.displacement_sensor import DisplacementSensor
from fenicsxconcrete.sensor_definition.strain_sensor import StrainSensor
from fenicsxconcrete.sensor_definition.stress_sensor import StressSensor
from fenicsxconcrete.util import QuadratureEvaluator, ureg

from amworkflow.simulation.experiment import ExperimentProcess, ExperimentStructure
Expand Down Expand Up @@ -97,7 +101,7 @@ def run(self, mesh_file: Path, out_xdmf: Path) -> None:
experiment = ExperimentStructure(self.pint_parameters)

if self.pint_parameters["material_type"].magnitude == "linear":
problem = LinearElasticity(
self.problem = LinearElasticity(
experiment,
self.pint_parameters,
pv_name=out_xdmf.stem,
Expand All @@ -107,15 +111,52 @@ def run(self, mesh_file: Path, out_xdmf: Path) -> None:
raise NotImplementedError(
"material type not yet implemented: ", self.material_type
)
problem.experiment.apply_displ_load(
self.pint_parameters["top_displacement"]

# set sensors
self.problem.add_sensor(ReactionForceSensor())
self.problem.add_sensor(StressSensor(where=experiment.sensor_location_middle_endge, name="stress sensor"))
self.problem.add_sensor(StrainSensor(where=experiment.sensor_location_middle_endge, name="strain sensor"))
self.problem.add_sensor(DisplacementSensor(where=experiment.sensor_location_corner_top, name="disp top"))



for i in range(0,self.pint_parameters["number_steps"].magnitude):
d_disp = 0 + (i+1) * self.pint_parameters["top_displacement"] / self.pint_parameters["number_steps"].magnitude
print('d_disp', d_disp)
self.problem.experiment.apply_displ_load(d_disp)
self.problem.solve()
self.problem.pv_plot()
print(
f"computed disp step: {i}, u_max={self.problem.fields.displacement.x.array[:].max()}, u_min={self.problem.fields.displacement.x.array[:].min()}"
)
print(
"force_x_y_z",
np.array(self.problem.sensors["ReactionForceSensor"].data),self.problem.sensors["ReactionForceSensor"].units
)
print("disp sensor", np.array(self.problem.sensors["disp top"].data))
print("stress sensor", np.array(self.problem.sensors["stress sensor"].data))
print("strain sensor", np.array(self.problem.sensors["strain sensor"].data))

if self.pint_parameters["bc_setting"].magnitude == "compr_disp_y":
force_y = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 1].max()
u_max = d_disp
print('equivalent modulus', force_y/(u_max*0.15))
elif self.pint_parameters["bc_setting"].magnitude == "compr_disp_x":
force_y = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 0].max()
u_max = d_disp
print('equivalent modulus', force_y/(u_max*0.15))

# save sensor output as csv file for postprocessing
sensor_dict = {}
sensor_dict["disp_x"] = np.array(self.problem.sensors["disp top"].data)[:, 0]
sensor_dict["disp_y"] = np.array(self.problem.sensors["disp top"].data)[:, 1]
sensor_dict["disp_z"] = np.array(self.problem.sensors["disp top"].data)[:, 2]
sensor_dict["force_x"] = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 0]
sensor_dict["force_y"] = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 1]
sensor_dict["force_z"] = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 2]
df_sensor = pd.DataFrame(sensor_dict)
df_sensor.to_csv(out_xdmf.parent / f"{out_xdmf.stem}.csv", index=False)

problem.solve()
problem.pv_plot()
print("problem solved")
print("max u", problem.fields.displacement.vector.array[:].max())
print("min u", problem.fields.displacement.vector.array[:].min())

elif self.pint_parameters["experiment_type"].magnitude == "process":
experiment = ExperimentProcess(self.pint_parameters)
Expand Down Expand Up @@ -247,4 +288,4 @@ def initial_path(self, time_layer, t_0: float = 0):
# set new path
self.problem.set_initial_path(q_path)

return
return
Loading

0 comments on commit c01553a

Please sign in to comment.