now including lower dirs

import sys

import numpy as np

import mechnet as mn

from multiprocessing import Pool

def get_parametercollection(seed):

    par = mn.ParameterCollection()

    nx = 8 # needs to == el**h

    ny = nx

    nz = 9 # need to be >=  h +2 + o_set (+2 bc bound cond)    

    n = nx*ny*nz

    par.set_N(n)

    par.set_Nx(nx)

    par.set_Ny(ny)

    par.set_Nz(nz)

    par.set_NxNy()

    #new parameters necessary for the structure

    par.set_xhierarchicalelementsize(2) # log_el(nx) = h <-> nx = el**h 

    par.set_yhierarchicalelementsize(2)

    par.set_gapzoffset(1) # == o_set

    par.set_thresholdrng(seed)

    par.set_structurerng(seed + 1_000_000)

    par.set_weibullparameter(1.5)

    par.set_scalingfactor(1.00)

    par.set_uniformupperdirichlet(1.0)

    return par

def run_a_simulation_fuse(seed):

    constructor = mn.cubic3D.Cubic3DFullConnectionConstructor()

    #new decorators to modify network structure:

    #random positioning of missing edges fitting to the amount that would be missing in a hierarchical network with

    #element sizes as determined by set_xhierarchicalelementsize() and set_yhierarchicalelementsize()

    #in the corresponding directions; arangement dependent on seed given by set_structurerng()

    #for the random positioning

    constructor = mn.cubic3D.DensityReferenceStructureDecorator(constructor)

    constructor = mn.cubic3D.ScalingWeibullThresholdDecorator(constructor)

    boundaries = mn.general.EmptyBoundariesConstructor()

    boundaries = mn.cubic3D.FixedLowerBoundaryDecorator(boundaries)

    boundaries = mn.cubic3D.UniformDisplacedUpperBoundaryDecorator(boundaries)

    parametercollection = get_parametercollection(seed)

    network = constructor.start_construction(parametercollection)

    bounded_network = boundaries.start_assign_boundaries(network)

    #switch to fuse simulation to make it easier to read the resulting stiffness matrix

    dict_of_simulation_parameters = {"niter":10_000_000, "externalV":1.0, "accuracy":12}

    builder = mn.sim.FuseDirichletBuilder()

    solver = mn.sim.PARDISOBindingsSolver()

    intrescal = mn.sim.StartInterimResultsCalculator() #InterimResultsCalculator stack pre-computes data needed for output and simulation step

    intrescal = mn.sim.FuseCurrentsVoltagesCalculator(intrescal)

    outgen = mn.sim.StartOutputGenerator()

    applier = mn.sim.HottestFuseBreakApplier()

    checker = mn.sim.NonInitialChecker(mn.sim.ConnectedPathChecker())

    simulation = mn.sim.Simulation(builder, solver, intrescal, outgen, applier, checker)

    dict_of_crs_edge_data = bounded_network.get_network_data() 

    simdata = simulation.start_simulation(dict_of_crs_edge_data, dict_of_simulation_parameters)

    simdata_dict = simdata.get_simdata_dict()

    network_descriptor_dict = bounded_network.get_construction_data()

    mn.datman.save_simulation_output(network_descriptor_dict, simdata_dict, simulation.simulationname, simulation.simulationdoc, suffix=str(seed), overwritemode = False)

    print("done with sim of seed", seed)

    return (network_descriptor_dict, simdata_dict, simulation.simulationname, simulation.simulationdoc)

if __name__=="__main__":

    run_a_simulation_fuse(1001)

    list_of_parameters =  []

    for seed in range(1000, 1010):

        list_of_parameters.append(seed) 

    with Pool(10) as p:

        results = p.map(run_a_simulation_fuse, list_of_parameters)

 #   for data in results:

 #       mn.datman.save_simulation_output(data[0], data[1], data[2], data[3], targetdirectory="./", suffix="", overwritemode=False)

    #print(data[1]["current_data"]["stiffnessmatrix"].toarray())

    print("\nRun successfull.")

 No newline at end of file

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import numpy as np

def x_position(node, nx, ny, nz):

    return (node%nx)

