pycad/py_src/gmsh.py

# -*- coding: utf-8 -*-

##############################################################################
#
# Copyright (c) 2003-2020 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
# Development from 2019 by School of Earth and Environmental Sciences
#
##############################################################################

from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2020 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

"""
mesh generation using gmsh

:var __author__: name of author
:var __copyright__: copyrights
:var __license__: licence agreement
:var __url__: url entry point on documentation
:var __version__: version
:var __date__: date of the version
"""

__author__="Lutz Gross, l.gross@uq.edu.au"

from . import design
import tempfile
import os
from .primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
from esys.escript import getMPIWorldMax, getMPIRankWorld, gmshGeo2Msh
from .transformations import DEG

class Design(design.AbstractDesign):
    """
    Design for gmsh.
    """
    DELAUNAY="Delauny"
    MESHADAPT="MeshAdapt"
    FRONTAL="Frontal"
    NETGEN="Frontal"
    TETGEN="Delauny"

    def __init__(self, dim=3, element_size=1., order=1, keep_files=False):
        """
        Initializes the gmsh design.

        :param dim: spatial dimension
        :param element_size: global element size
        :param order: element order
        :param keep_files: flag to keep work files
        """
        design.AbstractDesign.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
        self.__mshname_set = False
        self.__scriptname=""
        self.setScriptFileName()
        self.setOptions()
        self.setFileFormat(self.GMSH)

    def setScriptFileName(self, name=None):
        """
        Sets the filename for the gmsh input script. If no name is given a name
        with extension `geo` is generated.
        """
        if self.__scriptname:
            os.unlink(self.__scriptname)
        if name == None:
            self.__scriptname_set=False
            tmp_f_id=tempfile.mkstemp(suffix=".geo")
            self.__scriptname=tmp_f_id[1]
            os.close(tmp_f_id[0])
        else:
            self.__scriptname=name
            self.__scriptname_set=True

    def getScriptFileName(self):
        """
        Returns the name of the gmsh script file.
        """
        return self.__scriptname

    def setOptions(self, algorithm=None, optimize_quality=True, smoothing=1,
                   curvature_based_element_size=False, algorithm2D=None,
                   algorithm3D=None, generate_hexahedra=False,
                   random_factor=None):
        """
        Sets options for the mesh generator.

        :param algorithm: selects 2D meshing algorithm
        :type algorithm: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
        :param algorithm2D: must be equal to algorithm
        :type algorithm2D: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
        :param algorithm3D: selects 3D meshing algorithm
        :type algorithm3D: in self.DELAUNAY, self.FRONTAL
        :param curvature_based_element_size: switch for curvature based element size adaption
        :type curvature_based_element_size: ```bool```
        :param smoothing: number of smoothing steps
        :type smoothing: non-negative ```int```
        :param optimize_quality: switch for mesh quality optimization
        :type optimize_quality: ```bool```
        :param generate_hexahedra: switch for using quadrangles/hexahedra elements everywhere.
        :type generate_hexahedra: ```bool```
        :param random_factor: used in the 2D meshing algorithm (should be increased if RandomFactor * size(triangle)/size(model) approaches machine accuracy)
        :type random_factor: positive ```float```
        """
        if random_factor==None: random_factor=1.e-9
        if not random_factor > 0:
            raise ValueError("random_factor must be positive.")
        smoothing=int(smoothing)
        if not smoothing > 0:
            raise ValueError("smoothing must be positive.")

        if algorithm3D is None:
            algorithm3D=self.FRONTAL
        if algorithm is None:
            if algorithm2D is None:
                algorithm2D=self.MESHADAPT
        else:
            if not algorithm2D is None:
                if not algorithm == algorithm2D:
                    raise ValueError("argument algorithm (=%s) and algorithm2D (=%s) must have the same value if set."%(algorithm, algorithm2D))
            algorithm2D = algorithm
        if not algorithm2D in [ self.DELAUNAY, self.MESHADAPT, self.FRONTAL ]:
            raise ValueError("illegal 2D meshing algorithm %s."%algorithm2D)
        if not algorithm3D in [ self.DELAUNAY, self.FRONTAL ]:
            raise ValueError("illegal 3D meshing algorithm %s."%algorithm3D)

        self.__curvature_based_element_size=curvature_based_element_size
        self.__algo2D=algorithm2D
        self.__algo3D=algorithm3D
        self.__optimize_quality=optimize_quality
        self.__smoothing=smoothing
        self.__generate_hexahedra=generate_hexahedra
        self.__random_factor=random_factor

    def getOptions(self, name=None):
        """
        Returns the current options for the mesh generator.
        """
        if name is None:
            return {"optimize_quality" : self.__optimize_quality ,
                    "smoothing" : self.__smoothing,
                    "curvature_based_element_size" : self.__curvature_based_element_size,
                    "generate_hexahedra" : self.__generate_hexahedra,
                    "algorithm2D" : self.__algo2D,
                    "algorithm3D" : self.__algo3D ,
                    "random_factor" : self.__random_factor }
        else:
            return self.getOption()[name]

    def __del__(self):
        """
        Cleans up.
        """
        try:
            if not self.keepFiles():
                if not self.__scriptname_set: #i.e. it's a tempfile
                    os.unlink(self.getScriptFileName())
                if not self.__mshname_set: #i.e. it's a tempfile
                    os.unlink(self.getMeshFileName())
        except OSError:
            pass # The file might not have been created and there is nothing
                 # to do about a "failure" here anyway

    def getScriptHandler(self):
        """
        Returns a handler to the script file to generate the geometry.
        In the current implementation a script file name is returned.
        """
        if getMPIRankWorld() == 0:
            open(self.getScriptFileName(),"w").write(self.getScriptString())
        return self.getScriptFileName()

    def getMeshHandler(self):
        """
        Returns a handle to a mesh meshing the design. In the current
        implementation a mesh file name in gmsh format is returned.
        """

        verbosity = 3
        ret = gmshGeo2Msh(self.getScriptHandler(), self.getMeshFileName(),
                          self.getDim(), self.getElementOrder(), verbosity)
        if ret > 0:
            self.setKeepFilesOn() #no files to delete, so don't try to
            raise RuntimeError("Could not build mesh using gmsh.\n" + \
                    "Is gmsh available?")
        return self.getMeshFileName()

    def getScriptString(self):
        """
        Returns the gmsh script to generate the mesh.
        """
        h=self.getElementSize()
        out='// generated by esys.pycad\nGeneral.Terminal = 1;\nGeneral.ExpertMode = 1;\n'
        options=self.getOptions()
        if options["optimize_quality"]:
            out += "Mesh.Optimize = 1;\n"
        else:
            out += "Mesh.Optimize = 0;\n"

        if options["curvature_based_element_size"]:
            out += "Mesh.CharacteristicLengthFromCurvature = 1;\n"
        else:
            out += "Mesh.CharacteristicLengthFromCurvature = 0;\n"

        if options["generate_hexahedra"]:
            if self.getDim() == 2:
                out += "Mesh.SubdivisionAlgorithm = 1;\n"
            else:
                out += "Mesh.SubdivisionAlgorithm = 2;\n"
        else:
            out += "Mesh.SubdivisionAlgorithm = 0;\n"

        out += "Mesh.Smoothing = %d;\n"%options["smoothing"]
        out += "Mesh.RandomFactor = %.14e;\n"%options["random_factor"]
        if options["algorithm2D"] == self.MESHADAPT:
            out += "Mesh.Algorithm = 1; // = MeshAdapt\n"
        elif options["algorithm2D"] == self.DELAUNAY:
            out += "Mesh.Algorithm = 5; // = Delaunay\n"
        elif options["algorithm2D"] == self.FRONTAL:
            out += "Mesh.Algorithm = 6; // = Frontal\n"

        if options["algorithm3D"] == self.DELAUNAY:
            out += "Mesh.Algorithm3D = 1; // = Delaunay\n"
        elif options["algorithm3D"] == self.FRONTAL:
            out += "Mesh.Algorithm3D = 4; // = Frontal\n"

        for prim in self.getAllPrimitives():
            p=prim.getUnderlyingPrimitive()
            if isinstance(p, Point):
                c=p.getCoordinates()
                #out+="Point(%s) = {%f , %f, %f , %f };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
                out += "Point(%s) = {%.14e, %.14e, %.14e, %.14e};\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)

            elif isinstance(p, Spline):
                out += "Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)

            elif isinstance(p, BezierCurve):
                out += "Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)

            elif isinstance(p, BSpline):
                out += "BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)

            elif isinstance(p, Line):
                out += "Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)

            elif isinstance(p, Arc):
                out += "Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)

            elif isinstance(p, Ellipse):
                out += "Ellipse(%s) = {%s, %s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getPointOnMainAxis().getDirectedID(), p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)

            elif isinstance(p, CurveLoop):
                out += "Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))

            elif isinstance(p, RuledSurface):
                out += "Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)

            elif isinstance(p, PlaneSurface):
                line = self.__mkArgs(p.getHoles())
                if len(line) > 0:
                    out += "Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
                else:
                    out += "Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)

            elif isinstance(p, SurfaceLoop):
                out += "Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))

            elif isinstance(p, Volume):
                line = self.__mkArgs(p.getHoles())
                if len(line)>0:
                    out += "Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)+self.__mkTransfiniteVolume(p)
                else:
                    out += "Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())+self.__mkTransfiniteVolume(p)

            elif isinstance(p, PropertySet):
                if p.getNumItems() > 0:
                    dim=p.getDim()
                    line = "Physical "
                    if dim==0:
                        line += "Point"
                    elif dim==1:
                        line += "Line"
                    elif dim==2:
                        line += "Surface"
                    else:
                        line += "Volume"
                    out += line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"

            else:
                raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
        return out

    def __mkArgs(self, args, useAbs=False):
        line = ""
        for i in args:
            id = i.getDirectedID()
            if useAbs: id=abs(id)
            if len(line) > 0:
                line += ", %s"%id
            else:
                line = "%s"%id
        return line

    def __mkTransfiniteLine(self, p):
        s = p.getElementDistribution()
        if not s == None:
            if s[2]:
                out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
            else:
                out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
        else:
            out=""
        return out

    def __mkTransfiniteSurface(self, p):
        out = ""
        o = p.getRecombination()
        s = p.getTransfiniteMeshing()
        if not s == None:
            out2 = ""
            if not s[0] is None:
                for q in s[0]:
                    if len(out2)==0:
                        out2 = "%s"%q.getID()
                    else:
                        out2 = "%s,%s"%(out2, q.getID())
            if s[1] is None:
                out += "Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
            else:
                out += "Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
        if not o is None:
            out += "Recombine Surface {%s} = %f;\n"%(p.getID(), o/DEG)
        return out

    def __mkTransfiniteVolume(self, p):
        out=""
        s=p.getTransfiniteMeshing()
        if not s == None:
            if len(s)>0:
                out2=""
                for q in s[0]:
                    if len(out2)==0:
                        out2="%s"%q.getID()
                    else:
                        out2="%s,%s"%(out2,q.getID())
                out+="Transfinite Volume{%s} = {%s};\n"%(p.getID(),out2)
            else:
                out+="Transfinite Volume{%s};\n"%(p.getID(),)
        return out

