<jupyter_start><jupyter_text>**Q4 - MEMBRANE** ***<jupyter_code>%matplotlib inline<jupyter_output><empty_output><jupyter_text>**Step 1 : Create a Mesh**<jupyter_code>from pyoneer.mesh.rgrid import rgrid_Q4
import numpy as np

size = Lx, Ly = (10, 10)
shape = nx, ny = (100, 100)

gridparams = {
    'size' : size,
    'shape' : shape,
    'origo' : (0, 0),
    'start' : 0
    }

coordsQ4_, topoQ4 = rgrid_Q4(**gridparams)
coordsQ4 = np.zeros((coordsQ4_.shape[0], 3))
coordsQ4[:, :2] = coordsQ4_[:, :]<jupyter_output><empty_output><jupyter_text>**Step 2 : Define a Material**<jupyter_code>E = 12000
nu = 0.2
t = 0.25
A = np.array([[1, nu, 0], [nu, 1, 0], [0., 0, (1-nu)/2]]) * (t * E /(1-nu**2))<jupyter_output><empty_output><jupyter_text>**Step 3 : Set Boundary Conditions**<jupyter_code># fix points at x==0
cond = (coordsQ4[:, 0] <= 0.001) & (coordsQ4[:, 1] >= (Ly/2))
ebcinds = np.where(cond)[0]
fixity = np.zeros((coordsQ4.shape[0], 2), dtype=bool)
fixity[ebcinds, :] = True

# unit vertical load on points at (Lx, Ly)
Fz = -10
cond = (coordsQ4[:, 0] > (Lx-(1e-12))) & (coordsQ4[:, 1] > (Ly-(1e-12)))
nbcinds = np.where(cond)[0]
loads = np.zeros((coordsQ4.shape[0], 2))
loads[nbcinds, 1] = Fz<jupyter_output><empty_output><jupyter_text>**Step 4 : Assembly and Solution**<jupyter_code>from pyoneer.mech.fem.mesh import FemMesh
from pyoneer.mech.fem.structure import Structure
from pyoneer.mech.fem.elements.Q4 import Q4_M as Q4

Q4._selective_= False
meshQ4 = FemMesh(coords=coordsQ4, topo=topoQ4, model=A, fixity=fixity, \
    loads=loads, celltype=Q4)
structure = Structure(mesh=meshQ4)
structure.linsolve(permc_spec='NATURAL', sparsify=True)
structure.summary
dofsol = np.zeros_like(coordsQ4)
dofsol[:, :2] = structure.mesh.pointdata.dofsol.to_numpy()
structure.mesh.pointdata['x'] = coordsQ4 + dofsol
mb = structure.mesh.to_pv()
mb.GetBlock(0).point_data['dofsol'] = dofsol[:, 1]
#mb.plot(show_edges=False, notebook=False, cpos="xy")<jupyter_output><empty_output><jupyter_text>**Step 5 : Triangulate and Plot**<jupyter_code>from pyoneer.mesh.tr import Q4_to_T3
from pyoneer.mesh.trimesh import triangulation
from pyoneer.mesh.triplot import triplot
import matplotlib.pyplot as plt
from matplotlib import gridspec
import matplotlib as mpl
plt.style.use('dark_background')

points, triangles = Q4_to_T3(structure.mesh.coords(), topoQ4)
triobj = triangulation(points=points[:, :2], triangles=triangles)[-1]
fig = plt.figure(figsize=(10, 4))  # in inches
fig.patch.set_facecolor('black')
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1])
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1])
triplot(triobj, ax=ax1, fig=fig)
triplot(triobj, ax=ax2, fig=fig, 
        data=dofsol[:, 1], cmap='turbo', axis='off')<jupyter_output><empty_output><jupyter_text>**Step 6 : Optimize**<jupyter_code>%matplotlib qt
plt.style.use('classic')
from pyoneer.mech.topopt.SIMP_OC_FEM import OC_SIMP_COMP as OC
# iteration parameters
OC_params = {
    'p_start' : 1.0,  # SIMP penalty factor
    'p_stop' : 3.0,
    'p_inc' : 0.1,
    'p_step' : 5,
    'q' : 0.5,  # smoothing factor
    'vfrac' : 0.4,  # fraction of target volume over initial volume
    'dtol' : 0.1,  #  to control maximum change in the variables
    'r_min' : 1.1 * max(Lx/nx, Ly/ny),  # for the density filter
    'miniter' : 30,
    'maxiter' : 100
    }
from pyoneer.math import histogram
from matplotlib import gridspec
import matplotlib as mpl
from IPython.display import display, clear_output

