Fix Line and Plane constraint

[PatrickRoel]

When attempting to simulate a system where a node has a line or plane constraints the line_mask or plane_mask shape is set up incorrect.
When simulating a mesh a multiplication error arises between multiplication of mask *= line_mask.

the issue is in the following code snippet from ParticleSystem.py

        for p in self.__particles:
            if p.constraint_type == "plane":
                for i in range(3):
                    line_mask.append(True)
            if p.constraint_type == "plane":
                for i in range(3):
                    line_mask.append(True)
                constraint = p._Particle__constraint[0]
                if constraint[0] == 1:
                    plane_mask.append(False)
                    plane_mask.append(True)
                    plane_mask.append(True)
                elif constraint[1] == 1:
                    plane_mask.append(True)
                    plane_mask.append(False)
                    plane_mask.append(True)
                elif constraint[2] == 1:
                    plane_mask.append(True)
                    plane_mask.append(True)
                    plane_mask.append(False)
                else:
                    for i in range(3):
                        plane_mask.append(True)
                    for i in range(3):
                        plane_mask.append(True)

            elif p.constraint_type == "line":
                for i in range(3):
                    plane_mask.append(True)
            elif p.constraint_type == "line":
                for i in range(3):
                    plane_mask.append(True)
                constraint = p._Particle__constraint[0]
                if constraint[0] == 1:
                    line_mask.append(True)
                    line_mask.append(False)
                    line_mask.append(False)
                elif constraint[1] == 1:
                    line_mask.append(False)
                    line_mask.append(True)
                    line_mask.append(False)
                elif constraint[2] == 1:
                    line_mask.append(False)
                    line_mask.append(False)
                    line_mask.append(True)
                else:
                    for i in range(3):
                        line_mask.append(True)
                    for i in range(3):
                        line_mask.append(True)
            else:
                for i in range(3):
                    plane_mask.append(True)
                    line_mask.append(True)

        mask = np.outer(point_mask, [True, True, True]).flatten()
        mask *= plane_mask
        mask *= line_mask

on further investigation it seems that the following piece of code is run twice instead of once

for a line

            elif p.constraint_type == "line":
                for i in range(3):
                    plane_mask.append(True)

for a plane

            if p.constraint_type == "plane":
                for i in range(3):
                    line_mask.append(True)

removing this part of the code fixes the issue.

[PatrickRoel]

Example code using a "line" type constraint

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys, os

# sys.path.append(os.path.abspath(".."))
from PSS.particleSystem.ParticleSystem import ParticleSystem

# dictionary of required parameters
params = {
    # model parameters
    "k": 2000,  # [N/m]     spring stiffness
    "c": 20,  # [N s/m] damping coefficient
    "m_segment": 1,  # [kg] mass of each node
    # simulation settings
    "dt": 0.1,  # [s]       simulation timestep
    "t_steps": 10000,  # [-]      number of simulated time steps
    "abs_tol": 1e-50,  # [m/s]     absolute error tolerance iterative solver
    "rel_tol": 1e-5,  # [-]       relative error tolerance iterative solver
    "max_iter": 1e5,  # [-]       maximum number of iterations
    # physical parameters
    "g": 9.807,  # [m/s^2]   gravitational acceleration
}

params["l0"] = .25
params["n"] = 8


# Setting up the coordinates and connectivity matrix for the nodes
xyz_coordinates = np.empty((params["n"], 3))
initial_conditions = []
connections = []

for i in range(params["n"]):
    xyz_coordinates[i] = [i,0, 0]
    if i == 0:
        initial_conditions.append([xyz_coordinates[i], np.zeros(3), params["m_segment"], True, [0,0,0],"point"])
    elif i == params["n"]-1:
        initial_conditions.append([xyz_coordinates[i], np.zeros(3), params["m_segment"], True, [1,0,0],"line"])
        connections.append([i-1,i,params["k"], params["c"]])
    else:
        initial_conditions.append([xyz_coordinates[i], np.zeros(3), params["m_segment"], False, [0,0,0],"point"])
        connections.append([i-1,i,params["k"], params["c"]])


connections.append([0,params["n"]-1,params["k"]/1000, params["c"]])


# Setting up external forces
f_ext = np.empty((params["n"], 3))
for i in range(params["n"]):
    if i == 0 or i == params["n"]-1:
        f_ext[i] = [0, 0, 0]
    else:
        f_ext[i] = [0, .5, 0]

       
f_ext = f_ext.flatten()

# Plotting mesh with external forces
# Create a 2D plot
fig, ax = plt.subplots()

