> On 20 Mar 2019, at 08:01, Christophe Geuzaine <[email protected]> wrote:
>
>
> Hi Laurent,
>
> It's doable, but since you would need to connect this discrete surface with
> new CAD surfaces to build a volume, it would be easier if the input mesh had
> a topology, i.e. if you provided the surface + the bounding curves + points
> for the 4 corners.
>
> Then new surfaces and a volume could be constructed, connected to these
> curves/points, and all the standard meshing algorithms could be used.
>
I've played a bit with your file and you can actually do all of this directly
with the Gmsh api (here in Python):
import gmsh
import sys
gmsh.initialize(sys.argv)
gmsh.option.setNumber("General.Terminal", 1)
gmsh.option.setNumber("Mesh.Algorithm", 6) # frontal
gmsh.option.setNumber("Mesh.Algorithm3D", 10) # hxt
gmsh.open("AskerveinHill.stl")
# compute boundary
edges = gmsh.model.mesh.getElementEdgeNodes(2)
bnd = set()
for i in range(0, len(edges), 2):
e = (edges[i], edges[i+1])
ne = (edges[i+1], edges[i])
if(e in bnd):
bnd.remove(e)
elif(ne in bnd):
bnd.remove(ne)
else:
bnd.add(e)
# find corner nodes
nodes, coord, _ = gmsh.model.mesh.getNodes()
def findNode(x, y, coord):
for i in range(len(nodes)):
if(x == coord[3*i] and y == coord[3*i+1]):
return i
return -1
i1 = findNode(0, 0, coord)
i2 = findNode(2649, 0, coord)
i3 = findNode(2649, 2491, coord)
i4 = findNode(0, 2491, coord)
corners = {nodes[i1], nodes[i2], nodes[i3], nodes[i4]}
# create 4 point elements for the 4 corners
p1 = gmsh.model.geo.addPoint(coord[3*i1], coord[3*i1+1], coord[3*i1+2])
p2 = gmsh.model.geo.addPoint(coord[3*i2], coord[3*i2+1], coord[3*i2+2])
p3 = gmsh.model.geo.addPoint(coord[3*i3], coord[3*i3+1], coord[3*i3+2])
p4 = gmsh.model.geo.addPoint(coord[3*i4], coord[3*i4+1], coord[3*i4+2])
gmsh.model.geo.synchronize()
gmsh.model.mesh.setElementsByType(p1, 15, [], [nodes[i1]])
gmsh.model.mesh.setElementsByType(p2, 15, [], [nodes[i2]])
gmsh.model.mesh.setElementsByType(p3, 15, [], [nodes[i3]])
gmsh.model.mesh.setElementsByType(p4, 15, [], [nodes[i4]])
# compute ordered list of boundary line elements, for each side
input = list(bnd)
input_dict = dict(input)
elem = nodes[i1]
lines = []
for _ in range(len(input)):
if(elem in corners): lines.append([])
lines[-1].append((elem, input_dict[elem]))
elem = input_dict[elem]
# create line elements in the model
l1 = gmsh.model.addDiscreteEntity(1)
gmsh.model.mesh.setElementsByType(l1, 1, [], list(sum(lines[0], ())))
l2 = gmsh.model.addDiscreteEntity(1)
gmsh.model.mesh.setElementsByType(l2, 1, [], list(sum(lines[1], ())))
l3 = gmsh.model.addDiscreteEntity(1)
gmsh.model.mesh.setElementsByType(l3, 1, [], list(sum(lines[2], ())))
l4 = gmsh.model.addDiscreteEntity(1)
gmsh.model.mesh.setElementsByType(l4, 1, [], list(sum(lines[3], ())))
# reclassify mesh nodes in the correct entities
gmsh.model.mesh.reclassifyNodes()
# create 4 new points to create a box
z = 1000
p5 = gmsh.model.geo.addPoint(0, 0, z)
p6 = gmsh.model.geo.addPoint(2649, 0, z)
p7 = gmsh.model.geo.addPoint(2649, 2491, z)
p8 = gmsh.model.geo.addPoint(0, 2491, z)
# create new curves and surfaces
l5 = gmsh.model.geo.addLine(p1, p5)
l6 = gmsh.model.geo.addLine(p2, p6)
l7 = gmsh.model.geo.addLine(p3, p7)
l8 = gmsh.model.geo.addLine(p4, p8)
l9 = gmsh.model.geo.addLine(p5, p6)
l10 = gmsh.model.geo.addLine(p6, p7)
l11 = gmsh.model.geo.addLine(p7, p8)
l12 = gmsh.model.geo.addLine(p8, p5)
cl1 = gmsh.model.geo.addCurveLoop([l1, l6, l9, l5])
s1 = gmsh.model.geo.addPlaneSurface([cl1])
cl2 = gmsh.model.geo.addCurveLoop([l2, l6, l10, l7])
s2 = gmsh.model.geo.addPlaneSurface([cl2])
cl3 = gmsh.model.geo.addCurveLoop([l3, l7, l11, l8])
s3 = gmsh.model.geo.addPlaneSurface([cl3])
cl4 = gmsh.model.geo.addCurveLoop([l4, l8, l12, l5])
s4 = gmsh.model.geo.addPlaneSurface([cl4])
cl5 = gmsh.model.geo.addCurveLoop([l9, l10, l11, l12])
s5 = gmsh.model.geo.addPlaneSurface([cl5])
sl1 = gmsh.model.geo.addSurfaceLoop([1, s1, s2, s3, s4, s5])
v1 = gmsh.model.geo.addVolume([sl1])
gmsh.model.geo.synchronize()
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 50)
gmsh.fltk.run()
Christophe
> Christophe
>
>> On 19 Mar 2019, at 21:54, Laurent BRICTEUX <[email protected]>
>> wrote:
>>
>> Dear all,
>>
>> Is it possible to obtain a volume mesh starting from an stl terrain surface
>> (see attachment)?
>> I would like to make a CFD mesh consisting in a volume box for which the
>> bottom surface is a prescribed terrain.
>> The elements should be tets only.
>> Is there any example geo file that could help to solve my problem ?
>>
>> Thank you for your help and the nice piece of soft.
>>
>> Laurent
>>
>>
>>
>>
>> <AskerveinHill.stl.zip>_______________________________________________
>> gmsh mailing list
>> [email protected]
>> http://onelab.info/mailman/listinfo/gmsh
>
> —
> Prof. Christophe Geuzaine
> University of Liege, Electrical Engineering and Computer Science
> http://www.montefiore.ulg.ac.be/~geuzaine
>
>
>
>
> _______________________________________________
> gmsh mailing list
> [email protected]
> http://onelab.info/mailman/listinfo/gmsh
—
Prof. Christophe Geuzaine
University of Liege, Electrical Engineering and Computer Science
http://www.montefiore.ulg.ac.be/~geuzaine
_______________________________________________
gmsh mailing list
[email protected]
http://onelab.info/mailman/listinfo/gmsh
