Beam release and other modifications

src/compas_fea2/model/constraints.py

@@ -45,6 +45,21 @@ def registration(self, value: "Model") -> None:
     def __from_data__(cls, data, registry: Optional[Registry] = None, duplicate=True):
         raise NotImplementedError("LinearConnector does not support from_data method yet.")
 
+# ------------------------------------------------------------------------------
+# Coupling
+# ------------------------------------------------------------------------------
+
+class CouplingConstraint(_Constraint):
+    def __init__(self, x=True, y=True, z=True, xx=True, yy=True, zz=True, **kwargs):
+        super().__init__(**kwargs)
+
+        self.directions = ['x', 'y', 'z', 'xx', 'yy', 'zz']
+        self.x = x
+        self.y = y
+        self.z = z
+        self.xx = xx
+        self.yy = yy
+        self.zz = zz
 
 # ------------------------------------------------------------------------------
 # MPC
@@ -166,4 +181,4 @@ class TieConstraint(_SurfaceConstraint):
     Notes
     -----
     The tolerance (if used) is interpreted in the active unit system ("length").
-    """
+    """

src/compas_fea2/model/fields.py

@@ -9,10 +9,14 @@
 from compas_fea2.model.groups import NodesGroup
 from compas_fea2.model.ics import _InitialCondition
 from compas_fea2.model.nodes import Node
+from compas_fea2.model.groups import ElementsGroup
+from compas_fea2.model.elements import _Element1D
+from compas_fea2.model.releases import _BeamEndRelease
 
 if TYPE_CHECKING:
     from compas_fea2.model.bcs import _BoundaryCondition
     from compas_fea2.model.ics import _InitialCondition
+    from compas_fea2.model.parts import _Part
     from compas_fea2.model.model import Model
 
 
@@ -172,3 +176,90 @@ class InitialStressField(_InitialConditionField):
     def __init__(self, elements, condition, **kwargs):
         super().__init__(distribution=elements, conditions=condition, **kwargs)
         raise NotImplementedError("InitialStressField is not yet implemented.")
+
+
+# ------------------------------------------------------------------------------
+# BEAM END RELEASE FIELDS
+# ------------------------------------------------------------------------------
+
+class BeamReleaseField(FEAData):
+    """A BeamReleaseField is the spatial distribution of a specific set of beam end releases.
+
+    Parameters
+    ----------
+    release : `:class:compas_fea2.model._BeamEndRelease`
+        The release applied to the field.
+    elements : `:class:compas_fea2.model.ElementsGroup`
+        The elements composing the fields.
+    release_end : str 
+        The end to release is applied (start, end or both).
+    name : str, optional
+        Unique identifier for the field.
+
+    Attributes
+    ----------
+    release : `:class:compas_fea2.model._BeamEndRelease`
+        The release applied to the field.
+    elements : `:class:compas_fea2.model.ElementsGroup`
+        The elements composing the fields.
+    release_end : str
+        The end to release is applied (start, end, both)
+    name : str, optional
+        Unique identifier for the field.
+        Registered model to the field. None if the field has not been registered.
+    part :  `:class:compas_fea2.model._Part`
+        Part of the release field
+    """
+    def __init__(self, release, elements, end, name = None, **kwargs):
+        super().__init__(name, **kwargs)
+        if not isinstance(elements, ElementsGroup):
+            self._elements = ElementsGroup(elements)
+        else:
+            self._elements = elements
+        for element in self._elements:
+            if not isinstance(element, _Element1D):
+                raise ValueError('A BeamEndRelease can only applied on 1D elements.')
+        self._release = release
+
+
+        self._end = end
+        self._check_registration()
+
+    def _check_registration(self):
+        registrations = set([n.model for n in self._elements])
+        if len(registrations) != 1:
+            raise ValueError("All nodes in the distribution must be registered to the same model.")
+        if self._registration is None:
+            self._registration = registrations.pop()
+
+    @property
+    def release(self) -> "_BeamEndRelease":
+        if isinstance(self._release, _BeamEndRelease):
+            return self._release
+        else :
+            raise ValueError("Release is not a _BeamEndRelease.")
+
+    @property
+    def elements(self) -> "ElementsGroup":
+        return self._elements
+
+    @property
+    def end(self) -> str :
+        return self._end
+
+    @property
+    def elements_end_release(self):
+        """Return a list of tuples with the nodes and their assigned condition."""
+        return zip(self.elements, [self.end] * len(self.elements), [self.release] * len(self.elements))
+
+    @property
+    def part(self) -> "_Part":
+        if self._registration:
+            return self._registration
+        else:
+            raise ValueError("Register the ConditionField to a part first.")
+
+    @property
+    def model(self) -> "Model":
+        return self.part.model
+

src/compas_fea2/model/materials/steel.py

@@ -3,6 +3,7 @@
 from compas_fea2.base import Registry
 from compas_fea2.base import from_data
 from compas_fea2.units import units_io
