diff --git a/src/compas_fea2/UI/viewer/primitives.py b/src/compas_fea2/UI/viewer/primitives.py index 931a409ad..d81a68ed5 100644 --- a/src/compas_fea2/UI/viewer/primitives.py +++ b/src/compas_fea2/UI/viewer/primitives.py @@ -1,7 +1,6 @@ from compas.geometry import Box from compas.geometry import Circle from compas.geometry import Cone -from compas.geometry import Cylinder from compas.geometry import Frame from compas.geometry import Plane @@ -102,12 +101,15 @@ class RollerBCShape(_BCShape): The scale factor to apply when drawing. Default is 1. """ - def __init__(self, xyz, direction=[1, 0, 0], scale=1): + def __init__(self, xyz, direction=[0, 0, 1], scale=1): super(RollerBCShape, self).__init__(xyz, direction, scale) - self.height = 800 * self.scale - p = Plane([self.x, self.y, self.z / 2], [0, 1, 0]) - c = Circle(plane=p, radius=self.height / 2) - self.shape = Cylinder(circle=c, height=self.height) + self.height = 400 * self.scale + self.diameter = 400 * self.scale + # FIXME this is wrong because it should follow the normal + self.plane = Plane([self.x, self.y, self.z - self.height], direction) + self.circle = Circle(frame=Frame.from_plane(self.plane), radius=self.diameter / 2) + self.shape = Cone(radius=self.circle.radius, height=self.height, frame=Frame.from_plane(self.plane)) + class MomentShape(_BCShape): diff --git a/src/compas_fea2/model/__init__.py b/src/compas_fea2/model/__init__.py index aa2cd635d..e2ec38a5c 100644 --- a/src/compas_fea2/model/__init__.py +++ b/src/compas_fea2/model/__init__.py @@ -21,10 +21,12 @@ _Element3D, TetrahedronElement, HexahedronElement, + Face ) from .materials.material import ( _Material, ElasticIsotropic, + ThermalElasticIsotropic, ElasticOrthotropic, ElasticPlastic, Stiff, @@ -123,6 +125,8 @@ HardContactNoFriction, LinearContactFrictionPenalty, HardContactRough, + SurfaceConvection, + SurfaceRadiation ) __all__ = [ @@ -152,6 +156,7 @@ "ConcreteSmearedCrack", "ConcreteDamagedPlasticity", "ElasticIsotropic", + "ThermalElasticIsotropic", "Stiff", "ElasticOrthotropic", "ElasticPlastic", diff --git a/src/compas_fea2/model/bcs.py b/src/compas_fea2/model/bcs.py index ef5b6c221..695f135f9 100644 --- a/src/compas_fea2/model/bcs.py +++ b/src/compas_fea2/model/bcs.py @@ -1,4 +1,5 @@ from typing import Dict +from typing import Optional from compas_fea2.base import FEAData @@ -35,6 +36,8 @@ Restrain rotations around the y axis. zz : bool Restrain rotations around the z axis. +temp : float (False otherwise) + Imposed temperature for heat analysis. components : dict Dictionary with component-value pairs summarizing the boundary condition. axes : str @@ -56,6 +59,7 @@ def __init__(self, axes: str = "global", **kwargs): self._xx = False self._yy = False self._zz = False + self._temp = False @property def x(self) -> bool: @@ -80,6 +84,10 @@ def yy(self) -> bool: @property def zz(self) -> bool: return self._zz + + @property + def temp(self) -> bool: + return self._temp @property def axes(self) -> str: @@ -104,6 +112,7 @@ def __data__(self) -> dict: "xx": self._xx, "yy": self._yy, "zz": self._zz, + "temp": self._temp, } @classmethod @@ -302,3 +311,49 @@ def __init__(self, **kwargs): super().__init__(**kwargs) self._x = False self._z = False + +#=================================================================== +# HEAT ANALYSIS +#=================================================================== + +class _ThermalBoundaryCondition(FEAData): + """Base class for temperature boundary conditions.""" + + __doc__ += docs + + def __init__(self, temp: float = None, **kwargs): + super().__init__(**kwargs) + self._temp = None + + @property + def temp(self) -> Optional[float]: + return self._temp + + @property + def __data__(self) -> Dict[str, any]: + return { + "class": self.__class__.__base__.__name__, + "temp": self._temp, + } + + @classmethod + def __from_data__(cls, data: Dict[str, any]): + return cls( + temp=data.get("temp", None) + ) + +class ImposedTemperature(_ThermalBoundaryCondition): + """Imposed temperature conidtion for heat analysis. + + Additional Parameters +--------------------- +temperature : float + Value of imposed temperature applied + """ + + __doc__ += docs + + def __init__(self, temperature, **kwargs): + super().__init__(**kwargs) + self._temp = temperature + diff --git a/src/compas_fea2/model/elements.py b/src/compas_fea2/model/elements.py index 119f31400..821d38540 100644 --- a/src/compas_fea2/model/elements.py +++ b/src/compas_fea2/model/elements.py @@ -13,23 +13,24 @@ from compas.geometry import Polygon from compas.geometry import Polyhedron from compas.geometry import Vector +from compas.geometry import cross_vectors from compas.geometry import centroid_points from compas.geometry import distance_point_point from compas.itertools import pairwise - from compas_fea2.base import FEAData +from compas_fea2.model.materials.material import ElasticOrthotropic if TYPE_CHECKING: from compas_fea2.problem import Step + from compas_fea2.results import Result + from compas_fea2.results import ShellStressResult + from compas_fea2.results import SolidStressResult from .model import Model from .nodes import Node from .parts import _Part from .sections import _Section from .shapes import Shape - from compas_fea2.results import Result - from compas_fea2.results import ShellStressResult - from compas_fea2.results import SolidStressResult class _Element(FEAData): @@ -85,7 +86,7 @@ class _Element(FEAData): """ - def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, **kwargs): + def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, heat: bool = False, **kwargs): super().__init__(**kwargs) self._part_key = None self._nodes = self._check_nodes(nodes) @@ -98,6 +99,7 @@ def __init__(self, nodes: List["Node"], section: "_Section", implementation: Opt self._volume = None self._results_format = {} self._rigid = rigid + self._heat = heat self._reference_point = None self._shape = None @@ -202,6 +204,10 @@ def reference_point(self) -> "Point": def rigid(self) -> bool: return self._rigid + @property + def heat(self) -> bool: + return self._heat + @property def mass(self) -> float: return self.volume * self.section.material.density @@ -246,8 +252,8 @@ def __from_data__(cls, data): class _Element0D(_Element): """Element with 1 dimension.""" - def __init__(self, nodes: List["Node"], frame: Frame, implementation: Optional[str] = None, rigid: bool = False, **kwargs): - super().__init__(nodes, section=None, implementation=implementation, rigid=rigid, **kwargs) + def __init__(self, nodes: List["Node"], frame: Frame, implementation: Optional[str] = None, **kwargs): + super().__init__(nodes, section=None, implementation=implementation, rigid=False, heat=False, **kwargs) self._frame = frame self._ndim = 0 @@ -268,8 +274,8 @@ class SpringElement(_Element0D): """ - def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, **kwargs): - super().__init__(nodes, section=section, implementation=implementation, rigid=rigid, **kwargs) + def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, **kwargs): + super().__init__(nodes, section=section, implementation=implementation, rigid=False, **kwargs) class LinkElement(_Element0D): @@ -282,7 +288,7 @@ class LinkElement(_Element0D): """ def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, **kwargs): - super().__init__(nodes, section=section, implementation=implementation, rigid=rigid, **kwargs) + super().__init__(nodes, section=section, implementation=implementation, rigid=rigid, heat=False, **kwargs) class _Element1D(_Element): @@ -312,8 +318,10 @@ class _Element1D(_Element): The volume of the element. """ - def __init__(self, nodes: List["Node"], section: "_Section", frame: Optional[Frame] = None, implementation: Optional[str] = None, rigid: bool = False, **kwargs): - super().__init__(nodes, section, implementation=implementation, rigid=rigid, **kwargs) + def __init__( + self, nodes: List["Node"], section: "_Section", frame: Optional[Frame] = None, implementation: Optional[str] = None, rigid: bool = False, heat: bool = False, **kwargs + ): + super().__init__(nodes, section, implementation=implementation, rigid=rigid, heat=heat, **kwargs) if not frame: raise ValueError("Frame is required for 1D elements") self._frame = frame if isinstance(frame, Frame) else Frame(nodes[0].point, Vector(*frame), Vector.from_start_end(nodes[0].point, nodes[-1].point)) @@ -368,8 +376,8 @@ def plot_section(self): self.section.plot() def plot_stress_distribution(self, step: "Step", end: str = "end_1", nx: int = 100, ny: int = 100, *args, **kwargs): # noqa: F821 - """ Plot the stress distribution along the element. - + """Plot the stress distribution along the element. + Parameters ---------- step : :class:`compas_fea2.model.Step` @@ -424,8 +432,8 @@ def forces(self, step: "Step") -> "Result": return r.forces def moments(self, step: "_Step") -> "Result": # noqa: F821 - """ Get the moments result for the element. - + """Get the moments result for the element. + Parameters ---------- step : :class:`compas_fea2.model.Step` @@ -454,8 +462,8 @@ class BeamElement(_Element1D): class TrussElement(_Element1D): """A 1D element that resists axial loads.""" - def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, **kwargs): - super().__init__(nodes, section, frame=[1, 1, 1], implementation=implementation, rigid=rigid, **kwargs) + def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, heat: bool = False, **kwargs): + super().__init__(nodes, section, frame=[1, 1, 1], implementation=implementation, rigid=rigid, heat=heat, **kwargs) class StrutElement(TrussElement): @@ -466,6 +474,51 @@ class TieElement(TrussElement): """A truss element that resists axial tensile loads.""" +class Edge(FEAData): + def __init__(self, nodes: List["Node"], tag: str, element: Optional["_Element"] = None, **kwargs): + super().__init__(**kwargs) + self._nodes = nodes + self._tag = tag + self._line = Line(start=nodes[0].point, end=nodes[1].point) # TODO check when more than 3 nodes + self._registration = element # FIXME: not updated when copying parts + + @property + def nodes(self) -> List["Node"]: + return self._nodes + + @property + def tag(self) -> str: + return self._tag + + @property + def plane(self) -> Plane: + return self._plane + + @property + def element(self) -> Optional["_Element"]: + return self._registration + + @property + def part(self) -> "_Part": + return self.element.part + + @property + def model(self) -> "Model": + return self.element.model + + @property + def centroid(self) -> "Point": + return centroid_points([node.xyz for node in self.nodes]) + + @property + def nodes_key(self) -> List: + return [n._part_key for n in self.nodes] + + @property + def points(self) -> List["Point"]: + return [node.point for node in self.nodes] + + class Face(FEAData): """Element representing a face. @@ -587,29 +640,62 @@ def node_area(self) -> float: class _Element2D(_Element): """Element with 2 dimensions.""" - __doc__ += _Element.__doc__ + __doc__ = __doc__ or "" + __doc__ += _Element.__doc__ or "" __doc__ += """ Additional Parameters --------------------- + frame : list, optional. + Local Y axis in global coordinates (longitudinal axis). + This can be used to define the local orientation for anisotrop material. faces : [:class:`compas_fea2.model.elements.Face] The faces of the element. face_indices : dict Dictionary providing for each face the node indices. For example: {'s1': (0,1,2), ...} + edges : [:class:`compas_fea2.model.elements.Edge] + The edges of the element. + edge_indices : dict + Dictionary providing for each edge the node indices. For example: + {'e1': (0,1), ...} """ - def __init__(self, nodes: List["Node"], section: Optional["_Section"] = None, implementation: Optional[str] = None, rigid: bool = False, **kwargs): + def __init__( + self, + nodes: List["Node"], + section: Optional["_Section"] = None, + implementation: Optional[str] = None, + rigid: bool = False, + heat: bool = False, + frame: Optional[Frame] = None, + **kwargs, + ): super().__init__( nodes=nodes, section=section, implementation=implementation, rigid=rigid, + heat=heat, **kwargs, ) - self._faces = None self._face_indices = None + self._edges = None + self._edges_indices = None self._ndim = 2 + if frame: + if isinstance(frame, (Vector, List, Tuple)): + compas_polygon = Polygon(points=[node.point for node in nodes]) + # the third axis is built from the y-axis (frame input) and z-axis (normal of the element) + x_axis = cross_vectors(frame, compas_polygon.normal) + frame = Frame(nodes[0].point, Vector(*x_axis), Vector(*frame)) + self._frame = frame + else: + raise ValueError("The frame must be a Frame object or a Vector object.") + + else: + if isinstance(self.section.material, ElasticOrthotropic): + raise ValueError("For orthotropic or anisotropic materials, the frame must implemented to define the local orientation of the element.") @property def nodes(self) -> List["Node"]: @@ -620,6 +706,14 @@ def nodes(self, value: List["Node"]): self._nodes = self._check_nodes(value) self._faces = self._construct_faces(self._face_indices) + @property + def edge_indices(self) -> Optional[Dict[str, Tuple[int]]]: + return self._edgee_indices + + @property + def edges(self) -> Optional[List[Face]]: + return self._edges + @property def face_indices(self) -> Optional[Dict[str, Tuple[int]]]: return self._face_indices @@ -656,6 +750,22 @@ def _construct_faces(self, face_indices: Dict[str, Tuple[int]]) -> List[Face]: """ return [Face(nodes=itemgetter(*indices)(self.nodes), tag=name, element=self) for name, indices in face_indices.items()] + def _construct_edges(self, edge_indices: Dict[str, Tuple[int]]) -> List[Edge]: + """Construct the face-nodes dictionary. + + Parameters + ---------- + face_indices : dict + Dictionary providing for each face the node indices. For example: + {'s1': (0,1,2), ...} + + Returns + ------- + dict + Dictionary with edge names and the corresponding nodes. + """ + return [Edge(nodes=itemgetter(*indices)(self.nodes), tag=name, element=self) for name, indices in edge_indices.items()] + def stress_results(self, step: "Step") -> "Result": """Get the stress results for the element. @@ -682,17 +792,21 @@ class ShellElement(_Element2D): """ - def __init__(self, nodes: List["Node"], section: Optional["_Section"] = None, implementation: Optional[str] = None, rigid: bool = False, **kwargs): + def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid: bool = False, heat: bool = False, **kwargs): super().__init__( nodes=nodes, section=section, implementation=implementation, rigid=rigid, + heat=heat, **kwargs, ) self._face_indices = {"SPOS": tuple(range(len(nodes))), "SNEG": tuple(range(len(nodes)))[::-1]} self._faces = self._construct_faces(self._face_indices) + self._edges_indices = {f"s{i+1}": (i, i + 1) if i < len(nodes) - 1 else (i, 0) for i in range(len(nodes))} + # self._edge_indices[f"s{len(nodes)-1}"]= (len(nodes)-1, 0) + self._edges = self._construct_edges(self._edges_indices) @property def results_cls(self) -> "Result": @@ -724,11 +838,13 @@ class _Element3D(_Element): """ - def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, **kwargs): + def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, frame: Optional[List] = None, rigid=False, heat: bool = False, **kwargs): super().__init__( nodes=nodes, section=section, implementation=implementation, + rigid=rigid, + heat=heat, **kwargs, ) self._face_indices = None @@ -919,11 +1035,13 @@ class PentahedronElement(_Element3D): class HexahedronElement(_Element3D): """A Solid cuboid element with 6 faces (extruded rectangle).""" - def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, **kwargs): + def __init__(self, nodes: List["Node"], section: "_Section", implementation: Optional[str] = None, rigid=False, heat: bool = False,**kwargs): super().__init__( nodes=nodes, section=section, implementation=implementation, + rigid=rigid, + heat=heat, **kwargs, ) self._faces_indices = { diff --git a/src/compas_fea2/model/groups.py b/src/compas_fea2/model/groups.py index 22f89f48e..8a0b42ecd 100644 --- a/src/compas_fea2/model/groups.py +++ b/src/compas_fea2/model/groups.py @@ -445,6 +445,104 @@ def add_elements(self, elements): """ return self._add_members(elements) +class EdgesGroup(_Group): + """Base class elements edges groups. + + Parameters + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + edges : Set[:class:`compas_fea2.model.Edge`] + The Edges belonging to the group. + + Attributes + ---------- + name : str + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + edges : Set[:class:`compas_fea2.model.Edge`] + The Edges belonging to the group. + nodes : Set[:class:`compas_fea2.model.Node`] + The Nodes of the edges belonging to the group. + part : :class:`compas_fea2.model._Part` + The part where the group is registered, by default `None`. + model : :class:`compas_fea2.model.Model` + The model where the group is registered, by default `None`. + + Notes + ----- + EdgesGroups are registered to the same :class:`compas_fea2.model.Part` as the + elements of its edges. + + """ + + def __init__(self, edges, **kwargs): + super().__init__(members=edges, **kwargs) + + @property + def __data__(self): + data = super().__data__ + data.update({"edges": list(self.edges)}) + return data + + @classmethod + def __from_data__(cls, data): + obj = cls(edges=set(data["edges"])) + obj._registration = data["registration"] + return obj + + @property + def part(self): + return self._registration + + @property + def model(self): + return self.part._registration + + @property + def edges(self): + return self._members + + @property + def nodes(self): + nodes_set = set() + for edge in self.edges: + for node in edge.nodes: + nodes_set.add(node) + return nodes_set + + def add_edge(self, edge): + """ + Add a face to the group. + + Parameters + ---------- + face : :class:`compas_fea2.model.Face` + The face to add. + + Returns + ------- + :class:`compas_fea2.model.Face` + The face added. + """ + return self._add_member(edge) + + def add_edges(self, edges): + """ + Add multiple faces to the group. + + Parameters + ---------- + faces : [:class:`compas_fea2.model.Face`] + The faces to add. + + Returns + ------- + [:class:`compas_fea2.model.Face`] + The faces added. + """ + return self._add_members(edges) class FacesGroup(_Group): """Base class elements faces groups. diff --git a/src/compas_fea2/model/ics.py b/src/compas_fea2/model/ics.py index a92b08cf0..80e7e7e0e 100644 --- a/src/compas_fea2/model/ics.py +++ b/src/compas_fea2/model/ics.py @@ -75,8 +75,8 @@ def __from_data__(cls, data): @classmethod def from_file(cls, path, **kwargs): - cls._path = path - return cls(temperature=None, path=path, **kwargs) + return NotImplementedError( + "InitialTemperatureField is not implemented for the current backend. ") class InitialStressField(_InitialCondition): diff --git a/src/compas_fea2/model/interactions.py b/src/compas_fea2/model/interactions.py index c0f0d9b17..71baa2be6 100644 --- a/src/compas_fea2/model/interactions.py +++ b/src/compas_fea2/model/interactions.py @@ -195,3 +195,114 @@ class HardContactRough(Contact): def __init__(self, **kwargs) -> None: super().__init__(normal="HARD", tangent="ROUGH", **kwargs) + + +# ------------------------------------------------------------------------------ +# THERMAL INTERACTION +# ------------------------------------------------------------------------------ + +class ThermalInteraction(_Interaction): + """General thermal interaction of a part with the exterior. + + Parameters + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + temperature : float or [float] + Constant temperature value or temperature evolution. + + + Attributes + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + temperature : float or [float] + Constant temperature value or temperature evolution. + + """ + + def __init__(self, temperature, **kwargs): + super(_Interaction, self).__init__(**kwargs) + self._temperature = temperature + + @property + def temperature(self): + return self._temperature + +class SurfaceConvection(ThermalInteraction): + """Convection. + + Parameters + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + h : float + Convection coefficient. + surface : FaceGroup + Convection surface. + temperature : float or [float] + Constant temperature value or temperature evolution. + + Attributes + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + temperature : float or [float] + Constant temperature value or temperature evolution. + + """ + + def __init__(self, surface, h, temperature, **kwargs): + super().__init__(temperature=temperature, **kwargs) + self._h = h + self._surface = surface + + @property + def h(self): + return self._h + + @property + def surface(self): + return self._surface + +class SurfaceRadiation(ThermalInteraction): + """Radiation. + + Parameters + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + eps : float + Radiation coefficient. + surface : FaceGroup + Convection surface. + temperature : float or [float] + Constant temperature value or temperature evolution. + + Attributes + ---------- + name : str, optional + Uniqe identifier. If not provided it is automatically generated. Set a + name if you want a more human-readable input file. + temperature : float or [float] + Constant temperature value or temperature evolution. + + """ + + def __init__(self, surface, eps, temperature, **kwargs): + super().__init__(temperature=temperature, **kwargs) + self._eps = eps + self._surface = surface + + @property + def eps(self): + return self._eps + + @property + def surface(self): + return self._surface \ No newline at end of file diff --git a/src/compas_fea2/model/materials/material.py b/src/compas_fea2/model/materials/material.py index 4482f7154..7bf355493 100644 --- a/src/compas_fea2/model/materials/material.py +++ b/src/compas_fea2/model/materials/material.py @@ -156,12 +156,23 @@ def __init__( self.Ey = Ey self.Ez = Ez self.vxy = vxy + self.vyx = vxy*Ey/Ex self.vyz = vyz + self.vzy = vyz*Ez/Ey self.vzx = vzx + self.vxz = vzx*Ex/Ez self.Gxy = Gxy self.Gyz = Gyz self.Gzx = Gzx + #the equations below must be verified by an orthotropic elactic material + check_list = [self.Ex>0, self.Ey>0, self.Ez>0, self.Gxy>0, self.Gyz>0, self.Gzx>0, + self.vxy<(Ex/Ey)**0.5, vzx<(Ez/Ex)**0.5,vyz<(Ey/Ez)**0.5, + 1-self.vxy*self.vyx-self.vyz*self.vzy-self.vzx*self.vxz-2*self.vyx*self.vzy*self.vxz>0] + for check in check_list: + if not(check): + raise ValueError('The mechanical values do not respect the material stability criteria.') + @property def __data__(self): data = super().__data__ @@ -329,7 +340,7 @@ class ElasticPlastic(ElasticIsotropic): in the form of strain/stress value pairs. """ - def __init__(self, *, E: float, v: float, density: float, strain_stress: list[tuple[float, float]], expansion: Optional[float] = None, **kwargs): + def __init__(self, E: float, v: float, density: float, strain_stress: list[tuple[float, float]], expansion: float = None, **kwargs): super().__init__(E=E, v=v, density=density, expansion=expansion, **kwargs) self.strain_stress = strain_stress @@ -385,3 +396,74 @@ class UserMaterial(FEAData): def __init__(self, **kwargs): super().__init__(**kwargs) raise NotImplementedError("This class is not available for the selected backend plugin") + +# ============================================================================== +# Heat Material +# ============================================================================== + +class ThermalElasticIsotropic(ElasticIsotropic): + """Thermal isotropic material for heat analysis. + + Parameters + ---------- + k : float + Thermal conductivity. + c : float + Specific heat capacity. + + Attributes + ---------- + k : float + Thermal conductivity. + c : float + Specific heat capacity. + + """ + + + + def __init__(self, k: float, c:float, E: float, v: float, density: float, expansion: float = None, **kwargs): + super().__init__(E=E, v=v,density=density, expansion=expansion, **kwargs) + self._k = k + self._c = c + + @property + def __data__(self): + data = super().__data__ + data.update( + { + "k": self.k, + "c": self.c, + } + ) + return data + + @classmethod + def __from_data__(cls, data): + return cls(data["k"], data["c"],data["E"], data["v"], data["density"], data["expansion"]) + + def __str__(self) -> str: + return """ +Thermal ElasticIsotropic Material +------------------------- +name : {} +density : {} +expansion : {} + +E : {} +v : {} +G : {} + +k : {} +c : {} +""".format( + self.name, self.density, self.expansion, self.E, self.v, self.G, self.k, self.c + ) + + @property + def k(self) -> float: + return self._k + + @property + def c(self) -> float: + return self._c \ No newline at end of file diff --git a/src/compas_fea2/model/materials/timber.py b/src/compas_fea2/model/materials/timber.py index 9a9c1744e..a2a39d6ec 100644 --- a/src/compas_fea2/model/materials/timber.py +++ b/src/compas_fea2/model/materials/timber.py @@ -1,23 +1,111 @@ -from .material import _Material +from compas_fea2.units import UnitRegistry +from compas_fea2.units import units as u +from .material import ElasticOrthotropic -class Timber(_Material): - """Base class for Timber material""" +vLT_softwood = 0.3 +vTT_softwood = 0.4 +vLT_hardwood = 0.4 +vTT_hardwood = 0.5 - def __init__(self, *, density, **kwargs): - """ - Parameters - ---------- - density : float - Density of the timber material [kg/m^3]. - name : str, optional - Name of the material. - """ - super().__init__(density=density, **kwargs) + +class Timber(ElasticOrthotropic): + + """Base class for Timber material (elastic orthotropic behaviour). + The longitudinal axis (along the grain) is defined along the y-axis. + + Parameters + ---------- + fmk : float + Bending resistance. + ft0k : float + Characteristic tensile strength along the grain. + fc0k : float + Characteristic compressive strength along the grain. + ft90k : float + Characteristic tensile strength perpendicular to the grain. + fc90k : float + Characteristic tensile strength along the grain. + fvk : float + Shear resistance. + vLT : float + Value of Poisson ratio longitudinal/transverse. + vTT : float + Value of Poisson ratio transverse/transverse. + E0mean : float + Mean value of modulus of elasticity parallel. + E90mean : float + Mean value of modulus of elasticity perpendicular. + Gmean : float + Mean value of shear modulus. + density : float + Mean density of the timber material [kg/m^3]. + name : str, optional + Name of the material. + + Attributes : + ------------ + + fmk : float + Bending resistance. + ft0k : float + Characteristic tensile strength along the grain. + fc0k : float + Characteristic compressive strength along the grain. + ft90k : float + Characteristic tensile strength perpendicular to the grain. + fc90k : float + Characteristic tensile strength along the grain. + fvk : float + Shear resistance. + Ex : float + Young's modulus Ex in x direction (EN383 - E90 mean). + Ey : float + Young's modulus Ey in y direction (EN383 - E0 mean). + Ez : float + Young's modulus Ez in z direction (EN383 - E90 mean). + vxy : float + Poisson's ratio vxy in x-y directions. + vyz : float + Poisson's ratio vyz in y-z directions. + vzx : float + Poisson's ratio vzx in z-x directions. + Gxy : float + Shear modulus Gxy in x-y direction (EN383 - G mean). + Gyz : float + Shear modulus Gyz in y-z directions (EN383 - G mean). + Gzx : float + Shear modulus Gzx in z-x directions (EN383 - G mean). + + """ + __doc__+=ElasticOrthotropic.__doc__ + + def __init__(self, fmk, ft0k, fc0k, ft90k, fc90k, fvk, vLT, vTT, E0mean, E90mean, Gmean, densityk, density, **kwargs): + + super().__init__(Ex=E90mean, Ey=E0mean, Ez=E90mean, vxy=vLT*E90mean/E0mean, vyz=vLT, vzx=vTT , Gxy=Gmean, Gyz=Gmean, Gzx=Gmean, density=density, **kwargs) + self.fmk = fmk + self.ft0k = ft0k + self.fc0k = fc0k + self.ft90k = ft90k + self.fc90k =fc90k + self.fvk = fvk + self.densityk = densityk @property def __data__(self): return { + "fmk": self.fmk, + "ft0k": self.ft0k, + "fc0k": self.fc0k, + "ft90k": self.ft90k, + "fc90k": self.fc90k, + "fvk": self.fvk, + "vLT": self.vyz, + "vTT": self.vxy, + "E0mean": self.Ex, + "E90mean": self.Ey, + "Gmean": self.Gxy, + "densityk": self.densityk, "density": self.density, "name": self.name, } @@ -25,6 +113,400 @@ def __data__(self): @classmethod def __from_data__(cls, data): return cls( + fmk=data["fmk"], + ft0k=data["ft0k"], + fc0k=data["fc0k"], + ft90k=data["ft90k"], + fc90k=data["fc90k"], + fvk=data["fvk"], + vLT=data["vLT"], + vTT=data["vTT"], + E0mean=data["E0mean"], + E90mean=data["E90mean"], + Gmean=data["Gmean"], + densityk=data["densityk"], density=data["density"], name=data["name"], ) + + # --- Softwood Classes (C-classes) --- + @classmethod + def C14(cls, units=None, **kwargs): + """ + Softwood C14. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=14*units.MPa, ft0k=8*units.MPa, fc0k=16*units.MPa, ft90k=0.4*units.MPa, fc90k=2*units.MPa, + fvk=3*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=7*units.GPa, E90mean=0.23*units.GPa, Gmean=0.44*units.GPa, + densityk=290*units("kg/m**3"), density=350*units("kg/m**3"), **kwargs) + + @classmethod + def C16(cls, units=None, **kwargs): + """ + Softwood C16. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=16*units.MPa, ft0k=10*units.MPa, fc0k=17*units.MPa, ft90k=0.4*units.MPa, fc90k=2.2*units.MPa, + fvk=3.2*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=8*units.GPa, E90mean=0.27*units.GPa, Gmean=0.5*units.GPa, + densityk=310*units("kg/m**3"), density=370*units("kg/m**3"), **kwargs) + + @classmethod + def C18(cls, units=None, **kwargs): + """ + Softwood C18. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=18*units.MPa, ft0k=11*units.MPa, fc0k=18*units.MPa, ft90k=0.4*units.MPa, fc90k=2.2*units.MPa, + fvk=3.4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=9*units.GPa, E90mean=0.3*units.GPa, Gmean=0.56*units.GPa, + densityk=320*units("kg/m**3"), density=380*units("kg/m**3"), **kwargs) + + @classmethod + def C20(cls, units=None, **kwargs): + """ + Softwood C20. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=20*units.MPa, ft0k=12*units.MPa, fc0k=19*units.MPa, ft90k=0.4*units.MPa, fc90k=2.3*units.MPa, + fvk=3.6*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=9.5*units.GPa, E90mean=0.32*units.GPa, Gmean=0.59*units.GPa, + densityk=330*units("kg/m**3"), density=390*units("kg/m**3"), **kwargs) + + @classmethod + def C22(cls, units=None, **kwargs): + """ + Softwood C22. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=22*units.MPa, ft0k=13*units.MPa, fc0k=20*units.MPa, ft90k=0.4*units.MPa, fc90k=2.4*units.MPa, + fvk=3.8*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=10*units.GPa, E90mean=0.33*units.GPa, Gmean=0.63*units.GPa, + densityk=340*units("kg/m**3"), density=410*units("kg/m**3"), **kwargs) + + @classmethod + def C24(cls, units=None, **kwargs): + """ + Softwood C24. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=24*units.MPa, ft0k=14*units.MPa, fc0k=21*units.MPa, ft90k=0.4*units.MPa, fc90k=2.5*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=11*units.GPa, E90mean=0.37*units.GPa, Gmean=0.69*units.GPa, + densityk=350*units("kg/m**3"), density=420*units("kg/m**3"), **kwargs) + + @classmethod + def C27(cls, units=None, **kwargs): + """ + Softwood C27. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=27*units.MPa, ft0k=16*units.MPa, fc0k=22*units.MPa, ft90k=0.4*units.MPa, fc90k=2.6*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=11.5*units.GPa, E90mean=0.38*units.GPa, Gmean=0.72*units.GPa, + densityk=370*units("kg/m**3"), density=450*units("kg/m**3"), **kwargs) + + @classmethod + def C30(cls, units=None, **kwargs): + """ + Softwood C30. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=30*units.MPa, ft0k=18*units.MPa, fc0k=23*units.MPa, ft90k=0.4*units.MPa, fc90k=2.7*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=12*units.GPa, E90mean=0.4*units.GPa, Gmean=0.75*units.GPa, + densityk=380*units("kg/m**3"), density=460*units("kg/m**3"), **kwargs) + + @classmethod + def C35(cls, units=None, **kwargs): + """ + Softwood C35. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=35*units.MPa, ft0k=21*units.MPa, fc0k=25*units.MPa, ft90k=0.4*units.MPa, fc90k=2.8*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=13*units.GPa, E90mean=0.43*units.GPa, Gmean=0.81*units.GPa, + densityk=400*units("kg/m**3"), density=480*units("kg/m**3"), **kwargs) + + @classmethod + def C40(cls, units=None, **kwargs): + """ + Softwood C40. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=40*units.MPa, ft0k=24*units.MPa, fc0k=26*units.MPa, ft90k=0.4*units.MPa, fc90k=2.9*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=14*units.GPa, E90mean=0.47*units.GPa, Gmean=0.88*units.GPa, + densityk=420*units("kg/m**3"), density=500*units("kg/m**3"), **kwargs) + + @classmethod + def C45(cls, units=None, **kwargs): + """ + Softwood C45. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=45*units.MPa, ft0k=27*units.MPa, fc0k=27*units.MPa, ft90k=0.4*units.MPa, fc90k=3.1*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=15*units.GPa, E90mean=0.5*units.GPa, Gmean=0.94*units.GPa, + densityk=440*units("kg/m**3"), density=520*units("kg/m**3"), **kwargs) + + @classmethod + def C50(cls, units=None, **kwargs): + """ + Softwood C50. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled softwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=50*units.MPa, ft0k=30*units.MPa, fc0k=29*units.MPa, ft90k=0.4*units.MPa, fc90k=3.2*units.MPa, + fvk=4*units.MPa, vLT=vLT_softwood, vTT=vTT_softwood, E0mean=16*units.GPa, E90mean=0.53*units.GPa, Gmean=1*units.GPa, + densityk=460*units("kg/m**3"), density=550*units("kg/m**3"), **kwargs) + + # --- Hardwood Classes (D-classes) --- + @classmethod + def D18(cls, units=None, **kwargs): + """ + Hardwood D18. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=18*units.MPa, ft0k=11*units.MPa, fc0k=18*units.MPa, ft90k=0.6*units.MPa, fc90k=7.5*units.MPa, + fvk=3.4*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=9.5*units.GPa, E90mean=0.63*units.GPa, Gmean=0.59*units.GPa, + densityk=475*units("kg/m**3"), density=570*units("kg/m**3"), **kwargs) + + @classmethod + def D24(cls, units=None, **kwargs): + """ + Hardwood D24. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=24*units.MPa, ft0k=14*units.MPa, fc0k=21*units.MPa, ft90k=0.6*units.MPa, fc90k=7.8*units.MPa, + fvk=4*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=10*units.GPa, E90mean=0.67*units.GPa, Gmean=0.62*units.GPa, + densityk=485*units("kg/m**3"), density=580*units("kg/m**3"), **kwargs) + + @classmethod + def D30(cls, units=None, **kwargs): + """ + Hardwood D30. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=30*units.MPa, ft0k=18*units.MPa, fc0k=23*units.MPa, ft90k=0.6*units.MPa, fc90k=8*units.MPa, + fvk=4*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=11*units.GPa, E90mean=0.73*units.GPa, Gmean=0.69*units.GPa, + densityk=530*units("kg/m**3"), density=640*units("kg/m**3"), **kwargs) + + @classmethod + def D35(cls, units=None, **kwargs): + """ + Hardwood D35. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=35*units.MPa, ft0k=21*units.MPa, fc0k=25*units.MPa, ft90k=0.6*units.MPa, fc90k=8.1*units.MPa, + fvk=4*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=12*units.GPa, E90mean=0.8*units.GPa, Gmean=0.75*units.GPa, + densityk=540*units("kg/m**3"), density=650*units("kg/m**3"), **kwargs) + + @classmethod + def D40(cls, units=None, **kwargs): + """ + Hardwood D40. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=40*units.MPa, ft0k=24*units.MPa, fc0k=26*units.MPa, ft90k=0.6*units.MPa, fc90k=8.3*units.MPa, + fvk=4*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=13*units.GPa, E90mean=0.86*units.GPa, Gmean=0.81*units.GPa, + densityk=550*units("kg/m**3"), density=660*units("kg/m**3"), **kwargs) + + @classmethod + def D50(cls, units=None, **kwargs): + """ + Hardwood D50. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=50*units.MPa, ft0k=30*units.MPa, fc0k=29*units.MPa, ft90k=0.6*units.MPa, fc90k=9.3*units.MPa, + fvk=4*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=14*units.GPa, E90mean=0.93*units.GPa, Gmean=0.88*units.GPa, + densityk=620*units("kg/m**3"), density=750*units("kg/m**3"), **kwargs) + + @classmethod + def D60(cls, units=None, **kwargs): + """ + Hardwood D60. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=60*units.MPa, ft0k=36*units.MPa, fc0k=32*units.MPa, ft90k=0.6*units.MPa, fc90k=10.5*units.MPa, + fvk=4.5*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=17*units.GPa, E90mean=1.13*units.GPa, Gmean=1.06*units.GPa, + densityk=700*units("kg/m**3"), density=840*units("kg/m**3"), **kwargs) + + @classmethod + def D70(cls, units=None, **kwargs): + """ + Hardwood D70. + + Returns + ------- + :class:`compas_fea2.model.material.ElasticTimber` + The precompiled hardwood material. + """ + if not units: + units = u(system="SI_mm") + elif not isinstance(units, UnitRegistry): + units = u(system=units) + + return cls(fmk=70*units.MPa, ft0k=42*units.MPa, fc0k=34*units.MPa, ft90k=0.6*units.MPa, fc90k=13.5*units.MPa, + fvk=5*units.MPa, vLT=vLT_hardwood, vTT=vTT_hardwood, E0mean=20*units.GPa, E90mean=1.33*units.GPa, Gmean=1.25*units.GPa, + densityk=900*units("kg/m**3"), density=1080*units("kg/m**3"), **kwargs) \ No newline at end of file diff --git a/src/compas_fea2/model/model.py b/src/compas_fea2/model/model.py index d0e39d416..35f2b5bcd 100644 --- a/src/compas_fea2/model/model.py +++ b/src/compas_fea2/model/model.py @@ -8,8 +8,8 @@ from pathlib import Path from typing import Optional from typing import Set -from typing import Union from typing import Tuple +from typing import Union from compas.datastructures import Graph from compas.geometry import Box @@ -25,6 +25,7 @@ import compas_fea2 from compas_fea2.base import FEAData from compas_fea2.model.bcs import _BoundaryCondition +from compas_fea2.model.bcs import _ThermalBoundaryCondition from compas_fea2.model.connectors import Connector from compas_fea2.model.constraints import _Constraint from compas_fea2.model.elements import _Element @@ -34,6 +35,7 @@ from compas_fea2.model.groups import PartsGroup from compas_fea2.model.groups import _Group from compas_fea2.model.ics import _InitialCondition +from compas_fea2.model.interactions import ThermalInteraction from compas_fea2.model.interfaces import Interface from compas_fea2.model.materials.material import _Material from compas_fea2.model.nodes import Node @@ -101,6 +103,7 @@ def __init__(self, description: Optional[str] = None, author: Optional[str] = No self._materials: Set[_Material] = set() self._sections: Set[_Section] = set() self._bcs = {} + self._tbcs = {} self._ics = {} self._interfaces: Set[Interface] = set() self._connectors: Set[Connector] = set() @@ -108,6 +111,8 @@ def __init__(self, description: Optional[str] = None, author: Optional[str] = No self._partsgroups: Set[PartsGroup] = set() self._groups: Set[_Group] = set() self._problems: Set[Problem] = set() + self._thermalinteraction: Set[_Group] = set() + self._amplitudes = set() self._constants: dict = {"g": None} @@ -118,6 +123,7 @@ def __data__(self): "author": self.author, "parts": [part.__data__ for part in self.parts], "bcs": {bc.__data__: [node.__data__ for node in nodes] for bc, nodes in self.bcs.items()}, + "tbcs": {tbc.__data__: [node.__data__ for node in nodes] for tbc, nodes in self.tbcs.items()}, "ics": {ic.__data__: [node.__data__ for node in nodes] for ic, nodes in self.ics.items()}, "constraints": [constraint.__data__ for constraint in self.constraints], "partgroups": [group.__data__ for group in self.partgroups], @@ -151,6 +157,10 @@ def __from_data__(cls, data): for bc_data, nodes_data in data.get("bcs", {}).items(): model._bcs[bc_classes[bc_data["class"]].__from_data__(bc_data)] = [Node.__from_data__(node_data) for node_data in nodes_data] + bct_classes = {cls.__name__: cls for cls in _BoundaryCondition.__subclasses__()} + for bct_data, nodes_data in data.get("bcs", {}).items(): + model._tbcs[bct_classes[bct_data["class"]].__from_data__(bct_data)] = [Node.__from_data__(node_data) for node_data in nodes_data] + ic_classes = {cls.__name__: cls for cls in _InitialCondition.__subclasses__()} for ic_data, nodes_data in data.get("ics", {}).items(): model._ics[ic_classes[ic_data["class"]].__from_data__(ic_data)] = [Node.__from_data__(node_data) for node_data in nodes_data] @@ -211,6 +221,10 @@ def partgroups(self) -> Set[PartsGroup]: def bcs(self) -> dict: return self._bcs + @property + def tbcs(self) -> dict: + return self._tbcs + @property def ics(self) -> dict: return self._ics @@ -262,6 +276,14 @@ def interfaces(self) -> Set[Interface]: @property def interactions(self) -> Set[Interface]: return self.interfaces.group_by(lambda x: getattr(x, "behavior")) + + @property + def thermalinteraction(self): + return _Group(self._thermalinteraction) + + @property + def amplitudes(self): + return self._amplitudes @property def problems(self) -> Set[Problem]: @@ -1024,7 +1046,46 @@ def group_parts_where(self, attr: str, value: Union[str, int, float]) -> PartsGr # ========================================================================= # BCs methods # ========================================================================= + def add_tbcs(self, bct: _ThermalBoundaryCondition, nodes: Union[list[Node], NodesGroup], axes: str = "global") -> _BoundaryCondition: + """Add a :class:`compas_fea2.model._BoundaryCondition` to the model. + + Parameters + ---------- + bct : :class:`compas_fea2.model._ThermalBoundaryCondition` + nodes : list[:class:`compas_fea2.model.Node`] or :class:`compas_fea2.model.NodesGroup` + + Returns + ------- + :class:`compas_fea2.model._ThermalBoundaryCondition` + + """ + if isinstance(nodes, _Group): + nodes = nodes._members + + if isinstance(nodes, Node): + nodes = [nodes] + + if not isinstance(bct, _ThermalBoundaryCondition): + raise TypeError("{!r} is not a thermal Boundary Condition.".format(bct)) + + for node in nodes: + if not isinstance(node, Node): + raise TypeError("{!r} is not a Node.".format(node)) + if not node.part: + raise ValueError("{!r} is not registered to any part.".format(node)) + elif node.part not in self.parts: + raise ValueError("{!r} belongs to a part not registered to this model.".format(node)) + if isinstance(node.part, RigidPart): + if len(nodes) != 1 or not node.is_reference: + raise ValueError("For rigid parts bundary conditions can be assigned only to the reference point") + node._tbc = bct + + bct._key = len(self.bcs) + self._tbcs[bct] = set(nodes) + bct._registration = self + return bct + def add_bcs(self, bc: _BoundaryCondition, nodes: Union[list[Node], NodesGroup], axes: str = "global") -> _BoundaryCondition: """Add a :class:`compas_fea2.model._BoundaryCondition` to the model. @@ -1478,10 +1539,11 @@ def extract_interfaces(self, max_origin_distance=0.01, tol=1e-6) -> list[Tuple[P List of interfaces extracted from the model. """ - from compas.geometry import distance_point_plane_signed import itertools from typing import List + from compas.geometry import distance_point_plane_signed + # --- Nested helper function to check coplanarity --- def _are_two_planes_coplanar(pln1: Plane, pln2: Plane, tol: float) -> bool: """Checks if two COMPAS planes are geometrically the same.""" @@ -1505,6 +1567,28 @@ def _are_two_planes_coplanar(pln1: Plane, pln2: Plane, tol: float) -> bool: return coplanar_pairs + def add_thermalinteraction(self, interaction): + """Add a :class:`compas_fea2.model.ThermalInteraction` to the model. + + Parameters + ---------- + interface : :class:`compas_fea2.model.Interface` + The interface object to add to the model. + + Returns + ------- + :class:`compas_fea2.model.Interface` + + """ + if not isinstance(interaction, ThermalInteraction): + raise TypeError("{!r} is not a thermal interaction.".format(interaction)) + self._thermalinteraction.add(interaction) + if isinstance(interaction.surface, _Group): + self.add_group(interaction.surface) + if interaction.temperature.amplitude: + self._amplitudes.add(interaction.temperature.amplitude) + interaction._registration = self + # ============================================================================== # Summary # ============================================================================== diff --git a/src/compas_fea2/model/nodes.py b/src/compas_fea2/model/nodes.py index a5628eb83..923265fe3 100644 --- a/src/compas_fea2/model/nodes.py +++ b/src/compas_fea2/model/nodes.py @@ -86,6 +86,7 @@ def __init__(self, xyz: List[float], mass: Optional[float] = None, temperature: self._z = xyz[2] self._bc = None + self._tbc = None self._dof = {"x": True, "y": True, "z": True, "xx": True, "yy": True, "zz": True} self._mass = mass if isinstance(mass, list) else list([mass] * 6) @@ -236,6 +237,10 @@ def dof(self) -> Dict[str, bool]: @property def bc(self): return self._bc + + @property + def tbc(self): + return self._tbc @property def on_boundary(self) -> Optional[bool]: @@ -339,8 +344,10 @@ def reactions(self) -> Dict: def results_cls(self): from compas_fea2.results import DisplacementResult from compas_fea2.results import ReactionResult + from compas_fea2.results import TemperatureResult return { "u": DisplacementResult, "rf": ReactionResult, + "t" : TemperatureResult } diff --git a/src/compas_fea2/model/parts.py b/src/compas_fea2/model/parts.py index 271208188..66b366018 100644 --- a/src/compas_fea2/model/parts.py +++ b/src/compas_fea2/model/parts.py @@ -18,6 +18,8 @@ from compas.geometry import Frame from compas.geometry import Plane from compas.geometry import Point +from compas.geometry import Polyline +from compas.geometry import Line from compas.geometry import Line from compas.geometry import Scale from compas.geometry import Transformation @@ -42,6 +44,7 @@ from .elements import _Element1D from .elements import _Element2D from .elements import _Element3D +from .groups import EdgesGroup from .groups import ElementsGroup from .groups import FacesGroup from .groups import MaterialsGroup @@ -248,6 +251,10 @@ def points_sorted(self) -> List[Point]: def elements(self) -> ElementsGroup: return ElementsGroup(self._elements) + @property + def edges(self) -> EdgesGroup: + return EdgesGroup([edge for element in self.elements for edge in element.edges]) + @property def faces(self) -> FacesGroup: return FacesGroup([face for element in self.elements for face in element.faces]) @@ -663,7 +670,7 @@ def elements_by_dimension(self, dimension: int = 1) -> Iterable[_Element]: # ========================================================================= @classmethod - def from_compas_lines_discretized(cls, lines, targetlength, element_class, section, frame, name=None, **kwargs): + def from_compas_lines_discretized(cls, lines, targetlength, element_model, section, frame, name=None, **kwargs): """Generate a discretized model from a list of :class:`compas.geometry.Line`. Parameters @@ -677,8 +684,8 @@ def from_compas_lines_discretized(cls, lines, targetlength, element_class, secti The section to be assigned to the elements, by default None. element_model : str, optional Implementation model for the element, by default 'BeamElement'. - xaxis : list[float], optional - The x-axis direction, by default [0,1,0]. + frame : :class:`compas.geometry.Line` or list[float], optional + Local frame of the element or x-axis of the frame by default [0,1,0]. name : str, optional The name of the part, by default None (one is automatically generated). @@ -688,33 +695,13 @@ def from_compas_lines_discretized(cls, lines, targetlength, element_class, secti The part. """ - prt = cls(name=name) - - for line in lines: - # initialization of first point of first element - p1 = Point(line.start.x, line.start.y, line.start.z) - prt.add_nodes(prt.find_closest_nodes_to_point(p1, number_of_nodes=1) or [Node(xyz=[p1.x, p1.y, p1.z])]) - - # initialization of second point of first element - n = int(line.length // targetlength) if line.length>targetlength else 1 - v_discretized = Vector(line.vector[0] / n, line.vector[1] / n, line.vector[2] / n) - p2 = p1.translated(v_discretized) - - # loop to create the elements except the last one - for i in range(n-1): - prt.add_node(Node(xyz=[p2.x, p2.y, p2.z])) - nodes = (list(prt.find_closest_nodes_to_point(p1, number_of_nodes=1)) or [Node(xyz=[p1.x, p1.y, p1.z])]) + (list(prt.find_closest_nodes_to_point(p2, number_of_nodes=1) or [Node(xyz=[p2.x, p2.y, p2.z])])) - element = element_class(nodes=nodes, section=section, frame=frame) - prt.add_element(element) - - p1.translate(v_discretized) - p2.translate(v_discretized) - - #the last element might not respect the targetlength and is therefore created indepedantly - prt.add_node(Node(xyz=[line.end.x, line.end.y, line.end.z])) - nodes = (list(prt.find_closest_nodes_to_point(p1, number_of_nodes=1)) or [Node(xyz=[p1.x, p1.y, p1.z])]) + list(prt.find_closest_nodes_to_point(line.end, number_of_nodes=1) ) - element = element_class(nodes=nodes, section=section, frame=frame) - prt.add_element(element) + polyline = Polyline(points = [line.start for line in lines] + [lines[-1].end]) + dividedline_points = polyline.divide_by_length(targetlength, strict=False) + prt = Part.from_compas_lines(lines=[Line(start=dividedline_points[i], end=dividedline_points[i+1]) for i in range(len(dividedline_points)-1)], + section=section, + element_model=element_model, + frame=frame, + name=name) return prt @@ -724,7 +711,7 @@ def from_compas_lines( cls, lines: List["compas.geometry.Line"], element_model: str = "BeamElement", - xaxis: List[float] = [0, 1, 0], + frame: Frame or List[float] = [0, 1, 0], section: Optional["_Section"] = None, name: Optional[str] = None, **kwargs, @@ -737,7 +724,7 @@ def from_compas_lines( The lines to be converted. element_model : str, optional Implementation model for the element, by default 'BeamElement'. - xaxis : list[float], optional + frame : :class:`compas.geometry.Line` or list[float], optional The x-axis direction, by default [0,1,0]. section : :class:`compas_fea2.model.BeamSection`, optional The section to be assigned to the elements, by default None. @@ -755,7 +742,8 @@ def from_compas_lines( prt = cls(name=name) mass = kwargs.get("mass", None) for line in lines: - frame = Frame(line.start, xaxis, line.vector) + if not(isinstance(frame, Frame)): + frame = Frame(line.start, frame, line.vector) nodes = [] for p in [line.start, line.end]: @@ -796,12 +784,13 @@ def shell_from_compas_mesh(cls, mesh, section: ShellSection, name: Optional[str] """ implementation = kwargs.get("implementation", None) ndm = kwargs.get("ndm", None) + heat = kwargs.get("heat", False) part = cls(name=name, ndm=ndm) if ndm else cls(name=name) vertex_node = {vertex: part.add_node(Node(mesh.vertex_coordinates(vertex))) for vertex in mesh.vertices()} for face in mesh.faces(): nodes = [vertex_node[vertex] for vertex in mesh.face_vertices(face)] - element = ShellElement(nodes=nodes, section=section, implementation=implementation) + element = ShellElement(nodes=nodes, section=section, implementation=implementation, heat=heat, **kwargs) part.add_element(element) part._boundary_mesh = mesh @@ -866,6 +855,7 @@ def from_gmsh(cls, gmshModel, section: Optional[Union[SolidSection, ShellSection verbose = kwargs.get("verbose", False) rigid = kwargs.get("rigid", False) implementation = kwargs.get("implementation", None) + heat = kwargs.get("heat", False) for ntags in ntags_per_element: if kwargs.get("split", False): @@ -880,6 +870,8 @@ def from_gmsh(cls, gmshModel, section: Optional[Union[SolidSection, ShellSection section=section, rigid=rigid, implementation=implementation, + heat=heat + **kwargs ) ) @@ -891,18 +883,20 @@ def from_gmsh(cls, gmshModel, section: Optional[Union[SolidSection, ShellSection section=section, rigid=rigid, implementation=implementation, + heat=heat + **kwargs ) ) else: - part.add_element(TetrahedronElement(nodes=element_nodes, section=section, rigid=rigid)) + part.add_element(TetrahedronElement(nodes=element_nodes, section=section, rigid=rigid, heat=heat)) part.ndf = 3 # FIXME: move outside the loop elif ntags.size == 10: # C3D10 tetrahedral element - part.add_element(TetrahedronElement(nodes=element_nodes, section=section, rigid=rigid)) + part.add_element(TetrahedronElement(nodes=element_nodes, section=section, rigid=rigid, heat=heat)) part.ndf = 3 elif ntags.size == 8: - part.add_element(HexahedronElement(nodes=element_nodes, section=section, rigid=rigid)) + part.add_element(HexahedronElement(nodes=element_nodes, section=section, rigid=rigid, heat=heat)) else: raise NotImplementedError(f"Element with {ntags.size} nodes not supported") diff --git a/src/compas_fea2/problem/__init__.py b/src/compas_fea2/problem/__init__.py index 6cb54fa54..0294149f0 100644 --- a/src/compas_fea2/problem/__init__.py +++ b/src/compas_fea2/problem/__init__.py @@ -1,7 +1,7 @@ from .problem import Problem from .displacements import GeneralDisplacement from .loads import ( - Load, + VectorLoad, PrestressLoad, ConcentratedLoad, PressureLoad, @@ -10,6 +10,8 @@ HarmonicPointLoad, HarmonicPressureLoad, ThermalLoad, + HeatFluxLoad, + TemperatureLoad ) from .fields import ( @@ -19,9 +21,17 @@ PointLoadField, _PrescribedField, PrescribedTemperatureField, + HeatFluxField, + ConvectionField, + RadiationField, + TemperatureField, + SurfaceLoadField ) + from .combinations import LoadCombination +from.amplitudes import Amplitude + from .steps import ( Step, GeneralStep, @@ -34,13 +44,14 @@ DynamicStep, QuasiStaticStep, DirectCyclicStep, + HeatTransferStep ) __all__ = [ "Problem", "GeneralDisplacement", - "Load", + "VectorLoad", "PrestressLoad", "ConcentratedLoad", "PressureLoad", diff --git a/src/compas_fea2/problem/amplitudes.py b/src/compas_fea2/problem/amplitudes.py new file mode 100644 index 000000000..8a51906fa --- /dev/null +++ b/src/compas_fea2/problem/amplitudes.py @@ -0,0 +1,48 @@ +from compas_fea2.base import FEAData + +class Amplitude(FEAData): + """Amplitude object for creating a time-varying fonction that defines + how the magnitude of a load, boundary condition, or predifined field + changes throughout a step or the analysis. + + Parameters + ---------- + multipliers : list[float] + Multipliers list of values. + + times : list[float] + Corresponding time to the multipliers. + + Attributes + ---------- + q : float + Heat flux value of the load. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. """ + + def __init__(self, multipliers, times, name=None, **kwargs): + super().__init__(name, **kwargs) + self._multipliers = multipliers + self._times = times + if len(self._multipliers)!=len(self._times): + raise ValueError("The lists of values and times must have the same length.") + + @property + def multipliers(self): + return self._multipliers + + @property + def times(self): + return self._times + + @property + def multipliers_times(self): + """Return a list of tuples with the values and the assigned time.""" + return zip(self.multipliers, self.times) + + + @classmethod + def equally_spaced(cls, multipliers, first_value, fixed_interval, **kwargs): + times = [first_value+fixed_interval*i for i in range(len(fixed_interval))] + return cls(multipliers=multipliers, times=times, **kwargs) diff --git a/src/compas_fea2/problem/combinations.py b/src/compas_fea2/problem/combinations.py index cd41cf8b9..dff0b3bc9 100644 --- a/src/compas_fea2/problem/combinations.py +++ b/src/compas_fea2/problem/combinations.py @@ -4,8 +4,8 @@ from compas_fea2.base import FEAData if TYPE_CHECKING: - from compas_fea2.problem import Step from compas_fea2.problem import Problem + from compas_fea2.problem import Step diff --git a/src/compas_fea2/problem/fields.py b/src/compas_fea2/problem/fields.py index dbadea21d..1dc98aa22 100644 --- a/src/compas_fea2/problem/fields.py +++ b/src/compas_fea2/problem/fields.py @@ -2,6 +2,7 @@ from compas_fea2.base import FEAData from compas_fea2.problem.loads import GravityLoad +from compas_fea2.problem.loads import HeatFluxLoad # TODO implement __*__ magic method for combination @@ -53,10 +54,19 @@ def __init__( self.load_case = load_case self._registration = None + self._amplitudes = set() + for load in self._loads: + if load.amplitude: + self._amplitudes.add(load.amplitude) + @property def loads(self): return self._loads + @property + def amplitudes(self): + return self._amplitudes + @property def distribution(self): return self._distribution @@ -66,7 +76,7 @@ def step(self): if self._registration: return self._registration else: - raise ValueError("Register the Pattern to a Step first.") + raise ValueError("Register the LoadField to a Step first.") @property def problem(self): @@ -220,9 +230,9 @@ def __init__(self, **kwargs): class PrescribedTemperatureField(_PrescribedField): - """Temperature field""" + """Temperature field.""" - def __init__(self, temperature, **kwargs): + def __init__(self, temperature=None, **kwargs): super(PrescribedTemperatureField, self).__init__(**kwargs) self._t = temperature @@ -233,3 +243,100 @@ def temperature(self): @temperature.setter def temperature(self, value): self._t = value + + @classmethod + def from_file(cls, path, **kwargs): + cls._path = path + return cls(path=path, **kwargs) + +# ===================================================================== +# HEAT ANALYSIS +# ===================================================================== + +class TemperatureField(NodeLoadField): + def __init__(self, temperature, nodes, **kwargs): + nodes = nodes if isinstance(nodes, Iterable) else [nodes] + loads = temperature if isinstance(temperature, Iterable) else [temperature] * len(nodes) + super().__init__(loads=loads, nodes=nodes,**kwargs) + + +class SurfaceLoadField(LoadField): + def __init__(self, loads, surface, **kwargs): + super().__init__(loads=loads, distribution=surface, **kwargs) + + @property + def surface(self): + return self.distribution + + @property + def load_surface(self): + return zip(self.loads, self.surface) + + +class HeatFluxField(SurfaceLoadField): + """A distribution of a uniform convection flux on a set of face elements. + A non-uniform distribution is not yet implemented. + + Parameters + ---------- + heatflux : object + The heatflux to be applied. + face_elements : list + List of elements face where the surface heat flux is applied. + load_case : object, optional + The load case to which this pattern belongs. + """ + + def __init__(self, heatflux: HeatFluxLoad, surface, **kwargs): + surface = surface if isinstance(surface, Iterable) else [surface] + heatflux = heatflux if isinstance(heatflux, Iterable) else [heatflux] * len(surface) + super().__init__(loads=heatflux, surface=surface,**kwargs) + + @property + def heatflux(self): + return self.loads + +class ConvectionField(SurfaceLoadField): + """ + """ + + def __init__(self, h, temperature, surface, **kwargs): + surface = surface if isinstance(surface, Iterable) else [surface] + temperature = temperature if isinstance(temperature, Iterable) else [temperature] * len(surface) + super().__init__(loads=temperature, surface=surface,**kwargs) + self._h = h + + @property + def surface(self): + return self._distribution + + @property + def temperature(self): + return self._loads + + @property + def h(self): + return self._h + + +class RadiationField(SurfaceLoadField): + """ + """ + + def __init__(self, eps, temperature, surface, **kwargs): + surface = surface if isinstance(surface, Iterable) else [surface] + temperature = temperature if isinstance(temperature, Iterable) else [temperature] * len(surface) + super().__init__(loads=temperature, surface=surface,**kwargs) + self._eps = eps + + @property + def surface(self): + return self._distribution + + @property + def temperature(self): + return self._loads + + @property + def eps(self): + return self._eps diff --git a/src/compas_fea2/problem/loads.py b/src/compas_fea2/problem/loads.py index bbfdaf732..6fb3785a1 100644 --- a/src/compas_fea2/problem/loads.py +++ b/src/compas_fea2/problem/loads.py @@ -2,38 +2,154 @@ # TODO: make units independent using the utilities function - -class Load(FEAData): - """Initialises base Load object. +class _Load(FEAData): + """Initialises base _Load object. Parameters ---------- - name : str - Uniqe identifier. If not provided it is automatically generated. Set a - name if you want a more human-readable input file. - components : dict - Load components. - axes : str, optional - Load applied via 'local' or 'global' axes, by default 'global'. + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. Attributes ---------- + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + + field : :class:`compas_fea2.problem.LoadField` + Field associated with the load. + + step : :class:`compas_fea2.problem.Step` + Step associated with the load. + + problem : :class:`compas_fea2.problem.Problem` + Problem associated with the load. + + model : :class:`compas_fea2.model.Model` + Model associated with the load. + name : str Uniqe identifier. If not provided it is automatically generated. Set a name if you want a more human-readable input file. - components : dict - Load components. These differ according to each Load type - axes : str, optional - Load applied via 'local' or 'global' axes, by default 'global'. + """ + def __init__(self, amplitude = None, **kwargs): + super().__init__(**kwargs) + self._amplitude = amplitude + + @property + def amplitude(self): + return self._amplitude - Notes - ----- - Loads are registered to a :class:`compas_fea2.problem.Pattern`. + + @property + def field(self): + return self._registration + + @property + def step(self): + return self.field._registration + + @property + def problem(self): + return self.step._registration + + @property + def model(self): + return self.problem._registration + +class ScalarLoad(_Load): + """Scalar load object. + + Parameters + ---------- + scalar_load : float + Scalar value of the load. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + + Attributes + ---------- + scalar_load : float + Scalar value of the load + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + """ + def __init__(self, scalar_load, amplitude = None, **kwargs): + super().__init__(amplitude=amplitude, **kwargs) + if not(isinstance(scalar_load, (int,float))) : + raise ValueError("The scalar_load must be a float.") + self._scalar_load = scalar_load + + @property + def scalar_load(self): + return self._scalar_load + + +class VectorLoad(_Load): + """Vector load object. + + Parameters + ---------- + axes : str, "local" or "global" + The load is either defined in the local frame or the global one. + If not indicated, the global frame is considered. + + x : float + x-axis force value of the load. + + y : float + y-axis force value of the load. + + z : float + z-axis force value of the load. + + xx : float + Moment value of the load about the x-axis. + + yy : float + Moment value of the load about the y-axis. + + zz : float + Moment value of the load about the z-axis. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + + Attributes + ---------- + axes : str, "local" or "global" + The load is either defined in the local frame or the global one. + If not indicated, the global frame is considered. + + x : float + x-axis force value of the load. + + y : float + y-axis force value of the load. + + z : float + z-axis force value of the load. + + xx : float + Moment value of the load about the x-axis. + + yy : float + Moment value of the load about the y-axis. + + zz : float + Moment value of the load about the z-axis. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + + components : {str: float} + Dictionnary of the components of the load and values """ def __init__(self, x=None, y=None, z=None, xx=None, yy=None, zz=None, axes="global", **kwargs): - super(Load, self).__init__(**kwargs) + super(VectorLoad, self).__init__(**kwargs) self.axes = axes self.x = x self.y = y @@ -51,24 +167,7 @@ def components(self, value): for k, v in value: setattr(self, k, v) - @property - def pattern(self): - return self._registration - - @property - def step(self): - return self.pattern._registration - - @property - def problem(self): - return self.step._registration - - @property - def model(self): - return self.problem._registration - - -class ConcentratedLoad(Load): +class ConcentratedLoad(VectorLoad): """Concentrated forces and moments [units:N, Nm]. Parameters @@ -133,7 +232,7 @@ def __radd__(self, other): return self.__add__(other) -class PressureLoad(Load): +class PressureLoad(VectorLoad): """Distributed area force [e.g. units:N/m2] applied to element(s). Parameters @@ -167,7 +266,7 @@ def __init__(self, x=0, y=0, z=0, axes="local", **kwargs): raise NotImplementedError -class GravityLoad(Load): +class GravityLoad(VectorLoad): """Gravity load [units:N/m3] applied to element(s). Parameters @@ -234,7 +333,7 @@ def __rmul__(self, other): return self.__mul__(other) -class PrestressLoad(Load): +class PrestressLoad(VectorLoad): """Prestress load""" def __init__(self, components, axes="global", **kwargs): @@ -242,14 +341,14 @@ def __init__(self, components, axes="global", **kwargs): raise NotImplementedError -class ThermalLoad(Load): +class ThermalLoad(VectorLoad): """Thermal load""" def __init__(self, components, axes="global", **kwargs): super(ThermalLoad, self).__init__(components, axes, **kwargs) -class TributaryLoad(Load): +class TributaryLoad(VectorLoad): """Tributary load""" def __init__(self, components, axes="global", **kwargs): @@ -257,7 +356,7 @@ def __init__(self, components, axes="global", **kwargs): raise NotImplementedError -class HarmonicPointLoad(Load): +class HarmonicPointLoad(VectorLoad): """""" def __init__(self, components, axes="global", **kwargs): @@ -265,9 +364,68 @@ def __init__(self, components, axes="global", **kwargs): raise NotImplementedError -class HarmonicPressureLoad(Load): +class HarmonicPressureLoad(VectorLoad): """""" def __init__(self, components, axes="global", **kwargs): super(HarmonicPressureLoad, self).__init__(components, axes, **kwargs) raise NotImplementedError + + +#=================================================================== +# HEAT ANALYSIS +#=================================================================== + +class HeatFluxLoad(ScalarLoad): + """ Heat flux load for heat analysis. + + Parameters + ---------- + q : float + Heat flux value of the load. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + + Attributes + ---------- + q : float + Heat flux value of the load. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + """ + + def __init__(self, q, amplitude = None, **kwargs): + super().__init__(scalar_load=q, amplitude=amplitude, **kwargs) + + @property + def q(self): + return self._scalar_load + +class TemperatureLoad(ScalarLoad): + """ Temperature load for heat analysis + + Parameters + ---------- + temperature : float + Value of the temperature load. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + + Attributes + ---------- + temperature : float + Value of the temperature load. + + amplitude : :class:`compas_fea2.problem.Amplitude` + Amplitude associated to the load, optionnal. + """ + + def __init__(self, temperature, amplitude = None, **kwargs): + super().__init__(scalar_load=temperature, amplitude=amplitude, **kwargs) + + @property + def temperature(self): + return self._scalar_load diff --git a/src/compas_fea2/problem/problem.py b/src/compas_fea2/problem/problem.py index 819a6cc75..677aaf081 100644 --- a/src/compas_fea2/problem/problem.py +++ b/src/compas_fea2/problem/problem.py @@ -11,7 +11,8 @@ from compas_fea2.job.input_file import InputFile from compas_fea2.problem.steps import StaticStep from compas_fea2.problem.steps import Step -from compas_fea2.results.database import ResultsDatabase, SQLiteResultsDatabase +from compas_fea2.results.database import ResultsDatabase +from compas_fea2.results.database import SQLiteResultsDatabase class Problem(FEAData): diff --git a/src/compas_fea2/problem/steps/__init__.py b/src/compas_fea2/problem/steps/__init__.py index ee3afb841..41be43d8a 100644 --- a/src/compas_fea2/problem/steps/__init__.py +++ b/src/compas_fea2/problem/steps/__init__.py @@ -17,6 +17,8 @@ DirectCyclicStep, ) +from .heattransfer import HeatTransferStep + from .perturbations import ( _Perturbation, ModalAnalysis, @@ -42,4 +44,5 @@ "DynamicStep", "QuasiStaticStep", "DirectCyclicStep", + "HeatTransferStep" ] diff --git a/src/compas_fea2/problem/steps/heattransfer.py b/src/compas_fea2/problem/steps/heattransfer.py new file mode 100644 index 000000000..90d720df6 --- /dev/null +++ b/src/compas_fea2/problem/steps/heattransfer.py @@ -0,0 +1,160 @@ +from compas_fea2.problem.fields import HeatFluxField +from compas_fea2.problem.fields import TemperatureField + +from .step import Step + + +class HeatTransferStep(Step): + """HeatTransfer for use in a heat transfer analysis. + Specific for now to a transient heat transfer analysis as defined in Abaqus. + + Parameters + ---------- + max_increments : int + Max number of increments to perform during the case step. + (Typically 100 but you might have to increase it in highly non-linear + problems. This might increase the analysis time.). + initial_inc_size : float + Sets the the size of the increment for the first iteration. + (By default is equal to the total time, meaning that the software decrease + the size automatically.) + min_inc_size : float + Minimum increment size before stopping the analysis. + (By default is 1e-5, but you can set a smaller size for highly non-linear + problems. This might increase the analysis time.) + max_temp_delta : float + Maximum allowable temperature change per increment. + (By default is 10.) + max_emiss_delta : float + Maximum allowable emissivity change per increment. + time : float + Total time of the case step. Note that this not actual 'time', + but rather a proportionality factor. (By default is 1, meaning that the + analysis is complete when all the increments sum up to 1) + nlgeom : bool + if ``True`` nonlinear geometry effects are considered. + modify : bool + if ``True`` the loads applied in a previous step are substituted by the + ones defined in the present step, otherwise the loads are added. + + Attributes + ---------- + name : str + Automatically generated id. You can change the name if you want a more + human readable input file. + max_increments : int + Max number of increments to perform during the case step. + (Typically 100 but you might have to increase it in highly non-linear + problems. This might increase the analysis time.). + initial_inc_size : float + Sets the the size of the increment for the first iteration. + (By default is equal to the total time, meaning that the software decrease + the size automatically.) + min_inc_size : float + Minimum increment size before stopping the analysis. + (By default is 1e-5, but you can set a smaller size for highly non-linear + problems. This might increase the analysis time.) + max_temp_delta : float + Maximum allowable temperature change per increment. + (By default is 10.) + max_emiss_delta : float + Maximum allowable emissivity change per increment. + time : float + Total time of the case step. Note that this not actual 'time', + but rather a proportionality factor. (By default is 1, meaning that the + analysis is complete when all the increments sum up to 1) + nlgeom : bool + if ``True`` nonlinear geometry effects are considered. + modify : bool + if ``True`` the loads applied in a previous step are substituted by the + ones defined in the present step, otherwise the loads are added. + loads : dict + Dictionary of the loads assigned to each part in the model in the step. + displacements : dict + Dictionary of the displacements assigned to each part in the model in the step. + + """ + + def __init__( + self, + max_increments=100, + initial_inc_size=1, + min_inc_size=0.00001, + max_inc_size=1, + time=1, + max_temp_delta=10, + max_emiss_change=0.1, + nlgeom=False, + modify=True, + **kwargs, + ): + super().__init__( + **kwargs, + ) + + self._max_increments = max_increments + self._initial_inc_size = initial_inc_size + self._min_inc_size = min_inc_size + self._max_inc_size = max_inc_size + self._time = time + self.max_temp_delta = max_temp_delta + self.max_emiss_change = max_emiss_change + self._nlgeom = nlgeom + self._modify = modify + + + def __data__(self): + return { + "max_increments": self.max_increments, + "initial_inc_size": self.initial_inc_size, + "min_inc_size": self.min_inc_size, + "time": self.time, + "max_temp_delta": self.max_temp_delta, + "max_emiss_change": self.max_emiss_change, + "nlgeom": self.nlgeom, + "modify": self.modify, + # Add other attributes as needed + } + + @classmethod + def __from_data__(cls, data): + return cls( + max_increments=data["max_increments"], + initial_inc_size=data["initial_inc_size"], + min_inc_size=data["min_inc_size"], + max_temp_delta=data["max_temp_delta"], + max_emiss_change=data["max_emiss_change"], + time=data["time"], + nlgeom=data["nlgeom"], + modify=data["modify"], + # Add other attributes as needed + ) + + def add_load_field(self, field): + """Add a general :class:`compas_fea2.problem.patterns.Pattern` to the Step. + + Parameters + ---------- + load_pattern : :class:`compas_fea2.problem.patterns.Pattern` + The load pattern to add. + + Returns + ------- + :class:`compas_fea2.problem.patterns.Pattern` + + """ + from compas_fea2.problem.fields import LoadField + + if not isinstance(field, LoadField): + raise TypeError("{!r} is not a LoadPattern.".format(field)) + + if not(isinstance(field, (HeatFluxField, TemperatureField))): + raise ValueError("A non-thermal load can not be implemented in a heat analysis step.") + + self._load_fields.add(field) + self._load_cases.add(field.load_case) + field._registration = self + self.model.add_group(field.distribution) + for amplitude in field.amplitudes : + self.model.amplitudes.add(amplitude) + return field \ No newline at end of file diff --git a/src/compas_fea2/problem/steps/step.py b/src/compas_fea2/problem/steps/step.py index 91b084cfc..c4fc9c796 100644 --- a/src/compas_fea2/problem/steps/step.py +++ b/src/compas_fea2/problem/steps/step.py @@ -2,15 +2,23 @@ from compas.geometry import Point from compas.geometry import Vector +from compas.geometry import Polyline from compas.geometry import centroid_points_weighted from compas.geometry import sum_vectors from compas_fea2.base import FEAData from compas_fea2.UI import FEA2Viewer +from compas_fea2.model.nodes import Node +from compas_fea2.model.elements import _Element +from compas_fea2.model.groups import NodesGroup +from compas_fea2.model.nodes import Node from compas_fea2.problem.displacements import GeneralDisplacement from compas_fea2.problem.fields import DisplacementField from compas_fea2.problem.fields import NodeLoadField from compas_fea2.problem.fields import PointLoadField +from compas_fea2.problem.fields import _PrescribedField +from compas_fea2.results import TemperatureFieldResults +from compas_fea2.UI import FEA2Viewer # ============================================================================== # Base Steps @@ -61,6 +69,7 @@ def __init__(self, replace=False, **kwargs): self._load_fields = set() self._load_cases = set() self._combination = None + self._prescribed_fields = dict() @property def problem(self): @@ -85,6 +94,10 @@ def load_fields(self): @property def combination(self): return self._combination + + @property + def prescribed_fields(self): + return self._prescribed_fields @combination.setter def combination(self, combination): @@ -187,7 +200,7 @@ def reaction_field(self): @property def temperature_field(self): - raise NotImplementedError + return TemperatureFieldResults(self) @property def stress_field(self): @@ -649,6 +662,81 @@ def add_uniform_displacement_field(self, nodes, x=None, y=None, z=None, xx=None, displacement = GeneralDisplacement(x=x, y=y, z=z, xx=xx, yy=yy, zz=zz, axes=axes, **kwargs) return self.add_load_field(DisplacementField(displacement, nodes)) + # ============================================================================== + # Prescribed Fields + # ============================================================================== + def _add_prescribed_field(self, prescribed_field, group): + """Add a :class:`compas_fea2.model._InitialCondition` to the model. + + Parameters + ---------- + ic : :class:`compas_fea2.model._InitialCondition` + Initial condition object to add to the model. + group : :class:`compas_fea2.model._Group` + Group of Nodes/Elements where the initial condition is assigned. + + Returns + ------- + :class:`compas_fea2.model._InitialCondition` + + """ + # group.part.add_group(group) + + if not isinstance(prescribed_field, _PrescribedField): + raise TypeError("{!r} is not a InitialCondition.".format(ic)) + for member in group.members: + if not isinstance(member, (Node, _Element)): + raise TypeError("{!r} is not a Node or an Element.".format(member)) + if not member.part: + raise ValueError("{!r} is not registered to any part.".format(member)) + elif member.part not in self.model.parts: + raise ValueError("{!r} belongs to a part not registered to this model.".format(member)) + member._prescribed_field = prescribed_field + + prescribed_field._key = len(self._prescribed_fields) + self._prescribed_fields[prescribed_field] = group.members + prescribed_field._registration = self + + return prescribed_field + + def add_nodes_prescribed_field(self, prescribed_field, nodes): + """Add a :class:`compas_fea2.model._InitialCondition` to the model. + + Parameters + ---------- + ic : :class:`compas_fea2.model._InitialCondition` + nodes : list[:class:`compas_fea2.model.Node`] or :class:`compas_fea2.model.NodesGroup` + + Returns + ------- + :class:`compas_fea2.model._InitialCondition` + + """ + if not isinstance(nodes, NodesGroup): + raise TypeError("{} is not a group of nodes".format(nodes)) + self._add_prescribed_field(prescribed_field, nodes) + return prescribed_field + + def add_prescribed_field(self, field, *kwargs): + """Add a general :class:`compas_fea2.problem.patterns.Pattern` to the Step. + + Parameters + ---------- + + Returns + ------- + :class:`compas_fea2.problem.patterns.Pattern` + + """ + + if not isinstance(field, _PrescribedField): + raise TypeError("{!r} is not a _PrescribedField.".format(field)) + + self._prescribed_fields.add(field) + field._registration = self + return field + + # ============================================================================== # Combinations # ============================================================================== @@ -771,7 +859,6 @@ def show_deformed(self, opacity=1, show_bcs=1, scale_results=1, scale_model=1, s None """ - from compas_fea2.UI import FEA2Viewer viewer = FEA2Viewer(center=self.model.center, scale_model=scale_model) @@ -802,7 +889,6 @@ def show_displacements(self, fast=True, show_bcs=1, scale_model=1, show_loads=0. _description_, by default """ - from compas_fea2.UI import FEA2Viewer if not self.displacement_field: raise ValueError("No displacement field results available for this step") @@ -833,7 +919,6 @@ def show_reactions(self, fast=True, show_bcs=1, scale_model=1, show_loads=0.1, c scale_results : _type_, optional _description_, by default 1 """ - from compas_fea2.UI import FEA2Viewer if not self.reaction_field: raise ValueError("No reaction field results available for this step") @@ -848,7 +933,6 @@ def show_reactions(self, fast=True, show_bcs=1, scale_model=1, show_loads=0.1, c viewer.scene.clear() def show_stress(self, fast=True, show_bcs=1, scale_model=1, show_loads=0.1, component=None, show_vectors=1, show_contour=False, plane="mid", **kwargs): - from compas_fea2.UI import FEA2Viewer if not self.stress_field: raise ValueError("No reaction field results available for this step") @@ -878,7 +962,7 @@ def __data__(self): ) return data - def plot_deflection_along_line(self, line, step=None, n_divide=1000, plot_data=False): + def plot_deflection_along_line(self, line, n_divide=1000): """Plot the deflection along a compas line given as an input. This method can only be used on shell models. @@ -904,12 +988,9 @@ def plot_deflection_along_line(self, line, step=None, n_divide=1000, plot_data= #FIRST, the input line is discretized in n_divide points #----------------------------------------------------------- # TODO automatized the n_divide parameters with the mesh density - v_discretized=Vector(line.vector[0]/n_divide, line.vector[1]/n_divide, line.vector[2]/n_divide) - l_discretized=[line.start] - for i in range(n_divide-1): - l_discretized.append(l_discretized[i].translated(v_discretized)) - l_discretized.append(line.end) + length = line.length/n_divide + l_discretized=Polyline(points=[line.start, line.end]).divide_by_length(length) #-------------------------------------------------------------------------------------------------------------------- #SECOND, looking for the closest points of mesh to input line, according to their projection on the horizontal plan @@ -922,7 +1003,7 @@ def plot_deflection_along_line(self, line, step=None, n_divide=1000, plot_data= for node in nodes : element_XY_points.append(Point(node.xyz[0], node.xyz[1], 0)) #nodes of the shell are projected vertically - #determination of the closest nodes + #determination of the closest nodes of the projected mesh to the input line l_closestpoints=[] for point_line in l_discretized: tree=KDTree(element_XY_points) @@ -930,7 +1011,8 @@ def plot_deflection_along_line(self, line, step=None, n_divide=1000, plot_data= closest_3D_point=list(nodes)[closest_2D_point_index] l_closestpoints.append(closest_3D_point) - #Elimination of the points that have the same "closest node" + #The discretization of the input line might be more precised than the precision of the mesh + #The points of the line associated to the same mesh node are removed i=0 while i