def createQuadPart(res=None):
    """Create a quadrilateral domain from user input"""
    global x0,y0,x1,y1,x2,y2,x3,y3,nx,ny,eltype
    if model is not None:
        if ask('You have already merged the parts! I can not add new parts anymore.\nYou should first delete everything and recreate the parts.',['Delete','Cancel']) == 'Delete':
            deleteAll()
        else:
            return
    if res is None:
        res = askItems([
            I('Vertex 0',(x0,y0)),
            I('Vertex 1',(x1,y1)),
            I('Vertex 2',(x2,y2)),
            I('Vertex 3',(x3,y3)),
            I('nx',nx),
            I('ny',ny),
            I('eltype',eltype,itemtype='radio',choices=['quad','tri-u','tri-d']),
            ])
    if res:
        x0,y0 = res['Vertex 0']
        x1,y1 = res['Vertex 1']
        x2,y2 = res['Vertex 2']
        x3,y3 = res['Vertex 3']
        nx = res['nx']
        ny = res['ny']
        eltype = res['eltype']
        diag = {'quad':'', 'tri-u':'u', 'tri-d':'d'}[eltype]
        xold = rectangle(1,1).coords
        xnew = Coords([[x0,y0],[x1,y1],[x2,y2],[x3,y3]])
        M = rectangle(nx,ny,1.,1.,diag=diag).toMesh().isopar('quad4',xnew,xold)
        addPart(M)

def cone(r0,r1,h,t=360.,nr=1,nt=24,diag=None):
    """Constructs a Formex which is (a sector of) a
    circle / (truncated) cone / cylinder.

    r0,r1,h are the lower and upper radius and the height of the truncated
    cone. All can be positive, negative or zero.
    Special cases:
    r0 = r1 : cylinder
    h = 0 : (flat) circle
    r0 = 0 or r1 = 0 : untruncated cone

    Only a sector of the structure, with opening angle t, is modeled.
    The default results in a full circumference.

    The cone is modeled by nr elements in height direction and nt elements in
    circumferential direction.
    
    By default, the result is a 4-plex Formex whose elements are quadrilaterals
    (some of which may collapse into triangles).
    If diag='up' or diag = 'down', all quads are divided by an up directed
    diagonal and a plex-3 Formex results.
    """
    B = simple.rectangle(nt,nr,1.,1.,diag=diag) # grid with size 1x1
    B = B.map(lambda x,y,z:[x,y,r0-y*(r0-r1)]) # translate and tilt it
    B = B.scale([t,h,1.])             # scale to fit parameters
    return B.cylindrical(dir=[2,0,1]) # roll it into a cone

def run():
    clear()
    smoothwire()

    nx = 4
    ny = 3
    nz = 7

    delay(2)

    # A rectangular mesh
    M1 = simple.rectangle(nx,ny).toMesh().setProp(1)
    # Same mesh, rotated and translated
    M2 = M1.rotate(45,0).translate([1.,-1.,nz]).setProp(3)
    draw([M1,M2])

    # Leave out the first and the last two elements
    sel = arange(M1.nelems())[1:-2]
    m1 = M1.select(sel)
    m2 = M2.select(sel)
    clear()
    draw([m1,m2],view=None)

    # Connect both meshes to a hexaeder mesh
    m = m1.connect(m2,nz)
    clear()
    draw(m,color=red,view=None)