def y_position(node, nx, ny, nz):

    return ((node//nx)%ny)

def z_position(node,nx, ny, nz):

    return node//(nx*ny)

def get_edges_array(network_dict : dict, nxnynz : [int,int,int]):

    nx,ny,nz = nxnynz

    edges_array = []

    for key in network_dict.keys():

       node_1, node_2 = [int(key_string) for key_string in key.split('_')]

       node_1_coordinates = [x_position(node_1, nx, ny, nz), y_position(node_1, nx, ny, nz), z_position(node_1,  nx, ny, nz)] 

       node_2_coordinates = [x_position(node_2, nx, ny, nz), y_position(node_2, nx, ny, nz), z_position(node_2,  nx, ny, nz)] 

       line_entry = [*node_1_coordinates,*node_2_coordinates]

       edges_array.append(line_entry)

    return edges_array

def plot_edges_mechnet_network(network_dict : dict, nxnynz : [int,int,int], filename ="plot") -> None: 

    plt.rcParams.update({

    'font.size': 15,              # Default font size for all text

    'axes.labelsize': 15,         # Font size for axis labels

    'xtick.labelsize': 15,        # Font size for x-axis tick labels

    'ytick.labelsize': 15,        # Font size for y-axis tick labels

    'legend.fontsize': 15         # Font size for legend text

    })

    # Example data: Edges with values

    #edges = [

    #    (0, 0, 0, 1, 0, 0, 1.0),

    #    (1, 0, 0, 1, 1, 0, 1.0),

    #    (1, 1, 0, 0, 1, 0, 1.0),

    #    (0, 1, 0, 0, 1, 1, 4.0),

    #]

    edges = get_edges_array(network_dict, nxnynz)

    fig = plt.figure()

    ax = fig.add_subplot(111, projection='3d')

    # Plot each edge

    for edge in edges:

        x = [edge[0], edge[3]]

        y = [edge[1], edge[4]]

        z = [edge[2], edge[5]]

        ax.plot(x, y, z, color = 'r', linewidth=2)

    ax.set_xlabel(r'$L_{x}$')

    ax.set_xlim(0,16)

    ax.set_ylabel(r'$L_{y}$')

    ax.set_ylim(-8,8)

    ax.set_zlabel(r'$L_{z}$')

    ax.set_zlim(0,16)

    fig.savefig(f'{filename}.png')

    # Add a color bar

    plt.show()

 No newline at end of file

import mechnet as mn

import my_plot_edges

from my_plot_edges import x_position, z_position, y_position

NX = 16

NY = 2

NZ = 16

def get_parametercollection(seed):

    par = mn.ParameterCollection()

    nx = NX # needs to == el**h

    ny = NY

    nz = NZ # need to be >=  h +2 + o_set (+2 bc bound cond)    

    n = nx*ny*nz

    par.set_N(n)

    par.set_Nx(nx)

    par.set_Ny(ny)

    par.set_Nz(nz)

    par.set_NxNy()

    #new parameters necessary for the structure

    #par.set_xhierarchicalelementsize(2) # log_el(nx) = h <-> nx = el**h 

    #par.set_yhierarchicalelementsize(2)

    #par.set_gapzoffset(1) # == o_set

    par.set_thresholdrng(seed)

    par.set_structurerng(seed +1_000_000)

    par.set_weibullparameter(1.5)

    par.set_scalingfactor(1.00)

    par.set_uniformupperdirichlet(1.0)

    return par

def construct_system():

    constructor = mn.cubic3D.Cubic3DFullConnectionConstructor()

    #constructor = mn.cubic3D.ShuffledHierarchicalGapsDecorator(constructor)

    constructor = mn.cubic3D.ScalingWeibullThresholdDecorator(constructor)

    boundaries = mn.general.EmptyBoundariesConstructor()

    boundaries = mn.cubic3D.FixedLowerBoundaryDecorator(boundaries)

    boundaries = mn.cubic3D.UniformDisplacedUpperBoundaryDecorator(boundaries)

    parametercollection = get_parametercollection(1001)

    network = constructor.start_construction(parametercollection)

    net_dict = network.edge_dict

    new_net_dict = cut_heart(net_dict)

    network.edge_dict = new_net_dict

    return network

def cut_heart(net_dict : dict) -> dict:

    new_dict = {}

    nx, ny, nz = NX, NY, NZ

    for key, value in net_dict.items():

        node_1, node_2 = [int(key_string) for key_string in key.split('_')]

        node_1_coordinates = [x_position(node_1, nx, ny, nz), y_position(node_1, nx, ny, nz), z_position(node_1,  nx, ny, nz)] 

        node_2_coordinates = [x_position(node_2, nx, ny, nz), y_position(node_2, nx, ny, nz), z_position(node_2,  nx, ny, nz)] 

        if is_in_shape(x=node_1_coordinates[0], z=node_1_coordinates[2]) and is_in_shape(x=node_2_coordinates[0], z=node_2_coordinates[2]):

            new_dict[key] = value

    return new_dict

def is_in_shape(x, z):

    if x <= NX*0.5:

        if z > -x +(NX*0.5) and z < -x + (NX*1.25) : 

            ans = True

        else:

            ans =  False

    elif x > NX*0.5:

        if z > x - (NX*0.5) and z < x + (NX*0.25): 

            ans = True

        else: 

            ans = False

    return ans

if __name__ == "__main__":

    network = construct_system()

    edge_dict = network.edge_dict

    filename_plot = f"visualize_heart_{NX}x{NY}x{NZ}"

    my_plot_edges.plot_edges_mechnet_network(edge_dict,  [NX,NY,NZ], filename=filename_plot)

    print("cutoff prev")

 No newline at end of file

https://www.dropbox.com/scl/fo/d5by2bl05kj8bzti8n3kq/h?rlkey=7zr0v86xqz8mutigb6uogiwps&st=qmx4ndtw&dl=0

 No newline at end of file