added structured abtuse meshing for cuboid

_setup/geometries/geometry.py

@@ -7,13 +7,14 @@
 class Geometry:
   
   def __init__(self, geometry=None, meshSize=None, scaling=1.0 \
-    , structuredMesh=False):
+    , structuredMesh=False, structuredAbtuseMesh=False):
     
     self._geometry = geometry
     self.meshSize = meshSize
     self.scaling = scaling
 
     self.structuredMesh = structuredMesh
+    self.structuredAbtuseMesh = structuredAbtuseMesh
     
 #TOECDO    self.boundingBox = None
 
@@ -53,7 +54,7 @@ class Geometry:
       ngmesh = netgenMesh(3)
       ngmesh.Load( self._geometry )
     else:
-      if(self.structuredMesh):
+      if(self.structuredMesh or self.structuredAbtuseMesh):
         ngmesh = self._GenerateStructuredMesh()
       else:
         if( isinstance(self._geometry, _standardGeometry._StandardGeometry)):
@@ -80,15 +81,22 @@ class Geometry:
   def _GenerateStructuredMesh(self):
     from commix._setup.geometries.standardGeometries import cuboid
     from commix._setup.geometries.standardGeometries._structuredMeshing \
-      import StructuredMesh_Cuboid
+      import StructuredMesh_Cuboid, StructuredAbtuseMesh_Cuboid
     if(len(self.meshSize) != 3):
       raise ValueError("For structured mesh, provide meshSize as (hx, hy, hz)")
     if not(isinstance(self._geometry, cuboid.Cuboid)):
       raise NotImplementedError("Structured meshing only for Cuboid geo")
 
+    if self.structuredAbtuseMesh and self.structuredMesh:
+      raise ValueError("structuredAbtuseMesh=True and structuredMesh=True ?!")
+
     p1, p2 = self._geometry.GetXminXmax()
 
-    return StructuredMesh_Cuboid(p1, p2, self.meshSize)
+    if self.structuredMesh:
+      return StructuredMesh_Cuboid(p1, p2, self.meshSize)
+    else:
+      return StructuredAbtuseMesh_Cuboid(p1, p2, self.meshSize)
+
 
 
 #------------------------------------------------------------------------------#

_setup/geometries/standardGeometries/_structuredMeshing/__init__.py

@@ -2,5 +2,7 @@ print("Importing module", __name__)
 
 from commix._setup.geometries.standardGeometries._structuredMeshing.structuredMesh_Cuboid \
   import StructuredMesh_Cuboid
+from commix._setup.geometries.standardGeometries._structuredMeshing.structuredAbtuseMesh_Cuboid \
+  import StructuredAbtuseMesh_Cuboid
 
-__all__ = ["StructuredMesh_Cuboid"]
+__all__ = ["StructuredMesh_Cuboid", "StructuredAbtuseMesh_Cuboid"]

_setup/geometries/standardGeometries/_structuredMeshing/addStructuredAbtuseMesh_MicroCuboid.py

+
+
+################################################################################
+
+
+def AddStructuredAbtuseMesh_MicroCuboid(ngmesh, p, center):
+  from netgen.meshing import Element3D
+
+  ngmesh.Add(Element3D(1, [p[0, 0, 0], p[1, 0, 0], center, p[1, 1, 0]]))
+  ngmesh.Add(Element3D(1, [p[1, 0, 0], p[1, 0, 1], center, p[1, 1, 0]]))
+  ngmesh.Add(Element3D(1, [p[1, 0, 1], p[1, 1, 1], center, p[1, 1, 0]]))
+  ngmesh.Add(Element3D(1, [p[1, 1, 1], p[0, 1, 1], center, p[1, 1, 0]]))
+  ngmesh.Add(Element3D(1, [p[0, 1, 1], p[0, 1, 0], center, p[1, 1, 0]]))
+  ngmesh.Add(Element3D(1, [p[0, 1, 0], p[0, 0, 0], center, p[1, 1, 0]]))
+
+  ngmesh.Add(Element3D(1, [p[0, 0, 0], p[1, 0, 0], p[0, 0, 1], center]))
+  ngmesh.Add(Element3D(1, [p[1, 0, 0], p[1, 0, 1], p[0, 0, 1], center]))
+  ngmesh.Add(Element3D(1, [p[1, 0, 1], p[1, 1, 1], p[0, 0, 1], center]))
+  ngmesh.Add(Element3D(1, [p[1, 1, 1], p[0, 1, 1], p[0, 0, 1], center]))
+  ngmesh.Add(Element3D(1, [p[0, 1, 1], p[0, 1, 0], p[0, 0, 1], center]))
+  ngmesh.Add(Element3D(1, [p[0, 1, 0], p[0, 0, 0], p[0, 0, 1], center]))
+
+################################################################################

