Dear Miklos, Thanks a lot for getting back to me with this detailed response. I am looking to create a pipe geometry (with an annulus cross section) and possibly connections like T-joint, elbow joint etc in the future. I was particularly interested in using Gmsh because of their ability to generate boundary markers of curves, surfaces etc using Physical Groups. Do you have a suggestion on how I can do that in Firedrake? On a sidenote: I tried to create a hexahedral Q element as tensor product of Q quadilateral and Q interval. Firedrake said Q is supported but that this was not possible. I managed to create a Q element by using the CG family and then using GLL quadrature rule in dx(gllquad), similar to the Sum Factorization Jupyter Notebook example. What is the difference then between CG and Q family then? -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran@netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:72b8d1cc-99f4-4732-8629-c9716ea41ac1] ________________________________ From: Miklós Homolya <miklos.homolya@tum.de> Sent: Wednesday, May 5, 2021 5:26 AM To: Venketeswaran, Abhishek (FN) <Abhishek.Venketeswaran@netl.doe.gov>; firedrake@imperial.ac.uk <firedrake@imperial.ac.uk> Cc: Abhishek V <v.abhi6@gmail.com> Subject: [EXTERNAL] Re: [firedrake] KeyError using a Gmsh generated hexahedral mesh Dear Abhishek, As far as I know, fully unstructured hexahedral meshes, as generated by gmsh, are not supported yet in Firedrake. Firedrake can work with extruded quadrilateral meshes, whose cells are indeed hexahedra, but the mesh is structured in one dimension, while the base quadrilateral mesh can be fully unstructured. What kind of geometry do you need? Kind regards, Miklos "Venketeswaran, Abhishek (FN)" <Abhishek.Venketeswaran@netl.doe.gov> writes:
Hi, I successfully generated a hexahedral mesh using Gmsh. Please find the code in the .txt file attached here.
I then tried to use the mesh in Firedrake. Please find below the code: from firedrake import * mesh=Mesh('pipev1.msh')
I get a KeyError:6 from ufl.Cell module. Kindly help me.
Thanks
-- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran@netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:d5f4824f-4769-4311-b539-64bf8a6d05d2] # ------------------------------------------------------------------------------ # # Gmsh Mesh of a cylindrical pipe # # Dimensions: # # ------------------------------------------------------------------------------
# The Python API is entirely defined in the `gmsh.py' module (which contains the # full documentation of all the functions in the API): import gmsh import sys import numpy as np from math import isclose # Before using any functions in the Python API, Gmsh must be initialized: gmsh.initialize()
# Next we add a new model named "pipe1" (if gmsh.model.add() is not called a new # unnamed model will be created on the fly, if necessary): gmsh.model.add("pipev1")
# We start by defining some points and some lines. To make the code shorter we # can redefine a namespace: fG = gmsh.model.geo
# Pipe Dimensions Ri=7.0 # in Ro=8.0 # in ft2in=12.0 Lz=4.0*ft2in lin=0.5 # Target elements lout=0.5 nz=50 ############################################################################### # Adapted from t1.py # The first type of `elementary entity' in Gmsh is a `Point'. To create a point # with the built-in CAD kernel, the Python API function is # gmsh.model.geo.addPoint(): # - the first 3 arguments are the point coordinates (x, y, z) # - the next (optional) argument is the target mesh size close to the point # - the last (optional) argument is the point tag (a stricly positive integer # that uniquely identifies the point) # 5 points to define the circle origin and circumference rad=Ri cin_pts=[] cin_pts.append(fG.addPoint(0.0,0.0,0.0,lin)) cin_pts.append(fG.addPoint(rad,0.0,0.0,lin)) cin_pts.append(fG.addPoint(-rad,0.0,0.0,lin)) cin_pts.append(fG.addPoint(0.0,rad,0.0,lin)) cin_pts.append(fG.addPoint(0.0,-rad,0.0,lin))
# 5 points to define the outer circle rad=Ro cout_pts=[] cout_pts.append(fG.addPoint(0.0,0.0,0.0,lout)) cout_pts.append(fG.addPoint(rad,0.0,0.0,lout)) cout_pts.append(fG.addPoint(-rad,0.0,0.0,lout)) cout_pts.append(fG.addPoint(0.0,rad,0.0,lout)) cout_pts.append(fG.addPoint(0.0,-rad,0.0,lout))
# Curves are Gmsh's second type of elementery entities, and, amongst curves, # straight lines are the simplest. The API to create straight line segments with # the built-in kernel follows the same conventions: the first 2 arguments are # point tags (the start and end points of the line), and the last (optional one) # is the line tag. # Note that curve tags are separate from point tags - hence we can reuse tag `1' # for our first curve. And as a general rule, elementary entity tags in Gmsh # have to be unique per geometrical dimension.
# Define 4 arcs of the outer circle cout_arcs=[] cout_arcs.append(fG.addCircleArc(cout_pts[3],cout_pts[0] , cout_pts[1])) cout_arcs.append(fG.addCircleArc(cout_pts[1],cout_pts[0] , cout_pts[4])) cout_arcs.append(fG.addCircleArc(cout_pts[4],cout_pts[0] , cout_pts[2])) cout_arcs.append(fG.addCircleArc(cout_pts[2],cout_pts[0] , cout_pts[3])) #print(cout_arcs) # Define 4 arcs of the inner circle cin_arcs=[] cin_arcs.append(fG.addCircleArc(cin_pts[3],cin_pts[0] , cin_pts[1])) cin_arcs.append(fG.addCircleArc(cin_pts[1],cin_pts[0] , cin_pts[4])) cin_arcs.append(fG.addCircleArc(cin_pts[4],cin_pts[0] , cin_pts[2])) cin_arcs.append(fG.addCircleArc(cin_pts[2],cin_pts[0] , cin_pts[3])) #print(cin_arcs) # Define 4 lines b/w circles clines=[fG.addLine(cout_pts[index], cin_pts[index]) for index in range(1,5)] #print(clines)
# The third elementary entity is the surface. In order to define a simple # circular surface from the four curves defined above, a curve loop has first # to be defined. A curve loop is defined by an ordered list of connected curves, # a sign being associated with each curve (depending on the orientation of the # curve to form a loop). The API function to create curve loops takes a list # of integers as first argument, and the curve loop tag (which must be unique # amongst curve loops) as the second (optional) argument:
quads=[] quads.append(fG.addCurveLoop([cout_arcs[0],clines[0],-cin_arcs[0],-clines[2] ])) quads.append(fG.addCurveLoop([cout_arcs[1],clines[3],-cin_arcs[1],-clines[0] ]) ) quads.append(fG.addCurveLoop([cout_arcs[2],clines[1],-cin_arcs[2],-clines[3] ])) quads.append(fG.addCurveLoop([cout_arcs[3],clines[2],-cin_arcs[3],-clines[1] ]))
surfs=[] extsurfs=[] vols=[] extvec=[0.0,0.0,Lz] for quad in quads: circquad=fG.addPlaneSurface([quad]) cylquad=fG.extrude([(2,circquad)],extvec[0],extvec[1],extvec[2],[nz],[1],recombine=True) surfs.append(circquad) extsurfs.append(cylquad[0]) vols.append(cylquad[1]) fG.synchronize()
def avgdist(taglist,origin): dist=0.0 for e in taglist: nodeCoords = gmsh.model.getValue(e[0], e[1],[]) dist=dist+np.linalg.norm(nodeCoords-origin) return dist/len(taglist) coutext_arcs=[] cinext_arcs=[] for surf in extsurfs: topsurfcurves=gmsh.model.getBoundary([surf]) for c in topsurfcurves: bdypts=gmsh.model.getBoundary([c]) dist=avgdist(bdypts,extvec) if isclose(dist,Ro): coutext_arcs.append(c[1]) elif isclose(dist,Ri): cinext_arcs.append(c[1]) print(coutext_arcs) print(cinext_arcs) #entities = gmsh.model.getEntities() #for e in entities: # print("Entity " + str(e) + " of type " + gmsh.model.getType(e[0], e[1]))
# All the extrusion functions return a list of (dim,tag) tuples of extruded entities: the # "top" of the extruded surface (in `[0]'), the newly created volume (in # `[1]') and the lateral surfaces (in `[2]', `[3]', ...). #print(surfs) #print(vquads)
# Before they can be meshed (and, more generally, before they can be used by API # functions outside of the built-in CAD kernel functions), the CAD entities must # be synchronized with the Gmsh model, which will create the relevant Gmsh data # structures. This is achieved by the gmsh.model.geo.synchronize() API call for # the built-in CAD kernel. Synchronizations can be called at any time, but they # involve a non trivial amount of processing; so while you could synchronize the # internal CAD data after every CAD command, it is usually better to minimize # the number of synchronization points.
fG.synchronize()
# At this level, Gmsh knows everything to display the pipe and # to mesh it. An optional step is needed if we want to group elementary # geometrical entities into more meaningful groups, e.g. to define some # mathematical ("domain", "boundary"), functional ("left wing", "fuselage") or # material ("steel", "carbon") properties. # # Such groups are called "Physical Groups" in Gmsh. By default, if physical # groups are defined, Gmsh will export in output files only mesh elements that # belong to at least one physical group. (To force Gmsh to save all elements, # whether they belong to physical groups or not, set the `Mesh.SaveAll' option # to 1.) Physical groups are also identified by tags, i.e. stricly positive # integers, that should be unique per dimension (0D, 1D, 2D or 3D). Physical # groups can also be given names. # # Here we define a physical curve that groups the arcs of inner and outer # circles into individual groups : loop_cin=fG.addCurveLoop(cin_arcs) # Inner Circle loop_cout=fG.addCurveLoop(cout_arcs) # Outer Circle # Physical curve groups are started from 1000 # Physical surfaces are started from 2000 # Physical volumes are started from 3000 gmsh.model.addPhysicalGroup(1, cout_arcs,1001) gmsh.model.addPhysicalGroup(1, cin_arcs,1002) gmsh.model.addPhysicalGroup(1, coutext_arcs,1003) gmsh.model.addPhysicalGroup(1, cinext_arcs,1004)
surftags=surfs extsurftags=[t[1] for t in extsurfs] voltags=[v[1] for v in vols]
gmsh.model.addPhysicalGroup(2, surftags,2001) gmsh.model.addPhysicalGroup(2, extsurftags,2002) gmsh.model.addPhysicalGroup(3, voltags,3001)
#ps = gmsh.model.addPhysicalGroup(2, [1]) #gmsh.model.setPhysicalName(2, ps, "My surface")
################################################################################# # Adapted from t11.py # To generate quadrangles instead of triangles, we can simply add #gmsh.model.mesh.setRecombine(2, pl) # If we'd had several surfaces, we could have used the global option gmsh.option.setNumber("Mesh.RecombineAll", 1)
# The default recombination algorithm is called "Blossom": it uses a minimum # cost perfect matching algorithm to generate fully quadrilateral meshes from # triangulations. More details about the algorithm can be found in the # following paper: J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise, # A. Johnen and C. Geuzaine, "Blossom-Quad: a non-uniform quadrilateral mesh # generator using a minimum cost perfect matching algorithm", International # Journal for Numerical Methods in Engineering 89, pp. 1102-1119, 2012.
# For even better 2D (planar) quadrilateral meshes, you can try the # experimental "Frontal-Delaunay for quads" meshing algorithm, which is a # triangulation algorithm that enables to create right triangles almost # everywhere: J.-F. Remacle, F. Henrotte, T. Carrier-Baudouin, E. Bechet, # E. Marchandise, C. Geuzaine and T. Mouton. A frontal Delaunay quad mesh # generator using the L^inf norm. International Journal for Numerical Methods # in Engineering, 94, pp. 494-512, 2013. Uncomment the following line to try # the Frontal-Delaunay algorithms for quads: # gmsh.option.setNumber("Mesh.MshFileVersion", 2) gmsh.option.setNumber("Mesh.Algorithm", 8)
# You can also set the subdivision step alone, with # 0=none, 1=all quadrangles, 2=all hexahedra gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 2) gmsh.model.mesh.generate(3)
# Note that you could also apply the recombination algorithm and/or the # subdivision step explicitly after meshing, as follows: # # gmsh.model.mesh.generate(2) # gmsh.model.mesh.recombine() # gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) # gmsh.model.mesh.refine()
# ... and save it to disk gmsh.write("pipev1.msh") ########################################################################
# To visualize the model we can run the graphical user interface with # `gmsh.fltk.run()'. Here we run it only if the "-nopopup" is not provided in # the command line arguments: if '-nopopup' not in sys.argv: gmsh.fltk.run()
gmsh.finalize()
#ToDo # Check t6.py and use setTransfiniteCurve options to customize mesh # Look into Mesh Smoothing option in t6.py _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
******************************************************************** This message does not originate from a known Department of Energy email system. Use caution if this message contains attachments, links or requests for information. ********************************************************************
