Merge branch 'ngMeshingTools'

_setup/geometries/geometry.py

@@ -34,6 +34,8 @@ class Geometry:
      meshSize [float]
      scaling [float], e.g., 1e-9 (nm)
      structuedMesh=False
+  in B+D) force the maximum meshSize requirement by
+      forceNetgenMeshSize=True
 
   NOTE:
     The scaling parameter seems artificial. Indeed it is.
@@ -44,11 +46,12 @@ class Geometry:
     than Cuboid((0,0,0), (1e-9, 1e-9, 1e-9)) + scaling=1
   '''
   
-  def __init__(self, geometry=None, meshSize=None, scaling=1.0 \
+  def __init__(self, geometry=None, meshSize=None, forceNetgenMeshSize=False, scaling=1.0 \
     , structuredMesh=False, structuredAbtuseMesh=False):
     
     self._geometry = geometry
     self.meshSize = meshSize
+    self.forceNetgenMeshSize = forceNetgenMeshSize
     self.scaling = scaling
 
     self.structuredMesh = structuredMesh
@@ -68,6 +71,9 @@ class Geometry:
     from ngsolve import Mesh as ngsolveMesh
     from netgen.meshing import Mesh as netgenMesh
     from commics._setup.geometries.standardGeometries import _standardGeometry
+    from ngsolve import ngsglobals
+
+    ngsglobals.msg_level = 0
 
     if not(self._mesh == None):
       print("From Geometry.GetReady(): Nothing to do, mesh is ready.")
@@ -108,7 +114,20 @@ class Geometry:
         geo = CSGeometry()
         geo.Add(ngGeo)
         
-        ngmesh = geo.GenerateMesh( maxh = self.meshSize )
+        if self.forceNetgenMeshSize:
+          from commics._tools.misc import MaximumMeshSize
+          achievedMeshSize = False
+          j = 1
+          print("Netgen meshing ........ Desired meshSize =", self.meshSize) 
+          while not(achievedMeshSize):
+            ngmesh = geo.GenerateMesh( maxh = self.meshSize / j )
+            actualMeshSize = MaximumMeshSize(ngmesh)
+            achievedMeshSize = actualMeshSize <= self.meshSize
+            print("Itaration Nr.%i: Desired meshSize = %g --- Actual meshSize = %g" % (j, self.meshSize/j, actualMeshSize))
+            j += 1
+          print("Appropriate mesh was generated ........") 
+        else:
+          ngmesh = geo.GenerateMesh( maxh = self.meshSize )
       
     self._mesh = ngsolveMesh( ngmesh )
   
@@ -184,6 +203,7 @@ class Geometry:
       string += "Vertices: " + str(self._mesh.nv) + "\n"
       string += "Elements: " + str(self._mesh.ne) + "\n"
       string += "meshSize = " + str(self.meshSize) + "\n"
+      string += "forceNetgenMeshSize = " + str(self.forceNetgenMeshSize) + "\n"
       string += "scaling = " + str(self.scaling) +"\n"
       string += "structuredMesh = " + str(self.structuredMesh) +"\n"
 
@@ -216,7 +236,11 @@ class Geometry:
   
   def Save(self, filename):
     if not(self._mesh == None):
-      self._mesh.ngmesh.Save( filename )
+      ext = ".vol"
+      if(filename[-4:] == ext):
+        ext = ""
+      self._mesh.ngmesh.Save( filename + ext )
+
 # TOECDO
 #      if ( self.boundingBox.IsBoundingBox() ):
 #        tmp = filename.rfind("/")+1

_tools/misc/__init__.py

@@ -4,6 +4,7 @@ if __name__ == "__main__":
   print("Hello, this is the module. Explanation...")
 
 from commics._tools.misc.timeToString import TimeToString
+from commics._tools.misc.maximumMeshSize import MaximumMeshSize
 
 
-__all__ = ['TimeToString']
+__all__ = ['TimeToString', 'MaximumMeshSize']

_tools/misc/maximumMeshSize.py

+# return h_max of a ngmesh, object of class netgen.libngpy._meshing.Mesh
+
+
+#------------------------------------------------------------------------------#
+
+
+def MaximumMeshSize(ngmesh):
+  from math import sqrt
+
+  hmax = 0.0
+  for e in range(len(ngmesh.Elements3D())):
+    el = ngmesh.Elements3D()[e]
+    for i in range(4):
+      for j in range(i+1,4):
+        p1 = ngmesh.Points()[el.vertices[i]].p
+        p2 = ngmesh.Points()[el.vertices[j]].p
+  
+        h = sqrt(sum((p1[k]-p2[k])**2 for k in range(3)))
+
+        if h > hmax:
+          hmax = h
+
+  return hmax
+
+
+#------------------------------------------------------------------------------#