def createRectPart(res=None):
    """Create a rectangular domain from user input"""
    global x0,y0,x2,y2,nx,ny,eltype
    if model is not None:
        if ask('You have already merged the parts! I can not add new parts anymore.\nYou should first delete everything and recreate the parts.',['Delete','Cancel']) == 'Delete':
            deleteAll()
        else:
            return
    if res is None:
        res = askItems([
            I('x0',x0,tooltip='The x-value of one of the corners'),
            I('y0',y0),
            I('x2',x2),I('y2',y2),
            I('nx',nx),I('ny',ny),
            I('eltype',eltype,itemtype='radio',choices=['quad','tri-u','tri-d']),
            ])
    if res:
        globals().update(res)
        if x0 > x2:
            x0,x2 = x2,x0
        if y0 > y2:
            y0,y2 = y2,y0
        diag = {'quad':'', 'tri-u':'u', 'tri-d':'d'}[eltype]
        M = rectangle(nx,ny,x2-x0,y2-y0,diag=diag).toMesh().trl([x0,y0,0])
        addPart(M)

def createGrid():
    res = askItems([_I("name", "__auto__"), _I("nx", 3), _I("ny", 3), _I("b", 1), _I("h", 1)])
    if res:
        globals().update(res)
        # name = res['name']
        # nx = res['nx']
        # ny = res['ny']
        # b =
        S = TriSurface(simple.rectangle(nx, ny, b, h, diag="d"))
        export({name: S})
        selection.set([name])
        selection.draw()

def createGrid():
    res = askItems([('name','__auto__'),('nx',3),('ny',3),('b',1),('h',1)])
    if res:
        globals().update(res)
        #name = res['name']
        #nx = res['nx']
        #ny = res['ny']
        #b = 
        S = TriSurface(simple.rectangle(nx,ny,b,h,diag='d'))
        export({name:S})
        selection.set([name])
        selection.draw()

def showSuperEgg():
    """Draw a super Egg from set global parameters"""
    nx,ny = grid
    b,h = long_range[1]-long_range[0], lat_range[1]-lat_range[0]
    if grid_base.startswith('tri'):
        diag = grid_base[-1]
    else:
        diag = ''
    B = rectangle(nx,ny,b,h,diag=diag,bias=grid_bias).translate([long_range[0],lat_range[0],1.])
    if grid_skewness != 0.0:
        B = B.shear(0,1,grid_skewness)
    F = B.superSpherical(n=north_south,e=east_west,k=eggness).scale(scale)
    clear()
    draw(F,color=color)
    export({name:F})

def run():
    reset()
    smoothwire()
    transparent()
    lights(True)

    nx,ny = 20,10

    F = simple.rectangle(nx,ny)
    F = F.trl(-F.center()+[0.,0.,nx/2])
    draw(F)

    G = F.projectOnSphere(ny)
    draw(G,color=red)

    H = F.rotate(30).projectOnCylinder(ny)
    draw(H,color=blue)

def createGrid():
    """Create the grid from global parameters"""
    global B
    nx,ny = grid_size
    b,h = x_range[1]-x_range[0], y_range[1]-y_range[0]
    if grid_base.startswith('tri'):
        diag = grid_base[-1]
    else:
        diag = ''
    B = rectangle(nx,ny,b,h,diag=diag,bias=grid_bias).translate([x_range[0],y_range[0],1.])
    if grid_skewness != 0.0:
        B = B.shear(0,1,grid_skewness*b*ny/(h*nx))
    if x_clip:
        B = B.clip(B.test('any',dir=0,min=x_clip[0]+tol*b,max=x_clip[1]-tol*b))
    if y_clip:
        B = B.clip(B.test('any',dir=1,min=y_clip[0]+tol*h,max=y_clip[1]-tol*h))
    export({grid_name:B})

def createPart(res=None):
    """Create a rectangular domain from user input"""
    global x0,y0,x1,y1
    if model is not None:
        if ask('You have already merged the parts! I can not add new parts anymore.\nYou should first delete everything and recreate the parts.',['Delete','Cancel']) == 'Delete':
            deleteAll()
        else:
            return
    if res is None:
        res = askItems([('x0',x0),('y0',y0),
                        ('x1',x1),('y1',y1),
                        ('nx',nx),('ny',ny),
                        ('eltype',eltype,'select',['quad','tri-u','tri-d']),
                        ])
    if res:
        globals().update(res)
        if x0 > x1:
            x0,x1 = x1,x0
        if y0 > y1:
            y0,y1 = y1,y0
        diag = {'quad':'', 'tri-u':'u', 'tri-d':'d'}[eltype]
        F = rectangle(nx,ny,x1-x0,y1-y0,diag=diag).trl([x0,y0,0])
        addPart(F)
        drawParts()

