Add mean motion constraint and use rotation vector as default convention

docs/boundary_conditions.rst

@@ -105,9 +105,18 @@ ____________________________________________
    fedoo.MPC
 
 
-Avdanced BC and constraints
+Advanced BC and constraints
 ==============================
 
+Fedoo also provides higher-level boundary-condition objects for common
+distributed loads and constraint patterns. Distributed loads are converted to
+external Neumann forces, while the kinematic constraints below are based on
+MPC equations and can be added directly to ``pb.bc``.
+
+
+Distributed loads
+_________________
+
 Distributed loads are defined by dedicated assembly objects. The recommended
 way to apply them is to add the load assembly to the problem boundary
 conditions. If the assembly provides an ``as_neumann()`` method,
@@ -147,18 +156,118 @@ boundary-condition force.
     fedoo.constraint.Pressure
     fedoo.constraint.SurfaceForce
 
-Some advanced constraints based on multiple linearized MPC are available in
-fedoo. They can be created and added to the problem with the pb.bc.add method.
+
+Kinematic MPC constraints
+_________________________
+
+Rigid tie
+~~~~~~~~~
+
+The :py:class:`fedoo.constraint.RigidTie` constraint couples a set of nodes to
+a rigid-body motion. It creates global degrees of freedom for rigid
+translation and rotation, then eliminates the selected nodal displacement DOFs
+through MPC equations. Prescribing ``RigidDisp`` or ``RigidRot`` therefore
+drives the whole selected node set as a rigid object.
+
+This constraint is useful to impose a strict rigid connection, for instance to
+drive a face from a reference rigid point. It is more restrictive than
+:py:class:`fedoo.constraint.MeanMotion`, because all selected nodes follow the
+same rigid-body kinematics and local warping of the selected surface is
+suppressed.
 
 .. autosummary::
    :toctree: generated/
    :template: custom-class-template.rst
 
-   fedoo.constraint.PeriodicBC
    fedoo.constraint.RigidTie
    fedoo.constraint.RigidTie2D
 
 
+Mean motion control
+~~~~~~~~~~~~~~~~~~~
+
+The :py:class:`fedoo.constraint.MeanMotion` constraint is a less
+restrictive alternative to :py:class:`fedoo.constraint.RigidTie`. It extracts
+selected components of the best-fit motion of a node set or surface mesh in a
+mean sense, without forcing every selected node to follow a rigid motion. The
+selected surface can therefore deform locally while selected mean translation
+and rotation components are controlled.
+
+The ``components`` argument is required and defines which global DOFs are
+created. Accepted component names include ``"RotX"``, ``"MeanRotX"``,
+``"DispZ"`` and ``"MeanDispZ"``. The vector aliases ``"Rot"``/``"MeanRot"``
+and ``"Disp"``/``"MeanDisp"`` select all rotation or displacement components
+available in the current modeling dimension. In 2D, the rotation vector
+contains only ``MeanRotZ``.
+
+When a surface mesh is passed, nodal weights are computed from the element
+area, or from the element length in 2D. These weights define an area-weighted
+least-squares projection of the displacement field onto rigid-body modes.
+
+For displacement-only mean control, select only displacement components:
+
+.. code-block:: python
+
+    mean_disp = fd.constraint.MeanMotion(surface_mesh, components="DispZ")
+    pb.bc.add(mean_disp)
+    pb.bc.add("Dirichlet", "MeanDispZ", -delta)
+
+When ``finite_rotation`` is left to ``None`` (the default), rotational
+components use finite rotations automatically if geometrical nonlinearity is
+enabled on the problem. Otherwise, the rotation components are small-rotation
+vector components. The finite-rotation mode uses a rotation-vector convention
+and linearizes the constraint at each Newton iteration using the derivative of
+the rotation matrix with respect to the rotation vector.
+
+.. code-block:: python
+
+    face_motion = fd.constraint.MeanMotion(
+        surface_mesh,
+        components="RotZ",
+    )
+    pb.bc.add(face_motion)
+    pb.bc.add("Dirichlet", face_motion.node_by_variable["MeanRotZ"], "MeanRotZ", angle)
+
+.. code-block:: python
+
+    face_motion = fd.constraint.MeanMotion(
+        surface_mesh,
+        components=["RotZ", "DispZ"],
+    )
+    pb.bc.add(face_motion)
+    pb.bc.add("Dirichlet", "MeanRotZ", angle)
+    pb.bc.add("Dirichlet", "MeanDispZ", -delta)
+
+The reactions on the created global DOFs are generalized forces. For instance,
+the reaction associated with ``MeanRotZ`` is the conjugate moment around the
+mean-motion center.
+
+.. autosummary::
+   :toctree: generated/
+   :template: custom-class-template.rst
+
+   fedoo.constraint.MeanMotion
+
+
+Periodic boundary conditions
+____________________________
+
+The :py:class:`fedoo.constraint.PeriodicBC` constraint creates MPC equations
+between opposite boundaries so that their displacement fluctuations are
+periodic. It is mainly used for representative volume elements and
+homogenization problems, where the boundary kinematics must be compatible with
+a prescribed macroscopic strain or displacement gradient.
+
+The constraint can build the required node pairings from opposite node sets
+and then be added to ``pb.bc`` like the other MPC-based constraints.
+
+.. autosummary::
+   :toctree: generated/
+   :template: custom-class-template.rst
+
+   fedoo.constraint.PeriodicBC
+
+
 Contact
 ==============================
 

