<jupyter_start><jupyter_text>Malla básicaPara crear una malla, lo primero es crear todos sus nodossupongamos que queremos una malla que va desde 0 a 10 en x con 11 nodos y lo mismo en yEl total de nodos para una malla con $nx \times ny $ nodos es:$$ N_{nodos}=nx\times ny$$en este caso $ 11\times 11 =121$<jupyter_code>from numpy import *
import matplotlib.pyplot as plt
%matplotlib inline

xmin, xmax,ymin, ymax, nx, ny = 0, 10, 0, 10, 5, 5
x = linspace(xmin, xmax, nx)
y = linspace(ymin, ymax, ny)
xv, yv = meshgrid(x, y)
points = array([xv.flatten(),yv.flatten()]).T

plt.plot(points[:,0], points[:,1],'.' )
plt.show()
points[2+1*nx]<jupyter_output><empty_output><jupyter_text>El total de elementos está dado por:$$ n_{elementos}=(nx-1)\times(ny-1)\times 2 $$En este caso, serán $(11-1)\times(11-1)\times 2=200$<jupyter_code>nElem=(nx-1)*(ny-1)*2
T= zeros([nElem,3])                          # matriz para de conectividad cada fila tiene los tres nodos de un elemento
nElem<jupyter_output><empty_output><jupyter_text>Para crear la conectividad que define cada elemento![nodos.png](attachment:nodos.png)<jupyter_code>e=0;                                        # índice del elemento
for f in range(0, ny-1):                  # para cada fila excepto la última
    for nodo in range(0, nx-1):              # para cada nodo de la fila excepto el último
        i= nodo+nx*f                       # índice del nodo
        T[e,:]=[i, i+1 , i+nx+1]
        e+=1
        T[e,:]=[i, i+nx+1, i+nx]
        e+=1            
e                                             # verificamos el número total de elementos completados
plt.triplot(points[:,0], points[:,1], T.copy())<jupyter_output><empty_output><jupyter_text>Sin embargo en esa malla dos elementos (en verde) cuyos nodos están todos en el borde. Esto en muchos casos es indeseable, por ejemplo si se fijan condiciones de contorno en todos los nodos exteriores, esos nodos no n hay conexión con el resto de la malla. Por lo que hay 4 elementos que deben ser modificados, ![Captura%20de%20pantalla%20de%202020-05-19%2023-58-50.png](attachment:Captura%20de%20pantalla%20de%202020-05-19%2023-58-50.png)Eestos son los dos últimos creados en la primer fila de nodos$$ (nx-2)\times 2$$$$ (nx-2)\times 2 +1 $$y los dos primeros creados en la anteúltima fila de nodos$$ (nx-1)\times(ny-2)\times 2 $$$$ (nx-1)\times(ny-2)\times 2+1 $$<jupyter_code>(nx-2)*2, (nx-2)*2+1 ,(nx-1)*(ny-2)*2, (nx-1)*(ny-2)*2+1<jupyter_output><empty_output><jupyter_text>Para los dos dos útlimos nodos de la primer fila empezamos con el nodo  $ i=nx-2 $ y definimos la triangulación![Captura%20de%20pantalla%20de%202020-05-19%2023-53-59.png](attachment:Captura%20de%20pantalla%20de%202020-05-19%2023-53-59.png)Para los otros dos elementos empezamos con el nodo $ i=nx \times (ny-2) $ y definimos la conectividad de los elementos:<jupyter_code>i=nx-2       # verifico los nodos que busco
T[(nx-2)*2+1,:],T[(nx-2)*2+2,:], [i ,  i+1 , i+nx],[i+1, i+1+nx , i+nx]
i=nx-2
T[(nx-2)*2+1,:]=[i ,  i+1 , i+nx]
T[(nx-2)*2,:]=[i+1, i+1+nx , i+nx]

i=nx*(ny-2)
T[(nx-1)*(ny-2)*2+1,:] = [i ,  i+1 , i+nx]
T[(nx-1)*(ny-2)*2,:] = [i+1, i+1+nx , i+nx]

