Skip to content

Commit

Permalink
40 compression test simulation (#42)
Browse files Browse the repository at this point in the history 
* chaneg simulation to compression boundary conditions and add that in wall example

* fixed brocken test

* change some stuff

* fixed brocken test

* update

* add possibility of no infill in param wall

* change young#s modulus parameter

* saved .doit files in wall example

---------

Co-authored-by: aradermacher <[email protected]>
  • Loading branch information
@aradermacher
aradermacher and aradermacher authored Feb 20, 2024
1 parent 316189c commit d33d7cc
Show file tree
Hide file tree
Showing 7 changed files with 310 additions and 187 deletions.
15 changes: 15 additions & 0 deletions amworkflow/geometry/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,21 @@ def geometry_spawn(self) -> TopoDS_Shape:
)
points = [[0,0,0],[0,0,0]]

elif self.infill == "no":
W = self.width
L = self.length
t= self.line_width
points = [[0., 0., 0.],
[0., (W - t) / 2, 0.],
[(L - t), (W - t) / 2, 0.],
[(L - t), -(W - t) / 2, 0.],
[0., -(W - t) / 2, 0.]]

creator = composite_geometries.CreateWallByPoints(
points, th=self.line_width, height=self.height
)
shape = creator.Shape()

elif self.infill == "honeycomb":
points = self.honeycomb_infill(self.length, self.width, self.line_width)
creator = composite_geometries.CreateWallByPoints(
Expand Down
200 changes: 163 additions & 37 deletions amworkflow/simulation/experiment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from fenicsxconcrete.boundary_conditions.bcs import BoundaryConditions
from fenicsxconcrete.boundary_conditions.boundary import (
plane_at,
line_at,
point_at,
within_range,
)
Expand Down Expand Up @@ -175,6 +176,18 @@ def create_body_force_am(

return L

def compute_volume(self) -> float:
"""Computes the volume of the structure
Returns:
volume of the structure
"""
dx = ufl.Measure("dx", domain=self.mesh, metadata={"quadrature_degree": 2})
volume = df.fem.assemble_scalar(df.fem.form(1.0 * dx))

return volume


class ExperimentStructure(Experiment):
"""Experiment class for AM structure simulation in the end (after printing).
Expand Down Expand Up @@ -276,59 +289,142 @@ 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
# for each defined case please also define the bottom boundary for the reaction force sensor
# get mesh_points to define boundaries
mesh_points = self.mesh.geometry.x
if self.p["dim"] == 3:
self.min_x = mesh_points[:, 0].min()
self.max_x = mesh_points[:, 0].max()
self.min_y = mesh_points[:, 1].min()
self.max_y = mesh_points[:, 1].max()
self.min_z = mesh_points[:, 2].min()
self.max_z = mesh_points[:, 2].max()
else:
raise ValueError(f"wrong dimension: {self.p['dim']} is not implemented for problem setup")


if self.p["bc_setting"] == "fixed_y":
# loading displacement controlled in y direction at max_y surface, whereas y-z surface at min_y is full 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", self.min_y)
self.logger.debug("check points at y_min %s", mesh_points[mesh_points[:, 1] == self.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),
boundary=plane_at(self.min_y, 'y'),
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
self.logger.info("apply disp where y= %s", self.max_y)
self.logger.debug("check points at max_y: %s", mesh_points[mesh_points[:, 1] == self.max_y])
bc_generator.add_dirichlet_bc(
self.top_displacement,
boundary=plane_at(self.max_y, 'y'),
sub=1,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1,
)

elif 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 sligthly fixed
if self.p["dim"] == 3:
self.logger.info("fix at y= %s", self.min_y)
self.logger.debug("check points at y_min %s", mesh_points[mesh_points[:, 1] == self.min_y])
# fixed bottom in y
bc_generator.add_dirichlet_bc(
df.fem.Constant(
domain=self.mesh, c=0.0
),
boundary=plane_at(self.min_y, 'y'),
sub=1,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1, # surface
)
# one point in all dirc and one in z direction fixed
print('bc points: ', self.min_x, self.min_y, self.min_z, 'and', self.max_x, self.min_y, self.min_z)
bc_generator.add_dirichlet_bc(
np.array([0.0, 0.0, 0.0], dtype=ScalarType),
boundary=point_at([self.min_x, self.min_y, self.min_z]),
method="geometrical",
entity_dim=0, # point
)
bc_generator.add_dirichlet_bc(
df.fem.Constant(
domain=self.mesh, c=0.0
),
boundary=point_at([self.max_x, self.min_y, self.min_z]),
sub=2,
method="geometrical",
entity_dim=0, # point
)

# 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", self.max_y)
self.logger.debug("check points at max_y: %s", mesh_points[mesh_points[:, 1] == self.max_y])
bc_generator.add_dirichlet_bc(
self.top_displacement,
boundary=plane_at(max_y, 1),
boundary=plane_at(self.max_y, 'y'),
sub=1,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1,
)
# 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]
# ),
# method="geometrical",
# entity_dim=self.mesh.topology.dim - 2, # points
# )

# define displacement sensor position for this set-up
self.sensor_location_corner_top = [self.min_x, self.max_y, self.min_z] # corner point
self.sensor_location_middle_endge = [self.min_x, (self.max_y-self.min_y)/2., self.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 sligthly fixed
if self.p["dim"] == 3:
self.logger.info("fix at x=%s", self.min_x)
self.logger.debug("check points at x_min %s", mesh_points[mesh_points[:, 0] == self.min_x])
# fixed in x
bc_generator.add_dirichlet_bc(
df.fem.Constant(
domain=self.mesh, c=0.0
),
boundary=plane_at(self.min_x, 'x'),
sub=0,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1, # surface
)
#
print('bc points: ', self.min_x, self.min_y, self.min_z, 'and', self.min_x, self.max_y, self.min_z)
bc_generator.add_dirichlet_bc(
np.array([0.0, 0.0, 0.0], dtype=ScalarType),
boundary=point_at([self.min_x, self.min_y, self.min_z]),
method="geometrical",
entity_dim=0, # point
)
bc_generator.add_dirichlet_bc(
df.fem.Constant(
domain=self.mesh, c=0.0
),
boundary=point_at([self.min_x, self.max_y, self.min_z]),
sub=2,
method="geometrical",
entity_dim=0, # point
)

# top displacement at top (==max_x) in x direction
self.logger.info("apply disp where x=%s", self.max_x)
self.logger.debug("check points at max_x: %s", mesh_points[mesh_points[:, 0] == self.max_x])
bc_generator.add_dirichlet_bc(
self.top_displacement,
boundary=plane_at(self.max_x, 'x'),
sub=0,
method="geometrical",
entity_dim=self.mesh.topology.dim - 1,
)
# define displacement sensor position for this set-up
self.sensor_location_corner_top = [self.max_x, self.min_y, self.min_z] # corner point
self.sensor_location_middle_endge = [(self.max_x-self.min_x)/2, self.min_y, self.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 +440,36 @@ 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" or self.p["bc_setting"] == "fixed_y":
if self.p["dim"] == 3:
return plane_at(self.min_y, "y")

elif self.p["bc_setting"] == "compr_disp_x":
if self.p["dim"] == 3:
return plane_at(self.min_x, "x")

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

def compute_volume(self) -> float:
"""Computes the volume of the structure
Returns:
volume of the structure
"""
dx = ufl.Measure("dx", domain=self.mesh, metadata={"quadrature_degree": 2})
volume = df.fem.assemble_scalar(df.fem.form(1.0 * dx))

return volume


# def create_body_force(self, v: ufl.argument.Argument) -> ufl.form.Form:
# """defines body force
#
Expand All @@ -363,4 +489,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 @@ -95,9 +99,10 @@ def run(self, mesh_file: Path, out_xdmf: Path) -> None:
if self.pint_parameters["experiment_type"].magnitude == "structure":

experiment = ExperimentStructure(self.pint_parameters)
print('Volume of design', experiment.compute_volume())

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 +112,51 @@ 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())
if experiment.p["bc_setting"] == "compr_disp_x" or experiment.p["bc_setting"] == "compr_disp_y":
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
)

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())
## temporary solution to get equivalent modulus TO BE DELTED
if self.pint_parameters["bc_setting"].magnitude == "compr_disp_y":
force_y = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 1].min()
u_max = d_disp
print('equivalent modulus', force_y/(u_max*0.15)*10**(-6))
elif self.pint_parameters["bc_setting"].magnitude == "compr_disp_x":
force_y = np.array(self.problem.sensors["ReactionForceSensor"].data)[:, 0].min()
u_max = d_disp
print('equivalent modulus', force_y/(u_max*0.15)*10**(-6))

# save sensor output as csv file for postprocessing
if experiment.p["bc_setting"] == "compr_disp_x" or experiment.p["bc_setting"] == "compr_disp_y":
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)
print(df_sensor)
df_sensor.to_csv(out_xdmf.parent / f"{out_xdmf.stem}.csv", index=False)


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 d33d7cc

Please sign in to comment.