<jupyter_start><jupyter_text>VectoRelax IntroductionIn today's session we will focus on more advanced uses of scientific and computing libraries in Python. The main exercise will allow you to deform a mesh into a relaxed position, given an inverted gravitational force.  Aknowledgements*Dynamic relaxation* exercise was originally written by <NAME> and <NAME>, *NumPy/SciPy introduction* by <NAME>. Last edit on 09/2020Thanks to Ir. <NAME> for introducing the following reference:* <NAME>., <NAME>, <NAME>, and <NAME> (2014). Shell Structures for Architecture: Form Finding and Optimization. New York: Routledge. * Chapter 2: Review of Dynamic Relaxation with an extension to six degrees of freedom theory by <NAME>.  1. NumPy[Numpy](https://numpy.org/doc/1.17/reference/index.html) ([the user guide](https://numpy.org/doc/1.17/user/index.htmluser)) is the core library for scientific computing in Python. It provides very fast code for linear algebra operations, effectively replacing Matlab. Numpy is a **core Python library** upon which many others are built on (e.g. Pandas for data engineering, SciKit-learn for machine learning). So knowing the basics of Numpy will help you to faster understand the higher-level libraries and adopt solutions (look up Stack Overflow: as long as you know how to Google the solution will be there...) How does Numpy work?The base component in Numpy is an array. In short it works like an algebraic matrix: it can have one or more dimensions. A numpy array is a grid of values, **all of the same type**. It smartly stores data in one place in the computer memory, unlike other Python array-like objects. Those are scattered in memory and the code just has links to those locations. It's easy to use, but performance suffers a lot. SciPy : TODOSciPy is an open-source Python library used for scientific and technical computing. It has several sub-packages, such as on covering constants, spatial algorithms (e.g. distance, nearest neighbor functions), linear algebra methods and many more. It is built on top of NumPy, making it not only a useful, but also a powerful tool. More about it [here](https://www.tutorialspoint.com/scipy/index.htm).This library will not be used in the following exercise, but it will prove very useful if you need to do more complex linear algebra routines, such as integration. Let's create a NumPy array<jupyter_code>import numpy as np

# Create a rank 1 array (1D matrix)
a = np.array([1, 2, 3])   

print(a)        
print(type(a))            
print(a.shape)            
print(a[0], a[1], a[2]) 
a[0] = 5
print(a)<jupyter_output>[1 2 3]
<class 'numpy.ndarray'>
(3,)
1 2 3
[5 2 3]<jupyter_text>Arrays can be made in many other ways and remembering this functionality (or printing out a cheat sheet) will save you precious coding time. Also, you are **not** limited to 2 dimensions: all examples below can be extended by 1 dimension by adding an additional value after the comma.<jupyter_code># Creates an array of all zeros
a = np.zeros((2,2))  
print(a)              
                      
# Creates an array of all ones
b = np.ones((2,3))    
print('\n',b)              

# Creates an array of all of a specific number (5 in this case)
c = np.full((2,2), 5)  
print('\n',c)               
                       
# Creates an identity matrix of a specified size (3 in this case)
d = np.eye(3)         
print('\n',d)              
                      
# Creates a 2x3 array filled with random values in [0,5)
e = np.random.random((2,5))  
print('\n',e)<jupyter_output>[[0. 0.]
 [0. 0.]]

 [[1. 1. 1.]
 [1. 1. 1.]]

 [[5 5]
 [5 5]]

 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

 [[0.42040034 0.18807518 0.91282968 0.32322076 0.99204222]
 [0.00136697 0.49202091 0.80188799 0.84731008 0.10259576]]<jupyter_text>Exercise 1Create an identity matrix of size 3, then change all the ones to a different value each.<jupyter_code>array = np.eye(3)

array[0,0]=2
array[1,1]=4
array[2,2]=18

array<jupyter_output><empty_output><jupyter_text>Performance exampleTo compare the performance of pure Python and Numpy, follow the pseudocode instructions in the next cells.<jupyter_code># Leave it empty and use as example

# create a list, containing 1000**2 lists of [2,2,2]

my_list = []
for x in range(1000**2):
    my_list.append([2,2,2])
def python_test(my_list):
    # multiply each number by 3
    for i in range(len(my_list)):
        for x in range(len(my_list[i])):
            my_list[i][x]=my_list[i][x]*3

%time python_test(my_list)        
# create an array of 3x1000**2 containing only 2s
array = np.full((3,1000**2), 2) 
def numpy_test(array):
    # multiply each number in the array by 3
    array = array * 3
    
%time numpy_test(array)<jupyter_output>CPU times: user 8.11 ms, sys: 7.36 ms, total: 15.5 ms
Wall time: 13.8 ms<jupyter_text>2. Dynamic RelaxationThis exercise consists of 3 main steps:1. Reading in the data and defining the initial state of the mesh2. Relaxing the mesh3. Plotting the results<jupyter_code>from numpy import array, zeros, float64, append, sqrt, sign<jupyter_output><empty_output><jupyter_text>1. Reading in the data and defining the initial state of the meshDuring this step we need to read in the data, parse it to something we can manipulate, pre-process the geometeries and pass the results/ initial state of the mesh to the next process. The expected outputs of this step are:* points : List of points* edges : List of edges* vertices : List of vertices* xyzIS : Initial state of the system* bbmax, bbmin : Bounding box of the mesh The only **input** for the first step of the exercise in the location of the *mesh file (.obj)* you received together with this notebook.To begin with, let's locate the file and open it:<jupyter_code>filename = '../data/basic_dense_grid.obj'
#filename=filename.replace('\','//')
# opening the file
def open_obj_file(filename):
    with open(filename, 'r') as myfile:
        return myfile.read().split('\n')<jupyter_output><empty_output><jupyter_text>What is inside of data? (for a more elaborate description of the format, check out this [link](http://paulbourke.net/dataformats/obj/))<jupyter_code>data = open_obj_file(filename)
data<jupyter_output><empty_output><jupyter_text>How do we make this readable?<jupyter_code>def process_obj_file(data):
    # initiating the output
    xyz = []  # xyz coordinates
    points = []  # references to vertices
    lines = []  # pairs of references to vertices
    
    # parsing the obj file, filling the data lists
    for line in data:
        # split the string line by ' '
        parts = line.split(' ')

        # TODO: this seems unnecessary
        # checking if it contains data 
        if not parts:
            continue

        # assigning the 1st element to head and the rest to tail 
        head = parts[0]
        tail = parts[1:]

        # appending vertex cordinations in (xyz) and assigning an index to them in (points)
        if head == 'v':
            ftail = [float(x) for x in tail]
            xyz.append(ftail)
            points.append(len(xyz)-1)
        # iterating through the edges of the faces and appending them in (lines)
        elif head == 'f':
            ftail = [float(x) for x in tail]
            for i in range(len(ftail)):
                sp = ftail[i%len(ftail)]
                ep = ftail[(i+1)%len(ftail)]
                lines.append((int(sp)-1, int(ep)-1))
    return xyz, points, lines
xyz, points, lines = process_obj_file(data)
xyz<jupyter_output><empty_output><jupyter_text>How do we preprocess the data?<jupyter_code>def mesh_welding(xyz):
    
    # This is how we store the welded points
    vertices = []
    
    # initiating the bounding box max and min point
    bbmax = [0,0,0]
    bbmin = [0,0,0]

    # points to vertices (welding similar vertices)
    x2v = {}
    tol = 0.001 ** 2

    for i, x in enumerate(xyz):
        found = False
        for j, v in enumerate(vertices):
            if (x[0] - v[0]) ** 2 < tol \
                    and (x[1] - v[1]) ** 2 < tol \
                    and (x[2] - v[2]) ** 2 < tol:
                found = True
                x2v[i] = j
                break
        if not found:
            x2v[i] = len(vertices)
            vertices.append(x)

        #finding the bounding box
        for i in range(3):
            if x[i] > bbmax[i]: 
                bbmax[i] = x[i]
            if x[i] < bbmin[i]: 
                bbmin[i] = x[i]
    return bbmin, bbmax, vertices, x2v
bbmin, bbmax, vertices, x2v = mesh_welding(xyz)
vertices<jupyter_output><empty_output><jupyter_text>Anything we forgot?<jupyter_code># redefining the indexes
points[:] = [x2v[index] for index in points]

# redefining the edges by the refined indexes
edges = [(x2v[u[0]], x2v[u[1]]) for u in lines]

# reserving the Initial State of the system
import copy
xyzIS = copy.deepcopy(vertices)<jupyter_output><empty_output><jupyter_text>How does the initial state look like?<jupyter_code>xyzIS == vertices<jupyter_output><empty_output><jupyter_text>Simulation (DR) Input data:* vertices: a list of unique xyz coordinates* edges: a list of pairs of vertex indices* points: index list of fixed vertices* bbmax, bbmin : Bounding box coordinate list Output data:* xyzFS : Final state of the system<jupyter_code># fixed points here are set to be the boundary of the rectangle
fixed = []
for i, pnt in enumerate(iter(xyz)):
    # check if the points are no closer to 
    if abs(pnt[0] - bbmin[0]) < 0.01  \
        or abs(pnt[0] - bbmax[0]) < 0.01 \
        or abs(pnt[1] - bbmin[1]) < 0.01 \
        or abs(pnt[1] - bbmax[1]) < 0.01:
        fixed.append(i)
no_vert = len(vertices) 
no_ed = len(edges)

freePoints = list(set(range(no_vert)) - set(fixed))

xyz = array(vertices)

# p : Applied Load Component which in our case is zero
p = zeros((no_vert, 3), dtype=float64)

# W : Gravity Force
W = array([array([0,0,9.8]) for i in range(no_vert)])

# v : Velosity
v = zeros((no_vert, 3), dtype=float64)


# restlength : Rest Length of the springs
diff = array([(xyz[edges[i][1]] - xyz[edges[i][0]]) for i in range(no_ed)])

restlength = array([sqrt(diff[i,0]*diff[i,0] + diff[i,1]*diff[i,1] + diff[i,2]*diff[i,2]) for i in range(no_ed)])

# adamp : Acceleration damp controls the stability of the system
adamp = 10
IterationMax = 100

# K : Elasticity Constant in Hooke's law
K = 20

# dt : Time intervals
dt = 0.1
mass = 1
for k in range(IterationMax):
    
    # S : Stiffness Force
    S = zeros((no_vert, 3), dtype=float64)
    # R : Final Forces
    R = zeros((no_vert, 3), dtype=float64)
        
    # xyz0 : Initial Position
    xyz0 = xyz.copy()
    # vp : Initial Velocity
    vp = v.copy()
    
    # calculating the length of the edges
    di = array([(xyz[edges[i][1]] - xyz[edges[i][0]]) for i in range(no_ed)])
    dist = array([sqrt(di[i,0]*di[i,0] + di[i,1]*di[i,1] + di[i,2]*di[i,2]) for i in range(no_ed)])

    # Hooke's law : Force = length difference * Elasticity Constant (K)
    EdgeForce = (dist - restlength) * K
    # Edge Force : Decomposing the existing force in each edge into the 3 dimensions 
    EdgeForceS = zeros((no_ed, 3), dtype=float64)

    # iterating throughout the edges
    for i in range(no_ed):
        for j in range(3):
            # Separating the forces into X, Y, Z dimensions
            EdgeForceS[i,j] = (di[i,j] / dist[i]) * EdgeForce[i]

        # Adding the force of each edge to the coresponding nodes
        S[edges[i][0]] -= EdgeForceS[i]
        S[edges[i][1]] += EdgeForceS[i]
    
    R[freePoints] = p[freePoints] + W[freePoints] - S[freePoints]

    A = 1/(1 + adamp * dt * 0.5)
    B = (1 - adamp * dt * 0.5 ) / (1 + adamp * dt * 0.5)
    
    # updating the velocity
    v[freePoints] = A * dt * R[freePoints] / mass + B * vp[freePoints]
    # updating the position
    xyz[freePoints] = xyz0[freePoints] + dt * v[freePoints]


xyzFS = xyz<jupyter_output><empty_output><jupyter_text>3. Plotter<jupyter_code>import pyvista as pv
from compas.datastructures import Mesh as CompasMesh
import compas
import earthy_bridge as eb
basic_dense_grid_filepath = "../data/basic_dense_grid.obj"
basic_dense_grid_mesh = CompasMesh.from_obj(basic_dense_grid_filepath)
_, F = basic_dense_grid_mesh.to_vertices_and_faces()
relaxed_mesh = CompasMesh.from_vertices_and_faces(xyzFS, F)
eb.compas_to_pyvista(relaxed_mesh).plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10931174089068826, 'fraction_numerical': 0.028906882591093117, 'mean_word_length': 3.2443298969072165, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 61, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '305340', 'license': 'MIT', 'n_tokens_mistral': 4289, 'n_tokens_neox': 3960, 'n_words': 1445, 'path': 'notebooks/VectoRelax.ipynb', 'repo_name': 'shervinazadi/earthy_workshops', 'size': '12420'}