plt.triplot(points[:,0], points[:,1], T.copy())<jupyter_output><empty_output><jupyter_text>Encapsulamiento en una función<jupyter_code>from numpy import *
from matplotlib import pyplot as plt
#%matplotlib inline
%matplotlib notebook
%matplotlib notebook


def malla(xmin, xmax, nx, ymin, ymax, ny, correct_corner=True, plot=False, _3D=True):
    x = linspace(xmin, xmax, nx)
    y = linspace(ymin, ymax, ny)
    xv, yv = meshgrid(x, y)
    points = array([xv.flatten(),yv.flatten()]).T
    
    nElem=(nx-1)*(ny-1)*2
    T= zeros([nElem,3])                          # matriz para de conectividad cada fila tiene los tres nodos de un elemento

    e=0;                                        # índice del elemento
    for f in range(0, ny-1):                  # para cada fila excepto la última
        for nodo in range(0, nx-1):              # para cada nodo de la fila excepto el último
            i= nodo+nx*f                       # índice del nodo
            T[e,:]=[i, i+1 , i+nx+1]
            e+=1
            T[e,:]=[i, i+nx+1, i+nx]
            e+=1   
    
    if(correct_corner):                               # corrige las esquinas 
        i=nx-2                                        # últimos dos elementos de la primer fila
        T[(nx-2)*2+1,:]=[i ,  i+1 , i+nx]
        T[(nx-2)*2,:]=[i+1, i+1+nx , i+nx]

        i=nx*(ny-2)                                   # primeros dos elementos de la última fila
        T[(nx-1)*(ny-2)*2+1,:] = [i ,  i+1 , i+nx]
        T[(nx-1)*(ny-2)*2,:] = [i+1, i+1+nx , i+nx]
    
    if(plot):
        plt.triplot(points[:,0], points[:,1], T.copy())
    
    if(_3D):                                          # agrego una columna de 0 para z
        P3=zeros([len(points),3])
        P3[:,0:2]=points
        return P3, T
    
    return points, T                                  # devuelve los nodos y su conectividad
    
X,T= malla(0,10,11,0,10,11, plot=True)
import matplotlib.tri as mtri
from mpl_toolkits.mplot3d import Axes3D 
points,T= malla(0,10,11,0,10,11, _3D=True)
fig = plt.figure(figsize=plt.figaspect(0.5))

ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_trisurf(points[:,0], points[:,1], points[:,2], triangles=T)<jupyter_output><empty_output><jupyter_text>Identificando nodos del borde y esquinasEn muchos casos puede ser útil identificar aquellos nodos de los bordes y de las esquinas... En este caso los identificaremos con una lista de números llamada Creases:* nodos interiores =0* nodos de los bordes =1* nodos de las esquinas  =2<jupyter_code>creases=(X[:,0]==0) *1 + (X[:,0]==10 )*1+  (X[:,1]==0)*1+ (X[:,1]==10)*1 
# *1 convierte de boolan a int para sumar sin comparar any o all

creases
#cantidad de elementos internos, bordes, esquinas
count_nonzero(creases == 0),count_nonzero(creases == 1),count_nonzero(creases == 2)<jupyter_output><empty_output><jupyter_text>Convirtiendo a VTKUsando la librería PyVista<jupyter_code>def toPyVista(X,T):    # devuelve un objeto compatible de pyvista para facilitar las operaciones
    b = zeros((len(T),4))
    b[:,1:4]=T.copy()
    b[:,0]=3*ones([len(T)], dtype=int)
    faces=hstack(b)
    
    import pyvista as pv
    surf = pv.PolyData(X,faces)
    return surf

surf1=toPyVista(X,T)
surf1.save('ej.vtk')    # desde PyVista guarda en VTK<jupyter_output><empty_output><jupyter_text>DesplazamientoEn este caso se desplaza 10 unidades en x<jupyter_code>X_=X.copy()
X_[:,0]+=10
surf2=toPyVista(X_,T)
surf2.save('ej2.vtk')<jupyter_output><empty_output><jupyter_text>y se junta con la malla anterior<jupyter_code>surf=surf1+surf2
surf.plot(show_edges=True)
surf.save('doble.vtk')
P=surf.points
creases=(P[:,0]==0)*1+(P[:,0]==20)*1+(P[:,1]==0)*1+(P[:,1]==10)*1+ (P[:,0]==10)*1
#cantidad de elementos internos, bordes, esquinas
count_nonzero(creases == 0),count_nonzero(creases == 1),count_nonzero(creases == 2)