def createPart(res=None):
    """Create a rectangular domain from user input"""
    global x0, y0, x1, y1
    if model is not None:
        if (
            ask(
                "You have already merged the parts! I can not add new parts anymore.\nYou should first delete everything and recreate the parts.",
                ["Delete", "Cancel"],
            )
            == "Delete"
        ):
            deleteAll()
        else:
            return
    if res is None:
        res = askItems(
            [
                ("x0", x0),
                ("y0", y0),
                ("x1", x1),
                ("y1", y1),
                ("nx", nx),
                ("ny", ny),
                ("eltype", eltype, "select", ["quad", "tri-u", "tri-d"]),
            ]
        )
    if res:
        globals().update(res)
        if x0 > x1:
            x0, x1 = x1, x0
        if y0 > y1:
            y0, y1 = y1, y0
        diag = {"quad": "", "tri-u": "u", "tri-d": "d"}[eltype]
        F = rectangle(nx, ny, x1 - x0, y1 - y0, diag=diag).trl([x0, y0, 0])
        addPart(F)
        drawParts()

##

"""Projection

level = 'normal'
topics = ['geometry','surface']
techniques = ['colors']

"""

import simple
from gui.canvas import *

nx,ny = 20,10

F = simple.rectangle(nx,ny)
F = F.trl(-F.center()+[0.,0.,nx/2])
draw(F)

x = F.f.projectOnSphere(ny)
G = Formex(x)
draw(G,color=red)

x = F.f.rotate(30).projectOnCylinder(ny)
H = Formex(x)
draw(H,color=blue)

smooth()
n=200
for i in range (n):
    v = float(i)/(2*n)