fig = plt.figure(figsize=(10, 4))  # in inches
#fig.patch.set_facecolor('white')
gs = gridspec.GridSpec(2, 3, width_ratios=[6, 2, 4])

uz = structure.mesh.pointdata.dofsol.to_numpy()[nbcinds, 1].min()
nCell = structure.mesh.number_of_cells()
history = {'comp' : [], 'vol' : [], 'x' : []}
history['comp'].append(2*uz*Fz) 
history['x'].append(np.random.rand(nCell))
history['vol'].append(Lx*Ly)

nbins = 10
hist, bin_centers = histogram(history['x'][-1], nbins)
hist = hist.astype(float)
hist /= hist.max()
ax1 = fig.add_subplot(gs[:, 1])
ax1.set_facecolor('white')
bars = ax1.barh(bin_centers, hist, 0.5/nbins)
ax1.set_xlim(0, 1)
ax1.axis('off')

ax2 = fig.add_subplot(gs[0, 2])
ax2.axhline(y=history['comp'][-1], color="r", linestyle="--", lw=0.5)
#ax2.set_ylim(0, history['comp'][-1]*1.2)

ax3 = fig.add_subplot(gs[1, 2])
ax3.set_ylim(0, history['vol'][-1])
ax3.axhline(y=history['vol'][-1]*OC_params['vfrac'], 
            color="r", linestyle="--", lw=0.5)

cmap = plt.cm.winter  # define the colormap
# extract all colors from the .jet map
cmaplist = [cmap(i) for i in range(cmap.N)]
# force the first color entry to be grey
#cmaplist[0] = (.5, .5, .5, 1.0)
#cmaplist[-1] = (1., 0., 0., 1.0)
#cmaplist.reverse()
cmaplist[0] = (1., 1., 1., 1.0)
cmaplist[-1] = (1., 0., 0., 1.0)

# create the new map
cmap = mpl.colors.LinearSegmentedColormap.from_list(
    'Custom cmap', cmaplist, cmap.N)

# define the bins and normalize
bounds = np.linspace(0, 10, 11)
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

ax4 = fig.add_subplot(gs[:, 0])
points, triangles, edata = \
    Q4_to_T3(coordsQ4, topoQ4, edata=history['x'][-1])
triobj = triangulation(points=points[:, :2], triangles=triangles)[-1]
trifield = triplot(triobj, ax=ax4, data=edata, axis='off', 
                   lw=0.0, fig=fig, cmap=cmap)[0]


def callback_qt(i, comp, vol, dens):
    history['comp'].append(comp)
    history['vol'].append(vol)
    history['x'].append(dens)
    fig.canvas.manager.window.raise_()
    hist, _ = histogram(dens, nbins)
    hist = hist.astype(float)
    hist /= hist.max()
    for bar, h in zip(bars, hist):
        bar.set_width(h)
    ax2.plot(i, comp, marker='o', c='b', markersize='1')
    ax3.plot(i, vol, marker='*', c='g', markersize='1')
    ax2.set_xlim(0, i+1)
    ax3.set_xlim(0, i+1) 
    *_, edata = Q4_to_T3(coordsQ4, topoQ4, edata=dens)
    trifield.set_array(edata)
    fig.canvas.draw()
    fig.canvas.flush_events()
    
def callback_inline(i, comp, vol, dens):
    history['comp'].append(comp)
    history['vol'].append(vol)
    history['x'].append(dens)
    hist, _ = histogram(dens, nbins)
    hist = hist.astype(float)
    hist /= hist.max()
    
    ax1.cla()
    ax1.barh(bin_centers, hist, 0.5/nbins)
    ax2.plot(i, comp, marker='o', c='b', markersize='1')
    ax3.plot(i, vol, marker='*', c='g', markersize='1')
    ax4.cla()
    *_, edata = Q4_to_T3(coordsQ4, topoQ4, edata=dens)
    triplot(triobj, ax=ax4, data=edata, axis='off', lw=0.0, 
            fig=fig, cmap=cmap)[0]
    
    display(fig)
    clear_output(wait=True)
    plt.pause(0.1)
    
    
res, sum = OC(structure.mesh,  summary=True, callback=callback_qt, **OC_params)
ax2.set_xlim(0, len(history['comp']))
ax3.set_xlim(0, len(history['comp']))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1831845238095238, 'fraction_numerical': 0.038839285714285715, 'mean_word_length': 4.568351284175642, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '246934', 'license': 'MIT', 'n_tokens_mistral': 2953, 'n_tokens_neox': 2832, 'n_words': 549, 'path': 'tests/fem/test_Q4M.ipynb', 'repo_name': 'dewloosh/dewloosh-solid', 'size': '6720'}