SolidsPy scripts

block_uniform load.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Block with different loading and constraint conditions.

This examples runs with the developing version of SolidsPy (2.0).

@author: Nicolás Guarín-Zapata
@date: July 2020
"""
from time import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import spsolve
from solidspy.solids_GUI import solids_auto
from solidspy.preprocesor import rect_grid
from solidspy.assemutil import assembler, DME, loadasem
from solidspy.postprocesor import complete_disp, plot_node_field



# Setup
nx = 60
length = 30
width = 4
x, y, els = rect_grid(length, width, nx, nx//10)
nels = els.shape[0]
npts = x.shape[0]
nodes = np.zeros((npts, 3))
nodes[:, 1] = x
nodes[:, 2] = y

# Loads 
loads = np.zeros((npts, 3))
loads[:, 0] = range(npts)
loads[y == 2, 2] = -1.0
loads[(x == -15) * (y == 2), 2] = -0.5
loads[(x == 15) * (y == 2), 2] = -0.5

# Constraints 1
cons1 = np.zeros((npts, 2), dtype=int)
cons1[np.abs(x + 15) < 1e-6, 0:2] = -1

# Constraints 2
cons2 = np.zeros((npts, 2), dtype=int)
cons2[np.abs(x + 15) < 1e-6, 0:2] = -1
cons2[np.abs(x - 15) < 1e-6, 0:2] = -1

# Constraints 3
cons3 = np.zeros((npts, 2), dtype=int)
cons3[np.abs(y + 2) < 1e-6, 1] = -1
cons3[(np.abs(y + 2) < 1e-6) * (np.abs(x) < 1e-6), 0] = -1


# Materials
mats = np.array([[3e4, 1/3]])

# Data
data1 = {"nodes": nodes,
        "cons": cons1,
        "elements": els,
        "mats": mats,
        "loads": loads}

data2 = {"nodes": nodes,
        "cons": cons2,
        "elements": els,
        "mats": mats,
        "loads": loads}

data3 = {"nodes": nodes,
        "cons": cons3,
        "elements": els,
        "mats": mats,
        "loads": loads}

#%% Solution
# Run the simulation 1
disp1 = solids_auto(data1, compute_strains=False, plot_contours=False)
disp1_mag = np.linalg.norm(disp1, axis=1)

# Run the simulation 2
disp2 = solids_auto(data2, compute_strains=False, plot_contours=False)
disp2_mag = np.linalg.norm(disp2, axis=1)

# Run the simulation 2
disp3 = solids_auto(data3, compute_strains=False, plot_contours=False)
disp3_mag = np.linalg.norm(disp3, axis=1)


#%% Visualization
shape = nx//10 + 1, nx + 1
x.shape = shape
y.shape = shape

plt.subplot(3, 1, 1)
plt.contourf(x + np.reshape(disp1[:, 0], shape),
             y + np.reshape(disp1[:, 1], shape),
             np.reshape(disp1_mag, shape), cmap="magma", zorder=5,
             alpha=0.9)
plt.axis("image")
plt.grid()

plt.subplot(3, 1, 2)
plt.contourf(x + 10*np.reshape(disp2[:, 0], shape),
             y + 10*np.reshape(disp2[:, 1], shape),
             np.reshape(disp2_mag, shape), cmap="magma", zorder=5,
             alpha=0.9)
plt.axis("image")
plt.grid()

plt.subplot(3, 1, 3)
plt.contourf(x + 100*np.reshape(disp3[:, 0], shape),
             y + 100*np.reshape(disp3[:, 1], shape),
             np.reshape(disp3_mag, shape), cmap="magma", zorder=5,
             alpha=0.9)
plt.axis("image")
plt.ylim(-2, 2)
plt.grid()

# %%

solidspy_example.py

import numpy as np
from solidspy.solids_GUI import solids_auto
import matplotlib.pyplot as plt

### Define the data
nodes = np.array([
    [0, 0.00, 0.00],
    [1, 2.00, 0.00],
    [2, 2.00, 2.00],
    [3, 0.00, 2.00],
    [4, 1.00, 0.00],
    [5, 2.00, 1.00],
    [6, 1.00, 2.00],
    [7, 0.00, 1.00],
    [8, 1.00, 1.00]])

cons = np.array([
    [0, -1],
    [0, -1],
    [0,  0],
    [0,  0],
    [-1, -1],
    [0,  0],
    [0,  0],
    [0,  0],
    [0,  0]])

elements = np.array([
    [0, 1, 0, 0, 4, 8, 7],
    [1, 1, 0, 4, 1, 5, 8],
    [2, 1, 0, 7, 8, 6, 3],
    [3, 1, 0, 8, 5, 2, 6]])

mats = np.array([[1.0, 0.3]])

loads = np.array([
    [2, 0.0, 1.0],
    [3, 0.0, 1.0],
    [6, 0.0, 2.0]])

data = {"nodes": nodes,
        "cons": cons,
        "elements": elements,
        "mats": mats,
        "loads": loads}

### Run the simulation
disp = solids_auto(data)
plt.show()

solidspy_example2.py

import numpy as np
from solidspy.solids_GUI import solids_auto
from solidspy.postprocesor import plot_node_field
import matplotlib.pyplot as plt

### Define the data
nodes = np.array([
    [0, 0.00, 0.00],
    [1, 5.00, 0.00],
    [2, 10.00, 0.00],
    [3, 15.00, 0.00],
    [4, 0.00, 0.50],
    [5, 5.00, 0.50],
    [6, 10.00, 0.50],
    [7, 15.00, 0.50],
    [8, 0.00, 1.00],
    [9, 5.00, 1.00],
    [10, 10.00, 1.00],
    [11, 15.00, 1.00]])
nodes[:, 1:] *= 50

cons = np.array([
    [-1, -1],
    [0, 0],
    [0,  0],
    [0,  0],
    [-1, -1],
    [0,  0],
    [0,  0],
    [0,  0],
    [-1,  -1],
    [0,  0],
    [0,  0],
    [0,  0],])

# elements = np.array([
#     [0, 3, 0, 0, 5, 4],
#     [1, 3, 0, 0, 1, 5],
#     [2, 3, 0, 1, 6, 5],
#     [3, 3, 0, 1, 2, 6],
#     [4, 3, 0, 2, 7, 6],
#     [5, 3, 0, 2, 3, 7],
#     [6, 3, 0, 4, 9, 8],
#     [7, 3, 0, 4, 5, 9],
#     [8, 3, 0, 5, 10, 9],
#     [9, 3, 0, 5, 6, 10],
#     [10, 3, 0, 6, 11, 10],
#     [11, 3, 0, 6, 7, 11]])

elements = np.array([
    [0, 3, 0, 0, 5, 4],
    [1, 3, 0, 0, 1, 5],
    [2, 3, 0, 1, 6, 5],
    [3, 3, 0, 1, 2, 6],
    [4, 3, 0, 2, 7, 6],
    [5, 3, 0, 2, 3, 7],
    [6, 3, 0, 4, 5, 8],
    [7, 3, 0, 8, 5, 9],
    [8, 3, 0, 5, 6, 9],
    [9, 3, 0, 9, 6, 10],
    [10, 3, 0, 6, 7, 10],
    [11, 3, 0, 10, 7, 11]])

mats = np.array([[1e5, 0.45]])

loads = np.array([
    [3, 2500.0, 0.0],
    [7, 5000.0, 0.0],
    [11, 2500.0, 0.0]])

data = {"nodes": nodes,
        "cons": cons,
        "elements": elements,
        "mats": mats,
        "loads": loads}

#%% Run the simulation
disp = solids_auto(data, plot_contours=False)

#%% Visualization
# nodes_deform = nodes.copy()
# nodes_deform[:, 1:] += 10*disp
# # plot_node_field(disp[:, 1], nodes_deform, elements)
# plt.tricontourf(nodes_deform[:, 1], nodes_deform[:, 2],
#                 elements[:, 3:], disp[:, 1])
# plt.colorbar(shrink=0.8, ticks=[-0.02, 0, 0.02],
#              orientation="horizontal",
#              label="Vertical displacement")
# plt.triplot(nodes[:, 1], nodes[:, 2], elements[:, 3:],
#             color="black", alpha=0.4)
# plt.axis("image")
# plt.savefig("uniaxial_vert_disp.png", bbox_inches="tight",
#             dpi=100)
# plt.show()

solidspy_random_props.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Homogenization of a composite with random properties.

This examples runs with the developing version of SolidsPy (2.0).

@author: Nicolás Guarín-Zapata
@date: July 2020
"""
from time import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import spsolve
from solidspy.solids_GUI import solids_auto
from solidspy.preprocesor import rect_grid
from solidspy.assemutil import assembler, DME, loadasem
from solidspy.postprocesor import complete_disp