examples/01-simple/spherical_shell_compression.py

@@ -70,7 +70,7 @@
 #
 assembly = solid_assembly
 
-pb = fd.problem.NonLinear(assembly, nlgeom=False)
+pb = fd.problem.Linear(assembly)
 pb.bc.add(pressure_assembly)
 
 nodes = mesh.nodes
@@ -82,7 +82,7 @@
 pb.bc.add("Dirichlet", node_b, ["DispY", "DispZ"], 0)
 pb.bc.add("Dirichlet", node_c, "DispZ", 0)
 
-pb.nlsolve()
+pb.solve()
 
 ###############################################################################
 # Extract the results:

examples/02-constraints/mean_rigid_motion_torsion.py

@@ -0,0 +1,167 @@
+"""
+Mean motion torsion
+~~~~~~~~~~~~~~~~~~~
+
+Compare two ways of imposing a torsion on the right face of a cube:
+
+* ``RigidTie`` enforces a rigid motion of every node on the right face.
+* ``MeanMotion`` prescribes only the best-fit mean rotation of the face,
+  so local warping of the right face remains possible.
+"""
+
+import numpy as np
+import pyvista as pv
+import fedoo as fd
+
+
+fd.ModelingSpace("3D")
+
+NLGEOM = "UL"
+E = 200e3
+nu = 0.3
+L = 2.0
+ANGLE = np.pi / 3  # 0.10
+PRINT_INFO = 0
+
+
+def make_model(label):
+    mesh = fd.mesh.box_mesh(
+        nx=8,
+        ny=12,
+        nz=12,
+        x_min=0.0,
+        x_max=L,
+        y_min=-1.0,
+        y_max=1.0,
+        z_min=-1.0,
+        z_max=1.0,
+        elm_type="hex8",
+        name=f"torsion_{label}_mesh",
+    )
+
+    material = fd.constitutivelaw.ElasticIsotrop(E, nu, name=f"material_{label}")
+    wf = fd.weakform.StressEquilibrium(material, nlgeom=NLGEOM, name=f"wf_{label}")
+    assembly = fd.Assembly.create(wf, mesh, "hex8", name=f"assembly_{label}")
+
+    pb = fd.problem.NonLinear(assembly, nlgeom=NLGEOM, name=f"problem_{label}")
+    # pb.set_nr_criterion("Displacement", err0=1, tol=1e-3, max_subiter=20)
+    return mesh, assembly, pb
+
+
+def solve_rigid_tie_case():
+    mesh, assembly, pb = make_model("rigid")
+    left = mesh.find_nodes("X", mesh.bounding_box.xmin)
+    right = mesh.find_nodes("X", mesh.bounding_box.xmax)
+
+    pb.bc.add(fd.constraint.RigidTie(right))
+    pb.bc.add("Dirichlet", left, "Disp", 0.0)
+    pb.bc.add("Dirichlet", "RigidRotX", ANGLE)
+
+    pb.nlsolve(dt=0.1, tmax=1.0, update_dt=True, print_info=PRINT_INFO)
+
+    moment = pb.get_ext_forces("RigidRotX")[0]
+    disp = pb.get_dof_solution("Disp")
+    return pb, assembly, moment
+
+
+def solve_mean_motion_case():
+    mesh, assembly, pb = make_model("mean")
+    left = mesh.find_nodes("X", mesh.bounding_box.xmin)
+    right = mesh.find_nodes("X", mesh.bounding_box.xmax)
+    right_surface = fd.mesh.extract_surface(mesh, node_set=right, reduce_order=False)
+
+    mean_motion = fd.constraint.MeanMotion(
+        right_surface,
+        components=["Rot"],
+    )
+    pb.bc.add(mean_motion)
+    pb.bc.add("Dirichlet", left, "Disp", 0.0)
+    pb.bc.add("Dirichlet", "MeanRotX", ANGLE)
+    # pb.bc.add("Dirichlet", ["MeanDispX", "MeanDispY", "MeanDispZ"], [0.4, -1, 0])
+
+    pb.nlsolve(dt=0.1, tmax=1.0, update_dt=True, print_info=PRINT_INFO)
+
+    moment = pb.get_ext_forces("MeanRotX")[mean_motion.node_by_variable["MeanRotX"]]
+    disp = pb.get_dof_solution("Disp")
+    return pb, assembly, moment
+
+
+rigid_pb, rigid_assembly, rigid_moment = solve_rigid_tie_case()
+mean_pb, mean_assembly, mean_moment = solve_mean_motion_case()
+
+
+###############################################################################
+# Plot the two constraints side by side
+# -------------------------------------
+#
+# The deformed shapes are colored by the axial displacement ``DispX``.  The
+# rigid tie keeps the right face rigid, while the mean-motion constraint only
+# enforces the best-fit rotation and lets the face warp locally.
+
+rigid_res = rigid_pb.get_results(rigid_assembly, ["Disp", "Stress"], "Node")
+mean_res = mean_pb.get_results(mean_assembly, ["Disp", "Stress"], "Node")
+
+disp_x_rigid = rigid_res.get_data("Disp", component="X", data_type="Node")
+disp_x_mean = mean_res.get_data("Disp", component="X", data_type="Node")
+disp_x_lim = [
+    min(np.min(disp_x_rigid), np.min(disp_x_mean)),
+    max(np.max(disp_x_rigid), np.max(disp_x_mean)),
+]
+
+plotter = pv.Plotter(shape=(1, 2), window_size=(1400, 650))
+
+plotter.subplot(0, 0)
+rigid_res.plot(
+    "Disp",
+    component="X",
+    data_type="Node",
+    scale=1,
+    clim=disp_x_lim,
+    cmap="coolwarm",
+    show_edges=True,
+    show_scalar_bar=False,
+    lock_view=True,
+    plotter=plotter,
+    title="RigidTie",
+)
+plotter.hide_axes()
+plotter.add_text(
+    "Rigid face kinematics\n" f"Mx = {rigid_moment:.3e}\n",
+    position="lower_left",
+    font_size=10,
+)
+
+plotter.subplot(0, 1)
+mean_res.plot(
+    "Disp",
+    component="X",
+    data_type="Node",
+    scale=1,
+    clim=disp_x_lim,
+    cmap="coolwarm",
+    show_edges=True,
+    show_scalar_bar=False,
+    lock_view=True,
+    plotter=plotter,
+    title="MeanMotion",
+)
+plotter.hide_axes()
+plotter.add_text(
+    "Best-fit mean rotation\n" f"Mx = {mean_moment:.3e}\n",
+    position="lower_left",
+    font_size=10,
+)
+plotter.add_scalar_bar(
+    title="DispX",
+    vertical=True,
+    position_x=0.91,
+    position_y=0.18,
+    height=0.62,
+    width=0.04,
+    title_font_size=18,
+    label_font_size=14,
+)
+
+plotter.link_views()
+plotter.view_isometric()
+plotter.show()

fedoo/constraint/__init__.py

@@ -7,3 +7,4 @@
     SurfaceForce,
 )
 from .ipc_contact import IPCContact, IPCSelfContact
+from .mean_motion import MeanMotion