Merge pull request #255 from adtzlr/add-mesh-runout

Add Mesh-tool `mesh.runouts()`

felupe/__init__.py

@@ -104,6 +104,7 @@
     SolidBodyTensor,
     SolidBodyGravity,
 )
+
 try:
     from .mechanics import (
         MultiPointConstraint,

felupe/_field/_container.py

@@ -79,7 +79,7 @@ def extract(self, grad=True, sym=False, add_identity=True):
     def copy(self):
         "Return a copy of the field."
         return deepcopy(self)
-    
+
     def link(self, other_field):
         "Link value array of other field."
         for field, newfield in zip(self.fields, other_field.fields):

felupe/_field/_fields.py

@@ -47,7 +47,15 @@ class FieldsMixed(FieldContainer):
     r"""A mixed field based on a region and returns a :class:`FieldContainer`
     instance."""
 
-    def __init__(self, region, n=3, values=(0, 0, 1, 0), axisymmetric=False, offset=0, npoints=None):
+    def __init__(
+        self,
+        region,
+        n=3,
+        values=(0, 0, 1, 0),
+        axisymmetric=False,
+        offset=0,
+        npoints=None,
+    ):
         r"""Create a mixed field based on a region. The dual region is chosen
         automatically, i.e. for a :class:`RegionHexahedron` the dual region
         is :class:`RegionConstantHexahedron`. A total number of ``n`` fields
@@ -80,7 +88,7 @@ def __init__(self, region, n=3, values=(0, 0, 1, 0), axisymmetric=False, offset=
             RegionTetraMINI: RegionTetra,
             RegionTriangleMINI: RegionTriangle,
         }
-        
+
         kwargs = {}
         if offset > 0:
             kwargs["offset"] = offset

felupe/mechanics/__init__.py

@@ -2,10 +2,11 @@
 from ._solidbody_pressure import SolidBodyPressure
 from ._solidbody_tensor import SolidBodyTensor
 from ._solidbody_gravity import SolidBodyGravity
+
 try:
     from ._multipoint import (
         MultiPointConstraint,
         MultiPointContact,
     )
 except:
-    pass
+    pass

felupe/mechanics/_helpers.py

@@ -63,8 +63,8 @@ def __init__(self, stress=False, elasticity=False):
 
         if elasticity:
             self.elasticity = None
-        
+
     def update_statevars(self):
-        
+
         if self._statevars is not None:
             self.statevars = self._statevars

felupe/mechanics/_multipoint.py

@@ -43,10 +43,10 @@ def __init__(
         self.mask = ~np.array(skip, dtype=bool)[: self.mesh.dim]
         self.axes = np.arange(self.mesh.dim)[self.mask]
         self.multiplier = multiplier
-        
+
         self.results = Results(stress=False, elasticity=False)
         self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
-    
+
     def _vector(self, field=None, parallel=False, jit=False):
         "Calculate vector of residuals with RBE2 contributions."
         if field is not None:
@@ -98,7 +98,7 @@ def __init__(
         self.mask = ~np.array(skip, dtype=bool)[: self.mesh.dim]
         self.axes = np.arange(self.mesh.dim)[self.mask]
         self.multiplier = multiplier
-        
+
         self.results = Results(stress=False, elasticity=False)
         self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
 
@@ -122,7 +122,6 @@ def _vector(self, field=None):
         self.results.force = sparse.COO(r).reshape((-1, 1)).tocsr()
         return self.results.force
 
-
     def _matrix(self, field=None):
         "Calculate stiffness with RBE2 contributions."
         if field is not None:
@@ -152,4 +151,4 @@ def _matrix(self, field=None):
             )
             .tocsr()
         )
-        return self.results.stiffness
+        return self.results.stiffness

felupe/mechanics/_solidbody_gravity.py

@@ -41,37 +41,33 @@ def __init__(self, field, gravity, density):
         self.results = Results(stress=False, elasticity=False)
         self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
         self._form = IntegralFormMixed
-        
+
         self.results.gravity = np.array(gravity)
         self.results.density = density
 
-    def _vector(
-        self, field=None, parallel=False, jit=False
-    ):
+    def _vector(self, field=None, parallel=False, jit=False):
 
         if field is not None:
             self.field = field
-        
+
         # copy and take only the first (displacement) field of the container
         f = self.field.copy()
-        f.fields = f.fields[0: 1]
-            
+        f.fields = f.fields[0:1]
+
         self.results.force = self._form(
             fun=[self.results.density * self.results.gravity.reshape(-1, 1, 1)],
-            v=f, 
+            v=f,
             dV=self.field.region.dV,
             grad_v=[False],
         ).assemble(parallel=parallel, jit=jit)
-        
+
         if len(self.field) > 1:
             self.results.force.resize(np.sum(self.field.fieldsizes), 1)
 
         return self.results.force
 
-    def _matrix(
-        self, field=None, parallel=False, jit=False
-    ):
-
+    def _matrix(self, field=None, parallel=False, jit=False):
+
         from scipy.sparse import csr_matrix
 
         if field is not None:
@@ -80,4 +76,4 @@ def _matrix(
         n = np.sum(self.field.fieldsizes)
         self.results.stiffness = csr_matrix(([0], ([0], [0])), shape=(n, n))
 
-        return self.results.stiffness
+        return self.results.stiffness

felupe/mechanics/_solidbody_pressure.py

@@ -41,7 +41,7 @@ def __init__(self, field, pressure=None):
 
         self.results = Results()
         self.results.kinematics = self._extract(self.field)
-        
+
         if pressure is not None:
             self.results.pressure = pressure
         else:
@@ -75,7 +75,7 @@ def _vector(
 
         if pressure is not None:
             self.results.pressure = pressure
-        
+
         fun[0] *= self.results.pressure
 
         self.results.force = IntegralFormMixed(
@@ -104,7 +104,7 @@ def _matrix(
 
         if pressure is not None:
             self.results.pressure = pressure
-            
+
         fun[0] *= self.results.pressure
 
         self.results.stiffness = IntegralFormMixed(

felupe/mechanics/_solidbody_tensor.py

@@ -44,15 +44,17 @@ def __init__(self, umat, field, statevars=None):
 
         self.results = Results(stress=True, elasticity=True)
         self.results.kinematics = self._extract(self.field)
-        
+
         if statevars is not None:
             self.results.statevars = statevars
         else:
-            self.results.statevars = np.zeros((
-                *umat.x[-1].shape, 
-                field.region.quadrature.npoints, 
-                field.region.mesh.ncells,
-            ))
+            self.results.statevars = np.zeros(
+                (
+                    *umat.x[-1].shape,
+                    field.region.quadrature.npoints,
+                    field.region.mesh.ncells,
+                )
+            )
 
         self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
 
@@ -114,7 +116,7 @@ def _gradient(self, field=None, args=(), kwargs={}):
         if field is not None:
             self.field = field
             self.results.kinematics = self._extract(self.field)
-            
+
         function = self.umat.function(
             [*self.results.kinematics, self.results.statevars], *args, **kwargs
         )

felupe/mesh/__init__.py

@@ -6,6 +6,7 @@
     sweep,
     mirror,
     triangulate,
+    runouts,
 )
 from ._convert import (
     convert,

felupe/mesh/_tools.py

@@ -407,3 +407,68 @@ def triangulate(points, cells, cell_type, mode=3):
     cells_new = cells_new.reshape(-1, cells_new.shape[-1])
 
     return points, cells_new, cell_type_new
+
+
+@mesh_or_data
+def runouts(
+    points,
+    cells,
+    cell_type,
+    values=[0.1, 0.1],
+    centerpoint=[0, 0, 0],
+    axis=0,
+    exponent=5,
+):
+    """Add simple rubber-runouts for realistic rubber-metal structures.
+
+    Parameters
+    ----------
+    points : list or ndarray
+        Original point coordinates.
+    cells : list or ndarray
+        Original point-connectivity of cells.
+    cell_type : str
+        A string in VTK-convention that specifies the cell type.
+    values : list or ndarray, optional
+        Relative amount of runouts (per coordinate) perpendicular to the axis
+        (default is 10% per coordinate, i.e. [0.1, 0.1]).
+    centerpoint: list or ndarray, optional
+        Center-point coordinates (default is [0, 0, 0]).
+    axis: int or None, optional
+        Axis (default is 0).
+    exponent: int, optional
+        Positive exponent to control the shape of the runout. The higher
+        the exponent, the steeper the transition (default is 5).
+
+    Returns
+    -------
+    points : ndarray
+        Modified point coordinates.
+    cells : ndarray
+        Modified point-connectivity of cells.
+    cell_type : str or None
+        A string in VTK-convention that specifies the cell type.
+    """
+
+    dim = points.shape[1]
+    runout_along = {0: [1, 2], 1: [0, 2], 2: [0, 1]}
+
+    centerpoint = np.array(centerpoint, dtype=float)[:dim]
+    values = np.array(values, dtype=float)[:dim]
+
+    points_new = points - centerpoint
+    top = points[:, axis].max()
+    bottom = points[:, axis].min()
+
+    # check symmetry
+    if top == centerpoint[axis] or bottom == centerpoint[axis]:
+        half_height = top - bottom
+    else:
+        half_height = (top - bottom) / 2
+
+    for i, coord in enumerate(runout_along[axis][: dim - 1]):
+
+        factor = (abs(points_new[:, axis]) / half_height) ** exponent
+        points_new[:, coord] *= 1 + factor * values[i]
+
+    return points_new + centerpoint, cells, cell_type

felupe/region/_templates.py

@@ -60,7 +60,7 @@ def __init__(self, mesh, offset=0, npoints=None):
 
         element = ConstantQuad()
         quadrature = GaussLegendre(order=1, dim=2)
-        
+
         if npoints is not None:
             npts = npoints
         else:
@@ -109,7 +109,7 @@ def __init__(self, mesh, offset=0, npoints=None):
 
         element = ConstantHexahedron()
         quadrature = GaussLegendre(order=1, dim=3)
-        
+
         if npoints is not None:
             npts = npoints
             ncells = min(npoints, mesh.ncells)

felupe/tools/_newton.py

@@ -136,19 +136,19 @@ def solve(A, b, x, dof1, dof0, offsets=None, ext0=None, solver=spsolve):
 
 def check(dx, x, f, xtol, ftol, dof1=None, dof0=None, items=None, eps=1e-3):
     "Check result."
-    
+
     sumnorm = lambda x: np.sum(norm(x))
     xnorm = sumnorm(dx)
-    
+
     if dof1 is None:
         dof1 = slice(None)
-    
+
     if dof0 is None:
         dof0 = slice(0, 0)
-    
+
     fnorm = sumnorm(f[dof1]) / (eps + sumnorm(f[dof0]))
-    success = fnorm < ftol # and xnorm < xtol
-    
+    success = fnorm < ftol  # and xnorm < xtol
+
     if success and items is not None:
         for item in items:
             [item.results.update_statevars() for item in items]
@@ -158,7 +158,7 @@ def check(dx, x, f, xtol, ftol, dof1=None, dof0=None, items=None, eps=1e-3):
 
 def update(x, dx):
     "Update field."
-    #x += dx # in-place
+    # x += dx # in-place
     return x + dx
 
 
@@ -226,10 +226,10 @@ def newtonrhapson(
 
     if umat is not None:
         kwargs["umat"] = umat
-    
+
     if kwargs_solve is None:
         kwargs_solve = {}
-    
+
     if kwargs_check is None:
         kwargs_check = {}
 
@@ -277,7 +277,9 @@ def newtonrhapson(
             f = fun(x, *args, **kwargs)
 
         # check success of solution
-        xnorm, fnorm, success = check(dx, x, f, tol, tol, dof1, dof0, items, **kwargs_check)
+        xnorm, fnorm, success = check(
+            dx, x, f, tol, tol, dof1, dof0, items, **kwargs_check
+        )
 
         if verbose:
             print("|%2d | %1.3e | %1.3e |" % (1 + iteration, fnorm, xnorm))