diff --git a/.github/ISSUE_TEMPLATE/inquiry---need-help-for-my-project.md b/.github/ISSUE_TEMPLATE/inquiry---need-help-for-my-project.md index 20e02b22..b3ce0c72 100644 --- a/.github/ISSUE_TEMPLATE/inquiry---need-help-for-my-project.md +++ b/.github/ISSUE_TEMPLATE/inquiry---need-help-for-my-project.md @@ -12,7 +12,7 @@ A clear description of your problem. **To Reproduce** - URL for project repository -- Code implementation +- Code implementation (If your code is long, please use Gist and share the link.) **Expected behavior** A clear and concise description of what you expected to happen. diff --git a/examples/CatenaryCase/catenary.py b/examples/CatenaryCase/catenary.py new file mode 100644 index 00000000..e8f79671 --- /dev/null +++ b/examples/CatenaryCase/catenary.py @@ -0,0 +1,129 @@ +import numpy as np +from elastica import * + +from post_processing import ( + plot_video, + plot_catenary, +) + + +class CatenarySimulator(BaseSystemCollection, Constraints, Forcing, Damping, CallBacks): + pass + + +catenary_sim = CatenarySimulator() +final_time = 10 +damping_constant = 0.3 +time_step = 1e-4 +total_steps = int(final_time / time_step) +rendering_fps = 20 +step_skip = int(1.0 / (rendering_fps * time_step)) + +n_elem = 500 + +start = np.zeros((3,)) +direction = np.array([1.0, 0.0, 0.0]) +normal = np.array([0.0, 0.0, 1.0]) +binormal = np.cross(direction, normal) + +# catenary parameters +base_length = 1.0 +base_radius = 0.01 +base_area = np.pi * (base_radius ** 2) +volume = base_area * base_length +mass = 0.2 +density = mass / volume +E = 1e4 +poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) + +base_rod = CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus=E, + shear_modulus=shear_modulus, +) + +catenary_sim.append(base_rod) + +# add damping +catenary_sim.dampen(base_rod).using( + AnalyticalLinearDamper, + damping_constant=damping_constant, + time_step=time_step, +) + +# Add gravity +catenary_sim.add_forcing_to(base_rod).using( + GravityForces, acc_gravity=-9.80665 * normal +) + +# fix catenary ends +catenary_sim.constrain(base_rod).using( + FixedConstraint, constrained_position_idx=(0, -1), constrained_director_idx=(0, -1) +) + + +# Add call backs +class CatenaryCallBack(CallBackBaseClass): + """ + Call back function for continuum snake + """ + + def __init__(self, step_skip: int, callback_params: dict): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + + if current_step % self.every == 0: + + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["position"].append(system.position_collection.copy()) + self.callback_params["radius"].append(system.radius.copy()) + self.callback_params["internal_force"].append(system.internal_forces.copy()) + + return + + +pp_list = defaultdict(list) +catenary_sim.collect_diagnostics(base_rod).using( + CatenaryCallBack, step_skip=step_skip, callback_params=pp_list +) + + +catenary_sim.finalize() + + +timestepper = PositionVerlet() + +integrate(timestepper, catenary_sim, final_time, total_steps) +position = np.array(pp_list["position"]) +b = np.min(position[-1][2]) + +SAVE_VIDEO = True +if SAVE_VIDEO: + # plotting the videos + filename_video = "catenary.mp4" + plot_video( + pp_list, + video_name=filename_video, + fps=rendering_fps, + xlim=[0, base_length], + ylim=[-0.5 * base_length, 0.5 * base_length], + ) + +PLOT_RESULTS = True +if PLOT_RESULTS: + plot_catenary( + pp_list, + xlim=(0, base_length), + ylim=(b, 0.0), + ) diff --git a/examples/CatenaryCase/post_processing.py b/examples/CatenaryCase/post_processing.py new file mode 100644 index 00000000..ba6af0e3 --- /dev/null +++ b/examples/CatenaryCase/post_processing.py @@ -0,0 +1,76 @@ +import numpy as np +import matplotlib + +matplotlib.use("Agg") # Must be before importing matplotlib.pyplot or pylab! +from matplotlib import pyplot as plt +from tqdm import tqdm +import scipy as sci + + +def plot_video( + plot_params: dict, + video_name="video.mp4", + fps=15, + xlim=(0, 4), + ylim=(-1, 1), +): + import matplotlib.animation as manimation + + positions_over_time = np.array(plot_params["position"]) + + print("plot video") + FFMpegWriter = manimation.writers["ffmpeg"] + metadata = dict(title="Movie Test", artist="Matplotlib", comment="Movie support!") + writer = FFMpegWriter(fps=fps, metadata=metadata) + fig = plt.figure(figsize=(10, 8), frameon=True, dpi=150) + ax = fig.add_subplot(111) + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.set_xlabel("x [m]", fontsize=16) + ax.set_ylabel("y [m]", fontsize=16) + # plt.axis("equal") + with writer.saving(fig, video_name, dpi=150): + rod_lines_2d = ax.plot(positions_over_time[0][2], positions_over_time[0][0])[0] + for time in tqdm(range(1, len(plot_params["time"]))): + rod_lines_2d.set_xdata(positions_over_time[time][0]) + rod_lines_2d.set_ydata(positions_over_time[time][2]) + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + +def plot_catenary(plot_params: dict, xlim=(0, 1), ylim=(-0.5, 0.5), base_length=1.0): + """ + Catenary analytical solution from Routh, Edward John (1891). "Chapter X: On Strings". A Treatise on Analytical Statics. University Press. + """ + matplotlib.use("TkAgg") + position = np.array(plot_params["position"]) + lowest_point = np.min(position[-1][2]) + x_catenary = np.linspace(0, base_length, 100) + + def f_non_elastic_catenary(x): + return x * (1 - np.cosh(1 / (2 * x))) - lowest_point + + a = sci.optimize.fsolve(f_non_elastic_catenary, x0=1.0) # solve for a + y_catenary = a * np.cosh((x_catenary - 0.5) / a) - a * np.cosh(1 / (2 * a)) + plt.plot(position[-1][0], position[-1][2], label="Simulation", linewidth=3) + plt.plot( + x_catenary, + y_catenary, + label="Catenary (Analytical Solution)", + linewidth=3, + linestyle="dashed", + ) + plt.xlim(xlim) + plt.ylim(ylim) + plt.title("Catenary Final Shape") + plt.grid() + plt.legend() + plt.xlabel("x [m]", fontsize=16) + plt.ylabel("y [m]", fontsize=16) + plt.savefig("plot.png", dpi=300) + plt.show() diff --git a/examples/README.md b/examples/README.md index 489ea449..06f214be 100644 --- a/examples/README.md +++ b/examples/README.md @@ -26,17 +26,20 @@ Examples can serve as a starting template for customized usages. * __Purpose__: Physical convergence test of simple pendulum with flexible rod. * __Features__: CosseratRod, HingeBC, GravityForces * [ContinuumSnakeCase](./ContinuumSnakeCase) - * __Purpose__: Demonstrate simple case of modeling biological creature using PyElastica. The example use friction to create slithering snake, and optimize the speed using [CMA](https://github.com/CMA-ES/pycma). - * __Features__: CosseratRod, MuscleTorques, AnisotropicFrictionalPlane, Gravity, CMA Optimization - * [MuscularSnake](./MuscularSnake) - * __Purpose__: Example of [Parallel connection module](../elastica/experimental/connection_contact_joint/parallel_connection.py) and customized [Force module](./MuscularSnake/muscle_forces.py) to implement muscular snake. - * __Features__: MuscleForces(custom implemented) + * __Purpose__: Demonstrate simple case of modeling biological creature using PyElastica. The example uses friction to create slithering snake, and optimize the speed using [CMA](https://github.com/CMA-ES/pycma). + * __Features__: CosseratRod, MuscleTorques, RodPlaneContactWithAnisotropicFriction, Gravity, CMA Optimization +* [ContinuumSnakeWithLiftingWaveCase](./ContinuumSnakeWithLiftingWaveCase) + * __Purpose__: Demonstrate simple case of modeling biological creature using PyElastica. The example uses friction to create slithering snake with lift. + * __Features__: CosseratRod, MuscleTorquesLifting(custom implemented), SnakeRodPlaneContact(custom implemented), Gravity +* [MuscularSnake](./MuscularSnake) + * __Purpose__: Example of [Parallel connection module](../elastica/experimental/connection_contact_joint/parallel_connection.py) and customized [Force module](./MuscularSnake/muscle_forces.py) to implement muscular snake. + * __Features__: MuscleForces(custom implemented) * [ButterflyCase](./ButterflyCase) * __Purpose__: Demonstrate simple restoration with initial strain. * __Features__: CosseratRod * [FrictionValidationCases](./FrictionValidationCases) * __Purpose__: Physical validation of rolling and translational friction. - * __Features__: CosseratRod, UniformForces, AnisotropicFrictionalPlane + * __Features__: CosseratRod, UniformForces, RodPlaneContactWithAnisotropicFriction * [JointCases](./JointCases) * __Purpose__: Demonstrate various joint usage with Cosserat Rod. * __Features__: FreeJoint, FixedJoint, HingeJoint, OneEndFixedRod, EndpointForcesSinusoidal @@ -45,27 +48,27 @@ Examples can serve as a starting template for customized usages. * __Features__: Cylinder, Sphere * [RodRigidBodyContact](./RigidBodyCases/RodRigidBodyContact) * __Purpose__: Demonstrate contact between cylinder and rod, for different intial conditions. - * __Features__: Cylinder, CosseratRods, ExternalContact + * __Features__: Cylinder, CosseratRods, RodCylinderContact * [HelicalBucklingCase](./HelicalBucklingCase) * __Purpose__: Demonstrate helical buckling with extreme twisting boundary condition. * __Features__: HelicalBucklingBC * [ContinuumFlagellaCase](./ContinuumFlagellaCase) * __Purpose__: Demonstrate flagella modeling using PyElastica. * __Features__: SlenderBodyTheory, MuscleTorques, - * [MuscularFlagella](./MuscularFlagella) - * __Purpose__: Example of customizing [Joint module](./MuscularFlagella/connection_flagella.py) and [Force module](./MuscularFlagella/muscle_forces_flagella.py) to implement muscular flagella. - * __Features__: MuscleForces(custom implemented) +* [MuscularFlagella](./MuscularFlagella) + * __Purpose__: Example of customizing [Joint module](./MuscularFlagella/connection_flagella.py) and [Force module](./MuscularFlagella/muscle_forces_flagella.py) to implement muscular flagella. + * __Features__: MuscleForces(custom implemented) * [RodContactCase](./RodContactCase) * [RodRodContact](./RodContactCase/RodRodContact) * __Purpose__: Demonstrates contact between two rods, for different initial conditions. - * __Features__: CosseratRod, ExternalContact + * __Features__: CosseratRod, RodRodContact * [RodSelfContact](./RodContactCase/RodSelfContact) * [PlectonemesCase](./RodContactCase/RodSelfContact/PlectonemesCase) * __Purpose__: Demonstrates rod self contact with Plectoneme example, and how to use link-writhe-twist after simulation completed. - * __Features__: CosseratRod, SelonoidsBC, SelfContact, Link-Writhe-Twist + * __Features__: CosseratRod, SelonoidsBC, RodSelfContact, Link-Writhe-Twist * [SolenoidsCase](./RodContactCase/RodSelfContact/SolenoidsCase) * __Purpose__: Demonstrates rod self contact with Solenoid example, and how to use link-writhe-twist after simulation completed. - * __Features__: CosseratRod, SelonoidsBC, SelfContact, Link-Writhe-Twist + * __Features__: CosseratRod, SelonoidsBC, RodSelfContact, Link-Writhe-Twist * [BoundaryConditionsCases](./BoundaryConditionsCases) * __Purpose__: Demonstrate the usage of boundary conditions for constraining the movement of the system. * __Features__: GeneralConstraint, CosseratRod @@ -75,6 +78,9 @@ Examples can serve as a starting template for customized usages. * [RingRodCase](./RingRodCase) * __Purpose__: Demonstrate simulation of ring rod. * __Features__: RingCosseratRod, OneEndFixedRod, GravityForce +* [CatenaryCase](./CatenaryCase) + * __Purpose__: Demonstrate simulation of cosserat rod under gravity with fixed ends, compared with Catenary Analytical Solution from Routh, Edward John (1891). ["Chapter X: On Strings"](https://books.google.com/books?id=3N5JAAAAMAAJ&pg=PA315#v=onepage&q&f=false). A Treatise on Analytical Statics. University Press. + * __Features__: CosseratRod, FixedConstraint, GravityForce ## Functional Examples