m, nA, nB = (

len(myBCs.ground_structure),

myBCs.ground_structure[:, 0].astype(int),

myBCs.ground_structure[:, 1].astype(int),

)

l, X, Y = (

myBCs.ground_structure_length,

myBCs.nodes[nB, 0] - myBCs.nodes[nA, 0],

myBCs.nodes[nB, 1] - myBCs.nodes[nA, 1],

)

d0, d1, d2, d3 = (

myBCs.dofs[nA * 2],

myBCs.dofs[nA * 2 + 1],

myBCs.dofs[nB * 2],

myBCs.dofs[nB * 2 + 1],

) # Defined on boundary conditions. One line for each candidate

s = np.concatenate((-X / l * d0, -Y / l * d1, X / l * d2, Y / l * d3))

r = np.concatenate((nA * 2, nA * 2 + 1, nB * 2, nB * 2 + 1)) # dofs

c = np.concatenate((np.arange(m), np.arange(m), np.arange(m), np.arange(m)))

B = sp.coo_matrix((s, (r, c)), shape=(len(myBCs.dofs), m))

sp.coo_matrix.eliminate_zeros(B) # Eliminate explicit zeros

m, nA, nB = (

len(myBCs.ground_structure),

myBCs.ground_structure[:, 0].astype(int),

myBCs.ground_structure[:, 1].astype(int),

)

l, X, Y, Z = (

myBCs.ground_structure_length,

myBCs.nodes[nB, 0] - myBCs.nodes[nA, 0],

myBCs.nodes[nB, 1] - myBCs.nodes[nA, 1],

myBCs.nodes[nB, 2] - myBCs.nodes[nA, 2],

)

d0, d1, d2, d3, d4, d5 = (

myBCs.dofs[nA * 3],

myBCs.dofs[nA * 3 + 1],

myBCs.dofs[nA * 3 + 2],

myBCs.dofs[nB * 3],

myBCs.dofs[nB * 3 + 1],

myBCs.dofs[nB * 3 + 2],

) # Defined on boundary conditions. One line for each candidate

s = np.concatenate(

(-X / l * d0, -Y / l * d1, -Z / l * d2, X / l * d3, Y / l * d4, Z / l * d5)

)

r = np.concatenate(

(nA * 3, nA * 3 + 1, nA * 3 + 2, nB * 3, nB * 3 + 1, nB * 3 + 2)

) # dofs

c = np.tile(np.arange(m), 6)

B = sp.coo_matrix((s, (r, c)), shape=(len(myBCs.dofs), m))

sp.coo_matrix.eliminate_zeros(B) # Eliminate explicit zeros

BCs: BCS.MBB2D_Symm,

tol,

a,

q=False,

a_cell=None,

cell_mapping_vector=None,

delete_chain=False,

):