"""

import simple
from connectivity import reverseUniqueIndex
from plugins.mesh import *
        
def drawMesh(mesh,ncolor='blue',ecolor='red'):
    if ncolor:
        draw(mesh.coords,color=ncolor)
    if ecolor:
        draw(mesh,color=ecolor,bbox='last')

nx = 4
ny = 3
nz = 7
F = simple.rectangle(nx,ny).setProp(1)

c1,e1 = F.feModel()
c2 = c1.rotate(45,0).translate([1.,-1.,nz])

G = Formex(c2[e1]).setProp(3)
draw([F,G])


e1 = e1[1:-2]

m1 = Mesh(c1,e1)
m2 = Mesh(c2,e1)

clear()
drawMesh(m1)

def run():
    # GEOMETRICAL PARAMETERS FOR HE200B wide flange beam
    h = 200. #beam height
    b = 200. #flange width 
    tf = 15. #flange thickness
    tw = 9.  #body thickness
    l = 400. #beam length
    r = 18.  #filling radius

    # MESH PARAMETERS
    el = 20 #number of elements along the length
    etb = 2 #number of elements over half of the thickness of the body
    ehb = 5 #number of elements over half of the height of the body
    etf = 5 #number of elements over the thickness of the flange
    ewf = 8 #number of elements over half of the width of the flange
    er = 6  #number of elements in the circular segment

    Body = simple.rectangle(etb,ehb,tw/2.,h/2.-tf-r)
    Flange1 =  simple.rectangle(er/2,etf-etb,tw/2.+r,tf-tw/2.).translate([0.,h/2.-(tf-tw/2.),0.])
    Flange2 =  simple.rectangle(ewf,etf-etb,b/2.-r-tw/2.,tf-tw/2.).translate([tw/2.+r,h/2.-(tf-tw/2.),0.])
    Flange3 =  simple.rectangle(ewf,etb,b/2.-r-tw/2.,tw/2.).translate([tw/2.+r,h/2.-tf,0.])
    c1a = simple.line([0,h/2-tf-r,0],[0,h/2-tf+tw/2,0],er/2)
    c1b = simple.line([0,h/2-tf+tw/2,0],[tw/2+r,h/2-tf+tw/2,0],er/2)
    c1 = c1a + c1b
    c2 = simple.circle(90./er,0.,90.).reflect(0).scale(r).translate([tw/2+r,h/2-tf-r,0])
    Filled = simple.connectCurves(c2,c1,etb)
    Quarter = Body + Filled + Flange1 + Flange2 + Flange3
    Half = Quarter + Quarter.reflect(1).reverse()
    Full = Half + Half.reflect(0).reverse()
    Section = Full.toMesh()

    clear()
    draw(Section,color=red)
    #return

    #pause()

    method = ask("Choose extrude method:",['Cancel','Sweep','Connect','Extrude','ExtrudeQuadratic','Revolve','RevolveLoop'])

    import timer
    t = timer.Timer()
    if method == 'Sweep':
        L = simple.line([0,0,0],[0,0,l],el)
        x = concatenate([L.coords[:,0],L.coords[-1:,1]])
        path = curve.PolyLine(x)
        Beam = Section.sweep(path,normal=[0.,0.,1.],upvector=[0.,1.,0.])

    elif method == 'Connect':
        Section1 = Section.trl([0,0,l])
        Beam = Section.connect(Section1,el)

    elif method == 'Extrude':
        Beam = Section.extrude(el,step=l/el,dir=2)

    elif method == 'ExtrudeQuadratic':
        Section = Section.convert('quad9')
        Beam = Section.extrude(el,step=l/el,dir=2,degree=2)

    elif method == 'Revolve':
        Beam = Section.revolve(el,axis=1,angle=60.,around=[-l,0.,0.])

    elif method == 'RevolveLoop':
        Beam = Section.revolve(el,axis=1,angle=240.,around=[-l,0.,0.],loop=True)

    else:
        return

    print("Computing: %s seconds" % t.seconds())
    #print Beam.prop
    #print Beam.elems.shape

    t.reset()
    clear()
    #draw(Beam,color='red',linewidth=2)
    draw(Beam.getBorderMesh(),color='red',linewidth=2)
    print("Drawing: %s seconds" % t.seconds())
    export({'Beam':Beam})

view('myview1',True)

drawtimeout = 1
for i in range(19):
    createView('myview2',(i*10.,20.,0.))
    view('myview2',True)

# fly tru
if ack("Do you want to fly through the structure?"):
    totaltime = 10
    nsteps = 50
    # make sure bottom iz at y=0 and structure is centered in (x,z)
    F = F.centered()
    F = F.translate(1,-F.bbox()[0,1])
    clear()
    linewidth(1)
    draw(F)
    bb = F.bbox()
    # create a bottom plate
    B = simple.rectangle(1,1).swapAxes(1,2).centered().scale(F.sizes()[0]*1.5)
    smooth()
    draw(B)
    # Fly at reasonable height
    bb[0,1] = bb[1,1] = 170.
    ends = interpolate(Formex([[bb[0]]]),Formex([[bb[1]]]),[-0.5,0.6])
    path = connect([ends,ends],bias=[0,1]).divide(nsteps)
    linewidth(2)
    draw(path)
    steptime = float(totaltime)/nsteps
    flyAlong(path,sleeptime=steptime)

"""WedgeHex