1	gross	2935	# -- coding: utf-8 --
2	ksteube	1809
3	jfenwick	3981	##############################################################################
4	ksteube	1312	#
5	uqaeller	6939	# Copyright (c) 2003-2020 by The University of Queensland
6	jfenwick	3981	# http://www.uq.edu.au
7	ksteube	1312	#
8	ksteube	1809	# Primary Business: Queensland, Australia
9	jfenwick	6112	# Licensed under the Apache License, version 2.0
10			# http://www.apache.org/licenses/LICENSE-2.0
11	ksteube	1312	#
12	jfenwick	3981	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
14			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
15	uqaeller	6939	# Development from 2019 by School of Earth and Environmental Sciences
16	jfenwick	3981	#
17			##############################################################################
18	gross	932
19	sshaw	5706	from __future__ import print_function, division
20
21	uqaeller	6939	__copyright__="""Copyright (c) 2003-2020 by The University of Queensland
22	jfenwick	3981	http://www.uq.edu.au
23	ksteube	1809	Primary Business: Queensland, Australia"""
24	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
25			http://www.apache.org/licenses/LICENSE-2.0"""
26	jfenwick	2344	__url__="https://launchpad.net/escript-finley"
27	ksteube	1809
28	gross	932	"""
29			mesh generation using gmsh
30
31	jfenwick	2625	:var __author__: name of author
32			:var __copyright__: copyrights
33			:var __license__: licence agreement
34			:var __url__: url entry point on documentation
35			:var __version__: version
36			:var __date__: date of the version
37	gross	932	"""
38
39			__author__="Lutz Gross, l.gross@uq.edu.au"
40
41	jfenwick	3774	from . import design
42	gross	932	import tempfile
43			import os
44	jfenwick	3774	from .primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
45	caltinay	5397	from esys.escript import getMPIWorldMax, getMPIRankWorld, gmshGeo2Msh
46	jfenwick	3774	from .transformations import DEG
47	gross	932
48	sshaw	4559	class Design(design.AbstractDesign):
49	gross	932	"""
50	caltinay	2180	Design for gmsh.
51	gross	932	"""
52	gross	3002	DELAUNAY="Delauny"
53			MESHADAPT="MeshAdapt"
54			FRONTAL="Frontal"
55			NETGEN="Frontal"
56			TETGEN="Delauny"
57	phornby	2096
58	caltinay	5392	def __init__(self, dim=3, element_size=1., order=1, keep_files=False):
59			"""
60			Initializes the gmsh design.
61	gross	3002
62	caltinay	5392	:param dim: spatial dimension
63			:param element_size: global element size
64			:param order: element order
65			:param keep_files: flag to keep work files
66			"""
67			design.AbstractDesign.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
68			self.__mshname_set = False
69			self.__scriptname=""
70			self.setScriptFileName()
71			self.setOptions()
72			self.setFileFormat(self.GMSH)
73	gross	932
74	caltinay	5392	def setScriptFileName(self, name=None):
75			"""
76			Sets the filename for the gmsh input script. If no name is given a name
77			with extension `geo` is generated.
78			"""
79			if self.__scriptname:
80			os.unlink(self.__scriptname)
81			if name == None:
82			self.__scriptname_set=False
83			tmp_f_id=tempfile.mkstemp(suffix=".geo")
84			self.__scriptname=tmp_f_id[1]
85			os.close(tmp_f_id[0])
86			else:
87			self.__scriptname=name
88			self.__scriptname_set=True
89	caltinay	2180
90	gross	932	def getScriptFileName(self):
91	caltinay	5392	"""
92			Returns the name of the gmsh script file.
93			"""
94			return self.__scriptname
95	caltinay	2180
96	caltinay	5392	def setOptions(self, algorithm=None, optimize_quality=True, smoothing=1,
97			curvature_based_element_size=False, algorithm2D=None,
98			algorithm3D=None, generate_hexahedra=False,
99			random_factor=None):
100	gross	932	"""
101	caltinay	2180	Sets options for the mesh generator.
102	caltinay	5392
103			:param algorithm: selects 2D meshing algorithm
104			:type algorithm: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
105	sshaw	4559	:param algorithm2D: must be equal to algorithm
106	caltinay	5392	:type algorithm2D: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
107			:param algorithm3D: selects 3D meshing algorithm
108			:type algorithm3D: in self.DELAUNAY, self.FRONTAL
109	gross	3002	:param curvature_based_element_size: switch for curvature based element size adaption
110			:type curvature_based_element_size: ```bool```
111			:param smoothing: number of smoothing steps
112			:type smoothing: non-negative ```int```
113			:param optimize_quality: switch for mesh quality optimization
114			:type optimize_quality: ```bool```
115			:param generate_hexahedra: switch for using quadrangles/hexahedra elements everywhere.
116			:type generate_hexahedra: ```bool```
117	gross	4028	:param random_factor: used in the 2D meshing algorithm (should be increased if RandomFactor * size(triangle)/size(model) approaches machine accuracy)
118			:type random_factor: positive ```float```
119	caltinay	2180	"""
120	gross	4028	if random_factor==None: random_factor=1.e-9
121	caltinay	5392	if not random_factor > 0:
122			raise ValueError("random_factor must be positive.")
123	gross	4028	smoothing=int(smoothing)
124	caltinay	5392	if not smoothing > 0:
125			raise ValueError("smoothing must be positive.")
126	gross	4028
127	caltinay	5392	if algorithm3D is None:
128			algorithm3D=self.FRONTAL
129			if algorithm is None:
130			if algorithm2D is None:
131			algorithm2D=self.MESHADAPT
132	gross	3002	else:
133	caltinay	5392	if not algorithm2D is None:
134			if not algorithm == algorithm2D:
135			raise ValueError("argument algorithm (=%s) and algorithm2D (=%s) must have the same value if set."%(algorithm, algorithm2D))
136			algorithm2D = algorithm
137	gross	3002	if not algorithm2D in [ self.DELAUNAY, self.MESHADAPT, self.FRONTAL ]:
138	jfenwick	3774	raise ValueError("illegal 2D meshing algorithm %s."%algorithm2D)
139	gross	3002	if not algorithm3D in [ self.DELAUNAY, self.FRONTAL ]:
140	jfenwick	3774	raise ValueError("illegal 3D meshing algorithm %s."%algorithm3D)
141	caltinay	5392
142	gross	3002	self.__curvature_based_element_size=curvature_based_element_size
143			self.__algo2D=algorithm2D
144			self.__algo3D=algorithm3D
145	gross	932	self.__optimize_quality=optimize_quality
146			self.__smoothing=smoothing
147	gross	3002	self.__generate_hexahedra=generate_hexahedra
148	gross	4028	self.__random_factor=random_factor
149	caltinay	5392
150			def getOptions(self, name=None):
151	gross	3002	"""
152			Returns the current options for the mesh generator.
153			"""
154	caltinay	5392	if name is None:
155			return {"optimize_quality" : self.__optimize_quality ,
156			"smoothing" : self.__smoothing,
157			"curvature_based_element_size" : self.__curvature_based_element_size,
158			"generate_hexahedra" : self.__generate_hexahedra,
159			"algorithm2D" : self.__algo2D,
160			"algorithm3D" : self.__algo3D ,
161			"random_factor" : self.__random_factor }
162	gross	3002	else:
163	caltinay	5392	return self.getOption()[name]
164	caltinay	2180
165	gross	932	def __del__(self):
166			"""
167	caltinay	2180	Cleans up.
168	gross	932	"""
169	jfenwick	5330	try:
170	caltinay	5392	if not self.keepFiles():
171			if not self.__scriptname_set: #i.e. it's a tempfile
172			os.unlink(self.getScriptFileName())
173			if not self.__mshname_set: #i.e. it's a tempfile
174			os.unlink(self.getMeshFileName())
175	jfenwick	5330	except OSError:
176	sshaw	5348	pass # The file might not have been created and there is nothing
177			# to do about a "failure" here anyway
178	phornby	2096
179	gross	2620	def getScriptHandler(self):
180			"""
181			Returns a handler to the script file to generate the geometry.
182			In the current implementation a script file name is returned.
183			"""
184			if getMPIRankWorld() == 0:
185			open(self.getScriptFileName(),"w").write(self.getScriptString())
186			return self.getScriptFileName()
187	caltinay	2180
188	gross	932	def getMeshHandler(self):
189			"""
190	caltinay	2180	Returns a handle to a mesh meshing the design. In the current
191			implementation a mesh file name in gmsh format is returned.
192	gross	932	"""
193	caltinay	5392
194	caltinay	5397	verbosity = 3
195	caltinay	5398	ret = gmshGeo2Msh(self.getScriptHandler(), self.getMeshFileName(),
196	caltinay	5397	self.getDim(), self.getElementOrder(), verbosity)
197			if ret > 0:
198			self.setKeepFilesOn() #no files to delete, so don't try to
199			raise RuntimeError("Could not build mesh using gmsh.\n" + \
200			"Is gmsh available?")
201			return self.getMeshFileName()
202	gross	1045
203			def getScriptString(self):
204			"""
205	caltinay	2180	Returns the gmsh script to generate the mesh.
206	gross	1045	"""
207			h=self.getElementSize()
208	gross	3002	out='// generated by esys.pycad\nGeneral.Terminal = 1;\nGeneral.ExpertMode = 1;\n'
209			options=self.getOptions()
210			if options["optimize_quality"]:
211	caltinay	5392	out += "Mesh.Optimize = 1;\n"
212	gross	3002	else:
213	caltinay	5392	out += "Mesh.Optimize = 0;\n"
214
215	gross	3002	if options["curvature_based_element_size"]:
216	caltinay	5392	out += "Mesh.CharacteristicLengthFromCurvature = 1;\n"
217	gross	3002	else:
218	caltinay	5392	out += "Mesh.CharacteristicLengthFromCurvature = 0;\n"
219
220	gross	3002	if options["generate_hexahedra"]:
221	caltinay	5392	if self.getDim() == 2:
222			out += "Mesh.SubdivisionAlgorithm = 1;\n"
223			else:
224			out += "Mesh.SubdivisionAlgorithm = 2;\n"
225	gross	3002	else:
226	caltinay	5392	out += "Mesh.SubdivisionAlgorithm = 0;\n"
227
228			out += "Mesh.Smoothing = %d;\n"%options["smoothing"]
229			out += "Mesh.RandomFactor = %.14e;\n"%options["random_factor"]
230			if options["algorithm2D"] == self.MESHADAPT:
231			out += "Mesh.Algorithm = 1; // = MeshAdapt\n"
232			elif options["algorithm2D"] == self.DELAUNAY:
233			out += "Mesh.Algorithm = 5; // = Delaunay\n"
234			elif options["algorithm2D"] == self.FRONTAL:
235			out += "Mesh.Algorithm = 6; // = Frontal\n"
236
237			if options["algorithm3D"] == self.DELAUNAY:
238			out += "Mesh.Algorithm3D = 1; // = Delaunay\n"
239			elif options["algorithm3D"] == self.FRONTAL:
240			out += "Mesh.Algorithm3D = 4; // = Frontal\n"
241
242	gross	1045	for prim in self.getAllPrimitives():
243	caltinay	5392	p=prim.getUnderlyingPrimitive()
244			if isinstance(p, Point):
245			c=p.getCoordinates()
246			#out+="Point(%s) = {%f , %f, %f , %f };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
247			out += "Point(%s) = {%.14e, %.14e, %.14e, %.14e};\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
248	caltinay	2180
249	caltinay	5392	elif isinstance(p, Spline):
250			out += "Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
251	caltinay	2180
252	caltinay	5392	elif isinstance(p, BezierCurve):
253			out += "Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
254	gross	1045
255	caltinay	5392	elif isinstance(p, BSpline):
256			out += "BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
257	gross	1045
258	caltinay	5392	elif isinstance(p, Line):
259			out += "Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
260	gross	1045
261	caltinay	5392	elif isinstance(p, Arc):
262			out += "Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
263	caltinay	2180
264	caltinay	5392	elif isinstance(p, Ellipse):
265			out += "Ellipse(%s) = {%s, %s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getPointOnMainAxis().getDirectedID(), p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
266	ksteube	1312
267	caltinay	5392	elif isinstance(p, CurveLoop):
268			out += "Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
269	caltinay	2180
270	caltinay	5392	elif isinstance(p, RuledSurface):
271			out += "Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
272	caltinay	2180
273	caltinay	5392	elif isinstance(p, PlaneSurface):
274			line = self.__mkArgs(p.getHoles())
275			if len(line) > 0:
276			out += "Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
277			else:
278			out += "Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
279	caltinay	2180
280	caltinay	5392	elif isinstance(p, SurfaceLoop):
281			out += "Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
282	caltinay	2180
283	caltinay	5392	elif isinstance(p, Volume):
284			line = self.__mkArgs(p.getHoles())
285			if len(line)>0:
286			out += "Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)+self.__mkTransfiniteVolume(p)
287			else:
288			out += "Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())+self.__mkTransfiniteVolume(p)
289	gross	1045
290	caltinay	5392	elif isinstance(p, PropertySet):
291			if p.getNumItems() > 0:
292			dim=p.getDim()
293			line = "Physical "
294			if dim==0:
295			line += "Point"
296			elif dim==1:
297			line += "Line"
298			elif dim==2:
299			line += "Surface"
300			else:
301			line += "Volume"
302			out += line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
303	gross	1045
304	caltinay	5392	else:
305			raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
306	gross	1045	return out
307
308	caltinay	5392	def __mkArgs(self, args, useAbs=False):
309			line = ""
310	gross	1045	for i in args:
311	gross	2700	id = i.getDirectedID()
312			if useAbs: id=abs(id)
313	caltinay	5392	if len(line) > 0:
314			line += ", %s"%id
315	gross	1045	else:
316	caltinay	5392	line = "%s"%id
317	caltinay	2180	return line
318
319	caltinay	5392	def __mkTransfiniteLine(self, p):
320			s = p.getElementDistribution()
321			if not s == None:
322			if s[2]:
323			out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
324			else:
325			out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
326			else:
327			out=""
328			return out
329
330			def __mkTransfiniteSurface(self, p):
331			out = ""
332			o = p.getRecombination()
333			s = p.getTransfiniteMeshing()
334			if not s == None:
335			out2 = ""
336			if not s[0] is None:
337			for q in s[0]:
338			if len(out2)==0:
339			out2 = "%s"%q.getID()
340			else:
341			out2 = "%s,%s"%(out2, q.getID())
342			if s[1] is None:
343			out += "Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
344			else:
345			out += "Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
346			if not o is None:
347			out += "Recombine Surface {%s} = %f;\n"%(p.getID(), o/DEG)
348			return out
349
350			def __mkTransfiniteVolume(self, p):
351			out=""
352			s=p.getTransfiniteMeshing()
353			if not s == None:
354			if len(s)>0:
355			out2=""
356			for q in s[0]:
357			if len(out2)==0:
358			out2="%s"%q.getID()
359			else:
360			out2="%s,%s"%(out2,q.getID())
361			out+="Transfinite Volume{%s} = {%s};\n"%(p.getID(),out2)
362			else:
363			out+="Transfinite Volume{%s};\n"%(p.getID(),)
364			return out
365