import numpy as np
import vtk
from vtk.numpy_interface import dataset_adapter as dsa   # Para convertir de vtk a numpy

reader = vtk.vtkPolyDataReader()                         # SOURCE/READER 
reader.SetFileName('doble.vtk')                   
reader.Update()                                          
data_vtk=reader.GetOutput()
data = dsa.WrapDataObject(data_vtk)                      # convierte la salida a array compatible con numpy

N_points=data.GetNumberOfPoints()                        # número de nodos
Points=np.array(data.Points)                             # nodos

unq, index = np.unique(data.Points,                      # filtro valores uniq del array de puntos
                       axis=0,                           # compara filas enteras
                       return_index=True)                # índice de los no repetidos

master=np.zeros(N_points,dtype=int)
for i in range(0, N_points):                             # para cada nodo i
    if i in index:                                       # si es único
        master[i]=i                                      # su master es i
    else:                                                # si no es único
        for j in range(0,i):                             # comparo con cada nodo anterior  j
            if np.all(Points[i]==Points[j]):             # si es igual  al nodo j (el primero que aparece)
                master[i]=j                              # su master es j
                break                                    # detengo el bucle for
master

np.all(Points==Points[master])                            # si da True el algoritmo para los masters funciona bien
Master=dsa.numpyTovtkDataArray(master, name = "master" )   # convierto la lista master a un objeto vtkDataArray 
data_vtk.GetPointData().SetScalars(Master)                 # agrego en PointData al vtkArray llamado Master (NO ES EL data CONVERTIDO A NUMPY)

creases=dsa.numpyTovtkDataArray(creases,name = "CreasesId") # convierto la lista master a un objeto vtkDataArray 
data_vtk.GetPointData().AddArray(creases)                  # para agregar un segundo atributo uso AddArray()      ver http://vtk.1045678.n5.nabble.com/Multiple-fields-scalars-vectors-td1224224.html

#constraint=dsa.numpyTovtkDataArray(A, name= "ConstraintTags")
#data_vtk.GetPointData().AddArray(constraint)

writer = vtk.vtkPolyDataWriter()                           # creo el objeto PolyDataWriter
writer.SetFileName('doble.vtk')             # le pongo nombre al archivo de salida
writer.SetInputData(data_vtk)                              # le conecto los datos data_vtk
writer.SetFileType(1)                                      # opcional: set file type to ascii
print(writer.Write())                                             # escribo el archivo .vtk  Returns 1 on success and 0 on failure.
del(writer)<jupyter_output>1<jupyter_text>Rotación<jupyter_code>def rotZ(th):                # matriz de rotación Z con th radianes
    rot=array([                         
        [1,0,0],
        [0, cos(th), -sin(th)],
        [0, sin(th), cos(th)]
    ])
    return rot
surf=toPyVista(X,T)
for i in range (1,8):
    surf+=toPyVista(X@rotZ(i*3.0/4.0*pi),T)

surf.plot(show_edges=True)
surf.save("multi.vtk")<jupyter_output><empty_output><jupyter_text>Cilindro conforme Hollow Fiber<jupyter_code>#este es el ángulo para hace un lado de malla de cilindro ~=0.1
import numpy as np
def F(x):
    return (100*np.sin(x)*np.sin(x)+100-200*np.cos(x)+100*np.cos(x)*np.cos(x)-0.1)
import scipy.optimize
x = scipy.optimize.broyden1(F, [0.5], f_tol=1e-11)

(2*np.pi/x)%8   # entonces 6 * 8 =48 triángulos para dar la vuelta al cilindro => 2*pi/48 es el ángulo
#entonces 48*11=528 nodos
from numpy import *
X=zeros([11*48,3])