"""Use this method to reduce the size of the Bcs

reduced_BCs = copy.deepcopy(BCs)

keep_value = np.max(a) * tol

reduced_BCs.candidates = (

a >= keep_value

) # (q>0 & a>keep_value_t) + (q<0 & a>keep_value_c)

reduced_BCs.candidates_post_unstable_node_suppression = np.where(a >= keep_value)[

0

] # (q>0 & a>keep_value_t) + (q<0 & a>keep_value_c)

a_reduced = a[reduced_BCs.candidates]

# Ground structure (step 1)

reduced_BCs.ground_structure_not_reduced = (

BCs.ground_structure.copy()

) # Used for plotting

reduced_BCs.nodes_not_reduced = BCs.nodes.copy()

reduced_BCs.ground_structure = BCs.ground_structure[reduced_BCs.candidates].copy()

reduced_BCs.ground_structure_length = BCs.ground_structure_length[

reduced_BCs.candidates

].copy()

if a_cell is not None:

### Eliminates buckling chain

if delete_chain:

if BCs.is3D:

q_reduced = q[reduced_BCs.candidates]

# Identify if there is some chain buckling

unique, counts = np.unique(

reduced_BCs.ground_structure, return_counts=True

)

possible_chain_nodes = np.array(unique[np.where(counts == 2)]).ravel()

chain_nodes = []

for k in range(

possible_chain_nodes.shape[0]

): # needed to avoid adding BC nodes

candidates_id = (

reduced_BCs.ground_structure[:, 0] == possible_chain_nodes[k]

) | (reduced_BCs.ground_structure[:, 1] == possible_chain_nodes[k])

candidates = reduced_BCs.ground_structure[candidates_id, :]

diff = np.setxor1d(candidates[0, :], candidates[1, :])

start, end = np.min(diff), np.max(diff)

are_bars_in_compression = False

if np.all(q_reduced[candidates_id] < 0):

are_bars_in_compression = True

# Evaluate cross product

n1 = reduced_BCs.nodes[candidates[0, 0], :]

n2 = reduced_BCs.nodes[candidates[0, 1], :]

n3 = reduced_BCs.nodes[candidates[1, 0], :]

n4 = reduced_BCs.nodes[candidates[1, 1], :]

v1 = n2 - n1

v2 = n4 - n3

if np.allclose(

np.linalg.norm(np.cross(v2, v1), 2), 0, rtol=1e-2

): # cross prod = 0

if (possible_chain_nodes[k] > start) and (

possible_chain_nodes[k] < end

): # If not some BCs can become candidates

if (

are_bars_in_compression

): # We charge only nodes in a compression chain

chain_nodes.append(possible_chain_nodes[k])

merged_candidates = []

reduced_BCs.merged_candidates = []

# Merge chain nodes

for i, node in enumerate(

chain_nodes

): # need to use loops and not a vectorized routine

leng = 0

candidates_to_merge = reduced_BCs.ground_structure[

(reduced_BCs.ground_structure[:, 0] == node)

| (reduced_BCs.ground_structure[:, 1] == node),

:,

]

merged_candidates.append(

sorted([candidates_to_merge[candidates_to_merge != node]])

)

# Drop candidates_to_merge and add merged_candidates

delll = []

for k in range(candidates_to_merge.shape[0]):

dell = int(

np.where(

(

reduced_BCs.ground_structure

== candidates_to_merge[k, :]

).all(axis=1)

)[0]

)

delll.append(dell)

reduced_BCs.merged_candidates.append(delll)

for k in range(candidates_to_merge.shape[0]):

dell = np.where(

(

reduced_BCs.ground_structure

== candidates_to_merge[k, :]

).all(axis=1)

)

reduced_BCs.ground_structure = np.delete(

reduced_BCs.ground_structure, dell, axis=0

)

add_a = a_reduced[dell]

add_q = q_reduced[dell]

a_reduced = np.delete(a_reduced, dell, axis=0)

id_cand = np.nonzero(

reduced_BCs.candidates_post_unstable_node_suppression

)[0]

id_remove = id_cand[dell[0]]

reduced_BCs.candidates_post_unstable_node_suppression = (

np.delete(

reduced_BCs.candidates_post_unstable_node_suppression,

id_remove,

)

)

q_reduced = np.delete(q_reduced, dell, axis=0)

leng += reduced_BCs.ground_structure_length[dell]

reduced_BCs.ground_structure_length = np.delete(

reduced_BCs.ground_structure_length, dell, axis=0

)

add = np.array(merged_candidates[i]).reshape(1, 2)

if (

not np.equal(add, reduced_BCs.ground_structure)

.all(axis=1)

.any()

): # check before adding duplicates

reduced_BCs.ground_structure = np.vstack(

[reduced_BCs.ground_structure, add]

)

a_reduced = np.append(a_reduced, add_a)

reduced_BCs.candidates_post_unstable_node_suppression = (

np.append(

reduced_BCs.candidates_post_unstable_node_suppression,

True,

)

)

try:

q_reduced = np.vstack([q_reduced, add_q])

except:

q_reduced = np.append(q_reduced, add_q)

reduced_BCs.ground_structure_length = np.append(

reduced_BCs.ground_structure_length, leng

)

reduced_BCs.merged_candidates = np.array(reduced_BCs.merged_candidates)

else:

q_reduced = q[reduced_BCs.candidates]

# Identify if there is some chain buckling

unique, counts = np.unique(

reduced_BCs.ground_structure, return_counts=True

)

possible_chain_nodes = np.array(unique[np.where(counts == 2)]).ravel()

chain_nodes = []

for k in range(

possible_chain_nodes.shape[0]

): # needed to avoid adding BC nodes

cand = reduced_BCs.ground_structure[

(reduced_BCs.ground_structure[:, 0] == possible_chain_nodes[k])

| (

reduced_BCs.ground_structure[:, 1]

== possible_chain_nodes[k]

),

:,

]

diff = np.setxor1d(cand[0, :], cand[1, :])

start, end = np.min(diff), np.max(diff)

x0, y0 = reduced_BCs.nodes[cand[0, 0], :]

x1, y1 = reduced_BCs.nodes[cand[0, 1], :]

x2, y2 = reduced_BCs.nodes[cand[1, 0], :]

x3, y3 = reduced_BCs.nodes[cand[1, 1], :]

angle1 = np.arctan2(y1 - y0, x1 - x0)

angle2 = np.arctan2(y3 - y2, x3 - x2)

if np.allclose(np.abs(angle1), np.abs(angle2), rtol=1e-2):

if (possible_chain_nodes[k] > start) and (

possible_chain_nodes[k] < end

): # If not some BCs can become candidates

chain_nodes.append(possible_chain_nodes[k])

merged_candidates = []

reduced_BCs.merged_candidates = []

# Merge chain nodes

for i, node in enumerate(

chain_nodes

): # need to use loops and not a vectorized routine

leng = 0

candidates_to_merge = reduced_BCs.ground_structure[

(reduced_BCs.ground_structure[:, 0] == node)

| (reduced_BCs.ground_structure[:, 1] == node),

:,

]

merged_candidates.append(

sorted([candidates_to_merge[candidates_to_merge != node]])

)

# Drop candidates_to_merge and add merged_candidates

delll = []

for k in range(candidates_to_merge.shape[0]):

dell = int(

np.where(

(

reduced_BCs.ground_structure

== candidates_to_merge[k, :]

).all(axis=1)

)[0]

)

delll.append(dell)

reduced_BCs.merged_candidates.append(delll)

for k in range(candidates_to_merge.shape[0]):

dell = int(

np.where(

(

reduced_BCs.ground_structure

== candidates_to_merge[k, :]

).all(axis=1)

)[0]

)

reduced_BCs.ground_structure = np.delete(

reduced_BCs.ground_structure, dell, axis=0

)

add_a = a_reduced[dell]

add_q = q_reduced[dell]

a_reduced = np.delete(a_reduced, dell, axis=0)

try:

eq = int(

np.where(

(

BCs.ground_structure

== candidates_to_merge[k, :]

).all(axis=1)

)[0]

)

r = int(

np.where(

reduced_BCs.candidates_post_unstable_node_suppression

== eq

)[0]

)

reduced_BCs.candidates_post_unstable_node_suppression = np.delete(

reduced_BCs.candidates_post_unstable_node_suppression, r

)

delete = True

except:

delete = False

q_reduced = np.delete(q_reduced, dell, axis=0)

leng += reduced_BCs.ground_structure_length[dell]

reduced_BCs.ground_structure_length = np.delete(

reduced_BCs.ground_structure_length, dell, axis=0

)

add = np.array(merged_candidates[i]).reshape(1, 2)

if (

not np.equal(add, reduced_BCs.ground_structure)

.all(axis=1)

.any()

): # check before adding duplicates

reduced_BCs.ground_structure = np.vstack(

[reduced_BCs.ground_structure, add]

)

a_reduced = np.append(a_reduced, add_a)

if delete:

reduced_BCs.candidates_post_unstable_node_suppression = np.append(

reduced_BCs.candidates_post_unstable_node_suppression,

eq,

)

try:

q_reduced = np.vstack([q_reduced, add_q])

except:

q_reduced = np.append(q_reduced, add_q)

reduced_BCs.ground_structure_length = np.append(

reduced_BCs.ground_structure_length, leng

)

reduced_BCs.merged_candidates = np.array(reduced_BCs.merged_candidates)

# Nodes

# Identify the nodes to be removed from the array

reduced_BCs.nodes_candidates = reduced_BCs.ground_structure[:, [0, 1]].copy()

reduced_BCs.nodes_candidates = reduced_BCs.nodes_candidates.ravel().tolist()

# Sort and cancel duplicates

reduced_BCs.nodes_candidates = np.array(sorted(set(reduced_BCs.nodes_candidates)))

reduced_BCs.nodes = BCs.nodes[reduced_BCs.nodes_candidates, :].copy()

# Fixed DOFs

if BCs.is3D:

reduced_BCs.dofs_candidates = reduced_BCs.nodes_candidates.copy()

reduced_BCs.dofs_candidates = np.append(

reduced_BCs.dofs_candidates * 3, reduced_BCs.nodes_candidates * 3 + 1

)

reduced_BCs.dofs_candidates = np.sort(

np.append(reduced_BCs.dofs_candidates, reduced_BCs.nodes_candidates * 3 + 2)

)

else:

reduced_BCs.dofs_candidates = reduced_BCs.nodes_candidates.copy()

reduced_BCs.dofs_candidates = np.sort(

np.append(

reduced_BCs.dofs_candidates * 2, reduced_BCs.nodes_candidates * 2 + 1

)

)

# We need to translate the ID of the full BC to the reduced

dofs_list = np.arange(reduced_BCs.dofs_candidates.size)

reduced_BCs.dofs = BCs.dofs[reduced_BCs.dofs_candidates].copy()

reduced_BCs.fixed_dofs = np.array(np.where(reduced_BCs.dofs == 0)).ravel()

reduced_BCs.free_dofs = np.setdiff1d(dofs_list, reduced_BCs.fixed_dofs)

if BCs.is3D:

reduced_BCs.dofs_nodes_connectivity = reduced_BCs.dofs.reshape((-1, 3))

else:

reduced_BCs.dofs_nodes_connectivity = reduced_BCs.dofs.reshape((-1, 2))

# Force

tmp = np.zeros(BCs.dofs.shape[0], dtype="int")

tmp[BCs.force_dofs] = 1

tmp = tmp[reduced_BCs.dofs_candidates]

reduced_BCs.force_dofs = np.array(np.where(tmp == 1)).ravel()

reduced_BCs.R = BCs.R[reduced_BCs.dofs_candidates].copy()

reduced_BCs.R_free = reduced_BCs.R[reduced_BCs.free_dofs].copy()

# Ground structure (step 2)

# The node ID must be replaced as the nodes matrix has changed

replace_matrix = np.vstack(

[reduced_BCs.nodes_candidates, np.arange(reduced_BCs.nodes_candidates.size)]

).T # find (left col) and replace (right col) matrix used for the select conditions

condlist = replace_matrix[:, 0][:, None, None] == reduced_BCs.ground_structure

reduced_BCs.ground_structure_sel = np.select(condlist, replace_matrix[:, 1])

reduced_BCs.ground_structure[:, :2] = reduced_BCs.ground_structure_sel[:, :2]

if cell_mapping_vector is not None:

n_topologies_cell = len(set(cell_mapping_vector.ravel().tolist()))

else:

n_topologies_cell = 1

if a_cell is not None:

a_cell_reduced, reduced_BCs.keep_cell = reduce_A(a_cell, keep_value)

reduced_BCs.N_cell = np.zeros(n_topologies_cell, dtype="int")

for i in range(n_topologies_cell):

reduced_BCs.N_cell[i] = a_cell_reduced[i].size

reduced_BCs.N_cell_max = np.max(reduced_BCs.N_cell)

# add -1 values to the sections of the cells. This permits the Kronecher product

a_cell_reduced_pad = a_cell_reduced.copy()

for i in range(a_cell_reduced.size):

a_cell_reduced_pad[i] = np.pad(

a_cell_reduced_pad[i],

pad_width=(0, reduced_BCs.N_cell_max - a_cell_reduced_pad[i].size),

constant_values=-1,

)

# Create a bool matrix used to understand which sections are active

reduced_BCs.reduction_pattern = (a_cell_reduced_pad[0] >= 0).reshape(

1, -1

) # this array should be at least 2D

for i in range(a_cell_reduced.size - 1):

reduced_BCs.reduction_pattern = np.vstack(

(reduced_BCs.reduction_pattern, a_cell_reduced_pad[i + 1] >= 0)

)

# Create a bool matrix to calculate the lenght of every section

reduced_BCs.reduction_pattern_len = (a_cell[0] > keep_value).reshape(1, -1)

for i in range(a_cell_reduced.size - 1):

reduced_BCs.reduction_pattern_len = np.vstack(

(reduced_BCs.reduction_pattern_len, a_cell[i + 1] > keep_value)

)

# Create a matrix with the unreduced legnht of the members of every cell topology

reduced_BCs.l_cell_full = np.tile(

BCs.ground_structure_length[: a_cell[0].shape[0]], a_cell.shape[0]

).reshape((-1, a_cell[0].shape[0]))

# Update the list of the ground structure of the reduced problem

reduced_BCs.ground_stucture_list_cell = []

for i in range(n_topologies_cell):

reduced_BCs.ground_stucture_list_cell.append(

np.array(BCs.ground_stucture_list_cell)[a_cell[i] > keep_value].tolist()

)

reduced_BCs.ground_stucture_list_not_reduced = BCs.ground_structure[

:, [0, 1]

].tolist()

reduced_BCs.ground_stucture_list = reduced_BCs.ground_structure.tolist()

reduced_BCs.isReduced = True

if reduced_BCs.isFreeForm:

p = np.array(np.where(reduced_BCs.ground_structure[:, 2] == 1)).ravel()

a = np.array(np.where(reduced_BCs.ground_structure[:, 2] != 1)).ravel()

reduced_BCs.ground_structure_cellular = reduced_BCs.ground_structure[

p, :

].copy()

reduced_BCs.ground_structure_aperiodic = reduced_BCs.ground_structure[

a, :

].copy()

reduced_BCs.N = len(

reduced_BCs.ground_structure

) # Number of member of the Ground Structure

reduced_BCs.M = len(reduced_BCs.dofs) # Number of DOFs of the Groud Structure

reduced_BCs.N_cellular = len(reduced_BCs.ground_structure_cellular)

reduced_BCs.N_aperiodic = len(reduced_BCs.ground_structure_aperiodic)

reduced_BCs.N_cell = int(reduced_BCs.N_cell)

if a_cell is not None:

reduced_BCs.ground_structure_length_cell = BCs.ground_structure_length_cell[

keep_cell

].copy()

return reduced_BCs, a_reduced

if a_cell is not None:

return reduced_BCs, a_cell_reduced_pad

else:

a_reduced = a.copy()

keep = np.zeros(a_reduced.size, dtype="object_")

for i in range(a.size):

keep[i] = a_reduced[i] > thresold

a_reduced[i] = a_reduced[i][keep[i]]

# VL parameter

if np.any(cell_mapping_vector) == False:

cell_mapping_vector = np.ones(BC.ncel_str, dtype="int").ravel()

topologies_ID = set(cell_mapping_vector.ravel().tolist()) # unique

n_topologies_cell = len(topologies_ID) # How many different topologies are present

cell_mapping_matrix = np.zeros((BC.ncel_str, n_topologies_cell), dtype="int")

for i, item in enumerate(topologies_ID): # starting from the lower left corner

cell_mapping_matrix[np.where(cell_mapping_vector == item), i] += (

1 # create matrix to correctly assign the different cell topologies according to cell_mapping_matrix

)

if not BC.isReduced:

if np.any(n_topologies_cell) != False:

BC.N_cell = np.ones(n_topologies_cell, dtype="int") * BC.N_cell

NN = np.sum(

BC.N_cell

) # Number of elements in the cell times the number of different cells (area design variables)

myBCs, myBCs_reduced, N_design_var, area_id, force_id, U_id, E, a_init, a_fem

):

"""Method used to calculate the starting point of the design variables (x0) for a cellular case with or withour clusters and preinit"""

if myBCs.is3D:

B_old = calcB_3D(myBCs) # equilibrium matrix

else:

B_old = calcB(myBCs) # equilibrium matrix

M_old = myBCs.M # Number of DOFs of the Groud Structure

## Starting point

x0 = np.zeros(N_design_var) # Define an initial guess for the optimization

# Areas

x0[area_id] = a_init

area_phys = a_fem

area_phys[area_phys < np.max(a_init) * 1e-8] = np.max(a_init) * 1e-8

# Initial forces and displacements are calculated using FEM

U = np.zeros(M_old)

K = calcK(myBCs, B_old, E, area_phys)

keep = myBCs.free_dofs

K = K[keep, :][:, keep]

U[keep] = sp.linalg.spsolve(K, myBCs.R_free) # FEM analysis linear sistem

# Forces

F_old = area_phys * E / myBCs.ground_structure_length * (B_old.T @ U)

x0[force_id] = F_old[myBCs_reduced.candidates_post_unstable_node_suppression]

# Displacements

x0[U_id] = U[myBCs_reduced.dofs_candidates]

myBCs,

myBCs_reduced,

N_design_var,

area_id,

force_id,

U_id,

N,

M,

B,

E,

n_topologies_cell,

N_cell,

cell_mapping_matrix,

a_cell,

):

"""Method used to calculate the starting point of the design variables (x0) for a cellular case with or withour clusters and preinit"""

if myBCs.is3D:

B_old = calcB_3D(myBCs) # equilibrium matrix

else:

B_old = calcB(myBCs) # equilibrium matrix

M_old = myBCs.M # Number of DOFs of the Groud Structure

## Starting point

x0 = np.zeros(N_design_var) # Define an initial guess for the optimization

# Areas

a_cell_list = []

for i in range(n_topologies_cell):

a_cell_list = np.append(

a_cell_list, a_cell[i][a_cell[i] > -1e-5]

) # Eliminate -1 values used for the Kronecher product, a little neg value is accepted for numeric errors

a_cell

x0[area_id] = a_cell_list

area_phys = np.zeros(myBCs.ncel_str * myBCs.N_cell[0])

for i in range(n_topologies_cell): # starting from the lower left corner

a_cell_pad = np.ones(myBCs.N_cell[0]) * np.max(a_cell[0]) * 1e-8

a_cell_pad[myBCs_reduced.keep_cell[i]] = a_cell[i][

range(myBCs_reduced.N_cell[i])

]

area_phys += np.kron(cell_mapping_matrix[:, i], a_cell_pad)

# Initial forces and displacements are calculated using FEM

U = np.zeros(M_old)

K = calcK(myBCs, B_old, E, area_phys)

keep = myBCs.free_dofs

K = K[keep, :][:, keep]

U[keep] = sp.linalg.spsolve(K, myBCs.R_free) # FEM analysis linear sistem

# Forces

F_old = area_phys * E / myBCs.ground_structure_length * (B_old.T @ U)

x0[force_id] = F_old[myBCs_reduced.candidates]

# Displacements

x0[U_id] = U[myBCs_reduced.dofs_candidates]