_setup/geometries/standardGeometries/_structuredMeshing/structuredAbtuseMesh_Cuboid.py

+
+
+################################################################################
+
+
+def StructuredAbtuseMesh_Cuboid(Xmin, Xmax, h):
+  from commix._setup.geometries.standardGeometries._structuredMeshing.addStructuredAbtuseMesh_MicroCuboid \
+    import AddStructuredAbtuseMesh_MicroCuboid
+  from commix._setup.geometries.standardGeometries._structuredMeshing.addStructuredMesh_BoundaryFaceCuboid \
+    import AddStructuredMesh_BoundaryFaceCuboid
+
+  from netgen.meshing import MeshPoint, Mesh, Pnt, FaceDescriptor
+  import numpy as np
+  
+  Ncells = [0, 0, 0]
+  Nnodes = [0, 0, 0]
+  for d in range(3):
+    Ncells[d] = round((Xmax[d]-Xmin[d])/h[d])
+    Nnodes[d] = Ncells[d]+1
+  
+  ngmesh = Mesh()
+  ngmesh.dim = 3
+  
+  nodes = np.empty(Nnodes, dtype=object)
+  inside_nodes = np.empty(Ncells, dtype=object)
+  for nx in range(Nnodes[0]):
+    for ny in range(Nnodes[1]):
+      for nz in range(Nnodes[2]):
+        nodes[nx, ny, nz] = ngmesh.Add(MeshPoint( \
+          Pnt(Xmin[0]+nx*h[0], Xmin[1]+ny*h[1], Xmin[2]+nz*h[2])))
+
+  for nx in range(Ncells[0]):
+    for ny in range(Ncells[1]):
+      for nz in range(Ncells[2]):
+        inside_nodes[nx, ny, nz] = ngmesh.Add(MeshPoint( \
+          Pnt(Xmin[0]+(nx+0.5)*h[0], Xmin[1]+(ny+0.5)*h[1], Xmin[2]+(nz+0.5)*h[2])))
+  
+  
+  for cx in range(Ncells[0]):
+    for cy in range(Ncells[1]):
+      for cz in range(Ncells[2]):
+        AddStructuredAbtuseMesh_MicroCuboid(ngmesh, nodes[cx:cx+2, cy:cy+2, cz:cz+2], inside_nodes[cx, cy, cz])
+  
+  fd = ngmesh.Add (FaceDescriptor(bc=1,domin=1,domout=0,surfnr=1))
+  
+  AddStructuredMesh_BoundaryFaceCuboid(ngmesh, fd, nodes[0, :, :], "xmin")
+  AddStructuredMesh_BoundaryFaceCuboid(ngmesh, fd, nodes[-1, :, :], "xmax")
+  AddStructuredMesh_BoundaryFaceCuboid(ngmesh, fd, nodes[:, 0, :], "ymin")
+  AddStructuredMesh_BoundaryFaceCuboid(ngmesh, fd, nodes[:, -1, :], "ymax")
+  AddStructuredMesh_BoundaryFaceCuboid(ngmesh, fd, nodes[:, :, 0], "zmin")
+  AddStructuredMesh_BoundaryFaceCuboid(ngmesh, fd, nodes[:, :, -1], "zmax")
+  
+  # for unknown reasons, it seems this number (0) has to be the bc parameter in 
+  # FaceDescriptor, but -1.     ...
+  ngmesh.SetBCName(0, "bc_dirichlet")
+
+  return ngmesh
+
+
+################################################################################

tests/abtuseStructuredExperiment.py

+from commix import *
+
+h = 10 / 2
+
+geo = Geometry(geometry=Cuboid(120, 120, 10), meshSize=(h, h, h), scaling=1.0e-09, structuredAbtuseMesh=True)
+
+from math import pi
+par = Parameters(A=pi*2.0e-11, Ms=1.0e6, alpha=0.1, K=0.0 \
+  , D=pi*2.0e-3, dmCoupling="bulk", timeStepSize=0.1e-12, theta=1.0)
+
+par.maxwellCoupling = "none"
+
+par.m0 = (0, 0, 1)
+
+skyrmion = Integrator(scheme=TPS1, geometry=geo, parameters=par, numthreads=4)
+skyrmion.SetResultsFolder("ACOM_AbtuseStructured_Result")
+
+skyrmion.RecordMagnetization(record=True, maxwell=False, every=10e-12)
+
+skyrmion.Integrate(duration=1.36e-09)