+import compas_fea2.units as u
 
 from .material import ElasticIsotropic
 
@@ -51,7 +52,7 @@ class Steel(ElasticIsotropic):
 
     @units_io(types_in=("stress", None, "density", "stress", "stress", None), types_out=None)
     def __init__(self, *, E, v, density, fy, fu, eu, **kwargs):
-        super().__init__(E=E, v=v, density=density, **kwargs)
+        super().__init__(E=E*u.MPa, v=v, density=density*u.kg_per_m3, **kwargs)
 
         fu = fu or fy
 
@@ -70,8 +71,6 @@ def __init__(self, *, E, v, density, fy, fu, eu, **kwargs):
         self.fu = fu
         self.eu = eu
         self.ep = ep
-        self.E = E
-        self.v = v
         self.tension = {"f": f, "e": e}
         self.compression = {"f": fc, "e": ec}
 
@@ -139,4 +138,4 @@ def S355(cls):
         :class:`compas_fea2.model.material.Steel`
             The precompiled steel material.
         """
-        return cls(fy=355, fu=None, eu=20, E=210, v=0.3, density=7850, name=None)
+        return cls(fy=355, fu=None, eu=20, E=210000, v=0.3, density=7850, name=None)

src/compas_fea2/model/model.py

@@ -927,6 +927,9 @@ def _add_bc_type(self, bc_type: str, nodes: "Union[list[Node], NodesGroup]", fra
             "rollerY": "RollerBCY",
             "rollerZ": "RollerBCZ",
             "thermal": "ThermalBC",
+            "rollerXY": "RollerBCXY",
+            "rollerXZ": "RollerBCXZ",
+            "rollerYZ": "RollerBCYZ"
         }
         m = importlib.import_module("compas_fea2.model.bcs")
         bc = getattr(m, types[bc_type])(frame=frame)

src/compas_fea2/model/parts.py

@@ -53,9 +53,11 @@
 from .groups import MaterialsGroup
 from .groups import NodesGroup
 from .groups import SectionsGroup
+from .groups import FieldsGroup
 from .materials.material import _Material
 from .nodes import Node
 from .releases import _BeamEndRelease
+from .fields import BeamReleaseField
 from .sections import ShellSection
 from .sections import SolidSection
 from .sections import _Section
@@ -143,6 +145,8 @@ def __init__(self, **kwargs):
 
         self._reference_node: Optional[Node] = None
 
+        self._releases_fields = FieldsGroup(members=[], name=self.name+"_RELEASES_FIELDS")
+
     @property
     def __data__(self):
         data = super().__data__
@@ -207,7 +211,7 @@ def _build_graph(self):
     # =========================================================================
     @classmethod
     @units_io(types_in=(None, "length", None, None, None), types_out=None)
-    def from_compas_lines_discretized(cls, lines: List["Line"], targetlength: float, element_cls: type, section: "_Section2D", frame: "Union[Frame, List[float], Vector]", **kwargs):
+    def from_compas_lines_discretized(cls, lines: List["Line"], targetlength: float, element_cls: type, section: "_Section2D", orientation: "Union[Frame, List[float], Vector]", **kwargs):
         """Generate a discretized model from a list of :class:`compas.geometry.Line`.
 
         Parameters
@@ -221,7 +225,7 @@ def from_compas_lines_discretized(cls, lines: List["Line"], targetlength: float,
             The section to be assigned to the elements, by default None.
         element_model : str, optional
             Implementation model for the element, by default 'BeamElement'.
-        frame : :class:`compas.geometry.Vector` or list[float], optional
+        orientation : :class:`compas.geometry.Point` or list[float]
             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).
@@ -238,7 +242,7 @@ def from_compas_lines_discretized(cls, lines: List["Line"], targetlength: float,
             lines=[Line(start=dividedline_points[i], end=dividedline_points[i + 1]) for i in range(len(dividedline_points) - 1)],
             section=section,
             element_cls=element_cls,
-            frame=frame,
+            orientation=orientation,
             **kwargs,
         )
 