num=0
for j in range(0,11):
    for i in range(0,48):
        th = i*2*pi/48
        x, y, z = j, 10*cos(th), 10*sin(th)
        X[num]=[x,y,z]
        num+=1
X
from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt
%matplotlib notebook
%matplotlib notebook

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter3D(X[:,0],X[:,1],X[:,2], 'gray')
Nz=11
Nt=48
nElem = 2*(Nz-1)*Nt;
T= zeros([nElem,3])                          # matriz para de conectividad cada fila tiene los tres nodos de un elemento

e=0;                                         # índice del elemento
for k in range(0, (Nz-1)):                   # para cada anillo excepto el último
    for i in range(0,Nt-1):                # para cada nodo del anillo
        ii= (k)*Nt+i                     # el índice del nodo en la lista global de nodos
        T[e,:]=[ii, ii+1, Nt+ii+1]       # construye el elemento
        e+=1                             # paso al siguiente elemento
        
        T[e,:] = [ii, Nt+ii+1, Nt+ii]
        e = e+1
       
    T[e,:]= [(k+1)*Nt-1, k*Nt, (k+1)*Nt]
    e+=1                                 # paso al siguinte elemento
    
    T[e,:]=[(k+1)*Nt-1,  (k+1)*Nt,  (k+2)*Nt-1]
    e+=1
b = zeros((960,4))
b[:,1:4]=T
b[:,0]=3*ones([960], dtype=int)
faces=hstack(b)

import pyvista as pv
from itkwidgets import view



surf2 = pv.PolyData(X,faces)+surf
surf2.plot(show_edges=True)
surf2.save("hollow_fiver.vtk")<jupyter_output><empty_output><jupyter_text>Más transformaciones<jupyter_code>L0 = -0.5 * h                                  # para que el cilindro quede centrado

X = zeros([Nt*Nz,3])

ptId = 0
for i in range(0, Nz):                         # para cada anillo
    for j in range(0, Nt):                     # para cada punto del anillo
        theta = dt * (j + 0.5*mod(i,2))
        X[ptId,0] = R * cos(theta)
        X[ptId,1] = R * sin(theta)
        X[ptId,2] = (i*dz) + L0 
        ptId = ptId+1
th=-3
rot=array([                         # matriz de rotación
    [1,0,0],
    [0, cos(th), -sin(th)],
    [0, sin(th), cos(th)]
])
points=X@rot                        # matriz de puntos rotada

writeVtk(points, T, 'basico2')
import matplotlib.tri as mtri
from mpl_toolkits.mplot3d import Axes3D  

fig = plt.figure(figsize=plt.figaspect(0.5))

ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_trisurf(points[:,0], points[:,1], points[:,2], triangles=T)<jupyter_output><empty_output><jupyter_text>Transformación Bending<jupyter_code>import numpy as np

def bending(X ,A=1,B=1,C=0.5):
    X1, X2, X3=X[:,0], X[:,1], X[:,2]
    r=np.sqrt(2*A*X1)
    theta=B*X2
    #z=X3/(A*B)-B*C*X2
    
    X_=X.copy()              # si uso X_=X crea una referencia!!!
    X_[:,0]=r*np.cos(theta)
    X_[:,1]=r*np.sin(theta)
    
    return X_
X,T= malla(1,10,11,1,10,11, _3D=True)
points= bending(X, 1, 0.1, 1)
import matplotlib.tri as mtri
from mpl_toolkits.mplot3d import Axes3D  

fig = plt.figure(figsize=plt.figaspect(0.5))

ax = fig.add_subplot(1,1, 1, projection='3d')
ax.plot_trisurf(points[:,0], points[:,1], points[:,2], triangles=T)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13974998193511093, 'fraction_numerical': 0.040393091986415205, 'mean_word_length': 2.550538737814264, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '165016', 'license': 'MIT', 'n_tokens_mistral': 5466, 'n_tokens_neox': 5273, 'n_words': 1241, 'path': 'simple.ipynb', 'repo_name': 'ieferrari/python-3D-meshes', 'size': '13839'}