level = 'normal'
topics = ['mesh']
techniques = ['revolve','degenerate'] 
"""

import simple

clear()
smoothwire()

# create a 2D xy mesh
nx,ny = 6,2
G = simple.rectangle(1,1,1.,1.).replic2(nx,ny)
M = G.toMesh()
draw(M, color='red')
view('iso')

# create a 3D axial-symmetric mesh by REVOLVING
n,a = 8,45.
R = M.revolve(n,angle=a,axis=1,around=[1.,0.,0.])
sleep(2)
draw(R,color='yellow')

# reduce the degenerate elements to WEDGE6
sleep(2)
clear()
ML = R.splitDegenerate()
ML = [ Mi.setProp(i) for i,Mi in enumerate(ML) ]

def run():
    reset()
    clear()

    # Property numbers used
    pbol = 1  # Bol
    ptop = 2  # Top plate
    pbot = 3  # Bottom plate

    scale = 15.   # scale (grid unit in mm)

    # Create a solid sphere
    BolSurface = simple.sphere().scale(scale)
    try:
        # tetmesh may not be available
        Bol = BolSurface.tetmesh(quality=True).setProp(pbol)
    except:
        return
    draw(Bol)

    # Create top and bottom plates
    plate = simple.rectangle(4,4).toMesh().centered()
    topplate = plate.setProp(ptop).trl(2,1.).scale(scale)
    botplate = plate.setProp(pbot).trl(2,-1.).scale(scale)
    draw([topplate,botplate])

    # model is completely drawn, keep fixed bbox
    setDrawOptions({'bbox':'last','marksize':8})

    # Assemble the model
    M = Model(meshes=[Bol,topplate,botplate])

    # Create the property database
    P = PropertyDB()

    # In this simple example, we do not use a material/section database,
    # but define the data directly
    steel = {
        'name': 'steel',
        'young_modulus': 207000,
        'poisson_ratio': 0.3,
        'density': 7.85e-9,
        'plastic' : [
            (305.45,       0.),
            (306.52, 0.003507),
            (308.05, 0.008462),
            (310.96,  0.01784),
            (316.2, 0.018275),
            (367.5, 0.047015),
            (412.5, 0.093317),
            (448.11, 0.154839),
            (459.6, 0.180101),
            (494., 0.259978),
            (506.25, 0.297659),
            (497., 0.334071),
            (482.8, 0.348325),
            (422.5, 0.366015),
            (399.58,   0.3717),
            (1.,  0.37363),
            ],
        }
    solid_steel = {
        'name': 'solid_steel',
        'sectiontype': 'solid',
        'material': 'steel',   # Need material reference for Abaqus
        }
    steel_plate = {
        'name': 'solid_steel',
        'sectiontype': 'solid',
        'thickness': 3,
        'material': 'steel',   # Need material reference for Abaqus
        }

    # Set the element properties
    eset = dict([(p,where(M.prop==p)[0]) for p in [pbol,ptop,pbot]])

    # Bol is elasto/plastic
    P.elemProp(set=eset[pbol],name='Bol',eltype='C3D4',section=ElemSection(section=solid_steel,material=steel))

    # Top plate is rigid or elasto-plastic
    topplate_rigid = True
    if topplate_rigid:
        # Rigid bodies need a reference node.
        # We select the most central node, but any node would also work,
        # e.g. pbref = M.elems[1][0][0], the very first node in the group
        reftop = groupCentralPoint(M,1)
        print("Top plate refnode: %s" % reftop)
        draw(M.coords[reftop],color=green)
        P.elemProp(set=eset[ptop],name='TopPlate',eltype='R3D4',section=ElemSection(sectiontype='rigid',refnode=reftop))
    else:
        P.elemProp(set=eset[ptop],name='TopPlate',eltype='CPS4',section=ElemSection(section=steel_plate,material=steel))


    # Bottom plate is rigid or elasto-plastic
    refbot = groupCentralPoint(M,2)
    print("Bottom plate refnode: %s" % refbot)
    draw(M.coords[refbot],color=blue)
    P.elemProp(set=eset[pbot],name='BottomPlate',eltype='R3D4',section=ElemSection(sectiontype='rigid',refnode=refbot))

    # Set the boundary conditions
#.........这里部分代码省略.........

h = 200. #beam height
b = 200. #flange width 
tf = 15. #flange thickness
tw = 9.  #body thickness
l = 400. #beam length
r = 18.  #filling radius

# MESH PARAMETERS
el = 20 #number of elements along the length
etb = 2 #number of elements over half of the thickness of the body
ehb = 5 #number of elements over half of the height of the body
etf = 5 #number of elements over the thickness of the flange
ewf = 8 #number of elements over half of the width of the flange
er = 6  #number of elements in the circular segment

Body = simple.rectangle(etb,ehb,tw/2.,h/2.-tf-r)
Flange1 =  simple.rectangle(er/2,etf-etb,tw/2.+r,tf-tw/2.).translate([0.,h/2.-(tf-tw/2.),0.])
Flange2 =  simple.rectangle(ewf,etf-etb,b/2.-r-tw/2.,tf-tw/2.).translate([tw/2.+r,h/2.-(tf-tw/2.),0.])
Flange3 =  simple.rectangle(ewf,etb,b/2.-r-tw/2.,tw/2.).translate([tw/2.+r,h/2.-tf,0.])
c1a = simple.line([0,h/2-tf-r,0],[0,h/2-tf+tw/2,0],er/2)
c1b = simple.line([0,h/2-tf+tw/2,0],[tw/2+r,h/2-tf+tw/2,0],er/2)
c1 = c1a + c1b
c2 = simple.circle(90./er,0.,90.).reflect(0).scale(r).translate([tw/2+r,h/2-tf-r,0])
Filled = simple.connectCurves(c2,c1,etb)
Quarter = Body + Filled + Flange1 + Flange2 + Flange3
Half = Quarter + Quarter.reflect(1).reverse()
Full = Half + Half.reflect(0).reverse()
Section = Full.toMesh()

clear()
draw(Section,color=red)

def revolve_QuadMesh(n0, e0, nr, ang=None):
    """it takes a Quad mesh on xy and revolve it along the z axis nr times by ang. If ang==None, then it is calculated in order to fill 360 degrees with the nr revolutions."""
    if ang==None: ang=360./nr
    for i in range(int(nr)):
        n1=Formex(n0).rotate(-ang, 1)[:].reshape(-1, 3)
        n, e=connectMesh(n0,n1,e0)
        n0=n1.copy()
        parts.append(Formex(n[e], eltype='Hex8'))
    femodels = [part.feModel() for part in parts]
    nodes,elems = mergeModels(femodels)
    elems=concatenate([elems], 0).reshape(-1, 8)
    return nodes, elems
clear()
#create a 2D xy mesh
n=4
G=simple.rectangle(1,1,1.,1.).replic(n,1.,dir=1).replic(n,1.,dir=0)
draw(G, color='red')
view('front')
sleep(1)
#create a 3D axial-symmetric mesh by REVOLVING
n0, e0=G.feModel()
parts=[]



n, e=revolve_QuadMesh(n0, e0, nr=4, ang=20)
C=Formex(n[e], eltype='Hex8')

#check if there are Wedge elements in the global Hex mesh
w, h= detectHex2Wedge(C)
W=Formex(w, eltype='Wedge6')

"""

import simple
from mesh import Mesh

clear()
smoothwire()

nx = 4
ny = 3
nz = 7

delay(2)

# A rectangular mesh
M1 = simple.rectangle(nx,ny).toMesh().setProp(1)
# Same mesh, rotated and translated
M2 = M1.rotate(45,0).translate([1.,-1.,nz]).setProp(3)
draw([M1,M2])

# Leave out the first and the last two elements
sel = arange(M1.nelems())[1:-2]
m1 = M1.select(sel)
m2 = M2.select(sel)
clear()
draw([m1,m2],view=None)

# Connect both meshes to a hexaeder mesh
m = m1.connect(m2,nz)
clear()
draw(m,color=red,view=None)