@@ -249,7 +253,7 @@ def from_compas_lines(
         cls,
         lines: List["Line"],
         element_cls: type = BeamElement,
-        frame: "Union[Frame, List[float], Vector]" = [0, 1, 0],
+        orientation: "Union[Point, List[float]]" = [0, 1, 0],
         section: Optional["_Section"] = None,
         name: Optional[str] = None,
         **kwargs,
@@ -262,8 +266,8 @@ def from_compas_lines(
             The lines to be converted.
         element_cls : type, optional
             Implementation model for the element, by default BeamElement.
-        frame : :class:`compas.geometry.Line` or list[float], optional
-            The x-axis direction, by default [0,1,0].
+        orientation : :class:`compas.geometry.Point` or list[float], optional
+            Point determining the x-axis direction, by default [0,1,0].
         section : :class:`compas_fea2.model.Section1D`, optional
             The section to be assigned to the elements, by default None.
         name : str, optional
@@ -278,17 +282,16 @@ def from_compas_lines(
         prt = cls(name=name)
         mass = kwargs.get("mass", None)
         for line in lines:
-            if not (isinstance(frame, Frame)):
-                frame = Frame(line.start, frame, line.vector)
             nodes = []
             for p in [line.start, line.end]:
-                if g := prt.nodes.subgroup(condition=lambda node: node.point == p):
-                    nodes.append(list(g.nodes)[0])
-                else:
-                    nodes.append(Node(list(p), mass=mass))
-
-            prt.add_nodes(nodes)
-            element = element_cls(nodes=nodes, section=section, frame=frame)
+                if prt.nodes :
+                    if p not in prt.points:
+                        nodes.append(prt.add_node(Node(list(p), mass=mass)))
+                    else:
+                        nodes.append(list(prt.find_closest_nodes_to_point(p))[0])
+                else :
+                    nodes.append(prt.add_node(Node(list(p), mass=mass)))
+            element = element_cls(nodes=nodes, section=section, orientation=orientation)
             if not isinstance(element, _Element1D):
                 raise ValueError("Provide a 1D element")
             prt.add_element(element)
@@ -809,12 +812,10 @@ def faces(self) -> FacesGroup:
             raise ValueError("No elements in the part.")
         return FacesGroup([face for element in self.elements if element.faces is not None for face in element.faces])
 
-    # @property
-    # def releases(self) -> "ReleasesGroup | None":
-    #     """The releases of the part."""
-    #     for element in self.elements:
-    #         if hasattr(element, "releases") and element.releases is not None:
-    #             return ReleasesGroup(members=element.releases, name=f"{self.name}_releases_all")
+    @property
+    def releases_fields(self) -> "FieldsGroup | None":
+        """The releases of the part."""
+        return self._releases_fields
 
     @property
     def gkey_node(self) -> Dict[str, Node]:
@@ -2299,14 +2300,14 @@ def from_boundary_mesh(cls, boundary_mesh: "Mesh", section: Union["SolidSection"
     #                           Releases methods
     # =========================================================================
 
-    def add_beam_release(self, element: BeamElement, location: str, release: _BeamEndRelease) -> _BeamEndRelease:
+    def add_beams_releases_fields(self, elements: ElementsGroup, end: str, release: _BeamEndRelease) -> _BeamEndRelease:
         """Add a :class:`compas_fea2.model._BeamEndRelease` to an element in the part.
 
         Parameters
         ----------
-        element : :class:`compas_fea2.model.BeamElement`
+        elements : :class:`compas_fea2.model.ElementsGroup`
             The element to release.
-        location : str
+        end : str
             'start' or 'end'.
         release : :class:`compas_fea2.model._BeamEndRelease`
             Release type to apply.
@@ -2316,12 +2317,10 @@ def add_beam_release(self, element: BeamElement, location: str, release: _BeamEn
         :class:`compas_fea2.model._BeamEndRelease`
             The release applied to the element.
         """
-        raise NotImplementedError("Beam releases are not implemented in Part class. Use RigidPart instead.")
-        if not isinstance(release, _BeamEndRelease):
-            raise TypeError(f"{release!r} is not a beam release element.")
-        release.element = element
-        release.location = location
-        self._releases.add_member(release)
+
+        field = BeamReleaseField(elements=elements, end=end, release=release)
+        field._registration = self
+        self._releases_fields.add(field)
         return release
 
 

src/compas_fea2/problem/combinations.py

@@ -10,6 +10,7 @@
 from compas_fea2.base import FEAData
 from compas_fea2.base import from_data
 from compas_fea2.problem.groups import LoadsFieldGroup
+from compas_fea2.problem.fields import GravityLoadField
 
 if TYPE_CHECKING:
     from compas_fea2.model import Model
@@ -250,7 +251,11 @@ def select_factor(spec: FactorSpec, rank: int) -> float:
             factor = select_factor(spec, rank)
             if factor == 0.0:
                 continue
-            scaled_field = field if factor == 1.0 else factor * field
+            if isinstance(field, GravityLoadField):
+                field._g = factor * field._g
+                scaled_field = field
+            else:
+                scaled_field = field if factor == 1.0 else factor * field
             scaled_fields.append(scaled_field)
         return LoadsFieldGroup(members=scaled_fields)
 

src/compas_fea2/problem/fields.py

@@ -372,7 +372,7 @@ def __mul__(self, factor: float | int):
         if not isinstance(factor, (int, float)):
             return NotImplemented
         scaled = [self._scale_value(L, factor) for L in self._loads]
-        return self.__class__(scaled, self._distribution, load_case=self._load_case, combination_rank=self._combination_rank)
+        return self.__class__(scaled, distribution=self._distribution, load_case=self._load_case, combination_rank=self._combination_rank)
 
     def __rmul__(self, factor: float | int):
         return self.__mul__(factor)