diff --git a/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py b/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py
index 47893db..9bce4d7 100644
--- a/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py
+++ b/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py
@@ -1,362 +1,364 @@
import numpy as np
import h5py
import GooseFEM as gf
# ====================== create fictitious configuration + store to HDF5-file ======================
# create mesh object
mesh = gf.Mesh.Hex8.FineLayer(9,17,27)
# set node identifier for all boundary sets
Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int')
nodesFront = Z(mesh); nodesFront [mesh.nodesFront() ] = 1
nodesBack = Z(mesh); nodesBack [mesh.nodesBack() ] = 1
nodesLeft = Z(mesh); nodesLeft [mesh.nodesLeft() ] = 1
nodesRight = Z(mesh); nodesRight [mesh.nodesRight() ] = 1
nodesBottom = Z(mesh); nodesBottom [mesh.nodesBottom() ] = 1
nodesTop = Z(mesh); nodesTop [mesh.nodesTop() ] = 1
nodesFrontFace = Z(mesh); nodesFrontFace [mesh.nodesFrontFace() ] = 1
nodesBackFace = Z(mesh); nodesBackFace [mesh.nodesBackFace() ] = 1
nodesLeftFace = Z(mesh); nodesLeftFace [mesh.nodesLeftFace() ] = 1
nodesRightFace = Z(mesh); nodesRightFace [mesh.nodesRightFace() ] = 1
nodesBottomFace = Z(mesh); nodesBottomFace [mesh.nodesBottomFace() ] = 1
nodesTopFace = Z(mesh); nodesTopFace [mesh.nodesTopFace() ] = 1
nodesFrontBottomEdge = Z(mesh); nodesFrontBottomEdge [mesh.nodesFrontBottomEdge() ] = 1
nodesFrontTopEdge = Z(mesh); nodesFrontTopEdge [mesh.nodesFrontTopEdge() ] = 1
nodesFrontLeftEdge = Z(mesh); nodesFrontLeftEdge [mesh.nodesFrontLeftEdge() ] = 1
nodesFrontRightEdge = Z(mesh); nodesFrontRightEdge [mesh.nodesFrontRightEdge() ] = 1
nodesBackBottomEdge = Z(mesh); nodesBackBottomEdge [mesh.nodesBackBottomEdge() ] = 1
nodesBackTopEdge = Z(mesh); nodesBackTopEdge [mesh.nodesBackTopEdge() ] = 1
nodesBackLeftEdge = Z(mesh); nodesBackLeftEdge [mesh.nodesBackLeftEdge() ] = 1
nodesBackRightEdge = Z(mesh); nodesBackRightEdge [mesh.nodesBackRightEdge() ] = 1
nodesBottomLeftEdge = Z(mesh); nodesBottomLeftEdge [mesh.nodesBottomLeftEdge() ] = 1
nodesBottomRightEdge = Z(mesh); nodesBottomRightEdge [mesh.nodesBottomRightEdge() ] = 1
nodesTopLeftEdge = Z(mesh); nodesTopLeftEdge [mesh.nodesTopLeftEdge() ] = 1
nodesTopRightEdge = Z(mesh); nodesTopRightEdge [mesh.nodesTopRightEdge() ] = 1
nodesFrontBottomOpenEdge = Z(mesh); nodesFrontBottomOpenEdge [mesh.nodesFrontBottomOpenEdge() ] = 1
nodesFrontTopOpenEdge = Z(mesh); nodesFrontTopOpenEdge [mesh.nodesFrontTopOpenEdge() ] = 1
nodesFrontLeftOpenEdge = Z(mesh); nodesFrontLeftOpenEdge [mesh.nodesFrontLeftOpenEdge() ] = 1
nodesFrontRightOpenEdge = Z(mesh); nodesFrontRightOpenEdge [mesh.nodesFrontRightOpenEdge() ] = 1
nodesBackBottomOpenEdge = Z(mesh); nodesBackBottomOpenEdge [mesh.nodesBackBottomOpenEdge() ] = 1
nodesBackTopOpenEdge = Z(mesh); nodesBackTopOpenEdge [mesh.nodesBackTopOpenEdge() ] = 1
nodesBackLeftOpenEdge = Z(mesh); nodesBackLeftOpenEdge [mesh.nodesBackLeftOpenEdge() ] = 1
nodesBackRightOpenEdge = Z(mesh); nodesBackRightOpenEdge [mesh.nodesBackRightOpenEdge() ] = 1
nodesBottomLeftOpenEdge = Z(mesh); nodesBottomLeftOpenEdge [mesh.nodesBottomLeftOpenEdge() ] = 1
nodesBottomRightOpenEdge = Z(mesh); nodesBottomRightOpenEdge [mesh.nodesBottomRightOpenEdge() ] = 1
nodesTopLeftOpenEdge = Z(mesh); nodesTopLeftOpenEdge [mesh.nodesTopLeftOpenEdge() ] = 1
nodesTopRightOpenEdge = Z(mesh); nodesTopRightOpenEdge [mesh.nodesTopRightOpenEdge() ] = 1
nodesFrontBottomLeftCorner = Z(mesh); nodesFrontBottomLeftCorner [mesh.nodesFrontBottomLeftCorner() ] = 1
nodesFrontBottomRightCorner = Z(mesh); nodesFrontBottomRightCorner[mesh.nodesFrontBottomRightCorner()] = 1
nodesFrontTopLeftCorner = Z(mesh); nodesFrontTopLeftCorner [mesh.nodesFrontTopLeftCorner() ] = 1
nodesFrontTopRightCorner = Z(mesh); nodesFrontTopRightCorner [mesh.nodesFrontTopRightCorner() ] = 1
nodesBackBottomLeftCorner = Z(mesh); nodesBackBottomLeftCorner [mesh.nodesBackBottomLeftCorner() ] = 1
nodesBackBottomRightCorner = Z(mesh); nodesBackBottomRightCorner [mesh.nodesBackBottomRightCorner() ] = 1
nodesBackTopLeftCorner = Z(mesh); nodesBackTopLeftCorner [mesh.nodesBackTopLeftCorner() ] = 1
nodesBackTopRightCorner = Z(mesh); nodesBackTopRightCorner [mesh.nodesBackTopRightCorner() ] = 1
# DOFs after eliminating periodicity
dofsPeriodic = mesh.dofsPeriodic()[:,0]
# open data file
f = h5py.File('MeshHex8-FineLayer-nodes_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
+# write nodal positions, and a dummy connectivity
f.create_dataset('/coor' ,data=mesh.coor() )
f.create_dataset('/conn' ,data=np.arange(mesh.nnode()) )
f.create_dataset('/nodesFront' ,data=nodesFront )
f.create_dataset('/nodesBack' ,data=nodesBack )
f.create_dataset('/nodesLeft' ,data=nodesLeft )
f.create_dataset('/nodesRight' ,data=nodesRight )
f.create_dataset('/nodesBottom' ,data=nodesBottom )
f.create_dataset('/nodesTop' ,data=nodesTop )
f.create_dataset('/nodesFrontFace' ,data=nodesFrontFace )
f.create_dataset('/nodesBackFace' ,data=nodesBackFace )
f.create_dataset('/nodesLeftFace' ,data=nodesLeftFace )
f.create_dataset('/nodesRightFace' ,data=nodesRightFace )
f.create_dataset('/nodesBottomFace' ,data=nodesBottomFace )
f.create_dataset('/nodesTopFace' ,data=nodesTopFace )
f.create_dataset('/nodesFrontBottomEdge' ,data=nodesFrontBottomEdge )
f.create_dataset('/nodesFrontTopEdge' ,data=nodesFrontTopEdge )
f.create_dataset('/nodesFrontLeftEdge' ,data=nodesFrontLeftEdge )
f.create_dataset('/nodesFrontRightEdge' ,data=nodesFrontRightEdge )
f.create_dataset('/nodesBackBottomEdge' ,data=nodesBackBottomEdge )
f.create_dataset('/nodesBackTopEdge' ,data=nodesBackTopEdge )
f.create_dataset('/nodesBackLeftEdge' ,data=nodesBackLeftEdge )
f.create_dataset('/nodesBackRightEdge' ,data=nodesBackRightEdge )
f.create_dataset('/nodesBottomLeftEdge' ,data=nodesBottomLeftEdge )
f.create_dataset('/nodesBottomRightEdge' ,data=nodesBottomRightEdge )
f.create_dataset('/nodesTopLeftEdge' ,data=nodesTopLeftEdge )
f.create_dataset('/nodesTopRightEdge' ,data=nodesTopRightEdge )
f.create_dataset('/nodesFrontBottomOpenEdge' ,data=nodesFrontBottomOpenEdge )
f.create_dataset('/nodesFrontTopOpenEdge' ,data=nodesFrontTopOpenEdge )
f.create_dataset('/nodesFrontLeftOpenEdge' ,data=nodesFrontLeftOpenEdge )
f.create_dataset('/nodesFrontRightOpenEdge' ,data=nodesFrontRightOpenEdge )
f.create_dataset('/nodesBackBottomOpenEdge' ,data=nodesBackBottomOpenEdge )
f.create_dataset('/nodesBackTopOpenEdge' ,data=nodesBackTopOpenEdge )
f.create_dataset('/nodesBackLeftOpenEdge' ,data=nodesBackLeftOpenEdge )
f.create_dataset('/nodesBackRightOpenEdge' ,data=nodesBackRightOpenEdge )
f.create_dataset('/nodesBottomLeftOpenEdge' ,data=nodesBottomLeftOpenEdge )
f.create_dataset('/nodesBottomRightOpenEdge' ,data=nodesBottomRightOpenEdge )
f.create_dataset('/nodesTopLeftOpenEdge' ,data=nodesTopLeftOpenEdge )
f.create_dataset('/nodesTopRightOpenEdge' ,data=nodesTopRightOpenEdge )
f.create_dataset('/nodesFrontBottomLeftCorner' ,data=nodesFrontBottomLeftCorner )
f.create_dataset('/nodesFrontBottomRightCorner',data=nodesFrontBottomRightCorner)
f.create_dataset('/nodesFrontTopLeftCorner' ,data=nodesFrontTopLeftCorner )
f.create_dataset('/nodesFrontTopRightCorner' ,data=nodesFrontTopRightCorner )
f.create_dataset('/nodesBackBottomLeftCorner' ,data=nodesBackBottomLeftCorner )
f.create_dataset('/nodesBackBottomRightCorner' ,data=nodesBackBottomRightCorner )
f.create_dataset('/nodesBackTopLeftCorner' ,data=nodesBackTopLeftCorner )
f.create_dataset('/nodesBackTopRightCorner' ,data=nodesBackTopRightCorner )
f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic )
# ======================================== write XDMF-file =========================================
+# basic file-format
xmf = '''
MeshHex8-FineLayer-nodes_paraview.hdf5:/conn
MeshHex8-FineLayer-nodes_paraview.hdf5:/coor
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFront
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBack
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesLeft
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesRight
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottom
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTop
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontFace
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackFace
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesLeftFace
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesRightFace
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomFace
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopFace
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontLeftEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontRightEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackLeftEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackRightEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomLeftEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomRightEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopLeftEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopRightEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontLeftOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontRightOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackLeftOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackRightOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomLeftOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomRightOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopLeftOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopRightOpenEdge
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomLeftCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomRightCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopLeftCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopRightCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomLeftCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomRightCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopLeftCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopRightCorner
MeshHex8-FineLayer-nodes_paraview.hdf5:/dofsPeriodic
'''
+# write to file, fill mesh dimensions
open('MeshHex8-FineLayer-nodes_paraview.xmf','w').write(xmf.format(
nnode = mesh.nnode(),
))
diff --git a/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py b/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py
index c04f6a9..5198fbb 100644
--- a/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py
+++ b/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py
@@ -1,44 +1,56 @@
import numpy as np
import h5py
import GooseFEM as gf
# ====================== create fictitious configuration + store to HDF5-file ======================
# create mesh object
mesh = gf.Mesh.Hex8.FineLayer(9,17,27)
# open data file
f = h5py.File('MeshHex8-FineLayer_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
-f.create_dataset('/coor',data=mesh.coor())
-f.create_dataset('/conn',data=mesh.conn())
+# element set
+elementsMiddleLayer = np.zeros((mesh.nelem()),dtype='int')
+elementsMiddleLayer[mesh.elementsMiddleLayer()] = 1
+
+# write nodal coordinates and connectivity
+f.create_dataset('/coor' ,data=mesh.coor() )
+f.create_dataset('/conn' ,data=mesh.conn() )
+f.create_dataset('/elementsMiddleLayer',data=elementsMiddleLayer)
# ======================================== write XDMF-file =========================================
+# basic file-format
xmf = '''
MeshHex8-FineLayer_paraview.hdf5:/conn
MeshHex8-FineLayer_paraview.hdf5:/coor
+
+
+ MeshHex8-FineLayer_paraview.hdf5:/elementsMiddleLayer
+
+
'''
+# write to file, fill mesh dimensions
open('MeshHex8-FineLayer_paraview.xmf','w').write(xmf.format(
nnode = mesh.nnode(),
nelem = mesh.nelem(),
nne = mesh.nne(),
))
diff --git a/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py b/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py
index c03ef85..403e2f7 100644
--- a/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py
+++ b/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py
@@ -1,362 +1,364 @@
import numpy as np
import h5py
import GooseFEM as gf
# ====================== create fictitious configuration + store to HDF5-file ======================
# create mesh object
mesh = gf.Mesh.Hex8.Regular(6,8,12)
# set node identifier for all boundary sets
Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int')
nodesFront = Z(mesh); nodesFront [mesh.nodesFront() ] = 1
nodesBack = Z(mesh); nodesBack [mesh.nodesBack() ] = 1
nodesLeft = Z(mesh); nodesLeft [mesh.nodesLeft() ] = 1
nodesRight = Z(mesh); nodesRight [mesh.nodesRight() ] = 1
nodesBottom = Z(mesh); nodesBottom [mesh.nodesBottom() ] = 1
nodesTop = Z(mesh); nodesTop [mesh.nodesTop() ] = 1
nodesFrontFace = Z(mesh); nodesFrontFace [mesh.nodesFrontFace() ] = 1
nodesBackFace = Z(mesh); nodesBackFace [mesh.nodesBackFace() ] = 1
nodesLeftFace = Z(mesh); nodesLeftFace [mesh.nodesLeftFace() ] = 1
nodesRightFace = Z(mesh); nodesRightFace [mesh.nodesRightFace() ] = 1
nodesBottomFace = Z(mesh); nodesBottomFace [mesh.nodesBottomFace() ] = 1
nodesTopFace = Z(mesh); nodesTopFace [mesh.nodesTopFace() ] = 1
nodesFrontBottomEdge = Z(mesh); nodesFrontBottomEdge [mesh.nodesFrontBottomEdge() ] = 1
nodesFrontTopEdge = Z(mesh); nodesFrontTopEdge [mesh.nodesFrontTopEdge() ] = 1
nodesFrontLeftEdge = Z(mesh); nodesFrontLeftEdge [mesh.nodesFrontLeftEdge() ] = 1
nodesFrontRightEdge = Z(mesh); nodesFrontRightEdge [mesh.nodesFrontRightEdge() ] = 1
nodesBackBottomEdge = Z(mesh); nodesBackBottomEdge [mesh.nodesBackBottomEdge() ] = 1
nodesBackTopEdge = Z(mesh); nodesBackTopEdge [mesh.nodesBackTopEdge() ] = 1
nodesBackLeftEdge = Z(mesh); nodesBackLeftEdge [mesh.nodesBackLeftEdge() ] = 1
nodesBackRightEdge = Z(mesh); nodesBackRightEdge [mesh.nodesBackRightEdge() ] = 1
nodesBottomLeftEdge = Z(mesh); nodesBottomLeftEdge [mesh.nodesBottomLeftEdge() ] = 1
nodesBottomRightEdge = Z(mesh); nodesBottomRightEdge [mesh.nodesBottomRightEdge() ] = 1
nodesTopLeftEdge = Z(mesh); nodesTopLeftEdge [mesh.nodesTopLeftEdge() ] = 1
nodesTopRightEdge = Z(mesh); nodesTopRightEdge [mesh.nodesTopRightEdge() ] = 1
nodesFrontBottomOpenEdge = Z(mesh); nodesFrontBottomOpenEdge [mesh.nodesFrontBottomOpenEdge() ] = 1
nodesFrontTopOpenEdge = Z(mesh); nodesFrontTopOpenEdge [mesh.nodesFrontTopOpenEdge() ] = 1
nodesFrontLeftOpenEdge = Z(mesh); nodesFrontLeftOpenEdge [mesh.nodesFrontLeftOpenEdge() ] = 1
nodesFrontRightOpenEdge = Z(mesh); nodesFrontRightOpenEdge [mesh.nodesFrontRightOpenEdge() ] = 1
nodesBackBottomOpenEdge = Z(mesh); nodesBackBottomOpenEdge [mesh.nodesBackBottomOpenEdge() ] = 1
nodesBackTopOpenEdge = Z(mesh); nodesBackTopOpenEdge [mesh.nodesBackTopOpenEdge() ] = 1
nodesBackLeftOpenEdge = Z(mesh); nodesBackLeftOpenEdge [mesh.nodesBackLeftOpenEdge() ] = 1
nodesBackRightOpenEdge = Z(mesh); nodesBackRightOpenEdge [mesh.nodesBackRightOpenEdge() ] = 1
nodesBottomLeftOpenEdge = Z(mesh); nodesBottomLeftOpenEdge [mesh.nodesBottomLeftOpenEdge() ] = 1
nodesBottomRightOpenEdge = Z(mesh); nodesBottomRightOpenEdge [mesh.nodesBottomRightOpenEdge() ] = 1
nodesTopLeftOpenEdge = Z(mesh); nodesTopLeftOpenEdge [mesh.nodesTopLeftOpenEdge() ] = 1
nodesTopRightOpenEdge = Z(mesh); nodesTopRightOpenEdge [mesh.nodesTopRightOpenEdge() ] = 1
nodesFrontBottomLeftCorner = Z(mesh); nodesFrontBottomLeftCorner [mesh.nodesFrontBottomLeftCorner() ] = 1
nodesFrontBottomRightCorner = Z(mesh); nodesFrontBottomRightCorner[mesh.nodesFrontBottomRightCorner()] = 1
nodesFrontTopLeftCorner = Z(mesh); nodesFrontTopLeftCorner [mesh.nodesFrontTopLeftCorner() ] = 1
nodesFrontTopRightCorner = Z(mesh); nodesFrontTopRightCorner [mesh.nodesFrontTopRightCorner() ] = 1
nodesBackBottomLeftCorner = Z(mesh); nodesBackBottomLeftCorner [mesh.nodesBackBottomLeftCorner() ] = 1
nodesBackBottomRightCorner = Z(mesh); nodesBackBottomRightCorner [mesh.nodesBackBottomRightCorner() ] = 1
nodesBackTopLeftCorner = Z(mesh); nodesBackTopLeftCorner [mesh.nodesBackTopLeftCorner() ] = 1
nodesBackTopRightCorner = Z(mesh); nodesBackTopRightCorner [mesh.nodesBackTopRightCorner() ] = 1
# DOFs after eliminating periodicity
dofsPeriodic = mesh.dofsPeriodic()[:,0]
# open data file
f = h5py.File('MeshHex8-Regular-nodes_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
+# write nodal positions, and a dummy connectivity
f.create_dataset('/coor' ,data=mesh.coor() )
f.create_dataset('/conn' ,data=np.arange(mesh.nnode()) )
f.create_dataset('/nodesFront' ,data=nodesFront )
f.create_dataset('/nodesBack' ,data=nodesBack )
f.create_dataset('/nodesLeft' ,data=nodesLeft )
f.create_dataset('/nodesRight' ,data=nodesRight )
f.create_dataset('/nodesBottom' ,data=nodesBottom )
f.create_dataset('/nodesTop' ,data=nodesTop )
f.create_dataset('/nodesFrontFace' ,data=nodesFrontFace )
f.create_dataset('/nodesBackFace' ,data=nodesBackFace )
f.create_dataset('/nodesLeftFace' ,data=nodesLeftFace )
f.create_dataset('/nodesRightFace' ,data=nodesRightFace )
f.create_dataset('/nodesBottomFace' ,data=nodesBottomFace )
f.create_dataset('/nodesTopFace' ,data=nodesTopFace )
f.create_dataset('/nodesFrontBottomEdge' ,data=nodesFrontBottomEdge )
f.create_dataset('/nodesFrontTopEdge' ,data=nodesFrontTopEdge )
f.create_dataset('/nodesFrontLeftEdge' ,data=nodesFrontLeftEdge )
f.create_dataset('/nodesFrontRightEdge' ,data=nodesFrontRightEdge )
f.create_dataset('/nodesBackBottomEdge' ,data=nodesBackBottomEdge )
f.create_dataset('/nodesBackTopEdge' ,data=nodesBackTopEdge )
f.create_dataset('/nodesBackLeftEdge' ,data=nodesBackLeftEdge )
f.create_dataset('/nodesBackRightEdge' ,data=nodesBackRightEdge )
f.create_dataset('/nodesBottomLeftEdge' ,data=nodesBottomLeftEdge )
f.create_dataset('/nodesBottomRightEdge' ,data=nodesBottomRightEdge )
f.create_dataset('/nodesTopLeftEdge' ,data=nodesTopLeftEdge )
f.create_dataset('/nodesTopRightEdge' ,data=nodesTopRightEdge )
f.create_dataset('/nodesFrontBottomOpenEdge' ,data=nodesFrontBottomOpenEdge )
f.create_dataset('/nodesFrontTopOpenEdge' ,data=nodesFrontTopOpenEdge )
f.create_dataset('/nodesFrontLeftOpenEdge' ,data=nodesFrontLeftOpenEdge )
f.create_dataset('/nodesFrontRightOpenEdge' ,data=nodesFrontRightOpenEdge )
f.create_dataset('/nodesBackBottomOpenEdge' ,data=nodesBackBottomOpenEdge )
f.create_dataset('/nodesBackTopOpenEdge' ,data=nodesBackTopOpenEdge )
f.create_dataset('/nodesBackLeftOpenEdge' ,data=nodesBackLeftOpenEdge )
f.create_dataset('/nodesBackRightOpenEdge' ,data=nodesBackRightOpenEdge )
f.create_dataset('/nodesBottomLeftOpenEdge' ,data=nodesBottomLeftOpenEdge )
f.create_dataset('/nodesBottomRightOpenEdge' ,data=nodesBottomRightOpenEdge )
f.create_dataset('/nodesTopLeftOpenEdge' ,data=nodesTopLeftOpenEdge )
f.create_dataset('/nodesTopRightOpenEdge' ,data=nodesTopRightOpenEdge )
f.create_dataset('/nodesFrontBottomLeftCorner' ,data=nodesFrontBottomLeftCorner )
f.create_dataset('/nodesFrontBottomRightCorner',data=nodesFrontBottomRightCorner)
f.create_dataset('/nodesFrontTopLeftCorner' ,data=nodesFrontTopLeftCorner )
f.create_dataset('/nodesFrontTopRightCorner' ,data=nodesFrontTopRightCorner )
f.create_dataset('/nodesBackBottomLeftCorner' ,data=nodesBackBottomLeftCorner )
f.create_dataset('/nodesBackBottomRightCorner' ,data=nodesBackBottomRightCorner )
f.create_dataset('/nodesBackTopLeftCorner' ,data=nodesBackTopLeftCorner )
f.create_dataset('/nodesBackTopRightCorner' ,data=nodesBackTopRightCorner )
f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic )
# ======================================== write XDMF-file =========================================
+# basic file-format
xmf = '''
MeshHex8-Regular-nodes_paraview.hdf5:/conn
MeshHex8-Regular-nodes_paraview.hdf5:/coor
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFront
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBack
MeshHex8-Regular-nodes_paraview.hdf5:/nodesLeft
MeshHex8-Regular-nodes_paraview.hdf5:/nodesRight
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottom
MeshHex8-Regular-nodes_paraview.hdf5:/nodesTop
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontFace
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackFace
MeshHex8-Regular-nodes_paraview.hdf5:/nodesLeftFace
MeshHex8-Regular-nodes_paraview.hdf5:/nodesRightFace
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomFace
MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopFace
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontLeftEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontRightEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackLeftEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackRightEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomLeftEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomRightEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopLeftEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopRightEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontLeftOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontRightOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackLeftOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackRightOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomLeftOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomRightOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopLeftOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopRightOpenEdge
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomLeftCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomRightCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopLeftCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopRightCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomLeftCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomRightCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopLeftCorner
MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopRightCorner
MeshHex8-Regular-nodes_paraview.hdf5:/dofsPeriodic
'''
+# write to file, fill mesh dimensions
open('MeshHex8-Regular-nodes_paraview.xmf','w').write(xmf.format(
nnode = mesh.nnode(),
))
diff --git a/docs/examples/Mesh/MeshHex8-Regular_paraview.py b/docs/examples/Mesh/MeshHex8-Regular_paraview.py
index 3b32502..8a65f2d 100644
--- a/docs/examples/Mesh/MeshHex8-Regular_paraview.py
+++ b/docs/examples/Mesh/MeshHex8-Regular_paraview.py
@@ -1,44 +1,46 @@
import numpy as np
import h5py
import GooseFEM as gf
# ====================== create fictitious configuration + store to HDF5-file ======================
# create mesh object
mesh = gf.Mesh.Hex8.Regular(6,8,12)
# open data file
f = h5py.File('MeshHex8-Regular_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
+# write nodal coordinates and connectivity
f.create_dataset('/coor',data=mesh.coor())
f.create_dataset('/conn',data=mesh.conn())
# ======================================== write XDMF-file =========================================
+# basic file-format
xmf = '''
MeshHex8-Regular_paraview.hdf5:/conn
MeshHex8-Regular_paraview.hdf5:/coor
'''
+# write to file, fill mesh dimensions
open('MeshHex8-Regular_paraview.xmf','w').write(xmf.format(
nnode = mesh.nnode(),
nelem = mesh.nelem(),
nne = mesh.nne(),
))
diff --git a/docs/examples/Mesh/MeshQuad4-FineLayer-nodes_paraview.py b/docs/examples/Mesh/MeshQuad4-FineLayer-nodes_paraview.py
new file mode 100644
index 0000000..2398c14
--- /dev/null
+++ b/docs/examples/Mesh/MeshQuad4-FineLayer-nodes_paraview.py
@@ -0,0 +1,140 @@
+
+import numpy as np
+import h5py
+import GooseFEM as gf
+
+# ====================== create fictitious configuration + store to HDF5-file ======================
+
+# create mesh object
+mesh = gf.Mesh.Quad4.FineLayer(9,17)
+
+# set node identifier for all boundary sets
+Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int')
+nodesBottomEdge = Z(mesh); nodesBottomEdge [mesh.nodesBottomEdge() ] = 1
+nodesTopEdge = Z(mesh); nodesTopEdge [mesh.nodesTopEdge() ] = 1
+nodesLeftEdge = Z(mesh); nodesLeftEdge [mesh.nodesLeftEdge() ] = 1
+nodesRightEdge = Z(mesh); nodesRightEdge [mesh.nodesRightEdge() ] = 1
+nodesBottomOpenEdge = Z(mesh); nodesBottomOpenEdge [mesh.nodesBottomOpenEdge() ] = 1
+nodesTopOpenEdge = Z(mesh); nodesTopOpenEdge [mesh.nodesTopOpenEdge() ] = 1
+nodesLeftOpenEdge = Z(mesh); nodesLeftOpenEdge [mesh.nodesLeftOpenEdge() ] = 1
+nodesRightOpenEdge = Z(mesh); nodesRightOpenEdge [mesh.nodesRightOpenEdge() ] = 1
+nodesBottomLeftCorner = Z(mesh); nodesBottomLeftCorner [mesh.nodesBottomLeftCorner() ] = 1
+nodesBottomRightCorner = Z(mesh); nodesBottomRightCorner[mesh.nodesBottomRightCorner()] = 1
+nodesTopLeftCorner = Z(mesh); nodesTopLeftCorner [mesh.nodesTopLeftCorner() ] = 1
+nodesTopRightCorner = Z(mesh); nodesTopRightCorner [mesh.nodesTopRightCorner() ] = 1
+
+# DOFs after eliminating periodicity
+dofsPeriodic = mesh.dofsPeriodic()[:,0]
+
+# open data file
+f = h5py.File('MeshQuad4-FineLayer-nodes_paraview.hdf5','w')
+
+# write nodal positions, and a dummy connectivity
+f.create_dataset('/coor' ,data=mesh.coor() )
+f.create_dataset('/conn' ,data=np.arange(mesh.nnode()))
+f.create_dataset('/nodesBottomEdge' ,data=nodesBottomEdge )
+f.create_dataset('/nodesTopEdge' ,data=nodesTopEdge )
+f.create_dataset('/nodesLeftEdge' ,data=nodesLeftEdge )
+f.create_dataset('/nodesRightEdge' ,data=nodesRightEdge )
+f.create_dataset('/nodesBottomOpenEdge' ,data=nodesBottomOpenEdge )
+f.create_dataset('/nodesTopOpenEdge' ,data=nodesTopOpenEdge )
+f.create_dataset('/nodesLeftOpenEdge' ,data=nodesLeftOpenEdge )
+f.create_dataset('/nodesRightOpenEdge' ,data=nodesRightOpenEdge )
+f.create_dataset('/nodesBottomLeftCorner' ,data=nodesBottomLeftCorner )
+f.create_dataset('/nodesBottomRightCorner',data=nodesBottomRightCorner )
+f.create_dataset('/nodesTopLeftCorner' ,data=nodesTopLeftCorner )
+f.create_dataset('/nodesTopRightCorner' ,data=nodesTopRightCorner )
+f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic )
+
+# ======================================== write XDMF-file =========================================
+
+# basic file-format
+xmf = '''
+
+
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/conn
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/coor
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesLeftEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesRightEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomOpenEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopOpenEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesLeftOpenEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesRightOpenEdge
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomLeftCorner
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomRightCorner
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopLeftCorner
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopRightCorner
+
+
+
+
+ MeshQuad4-FineLayer-nodes_paraview.hdf5:/dofsPeriodic
+
+
+
+
+
+'''
+
+# write to file, fill mesh dimensions
+open('MeshQuad4-FineLayer-nodes_paraview.xmf','w').write(xmf.format(
+ nnode = mesh.nnode(),
+))
diff --git a/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py b/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py
index 3a027ef..4ec21e1 100644
--- a/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py
+++ b/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py
@@ -1,115 +1,56 @@
-# ========================================= load libraries =========================================
-
-import GooseFEM as gf
import numpy as np
import h5py
+import GooseFEM as gf
-# ========================================== define mesh ===========================================
-
-mesh = gf.Mesh.Quad4.FineLayer(36,36)
-coor = mesh.coor()
-conn = mesh.conn()
-
-# ======================================== write HDF5=file =========================================
+# ====================== create fictitious configuration + store to HDF5-file ======================
-# set file-name base
-name = 'MeshQuad4-FineLayer_paraview'
+# create mesh object
+mesh = gf.Mesh.Quad4.FineLayer(9,17)
# open data file
-f = h5py.File('{name:s}.hdf5'.format(name=name),'w')
+f = h5py.File('MeshQuad4-FineLayer_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
-f.create_dataset('/coor',data=coor)
-f.create_dataset('/conn',data=conn)
+# element set
+elementsMiddleLayer = np.zeros((mesh.nelem()),dtype='int')
+elementsMiddleLayer[mesh.elementsMiddleLayer()] = 1
-# create a sample deformation: simple shear
-for inc,gamma in enumerate(np.linspace(0,1,100)):
-
- # - initialize displacement (should be 3-d for ParaView)
- U = np.zeros((mesh.nnode(),3),dtype='float64')
- # - set
- U[:,0] += gamma * coor[:,1]
- # - store
- f.create_dataset('/disp/%d'%inc,data=U)
+# write nodal coordinates and connectivity
+f.create_dataset('/coor' ,data=mesh.coor() )
+f.create_dataset('/conn' ,data=mesh.conn() )
+f.create_dataset('/elementsMiddleLayer',data=elementsMiddleLayer)
# ======================================== write XDMF-file =========================================
-# --------------------------------- format of the main structure ----------------------------------
-
+# basic file-format
xmf = '''
-
-{series:s}
-
+
+
+
+ MeshQuad4-FineLayer_paraview.hdf5:/conn
+
+
+
+
+ MeshQuad4-FineLayer_paraview.hdf5:/coor
+
+
+
+
+ MeshQuad4-FineLayer_paraview.hdf5:/elementsMiddleLayer
+
+
+
'''
-# ----------------------------- format of an increment in time-series ------------------------------
-
-grid = '''
-
-
-
- {name:s}.hdf5:/conn
-
-
-
-
- {name:s}.hdf5:/coor
-
-
-
-
- {name:s}.hdf5:/disp/{inc:d}
-
-
-
-'''
-
-# ------------------------------------------- write file -------------------------------------------
-
-# initialize string that will contain the full time series
-txt = ''
-
-# loop over all increments, append the time series
-for inc in range(100):
- txt += grid.format(
- name = name,
- inc = inc,
- nnode = mesh.nnode(),
- ndim = mesh.ndim(),
- nelem = mesh.nelem(),
- nne = mesh.nne()
- )
-
-# write xmf-file, fix the indentation
-fname = '{name:s}.xmf'.format(name=name)
-open(fname,'w').write(xmf.format(series=' '+txt.replace('\n','\n ')))
-
-# ======================================== print node=sets =========================================
-
-np.set_printoptions(linewidth=200)
-
-print('nodesBottomEdge =', mesh.nodesBottomEdge())
-print('nodesTopEdge =', mesh.nodesTopEdge())
-print('nodesLeftEdge =', mesh.nodesLeftEdge())
-print('nodesRightEdge =', mesh.nodesRightEdge())
-print('nodesLeftBottomCorner =', mesh.nodesLeftBottomCorner())
-print('nodesRightBottomCorner =', mesh.nodesRightBottomCorner())
-print('nodesRightTopCorner =', mesh.nodesRightTopCorner())
-print('nodesLeftTopCorner =', mesh.nodesLeftTopCorner())
-
-print('nodesPeriodic =')
-print(mesh.nodesPeriodic())
-print('nodesOrigin =')
-print(mesh.nodesOrigin())
-
-# check uniqueness of tied nodes
-ties = mesh.nodesPeriodic()[:,1]
-
-if len(np.unique(ties)) != len(ties):
- raise IOError('Nodal ties not unique')
+# write to file, fill mesh dimensions
+open('MeshQuad4-FineLayer_paraview.xmf','w').write(xmf.format(
+ nnode = mesh.nnode(),
+ nelem = mesh.nelem(),
+ nne = mesh.nne(),
+))
diff --git a/docs/examples/Mesh/MeshQuad4-Regular-nodes_paraview.py b/docs/examples/Mesh/MeshQuad4-Regular-nodes_paraview.py
new file mode 100644
index 0000000..ed1e6d6
--- /dev/null
+++ b/docs/examples/Mesh/MeshQuad4-Regular-nodes_paraview.py
@@ -0,0 +1,140 @@
+
+import numpy as np
+import h5py
+import GooseFEM as gf
+
+# ====================== create fictitious configuration + store to HDF5-file ======================
+
+# create mesh object
+mesh = gf.Mesh.Quad4.Regular(9,11)
+
+# set node identifier for all boundary sets
+Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int')
+nodesBottomEdge = Z(mesh); nodesBottomEdge [mesh.nodesBottomEdge() ] = 1
+nodesTopEdge = Z(mesh); nodesTopEdge [mesh.nodesTopEdge() ] = 1
+nodesLeftEdge = Z(mesh); nodesLeftEdge [mesh.nodesLeftEdge() ] = 1
+nodesRightEdge = Z(mesh); nodesRightEdge [mesh.nodesRightEdge() ] = 1
+nodesBottomOpenEdge = Z(mesh); nodesBottomOpenEdge [mesh.nodesBottomOpenEdge() ] = 1
+nodesTopOpenEdge = Z(mesh); nodesTopOpenEdge [mesh.nodesTopOpenEdge() ] = 1
+nodesLeftOpenEdge = Z(mesh); nodesLeftOpenEdge [mesh.nodesLeftOpenEdge() ] = 1
+nodesRightOpenEdge = Z(mesh); nodesRightOpenEdge [mesh.nodesRightOpenEdge() ] = 1
+nodesBottomLeftCorner = Z(mesh); nodesBottomLeftCorner [mesh.nodesBottomLeftCorner() ] = 1
+nodesBottomRightCorner = Z(mesh); nodesBottomRightCorner[mesh.nodesBottomRightCorner()] = 1
+nodesTopLeftCorner = Z(mesh); nodesTopLeftCorner [mesh.nodesTopLeftCorner() ] = 1
+nodesTopRightCorner = Z(mesh); nodesTopRightCorner [mesh.nodesTopRightCorner() ] = 1
+
+# DOFs after eliminating periodicity
+dofsPeriodic = mesh.dofsPeriodic()[:,0]
+
+# open data file
+f = h5py.File('MeshQuad4-Regular-nodes_paraview.hdf5','w')
+
+# write nodal positions, and a dummy connectivity
+f.create_dataset('/coor' ,data=mesh.coor() )
+f.create_dataset('/conn' ,data=np.arange(mesh.nnode()))
+f.create_dataset('/nodesBottomEdge' ,data=nodesBottomEdge )
+f.create_dataset('/nodesTopEdge' ,data=nodesTopEdge )
+f.create_dataset('/nodesLeftEdge' ,data=nodesLeftEdge )
+f.create_dataset('/nodesRightEdge' ,data=nodesRightEdge )
+f.create_dataset('/nodesBottomOpenEdge' ,data=nodesBottomOpenEdge )
+f.create_dataset('/nodesTopOpenEdge' ,data=nodesTopOpenEdge )
+f.create_dataset('/nodesLeftOpenEdge' ,data=nodesLeftOpenEdge )
+f.create_dataset('/nodesRightOpenEdge' ,data=nodesRightOpenEdge )
+f.create_dataset('/nodesBottomLeftCorner' ,data=nodesBottomLeftCorner )
+f.create_dataset('/nodesBottomRightCorner',data=nodesBottomRightCorner )
+f.create_dataset('/nodesTopLeftCorner' ,data=nodesTopLeftCorner )
+f.create_dataset('/nodesTopRightCorner' ,data=nodesTopRightCorner )
+f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic )
+
+# ======================================== write XDMF-file =========================================
+
+# basic file-format
+xmf = '''
+
+
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/conn
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/coor
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesLeftEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesRightEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomOpenEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopOpenEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesLeftOpenEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesRightOpenEdge
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomLeftCorner
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomRightCorner
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopLeftCorner
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopRightCorner
+
+
+
+
+ MeshQuad4-Regular-nodes_paraview.hdf5:/dofsPeriodic
+
+
+
+
+
+'''
+
+# write to file, fill mesh dimensions
+open('MeshQuad4-Regular-nodes_paraview.xmf','w').write(xmf.format(
+ nnode = mesh.nnode(),
+))
diff --git a/docs/examples/Mesh/MeshQuad4-Regular_paraview.py b/docs/examples/Mesh/MeshQuad4-Regular_paraview.py
index 78ba252..3c58990 100644
--- a/docs/examples/Mesh/MeshQuad4-Regular_paraview.py
+++ b/docs/examples/Mesh/MeshQuad4-Regular_paraview.py
@@ -1,115 +1,46 @@
-# ========================================= load libraries =========================================
-
-import GooseFEM as gf
import numpy as np
import h5py
+import GooseFEM as gf
-# ========================================== define mesh ===========================================
-
-mesh = gf.Mesh.Quad4.Regular(4,5)
-coor = mesh.coor()
-conn = mesh.conn()
-
-# ======================================== write HDF5=file =========================================
+# ====================== create fictitious configuration + store to HDF5-file ======================
-# set file-name base
-name = 'MeshQuad4-Regular_paraview'
+# create mesh object
+mesh = gf.Mesh.Quad4.Regular(9,11)
# open data file
-f = h5py.File('{name:s}.hdf5'.format(name=name),'w')
+f = h5py.File('MeshQuad4-Regular_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
-f.create_dataset('/coor',data=coor)
-f.create_dataset('/conn',data=conn)
-
-# create a sample deformation: simple shear
-for inc,gamma in enumerate(np.linspace(0,1,100)):
-
- # - initialize displacement (should be 3-d for ParaView)
- U = np.zeros((mesh.nnode(),3),dtype='float64')
- # - set
- U[:,0] += gamma * coor[:,1]
- # - store
- f.create_dataset('/disp/%d'%inc,data=U)
+# write nodal coordinates and connectivity
+f.create_dataset('/coor',data=mesh.coor())
+f.create_dataset('/conn',data=mesh.conn())
# ======================================== write XDMF-file =========================================
-# --------------------------------- format of the main structure ----------------------------------
-
+# basic file-format
xmf = '''
-
-{series:s}
-
+
+
+
+ MeshQuad4-Regular_paraview.hdf5:/conn
+
+
+
+
+ MeshQuad4-Regular_paraview.hdf5:/coor
+
+
+
'''
-# ----------------------------- format of an increment in time-series ------------------------------
-
-grid = '''
-
-
-
- {name:s}.hdf5:/conn
-
-
-
-
- {name:s}.hdf5:/coor
-
-
-
-
- {name:s}.hdf5:/disp/{inc:d}
-
-
-
-'''
-
-# ------------------------------------------- write file -------------------------------------------
-
-# initialize string that will contain the full time series
-txt = ''
-
-# loop over all increments, append the time series
-for inc in range(100):
- txt += grid.format(
- name = name,
- inc = inc,
- nnode = mesh.nnode(),
- ndim = mesh.ndim(),
- nelem = mesh.nelem(),
- nne = mesh.nne()
- )
-
-# write xmf-file, fix the indentation
-fname = '{name:s}.xmf'.format(name=name)
-open(fname,'w').write(xmf.format(series=' '+txt.replace('\n','\n ')))
-
-# ======================================== print node=sets =========================================
-
-np.set_printoptions(linewidth=200)
-
-print('nodesBottomEdge =', mesh.nodesBottomEdge())
-print('nodesTopEdge =', mesh.nodesTopEdge())
-print('nodesLeftEdge =', mesh.nodesLeftEdge())
-print('nodesRightEdge =', mesh.nodesRightEdge())
-print('nodesLeftBottomCorner =', mesh.nodesLeftBottomCorner())
-print('nodesRightBottomCorner =', mesh.nodesRightBottomCorner())
-print('nodesRightTopCorner =', mesh.nodesRightTopCorner())
-print('nodesLeftTopCorner =', mesh.nodesLeftTopCorner())
-
-print('nodesPeriodic =')
-print(mesh.nodesPeriodic())
-print('nodesOrigin =')
-print(mesh.nodesOrigin())
-
-# check uniqueness of tied nodes
-ties = mesh.nodesPeriodic()[:,1]
-
-if len(np.unique(ties)) != len(ties):
- raise IOError('Nodal ties not unique')
+# write to file, fill mesh dimensions
+open('MeshQuad4-Regular_paraview.xmf','w').write(xmf.format(
+ nnode = mesh.nnode(),
+ nelem = mesh.nelem(),
+ nne = mesh.nne(),
+))
diff --git a/docs/examples/Mesh/MeshTri3-Regular-nodes_paraview.py b/docs/examples/Mesh/MeshTri3-Regular-nodes_paraview.py
new file mode 100644
index 0000000..b320616
--- /dev/null
+++ b/docs/examples/Mesh/MeshTri3-Regular-nodes_paraview.py
@@ -0,0 +1,140 @@
+
+import numpy as np
+import h5py
+import GooseFEM as gf
+
+# ====================== create fictitious configuration + store to HDF5-file ======================
+
+# create mesh object
+mesh = gf.Mesh.Tri3.Regular(9,11)
+
+# set node identifier for all boundary sets
+Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int')
+nodesBottomEdge = Z(mesh); nodesBottomEdge [mesh.nodesBottomEdge() ] = 1
+nodesTopEdge = Z(mesh); nodesTopEdge [mesh.nodesTopEdge() ] = 1
+nodesLeftEdge = Z(mesh); nodesLeftEdge [mesh.nodesLeftEdge() ] = 1
+nodesRightEdge = Z(mesh); nodesRightEdge [mesh.nodesRightEdge() ] = 1
+nodesBottomOpenEdge = Z(mesh); nodesBottomOpenEdge [mesh.nodesBottomOpenEdge() ] = 1
+nodesTopOpenEdge = Z(mesh); nodesTopOpenEdge [mesh.nodesTopOpenEdge() ] = 1
+nodesLeftOpenEdge = Z(mesh); nodesLeftOpenEdge [mesh.nodesLeftOpenEdge() ] = 1
+nodesRightOpenEdge = Z(mesh); nodesRightOpenEdge [mesh.nodesRightOpenEdge() ] = 1
+nodesBottomLeftCorner = Z(mesh); nodesBottomLeftCorner [mesh.nodesBottomLeftCorner() ] = 1
+nodesBottomRightCorner = Z(mesh); nodesBottomRightCorner[mesh.nodesBottomRightCorner()] = 1
+nodesTopLeftCorner = Z(mesh); nodesTopLeftCorner [mesh.nodesTopLeftCorner() ] = 1
+nodesTopRightCorner = Z(mesh); nodesTopRightCorner [mesh.nodesTopRightCorner() ] = 1
+
+# DOFs after eliminating periodicity
+dofsPeriodic = mesh.dofsPeriodic()[:,0]
+
+# open data file
+f = h5py.File('MeshTri3-Regular-nodes_paraview.hdf5','w')
+
+# write nodal positions, and a dummy connectivity
+f.create_dataset('/coor' ,data=mesh.coor() )
+f.create_dataset('/conn' ,data=np.arange(mesh.nnode()))
+f.create_dataset('/nodesBottomEdge' ,data=nodesBottomEdge )
+f.create_dataset('/nodesTopEdge' ,data=nodesTopEdge )
+f.create_dataset('/nodesLeftEdge' ,data=nodesLeftEdge )
+f.create_dataset('/nodesRightEdge' ,data=nodesRightEdge )
+f.create_dataset('/nodesBottomOpenEdge' ,data=nodesBottomOpenEdge )
+f.create_dataset('/nodesTopOpenEdge' ,data=nodesTopOpenEdge )
+f.create_dataset('/nodesLeftOpenEdge' ,data=nodesLeftOpenEdge )
+f.create_dataset('/nodesRightOpenEdge' ,data=nodesRightOpenEdge )
+f.create_dataset('/nodesBottomLeftCorner' ,data=nodesBottomLeftCorner )
+f.create_dataset('/nodesBottomRightCorner',data=nodesBottomRightCorner )
+f.create_dataset('/nodesTopLeftCorner' ,data=nodesTopLeftCorner )
+f.create_dataset('/nodesTopRightCorner' ,data=nodesTopRightCorner )
+f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic )
+
+# ======================================== write XDMF-file =========================================
+
+# basic file-format
+xmf = '''
+
+
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/conn
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/coor
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesLeftEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesRightEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomOpenEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopOpenEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesLeftOpenEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesRightOpenEdge
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomLeftCorner
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomRightCorner
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopLeftCorner
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopRightCorner
+
+
+
+
+ MeshTri3-Regular-nodes_paraview.hdf5:/dofsPeriodic
+
+
+
+
+
+'''
+
+# write to file, fill mesh dimensions
+open('MeshTri3-Regular-nodes_paraview.xmf','w').write(xmf.format(
+ nnode = mesh.nnode(),
+))
diff --git a/docs/examples/Mesh/MeshTri3-Regular_paraview.py b/docs/examples/Mesh/MeshTri3-Regular_paraview.py
index 1e95a2f..66cc261 100644
--- a/docs/examples/Mesh/MeshTri3-Regular_paraview.py
+++ b/docs/examples/Mesh/MeshTri3-Regular_paraview.py
@@ -1,115 +1,46 @@
-# ========================================= load libraries =========================================
-
-import GooseFEM as gf
import numpy as np
import h5py
+import GooseFEM as gf
-# ========================================== define mesh ===========================================
-
-mesh = gf.Mesh.Tri3.Regular(4,5)
-coor = mesh.coor()
-conn = mesh.conn()
-
-# ======================================== write HDF5=file =========================================
+# ====================== create fictitious configuration + store to HDF5-file ======================
-# set file-name base
-name = 'MeshTri3-Regular_paraview'
+# create mesh object
+mesh = gf.Mesh.Tri3.Regular(9,11)
# open data file
-f = h5py.File('{name:s}.hdf5'.format(name=name),'w')
+f = h5py.File('MeshTri3-Regular_paraview.hdf5','w')
-# write particle positions, and a dummy connectivity
-f.create_dataset('/coor',data=coor)
-f.create_dataset('/conn',data=conn)
-
-# create a sample deformation: simple shear
-for inc,gamma in enumerate(np.linspace(0,1,100)):
-
- # - initialize displacement (should be 3-d for ParaView)
- U = np.zeros((mesh.nnode(),3),dtype='float64')
- # - set
- U[:,0] += gamma * coor[:,1]
- # - store
- f.create_dataset('/disp/%d'%inc,data=U)
+# write nodal coordinates and connectivity
+f.create_dataset('/coor',data=mesh.coor())
+f.create_dataset('/conn',data=mesh.conn())
# ======================================== write XDMF-file =========================================
-# --------------------------------- format of the main structure ----------------------------------
-
+# basic file-format
xmf = '''
-
-{series:s}
-
+
+
+
+ MeshTri3-Regular_paraview.hdf5:/conn
+
+
+
+
+ MeshTri3-Regular_paraview.hdf5:/coor
+
+
+
'''
-# ----------------------------- format of an increment in time-series ------------------------------
-
-grid = '''
-
-
-
- {name:s}.hdf5:/conn
-
-
-
-
- {name:s}.hdf5:/coor
-
-
-
-
- {name:s}.hdf5:/disp/{inc:d}
-
-
-
-'''
-
-# ------------------------------------------- write file -------------------------------------------
-
-# initialize string that will contain the full time series
-txt = ''
-
-# loop over all increments, append the time series
-for inc in range(100):
- txt += grid.format(
- name = name,
- inc = inc,
- nnode = mesh.nnode(),
- ndim = mesh.ndim(),
- nelem = mesh.nelem(),
- nne = mesh.nne()
- )
-
-# write xmf-file, fix the indentation
-fname = '{name:s}.xmf'.format(name=name)
-open(fname,'w').write(xmf.format(series=' '+txt.replace('\n','\n ')))
-
-# ======================================== print node=sets =========================================
-
-np.set_printoptions(linewidth=200)
-
-print('nodesBottomEdge =', mesh.nodesBottomEdge())
-print('nodesTopEdge =', mesh.nodesTopEdge())
-print('nodesLeftEdge =', mesh.nodesLeftEdge())
-print('nodesRightEdge =', mesh.nodesRightEdge())
-print('nodesLeftBottomCorner =', mesh.nodesLeftBottomCorner())
-print('nodesRightBottomCorner =', mesh.nodesRightBottomCorner())
-print('nodesRightTopCorner =', mesh.nodesRightTopCorner())
-print('nodesLeftTopCorner =', mesh.nodesLeftTopCorner())
-
-print('nodesPeriodic =')
-print(mesh.nodesPeriodic())
-print('nodesOrigin =')
-print(mesh.nodesOrigin())
-
-# check uniqueness of tied nodes
-ties = mesh.nodesPeriodic()[:,1]
-
-if len(np.unique(ties)) != len(ties):
- raise IOError('Nodal ties not unique')
+# write to file, fill mesh dimensions
+open('MeshTri3-Regular_paraview.xmf','w').write(xmf.format(
+ nnode = mesh.nnode(),
+ nelem = mesh.nelem(),
+ nne = mesh.nne(),
+))
diff --git a/python/python_interface.cpp b/python/python_interface.cpp
index 3bc3025..a01ceea 100644
--- a/python/python_interface.cpp
+++ b/python/python_interface.cpp
@@ -1,492 +1,505 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#include
#include
#include
#include "../src/GooseFEM/GooseFEM.h"
// alias for short-hand notation below
namespace py = pybind11;
PYBIND11_MODULE(GooseFEM, m) {
// =================================================================================================
m.doc() = "Some simple finite element meshes and operations";
// define submodules "mXXX"
py::module mMesh = m .def_submodule("Mesh" , "Generic mesh routines" );
py::module mMeshTri3 = mMesh.def_submodule("Tri3" , "Linear triangular elements (2D)" );
py::module mMeshQuad4 = mMesh.def_submodule("Quad4" , "Linear quadrilateral elements (2D)" );
py::module mMeshHex8 = mMesh.def_submodule("Hex8" , "Linear hexahedron (brick) elements (3D)");
// ======================================= GooseFEM/Mesh.h ========================================
mMesh.def("elem2node",
&GooseFEM::Mesh::elem2node,
"Elements connect to each node: [ number of elements , element numbers ]",
py::arg("conn")
);
mMesh.def("dofs",
&GooseFEM::Mesh::dofs,
"List with DOF-numbers (in sequential order)",
py::arg("nnode"),
py::arg("ndim")
);
mMesh.def("renumber",
py::overload_cast(&GooseFEM::Mesh::renumber),
"Renumber DOF-list to use the lowest possible index",
py::arg("dofs")
);
mMesh.def("renumber",
py::overload_cast(&GooseFEM::Mesh::renumber),
"Renumber DOF-list to begin or end with 'idx'",
py::arg("dofs"),
py::arg("idx"),
py::arg("location")="end"
);
// ====================================== GooseFEM/MeshHex8.h ======================================
py::class_(mMeshHex8, "Regular")
// constructor
.def(
py::init(),
"mesh with nx*ny*nz 'pixels' and edge size h",
py::arg("nx"),
py::arg("ny"),
py::arg("nz"),
py::arg("h")=1.
)
// sizes
.def("nelem" , &GooseFEM::Mesh::Hex8::Regular::nelem )
.def("nnode" , &GooseFEM::Mesh::Hex8::Regular::nnode )
.def("nne" , &GooseFEM::Mesh::Hex8::Regular::nne )
.def("ndim" , &GooseFEM::Mesh::Hex8::Regular::ndim )
// mesh
.def("coor" , &GooseFEM::Mesh::Hex8::Regular::coor )
.def("conn" , &GooseFEM::Mesh::Hex8::Regular::conn )
// boundary nodes: planes
.def("nodesFront" , &GooseFEM::Mesh::Hex8::Regular::nodesFront )
.def("nodesBack" , &GooseFEM::Mesh::Hex8::Regular::nodesBack )
.def("nodesLeft" , &GooseFEM::Mesh::Hex8::Regular::nodesLeft )
.def("nodesRight" , &GooseFEM::Mesh::Hex8::Regular::nodesRight )
.def("nodesBottom" , &GooseFEM::Mesh::Hex8::Regular::nodesBottom )
.def("nodesTop" , &GooseFEM::Mesh::Hex8::Regular::nodesTop )
// boundary nodes: faces
.def("nodesFrontFace" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontFace )
.def("nodesBackFace" , &GooseFEM::Mesh::Hex8::Regular::nodesBackFace )
.def("nodesLeftFace" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFace )
.def("nodesRightFace" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFace )
.def("nodesBottomFace" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFace )
.def("nodesTopFace" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFace )
// boundary nodes: edges
.def("nodesFrontBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomEdge )
.def("nodesFrontTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopEdge )
.def("nodesFrontLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftEdge )
.def("nodesFrontRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightEdge )
.def("nodesBackBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomEdge )
.def("nodesBackTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopEdge )
.def("nodesBackLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftEdge )
.def("nodesBackRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightEdge )
.def("nodesBottomLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftEdge )
.def("nodesBottomRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightEdge )
.def("nodesTopLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftEdge )
.def("nodesTopRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightEdge )
// boundary nodes: faces (aliases)
.def("nodesBottomFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontEdge )
.def("nodesBottomBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackEdge )
.def("nodesTopFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontEdge )
.def("nodesTopBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackEdge )
.def("nodesLeftBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomEdge )
.def("nodesLeftFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontEdge )
.def("nodesLeftBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackEdge )
.def("nodesLeftTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopEdge )
.def("nodesRightBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomEdge )
.def("nodesRightTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopEdge )
.def("nodesRightFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontEdge )
.def("nodesRightBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackEdge )
// boundary nodes: edges, without corners
.def("nodesFrontBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomOpenEdge )
.def("nodesFrontTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopOpenEdge )
.def("nodesFrontLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftOpenEdge )
.def("nodesFrontRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightOpenEdge )
.def("nodesBackBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomOpenEdge )
.def("nodesBackTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopOpenEdge )
.def("nodesBackLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftOpenEdge )
.def("nodesBackRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightOpenEdge )
.def("nodesBottomLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftOpenEdge )
.def("nodesBottomRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightOpenEdge )
.def("nodesTopLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftOpenEdge )
.def("nodesTopRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightOpenEdge )
// boundary nodes: edges, without corners (aliases)
.def("nodesBottomFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontOpenEdge )
.def("nodesBottomBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackOpenEdge )
.def("nodesTopFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontOpenEdge )
.def("nodesTopBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackOpenEdge )
.def("nodesLeftBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomOpenEdge )
.def("nodesLeftFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontOpenEdge )
.def("nodesLeftBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackOpenEdge )
.def("nodesLeftTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopOpenEdge )
.def("nodesRightBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomOpenEdge )
.def("nodesRightTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopOpenEdge )
.def("nodesRightFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontOpenEdge )
.def("nodesRightBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackOpenEdge )
// boundary nodes: corners
.def("nodesFrontBottomLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomLeftCorner )
.def("nodesFrontBottomRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomRightCorner)
.def("nodesFrontTopLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopLeftCorner )
.def("nodesFrontTopRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopRightCorner )
.def("nodesBackBottomLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomLeftCorner )
.def("nodesBackBottomRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomRightCorner )
.def("nodesBackTopLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopLeftCorner )
.def("nodesBackTopRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopRightCorner )
// boundary nodes: corners (aliases)
.def("nodesFrontLeftBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftBottomCorner )
.def("nodesBottomFrontLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontLeftCorner )
.def("nodesBottomLeftFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftFrontCorner )
.def("nodesLeftFrontBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontBottomCorner )
.def("nodesLeftBottomFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomFrontCorner )
.def("nodesFrontRightBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightBottomCorner)
.def("nodesBottomFrontRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontRightCorner)
.def("nodesBottomRightFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightFrontCorner)
.def("nodesRightFrontBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontBottomCorner)
.def("nodesRightBottomFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomFrontCorner)
.def("nodesFrontLeftTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftTopCorner )
.def("nodesTopFrontLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontLeftCorner )
.def("nodesTopLeftFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftFrontCorner )
.def("nodesLeftFrontTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontTopCorner )
.def("nodesLeftTopFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopFrontCorner )
.def("nodesFrontRightTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightTopCorner )
.def("nodesTopFrontRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontRightCorner )
.def("nodesTopRightFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightFrontCorner )
.def("nodesRightFrontTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontTopCorner )
.def("nodesRightTopFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopFrontCorner )
.def("nodesBackLeftBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftBottomCorner )
.def("nodesBottomBackLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackLeftCorner )
.def("nodesBottomLeftBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftBackCorner )
.def("nodesLeftBackBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackBottomCorner )
.def("nodesLeftBottomBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomBackCorner )
.def("nodesBackRightBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightBottomCorner )
.def("nodesBottomBackRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackRightCorner )
.def("nodesBottomRightBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightBackCorner )
.def("nodesRightBackBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackBottomCorner )
.def("nodesRightBottomBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomBackCorner )
.def("nodesBackLeftTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftTopCorner )
.def("nodesTopBackLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackLeftCorner )
.def("nodesTopLeftBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftBackCorner )
.def("nodesLeftBackTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackTopCorner )
.def("nodesLeftTopBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopBackCorner )
.def("nodesBackRightTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightTopCorner )
.def("nodesTopBackRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackRightCorner )
.def("nodesTopRightBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightBackCorner )
.def("nodesRightBackTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackTopCorner )
.def("nodesRightTopBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopBackCorner )
// periodicity
.def("nodesPeriodic" , &GooseFEM::Mesh::Hex8::Regular::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Hex8::Regular::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Hex8::Regular::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Hex8::Regular::dofsPeriodic )
// print to screen
.def("__repr__",
[](const GooseFEM::Mesh::Hex8::Regular &a){ return ""; }
);
// -------------------------------------------------------------------------------------------------
py::class_(mMeshHex8, "FineLayer")
// constructor
.def(
py::init(),
"mesh with nx*ny*nz 'pixels' and edge size h",
py::arg("nx"),
py::arg("ny"),
py::arg("nz"),
py::arg("h")=1.,
py::arg("nfine")=1
)
// sizes
.def("nelem" , &GooseFEM::Mesh::Hex8::FineLayer::nelem )
.def("nnode" , &GooseFEM::Mesh::Hex8::FineLayer::nnode )
.def("nne" , &GooseFEM::Mesh::Hex8::FineLayer::nne )
.def("ndim" , &GooseFEM::Mesh::Hex8::FineLayer::ndim )
.def("shape" , &GooseFEM::Mesh::Hex8::FineLayer::shape )
// mesh
.def("coor" , &GooseFEM::Mesh::Hex8::FineLayer::coor )
.def("conn" , &GooseFEM::Mesh::Hex8::FineLayer::conn )
+ // element sets
+ .def("elementsMiddleLayer" , &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer )
// boundary nodes: planes
.def("nodesFront" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFront )
.def("nodesBack" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBack )
.def("nodesLeft" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeft )
.def("nodesRight" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRight )
.def("nodesBottom" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottom )
.def("nodesTop" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTop )
// boundary nodes: faces
.def("nodesFrontFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontFace )
.def("nodesBackFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackFace )
.def("nodesLeftFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFace )
.def("nodesRightFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFace )
.def("nodesBottomFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFace )
.def("nodesTopFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFace )
// boundary nodes: edges
.def("nodesFrontBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomEdge )
.def("nodesFrontTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopEdge )
.def("nodesFrontLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftEdge )
.def("nodesFrontRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightEdge )
.def("nodesBackBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomEdge )
.def("nodesBackTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopEdge )
.def("nodesBackLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftEdge )
.def("nodesBackRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightEdge )
.def("nodesBottomLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftEdge )
.def("nodesBottomRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightEdge )
.def("nodesTopLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftEdge )
.def("nodesTopRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightEdge )
// boundary nodes: faces (aliases)
.def("nodesBottomFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontEdge )
.def("nodesBottomBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackEdge )
.def("nodesTopFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontEdge )
.def("nodesTopBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackEdge )
.def("nodesLeftBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomEdge )
.def("nodesLeftFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontEdge )
.def("nodesLeftBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackEdge )
.def("nodesLeftTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopEdge )
.def("nodesRightBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomEdge )
.def("nodesRightTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopEdge )
.def("nodesRightFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontEdge )
.def("nodesRightBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackEdge )
// boundary nodes: edges, without corners
.def("nodesFrontBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomOpenEdge )
.def("nodesFrontTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopOpenEdge )
.def("nodesFrontLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftOpenEdge )
.def("nodesFrontRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightOpenEdge )
.def("nodesBackBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomOpenEdge )
.def("nodesBackTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopOpenEdge )
.def("nodesBackLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftOpenEdge )
.def("nodesBackRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightOpenEdge )
.def("nodesBottomLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftOpenEdge )
.def("nodesBottomRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightOpenEdge )
.def("nodesTopLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftOpenEdge )
.def("nodesTopRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightOpenEdge )
// boundary nodes: edges, without corners (aliases)
.def("nodesBottomFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontOpenEdge )
.def("nodesBottomBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackOpenEdge )
.def("nodesTopFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontOpenEdge )
.def("nodesTopBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackOpenEdge )
.def("nodesLeftBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomOpenEdge )
.def("nodesLeftFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontOpenEdge )
.def("nodesLeftBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackOpenEdge )
.def("nodesLeftTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopOpenEdge )
.def("nodesRightBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomOpenEdge )
.def("nodesRightTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopOpenEdge )
.def("nodesRightFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontOpenEdge )
.def("nodesRightBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackOpenEdge )
// boundary nodes: corners
.def("nodesFrontBottomLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomLeftCorner )
.def("nodesFrontBottomRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomRightCorner)
.def("nodesFrontTopLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopLeftCorner )
.def("nodesFrontTopRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopRightCorner )
.def("nodesBackBottomLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomLeftCorner )
.def("nodesBackBottomRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomRightCorner )
.def("nodesBackTopLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopLeftCorner )
.def("nodesBackTopRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopRightCorner )
// boundary nodes: corners (aliases)
.def("nodesFrontLeftBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftBottomCorner )
.def("nodesBottomFrontLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontLeftCorner )
.def("nodesBottomLeftFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftFrontCorner )
.def("nodesLeftFrontBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontBottomCorner )
.def("nodesLeftBottomFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomFrontCorner )
.def("nodesFrontRightBottomCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightBottomCorner)
.def("nodesBottomFrontRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontRightCorner)
.def("nodesBottomRightFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightFrontCorner)
.def("nodesRightFrontBottomCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontBottomCorner)
.def("nodesRightBottomFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomFrontCorner)
.def("nodesFrontLeftTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftTopCorner )
.def("nodesTopFrontLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontLeftCorner )
.def("nodesTopLeftFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftFrontCorner )
.def("nodesLeftFrontTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontTopCorner )
.def("nodesLeftTopFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopFrontCorner )
.def("nodesFrontRightTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightTopCorner )
.def("nodesTopFrontRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontRightCorner )
.def("nodesTopRightFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightFrontCorner )
.def("nodesRightFrontTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontTopCorner )
.def("nodesRightTopFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopFrontCorner )
.def("nodesBackLeftBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftBottomCorner )
.def("nodesBottomBackLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackLeftCorner )
.def("nodesBottomLeftBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftBackCorner )
.def("nodesLeftBackBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackBottomCorner )
.def("nodesLeftBottomBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomBackCorner )
.def("nodesBackRightBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightBottomCorner )
.def("nodesBottomBackRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackRightCorner )
.def("nodesBottomRightBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightBackCorner )
.def("nodesRightBackBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackBottomCorner )
.def("nodesRightBottomBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomBackCorner )
.def("nodesBackLeftTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftTopCorner )
.def("nodesTopBackLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackLeftCorner )
.def("nodesTopLeftBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftBackCorner )
.def("nodesLeftBackTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackTopCorner )
.def("nodesLeftTopBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopBackCorner )
.def("nodesBackRightTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightTopCorner )
.def("nodesTopBackRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackRightCorner )
.def("nodesTopRightBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightBackCorner )
.def("nodesRightBackTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackTopCorner )
.def("nodesRightTopBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopBackCorner )
// periodicity
.def("nodesPeriodic" , &GooseFEM::Mesh::Hex8::FineLayer::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Hex8::FineLayer::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Hex8::FineLayer::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Hex8::FineLayer::dofsPeriodic )
// print to screen
.def("__repr__",
[](const GooseFEM::Mesh::Hex8::FineLayer &a){ return ""; }
);
// ===================================== GooseFEM/MeshQuad4.h =====================================
py::class_(mMeshQuad4, "Regular")
.def(
py::init(),
"Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'",
py::arg("nx"),
py::arg("ny"),
py::arg("h")=1.
)
.def("coor" , &GooseFEM::Mesh::Quad4::Regular::coor )
.def("conn" , &GooseFEM::Mesh::Quad4::Regular::conn )
.def("nelem" , &GooseFEM::Mesh::Quad4::Regular::nelem )
.def("nnode" , &GooseFEM::Mesh::Quad4::Regular::nnode )
.def("nne" , &GooseFEM::Mesh::Quad4::Regular::nne )
.def("ndim" , &GooseFEM::Mesh::Quad4::Regular::ndim )
.def("nodesBottomEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge )
.def("nodesTopEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge )
.def("nodesLeftEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge )
.def("nodesRightEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge )
+ .def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge )
+ .def("nodesTopOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge )
+ .def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge )
+ .def("nodesRightOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge )
.def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner )
.def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner)
.def("nodesTopLeftCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner )
.def("nodesTopRightCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner )
.def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner )
.def("nodesLeftTopCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner )
.def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner)
.def("nodesRightTopCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner )
.def("nodesPeriodic" , &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Quad4::Regular::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Quad4::Regular::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic )
.def("__repr__",
[](const GooseFEM::Mesh::Quad4::Regular &a){ return ""; }
);
// -------------------------------------------------------------------------------------------------
py::class_(mMeshQuad4, "FineLayer")
.def(
- py::init(),
+ py::init(),
"FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical direction",
py::arg("nx"),
py::arg("ny"),
py::arg("h")=1.,
- py::arg("nfine")=1,
- py::arg("nskip")=0
+ py::arg("nfine")=1
)
.def("shape" , &GooseFEM::Mesh::Quad4::FineLayer::shape )
.def("coor" , &GooseFEM::Mesh::Quad4::FineLayer::coor )
.def("conn" , &GooseFEM::Mesh::Quad4::FineLayer::conn )
.def("nelem" , &GooseFEM::Mesh::Quad4::FineLayer::nelem )
.def("nnode" , &GooseFEM::Mesh::Quad4::FineLayer::nnode )
.def("nne" , &GooseFEM::Mesh::Quad4::FineLayer::nne )
.def("ndim" , &GooseFEM::Mesh::Quad4::FineLayer::ndim )
.def("elementsMiddleLayer" , &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer )
.def("nodesBottomEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge )
.def("nodesTopEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge )
.def("nodesLeftEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge )
.def("nodesRightEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge )
+ .def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge )
+ .def("nodesTopOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge )
+ .def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge )
+ .def("nodesRightOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge )
.def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner )
.def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner)
.def("nodesTopLeftCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner )
.def("nodesTopRightCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner )
.def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner )
.def("nodesLeftTopCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner )
.def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner)
.def("nodesRightTopCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner )
.def("nodesPeriodic" , &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Quad4::FineLayer::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic )
.def("__repr__",
[](const GooseFEM::Mesh::Quad4::FineLayer &a){ return ""; }
);
// ===================================== GooseFEM/MeshTri3.h ======================================
py::class_(mMeshTri3, "Regular")
.def(
py::init(),
"Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'",
py::arg("nx"),
py::arg("ny"),
py::arg("h")=1.
)
.def("coor" , &GooseFEM::Mesh::Tri3::Regular::coor )
.def("conn" , &GooseFEM::Mesh::Tri3::Regular::conn )
.def("nelem" , &GooseFEM::Mesh::Tri3::Regular::nelem )
.def("nnode" , &GooseFEM::Mesh::Tri3::Regular::nnode )
.def("nne" , &GooseFEM::Mesh::Tri3::Regular::nne )
.def("ndim" , &GooseFEM::Mesh::Tri3::Regular::ndim )
.def("nodesBottomEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomEdge )
.def("nodesTopEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesTopEdge )
.def("nodesLeftEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftEdge )
.def("nodesRightEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesRightEdge )
+ .def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomOpenEdge )
+ .def("nodesTopOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesTopOpenEdge )
+ .def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftOpenEdge )
+ .def("nodesRightOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesRightOpenEdge )
.def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomLeftCorner )
.def("nodesBottomRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomRightCorner)
.def("nodesTopLeftCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesTopLeftCorner )
.def("nodesTopRightCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesTopRightCorner )
.def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftBottomCorner )
.def("nodesLeftTopCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftTopCorner )
.def("nodesRightBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightBottomCorner)
.def("nodesRightTopCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesRightTopCorner )
.def("nodesPeriodic" , &GooseFEM::Mesh::Tri3::Regular::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Tri3::Regular::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Tri3::Regular::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Tri3::Regular::dofsPeriodic )
.def("__repr__",
[](const GooseFEM::Mesh::Tri3::Regular &a){ return ""; }
);
// -------------------------------------------------------------------------------------------------
mMeshTri3.def("getOrientation",&GooseFEM::Mesh::Tri3::getOrientation,
"Get the orientation of each element",
py::arg("coor"),
py::arg("conn")
);
// -------------------------------------------------------------------------------------------------
mMeshTri3.def("retriangulate",&GooseFEM::Mesh::Tri3::retriangulate,
"Re-triangulate existing mesh",
py::arg("coor"),
py::arg("conn"),
py::arg("orientation")=-1
);
// =================================================================================================
}
diff --git a/src/GooseFEM/MeshHex8.cpp b/src/GooseFEM/MeshHex8.cpp
index f27a375..9d181f1 100644
--- a/src/GooseFEM/MeshHex8.cpp
+++ b/src/GooseFEM/MeshHex8.cpp
@@ -1,2608 +1,2625 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHHEX8_CPP
#define GOOSEFEM_MESHHEX8_CPP
// -------------------------------------------------------------------------------------------------
#include "MeshHex8.h"
// ===================================== GooseFEM::Mesh::Hex8 ======================================
namespace GooseFEM {
namespace Mesh {
namespace Hex8 {
// ===================================== CLASS - REGULAR MESH ======================================
// ------------------------------------------ constructor ------------------------------------------
inline Regular::Regular(size_t nelx, size_t nely, size_t nelz, double h):
m_h(h), m_nelx(nelx), m_nely(nely), m_nelz(nelz)
{
assert( m_nelx >= 1 );
assert( m_nely >= 1 );
assert( m_nelz >= 1 );
m_nnode = (m_nelx+1) * (m_nely+1) * (m_nelz+1);
m_nelem = m_nelx * m_nely * m_nelz ;
}
// -------------------------------------- number of elements ---------------------------------------
inline size_t Regular::nelem()
{
return m_nelem;
}
// ---------------------------------------- number of nodes ----------------------------------------
inline size_t Regular::nnode()
{
return m_nnode;
}
// ---------------------------------- number of nodes per element ----------------------------------
inline size_t Regular::nne()
{
return m_nne;
}
// ------------------------------------- number of dimensions --------------------------------------
inline size_t Regular::ndim()
{
return m_ndim;
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline MatD Regular::coor()
{
MatD out(m_nnode,m_ndim);
ColD x = ColD::LinSpaced(m_nelx+1, 0.0, m_h*static_cast(m_nelx));
ColD y = ColD::LinSpaced(m_nely+1, 0.0, m_h*static_cast(m_nely));
ColD z = ColD::LinSpaced(m_nelz+1, 0.0, m_h*static_cast(m_nelz));
size_t inode = 0;
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) {
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) {
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy);
out(inode,2) = z(iz);
++inode;
}
}
}
return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline MatS Regular::conn()
{
MatS out(m_nelem,m_nne);
size_t ielem = 0;
for ( size_t iz = 0 ; iz < m_nelz ; ++iz ) {
for ( size_t iy = 0 ; iy < m_nely ; ++iy ) {
for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) {
- out(ielem,0) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+0);
- out(ielem,1) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+1);
- out(ielem,3) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+0);
- out(ielem,2) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+1);
- out(ielem,4) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+0);
- out(ielem,5) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+1);
- out(ielem,7) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+0);
- out(ielem,6) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+1);
+ out(ielem,0) = (iy )*(m_nelx+1) + (ix ) + (iz )*(m_nely+1)*(m_nelx+1);
+ out(ielem,1) = (iy )*(m_nelx+1) + (ix+1) + (iz )*(m_nely+1)*(m_nelx+1);
+ out(ielem,3) = (iy+1)*(m_nelx+1) + (ix ) + (iz )*(m_nely+1)*(m_nelx+1);
+ out(ielem,2) = (iy+1)*(m_nelx+1) + (ix+1) + (iz )*(m_nely+1)*(m_nelx+1);
+ out(ielem,4) = (iy )*(m_nelx+1) + (ix ) + (iz+1)*(m_nely+1)*(m_nelx+1);
+ out(ielem,5) = (iy )*(m_nelx+1) + (ix+1) + (iz+1)*(m_nely+1)*(m_nelx+1);
+ out(ielem,7) = (iy+1)*(m_nelx+1) + (ix ) + (iz+1)*(m_nely+1)*(m_nelx+1);
+ out(ielem,6) = (iy+1)*(m_nelx+1) + (ix+1) + (iz+1)*(m_nely+1)*(m_nelx+1);
++ielem;
}
}
}
return out;
}
// ------------------------------ node-numbers along the front plane -------------------------------
inline ColS Regular::nodesFront()
{
ColS out((m_nelx+1)*(m_nely+1));
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iy*(m_nelx+1)+ix) = iy*(m_nelx+1) + ix;
return out;
}
// ------------------------------- node-numbers along the back plane -------------------------------
inline ColS Regular::nodesBack()
{
ColS out((m_nelx+1)*(m_nely+1));
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iy*(m_nelx+1)+ix) = iy*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// ------------------------------- node-numbers along the left plane -------------------------------
inline ColS Regular::nodesLeft()
{
ColS out((m_nely+1)*(m_nelz+1));
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iz*(m_nely+1)+iy) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------ node-numbers along the right plane -------------------------------
inline ColS Regular::nodesRight()
{
ColS out((m_nely+1)*(m_nelz+1));
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iz*(m_nely+1)+iy) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ------------------------------ node-numbers along the bottom plane ------------------------------
inline ColS Regular::nodesBottom()
{
ColS out((m_nelx+1)*(m_nelz+1));
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iz*(m_nelx+1)+ix) = ix + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------- node-numbers along the top plane --------------------------------
inline ColS Regular::nodesTop()
{
ColS out((m_nelx+1)*(m_nelz+1));
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iz*(m_nelx+1)+ix) = ix + m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------- node-numbers along the front face -------------------------------
inline ColS Regular::nodesFrontFace()
{
ColS out((m_nelx-1)*(m_nely-1));
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out((iy-1)*(m_nelx-1)+(ix-1)) = iy*(m_nelx+1) + ix;
return out;
}
// ------------------------------- node-numbers along the back face --------------------------------
inline ColS Regular::nodesBackFace()
{
ColS out((m_nelx-1)*(m_nely-1));
for ( size_t iy = 1 ; iy < m_nely ; ++iy ) {
for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) {
out((iy-1)*(m_nelx-1)+(ix-1)) = iy*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
}
}
return out;
}
// ------------------------------- node-numbers along the left face --------------------------------
inline ColS Regular::nodesLeftFace()
{
ColS out((m_nely-1)*(m_nelz-1));
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out((iz-1)*(m_nely-1)+(iy-1)) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------- node-numbers along the right face -------------------------------
inline ColS Regular::nodesRightFace()
{
ColS out((m_nely-1)*(m_nelz-1));
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out((iz-1)*(m_nely-1)+(iy-1)) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ------------------------------ node-numbers along the bottom face -------------------------------
inline ColS Regular::nodesBottomFace()
{
ColS out((m_nelx-1)*(m_nelz-1));
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out((iz-1)*(m_nelx-1)+(ix-1)) = ix + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// -------------------------------- node-numbers along the top face --------------------------------
inline ColS Regular::nodesTopFace()
{
ColS out((m_nelx-1)*(m_nelz-1));
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out((iz-1)*(m_nelx-1)+(ix-1)) = ix + m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// --------------------------- node-numbers along the front-bottom edge ----------------------------
inline ColS Regular::nodesFrontBottomEdge()
{
ColS out(m_nelx+1);
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = ix;
return out;
}
// ----------------------------- node-numbers along the front-top edge -----------------------------
inline ColS Regular::nodesFrontTopEdge()
{
ColS out(m_nelx+1);
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
- out(ix) = m_nely*(m_nelx+1) + ix;
+ out(ix) = ix + m_nely*(m_nelx+1);
return out;
}
// ---------------------------- node-numbers along the front-left edge -----------------------------
inline ColS Regular::nodesFrontLeftEdge()
{
ColS out(m_nely+1);
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1);
return out;
}
// ---------------------------- node-numbers along the front-right edge ----------------------------
inline ColS Regular::nodesFrontRightEdge()
{
ColS out(m_nely+1);
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1) + m_nelx;
return out;
}
// ---------------------------- node-numbers along the back-bottom edge ----------------------------
inline ColS Regular::nodesBackBottomEdge()
{
ColS out(m_nelx+1);
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// ----------------------------- node-numbers along the back-top edge ------------------------------
inline ColS Regular::nodesBackTopEdge()
{
ColS out(m_nelx+1);
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = m_nely*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// ----------------------------- node-numbers along the back-left edge -----------------------------
inline ColS Regular::nodesBackLeftEdge()
{
ColS out(m_nely+1);
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1);
return out;
}
// ---------------------------- node-numbers along the back-right edge -----------------------------
inline ColS Regular::nodesBackRightEdge()
{
ColS out(m_nely+1);
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ---------------------------- node-numbers along the bottom-left edge ----------------------------
inline ColS Regular::nodesBottomLeftEdge()
{
ColS out(m_nelz+1);
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = iz*(m_nelx+1)*(m_nely+1);
return out;
}
// --------------------------- node-numbers along the bottom-right edge ----------------------------
inline ColS Regular::nodesBottomRightEdge()
{
ColS out(m_nelz+1);
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ----------------------------- node-numbers along the top-left edge ------------------------------
inline ColS Regular::nodesTopLeftEdge()
{
ColS out(m_nelz+1);
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ----------------------------- node-numbers along the top-right edge -----------------------------
inline ColS Regular::nodesTopRightEdge()
{
ColS out(m_nelz+1);
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline ColS Regular::nodesBottomFrontEdge() { return nodesFrontBottomEdge(); }
inline ColS Regular::nodesBottomBackEdge() { return nodesBackBottomEdge(); }
inline ColS Regular::nodesTopFrontEdge() { return nodesFrontTopEdge(); }
inline ColS Regular::nodesTopBackEdge() { return nodesBackTopEdge(); }
inline ColS Regular::nodesLeftBottomEdge() { return nodesBottomLeftEdge(); }
inline ColS Regular::nodesLeftFrontEdge() { return nodesFrontLeftEdge(); }
inline ColS Regular::nodesLeftBackEdge() { return nodesBackLeftEdge(); }
inline ColS Regular::nodesLeftTopEdge() { return nodesTopLeftEdge(); }
inline ColS Regular::nodesRightBottomEdge() { return nodesBottomRightEdge(); }
inline ColS Regular::nodesRightTopEdge() { return nodesTopRightEdge(); }
inline ColS Regular::nodesRightFrontEdge() { return nodesFrontRightEdge(); }
inline ColS Regular::nodesRightBackEdge() { return nodesBackRightEdge(); }
// ------------------- node-numbers along the front-bottom edge, without corners -------------------
inline ColS Regular::nodesFrontBottomOpenEdge()
{
ColS out(m_nelx-1);
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = ix;
return out;
}
// -------------------- node-numbers along the front-top edge, without corners ---------------------
inline ColS Regular::nodesFrontTopOpenEdge()
{
ColS out(m_nelx-1);
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
- out(ix-1) = m_nely*(m_nelx+1) + ix;
+ out(ix-1) = ix + m_nely*(m_nelx+1);
return out;
}
// -------------------- node-numbers along the front-left edge, without corners --------------------
inline ColS Regular::nodesFrontLeftOpenEdge()
{
ColS out(m_nely-1);
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1);
return out;
}
// ------------------- node-numbers along the front-right edge, without corners --------------------
inline ColS Regular::nodesFrontRightOpenEdge()
{
ColS out(m_nely-1);
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelx;
return out;
}
// ------------------- node-numbers along the back-bottom edge, without corners --------------------
inline ColS Regular::nodesBackBottomOpenEdge()
{
ColS out(m_nelx-1);
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// --------------------- node-numbers along the back-top edge, without corners ---------------------
inline ColS Regular::nodesBackTopOpenEdge()
{
ColS out(m_nelx-1);
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = m_nely*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// -------------------- node-numbers along the back-left edge, without corners ---------------------
inline ColS Regular::nodesBackLeftOpenEdge()
{
ColS out(m_nely-1);
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1);
return out;
}
// -------------------- node-numbers along the back-right edge, without corners --------------------
inline ColS Regular::nodesBackRightOpenEdge()
{
ColS out(m_nely-1);
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ------------------- node-numbers along the bottom-left edge, without corners --------------------
inline ColS Regular::nodesBottomLeftOpenEdge()
{
ColS out(m_nelz-1);
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------- node-numbers along the bottom-right edge, without corners -------------------
inline ColS Regular::nodesBottomRightOpenEdge()
{
ColS out(m_nelz-1);
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// --------------------- node-numbers along the top-left edge, without corners ---------------------
inline ColS Regular::nodesTopLeftOpenEdge()
{
ColS out(m_nelz-1);
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// -------------------- node-numbers along the top-right edge, without corners ---------------------
inline ColS Regular::nodesTopRightOpenEdge()
{
ColS out(m_nelz-1);
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline ColS Regular::nodesBottomFrontOpenEdge() { return nodesFrontBottomOpenEdge(); }
inline ColS Regular::nodesBottomBackOpenEdge() { return nodesBackBottomOpenEdge(); }
inline ColS Regular::nodesTopFrontOpenEdge() { return nodesFrontTopOpenEdge(); }
inline ColS Regular::nodesTopBackOpenEdge() { return nodesBackTopOpenEdge(); }
inline ColS Regular::nodesLeftBottomOpenEdge() { return nodesBottomLeftOpenEdge(); }
inline ColS Regular::nodesLeftFrontOpenEdge() { return nodesFrontLeftOpenEdge(); }
inline ColS Regular::nodesLeftBackOpenEdge() { return nodesBackLeftOpenEdge(); }
inline ColS Regular::nodesLeftTopOpenEdge() { return nodesTopLeftOpenEdge(); }
inline ColS Regular::nodesRightBottomOpenEdge() { return nodesBottomRightOpenEdge(); }
inline ColS Regular::nodesRightTopOpenEdge() { return nodesTopRightOpenEdge(); }
inline ColS Regular::nodesRightFrontOpenEdge() { return nodesFrontRightOpenEdge(); }
inline ColS Regular::nodesRightBackOpenEdge() { return nodesBackRightOpenEdge(); }
// -------------------------- node-number of the front-bottom-left corner --------------------------
inline size_t Regular::nodesFrontBottomLeftCorner()
{
return 0;
}
// ------------------------- node-number of the front-bottom-right corner --------------------------
inline size_t Regular::nodesFrontBottomRightCorner()
{
return m_nelx;
}
// --------------------------- node-number of the front-top-left corner ----------------------------
inline size_t Regular::nodesFrontTopLeftCorner()
{
return m_nely*(m_nelx+1);
}
// --------------------------- node-number of the front-top-right corner ---------------------------
inline size_t Regular::nodesFrontTopRightCorner()
{
return m_nely*(m_nelx+1) + m_nelx;
}
// -------------------------- node-number of the back-bottom-left corner ---------------------------
inline size_t Regular::nodesBackBottomLeftCorner()
{
return m_nelz*(m_nely+1)*(m_nelx+1);
}
// -------------------------- node-number of the back-bottom-right corner --------------------------
inline size_t Regular::nodesBackBottomRightCorner()
{
return m_nelx + m_nelz*(m_nely+1)*(m_nelx+1);
}
// ---------------------------- node-number of the back-top-left corner ----------------------------
inline size_t Regular::nodesBackTopLeftCorner()
{
return m_nely*(m_nelx+1) + m_nelz*(m_nely+1)*(m_nelx+1);
}
// --------------------------- node-number of the back-top-right corner ----------------------------
inline size_t Regular::nodesBackTopRightCorner()
{
return m_nely*(m_nelx+1) + m_nelx + m_nelz*(m_nely+1)*(m_nelx+1);
}
// -------------------------------------------- aliases --------------------------------------------
inline size_t Regular::nodesFrontLeftBottomCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesBottomFrontLeftCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesBottomLeftFrontCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesLeftFrontBottomCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesLeftBottomFrontCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesFrontRightBottomCorner() { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesBottomFrontRightCorner() { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesBottomRightFrontCorner() { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesRightFrontBottomCorner() { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesRightBottomFrontCorner() { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesFrontLeftTopCorner() { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesTopFrontLeftCorner() { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesTopLeftFrontCorner() { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesLeftFrontTopCorner() { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesLeftTopFrontCorner() { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesFrontRightTopCorner() { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesTopFrontRightCorner() { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesTopRightFrontCorner() { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesRightFrontTopCorner() { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesRightTopFrontCorner() { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesBackLeftBottomCorner() { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesBottomBackLeftCorner() { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesBottomLeftBackCorner() { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesLeftBackBottomCorner() { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesLeftBottomBackCorner() { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesBackRightBottomCorner() { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesBottomBackRightCorner() { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesBottomRightBackCorner() { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesRightBackBottomCorner() { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesRightBottomBackCorner() { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesBackLeftTopCorner() { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesTopBackLeftCorner() { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesTopLeftBackCorner() { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesLeftBackTopCorner() { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesLeftTopBackCorner() { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesBackRightTopCorner() { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesTopBackRightCorner() { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesTopRightBackCorner() { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesRightBackTopCorner() { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesRightTopBackCorner() { return nodesBackTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline MatS Regular::nodesPeriodic()
{
// faces
ColS fro = nodesFrontFace();
ColS bck = nodesBackFace();
ColS lft = nodesLeftFace();
ColS rgt = nodesRightFace();
ColS bot = nodesBottomFace();
ColS top = nodesTopFace();
// edges
ColS froBot = nodesFrontBottomOpenEdge();
ColS froTop = nodesFrontTopOpenEdge();
ColS froLft = nodesFrontLeftOpenEdge();
ColS froRgt = nodesFrontRightOpenEdge();
ColS bckBot = nodesBackBottomOpenEdge();
ColS bckTop = nodesBackTopOpenEdge();
ColS bckLft = nodesBackLeftOpenEdge();
ColS bckRgt = nodesBackRightOpenEdge();
ColS botLft = nodesBottomLeftOpenEdge();
ColS botRgt = nodesBottomRightOpenEdge();
ColS topLft = nodesTopLeftOpenEdge();
ColS topRgt = nodesTopRightOpenEdge();
// allocate nodal ties
// - number of tying per category
size_t tface = fro.size() + lft.size() + bot.size();
size_t tedge = 3*froBot.size() + 3*froLft.size() + 3*botLft.size();
size_t tnode = 7;
// - allocate
MatS out(tface+tedge+tnode, 2);
// counter
size_t i = 0;
// tie all corners to one corner
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopLeftCorner(); ++i;
// tie all corresponding edges to each other (exclude corners)
for ( auto j = 0 ; j(nodePer.rows());
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for ( size_t i = 0 ; i < nper ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
out(nodePer(i,1),j) = out(nodePer(i,0),j);
// renumber "out" to be sequential
return GooseFEM::Mesh::renumber(out);
}
// ==================================== CLASS - FINELAYER MESH =====================================
// ------------------------------------------ constructor ------------------------------------------
inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine):
m_h(h), m_nelx(nelx), m_nelz(nelz)
{
// basic assumptions
assert( nelx >= 1 );
assert( nely >= 1 );
assert( nelz >= 1 );
// store basic info
m_Lx = m_h * static_cast(nelx);
m_Lz = m_h * static_cast(nelz);
// compute element size in y-direction (use symmetry, compute upper half)
// ----------------------------------------------------------------------
// temporary variables
size_t nmin, ntot;
ColS nhx(nely), nhy(nely), nhz(nely);
ColI refine(nely);
// minimum height in y-direction (half of the height because of symmetry)
if ( nely % 2 == 0 ) nmin = nely /2;
else nmin = (nely +1)/2;
// minimum number of fine layers in y-direction (minimum 1, middle layer part of this half)
if ( nfine % 2 == 0 ) nfine = nfine /2 + 1;
else nfine = (nfine+1)/2;
if ( nfine < 1 ) nfine = 1;
if ( nfine > nmin ) nfine = nmin;
// initialize to state with only fine elements
nhx .setOnes();
nhy .setOnes();
nhz .setOnes();
refine.setConstant(-1);
// loop over element layers in y-direction, try to coarsen using these rules:
// (1) element size in y-direction <= distance to origin in y-direction
// (2) element size in x-(z-)direction should fit the total number of elements in x-(z-)direction
// (3) a certain number of layers have the minimum size "1" (are fine)
for ( size_t iy = nfine ; ; )
{
// initialize current size in y-direction
if ( iy == nfine ) ntot = nfine;
// check to stop
if ( iy >= nely or ntot >= nmin ) { nely = iy; break; }
// rules (1,2) satisfied: coarsen in x-direction (and z-direction)
if ( 3*nhy(iy) <= ntot and nelx%(3*nhx(iy)) == 0 )
{
// - process refinement in x-direction
refine (iy ) = 0;
nhy (iy ) *= 2;
nhy.segment(iy+1,nely-iy-1) *= 3;
nhx.segment(iy ,nely-iy ) *= 3;
// - rule (2) satisfied: coarsen next element layer in z-direction
if ( iy+1 < nely and ntot+nhy(iy) < nmin )
{
if ( nelz%(3*nhz(iy+1)) == 0 )
{
// - update the number of elements in y-direction
ntot += nhy(iy);
// - proceed to next element layer in y-direction
++iy;
// - process refinement in z-direction
refine (iy ) = 2;
nhy (iy ) = nhy(iy-1);
nhz.segment(iy,nely-iy) *= 3;
}
}
}
// rules (1,2) satisfied: coarse in z-direction
else if ( 3*nhy(iy) <= ntot and nelz%(3*nhz(iy)) == 0 )
{
// - process refinement in z-direction
refine (iy ) = 2;
nhy (iy ) *= 2;
nhy.segment(iy+1,nely-iy-1) *= 3;
nhz.segment(iy ,nely-iy ) *= 3;
}
// update the number of elements in y-direction
ntot += nhy(iy);
// proceed to next element layer in y-direction
++iy;
// check to stop
if ( iy >= nely or ntot >= nmin ) { nely = iy; break; }
}
// symmetrize, compute full information
// ------------------------------------
// allocate mesh constructor parameters
m_nhx .conservativeResize(nely*2-1);
m_nhy .conservativeResize(nely*2-1);
m_nhz .conservativeResize(nely*2-1);
m_refine .conservativeResize(nely*2-1);
m_nelx .conservativeResize(nely*2-1);
m_nelz .conservativeResize(nely*2-1);
m_nnd .conservativeResize(nely*2 );
m_startElem.conservativeResize(nely*2-1);
m_startNode.conservativeResize(nely*2 );
// fill
// - lower half
for ( size_t iy = 0 ; iy < nely ; ++iy )
{
m_nhx (iy) = nhx (nely-iy-1);
m_nhy (iy) = nhy (nely-iy-1);
m_nhz (iy) = nhz (nely-iy-1);
m_refine(iy) = refine(nely-iy-1);
}
// - upper half
for ( size_t iy = 0 ; iy < nely-1 ; ++iy )
{
m_nhx (iy+nely) = nhx (iy+1);
m_nhy (iy+nely) = nhy (iy+1);
m_nhz (iy+nely) = nhz (iy+1);
m_refine(iy+nely) = refine(iy+1);
}
// update size
nely = m_nhx.size();
// compute the number of elements per element layer in y-direction
for ( size_t iy = 0 ; iy < nely ; ++iy )
{
m_nelx(iy) = nelx / m_nhx(iy);
m_nelz(iy) = nelz / m_nhz(iy);
}
// compute the number of nodes per node layer in y-direction
// - bottom half
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
m_nnd(iy) = (m_nelx(iy)+1) * (m_nelz(iy)+1);
// - top half
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
m_nnd(iy+1) = (m_nelx(iy)+1) * (m_nelz(iy)+1);
// compute mesh dimensions
// -----------------------
// initialize
m_nnode = 0;
m_nelem = 0;
m_startNode(0) = 0;
// loop over element layers (bottom -> middle, elements become finer)
for ( size_t i = 0 ; i < (nely-1)/2 ; ++i )
{
// - store the first element of the layer
m_startElem(i) = m_nelem;
// - add the nodes of this layer
if ( m_refine(i) == 0 ) { m_nnode += (3*m_nelx(i)+1) * ( m_nelz(i)+1); }
else if ( m_refine(i) == 2 ) { m_nnode += ( m_nelx(i)+1) * (3*m_nelz(i)+1); }
else { m_nnode += ( m_nelx(i)+1) * ( m_nelz(i)+1); }
// - add the elements of this layer
if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ) * ( m_nelz(i) ); }
else if ( m_refine(i) == 2 ) { m_nelem += ( m_nelx(i) ) * (4*m_nelz(i) ); }
else { m_nelem += ( m_nelx(i) ) * ( m_nelz(i) ); }
// - store the starting node of the next layer
m_startNode(i+1) = m_nnode;
}
// loop over element layers (middle -> top, elements become coarser)
for ( size_t i = (nely-1)/2 ; i < nely ; ++i )
{
// - store the first element of the layer
m_startElem(i) = m_nelem;
// - add the nodes of this layer
if ( m_refine(i) == 0 ) { m_nnode += (5*m_nelx(i)+1) * ( m_nelz(i)+1); }
else if ( m_refine(i) == 2 ) { m_nnode += ( m_nelx(i)+1) * (5*m_nelz(i)+1); }
else { m_nnode += ( m_nelx(i)+1) * ( m_nelz(i)+1); }
// - add the elements of this layer
if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ) * ( m_nelz(i) ); }
else if ( m_refine(i) == 2 ) { m_nelem += ( m_nelx(i) ) * (4*m_nelz(i) ); }
else { m_nelem += ( m_nelx(i) ) * ( m_nelz(i) ); }
// - store the starting node of the next layer
m_startNode(i+1) = m_nnode;
}
// - add the top row of nodes
m_nnode += (m_nelx(nely-1)+1) * (m_nelz(nely-1)+1);
}
// -------------------------------------- number of elements ---------------------------------------
inline size_t FineLayer::nelem()
{
return m_nelem;
}
// ---------------------------------------- number of nodes ----------------------------------------
inline size_t FineLayer::nnode()
{
return m_nnode;
}
// ---------------------------------- number of nodes per element ----------------------------------
inline size_t FineLayer::nne()
{
return m_nne;
}
// ------------------------------------- number of dimensions --------------------------------------
inline size_t FineLayer::ndim()
{
return m_ndim;
}
// ---------------------------- actual number of nodes in one direction ----------------------------
inline size_t FineLayer::shape(size_t i)
{
assert( i >= 0 and i <= 2 );
if ( i == 0 ) return m_nelx.maxCoeff();
else if ( i == 2 ) return m_nelz.maxCoeff();
else return m_nhy .sum();
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline MatD FineLayer::coor()
{
// allocate output
MatD out(m_nnode, m_ndim);
// current node, number of element layers
size_t inode = 0;
size_t nely = static_cast(m_nhy.size());
// y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
// - allocate
ColD y(nely+1);
// - initialize
y(0) = 0.0;
// - compute
for ( size_t iy = 1 ; iy < nely+1 ; ++iy )
y(iy) = y(iy-1) + m_nhy(iy-1) * m_h;
// loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes
// --------------------------------------------------------------------------------------------
for ( size_t iy = 0 ; ; ++iy )
{
// get positions along the x- and z-axis
ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx);
ColD z = ColD::LinSpaced(m_nelz(iy)+1, 0.0, m_Lz);
// add nodes of the bottom layer of this element
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy);
out(inode,2) = z(iz);
++inode;
}
}
// stop at middle layer
if ( iy == (nely-1)/2 )
break;
// add extra nodes of the intermediate layer, for refinement in x-direction
if ( m_refine(iy) == 0 )
{
// - get position offset in x- and y-direction
double dx = m_h * static_cast(m_nhx(iy)/3);
double dy = m_h * static_cast(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
for ( size_t j = 0 ; j < 2 ; ++j ) {
out(inode,0) = x(ix) + dx * static_cast(j+1);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz);
++inode;
}
}
}
}
// add extra nodes of the intermediate layer, for refinement in z-direction
else if ( m_refine(iy) == 2 )
{
// - get position offset in y- and z-direction
double dz = m_h * static_cast(m_nhz(iy)/3);
double dy = m_h * static_cast(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t j = 0 ; j < 2 ; ++j ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz) + dz * static_cast(j+1);
++inode;
}
}
}
}
}
// loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes
// --------------------------------------------------------------------------------------
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// get positions along the x- and z-axis
ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx);
ColD z = ColD::LinSpaced(m_nelz(iy)+1, 0.0, m_Lz);
// add extra nodes of the intermediate layer, for refinement in x-direction
if ( m_refine(iy) == 0 )
{
// - get position offset in x- and y-direction
double dx = m_h * static_cast(m_nhx(iy)/3);
double dy = m_h * static_cast(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
for ( size_t j = 0 ; j < 2 ; ++j ) {
out(inode,0) = x(ix) + dx * static_cast(j+1);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz);
++inode;
}
}
}
}
// add extra nodes of the intermediate layer, for refinement in z-direction
else if ( m_refine(iy) == 2 )
{
// - get position offset in y- and z-direction
double dz = m_h * static_cast(m_nhz(iy)/3);
double dy = m_h * static_cast(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t j = 0 ; j < 2 ; ++j ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz) + dz * static_cast(j+1);
++inode;
}
}
}
}
// add nodes of the top layer of this element
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix );
out(inode,1) = y(iy+1);
out(inode,2) = z(iz );
++inode;
}
}
}
return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline MatS FineLayer::conn()
{
// allocate output
MatS out(m_nelem, m_nne);
// current element, number of element layers, starting nodes of each node layer
size_t ielem = 0;
size_t nely = static_cast(m_nhy.size());
size_t bot,mid,top;
// loop over all element layers
for ( size_t iy = 0 ; iy < nely ; ++iy )
{
// - get: starting nodes of bottom(, middle) and top layer
bot = m_startNode(iy );
mid = m_startNode(iy ) + m_nnd(iy);
top = m_startNode(iy+1);
// - define connectivity: no coarsening/refinement
if ( m_refine(iy) == -1 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
out(ielem,0) = bot + (ix ) + (iz ) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (iz ) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (iz+1) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (iz+1) * (m_nelx(iy)+1);
ielem++;
}
}
}
// - define connectivity: refinement along the x-direction (below the middle layer)
else if ( m_refine(iy) == 0 and iy <= (nely-1)/2 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- bottom element
out(ielem,0) = bot + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,1) = bot + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,2) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,3) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,5) = bot + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,6) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,7) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- top-right element
out(ielem,0) = bot + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,1) = top + (3*ix+3) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = top + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,3) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,5) = top + (3*ix+3) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = top + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,7) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- top-center element
out(ielem,0) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,1) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,2) = top + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,3) = top + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,4) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,5) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,6) = top + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,7) = top + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
ielem++;
// -- top-left element
out(ielem,0) = bot + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,1) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,2) = top + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,3) = top + (3*ix ) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,4) = bot + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,5) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,6) = top + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,7) = top + (3*ix ) + (iz+1) * (3*m_nelx(iy)+1);
ielem++;
}
}
}
// - define connectivity: coarsening along the x-direction (above the middle layer)
else if ( m_refine(iy) == 0 and iy > (nely-1)/2 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- lower-left element
out(ielem,0) = bot + (3*ix ) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,1) = bot + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,3) = top + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,4) = bot + (3*ix ) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,5) = bot + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,7) = top + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
ielem++;
// -- lower-center element
out(ielem,0) = bot + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,1) = bot + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,3) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,5) = bot + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,7) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- lower-right element
out(ielem,0) = bot + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,1) = bot + (3*ix+3) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = top + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,3) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,5) = bot + (3*ix+3) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = top + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,7) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- upper element
out(ielem,0) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,1) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,2) = top + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,3) = top + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,4) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,5) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,6) = top + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,7) = top + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
ielem++;
}
}
}
// - define connectivity: refinement along the z-direction (below the middle layer)
else if ( m_refine(iy) == 2 and iy <= (nely-1)/2 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- bottom element
out(ielem,0) = bot + (ix ) + iz * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + iz * (m_nelx(iy)+1);
out(ielem,2) = bot + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,3) = bot + (ix ) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + 2*iz * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + 2*iz * (m_nelx(iy)+1);
out(ielem,6) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,7) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
ielem++;
// -- top-back element
out(ielem,0) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,1) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (3*iz+3) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (3*iz+3) * (m_nelx(iy)+1);
ielem++;
// -- top-center element
out(ielem,0) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,1) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
ielem++;
// -- top-front element
out(ielem,0) = bot + (ix ) + ( iz ) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + ( iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
ielem++;
}
}
}
// - define connectivity: coarsening along the z-direction (above the middle layer)
else if ( m_refine(iy) == 2 and iy > (nely-1)/2 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- bottom-front element
out(ielem,0) = bot + (ix ) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + ( iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + ( iz ) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,7) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
ielem++;
// -- bottom-center element
out(ielem,0) = bot + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,2) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,3) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,6) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,7) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
ielem++;
// -- bottom-back element
out(ielem,0) = bot + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,2) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,3) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (3*iz+3) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (3*iz+3) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + ( iz+1) * (m_nelx(iy)+1);
ielem++;
// -- top element
out(ielem,0) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,1) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + ( iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + ( iz ) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + ( iz+1) * (m_nelx(iy)+1);
ielem++;
}
}
}
}
return out;
}
+// ------------------------------ element numbers of the middle layer ------------------------------
+
+inline ColS FineLayer::elementsMiddleLayer()
+{
+ // number of element layers in y-direction, the index of the middle layer
+ size_t nely = static_cast(m_nhy.size());
+ size_t iy = (nely-1)/2;
+
+ ColS out(m_nelx(iy)*m_nelz(iy));
+
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix )
+ for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz )
+ out(ix+iz*m_nelx(iy)) = m_startElem(iy) + ix + iz*m_nelx(iy);
+
+ return out;
+}
+
// ------------------------------ node-numbers along the front plane -------------------------------
inline ColS FineLayer::nodesFront()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy) + ix;
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy+1) + ix;
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the back plane -------------------------------
inline ColS FineLayer::nodesBack()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy+1) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the left plane -------------------------------
inline ColS FineLayer::nodesLeft()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1);
++j;
}
}
return out;
}
// ------------------------------ node-numbers along the right plane -------------------------------
inline ColS FineLayer::nodesRight()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
return out;
}
// ------------------------------ node-numbers along the bottom plane ------------------------------
inline ColS FineLayer::nodesBottom()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// allocate node list
ColS out(m_nnd(nely));
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) {
for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz ) {
out(j) = m_startNode(0) + ix + iz * (m_nelx(0)+1);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the top plane --------------------------------
inline ColS FineLayer::nodesTop()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// allocate node list
ColS out(m_nnd(nely));
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) {
for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz ) {
out(j) = m_startNode(nely) + ix + iz * (m_nelx(nely-1)+1);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the front face -------------------------------
inline ColS FineLayer::nodesFrontFace()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix;
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy+1) + ix;
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the back face --------------------------------
inline ColS FineLayer::nodesBackFace()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy+1) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the left face --------------------------------
inline ColS FineLayer::nodesLeftFace()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the right face -------------------------------
inline ColS FineLayer::nodesRightFace()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// allocate node-list
ColS out(n);
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
return out;
}
// ------------------------------ node-numbers along the bottom face -------------------------------
inline ColS FineLayer::nodesBottomFace()
{
// allocate node list
ColS out((m_nelx(0)-1)*(m_nelz(0)-1));
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) {
for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz ) {
out(j) = m_startNode(0) + ix + iz * (m_nelx(0)+1);
++j;
}
}
return out;
}
// -------------------------------- node-numbers along the top face --------------------------------
inline ColS FineLayer::nodesTopFace()
{
// number of element layers in y-direction
size_t nely = static_cast(m_nhy.size());
// allocate node list
ColS out((m_nelx(nely-1)-1)*(m_nelz(nely-1)-1));
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) {
for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz ) {
out(j) = m_startNode(nely) + ix + iz * (m_nelx(nely-1)+1);
++j;
}
}
return out;
}
// --------------------------- node-numbers along the front-bottom edge ----------------------------
inline ColS FineLayer::nodesFrontBottomEdge()
{
ColS out(m_nelx(0)+1);
for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix )
out(ix) = m_startNode(0) + ix;
return out;
}
// ----------------------------- node-numbers along the front-top edge -----------------------------
inline ColS FineLayer::nodesFrontTopEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelx(nely-1)+1);
for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix )
out(ix) = m_startNode(nely) + ix;
return out;
}
// ---------------------------- node-numbers along the front-left edge -----------------------------
inline ColS FineLayer::nodesFrontLeftEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely+1);
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
out(iy) = m_startNode(iy);
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
out(iy+1) = m_startNode(iy+1);
return out;
}
// ---------------------------- node-numbers along the front-right edge ----------------------------
inline ColS FineLayer::nodesFrontRightEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely+1);
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
out(iy) = m_startNode(iy) + m_nelx(iy);
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
out(iy+1) = m_startNode(iy+1) + m_nelx(iy);
return out;
}
// ---------------------------- node-numbers along the back-bottom edge ----------------------------
inline ColS FineLayer::nodesBackBottomEdge()
{
ColS out(m_nelx(0)+1);
for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix )
out(ix) = m_startNode(0) + ix + (m_nelx(0)+1)*(m_nelz(0));
return out;
}
// ----------------------------- node-numbers along the back-top edge ------------------------------
inline ColS FineLayer::nodesBackTopEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelx(nely-1)+1);
for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix )
out(ix) = m_startNode(nely) + ix + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
return out;
}
// ----------------------------- node-numbers along the back-left edge -----------------------------
inline ColS FineLayer::nodesBackLeftEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely+1);
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
out(iy) = m_startNode(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
out(iy+1) = m_startNode(iy+1) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// ---------------------------- node-numbers along the back-right edge -----------------------------
inline ColS FineLayer::nodesBackRightEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely+1);
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
out(iy) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
out(iy+1) = m_startNode(iy+1) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// ---------------------------- node-numbers along the bottom-left edge ----------------------------
inline ColS FineLayer::nodesBottomLeftEdge()
{
ColS out(m_nelz(0)+1);
for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz )
out(iz) = m_startNode(0) + iz * (m_nelx(0)+1);
return out;
}
// --------------------------- node-numbers along the bottom-right edge ----------------------------
inline ColS FineLayer::nodesBottomRightEdge()
{
ColS out(m_nelz(0)+1);
for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz )
out(iz) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0)+1);
return out;
}
// ----------------------------- node-numbers along the top-left edge ------------------------------
inline ColS FineLayer::nodesTopLeftEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelz(nely-1)+1);
for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz )
out(iz) = m_startNode(nely) + iz * (m_nelx(nely-1)+1);
return out;
}
// ----------------------------- node-numbers along the top-right edge -----------------------------
inline ColS FineLayer::nodesTopRightEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelz(nely-1)+1);
for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz )
out(iz) = m_startNode(nely) + m_nelx(nely-1) + iz * (m_nelx(nely-1)+1);
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline ColS FineLayer::nodesBottomFrontEdge() { return nodesFrontBottomEdge(); }
inline ColS FineLayer::nodesBottomBackEdge() { return nodesBackBottomEdge(); }
inline ColS FineLayer::nodesTopFrontEdge() { return nodesFrontTopEdge(); }
inline ColS FineLayer::nodesTopBackEdge() { return nodesBackTopEdge(); }
inline ColS FineLayer::nodesLeftBottomEdge() { return nodesBottomLeftEdge(); }
inline ColS FineLayer::nodesLeftFrontEdge() { return nodesFrontLeftEdge(); }
inline ColS FineLayer::nodesLeftBackEdge() { return nodesBackLeftEdge(); }
inline ColS FineLayer::nodesLeftTopEdge() { return nodesTopLeftEdge(); }
inline ColS FineLayer::nodesRightBottomEdge() { return nodesBottomRightEdge(); }
inline ColS FineLayer::nodesRightTopEdge() { return nodesTopRightEdge(); }
inline ColS FineLayer::nodesRightFrontEdge() { return nodesFrontRightEdge(); }
inline ColS FineLayer::nodesRightBackEdge() { return nodesBackRightEdge(); }
// ------------------- node-numbers along the front-bottom edge, without corners -------------------
inline ColS FineLayer::nodesFrontBottomOpenEdge()
{
ColS out(m_nelx(0)-1);
for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix )
out(ix-1) = m_startNode(0) + ix;
return out;
}
// -------------------- node-numbers along the front-top edge, without corners ---------------------
inline ColS FineLayer::nodesFrontTopOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelx(nely-1)-1);
for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix )
out(ix-1) = m_startNode(nely) + ix;
return out;
}
// -------------------- node-numbers along the front-left edge, without corners --------------------
inline ColS FineLayer::nodesFrontLeftOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely-1);
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
out(iy-1) = m_startNode(iy);
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
out(iy) = m_startNode(iy+1);
return out;
}
// ------------------- node-numbers along the front-right edge, without corners --------------------
inline ColS FineLayer::nodesFrontRightOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely-1);
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
out(iy-1) = m_startNode(iy) + m_nelx(iy);
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
out(iy) = m_startNode(iy+1) + m_nelx(iy);
return out;
}
// ------------------- node-numbers along the back-bottom edge, without corners --------------------
inline ColS FineLayer::nodesBackBottomOpenEdge()
{
ColS out(m_nelx(0)-1);
for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix )
out(ix-1) = m_startNode(0) + ix + (m_nelx(0)+1)*(m_nelz(0));
return out;
}
// --------------------- node-numbers along the back-top edge, without corners ---------------------
inline ColS FineLayer::nodesBackTopOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelx(nely-1)-1);
for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix )
out(ix-1) = m_startNode(nely) + ix + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
return out;
}
// -------------------- node-numbers along the back-left edge, without corners ---------------------
inline ColS FineLayer::nodesBackLeftOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely-1);
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
out(iy-1) = m_startNode(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
out(iy) = m_startNode(iy+1) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// -------------------- node-numbers along the back-right edge, without corners --------------------
inline ColS FineLayer::nodesBackRightOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(nely-1);
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
out(iy-1) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
out(iy) = m_startNode(iy+1) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// ------------------- node-numbers along the bottom-left edge, without corners --------------------
inline ColS FineLayer::nodesBottomLeftOpenEdge()
{
ColS out(m_nelz(0)-1);
for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz )
out(iz-1) = m_startNode(0) + iz * (m_nelx(0)+1);
return out;
}
// ------------------- node-numbers along the bottom-right edge, without corners -------------------
inline ColS FineLayer::nodesBottomRightOpenEdge()
{
ColS out(m_nelz(0)-1);
for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz )
out(iz-1) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0)+1);
return out;
}
// --------------------- node-numbers along the top-left edge, without corners ---------------------
inline ColS FineLayer::nodesTopLeftOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelz(nely-1)-1);
for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz )
out(iz-1) = m_startNode(nely) + iz * (m_nelx(nely-1)+1);
return out;
}
// -------------------- node-numbers along the top-right edge, without corners ---------------------
inline ColS FineLayer::nodesTopRightOpenEdge()
{
size_t nely = static_cast(m_nhy.size());
ColS out(m_nelz(nely-1)-1);
for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz )
out(iz-1) = m_startNode(nely) + m_nelx(nely-1) + iz * (m_nelx(nely-1)+1);
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline ColS FineLayer::nodesBottomFrontOpenEdge() { return nodesFrontBottomOpenEdge(); }
inline ColS FineLayer::nodesBottomBackOpenEdge() { return nodesBackBottomOpenEdge(); }
inline ColS FineLayer::nodesTopFrontOpenEdge() { return nodesFrontTopOpenEdge(); }
inline ColS FineLayer::nodesTopBackOpenEdge() { return nodesBackTopOpenEdge(); }
inline ColS FineLayer::nodesLeftBottomOpenEdge() { return nodesBottomLeftOpenEdge(); }
inline ColS FineLayer::nodesLeftFrontOpenEdge() { return nodesFrontLeftOpenEdge(); }
inline ColS FineLayer::nodesLeftBackOpenEdge() { return nodesBackLeftOpenEdge(); }
inline ColS FineLayer::nodesLeftTopOpenEdge() { return nodesTopLeftOpenEdge(); }
inline ColS FineLayer::nodesRightBottomOpenEdge() { return nodesBottomRightOpenEdge(); }
inline ColS FineLayer::nodesRightTopOpenEdge() { return nodesTopRightOpenEdge(); }
inline ColS FineLayer::nodesRightFrontOpenEdge() { return nodesFrontRightOpenEdge(); }
inline ColS FineLayer::nodesRightBackOpenEdge() { return nodesBackRightOpenEdge(); }
// -------------------------- node-number of the front-bottom-left corner --------------------------
inline size_t FineLayer::nodesFrontBottomLeftCorner()
{
return m_startNode(0);
}
// ------------------------- node-number of the front-bottom-right corner --------------------------
inline size_t FineLayer::nodesFrontBottomRightCorner()
{
return m_startNode(0) + m_nelx(0);
}
// --------------------------- node-number of the front-top-left corner ----------------------------
inline size_t FineLayer::nodesFrontTopLeftCorner()
{
size_t nely = static_cast(m_nhy.size());
return m_startNode(nely);
}
// --------------------------- node-number of the front-top-right corner ---------------------------
inline size_t FineLayer::nodesFrontTopRightCorner()
{
size_t nely = static_cast(m_nhy.size());
return m_startNode(nely) + m_nelx(nely-1);
}
// -------------------------- node-number of the back-bottom-left corner ---------------------------
inline size_t FineLayer::nodesBackBottomLeftCorner()
{
return m_startNode(0) + (m_nelx(0)+1)*(m_nelz(0));
}
// -------------------------- node-number of the back-bottom-right corner --------------------------
inline size_t FineLayer::nodesBackBottomRightCorner()
{
return m_startNode(0) + m_nelx(0) + (m_nelx(0)+1)*(m_nelz(0));
}
// ---------------------------- node-number of the back-top-left corner ----------------------------
inline size_t FineLayer::nodesBackTopLeftCorner()
{
size_t nely = static_cast(m_nhy.size());
return m_startNode(nely) + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
}
// --------------------------- node-number of the back-top-right corner ----------------------------
inline size_t FineLayer::nodesBackTopRightCorner()
{
size_t nely = static_cast(m_nhy.size());
return m_startNode(nely) + m_nelx(nely-1) + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
}
// -------------------------------------------- aliases --------------------------------------------
inline size_t FineLayer::nodesFrontLeftBottomCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomFrontLeftCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomLeftFrontCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftFrontBottomCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftBottomFrontCorner() { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesFrontRightBottomCorner() { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesBottomFrontRightCorner() { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesBottomRightFrontCorner() { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesRightFrontBottomCorner() { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesRightBottomFrontCorner() { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesFrontLeftTopCorner() { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesTopFrontLeftCorner() { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesTopLeftFrontCorner() { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesLeftFrontTopCorner() { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesLeftTopFrontCorner() { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesFrontRightTopCorner() { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesTopFrontRightCorner() { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesTopRightFrontCorner() { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesRightFrontTopCorner() { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesRightTopFrontCorner() { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesBackLeftBottomCorner() { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomBackLeftCorner() { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomLeftBackCorner() { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftBackBottomCorner() { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftBottomBackCorner() { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesBackRightBottomCorner() { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesBottomBackRightCorner() { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesBottomRightBackCorner() { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesRightBackBottomCorner() { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesRightBottomBackCorner() { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesBackLeftTopCorner() { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesTopBackLeftCorner() { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesTopLeftBackCorner() { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesLeftBackTopCorner() { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesLeftTopBackCorner() { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesBackRightTopCorner() { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesTopBackRightCorner() { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesTopRightBackCorner() { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesRightBackTopCorner() { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesRightTopBackCorner() { return nodesBackTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline MatS FineLayer::nodesPeriodic()
{
// faces
ColS fro = nodesFrontFace();
ColS bck = nodesBackFace();
ColS lft = nodesLeftFace();
ColS rgt = nodesRightFace();
ColS bot = nodesBottomFace();
ColS top = nodesTopFace();
// edges
ColS froBot = nodesFrontBottomOpenEdge();
ColS froTop = nodesFrontTopOpenEdge();
ColS froLft = nodesFrontLeftOpenEdge();
ColS froRgt = nodesFrontRightOpenEdge();
ColS bckBot = nodesBackBottomOpenEdge();
ColS bckTop = nodesBackTopOpenEdge();
ColS bckLft = nodesBackLeftOpenEdge();
ColS bckRgt = nodesBackRightOpenEdge();
ColS botLft = nodesBottomLeftOpenEdge();
ColS botRgt = nodesBottomRightOpenEdge();
ColS topLft = nodesTopLeftOpenEdge();
ColS topRgt = nodesTopRightOpenEdge();
// allocate nodal ties
// - number of tying per category
size_t tface = fro.size() + lft.size() + bot.size();
size_t tedge = 3*froBot.size() + 3*froLft.size() + 3*botLft.size();
size_t tnode = 7;
// - allocate
MatS out(tface+tedge+tnode, 2);
// counter
size_t i = 0;
// tie all corners to one corner
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopLeftCorner(); ++i;
// tie all corresponding edges to each other (exclude corners)
for ( auto j = 0 ; j(nodePer.rows());
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for ( size_t i = 0 ; i < nper ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
out(nodePer(i,1),j) = out(nodePer(i,0),j);
// renumber "out" to be sequential
return GooseFEM::Mesh::renumber(out);
}
// =================================================================================================
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/src/GooseFEM/MeshHex8.h b/src/GooseFEM/MeshHex8.h
index 2a7b3f6..8290448 100644
--- a/src/GooseFEM/MeshHex8.h
+++ b/src/GooseFEM/MeshHex8.h
@@ -1,327 +1,329 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHHEX8_H
#define GOOSEFEM_MESHHEX8_H
// -------------------------------------------------------------------------------------------------
#include "Mesh.h"
// ===================================== GooseFEM::Mesh::Hex8 ======================================
namespace GooseFEM {
namespace Mesh {
namespace Hex8 {
// ========================================== REGULAR MESH =========================================
class Regular
{
private:
double m_h; // elementary element edge-size (in all directions)
- size_t m_nelx; // number of element x-direction (length == "m_nelx * m_h")
- size_t m_nely; // number of element y-direction (length == "m_nely * m_h")
- size_t m_nelz; // number of element z-direction (length == "m_nely * m_h")
+ size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h")
+ size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h")
+ size_t m_nelz; // number of elements in z-direction (length == "m_nely * m_h")
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
size_t m_nne=8; // number of nodes-per-element
size_t m_ndim=3; // number of dimensions
public:
- // mesh with "nelx*nely*nelz" 'pixels' of edge size "h"
+ // mesh with "nelx*nely*nelz" 'elements' of edge size "h"
Regular(size_t nelx, size_t nely, size_t nelz, double h=1.);
// sizes
size_t nelem(); // number of elements
size_t nnode(); // number of nodes
size_t nne(); // number of nodes-per-element
size_t ndim(); // number of dimensions
// mesh
MatD coor(); // nodal positions [nnode ,ndim]
MatS conn(); // connectivity [nelem ,nne ]
// boundary nodes: planes
ColS nodesFront(); // node-numbers along the front plane
ColS nodesBack(); // node-numbers along the back plane
ColS nodesLeft(); // node-numbers along the left plane
ColS nodesRight(); // node-numbers along the right plane
ColS nodesBottom(); // node-numbers along the bottom plane
ColS nodesTop(); // node-numbers along the top plane
// boundary nodes: faces
ColS nodesFrontFace(); // node-numbers along the front face
ColS nodesBackFace(); // node-numbers along the back face
ColS nodesLeftFace(); // node-numbers along the left face
ColS nodesRightFace(); // node-numbers along the right face
ColS nodesBottomFace(); // node-numbers along the bottom face
ColS nodesTopFace(); // node-numbers along the top face
// boundary nodes: edges
ColS nodesFrontBottomEdge(); // node-numbers along the front - bottom edge
ColS nodesFrontTopEdge(); // node-numbers along the front - top edge
ColS nodesFrontLeftEdge(); // node-numbers along the front - left edge
ColS nodesFrontRightEdge(); // node-numbers along the front - right edge
ColS nodesBackBottomEdge(); // node-numbers along the back - bottom edge
ColS nodesBackTopEdge(); // node-numbers along the back - top edge
ColS nodesBackLeftEdge(); // node-numbers along the back - left edge
ColS nodesBackRightEdge(); // node-numbers along the back - right edge
ColS nodesBottomLeftEdge(); // node-numbers along the bottom - left edge
ColS nodesBottomRightEdge(); // node-numbers along the bottom - right edge
ColS nodesTopLeftEdge(); // node-numbers along the top - left edge
ColS nodesTopRightEdge(); // node-numbers along the top - right edge
// boundary nodes: edges (aliases)
ColS nodesBottomFrontEdge(); // alias, see above: nodesFrontBottomEdge
ColS nodesBottomBackEdge(); // alias, see above: nodesBackBottomEdge
ColS nodesTopFrontEdge(); // alias, see above: nodesFrontTopEdge
ColS nodesTopBackEdge(); // alias, see above: nodesBackTopEdge
ColS nodesLeftBottomEdge(); // alias, see above: nodesBottomLeftEdge
ColS nodesLeftFrontEdge(); // alias, see above: nodesFrontLeftEdge
ColS nodesLeftBackEdge(); // alias, see above: nodesBackLeftEdge
ColS nodesLeftTopEdge(); // alias, see above: nodesTopLeftEdge
ColS nodesRightBottomEdge(); // alias, see above: nodesBottomRightEdge
ColS nodesRightTopEdge(); // alias, see above: nodesTopRightEdge
ColS nodesRightFrontEdge(); // alias, see above: nodesFrontRightEdge
ColS nodesRightBackEdge(); // alias, see above: nodesBackRightEdge
// boundary nodes: edges, without corners
ColS nodesFrontBottomOpenEdge(); // node-numbers along the front - bottom edge
ColS nodesFrontTopOpenEdge(); // node-numbers along the front - top edge
ColS nodesFrontLeftOpenEdge(); // node-numbers along the front - left edge
ColS nodesFrontRightOpenEdge(); // node-numbers along the front - right edge
ColS nodesBackBottomOpenEdge(); // node-numbers along the back - bottom edge
ColS nodesBackTopOpenEdge(); // node-numbers along the back - top edge
ColS nodesBackLeftOpenEdge(); // node-numbers along the back - left edge
ColS nodesBackRightOpenEdge(); // node-numbers along the back - right edge
ColS nodesBottomLeftOpenEdge(); // node-numbers along the bottom - left edge
ColS nodesBottomRightOpenEdge(); // node-numbers along the bottom - right edge
ColS nodesTopLeftOpenEdge(); // node-numbers along the top - left edge
ColS nodesTopRightOpenEdge(); // node-numbers along the top - right edge
// boundary nodes: edges, without corners (aliases)
ColS nodesBottomFrontOpenEdge(); // alias, see above: nodesFrontBottomOpenEdge
ColS nodesBottomBackOpenEdge(); // alias, see above: nodesBackBottomOpenEdge
ColS nodesTopFrontOpenEdge(); // alias, see above: nodesFrontTopOpenEdge
ColS nodesTopBackOpenEdge(); // alias, see above: nodesBackTopOpenEdge
ColS nodesLeftBottomOpenEdge(); // alias, see above: nodesBottomLeftOpenEdge
ColS nodesLeftFrontOpenEdge(); // alias, see above: nodesFrontLeftOpenEdge
ColS nodesLeftBackOpenEdge(); // alias, see above: nodesBackLeftOpenEdge
ColS nodesLeftTopOpenEdge(); // alias, see above: nodesTopLeftOpenEdge
ColS nodesRightBottomOpenEdge(); // alias, see above: nodesBottomRightOpenEdge
ColS nodesRightTopOpenEdge(); // alias, see above: nodesTopRightOpenEdge
ColS nodesRightFrontOpenEdge(); // alias, see above: nodesFrontRightOpenEdge
ColS nodesRightBackOpenEdge(); // alias, see above: nodesBackRightOpenEdge
// boundary nodes: corners
size_t nodesFrontBottomLeftCorner(); // node-number of the front - bottom - left corner
size_t nodesFrontBottomRightCorner(); // node-number of the front - bottom - right corner
size_t nodesFrontTopLeftCorner(); // node-number of the front - top - left corner
size_t nodesFrontTopRightCorner(); // node-number of the front - top - right corner
size_t nodesBackBottomLeftCorner(); // node-number of the back - bottom - left corner
size_t nodesBackBottomRightCorner(); // node-number of the back - bottom - right corner
size_t nodesBackTopLeftCorner(); // node-number of the back - top - left corner
size_t nodesBackTopRightCorner(); // node-number of the back - top - right corner
// boundary nodes: corners (aliases)
size_t nodesFrontLeftBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesBottomFrontLeftCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesBottomLeftFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesLeftFrontBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesLeftBottomFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesFrontRightBottomCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesBottomFrontRightCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesBottomRightFrontCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesRightFrontBottomCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesRightBottomFrontCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesFrontLeftTopCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesTopFrontLeftCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesTopLeftFrontCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesLeftFrontTopCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesLeftTopFrontCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesFrontRightTopCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesTopFrontRightCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesTopRightFrontCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesRightFrontTopCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesRightTopFrontCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesBackLeftBottomCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesBottomBackLeftCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesBottomLeftBackCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesLeftBackBottomCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesLeftBottomBackCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesBackRightBottomCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesBottomBackRightCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesBottomRightBackCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesRightBackBottomCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesRightBottomBackCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesBackLeftTopCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesTopBackLeftCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesTopLeftBackCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesLeftBackTopCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesLeftTopBackCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesBackRightTopCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesTopBackRightCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesTopRightBackCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesRightBackTopCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesRightTopBackCorner(); // alias, see above: nodesBackTopRightCorner
// periodicity
MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent)
size_t nodesOrigin(); // front-bottom-left node, used as reference for periodicity
MatS dofs(); // DOF-numbers for each component of each node (sequential)
MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated
};
// ====================== MESH WITH A FINE LAYER THAT EXPONENTIALLY COARSENS =======================
class FineLayer
{
private:
double m_h; // elementary element edge-size (in all directions)
double m_Lx, m_Lz; // mesh size in "x" and "z"
ColS m_nelx, m_nelz; // number of elements in "x" and "z" (per el.layer in "y")
ColS m_nnd; // total number of nodes in the main node layer (per nd.layer in "y")
ColS m_nhx, m_nhy, m_nhz; // element size in each direction (per el.layer in "y")
- ColI m_refine; // refine direction (-1: no refine, 0: "x", ...) (per el.layer in "y")
+ ColI m_refine; // refine direction (-1:no refine, 0:"x", 2:"z") (per el.layer in "y")
ColS m_startElem; // start element (per el.layer in "y")
ColS m_startNode; // start node (per nd.layer in "y")
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
size_t m_nne=8; // number of nodes-per-element
size_t m_ndim=3; // number of dimensions
public:
- // mesh with "nelx*nely*nelz" 'pixels' of edge size "h"; elements are coarsened in "y"-direction
+ // mesh with "nelx*nely*nelz" elements of edge size "h"; elements are coarsened in "y"-direction
FineLayer(size_t nelx, size_t nely, size_t nelz, double h=1., size_t nfine=1);
// sizes
size_t nelem(); // number of elements
size_t nnode(); // number of nodes
size_t nne(); // number of nodes-per-element
size_t ndim(); // number of dimensions
size_t shape(size_t i); // actual shape in a certain direction
// mesh
MatD coor(); // nodal positions [nnode ,ndim]
MatS conn(); // connectivity [nelem ,nne ]
+ // element sets
+ ColS elementsMiddleLayer(); // elements in the middle, fine, layer
// boundary nodes: planes
ColS nodesFront(); // node-numbers along the front plane
ColS nodesBack(); // node-numbers along the back plane
ColS nodesLeft(); // node-numbers along the left plane
ColS nodesRight(); // node-numbers along the right plane
ColS nodesBottom(); // node-numbers along the bottom plane
ColS nodesTop(); // node-numbers along the top plane
// boundary nodes: faces
ColS nodesFrontFace(); // node-numbers along the front face
ColS nodesBackFace(); // node-numbers along the back face
ColS nodesLeftFace(); // node-numbers along the left face
ColS nodesRightFace(); // node-numbers along the right face
ColS nodesBottomFace(); // node-numbers along the bottom face
ColS nodesTopFace(); // node-numbers along the top face
// boundary nodes: edges
ColS nodesFrontBottomEdge(); // node-numbers along the front - bottom edge
ColS nodesFrontTopEdge(); // node-numbers along the front - top edge
ColS nodesFrontLeftEdge(); // node-numbers along the front - left edge
ColS nodesFrontRightEdge(); // node-numbers along the front - right edge
ColS nodesBackBottomEdge(); // node-numbers along the back - bottom edge
ColS nodesBackTopEdge(); // node-numbers along the back - top edge
ColS nodesBackLeftEdge(); // node-numbers along the back - left edge
ColS nodesBackRightEdge(); // node-numbers along the back - right edge
ColS nodesBottomLeftEdge(); // node-numbers along the bottom - left edge
ColS nodesBottomRightEdge(); // node-numbers along the bottom - right edge
ColS nodesTopLeftEdge(); // node-numbers along the top - left edge
ColS nodesTopRightEdge(); // node-numbers along the top - right edge
// boundary nodes: edges (aliases)
ColS nodesBottomFrontEdge(); // alias, see above: nodesFrontBottomEdge
ColS nodesBottomBackEdge(); // alias, see above: nodesBackBottomEdge
ColS nodesTopFrontEdge(); // alias, see above: nodesFrontTopEdge
ColS nodesTopBackEdge(); // alias, see above: nodesBackTopEdge
ColS nodesLeftBottomEdge(); // alias, see above: nodesBottomLeftEdge
ColS nodesLeftFrontEdge(); // alias, see above: nodesFrontLeftEdge
ColS nodesLeftBackEdge(); // alias, see above: nodesBackLeftEdge
ColS nodesLeftTopEdge(); // alias, see above: nodesTopLeftEdge
ColS nodesRightBottomEdge(); // alias, see above: nodesBottomRightEdge
ColS nodesRightTopEdge(); // alias, see above: nodesTopRightEdge
ColS nodesRightFrontEdge(); // alias, see above: nodesFrontRightEdge
ColS nodesRightBackEdge(); // alias, see above: nodesBackRightEdge
// boundary nodes: edges, without corners
ColS nodesFrontBottomOpenEdge(); // node-numbers along the front - bottom edge
ColS nodesFrontTopOpenEdge(); // node-numbers along the front - top edge
ColS nodesFrontLeftOpenEdge(); // node-numbers along the front - left edge
ColS nodesFrontRightOpenEdge(); // node-numbers along the front - right edge
ColS nodesBackBottomOpenEdge(); // node-numbers along the back - bottom edge
ColS nodesBackTopOpenEdge(); // node-numbers along the back - top edge
ColS nodesBackLeftOpenEdge(); // node-numbers along the back - left edge
ColS nodesBackRightOpenEdge(); // node-numbers along the back - right edge
ColS nodesBottomLeftOpenEdge(); // node-numbers along the bottom - left edge
ColS nodesBottomRightOpenEdge(); // node-numbers along the bottom - right edge
ColS nodesTopLeftOpenEdge(); // node-numbers along the top - left edge
ColS nodesTopRightOpenEdge(); // node-numbers along the top - right edge
// boundary nodes: edges, without corners (aliases)
ColS nodesBottomFrontOpenEdge(); // alias, see above: nodesFrontBottomOpenEdge
ColS nodesBottomBackOpenEdge(); // alias, see above: nodesBackBottomOpenEdge
ColS nodesTopFrontOpenEdge(); // alias, see above: nodesFrontTopOpenEdge
ColS nodesTopBackOpenEdge(); // alias, see above: nodesBackTopOpenEdge
ColS nodesLeftBottomOpenEdge(); // alias, see above: nodesBottomLeftOpenEdge
ColS nodesLeftFrontOpenEdge(); // alias, see above: nodesFrontLeftOpenEdge
ColS nodesLeftBackOpenEdge(); // alias, see above: nodesBackLeftOpenEdge
ColS nodesLeftTopOpenEdge(); // alias, see above: nodesTopLeftOpenEdge
ColS nodesRightBottomOpenEdge(); // alias, see above: nodesBottomRightOpenEdge
ColS nodesRightTopOpenEdge(); // alias, see above: nodesTopRightOpenEdge
ColS nodesRightFrontOpenEdge(); // alias, see above: nodesFrontRightOpenEdge
ColS nodesRightBackOpenEdge(); // alias, see above: nodesBackRightOpenEdge
// boundary nodes: corners
size_t nodesFrontBottomLeftCorner(); // node-number of the front - bottom - left corner
size_t nodesFrontBottomRightCorner(); // node-number of the front - bottom - right corner
size_t nodesFrontTopLeftCorner(); // node-number of the front - top - left corner
size_t nodesFrontTopRightCorner(); // node-number of the front - top - right corner
size_t nodesBackBottomLeftCorner(); // node-number of the back - bottom - left corner
size_t nodesBackBottomRightCorner(); // node-number of the back - bottom - right corner
size_t nodesBackTopLeftCorner(); // node-number of the back - top - left corner
size_t nodesBackTopRightCorner(); // node-number of the back - top - right corner
// boundary nodes: corners (aliases)
size_t nodesFrontLeftBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesBottomFrontLeftCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesBottomLeftFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesLeftFrontBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesLeftBottomFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner
size_t nodesFrontRightBottomCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesBottomFrontRightCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesBottomRightFrontCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesRightFrontBottomCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesRightBottomFrontCorner(); // alias, see above: nodesFrontBottomRightCorner
size_t nodesFrontLeftTopCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesTopFrontLeftCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesTopLeftFrontCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesLeftFrontTopCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesLeftTopFrontCorner(); // alias, see above: nodesFrontTopLeftCorner
size_t nodesFrontRightTopCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesTopFrontRightCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesTopRightFrontCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesRightFrontTopCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesRightTopFrontCorner(); // alias, see above: nodesFrontTopRightCorner
size_t nodesBackLeftBottomCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesBottomBackLeftCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesBottomLeftBackCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesLeftBackBottomCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesLeftBottomBackCorner(); // alias, see above: nodesBackBottomLeftCorner
size_t nodesBackRightBottomCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesBottomBackRightCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesBottomRightBackCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesRightBackBottomCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesRightBottomBackCorner(); // alias, see above: nodesBackBottomRightCorner
size_t nodesBackLeftTopCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesTopBackLeftCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesTopLeftBackCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesLeftBackTopCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesLeftTopBackCorner(); // alias, see above: nodesBackTopLeftCorner
size_t nodesBackRightTopCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesTopBackRightCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesTopRightBackCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesRightBackTopCorner(); // alias, see above: nodesBackTopRightCorner
size_t nodesRightTopBackCorner(); // alias, see above: nodesBackTopRightCorner
// periodicity
MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent)
size_t nodesOrigin(); // front-bottom-left node, used as reference for periodicity
MatS dofs(); // DOF-numbers for each component of each node (sequential)
MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated
};
// =================================================================================================
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/src/GooseFEM/MeshQuad4.cpp b/src/GooseFEM/MeshQuad4.cpp
index 9e23c3b..58e52e5 100644
--- a/src/GooseFEM/MeshQuad4.cpp
+++ b/src/GooseFEM/MeshQuad4.cpp
@@ -1,876 +1,922 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHQUAD4_CPP
#define GOOSEFEM_MESHQUAD4_CPP
// -------------------------------------------------------------------------------------------------
#include "MeshQuad4.h"
// ===================================== GooseFEM::Mesh::Quad4 =====================================
namespace GooseFEM {
namespace Mesh {
namespace Quad4 {
// ===================================== CLASS - REGULAR MESH ======================================
// ------------------------------------------ constructor ------------------------------------------
-inline Regular::Regular(size_t nx, size_t ny, double h):
-m_nx(nx), m_ny(ny), m_h(h)
+inline Regular::Regular(size_t nelx, size_t nely, double h):
+m_h(h), m_nelx(nelx), m_nely(nely)
{
- assert( m_nx >= 1 );
- assert( m_ny >= 1 );
+ assert( m_nelx >= 1 );
+ assert( m_nely >= 1 );
- m_nnode = (m_nx+1) * (m_ny+1);
- m_nelem = m_nx * m_ny ;
+ m_nnode = (m_nelx+1) * (m_nely+1);
+ m_nelem = m_nelx * m_nely ;
}
// -------------------------------------- number of elements ---------------------------------------
inline size_t Regular::nelem()
{
return m_nelem;
}
// ---------------------------------------- number of nodes ----------------------------------------
inline size_t Regular::nnode()
{
return m_nnode;
}
// ---------------------------------- number of nodes per element ----------------------------------
inline size_t Regular::nne()
{
return m_nne;
}
// ------------------------------------- number of dimensions --------------------------------------
inline size_t Regular::ndim()
{
return m_ndim;
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline MatD Regular::coor()
{
- MatD coor(m_nnode,m_ndim);
+ MatD out(m_nnode, m_ndim);
- ColD x = ColD::LinSpaced(m_nx+1, 0.0, m_h*static_cast(m_nx));
- ColD y = ColD::LinSpaced(m_ny+1, 0.0, m_h*static_cast(m_ny));
+ ColD x = ColD::LinSpaced(m_nelx+1, 0.0, m_h*static_cast(m_nelx));
+ ColD y = ColD::LinSpaced(m_nely+1, 0.0, m_h*static_cast(m_nely));
size_t inode = 0;
- for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) {
- for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) {
- coor(inode,0) = x(ix);
- coor(inode,1) = y(iy);
+ for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) {
+ for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) {
+ out(inode,0) = x(ix);
+ out(inode,1) = y(iy);
++inode;
}
}
- return coor;
+ return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline MatS Regular::conn()
{
- MatS conn(m_nelem,m_nne);
+ MatS out(m_nelem,m_nne);
size_t ielem = 0;
- for ( size_t iy = 0 ; iy < m_ny ; ++iy ) {
- for ( size_t ix = 0 ; ix < m_nx ; ++ix ) {
- conn(ielem,0) = (iy+0)*(m_nx+1) + (ix+0);
- conn(ielem,1) = (iy+0)*(m_nx+1) + (ix+1);
- conn(ielem,3) = (iy+1)*(m_nx+1) + (ix+0);
- conn(ielem,2) = (iy+1)*(m_nx+1) + (ix+1);
+ for ( size_t iy = 0 ; iy < m_nely ; ++iy ) {
+ for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) {
+ out(ielem,0) = (iy )*(m_nelx+1) + (ix );
+ out(ielem,1) = (iy )*(m_nelx+1) + (ix+1);
+ out(ielem,3) = (iy+1)*(m_nelx+1) + (ix );
+ out(ielem,2) = (iy+1)*(m_nelx+1) + (ix+1);
++ielem;
}
}
- return conn;
+ return out;
}
// ------------------------------ node-numbers along the bottom edge -------------------------------
inline ColS Regular::nodesBottomEdge()
{
- ColS nodes(m_nx+1);
+ ColS out(m_nelx+1);
- for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix )
- nodes(ix) = ix;
+ for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
+ out(ix) = ix;
- return nodes;
+ return out;
}
// -------------------------------- node-numbers along the top edge --------------------------------
inline ColS Regular::nodesTopEdge()
{
- ColS nodes(m_nx+1);
+ ColS out(m_nelx+1);
- for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix )
- nodes(ix) = ix + m_ny*(m_nx+1);
+ for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
+ out(ix) = ix + m_nely*(m_nelx+1);
- return nodes;
+ return out;
}
-// ----------------------------- node-numbers along the node left edge -----------------------------
+// ------------------------------- node-numbers along the left edge --------------------------------
inline ColS Regular::nodesLeftEdge()
{
- ColS nodes(m_ny+1);
+ ColS out(m_nely+1);
- for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy )
- nodes(iy) = iy*(m_nx+1);
+ for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
+ out(iy) = iy*(m_nelx+1);
- return nodes;
+ return out;
}
// ------------------------------- node-numbers along the right edge -------------------------------
inline ColS Regular::nodesRightEdge()
{
- ColS nodes(m_ny+1);
+ ColS out(m_nely+1);
+
+ for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
+ out(iy) = iy*(m_nelx+1) + m_nelx;
+
+ return out;
+}
+
+// ---------------------- node-numbers along the bottom edge, without corners ----------------------
+
+inline ColS Regular::nodesBottomOpenEdge()
+{
+ ColS out(m_nelx-1);
+
+ for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
+ out(ix-1) = ix;
+
+ return out;
+}
+
+// ----------------------- node-numbers along the top edge, without corners ------------------------
+
+inline ColS Regular::nodesTopOpenEdge()
+{
+ ColS out(m_nelx-1);
+
+ for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
+ out(ix-1) = ix + m_nely*(m_nelx+1);
+
+ return out;
+}
+
+// ----------------------- node-numbers along the left edge, without corners -----------------------
+
+inline ColS Regular::nodesLeftOpenEdge()
+{
+ ColS out(m_nely-1);
- for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy )
- nodes(iy) = iy*(m_nx+1) + m_nx;
+ for ( size_t iy = 1 ; iy < m_nely ; ++iy )
+ out(iy-1) = iy*(m_nelx+1);
- return nodes;
+ return out;
+}
+
+// ---------------------- node-numbers along the right edge, without corners -----------------------
+
+inline ColS Regular::nodesRightOpenEdge()
+{
+ ColS out(m_nely-1);
+
+ for ( size_t iy = 1 ; iy < m_nely ; ++iy )
+ out(iy-1) = iy*(m_nelx+1) + m_nelx;
+
+ return out;
}
// ----------------------------- node-number of the bottom-left corner -----------------------------
inline size_t Regular::nodesBottomLeftCorner()
{
return 0;
}
// ---------------------------- node-number of the bottom-right corner -----------------------------
inline size_t Regular::nodesBottomRightCorner()
{
- return m_nx;
+ return m_nelx;
}
// ------------------------------ node-number of the top-left corner -------------------------------
inline size_t Regular::nodesTopLeftCorner()
{
- return m_ny*(m_nx+1);
+ return m_nely*(m_nelx+1);
}
// ------------------------------ node-number of the top-right corner ------------------------------
inline size_t Regular::nodesTopRightCorner()
{
- return (m_ny+1)*(m_nx+1) - 1;
+ return m_nely*(m_nelx+1) + m_nelx;
}
// ----------------------------- node-number of the corners (aliases) ------------------------------
inline size_t Regular::nodesLeftBottomCorner() { return nodesBottomLeftCorner(); }
inline size_t Regular::nodesLeftTopCorner() { return nodesTopLeftCorner(); }
inline size_t Regular::nodesRightBottomCorner() { return nodesBottomRightCorner(); }
inline size_t Regular::nodesRightTopCorner() { return nodesTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline MatS Regular::nodesPeriodic()
{
- // edges
- ColS bot = nodesBottomEdge();
- ColS top = nodesTopEdge();
- ColS rgt = nodesRightEdge();
- ColS lft = nodesLeftEdge();
+ // edges (without corners)
+ ColS bot = nodesBottomOpenEdge();
+ ColS top = nodesTopOpenEdge();
+ ColS lft = nodesLeftOpenEdge();
+ ColS rgt = nodesRightOpenEdge();
// allocate nodal ties
- MatS nodes(bot.size()-2 + lft.size()-2 + 3, 2);
+ // - number of tying per category
+ size_t tedge = bot.size() + lft.size();
+ size_t tnode = 3;
+ // - allocate
+ MatS out(tedge+tnode, 2);
// counter
size_t i = 0;
// tie all corners to one corner
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesBottomRightCorner(); ++i;
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopRightCorner(); ++i;
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopLeftCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i;
- // tie edges to each other (exclude corners)
- for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; }
- for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; }
+ // tie all corresponding edges to each other
+ for ( auto j = 0 ; j(nodePer.rows());
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for ( size_t i = 0 ; i < nper ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
out(nodePer(i,1),j) = out(nodePer(i,0),j);
// renumber "out" to be sequential
return GooseFEM::Mesh::renumber(out);
}
// ==================================== CLASS - FINELAYER MESH =====================================
// ------------------------------------------ constructor ------------------------------------------
-inline FineLayer::FineLayer(size_t nx, size_t ny, double h, size_t nfine, size_t nskip):
-m_h(h), m_nx(nx)
+inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine):
+m_h(h), m_nelx(nelx)
{
- assert( nx >= 1 );
- assert( ny >= 1 );
+ // basic assumptions
+ assert( nelx >= 1 );
+ assert( nely >= 1 );
+
+ // store basic info
+ m_Lx = m_h * static_cast(nelx);
- // local counters
- size_t N;
+ // compute element size in y-direction (use symmetry, compute upper half)
+ // ----------------------------------------------------------------------
- // compute the element size : based on symmetric half
- // --------------------------------------------------
+ // temporary variables
+ size_t nmin, ntot;
+ ColS nhx(nely), nhy(nely);
+ ColI refine(nely);
- // convert height to half of the height
- if ( ny % 2 == 0 ) ny = ny /2;
- else ny = (ny+1)/2;
+ // minimum height in y-direction (half of the height because of symmetry)
+ if ( nely % 2 == 0 ) nmin = nely /2;
+ else nmin = (nely +1)/2;
- // convert to half the number of fine layer (minimum 1)
+ // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half)
if ( nfine % 2 == 0 ) nfine = nfine /2 + 1;
else nfine = (nfine+1)/2;
if ( nfine < 1 ) nfine = 1;
+ if ( nfine > nmin ) nfine = nmin;
+
+ // initialize to state with only fine elements
+ nhx .setOnes();
+ nhy .setOnes();
+ refine.setConstant(-1);
+
+ // loop over element layers in y-direction, try to coarsen using these rules:
+ // (1) element size in y-direction <= distance to origin in y-direction
+ // (2) element size in x-direction should fit the total number of elements in x-direction
+ // (3) a certain number of layers have the minimum size "1" (are fine)
+ for ( size_t iy = nfine ; ; )
+ {
+ // initialize current size in y-direction
+ if ( iy == nfine ) ntot = nfine;
+ // check to stop
+ if ( iy >= nely or ntot >= nmin ) { nely = iy; break; }
- // check the number of fine layers from the center
- assert( nfine <= ny );
-
- // define arrays to determine to coarsen
- ColS n (ny+1); // size of the element (number of times "h"): even sizes -> coarsen
- ColS nsum(ny+1); // current 'size' of the mesh in y-direction
-
- // initialize all elements of size "1"
- n.setOnes();
+ // rules (1,2) satisfied: coarsen in x-direction
+ if ( 3*nhy(iy) <= ntot and nelx%(3*nhx(iy)) == 0 )
+ {
+ refine (iy ) = 0;
+ nhy (iy ) *= 2;
+ nhy.segment(iy+1,nely-iy-1) *= 3;
+ nhx.segment(iy ,nely-iy ) *= 3;
+ }
- // size of the mesh in y-direction after "i" elements
- // - initialize
- nsum(0) = 1;
- // - compute based on the current value of "n"
- for ( size_t i = 1 ; i < ny+1 ; ++i )
- nsum(i) = nsum(i-1) + n(i-1);
-
- // loop in vertical direction and check to coarsen; rules:
- // * the size of the element cannot be greater than the distance
- // * the size of the element should match the horizontal direction
- // * skip a certain number of layers that have the minimum size "1" (are fine)
- // - initialize next element size
- size_t next = 1;
- // - loop to coarsen
- for ( size_t i = nfine ; i < ny+1 ; ++i )
- {
- // -- set element size, based on whether to coarsen or not
- if ( 3*n(i-1) <= nsum(i-1) and m_nx % ( 3*next ) == 0 ) { n(i) = next * 2; next *= 3; }
- else { n(i) = next; }
- // -- compute the total number of elements
- nsum(i) = nsum(i-1) + n(i);
+ // update the number of elements in y-direction
+ ntot += nhy(iy);
+ // proceed to next element layer in y-direction
+ ++iy;
+ // check to stop
+ if ( iy >= nely or ntot >= nmin ) { nely = iy; break; }
}
- // truncate to size "N" such that the height does not exceed "ny"
- for ( N = 0 ; N < ny+1 ; ++N )
- if ( nsum(N) > ny+1 )
- break;
-
- // fiddle a little bit to approximate the size as best as possible
- if ( N > 1 and N < ny )
+ // symmetrize, compute full information
+ // ------------------------------------
+
+ // allocate mesh constructor parameters
+ m_nhx .conservativeResize(nely*2-1);
+ m_nhy .conservativeResize(nely*2-1);
+ m_refine .conservativeResize(nely*2-1);
+ m_nelx .conservativeResize(nely*2-1);
+ m_nnd .conservativeResize(nely*2 );
+ m_startElem.conservativeResize(nely*2-1);
+ m_startNode.conservativeResize(nely*2 );
+
+ // fill
+ // - lower half
+ for ( size_t iy = 0 ; iy < nely ; ++iy )
{
- // - set all sequential steps to last value, make sure not to use a coarsening step
- if ( n(N-1) % 2 == 0 ) for ( size_t i = N ; i < ny+1 ; ++i ) n(i) = n(i-1)/2*3;
- else for ( size_t i = N ; i < ny+1 ; ++i ) n(i) = n(i-1);
- // - fix total size
- for ( size_t i = N ; i < ny+1 ; ++i ) nsum(i) = nsum(i-1) + n(i);
- // - truncate such that the height does not exceed "ny"
- for ( N = 0 ; N < ny ; ++N )
- if ( nsum(N) > ny+1 )
- break;
- // - check to extend one
- if ( nsum(N) - ny < ny - nsum(N-1) )
- ++N;
+ m_nhx (iy) = nhx (nely-iy-1);
+ m_nhy (iy) = nhy (nely-iy-1);
+ m_refine(iy) = refine(nely-iy-1);
}
-
- // skip the last "nskip" coarsening steps
- if ( nskip > 0 )
+ // - upper half
+ for ( size_t iy = 0 ; iy < nely-1 ; ++iy )
{
- // - counter : number of steps skipped
- size_t iskip = 0;
- // - loop to skip
- for ( size_t i = N ; i-- > 1 ; )
- {
- // -- check -> quit if the number of steps is reached
- if ( iskip >= nskip ) break;
- // -- check -> quit if the finest layer has been reached
- if ( n(i) == 1 ) break;
- // -- if coarsening step found -> remove
- if ( n(i) % 2 == 0 )
- {
- size_t k = n(i) / 2;
- for ( size_t j = i ; j < ny+1 ; ++j ) n(j) = k;
- iskip++;
- }
- }
- // - update the height after "i" elements
- for ( size_t i = 1 ; i < ny+1 ; ++i )
- nsum(i) = nsum(i-1) + n(i);
- // - truncate such that the height does not exceed "ny"
- for ( N = 0 ; N < ny+1 ; ++N )
- if ( nsum(N) > ny+1 )
- break;
+ m_nhx (iy+nely) = nhx (iy+1);
+ m_nhy (iy+nely) = nhy (iy+1);
+ m_refine(iy+nely) = refine(iy+1);
}
- // truncate
- n.conservativeResize(N);
+ // update size
+ nely = m_nhx.size();
- // compute mesh dimensions : based on the entire mesh
- // --------------------------------------------------
+ // compute the number of elements per element layer in y-direction
+ for ( size_t iy = 0 ; iy < nely ; ++iy )
+ m_nelx(iy) = nelx / m_nhx(iy);
- // symmetrize
- // - get size
- N = static_cast(n.size());
- // - allocate
- m_nh.conservativeResize(2*N-1);
- // - store symmetrized
- for ( size_t i = 0 ; i < N ; i++ ) m_nh( i ) = n(N-1-i);
- for ( size_t i = 1 ; i < N ; i++ ) m_nh(N+i-1) = n( i);
+ // compute the number of nodes per node layer in y-direction
+ for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) m_nnd(iy ) = m_nelx(iy)+1;
+ for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) m_nnd(iy+1) = m_nelx(iy)+1;
- // allocate counters
- m_startElem.conservativeResize(m_nh.size() );
- m_startNode.conservativeResize(m_nh.size()+1);
+ // compute mesh dimensions
+ // -----------------------
- // compute the dimensions of the mesh
- // - initialize
+ // initialize
m_nnode = 0;
m_nelem = 0;
- N = static_cast(m_nh.size());
m_startNode(0) = 0;
- // loop from bottom to middle : elements become finer
- for ( size_t i = 0 ; i < (N-1)/2 ; ++i )
+ // loop over element layers (bottom -> middle, elements become finer)
+ for ( size_t i = 0 ; i < (nely-1)/2 ; ++i )
{
- // - store the first element of the row
+ // - store the first element of the layer
m_startElem(i) = m_nelem;
- // - add the nodes of this row
- if ( m_nh(i)%2 == 0 ) { m_nnode += 3 * m_nx/(3*m_nh(i)/2) + 1; }
- else { m_nnode += m_nx/ m_nh(i) + 1; }
- // - add the elements of this row
- if ( m_nh(i)%2 == 0 ) { m_nelem += 4 * m_nx/(3*m_nh(i)/2); }
- else { m_nelem += m_nx/ m_nh(i) ; }
- // - store the starting node of the next row
+ // - add the nodes of this layer
+ if ( m_refine(i) == 0 ) { m_nnode += (3*m_nelx(i)+1); }
+ else { m_nnode += ( m_nelx(i)+1); }
+ // - add the elements of this layer
+ if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ); }
+ else { m_nelem += ( m_nelx(i) ); }
+ // - store the starting node of the next layer
m_startNode(i+1) = m_nnode;
}
- // loop from middle to top : elements become coarser
- for ( size_t i = (N-1)/2 ; i < N ; ++i )
+ // loop over element layers (middle -> top, elements become coarser)
+ for ( size_t i = (nely-1)/2 ; i < nely ; ++i )
{
- // - store the first element of the row
+ // - store the first element of the layer
m_startElem(i) = m_nelem;
- // - add the nodes of this row
- if ( m_nh(i)%2 == 0 ) { m_nnode += 1 + m_nx/(m_nh(i)/2) + 2 * m_nx/(3*m_nh(i)/2); }
- else { m_nnode += 1 + m_nx/ m_nh(i); }
- // - add the elements of this row
- if ( m_nh(i)%2 == 0 ) { m_nelem += 4 * m_nx/(3*m_nh(i)/2); }
- else { m_nelem += m_nx/ m_nh(i) ; }
- // - store the starting node of the next row
+ // - add the nodes of this layer
+ if ( m_refine(i) == 0 ) { m_nnode += (5*m_nelx(i)+1); }
+ else { m_nnode += ( m_nelx(i)+1); }
+ // - add the elements of this layer
+ if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ); }
+ else { m_nelem += ( m_nelx(i) ); }
+ // - store the starting node of the next layer
m_startNode(i+1) = m_nnode;
}
-
- // add the top row of nodes to the number of rows (starting node already stored)
- if ( m_nh(N-1)%2 == 0 ) { m_nnode += 1 + m_nx/(3*m_nh(N-1)/2); }
- else { m_nnode += 1 + m_nx/ m_nh(N-1); }
+ // - add the top row of nodes
+ m_nnode += m_nelx(nely-1)+1;
}
// -------------------------------------- number of elements ---------------------------------------
inline size_t FineLayer::nelem()
{
return m_nelem;
}
// ---------------------------------------- number of nodes ----------------------------------------
inline size_t FineLayer::nnode()
{
return m_nnode;
}
// ---------------------------------- number of nodes per element ----------------------------------
inline size_t FineLayer::nne()
{
return m_nne;
}
// ------------------------------------- number of dimensions --------------------------------------
inline size_t FineLayer::ndim()
{
return m_ndim;
}
// ---------------------------- actual number of nodes in one direction ----------------------------
inline size_t FineLayer::shape(size_t i)
{
- // check the index
- assert( i >= 0 and i < 2 );
-
- // horizontal direction
- if ( i == 0 )
- return m_nx;
+ assert( i >= 0 and i <= 1 );
- // vertical direction
- // - zero-initialize
- size_t n = 0;
- // - cumulative sum of the size per row
- for ( size_t i = 0 ; i < static_cast(m_nh.size()) ; ++i )
- n += m_nh(i);
+ if ( i == 0 ) return m_nelx.maxCoeff();
+ else return m_nhy .sum();
- return n;
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline MatD FineLayer::coor()
{
// allocate output
MatD out(m_nnode, m_ndim);
- // initialize position in horizontal direction (of the finest nodes)
- ColD x = ColD::LinSpaced(m_nx+1, 0.0, m_h*static_cast(m_nx));
-
- // zero-initialize current node and height: loop from bottom to top; number of rows
- double h = 0;
+ // current node, number of element layers
size_t inode = 0;
- size_t N = static_cast(m_nh.size());
+ size_t nely = static_cast(m_nhy.size());
+
+ // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
+ // - allocate
+ ColD y(nely+1);
+ // - initialize
+ y(0) = 0.0;
+ // - compute
+ for ( size_t iy = 1 ; iy < nely+1 ; ++iy )
+ y(iy) = y(iy-1) + m_nhy(iy-1) * m_h;
+
+ // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes
+ // --------------------------------------------------------------------------------------------
- // lower half
- for ( size_t i = 0 ; i < (N-1)/2 ; ++i )
+ for ( size_t iy = 0 ; ; ++iy )
{
- // - refinement row
- if ( m_nh(i)%2 == 0 )
- {
- // -- bottom row: coarse nodes
- for ( size_t j = 0 ; j < m_nx+1 ; j += 3*m_nh(i)/2 )
- {
- out(inode,0) = x(j); out(inode,1) = h*m_h; inode++;
- }
- h += static_cast(m_nh(i)/2);
- // -- middle row: part the fine nodes
- for ( size_t j = 0 ; j < m_nx ; j += 3*m_nh(i)/2 )
- {
- out(inode,0) = x(j+m_nh(i)/2); out(inode,1) = h*m_h; inode++;
- out(inode,0) = x(j+m_nh(i) ); out(inode,1) = h*m_h; inode++;
- }
- h += static_cast(m_nh(i)/2);
+ // get positions along the x- and z-axis
+ ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx);
+
+ // add nodes of the bottom layer of this element
+ for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
+ out(inode,0) = x(ix);
+ out(inode,1) = y(iy);
+ ++inode;
}
- // - normal row
- else
+
+ // stop at middle layer
+ if ( iy == (nely-1)/2 )
+ break;
+
+ // add extra nodes of the intermediate layer, for refinement in x-direction
+ if ( m_refine(iy) == 0 )
{
- for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i) )
- {
- out(inode,0) = x(j); out(inode,1) = h*m_h; inode++;
+ // - get position offset in x- and y-direction
+ double dx = m_h * static_cast(m_nhx(iy)/3);
+ double dy = m_h * static_cast(m_nhy(iy)/2);
+ // - add nodes of the intermediate layer
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
+ for ( size_t j = 0 ; j < 2 ; ++j ) {
+ out(inode,0) = x(ix) + dx * static_cast(j+1);
+ out(inode,1) = y(iy) + dy;
+ ++inode;
+ }
}
- h += static_cast(m_nh(i));
}
}
- // top half
- for ( size_t i = (N-1)/2 ; i < N ; ++i )
+ // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes
+ // --------------------------------------------------------------------------------------
+
+ for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
- // - coarsening row
- if ( m_nh(i)%2 == 0 )
- {
- // -- bottom row: fine nodes
- for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i)/2 )
- {
- out(inode,0) = x(j); out(inode,1) = h*m_h; inode++;
- }
- h += static_cast(m_nh(i)/2);
- // -- middle row: part the fine nodes
- for ( size_t j = 0 ; j < m_nx ; j += 3*m_nh(i)/2 )
- {
- out(inode,0) = x(j+m_nh(i)/2); out(inode,1) = h*m_h; inode++;
- out(inode,0) = x(j+m_nh(i) ); out(inode,1) = h*m_h; inode++;
- }
- h += static_cast(m_nh(i)/2);
- }
- // - normal row
- else
+ // get positions along the x- and z-axis
+ ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx);
+
+ // add extra nodes of the intermediate layer, for refinement in x-direction
+ if ( m_refine(iy) == 0 )
{
- for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i) )
- {
- out(inode,0) = x(j); out(inode,1) = h*m_h; inode++;
+ // - get position offset in x- and y-direction
+ double dx = m_h * static_cast(m_nhx(iy)/3);
+ double dy = m_h * static_cast(m_nhy(iy)/2);
+ // - add nodes of the intermediate layer
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
+ for ( size_t j = 0 ; j < 2 ; ++j ) {
+ out(inode,0) = x(ix) + dx * static_cast(j+1);
+ out(inode,1) = y(iy) + dy;
+ ++inode;
+ }
}
- h += static_cast(m_nh(i));
}
- }
- // top row
- if ( m_nh(N-1)%2 == 0 )
- {
- for ( size_t j = 0 ; j < m_nx+1 ; j += 3*m_nh(N-1)/2 )
- {
- out(inode,0) = x(j); out(inode,1) = h*m_h; inode++;
- }
- }
- else
- {
- for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(N-1) )
- {
- out(inode,0) = x(j); out(inode,1) = h*m_h; inode++;
+ // add nodes of the top layer of this element
+ for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
+ out(inode,0) = x(ix );
+ out(inode,1) = y(iy+1);
+ ++inode;
}
}
return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline MatS FineLayer::conn()
{
// allocate output
MatS out(m_nelem, m_nne);
- // current element, number of rows, beginning nodes
+ // current element, number of element layers, starting nodes of each node layer
size_t ielem = 0;
- size_t N = static_cast(m_nh.size());
+ size_t nely = static_cast(m_nhy.size());
size_t bot,mid,top;
- // bottom half
- for ( size_t i = 0 ; i < (N-1)/2 ; ++i )
+ // loop over all element layers
+ for ( size_t iy = 0 ; iy < nely ; ++iy )
{
- // - starting nodes of bottom and top layer (and middle layer, only relevant for even shape)
- if ( true ) bot = m_startNode(i );
- if ( m_nh(i)%2 == 0 ) mid = m_startNode(i )+m_nx/(3*m_nh(i)/2)+1;
- if ( true ) top = m_startNode(i+1);
+ // - get: starting nodes of bottom(, middle) and top layer
+ bot = m_startNode(iy );
+ mid = m_startNode(iy ) + m_nnd(iy);
+ top = m_startNode(iy+1);
- // - even sized shape: refinement layer
- if ( m_nh(i)%2 == 0 )
+ // - define connectivity: no coarsening/refinement
+ if ( m_refine(iy) == -1 )
{
- for ( size_t j = 0 ; j < m_nx/(3*m_nh(i)/2) ; ++j )
- {
- // -- lower
- out(ielem,0) = j +bot;
- out(ielem,1) = j +1+bot;
- out(ielem,2) = j*2+1+mid;
- out(ielem,3) = j*2 +mid;
- ielem++;
- // -- upper right
- out(ielem,0) = j +1+bot;
- out(ielem,1) = j*3+3+top;
- out(ielem,2) = j*3+2+top;
- out(ielem,3) = j*2+1+mid;
- ielem++;
- // -- upper middle
- out(ielem,0) = j*2 +mid;
- out(ielem,1) = j*2+1+mid;
- out(ielem,2) = j*3+2+top;
- out(ielem,3) = j*3+1+top;
- ielem++;
- // -- upper left
- out(ielem,0) = j +bot;
- out(ielem,1) = j*2 +mid;
- out(ielem,2) = j*3+1+top;
- out(ielem,3) = j*3 +top;
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
+ out(ielem,0) = bot + (ix );
+ out(ielem,1) = bot + (ix+1);
+ out(ielem,2) = top + (ix+1);
+ out(ielem,3) = top + (ix );
ielem++;
}
}
- // - odd sized shape: normal layer
- else
- {
- for ( size_t j = 0 ; j < m_nx/m_nh(i) ; ++j )
- {
- out(ielem,0) = j +bot;
- out(ielem,1) = j+1+bot;
- out(ielem,2) = j+1+top;
- out(ielem,3) = j +top;
- ielem++;
- }
- }
- }
- // top half
- for ( size_t i = (N-1)/2 ; i < N ; ++i )
- {
- // - starting nodes of bottom and top layer (and middle layer, only relevant for even shape)
- if ( true ) bot = m_startNode(i );
- if ( m_nh(i)%2 == 0 ) mid = m_startNode(i )+m_nx/(m_nh(i)/2)+1;
- if ( true ) top = m_startNode(i+1);
-
- // - even sized shape: refinement layer
- if ( m_nh(i)%2 == 0 )
+ // - define connectivity: refinement along the x-direction (below the middle layer)
+ else if ( m_refine(iy) == 0 and iy <= (nely-1)/2 )
{
- for ( size_t j = 0 ; j < m_nx/(3*m_nh(i)/2) ; ++j )
- {
- // -- lower left
- out(ielem,0) = j*3 +bot;
- out(ielem,1) = j*3+1+bot;
- out(ielem,2) = j*2 +mid;
- out(ielem,3) = j +top;
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
+ // -- bottom element
+ out(ielem,0) = bot + ( ix );
+ out(ielem,1) = bot + ( ix+1);
+ out(ielem,2) = mid + (2*ix+1);
+ out(ielem,3) = mid + (2*ix );
ielem++;
- // -- lower middle
- out(ielem,0) = j*3+1+bot;
- out(ielem,1) = j*3+2+bot;
- out(ielem,2) = j*2+1+mid;
- out(ielem,3) = j*2 +mid;
+ // -- top-right element
+ out(ielem,0) = bot + ( ix+1);
+ out(ielem,1) = top + (3*ix+3);
+ out(ielem,2) = top + (3*ix+2);
+ out(ielem,3) = mid + (2*ix+1);
ielem++;
- // -- lower right
- out(ielem,0) = j*3+2+bot;
- out(ielem,1) = j*3+3+bot;
- out(ielem,2) = j +1+top;
- out(ielem,3) = j*2+1+mid;
+ // -- top-center element
+ out(ielem,0) = mid + (2*ix );
+ out(ielem,1) = mid + (2*ix+1);
+ out(ielem,2) = top + (3*ix+2);
+ out(ielem,3) = top + (3*ix+1);
ielem++;
- // -- upper
- out(ielem,0) = j +top;
- out(ielem,1) = j*2 +mid;
- out(ielem,2) = j*2+1+mid;
- out(ielem,3) = j +1+top;
+ // -- top-left element
+ out(ielem,0) = bot + ( ix );
+ out(ielem,1) = mid + (2*ix );
+ out(ielem,2) = top + (3*ix+1);
+ out(ielem,3) = top + (3*ix );
ielem++;
}
}
- // - odd sized shape: normal layer
- else
+
+ // - define connectivity: coarsening along the x-direction (above the middle layer)
+ else if ( m_refine(iy) == 0 and iy > (nely-1)/2 )
{
- for ( size_t j = 0 ; j < m_nx/m_nh(i) ; ++j )
- {
- out(ielem,0) = j +bot;
- out(ielem,1) = j+1+bot;
- out(ielem,2) = j+1+top;
- out(ielem,3) = j +top;
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
+ // -- lower-left element
+ out(ielem,0) = bot + (3*ix );
+ out(ielem,1) = bot + (3*ix+1);
+ out(ielem,2) = mid + (2*ix );
+ out(ielem,3) = top + ( ix );
+ ielem++;
+ // -- lower-center element
+ out(ielem,0) = bot + (3*ix+1);
+ out(ielem,1) = bot + (3*ix+2);
+ out(ielem,2) = mid + (2*ix+1);
+ out(ielem,3) = mid + (2*ix );
+ ielem++;
+ // -- lower-right element
+ out(ielem,0) = bot + (3*ix+2);
+ out(ielem,1) = bot + (3*ix+3);
+ out(ielem,2) = top + ( ix+1);
+ out(ielem,3) = mid + (2*ix+1);
+ ielem++;
+ // -- upper element
+ out(ielem,0) = mid + (2*ix );
+ out(ielem,1) = mid + (2*ix+1);
+ out(ielem,2) = top + ( ix+1);
+ out(ielem,3) = top + ( ix );
ielem++;
}
}
}
return out;
}
// ------------------------------ element numbers of the middle layer ------------------------------
inline ColS FineLayer::elementsMiddleLayer()
{
- ColS out(m_nx);
+ // number of element layers in y-direction, the index of the middle layer
+ size_t nely = static_cast(m_nhy.size());
+ size_t iy = (nely-1)/2;
- size_t j = m_startElem((m_startElem.size()-1)/2);
+ ColS out(m_nelx(iy));
- for ( size_t i = 0 ; i < m_nx ; ++i ) out(i) = i+j;
+ for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix )
+ out(ix) = m_startElem(iy) + ix;
return out;
}
// ------------------------------ node-numbers along the bottom edge -------------------------------
inline ColS FineLayer::nodesBottomEdge()
{
- ColS nodes;
-
- if ( m_nh(0)%2 == 0 ) nodes.conservativeResize(m_nx/(3*m_nh(0)/2)+1);
- else nodes.conservativeResize(m_nx/ m_nh(0) +1);
+ ColS out(m_nelx(0)+1);
- for ( size_t i = 0 ; i < static_cast(nodes.size()) ; ++i ) nodes(i) = i;
+ for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix )
+ out(ix) = m_startNode(0) + ix;
- return nodes;
+ return out;
}
// -------------------------------- node-numbers along the top edge --------------------------------
inline ColS FineLayer::nodesTopEdge()
{
- ColS nodes;
+ size_t nely = static_cast(m_nhy.size());
- if ( m_nh(0)%2 == 0 ) nodes.conservativeResize(m_nx/(3*m_nh(0)/2)+1);
- else nodes.conservativeResize(m_nx/ m_nh(0) +1);
+ ColS out(m_nelx(nely-1)+1);
- size_t j = m_startNode(m_startNode.size()-1);
+ for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix )
+ out(ix) = m_startNode(nely) + ix;
- for ( size_t i = 0 ; i < static_cast(nodes.size()) ; ++i ) nodes(i) = i+j;
-
- return nodes;
+ return out;
}
-// ----------------------------- node-numbers along the node left edge -----------------------------
+// ------------------------------- node-numbers along the left edge --------------------------------
inline ColS FineLayer::nodesLeftEdge()
{
- return m_startNode;
+ size_t nely = static_cast(m_nhy.size());
+
+ ColS out(nely+1);
+
+ for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
+ out(iy) = m_startNode(iy);
+
+ for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
+ out(iy+1) = m_startNode(iy+1);
+
+ return out;
}
// ------------------------------- node-numbers along the right edge -------------------------------
inline ColS FineLayer::nodesRightEdge()
{
- size_t N = static_cast(m_nh.size());
- ColS nodes(N+1);
+ size_t nely = static_cast(m_nhy.size());
- // bottom half
- for ( size_t i = 0 ; i < (N-1)/2 ; ++i )
- {
- if ( m_nh(i)%2 == 0 ) nodes(i) = m_startNode(i) + m_nx/(3*m_nh(i)/2);
- else nodes(i) = m_startNode(i) + m_nx/ m_nh(i) ;
- }
+ ColS out(nely+1);
- // top half
- for ( size_t i = (N-1)/2 ; i < N ; ++i )
- {
- if ( m_nh(i)%2 == 0 ) nodes(i) = m_startNode(i) + m_nx/(m_nh(i)/2);
- else nodes(i) = m_startNode(i) + m_nx/ m_nh(i) ;
- }
+ for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
+ out(iy) = m_startNode(iy) + m_nelx(iy);
+
+ for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
+ out(iy+1) = m_startNode(iy+1) + m_nelx(iy);
+
+ return out;
+}
+
+// ---------------------- node-numbers along the bottom edge, without corners ----------------------
- // last node
- nodes(N) = m_nnode-1;
+inline ColS FineLayer::nodesBottomOpenEdge()
+{
+ ColS out(m_nelx(0)-1);
+
+ for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix )
+ out(ix-1) = m_startNode(0) + ix;
+
+ return out;
+}
+
+// ----------------------- node-numbers along the top edge, without corners ------------------------
+
+inline ColS FineLayer::nodesTopOpenEdge()
+{
+ size_t nely = static_cast(m_nhy.size());
+
+ ColS out(m_nelx(nely-1)-1);
+
+ for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix )
+ out(ix-1) = m_startNode(nely) + ix;
+
+ return out;
+}
+
+// ----------------------- node-numbers along the left edge, without corners -----------------------
+
+inline ColS FineLayer::nodesLeftOpenEdge()
+{
+ size_t nely = static_cast(m_nhy.size());
+
+ ColS out(nely-1);
+
+ for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
+ out(iy-1) = m_startNode(iy);
+
+ for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
+ out(iy) = m_startNode(iy+1);
- return nodes;
+ return out;
+}
+
+// ---------------------- node-numbers along the right edge, without corners -----------------------
+
+inline ColS FineLayer::nodesRightOpenEdge()
+{
+ size_t nely = static_cast(m_nhy.size());
+
+ ColS out(nely-1);
+
+ for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
+ out(iy-1) = m_startNode(iy) + m_nelx(iy);
+
+ for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
+ out(iy) = m_startNode(iy+1) + m_nelx(iy);
+
+ return out;
}
// ----------------------------- node-number of the bottom-left corner -----------------------------
inline size_t FineLayer::nodesBottomLeftCorner()
{
- ColS nodes = nodesBottomEdge();
-
- return nodes(0);
+ return m_startNode(0);
}
// ---------------------------- node-number of the bottom-right corner -----------------------------
inline size_t FineLayer::nodesBottomRightCorner()
{
- ColS nodes = nodesBottomEdge();
-
- return nodes(nodes.size()-1);
+ return m_startNode(0) + m_nelx(0);
}
// ------------------------------ node-number of the top-left corner -------------------------------
inline size_t FineLayer::nodesTopLeftCorner()
{
- ColS nodes = nodesTopEdge();
+ size_t nely = static_cast(m_nhy.size());
- return nodes(0);
+ return m_startNode(nely);
}
// ------------------------------ node-number of the top-right corner ------------------------------
inline size_t FineLayer::nodesTopRightCorner()
{
- ColS nodes = nodesTopEdge();
+ size_t nely = static_cast(m_nhy.size());
- return nodes(nodes.size()-1);
+ return m_startNode(nely) + m_nelx(nely-1);
}
-// ----------------------------- node-number of the corners (aliases) ------------------------------
+// -------------------------------------------- aliases --------------------------------------------
inline size_t FineLayer::nodesLeftBottomCorner() { return nodesBottomLeftCorner(); }
-inline size_t FineLayer::nodesLeftTopCorner() { return nodesTopLeftCorner(); }
inline size_t FineLayer::nodesRightBottomCorner() { return nodesBottomRightCorner(); }
+inline size_t FineLayer::nodesLeftTopCorner() { return nodesTopLeftCorner(); }
inline size_t FineLayer::nodesRightTopCorner() { return nodesTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline MatS FineLayer::nodesPeriodic()
{
- // edges
- ColS bot = nodesBottomEdge();
- ColS top = nodesTopEdge();
- ColS rgt = nodesRightEdge();
- ColS lft = nodesLeftEdge();
+ // edges (without corners)
+ ColS bot = nodesBottomOpenEdge();
+ ColS top = nodesTopOpenEdge();
+ ColS lft = nodesLeftOpenEdge();
+ ColS rgt = nodesRightOpenEdge();
// allocate nodal ties
- MatS nodes(bot.size()-2 + lft.size()-2 + 3, 2);
+ // - number of tying per category
+ size_t tedge = bot.size() + lft.size();
+ size_t tnode = 3;
+ // - allocate
+ MatS out(tedge+tnode, 2);
// counter
size_t i = 0;
// tie all corners to one corner
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesBottomRightCorner(); ++i;
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopRightCorner(); ++i;
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopLeftCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i;
- // tie edges to each other (exclude corners)
- for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; }
- for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; }
+ // tie all corresponding edges to each other
+ for ( auto j = 0 ; j(nodePer.rows());
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for ( size_t i = 0 ; i < nper ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
out(nodePer(i,1),j) = out(nodePer(i,0),j);
// renumber "out" to be sequential
return GooseFEM::Mesh::renumber(out);
}
// =================================================================================================
}}} // namespace ...
// =================================================================================================
#endif
-
-
diff --git a/src/GooseFEM/MeshQuad4.h b/src/GooseFEM/MeshQuad4.h
index 35a776f..5960947 100644
--- a/src/GooseFEM/MeshQuad4.h
+++ b/src/GooseFEM/MeshQuad4.h
@@ -1,122 +1,136 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHQUAD4_H
#define GOOSEFEM_MESHQUAD4_H
// -------------------------------------------------------------------------------------------------
#include "GooseFEM.h"
// ===================================== GooseFEM::Mesh::Quad4 =====================================
namespace GooseFEM {
namespace Mesh {
namespace Quad4 {
// ========================================== REGULAR MESH =========================================
class Regular
{
private:
- size_t m_nx; // number of 'pixels' horizontal direction (length == "m_nx * m_h")
- size_t m_ny; // number of 'pixels' vertical direction (length == "m_ny * m_h")
- double m_h; // size of the element edge (equal in both directions)
+ double m_h; // elementary element edge-size (in both directions)
+ size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h")
+ size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h")
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
size_t m_nne=4; // number of nodes-per-element
size_t m_ndim=2; // number of dimensions
public:
- // mesh with "nx*ny" 'pixels' and edge size "h"
- Regular(size_t nx, size_t ny, double h=1.);
+ // mesh with "nelx*nely" 'elements' of edge size "h"
+ Regular(size_t nelx, size_t nely, double h=1.);
// sizes
size_t nelem(); // number of elements
size_t nnode(); // number of nodes
size_t nne(); // number of nodes-per-element
size_t ndim(); // number of dimensions
// mesh
- MatD coor(); // nodal positions [ nnode , ndim ]
- MatS conn(); // connectivity [ nelem , nne ]
+ MatD coor(); // nodal positions [nnode ,ndim]
+ MatS conn(); // connectivity [nelem ,nne ]
// boundary nodes: edges
- ColS nodesBottomEdge(); // nodes along the bottom edge
- ColS nodesTopEdge(); // nodes along the top edge
- ColS nodesLeftEdge(); // nodes along the left edge
- ColS nodesRightEdge(); // nodes along the right edge
+ ColS nodesBottomEdge(); // node-numbers along the bottom edge
+ ColS nodesTopEdge(); // node-numbers along the top edge
+ ColS nodesLeftEdge(); // node-numbers along the left edge
+ ColS nodesRightEdge(); // node-numbers along the right edge
+ // boundary nodes: edges, without corners
+ ColS nodesBottomOpenEdge(); // node-numbers along the bottom edge
+ ColS nodesTopOpenEdge(); // node-numbers along the top edge
+ ColS nodesLeftOpenEdge(); // node-numbers along the left edge
+ ColS nodesRightOpenEdge(); // node-numbers along the right edge
// boundary nodes: corners
- size_t nodesBottomLeftCorner(); // bottom - left corner node
- size_t nodesBottomRightCorner(); // bottom - right corner node
- size_t nodesTopLeftCorner(); // top - left corner node
- size_t nodesTopRightCorner(); // top - right corner node
+ size_t nodesBottomLeftCorner(); // node-number of the bottom - left corner
+ size_t nodesBottomRightCorner(); // node-number of the bottom - right corner
+ size_t nodesTopLeftCorner(); // node-number of the top - left corner
+ size_t nodesTopRightCorner(); // node-number of the top - right corner
// boundary nodes: corners (aliases)
- size_t nodesLeftBottomCorner(); // bottom - left corner node
- size_t nodesLeftTopCorner(); // top - left corner node
- size_t nodesRightBottomCorner(); // bottom - right corner node
- size_t nodesRightTopCorner(); // top - right corner node
+ size_t nodesLeftBottomCorner(); // alias, see above: nodesBottomLeftCorner
+ size_t nodesLeftTopCorner(); // alias, see above: nodesBottomRightCorner
+ size_t nodesRightBottomCorner(); // alias, see above: nodesTopLeftCorner
+ size_t nodesRightTopCorner(); // alias, see above: nodesTopRightCorner
// periodicity
MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent)
- size_t nodesOrigin(); // bottom-left node, to be used as reference for periodicity
+ size_t nodesOrigin(); // bottom-left node, used as reference for periodicity
MatS dofs(); // DOF-numbers for each component of each node (sequential)
- MatS dofsPeriodic(); // DOF-numbers for each component of each node (sequential)
+ MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated
};
// ====================== MESH WITH A FINE LAYER THAT EXPONENTIALLY COARSENS =======================
class FineLayer
{
private:
- double m_h; // elementary element size (middle-layer in x-direction has "L = m_nx * m_h")
- size_t m_nx; // number of elements in x-direction
- ColS m_nh; // element size in y-direction of each layer
- ColS m_startNode; // start node of each layer
- ColS m_startElem; // start element of each layer
+ double m_h; // elementary element edge-size (in all directions)
+ double m_Lx; // mesh size in "x"
+ ColS m_nelx; // number of elements in "x" (per el.layer in "y")
+ ColS m_nnd; // total number of nodes in the main node layer (per nd.layer in "y")
+ ColS m_nhx, m_nhy; // element size in each direction (per el.layer in "y")
+ ColI m_refine; // refine direction (-1:no refine, 0:"x") (per el.layer in "y")
+ ColS m_startElem; // start element (per el.layer in "y")
+ ColS m_startNode; // start node (per nd.layer in "y")
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
size_t m_nne=4; // number of nodes-per-element
size_t m_ndim=2; // number of dimensions
public:
- // mesh with "nx*ny" 'pixels' and edge size "h"; the elements in y-direction are coarsened
- FineLayer(size_t nx, size_t ny, double h=1., size_t nfine=1, size_t nskip=0);
+ // mesh with "nelx*nely" elements of edge size "h"; elements are coarsened in "y"-direction
+ FineLayer(size_t nelx, size_t nely, double h=1., size_t nfine=1);
// sizes
size_t nelem(); // number of elements
size_t nnode(); // number of nodes
size_t nne(); // number of nodes-per-element
size_t ndim(); // number of dimensions
- size_t shape(size_t i); // actual shape in horizontal and vertical direction
+ size_t shape(size_t i); // actual shape in a certain direction
// mesh
- MatD coor(); // nodal positions [ nnode , ndim ]
- MatS conn(); // connectivity [ nelem , nne ]
- // boundary nodes: edges
+ MatD coor(); // nodal positions [nnode ,ndim]
+ MatS conn(); // connectivity [nelem ,nne ]
+ // element sets
ColS elementsMiddleLayer(); // elements in the middle, fine, layer
- ColS nodesBottomEdge(); // nodes along the bottom edge
- ColS nodesTopEdge(); // nodes along the top edge
- ColS nodesLeftEdge(); // nodes along the left edge
- ColS nodesRightEdge(); // nodes along the right edge
+ // boundary nodes: edges
+ ColS nodesBottomEdge(); // node-numbers along the bottom edge
+ ColS nodesTopEdge(); // node-numbers along the top edge
+ ColS nodesLeftEdge(); // node-numbers along the left edge
+ ColS nodesRightEdge(); // node-numbers along the right edge
+ // boundary nodes: edges, without corners
+ ColS nodesBottomOpenEdge(); // node-numbers along the bottom edge
+ ColS nodesTopOpenEdge(); // node-numbers along the top edge
+ ColS nodesLeftOpenEdge(); // node-numbers along the left edge
+ ColS nodesRightOpenEdge(); // node-numbers along the right edge
// boundary nodes: corners
- size_t nodesBottomLeftCorner(); // bottom - left corner node
- size_t nodesBottomRightCorner(); // bottom - right corner node
- size_t nodesTopLeftCorner(); // top - left corner node
- size_t nodesTopRightCorner(); // top - right corner node
+ size_t nodesBottomLeftCorner(); // node-number of the bottom - left corner
+ size_t nodesBottomRightCorner(); // node-number of the bottom - right corner
+ size_t nodesTopLeftCorner(); // node-number of the top - left corner
+ size_t nodesTopRightCorner(); // node-number of the top - right corner
// boundary nodes: corners (aliases)
- size_t nodesLeftBottomCorner(); // bottom - left corner node
- size_t nodesLeftTopCorner(); // top - left corner node
- size_t nodesRightBottomCorner(); // bottom - right corner node
- size_t nodesRightTopCorner(); // top - right corner node
+ size_t nodesLeftBottomCorner(); // alias, see above: nodesBottomLeftCorner
+ size_t nodesLeftTopCorner(); // alias, see above: nodesBottomRightCorner
+ size_t nodesRightBottomCorner(); // alias, see above: nodesTopLeftCorner
+ size_t nodesRightTopCorner(); // alias, see above: nodesTopRightCorner
// periodicity
MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent)
- size_t nodesOrigin(); // bottom-left node, to be used as reference for periodicity
+ size_t nodesOrigin(); // bottom-left node, used as reference for periodicity
MatS dofs(); // DOF-numbers for each component of each node (sequential)
- MatS dofsPeriodic(); // DOF-numbers for each component of each node (sequential)
+ MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated
};
// =================================================================================================
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/src/GooseFEM/MeshTri3.cpp b/src/GooseFEM/MeshTri3.cpp
index 4206e23..bb0840d 100644
--- a/src/GooseFEM/MeshTri3.cpp
+++ b/src/GooseFEM/MeshTri3.cpp
@@ -1,553 +1,606 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHTRI3_CPP
#define GOOSEFEM_MESHTRI3_CPP
// -------------------------------------------------------------------------------------------------
#include "MeshTri3.h"
// ===================================== GooseFEM::Mesh::Tri3 =====================================
namespace GooseFEM {
namespace Mesh {
namespace Tri3 {
// ===================================== CLASS - REGULAR MESH ======================================
// ------------------------------------------ constructor ------------------------------------------
-inline Regular::Regular(size_t nx, size_t ny, double h): m_nx(nx), m_ny(ny), m_h(h)
+inline Regular::Regular(size_t nelx, size_t nely, double h):
+m_h(h), m_nelx(nelx), m_nely(nely)
{
- assert( m_nx >= 1 );
- assert( m_ny >= 1 );
+ assert( m_nelx >= 1 );
+ assert( m_nely >= 1 );
- m_nnode = (m_nx+1) * (m_ny+1) ;
- m_nelem = m_nx * m_ny * 2;
+ m_nnode = (m_nelx+1) * (m_nely+1) ;
+ m_nelem = m_nelx * m_nely * 2;
}
// -------------------------------------- number of elements ---------------------------------------
inline size_t Regular::nelem()
{
return m_nelem;
}
// ---------------------------------------- number of nodes ----------------------------------------
inline size_t Regular::nnode()
{
return m_nnode;
}
// ---------------------------------- number of nodes per element ----------------------------------
inline size_t Regular::nne()
{
return m_nne;
}
// ------------------------------------- number of dimensions --------------------------------------
inline size_t Regular::ndim()
{
return m_ndim;
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline MatD Regular::coor()
{
- MatD coor( m_nnode , m_ndim );
+ MatD out(m_nnode , m_ndim);
- ColD x = ColD::LinSpaced( m_nx+1 , 0.0 , m_h * static_cast(m_nx) );
- ColD y = ColD::LinSpaced( m_ny+1 , 0.0 , m_h * static_cast(m_ny) );
+ ColD x = ColD::LinSpaced( m_nelx+1 , 0.0 , m_h * static_cast(m_nelx) );
+ ColD y = ColD::LinSpaced( m_nely+1 , 0.0 , m_h * static_cast(m_nely) );
size_t inode = 0;
- for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) {
- for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) {
- coor(inode,0) = x(ix);
- coor(inode,1) = y(iy);
+ for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) {
+ for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) {
+ out(inode,0) = x(ix);
+ out(inode,1) = y(iy);
++inode;
}
}
- return coor;
+ return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline MatS Regular::conn()
{
- MatS conn( m_nelem , m_nne );
+ MatS out(m_nelem , m_nne);
size_t ielem = 0;
- for ( size_t iy = 0 ; iy < m_ny ; ++iy ) {
- for ( size_t ix = 0 ; ix < m_nx ; ++ix ) {
- conn(ielem,0) = (iy+0)*(m_nx+1) + (ix+0);
- conn(ielem,1) = (iy+0)*(m_nx+1) + (ix+1);
- conn(ielem,2) = (iy+1)*(m_nx+1) + (ix+0);
+ for ( size_t iy = 0 ; iy < m_nely ; ++iy ) {
+ for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) {
+ out(ielem,0) = (iy )*(m_nelx+1) + (ix );
+ out(ielem,1) = (iy )*(m_nelx+1) + (ix+1);
+ out(ielem,2) = (iy+1)*(m_nelx+1) + (ix );
++ielem;
- conn(ielem,0) = (iy+0)*(m_nx+1) + (ix+1);
- conn(ielem,1) = (iy+1)*(m_nx+1) + (ix+1);
- conn(ielem,2) = (iy+1)*(m_nx+1) + (ix+0);
+ out(ielem,0) = (iy )*(m_nelx+1) + (ix+1);
+ out(ielem,1) = (iy+1)*(m_nelx+1) + (ix+1);
+ out(ielem,2) = (iy+1)*(m_nelx+1) + (ix );
++ielem;
}
}
- return conn;
+ return out;
}
// ------------------------------ node-numbers along the bottom edge -------------------------------
inline ColS Regular::nodesBottomEdge()
{
- ColS nodes(m_nx+1);
+ ColS out(m_nelx+1);
- for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix )
- nodes(ix) = ix;
+ for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
+ out(ix) = ix;
- return nodes;
+ return out;
}
// -------------------------------- node-numbers along the top edge --------------------------------
inline ColS Regular::nodesTopEdge()
{
- ColS nodes(m_nx+1);
+ ColS out(m_nelx+1);
- for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix )
- nodes(ix) = ix + m_ny*(m_nx+1);
+ for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
+ out(ix) = ix + m_nely*(m_nelx+1);
- return nodes;
+ return out;
}
-// ----------------------------- node-numbers along the node left edge -----------------------------
+// ------------------------------- node-numbers along the left edge --------------------------------
inline ColS Regular::nodesLeftEdge()
{
- ColS nodes(m_ny+1);
+ ColS out(m_nely+1);
- for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy )
- nodes(iy) = iy*(m_nx+1);
+ for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
+ out(iy) = iy*(m_nelx+1);
- return nodes;
+ return out;
}
// ------------------------------- node-numbers along the right edge -------------------------------
inline ColS Regular::nodesRightEdge()
{
- ColS nodes(m_ny+1);
+ ColS out(m_nely+1);
+
+ for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
+ out(iy) = iy*(m_nelx+1) + m_nelx;
+
+ return out;
+}
+
+// ---------------------- node-numbers along the bottom edge, without corners ----------------------
+
+inline ColS Regular::nodesBottomOpenEdge()
+{
+ ColS out(m_nelx-1);
- for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy )
- nodes(iy) = iy*(m_nx+1) + m_nx;
+ for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
+ out(ix-1) = ix;
- return nodes;
+ return out;
+}
+
+// ----------------------- node-numbers along the top edge, without corners ------------------------
+
+inline ColS Regular::nodesTopOpenEdge()
+{
+ ColS out(m_nelx-1);
+
+ for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
+ out(ix-1) = ix + m_nely*(m_nelx+1);
+
+ return out;
+}
+
+// ----------------------- node-numbers along the left edge, without corners -----------------------
+
+inline ColS Regular::nodesLeftOpenEdge()
+{
+ ColS out(m_nely-1);
+
+ for ( size_t iy = 1 ; iy < m_nely ; ++iy )
+ out(iy-1) = iy*(m_nelx+1);
+
+ return out;
+}
+
+// ---------------------- node-numbers along the right edge, without corners -----------------------
+
+inline ColS Regular::nodesRightOpenEdge()
+{
+ ColS out(m_nely-1);
+
+ for ( size_t iy = 1 ; iy < m_nely ; ++iy )
+ out(iy-1) = iy*(m_nelx+1) + m_nelx;
+
+ return out;
}
// ----------------------------- node-number of the bottom-left corner -----------------------------
inline size_t Regular::nodesBottomLeftCorner()
{
return 0;
}
// ---------------------------- node-number of the bottom-right corner -----------------------------
inline size_t Regular::nodesBottomRightCorner()
{
- return m_nx;
+ return m_nelx;
}
// ------------------------------ node-number of the top-left corner -------------------------------
inline size_t Regular::nodesTopLeftCorner()
{
- return m_ny*(m_nx+1);
+ return m_nely*(m_nelx+1);
}
// ------------------------------ node-number of the top-right corner ------------------------------
inline size_t Regular::nodesTopRightCorner()
{
- return (m_ny+1)*(m_nx+1) - 1;
+ return m_nely*(m_nelx+1) + m_nelx;
}
// ----------------------------- node-number of the corners (aliases) ------------------------------
inline size_t Regular::nodesLeftBottomCorner() { return nodesBottomLeftCorner(); }
inline size_t Regular::nodesLeftTopCorner() { return nodesTopLeftCorner(); }
inline size_t Regular::nodesRightBottomCorner() { return nodesBottomRightCorner(); }
inline size_t Regular::nodesRightTopCorner() { return nodesTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline MatS Regular::nodesPeriodic()
{
- // edges
- ColS bot = nodesBottomEdge();
- ColS top = nodesTopEdge();
- ColS rgt = nodesRightEdge();
- ColS lft = nodesLeftEdge();
+ // edges (without corners)
+ ColS bot = nodesBottomOpenEdge();
+ ColS top = nodesTopOpenEdge();
+ ColS lft = nodesLeftOpenEdge();
+ ColS rgt = nodesRightOpenEdge();
// allocate nodal ties
- MatS nodes(bot.size()-2 + lft.size()-2 + 3, 2);
+ // - number of tying per category
+ size_t tedge = bot.size() + lft.size();
+ size_t tnode = 3;
+ // - allocate
+ MatS out(tedge+tnode, 2);
// counter
size_t i = 0;
// tie all corners to one corner
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesBottomRightCorner(); ++i;
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopRightCorner(); ++i;
- nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopLeftCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i;
+ out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i;
- // tie edges to each other (exclude corners)
- for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; }
- for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; }
+ // tie all corresponding edges to each other
+ for ( auto j = 0 ; j(nodePer.rows());
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for ( size_t i = 0 ; i < nper ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
out(nodePer(i,1),j) = out(nodePer(i,0),j);
// renumber "out" to be sequential
return GooseFEM::Mesh::renumber(out);
}
// ========================================= MESH ANALYSIS =========================================
// ------------------------------ get the orientation of each element ------------------------------
inline ColI getOrientation(const MatD &coor, const MatS &conn)
{
assert( conn.cols() == 3 );
assert( coor.cols() == 2 );
Eigen::Vector2d v1,v2;
double k;
size_t nelem = static_cast( conn.rows() );
ColI out( nelem );
for ( size_t ielem = 0 ; ielem < nelem ; ++ielem )
{
v1 = coor.row( conn(ielem,0) ) - coor.row( conn(ielem,1) );
v2 = coor.row( conn(ielem,2) ) - coor.row( conn(ielem,1) );
k = v1(0) * v2(1) - v2(0) * v1(1);
if ( k < 0 ) out( ielem ) = -1;
else out( ielem ) = +1;
}
return out;
}
// ------------------------------ set the orientation of each element ------------------------------
inline MatS setOrientation(const MatD &coor, const MatS &conn, int orientation)
{
assert( conn.cols() == 3 );
assert( coor.cols() == 2 );
assert( orientation == -1 || orientation == +1 );
ColI val = getOrientation(coor, conn);
return setOrientation( coor, conn, val, orientation );
}
// -------------------- set the orientation of each element to a certain value ---------------------
inline MatS setOrientation(const MatD &coor, const MatS &conn, const ColI &val, int orientation)
{
assert( conn.cols() == 3 );
assert( coor.cols() == 2 );
assert( conn.rows() == val.size() );
assert( orientation == -1 || orientation == +1 );
// avoid compiler warning
UNUSED(coor);
size_t nelem = static_cast(conn.rows());
MatS out = conn;
for ( size_t ielem = 0 ; ielem < nelem ; ++ielem )
if ( ( orientation == -1 and val(ielem) > 0 ) or ( orientation == +1 and val(ielem) < 0 ) )
std::swap( out(ielem,2) , out(ielem,1) );
return out;
}
// ===================================== RE-TRIANGULATION API ======================================
inline MatS retriangulate(const MatD &coor, const MatS &conn, int orientation)
{
// get the orientation of all elements
ColI dir = getOrientation( coor, conn );
// check the orientation
bool eq = std::abs( dir.sum() ) == conn.rows();
// new connectivity
MatS out;
// perform re-triangulation
// - use "TriUpdate"
if ( eq )
{
TriUpdate obj(coor,conn);
obj.eval();
out = obj.conn();
}
// - using TriCompute
else
{
throw std::runtime_error("Work-in-progress, has to be re-triangulated using 'TriCompute'");
}
return setOrientation(coor,out,orientation);
}
// ================ RE-TRIANGULATION SUPPORT CLASS - UPDATE EXISTING TRIANGULATION =================
// ------------------------------------------ constructor ------------------------------------------
inline TriUpdate::TriUpdate(const MatD &coor, const MatS &conn): m_conn(conn), m_coor(coor)
{
assert( conn.cols() == 3 );
assert( coor.cols() == 2 );
// store shapes
m_nnode = static_cast( coor.rows() );
m_ndim = static_cast( coor.cols() );
m_nelem = static_cast( conn.rows() );
m_nne = static_cast( conn.cols() );
// resize internal arrays
m_elem.resize(2);
m_node.resize(4);
// set default to out-of-bounds, to make clear that nothing happened yet
m_elem.setConstant(m_nelem);
m_node.setConstant(m_nnode);
edge();
}
// -------------------------- compute neighbors per edge of all elements ---------------------------
inline void TriUpdate::edge()
{
m_edge.resize( m_nelem , m_nne );
m_edge.setConstant( m_nelem ); // signal that nothing has been set
std::vector idx = {0,1,2}; // lists to convert connectivity -> edges
std::vector jdx = {1,2,0}; // lists to convert connectivity -> edges
std::vector edge;
edge.reserve( m_nelem*idx.size() );
for ( size_t e = 0 ; e < m_nelem ; ++e )
for ( size_t i = 0 ; i < m_nne ; ++i )
edge.push_back( Edge( m_conn(e,idx[i]), m_conn(e,jdx[i]) , e , i , true ) );
std::sort( edge.begin() , edge.end() , Edge_sort );
for ( size_t i = 0 ; i < edge.size()-1 ; ++i )
{
if ( edge[i].n1 == edge[i+1].n1 and edge[i].n2 == edge[i+1].n2 )
{
m_edge( edge[i ].elem , edge[i ].edge ) = edge[i+1].elem;
m_edge( edge[i+1].elem , edge[i+1].edge ) = edge[i ].elem;
}
}
}
// ---------------------------- update edges around renumbered elements ----------------------------
inline void TriUpdate::chedge(size_t edge, size_t old_elem, size_t new_elem)
{
size_t m;
size_t neigh = m_edge( old_elem , edge );
if ( neigh >= m_nelem ) return;
for ( m = 0 ; m < m_nne ; ++m )
if ( m_edge( neigh , m ) == old_elem )
break;
m_edge( neigh , m ) = new_elem;
}
// --------------------------------- re-triangulate the full mesh ----------------------------------
inline bool TriUpdate::eval()
{
bool change = false;
while ( increment() ) { change = true; }
return change;
}
// ----------------------------------- one re-triangulation step -----------------------------------
inline bool TriUpdate::increment()
{
size_t ielem,jelem,iedge,jedge;
double phi1,phi2;
ColS c(4);
ColS n(4);
// loop over all elements
for ( ielem = 0 ; ielem < m_nelem ; ++ielem )
{
// loop over all edges
for ( iedge = 0 ; iedge < m_nne ; ++iedge )
{
// only proceed if the edge is shared with another element
if ( m_edge(ielem,iedge) >= m_nelem )
continue;
// read "jelem"
jelem = m_edge(ielem,iedge);
// find the edge involved for "jelem"
for ( jedge=0; jedge M_PI )
{
// update connectivity
m_conn(ielem,0) = c(0); m_conn(ielem,1) = c(3); m_conn(ielem,2) = c(2);
m_conn(jelem,0) = c(0); m_conn(jelem,1) = c(1); m_conn(jelem,2) = c(3);
// change list with neighbors for the elements around (only two neighbors change)
if ( iedge==0 ) { chedge(2,ielem,jelem); }
else if ( iedge==1 ) { chedge(0,ielem,jelem); }
else if ( iedge==2 ) { chedge(1,ielem,jelem); }
if ( jedge==0 ) { chedge(2,jelem,ielem); }
else if ( jedge==1 ) { chedge(0,jelem,ielem); }
else if ( jedge==2 ) { chedge(1,jelem,ielem); }
// convert to four static nodes
if ( iedge==0 ) { n(0)=m_edge(ielem,2); n(3)=m_edge(ielem,1); }
else if ( iedge==1 ) { n(0)=m_edge(ielem,0); n(3)=m_edge(ielem,2); }
else if ( iedge==2 ) { n(0)=m_edge(ielem,1); n(3)=m_edge(ielem,0); }
if ( jedge==0 ) { n(1)=m_edge(jelem,1); n(2)=m_edge(jelem,2); }
else if ( jedge==1 ) { n(1)=m_edge(jelem,2); n(2)=m_edge(jelem,0); }
else if ( jedge==2 ) { n(1)=m_edge(jelem,0); n(2)=m_edge(jelem,1); }
// store the neighbors for the changed elements
m_edge(ielem,0) = jelem; m_edge(jelem,0) = n(0) ;
m_edge(ielem,1) = n(2) ; m_edge(jelem,1) = n(1) ;
m_edge(ielem,2) = n(3) ; m_edge(jelem,2) = ielem;
// store information for transfer algorithm
m_node = c;
m_elem(0) = ielem;
m_elem(1) = jelem;
return true;
}
}
}
return false;
}
// ================================ SUPPORT CLASS - EDGE DEFINITION ================================
// ------------------------------------------ constructor ------------------------------------------
Edge::Edge(size_t i, size_t j, size_t el, size_t ed, bool sort):
n1(i), n2(j), elem(el), edge(ed)
{
if ( sort && n1>n2 )
std::swap(n1,n2);
}
// --------------------------------------- compare two edges ---------------------------------------
inline bool Edge_cmp(Edge a, Edge b)
{
if ( a.n1 == b.n1 and a.n2 == b.n2 )
return true;
return false;
}
// ----------------------- sort edges by comparing the first and second node -----------------------
inline bool Edge_sort(Edge a, Edge b)
{
if ( a.n1 < b.n1 or a.n2 < b.n2 )
return true;
return false;
}
// =================================================================================================
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/src/GooseFEM/MeshTri3.h b/src/GooseFEM/MeshTri3.h
index 245d8e0..baffb16 100644
--- a/src/GooseFEM/MeshTri3.h
+++ b/src/GooseFEM/MeshTri3.h
@@ -1,141 +1,146 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHTRI3_H
#define GOOSEFEM_MESHTRI3_H
// -------------------------------------------------------------------------------------------------
#include "Mesh.h"
// ===================================== GooseFEM::Mesh::Tri3 =====================================
namespace GooseFEM {
namespace Mesh {
namespace Tri3 {
// ========================================== REGULAR MESH =========================================
class Regular
{
private:
- size_t m_nx; // number of 'pixels' horizontal direction (length == "m_nx * m_h")
- size_t m_ny; // number of 'pixels' vertical direction (length == "m_ny * m_h")
- double m_h; // size of the element edge (equal in both directions)
- size_t m_nelem; // number of elements
- size_t m_nnode; // number of nodes
- size_t m_nne=3; // number of nodes-per-element
- size_t m_ndim=2; // number of dimensions
+ double m_h; // elementary element edge-size (in both directions)
+ size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h")
+ size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h")
+ size_t m_nelem; // number of elements
+ size_t m_nnode; // number of nodes
+ size_t m_nne=3; // number of nodes-per-element
+ size_t m_ndim=2; // number of dimensions
public:
- // mesh with "nx" pixels in horizontal direction, "ny" in vertical direction and "h" the edge size
- Regular(size_t nx, size_t ny, double h=1.);
+ // mesh with "2*nelx*nely" 'elements' of edge size "h"
+ Regular(size_t nelx, size_t nely, double h=1.);
// sizes
- size_t nelem(); // number of elements
- size_t nnode(); // number of nodes
- size_t nne(); // number of nodes-per-element
- size_t ndim(); // number of dimensions
+ size_t nelem(); // number of elements
+ size_t nnode(); // number of nodes
+ size_t nne(); // number of nodes-per-element
+ size_t ndim(); // number of dimensions
// mesh
- MatD coor(); // nodal positions [ nnode , ndim ]
- MatS conn(); // connectivity [ nelem , nne ]
+ MatD coor(); // nodal positions [nnode ,ndim]
+ MatS conn(); // connectivity [nelem ,nne ]
// boundary nodes: edges
- ColS nodesBottomEdge(); // nodes along the bottom edge
- ColS nodesTopEdge(); // nodes along the top edge
- ColS nodesLeftEdge(); // nodes along the left edge
- ColS nodesRightEdge(); // nodes along the right edge
+ ColS nodesBottomEdge(); // node-numbers along the bottom edge
+ ColS nodesTopEdge(); // node-numbers along the top edge
+ ColS nodesLeftEdge(); // node-numbers along the left edge
+ ColS nodesRightEdge(); // node-numbers along the right edge
+ // boundary nodes: edges, without corners
+ ColS nodesBottomOpenEdge(); // node-numbers along the bottom edge
+ ColS nodesTopOpenEdge(); // node-numbers along the top edge
+ ColS nodesLeftOpenEdge(); // node-numbers along the left edge
+ ColS nodesRightOpenEdge(); // node-numbers along the right edge
// boundary nodes: corners
- size_t nodesBottomLeftCorner(); // bottom - left corner node
- size_t nodesBottomRightCorner(); // bottom - right corner node
- size_t nodesTopLeftCorner(); // top - left corner node
- size_t nodesTopRightCorner(); // top - right corner node
+ size_t nodesBottomLeftCorner(); // node-number of the bottom - left corner
+ size_t nodesBottomRightCorner(); // node-number of the bottom - right corner
+ size_t nodesTopLeftCorner(); // node-number of the top - left corner
+ size_t nodesTopRightCorner(); // node-number of the top - right corner
// boundary nodes: corners (aliases)
- size_t nodesLeftBottomCorner(); // bottom - left corner node
- size_t nodesLeftTopCorner(); // top - left corner node
- size_t nodesRightBottomCorner(); // bottom - right corner node
- size_t nodesRightTopCorner(); // top - right corner node
+ size_t nodesLeftBottomCorner(); // alias, see above: nodesBottomLeftCorner
+ size_t nodesLeftTopCorner(); // alias, see above: nodesBottomRightCorner
+ size_t nodesRightBottomCorner(); // alias, see above: nodesTopLeftCorner
+ size_t nodesRightTopCorner(); // alias, see above: nodesTopRightCorner
// periodicity
- MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent)
- size_t nodesOrigin(); // bottom-left node, to be used as reference for periodicity
- MatS dofs(); // DOF-numbers for each component of each node (sequential)
- MatS dofsPeriodic(); // DOF-numbers for each component of each node (sequential)
+ MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent)
+ size_t nodesOrigin(); // bottom-left node, used as reference for periodicity
+ MatS dofs(); // DOF-numbers for each component of each node (sequential)
+ MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated
};
// ========================================= MESH ANALYSIS =========================================
// read / set the orientation (-1 / +1) of all triangles
inline ColI getOrientation(const MatD &coor, const MatS &conn );
inline MatS setOrientation(const MatD &coor, const MatS &conn, int orientation=-1);
inline MatS setOrientation(const MatD &coor, const MatS &conn, const ColI &val, int orientation=-1);
// ======================================= RE-TRIANGULATION ========================================
// simple interface to compute the full re-triangulation; it uses, depending on the input mesh:
// (1) the minimal evasive "TriUpdate"
// (2) the more rigorous "TriCompute"
inline MatS retriangulate(const MatD &coor, const MatS &conn, int orientation=-1);
// ------------------------- support class - update existing triangulation -------------------------
// minimal evasive re-triangulation which only flips edges of the existing connectivity
class TriUpdate
{
private:
MatS m_edge; // the element that neighbors along each edge (m_nelem: no neighbor)
MatS m_conn; // connectivity (updated)
MatD m_coor; // nodal positions (does not change)
size_t m_nelem; // #elements
size_t m_nnode; // #nodes
size_t m_nne; // #nodes-per-element
size_t m_ndim; // #dimensions
ColS m_elem; // the two elements involved in the last element change (see below)
ColS m_node; // the four nodes involved in the last element change (see below)
// old: m_elem(0) = [ m_node(0) , m_node(1) , m_node(2) ]
// m_elem(1) = [ m_node(1) , m_node(3) , m_node(2) ]
// new: m_elem(0) = [ m_node(0) , m_node(3) , m_node(2) ]
// m_elem(1) = [ m_node(0) , m_node(1) , m_node(3) ]
// compute neighbors per edge of all elements
void edge();
// update edges around renumbered elements
void chedge(size_t edge, size_t old_elem, size_t new_elem);
public:
TriUpdate(){};
TriUpdate(const MatD &coor, const MatS &conn);
bool eval (); // re-triangulate the full mesh (returns "true" if changed)
bool increment(); // one re-triangulation step (returns "true" if changed)
MatS conn () { return m_conn; }; // return (new) connectivity
MatS ch_elem () { return m_elem; }; // return element involved in last element change
MatS ch_node () { return m_node; }; // return nodes involved in last element change
};
// ================================ SUPPORT CLASS - EDGE DEFINITION ================================
// support class to allow the storage of a list of edges
class Edge {
public:
size_t n1 ; // node 1 (edge from node 1-2)
size_t n2 ; // node 2 (edge from node 1-2)
size_t elem; // element to which the edge belong
size_t edge; // edge index within the element (e.g. edge==1 -> n1=conn(0,elem), n2=conn(1,elem))
Edge(){};
Edge(size_t i, size_t j, size_t el, size_t ed, bool sort=false);
};
// -------------------------------------------------------------------------------------------------
// compare edges
inline bool Edge_cmp (Edge a, Edge b); // check equality
inline bool Edge_sort(Edge a, Edge b); // check if "a" is smaller than "b" in terms of node-numbers
// =================================================================================================
}}} // namespace ...
#endif