np.random.seed(3)
G1_homog = []
G2_homog = []
nels_x = list(range(10, 110, 10))
for nx in nels_x:
    # Setup
    x, y, els = rect_grid(1, 1, nx, nx)
    nels = els.shape[0]
    npts = x.shape[0]
    nodes = np.zeros((npts, 3))
    nodes[:, 1] = x
    nodes[:, 2] = y

    # Constraints 1
    cons1 = np.zeros((npts, 2), dtype=int)
    cons1[np.abs(x + 0.5) < 1e-6, 0] = -1
    cons1[np.abs(y + 0.5) < 1e-6, 1] = -1

    # Loads 1
    loads1 = np.zeros((npts, 3))
    loads1[:, 0] = range(npts)
    loads1[x == 0.5, 1] = 1.0
    loads1[(x == -0.5) * (y == 0.5), 1] = 0.5
    loads1[(x == 0.5) * (y == 0.5), 1] = 0.5
    loads1[y == 0.5, 2] = -1.0
    loads1[(x == 0.5) * (y == 0.5), 2] = -0.5
    loads1[(x == 0.5) * (y == -0.5), 2] = -0.5

    # Constraints 2
    cons2 = np.zeros((npts, 2), dtype=int)
    cons2[(y == -0.5), 0] = -1
    cons2[np.abs(y + 0.5) < 1e-6, 1] = -1

    # Loads 2
    loads2 = np.zeros((npts, 3))
    loads2[:, 0] = range(npts)
    loads2[y == 0.5, 1] = -1.0
    loads2[(x == 0.5) * (y == 0.5), 1] = -0.5
    loads2[(x == 0.5) * (y == -0.5), 1] = -0.5

    # Materials
    G1_acu = 0
    G2_acu = 0
    nsamples = 5
    for cont in range(nsamples):
        mats = np.ones((nels, 2))
        mats[:, 0] = np.random.normal(loc=8/3, scale=0.5, size=nels)
        mats[:, 1] = 1/3
        els[:, 2] = range(nels)

        # Data
        data1 = {"nodes": nodes,
                "cons": cons1,
                "elements": els,
                "mats": mats,
                "loads": loads1}
        data2 = {"nodes": nodes,
                "cons": cons2,
                "elements": els,
                "mats": mats,
                "loads": loads2}

        # Run the simulation 1
        disp, strain, stress = solids_auto(data1, compute_strains=True,
                                                plot_contours=False)
        Sxy = 0.5*np.abs(stress[:, 0] - stress[:, 1])
        Exy = np.abs(strain[:, 0] - strain[:, 1])
        G = np.mean(Sxy)/np.mean(Exy)
        G1_acu += G
        # Run the simulation 2
        disp, strain, stress = solids_auto(data2, compute_strains=True,
                                                plot_contours=False)
        Sxy = 0.5*np.abs(stress[:, 0] - stress[:, 1])
        Exy = np.abs(strain[:, 0] - strain[:, 1])
        G = np.mean(Sxy)/np.mean(Exy)
        G2_acu += G
    G1_homog.append(G1_acu/nsamples)
    G2_homog.append(G2_acu/nsamples)


#%% Plotting
plt.plot(nels_x, G1_homog)
plt.plot(nels_x, G2_homog)
plt.xlabel("Number of elements per side")
plt.ylabel(r"Shear modulus")
plt.legend(["$G_1$", "$G_2$"])