diff --git a/docs/details/figures/MeshTri3/ReTri/ReTri.py b/docs/details/figures/MeshTri3/ReTri/ReTri.py
deleted file mode 100644
index 7fa4d06..0000000
--- a/docs/details/figures/MeshTri3/ReTri/ReTri.py
+++ /dev/null
@@ -1,131 +0,0 @@
-
-import numpy as np
-import GooseFEM
-
-# define a mesh, and extract the relevant data
-mesh = GooseFEM.Mesh.Tri3.Regular(10,10)
-coor = mesh.coor()
-conn = mesh.conn()
-nnode = coor.shape[0]
-nelem = conn.shape[0]
-per = mesh.nodesPeriodic()
-
-# define a random perturbation of the nodal coordinates
-delta = np.zeros(coor.shape)
-delta[ 5,0] = .08
-delta[ 5,1] = .08
-delta[35,0] = -.08
-delta[35,1] = -.08
-delta[55,0] = .08
-delta[55,1] = -.08
-# make periodic
-delta[per[:,1],0] = delta[per[:,0],0]
-delta[per[:,1],1] = delta[per[:,0],1]
-
-# copy the coordinates and connectivity (for plotting)
-newcoor = np.array(coor,copy=True)
-newconn = np.array(conn,copy=True)
-
-# add the perturbation incrementally, compute the new connectivity each time
-newcoor += delta*.5
-newconn = GooseFEM.Mesh.Tri3.retriangulate(newcoor,newconn)
-
-newcoor += delta*.01
-newconn = GooseFEM.Mesh.Tri3.retriangulate(newcoor,newconn)
-
-newcoor += delta*.01
-newconn = GooseFEM.Mesh.Tri3.retriangulate(newcoor,newconn)
-
-newcoor += delta*.4
-newconn = GooseFEM.Mesh.Tri3.retriangulate(newcoor,newconn)
-
-# ==================================================================================================
-
-# define colors to clearly show the difference (and element overlaps)
-orig = np.zeros(( nelem ))
-diff = np.zeros(( nelem ))
-# find the differences
-for i in range(nelem):
- if not np.all( conn[i,:] == newconn[i,:] ):
- diff[i] = 1
-
-# define shifts to show periodicity
-DX = np.zeros( coor.shape )
-DY = np.zeros( coor.shape )
-DX[:,0] = +1.
-DY[:,1] = +1.
-
-# --------------------------------------------------------------------------------------------------
-
-import matplotlib.pylab as plt
-import goosempl as gplt
-
-plt.style.use(['goose','goose-latex'])
-
-fig = plt.figure(figsize=(18,6))
-fig.set_tight_layout(True)
-
-# --------------------------------------------------------------------------------------------------
-
-ax = fig.add_subplot(1,3,1)
-
-gplt.patch(coor=coor,conn=conn,cindex=orig,clim=[-0.5,1.5],cmap='RdBu_r',alpha=.5)
-
-gplt.patch(coor=coor-DX ,conn=conn)
-gplt.patch(coor=coor+DX ,conn=conn)
-gplt.patch(coor=coor+DX-DY,conn=conn)
-gplt.patch(coor=coor+DX+DY,conn=conn)
-gplt.patch(coor=coor-DX-DY,conn=conn)
-gplt.patch(coor=coor-DX+DY,conn=conn)
-gplt.patch(coor=coor -DY,conn=conn)
-gplt.patch(coor=coor +DY,conn=conn)
-
-plt.xlim([ -0.15 , 1.2 ])
-plt.ylim([ -0.15 , 1.2 ])
-
-ax.set_aspect('equal')
-
-# --------------------------------------------------------------------------------------------------
-
-ax = fig.add_subplot(1,3,2)
-
-gplt.patch(coor=newcoor,conn=conn,cindex=orig,clim=[-0.5,1.5],cmap='RdBu_r',alpha=.5)
-
-gplt.patch(coor=newcoor-DX ,conn=conn)
-gplt.patch(coor=newcoor+DX ,conn=conn)
-gplt.patch(coor=newcoor+DX-DY,conn=conn)
-gplt.patch(coor=newcoor+DX+DY,conn=conn)
-gplt.patch(coor=newcoor-DX-DY,conn=conn)
-gplt.patch(coor=newcoor-DX+DY,conn=conn)
-gplt.patch(coor=newcoor -DY,conn=conn)
-gplt.patch(coor=newcoor +DY,conn=conn)
-
-plt.xlim([ -0.15 , 1.2 ])
-plt.ylim([ -0.15 , 1.2 ])
-
-ax.set_aspect('equal')
-
-# --------------------------------------------------------------------------------------------------
-
-ax = fig.add_subplot(1,3,3)
-
-gplt.patch(coor=newcoor,conn=newconn,cindex=diff,clim=[-0.5,1.5],cmap='RdBu_r',alpha=.5)
-
-gplt.patch(coor=newcoor-DX ,conn=newconn)
-gplt.patch(coor=newcoor+DX ,conn=newconn)
-gplt.patch(coor=newcoor+DX-DY,conn=newconn)
-gplt.patch(coor=newcoor+DX+DY,conn=newconn)
-gplt.patch(coor=newcoor-DX-DY,conn=newconn)
-gplt.patch(coor=newcoor-DX+DY,conn=newconn)
-gplt.patch(coor=newcoor -DY,conn=newconn)
-gplt.patch(coor=newcoor +DY,conn=newconn)
-
-plt.xlim([ -0.15 , 1.2 ])
-plt.ylim([ -0.15 , 1.2 ])
-
-ax.set_aspect('equal')
-
-# --------------------------------------------------------------------------------------------------
-
-plt.savefig('ReTri.svg')
-plt.show()
diff --git a/docs/details/figures/MeshTri3/ReTri/ReTri.svg b/docs/details/figures/MeshTri3/ReTri/ReTri.svg
deleted file mode 100644
index 9ef4353..0000000
--- a/docs/details/figures/MeshTri3/ReTri/ReTri.svg
+++ /dev/null
@@ -1,28197 +0,0 @@
-
-
-
-
diff --git a/docs/details/figures/MeshTri3/ReTri/ReTri_random.py b/docs/details/figures/MeshTri3/ReTri/ReTri_random.py
deleted file mode 100644
index d2c8f46..0000000
--- a/docs/details/figures/MeshTri3/ReTri/ReTri_random.py
+++ /dev/null
@@ -1,119 +0,0 @@
-
-import numpy as np
-import GooseFEM
-
-# define a mesh, and extract the relevant data
-mesh = GooseFEM.Mesh.Tri3.Regular(10,10)
-coor = mesh.coor()
-conn = mesh.conn()
-nnode = coor.shape[0]
-nelem = conn.shape[0]
-per = mesh.nodesPeriodic()
-
-# define a random perturbation of the nodal coordinates
-delta = np.random.random(coor.shape)
-delta[:,0] *= .1
-delta[:,1] *= .1
-# make periodic
-delta[per[:,1],0] = delta[per[:,0],0]
-delta[per[:,1],1] = delta[per[:,0],1]
-
-# copy the coordinates and connectivity (for plotting)
-newcoor = np.array(coor,copy=True)
-newconn = np.array(conn,copy=True)
-
-# add the perturbation incrementally, compute the new connectivity each time
-for i in range(4):
- newcoor += delta/4.
- newconn = GooseFEM.Mesh.Tri3.retriangulate(newcoor,newconn)
-
-# ==================================================================================================
-
-# define colors to clearly show the difference (and element overlaps)
-orig = np.zeros(( nelem ))
-diff = np.zeros(( nelem ))
-# find the differences
-for i in range(nelem):
- if not np.all( conn[i,:] == newconn[i,:] ):
- diff[i] = 1
-
-# define shifts to show periodicity
-DX = np.zeros( coor.shape )
-DY = np.zeros( coor.shape )
-DX[:,0] = +1.
-DY[:,1] = +1.
-
-# --------------------------------------------------------------------------------------------------
-
-import matplotlib.pylab as plt
-import goosempl as gplt
-
-plt.style.use(['goose','goose-latex'])
-
-fig = plt.figure(figsize=(18,6))
-fig.set_tight_layout(True)
-
-# --------------------------------------------------------------------------------------------------
-
-ax = fig.add_subplot(1,3,1)
-
-gplt.patch(coor=coor,conn=conn,cindex=orig,clim=[-0.5,1.5],cmap='RdBu_r',alpha=.5)
-
-gplt.patch(coor=coor-DX ,conn=conn)
-gplt.patch(coor=coor+DX ,conn=conn)
-gplt.patch(coor=coor+DX-DY,conn=conn)
-gplt.patch(coor=coor+DX+DY,conn=conn)
-gplt.patch(coor=coor-DX-DY,conn=conn)
-gplt.patch(coor=coor-DX+DY,conn=conn)
-gplt.patch(coor=coor -DY,conn=conn)
-gplt.patch(coor=coor +DY,conn=conn)
-
-plt.xlim([ -0.15 , 1.2 ])
-plt.ylim([ -0.15 , 1.2 ])
-
-ax.set_aspect('equal')
-
-# --------------------------------------------------------------------------------------------------
-
-ax = fig.add_subplot(1,3,2)
-
-gplt.patch(coor=newcoor,conn=conn,cindex=orig,clim=[-0.5,1.5],cmap='RdBu_r',alpha=.5)
-
-gplt.patch(coor=newcoor-DX ,conn=conn)
-gplt.patch(coor=newcoor+DX ,conn=conn)
-gplt.patch(coor=newcoor+DX-DY,conn=conn)
-gplt.patch(coor=newcoor+DX+DY,conn=conn)
-gplt.patch(coor=newcoor-DX-DY,conn=conn)
-gplt.patch(coor=newcoor-DX+DY,conn=conn)
-gplt.patch(coor=newcoor -DY,conn=conn)
-gplt.patch(coor=newcoor +DY,conn=conn)
-
-plt.xlim([ -0.15 , 1.2 ])
-plt.ylim([ -0.15 , 1.2 ])
-
-ax.set_aspect('equal')
-
-# --------------------------------------------------------------------------------------------------
-
-ax = fig.add_subplot(1,3,3)
-
-gplt.patch(coor=newcoor,conn=newconn,cindex=diff,clim=[-0.5,1.5],cmap='RdBu_r',alpha=.5)
-
-gplt.patch(coor=newcoor-DX ,conn=newconn)
-gplt.patch(coor=newcoor+DX ,conn=newconn)
-gplt.patch(coor=newcoor+DX-DY,conn=newconn)
-gplt.patch(coor=newcoor+DX+DY,conn=newconn)
-gplt.patch(coor=newcoor-DX-DY,conn=newconn)
-gplt.patch(coor=newcoor-DX+DY,conn=newconn)
-gplt.patch(coor=newcoor -DY,conn=newconn)
-gplt.patch(coor=newcoor +DY,conn=newconn)
-
-plt.xlim([ -0.15 , 1.2 ])
-plt.ylim([ -0.15 , 1.2 ])
-
-ax.set_aspect('equal')
-
-# --------------------------------------------------------------------------------------------------
-
-plt.savefig('ReTri_random.svg')
-plt.show()
diff --git a/docs/details/figures/MeshTri3/ReTri/ReTri_random.svg b/docs/details/figures/MeshTri3/ReTri/ReTri_random.svg
deleted file mode 100644
index 14284a4..0000000
--- a/docs/details/figures/MeshTri3/ReTri/ReTri_random.svg
+++ /dev/null
@@ -1,28197 +0,0 @@
-
-
-
-
diff --git a/docs/details/figures/MeshTri3/ReTri/element_change.svg b/docs/details/figures/MeshTri3/ReTri/element_change.svg
deleted file mode 100644
index ddd21f0..0000000
--- a/docs/details/figures/MeshTri3/ReTri/element_change.svg
+++ /dev/null
@@ -1,167 +0,0 @@
-
-
-
-
diff --git a/include/GooseFEM/MeshTri3.h b/include/GooseFEM/MeshTri3.h
index 2a2014f..5388889 100644
--- a/include/GooseFEM/MeshTri3.h
+++ b/include/GooseFEM/MeshTri3.h
@@ -1,162 +1,101 @@
/* =================================================================================================
(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 "config.h"
// =================================================================================================
namespace GooseFEM {
namespace Mesh {
namespace Tri3 {
// -------------------------------------------------------------------------------------------------
class Regular
{
private:
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
static const size_t m_nne=3; // number of nodes-per-element
static const size_t m_ndim=2; // number of dimensions
public:
// mesh with "2*nelx*nely" 'elements' of edge size "h"
Regular(size_t nelx, size_t nely, double h=1.);
// sizes
size_t nelem() const; // number of elements
size_t nnode() const; // number of nodes
size_t nne() const; // number of nodes-per-element
size_t ndim() const; // number of dimensions
// type
ElementType getElementType() const;
// mesh
xt::xtensor coor() const; // nodal positions [nnode ,ndim]
xt::xtensor conn() const; // connectivity [nelem ,nne ]
// boundary nodes: edges
xt::xtensor nodesBottomEdge() const; // node-numbers along the bottom edge
xt::xtensor nodesTopEdge() const; // node-numbers along the top edge
xt::xtensor nodesLeftEdge() const; // node-numbers along the left edge
xt::xtensor nodesRightEdge() const; // node-numbers along the right edge
// boundary nodes: edges, without corners
xt::xtensor nodesBottomOpenEdge() const; // node-numbers along the bottom edge
xt::xtensor nodesTopOpenEdge() const; // node-numbers along the top edge
xt::xtensor nodesLeftOpenEdge() const; // node-numbers along the left edge
xt::xtensor nodesRightOpenEdge() const; // node-numbers along the right edge
// boundary nodes: corners
size_t nodesBottomLeftCorner() const; // node-number of the bottom - left corner
size_t nodesBottomRightCorner() const; // node-number of the bottom - right corner
size_t nodesTopLeftCorner() const; // node-number of the top - left corner
size_t nodesTopRightCorner() const; // node-number of the top - right corner
// boundary nodes: corners (aliases)
size_t nodesLeftBottomCorner() const; // alias, see above: nodesBottomLeftCorner
size_t nodesLeftTopCorner() const; // alias, see above: nodesBottomRightCorner
size_t nodesRightBottomCorner() const; // alias, see above: nodesTopLeftCorner
size_t nodesRightTopCorner() const; // alias, see above: nodesTopRightCorner
// periodicity
xt::xtensor nodesPeriodic() const; // periodic node pairs [:,2]: (independent, dependent)
size_t nodesOrigin() const; // bottom-left node, used as reference for periodicity
xt::xtensor dofs() const; // DOF-numbers for each component of each node (sequential)
xt::xtensor dofsPeriodic() const; // ,, for the case that the periodicity if fully eliminated
};
// -------------------------------------------------------------------------------------------------
-// read / set the orientation (-1 / +1) of all triangles
-inline xt::xtensor getOrientation(const xt::xtensor &coor, const xt::xtensor &conn );
-inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, int orientation=-1);
-inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, const xt::xtensor &val, int orientation=-1);
-
-// -------------------------------------------------------------------------------------------------
-
-// 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 xt::xtensor retriangulate(const xt::xtensor &coor, const xt::xtensor &conn, int orientation=-1);
-
-// =================================================================================================
-
-namespace Private {
-
-// -------------------------------------------------------------------------------------------------
-
-// minimal evasive re-triangulation which only flips edges of the existing connectivity
-class TriUpdate
-{
-private:
- size_t m_nelem; // #elements
- size_t m_nnode; // #nodes
- size_t m_nne; // #nodes-per-element
- size_t m_ndim; // #dimensions
- xt::xtensor m_edge; // the element that neighbors along each edge (m_nelem: no neighbor)
- xt::xtensor m_conn; // connectivity (updated)
- xt::xtensor m_coor; // nodal positions (does not change)
- xt::xtensor_fixed> m_elem; // the two elements in the last element change
- xt::xtensor_fixed> m_node; // the four nodes in the last element change
- // 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() = default;
- TriUpdate(const xt::xtensor &coor, const xt::xtensor &conn);
-
- bool eval (); // re-triangulate the full mesh (returns "true" if changed)
- bool increment(); // one re-triangulation step (returns "true" if changed)
- xt::xtensor conn() { return m_conn; }; // (new) connectivity
- xt::xtensor ch_elem() { return m_elem; }; // element involved in last element change
- xt::xtensor ch_node() { return m_node; }; // nodes involved in last element change
-};
-
-// -------------------------------------------------------------------------------------------------
-
-// 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() = default;
- 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 Private
+// read / set the orientation (-1/+1) of all triangles
+inline xt::xtensor getOrientation(
+ const xt::xtensor &coor, // nodal coordinates
+ const xt::xtensor &conn); // connectivity
+
+inline xt::xtensor setOrientation(
+ const xt::xtensor &coor, // nodal coordinates
+ const xt::xtensor &conn, // connectivity
+ int orientation=-1); // target orientation (-1/+1)
+
+inline xt::xtensor setOrientation(
+ const xt::xtensor &coor, // nodal coordinates
+ const xt::xtensor &conn, // connectivity
+ const xt::xtensor &val, // current orientations (== output of "getOrientation")
+ int orientation=-1); // target orientation (-1/+1)
// =================================================================================================
}}} // namespace ...
// =================================================================================================
#include "MeshTri3.hpp"
// =================================================================================================
#endif
diff --git a/include/GooseFEM/MeshTri3.hpp b/include/GooseFEM/MeshTri3.hpp
index 00d260f..e0c9277 100644
--- a/include/GooseFEM/MeshTri3.hpp
+++ b/include/GooseFEM/MeshTri3.hpp
@@ -1,614 +1,392 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHTRI3_HPP
#define GOOSEFEM_MESHTRI3_HPP
// -------------------------------------------------------------------------------------------------
#include "MeshTri3.h"
// =================================================================================================
namespace GooseFEM {
namespace Mesh {
namespace Tri3 {
// -------------------------------------------------------------------------------------------------
inline Regular::Regular(size_t nelx, size_t nely, double h):
m_h(h), m_nelx(nelx), m_nely(nely)
{
GOOSEFEM_ASSERT(m_nelx >= 1ul);
GOOSEFEM_ASSERT(m_nely >= 1ul);
m_nnode = (m_nelx+1) * (m_nely+1);
m_nelem = m_nelx * m_nely * 2;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nelem() const
{
return m_nelem;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nnode() const
{
return m_nnode;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nne() const
{
return m_nne;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::ndim() const
{
return m_ndim;
}
// -------------------------------------------------------------------------------------------------
inline ElementType Regular::getElementType() const
{
return ElementType::Tri3;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::coor() const
{
xt::xtensor out = xt::empty({m_nnode, m_ndim});
xt::xtensor x = xt::linspace(0.0, m_h*static_cast(m_nelx), m_nelx+1);
xt::xtensor y = xt::linspace(0.0, m_h*static_cast(m_nely), m_nely+1);
size_t inode = 0;
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 out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::conn() const
{
xt::xtensor out = xt::empty({m_nelem,m_nne});
size_t ielem = 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;
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 out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesBottomEdge() const
{
xt::xtensor out = xt::empty({m_nelx+1});
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = ix;
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesTopEdge() const
{
xt::xtensor out = xt::empty({m_nelx+1});
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = ix + m_nely*(m_nelx+1);
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesLeftEdge() const
{
xt::xtensor out = xt::empty({m_nely+1});
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1);
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesRightEdge() const
{
xt::xtensor out = xt::empty({m_nely+1});
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1) + m_nelx;
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesBottomOpenEdge() const
{
xt::xtensor out = xt::empty({m_nelx-1});
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = ix;
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesTopOpenEdge() const
{
xt::xtensor out = xt::empty({m_nelx-1});
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = ix + m_nely*(m_nelx+1);
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesLeftOpenEdge() const
{
xt::xtensor out = xt::empty({m_nely-1});
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1);
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesRightOpenEdge() const
{
xt::xtensor out = xt::empty({m_nely-1});
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelx;
return out;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nodesBottomLeftCorner() const
{
return 0;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nodesBottomRightCorner() const
{
return m_nelx;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nodesTopLeftCorner() const
{
return m_nely*(m_nelx+1);
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nodesTopRightCorner() const
{
return m_nely*(m_nelx+1) + m_nelx;
}
// -------------------------------------------------------------------------------------------------
inline size_t Regular::nodesLeftBottomCorner() const
{ return nodesBottomLeftCorner(); }
inline size_t Regular::nodesLeftTopCorner() const
{ return nodesTopLeftCorner(); }
inline size_t Regular::nodesRightBottomCorner() const
{ return nodesBottomRightCorner(); }
inline size_t Regular::nodesRightTopCorner() const
{ return nodesTopRightCorner(); }
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::nodesPeriodic() const
{
// edges (without corners)
xt::xtensor bot = nodesBottomOpenEdge();
xt::xtensor top = nodesTopOpenEdge();
xt::xtensor lft = nodesLeftOpenEdge();
xt::xtensor rgt = nodesRightOpenEdge();
// allocate nodal ties
// - number of tying per category
size_t tedge = bot.size() + lft.size();
size_t tnode = 3;
// - allocate
xt::xtensor out = xt::empty({tedge+tnode, std::size_t(2)});
// counter
size_t i = 0;
// tie all corners to one corner
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 all corresponding edges to each other
for ( size_t j = 0 ; j Regular::dofs() const
{
return GooseFEM::Mesh::dofs(m_nnode,m_ndim);
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor Regular::dofsPeriodic() const
{
// DOF-numbers for each component of each node (sequential)
xt::xtensor out = GooseFEM::Mesh::dofs(m_nnode,m_ndim);
// periodic node-pairs
xt::xtensor nodePer = nodesPeriodic();
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for ( size_t i = 0 ; i < nodePer.shape()[0] ; ++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);
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor getOrientation(const xt::xtensor &coor, const xt::xtensor &conn)
+inline xt::xtensor getOrientation(
+ const xt::xtensor &coor,
+ const xt::xtensor &conn)
{
GOOSEFEM_ASSERT(conn.shape()[1] == 3ul);
GOOSEFEM_ASSERT(coor.shape()[1] == 2ul);
double k;
size_t nelem = conn.shape()[0];
xt::xtensor out = xt::empty({nelem});
for ( size_t ielem = 0 ; ielem < nelem ; ++ielem )
{
auto v1 = xt::view(coor, conn(ielem,0), xt::all()) - xt::view(coor, conn(ielem,1), xt::all());
auto v2 = xt::view(coor, conn(ielem,2), xt::all()) - xt::view(coor, conn(ielem,1), xt::all());
k = v1(0) * v2(1) - v2(0) * v1(1);
if ( k < 0 ) out(ielem) = -1;
else out(ielem) = +1;
}
return out;
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, int orientation)
+inline xt::xtensor setOrientation(
+ const xt::xtensor &coor,
+ const xt::xtensor &conn,
+ int orientation)
{
GOOSEFEM_ASSERT(conn.shape()[1] == 3ul);
GOOSEFEM_ASSERT(coor.shape()[1] == 2ul);
GOOSEFEM_ASSERT(orientation == -1 || orientation == +1);
xt::xtensor val = getOrientation(coor, conn);
return setOrientation(coor, conn, val, orientation);
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, const xt::xtensor &val, int orientation)
+inline xt::xtensor setOrientation(
+ const xt::xtensor &coor,
+ const xt::xtensor &conn,
+ const xt::xtensor &val,
+ int orientation)
{
GOOSEFEM_ASSERT(conn.shape()[1] == 3ul);
GOOSEFEM_ASSERT(coor.shape()[1] == 2ul);
GOOSEFEM_ASSERT(conn.shape()[0] == val.size());
GOOSEFEM_ASSERT(orientation == -1 || orientation == +1);
// avoid compiler warning
UNUSED(coor);
size_t nelem = conn.shape()[0];
xt::xtensor 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) );
+ for (size_t ielem = 0; ielem < nelem; ++ielem)
+ if ((orientation == -1 && val(ielem) > 0) || (orientation == +1 && val(ielem) < 0))
+ std::swap(out(ielem,2) , out(ielem,1));
return out;
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor retriangulate(const xt::xtensor &coor, const xt::xtensor &conn, int orientation)
-{
- // get the orientation of all elements
- xt::xtensor dir = getOrientation(coor, conn);
- // check the orientation
- bool eq = static_cast(std::abs(xt::sum(dir)[0])) == conn.shape()[0];
-
- // new connectivity
- xt::xtensor out;
-
- // perform re-triangulation
- // - use "TriUpdate"
- if ( eq )
- {
- Private::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);
-}
-
-// =================================================================================================
-
-namespace Private {
-
-// -------------------------------------------------------------------------------------------------
-
-inline TriUpdate::TriUpdate(const xt::xtensor &coor, const xt::xtensor &conn): m_conn(conn), m_coor(coor)
-{
- GOOSEFEM_ASSERT(conn.shape()[1] == 3ul);
- GOOSEFEM_ASSERT(coor.shape()[1] == 2ul);
-
- // store shapes
- m_nnode = coor.shape()[0];
- m_ndim = coor.shape()[1];
- m_nelem = conn.shape()[0];
- m_nne = conn.shape()[1];
-
- // set default to out-of-bounds, to make clear that nothing happened yet
- m_elem = m_nelem * xt::ones({2});
- m_node = m_nnode * xt::ones({4});
-
- edge();
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline void TriUpdate::edge()
-{
- // signal that nothing has been set
- m_edge = m_nelem * xt::ones({m_nelem , m_nne});
-
- 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;
- }
- }
-}
-
-// -------------------------------------------------------------------------------------------------
-
-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;
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline bool TriUpdate::eval()
-{
- bool change = false;
-
- while ( increment() ) { change = true; }
-
- return change;
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline bool TriUpdate::increment()
-{
- size_t ielem,jelem,iedge,jedge;
- double phi1,phi2;
-
- xt::xtensor_fixed> c = xt::empty({4});
- xt::xtensor_fixed> n = xt::empty({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;
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline 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);
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline bool Edge_cmp(Edge a, Edge b)
-{
- if ( a.n1 == b.n1 and a.n2 == b.n2 )
- return true;
-
- return false;
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline bool Edge_sort(Edge a, Edge b)
-{
- if ( a.n1 < b.n1 or a.n2 < b.n2 )
- return true;
-
- return false;
-}
-
-// -------------------------------------------------------------------------------------------------
-
-} // namespace Private
-
-// -------------------------------------------------------------------------------------------------
-
}}} // namespace ...
// =================================================================================================
#endif