# Plot the nodes
for i, node in enumerate(initial_conditions):
    if node[3]:  # Fixed node
        if node[5] == "line":
            ax.scatter(node[0][0], node[0][1], color="purple", marker="o", label="Line Fixed Node" if i == params["n"]-1 else "")

        else:
            ax.scatter(node[0][0], node[0][1], color="red", marker="o", label="Fixed Node" if i == 0 else "")
        
    else:  # Free node
        ax.scatter(node[0][0], node[0][1], color="blue", marker="o", label="Free Node" if i == 1 else "")

    # Plot the external forces as arrows
    ax.quiver(
        node[0][0], node[0][1], f_ext[3 * i], f_ext[3 * i + 1], angles="xy", scale_units="xy", scale=.1, color="green"
    )

# Plot the connections between nodes
for i,connection in enumerate(connections):
    line = np.column_stack(
        [initial_conditions[connection[0]][0][:2], initial_conditions[connection[1]][0][:2]]
    )
    if i == params["n"]-1:
        ax.plot(line[0], line[1], color="black",linestyle="--")
    else:
        ax.plot(line[0], line[1], color="black")

# Set axis limits
ax.set_xlim(-1, 8)
ax.set_ylim(-10, 5)

# Add labels and title
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("Initial State")

# Add legend
ax.legend()

# Show the plot
plt.show()

# Now we can setup the particle system and simulation
PS = ParticleSystem(connections, initial_conditions, params,init_surface=False)

t_vector = np.linspace(
    params["dt"], params["t_steps"] * params["dt"], params["t_steps"]
)
final_step = 0
E_kin = []
f_int = []

f_external = f_ext
# print(
#     f"f_external: shape: {np.shape(f_external)}, average: {np.mean(f_external)}, min: {np.min(f_external)}, max: {np.max(f_external)}"
# )

# And run the simulation
for step in t_vector:
    PS.kin_damp_sim(f_ext)

    final_step = step
    (
        x,
        v,
    ) = PS.x_v_current
    E_kin.append(np.linalg.norm(v * v))
    f_int.append(np.linalg.norm(PS.f_int))

    converged = False
    if step > 10:
        if np.max(E_kin[-10:-1]) <= 1e-29:
            converged = True
    if converged and step > 1:
        print("Kinetic damping PS converged", step)
        break

particles = PS.particles
particle = particles[1]

final_positions = [
    [particle.x, particle.v, particle.m, particle.fixed, particle.constraint_type] for particle in PS.particles
]

# Plotting final results in 2D
fig, ax = plt.subplots()

# Plot the nodes
for i, node in enumerate(final_positions):
    if node[3]:  # Fixed node
        if node[4] == "line":
            ax.scatter(node[0][0], node[0][1], color="purple", marker="o", label="Line Fixed Node" if i == params["n"]-1 else "")
        else:
            ax.scatter(node[0][0], node[0][1], color="red", marker="o", label="Fixed Node" if i == 0 else "")
    else:  # Free node
        ax.scatter(node[0][0], node[0][1], color="blue", marker="o", label="Free Node" if i == 1 else "")

    # Plot the external forces as arrows
    ax.quiver(
        node[0][0], node[0][1], f_ext[3 * i], f_ext[3 * i + 1], angles="xy", scale_units="xy", scale=.1, color="green"
        )

# Plot the connections between nodes
for i,connection in enumerate(connections):
    line = np.column_stack(
        [final_positions[connection[0]][0][:2], final_positions[connection[1]][0][:2]]
    )
    if i == params["n"]-1:
        ax.plot(line[0], line[1], color="black",linestyle="--")
    else:
        ax.plot(line[0], line[1], color="black")

# Set axis limits
ax.set_xlim(-1, 8)
ax.set_ylim(-10, 5)

# Add labels and title
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("Final State")

# Add legend
ax.legend()

# Show the plot
plt.show()


fig, ax = plt.subplots()

# Plot the nodes
for i, node in enumerate(initial_conditions):
    if node[3]:  # Fixed node
        if node[5] == "line":
            ax.scatter(node[0][0], node[0][1], color="purple", marker="o")
        else:
            ax.scatter(node[0][0], node[0][1], color="red", marker="o")
    else:  # Free node
        ax.scatter(node[0][0], node[0][1], color="blue", marker="o")
        
for i, node in enumerate(final_positions):
    if node[3]:  # Fixed node
        if node[4] == "line":
            ax.scatter(node[0][0], node[0][1], color="purple", marker="o")
        else:
            ax.scatter(node[0][0], node[0][1], color="red", marker="o")
    else:  # Free node
        ax.scatter(node[0][0], node[0][1], color="blue", marker="o")
        

# Plot the connections between nodes
for i, connection in enumerate(connections):
    line = np.column_stack(
        [initial_conditions[connection[0]][0][:2], initial_conditions[connection[1]][0][:2]]
    )
    ax.plot(line[0], line[1], color="black",label="Initial state" if i == 0 else "")
    line = np.column_stack(
        [final_positions[connection[0]][0][:2], final_positions[connection[1]][0][:2]]
    )
    ax.plot(line[0], line[1], color="red",label="Final state" if i == 0 else "")
    
plt.legend()
        
    
plt.show()

[PatrickRoel]

pull request #27