diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h
index 5f99eff..9927c2a 100644
--- a/include/GooseFEM/Mesh.h
+++ b/include/GooseFEM/Mesh.h
@@ -1,114 +1,128 @@
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_MESH_H
#define GOOSEFEM_MESH_H
#include "config.h"
namespace GooseFEM {
namespace Mesh {
enum class ElementType {
Quad4,
Hex8,
Tri3 };
inline ElementType defaultElementType(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn);
// Renumber to lowest possible index. For example [0,3,4,2] -> [0,2,3,1]
class Renumber {
public:
// constructors
Renumber() = default;
Renumber(const xt::xarray<size_t>& dofs);
// get renumbered DOFs (same as "Renumber::apply(dofs)")
xt::xtensor<size_t, 2> get(const xt::xtensor<size_t, 2>& dofs) const;
// apply renumbering to other set
template <class T>
T apply(const T& list) const;
// get the list needed to renumber, e.g.:
// dofs_renumbered(i,j) = index(dofs(i,j))
xt::xtensor<size_t, 1> index() const;
private:
xt::xtensor<size_t, 1> m_renum;
};
// Reorder to lowest possible index, in specific order.
//
// For example for "Reorder({iiu,iip})" after reordering:
//
// iiu = xt::range<size_t>(nnu);
// iip = xt::range<size_t>(nnp) + nnu;
class Reorder {
public:
// constructors
Reorder() = default;
Reorder(const std::initializer_list<xt::xtensor<size_t, 1>> args);
// get reordered DOFs (same as "Reorder::apply(dofs)")
xt::xtensor<size_t, 2> get(const xt::xtensor<size_t, 2>& dofs) const;
// apply renumbering to other set
template <class T>
T apply(const T& list) const;
// get the list needed to reorder, e.g.:
// dofs_reordered(i,j) = index(dofs(i,j))
xt::xtensor<size_t, 1> index() const;
private:
xt::xtensor<size_t, 1> m_renum;
};
// list with DOF-numbers in sequential order
inline xt::xtensor<size_t, 2> dofs(size_t nnode, size_t ndim);
// renumber to lowest possible index (see "GooseFEM::Mesh::Renumber")
inline xt::xtensor<size_t, 2> renumber(const xt::xtensor<size_t, 2>& dofs);
// number of elements connected to each node
inline xt::xtensor<size_t, 1> coordination(const xt::xtensor<size_t, 2>& conn);
// elements connected to each node
inline std::vector<std::vector<size_t>> elem2node(
const xt::xtensor<size_t, 2>& conn,
bool sorted=true); // ensure the output to be sorted
// return size of each element edge
inline xt::xtensor<double, 2> edgesize(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn,
ElementType type);
inline xt::xtensor<double, 2> edgesize(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn); // extract element-type based on shape of "conn"
// Coordinates of the center of each element
inline xt::xtensor<double, 2> centers(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn,
ElementType type);
inline xt::xtensor<double, 2> centers(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn); // extract element-type based on shape of "conn"
+// Convert an element-map to a node-map.
+// Input: new_elvar = elvar[elem_map]
+// Return: new_nodevar = nodevar[node_map]
+inline xt::xtensor<size_t, 1> elemmap2nodemap(
+ const xt::xtensor<size_t, 1>& elem_map,
+ const xt::xtensor<double, 2>& coor,
+ const xt::xtensor<size_t, 2>& conn); // extract element-type based on shape of "conn"
+
+inline xt::xtensor<size_t, 1> elemmap2nodemap(
+ const xt::xtensor<size_t, 1>& elem_map,
+ const xt::xtensor<double, 2>& coor,
+ const xt::xtensor<size_t, 2>& conn,
+ ElementType type);
+
} // namespace Mesh
} // namespace GooseFEM
#include "Mesh.hpp"
#endif
diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp
index 43c516f..407acaf 100644
--- a/include/GooseFEM/Mesh.hpp
+++ b/include/GooseFEM/Mesh.hpp
@@ -1,247 +1,284 @@
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_MESH_HPP
#define GOOSEFEM_MESH_HPP
#include "Mesh.h"
namespace GooseFEM {
namespace Mesh {
inline ElementType defaultElementType(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn)
{
if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) {
return ElementType::Tri3;
}
if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) {
return ElementType::Quad4;
}
if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) {
return ElementType::Hex8;
}
throw std::runtime_error("Element-type not implemented");
}
inline Renumber::Renumber(const xt::xarray<size_t>& dofs)
{
size_t n = xt::amax(dofs)() + 1;
size_t i = 0;
xt::xtensor<size_t, 1> unique = xt::unique(dofs);
m_renum = xt::empty<size_t>({n});
for (auto& j : unique) {
m_renum(j) = i;
++i;
}
}
// ret(i,j) = renum(list(i,j))
template <class T>
T Renumber::apply(const T& list) const
{
T ret = T::from_shape(list.shape());
auto jt = ret.begin();
for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
*jt = m_renum(*it);
}
return ret;
}
inline xt::xtensor<size_t, 2> Renumber::get(const xt::xtensor<size_t, 2>& dofs) const
{
return this->apply(dofs);
}
inline xt::xtensor<size_t, 1> Renumber::index() const
{
return m_renum;
}
inline Reorder::Reorder(const std::initializer_list<xt::xtensor<size_t, 1>> args)
{
size_t n = 0;
size_t i = 0;
for (auto& arg : args) {
if (arg.size() == 0) {
continue;
}
n = std::max(n, xt::amax(arg)() + 1);
}
#ifdef GOOSEFEM_ENABLE_ASSERT
for (auto& arg : args) {
GOOSEFEM_ASSERT(xt::unique(arg) == xt::sort(arg));
}
#endif
m_renum = xt::empty<size_t>({n});
for (auto& arg : args) {
for (auto& j : arg) {
m_renum(j) = i;
++i;
}
}
}
inline xt::xtensor<size_t, 2> Reorder::get(const xt::xtensor<size_t, 2>& dofs) const
{
return this->apply(dofs);
}
inline xt::xtensor<size_t, 1> Reorder::index() const
{
return m_renum;
}
// apply renumbering, e.g. for a matrix:
//
// ret(i,j) = renum(list(i,j))
template <class T>
T Reorder::apply(const T& list) const
{
T ret = T::from_shape(list.shape());
auto jt = ret.begin();
for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
*jt = m_renum(*it);
}
return ret;
}
inline xt::xtensor<size_t, 2> renumber(const xt::xtensor<size_t, 2>& dofs)
{
return Renumber(dofs).get(dofs);
}
inline xt::xtensor<size_t, 2> dofs(size_t nnode, size_t ndim)
{
return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
}
inline xt::xtensor<size_t, 1> coordination(const xt::xtensor<size_t, 2>& conn)
{
size_t nnode = xt::amax(conn)() + 1;
xt::xtensor<size_t, 1> N = xt::zeros<size_t>({nnode});
for (auto it = conn.begin(); it != conn.end(); ++it) {
N(*it) += 1;
}
return N;
}
inline std::vector<std::vector<size_t>> elem2node(const xt::xtensor<size_t, 2>& conn, bool sorted)
{
auto N = coordination(conn);
auto nnode = N.size();
std::vector<std::vector<size_t>> ret(nnode);
for (size_t i = 0; i < nnode; ++i) {
ret[i].reserve(N(i));
}
for (size_t e = 0; e < conn.shape(0); ++e) {
for (size_t m = 0; m < conn.shape(1); ++m) {
ret[conn(e, m)].push_back(e);
}
}
if (sorted) {
for (auto& row : ret) {
std::sort(row.begin(), row.end());
}
}
return ret;
}
inline xt::xtensor<double, 2> edgesize(
const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn, ElementType type)
{
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
if (type == ElementType::Quad4) {
GOOSEFEM_ASSERT(coor.shape(1) == 2ul);
GOOSEFEM_ASSERT(conn.shape(1) == 4ul);
xt::xtensor<size_t, 1> n0 = xt::view(conn, xt::all(), 0);
xt::xtensor<size_t, 1> n1 = xt::view(conn, xt::all(), 1);
xt::xtensor<size_t, 1> n2 = xt::view(conn, xt::all(), 2);
xt::xtensor<size_t, 1> n3 = xt::view(conn, xt::all(), 3);
xt::xtensor<double, 1> x0 = xt::view(coor, xt::keep(n0), 0);
xt::xtensor<double, 1> x1 = xt::view(coor, xt::keep(n1), 0);
xt::xtensor<double, 1> x2 = xt::view(coor, xt::keep(n2), 0);
xt::xtensor<double, 1> x3 = xt::view(coor, xt::keep(n3), 0);
xt::xtensor<double, 1> y0 = xt::view(coor, xt::keep(n0), 1);
xt::xtensor<double, 1> y1 = xt::view(coor, xt::keep(n1), 1);
xt::xtensor<double, 1> y2 = xt::view(coor, xt::keep(n2), 1);
xt::xtensor<double, 1> y3 = xt::view(coor, xt::keep(n3), 1);
xt::xtensor<double, 2> ret = xt::empty<double>(conn.shape());
xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0));
xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0));
xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0));
xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0));
return ret;
}
throw std::runtime_error("Element-type not implemented");
}
inline xt::xtensor<double, 2> edgesize(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn)
{
return edgesize(coor, conn, defaultElementType(coor, conn));
}
inline xt::xtensor<double, 2> centers(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn,
ElementType type)
{
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
xt::xtensor<double, 2> ret = xt::zeros<double>({conn.shape(0), coor.shape(1)});
if (type == ElementType::Quad4) {
GOOSEFEM_ASSERT(coor.shape(1) == 2);
GOOSEFEM_ASSERT(conn.shape(1) == 4);
for (size_t i = 0; i < 4; ++i) {
auto n = xt::view(conn, xt::all(), i);
ret += xt::view(coor, xt::keep(n), xt::all());
}
ret /= 4.0;
return ret;
}
throw std::runtime_error("Element-type not implemented");
}
inline xt::xtensor<double, 2> centers(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn)
{
return centers(coor, conn, defaultElementType(coor, conn));
}
+inline xt::xtensor<size_t, 1> elemmap2nodemap(
+ const xt::xtensor<size_t, 1>& elem_map,
+ const xt::xtensor<double, 2>& coor,
+ const xt::xtensor<size_t, 2>& conn,
+ ElementType type)
+{
+ GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
+ GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0));
+ size_t N = coor.shape(0);
+
+ xt::xtensor<size_t, 1> ret = N * xt::ones<size_t>({N});
+
+ if (type == ElementType::Quad4) {
+ GOOSEFEM_ASSERT(coor.shape(1) == 2);
+ GOOSEFEM_ASSERT(conn.shape(1) == 4);
+
+ for (size_t i = 0; i < 4; ++i) {
+ xt::xtensor<size_t, 1> t = N * xt::ones<size_t>({N});
+ auto old_nd = xt::view(conn, xt::all(), i);
+ auto new_nd = xt::view(conn, xt::keep(elem_map), i);
+ xt::view(t, xt::keep(old_nd)) = new_nd;
+ ret = xt::where(xt::equal(ret, N), t, ret);
+ }
+
+ return ret;
+ }
+
+ throw std::runtime_error("Element-type not implemented");
+}
+
+inline xt::xtensor<size_t, 1> elemmap2nodemap(
+ const xt::xtensor<size_t, 1>& elem_map,
+ const xt::xtensor<double, 2>& coor,
+ const xt::xtensor<size_t, 2>& conn)
+{
+ return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn));
+}
} // namespace Mesh
} // namespace GooseFEM
#endif
diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h
index 1494415..52a614a 100644
--- a/include/GooseFEM/MeshQuad4.h
+++ b/include/GooseFEM/MeshQuad4.h
@@ -1,258 +1,258 @@
/*
(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 "config.h"
namespace GooseFEM {
namespace Mesh {
namespace Quad4 {
namespace Map {
class FineLayer2Regular;
} // namespace Map
class Regular {
public:
Regular() = default;
Regular(size_t nelx, size_t nely, double h = 1.0);
// size
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
size_t nelx() const; // number of elements in x-direction
size_t nely() const; // number of elements in y-direction
double h() const; // edge size
// type
ElementType getElementType() const;
// mesh
xt::xtensor<double, 2> coor() const; // nodal positions [nnode, ndim]
xt::xtensor<size_t, 2> conn() const; // connectivity [nelem, nne]
// boundary nodes: edges
xt::xtensor<size_t, 1> nodesBottomEdge() const;
xt::xtensor<size_t, 1> nodesTopEdge() const;
xt::xtensor<size_t, 1> nodesLeftEdge() const;
xt::xtensor<size_t, 1> nodesRightEdge() const;
// boundary nodes: edges, without corners
xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
xt::xtensor<size_t, 1> nodesRightOpenEdge() const;
// boundary nodes: corners (including aliases)
size_t nodesBottomLeftCorner() const;
size_t nodesBottomRightCorner() const;
size_t nodesTopLeftCorner() const;
size_t nodesTopRightCorner() const;
size_t nodesLeftBottomCorner() const;
size_t nodesLeftTopCorner() const;
size_t nodesRightBottomCorner() const;
size_t nodesRightTopCorner() const;
// DOF-numbers for each component of each node (sequential)
xt::xtensor<size_t, 2> dofs() const;
// DOF-numbers for the case that the periodicity if fully eliminated
xt::xtensor<size_t, 2> dofsPeriodic() const;
// periodic node pairs [:,2]: (independent, dependent)
xt::xtensor<size_t, 2> nodesPeriodic() const;
// front-bottom-left node, used as reference for periodicity
size_t nodesOrigin() const;
// element numbers as matrix
xt::xtensor<size_t, 2> elementMatrix() const;
private:
double m_h; // elementary element edge-size (in all 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 = 4; // number of nodes-per-element
static const size_t m_ndim = 2; // number of dimensions
};
class FineLayer {
public:
FineLayer() = default;
FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1);
// Reconstruct class for given coordinates / connectivity
FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
// size
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
size_t nelx() const; // number of elements in x-direction
size_t nely() const; // number of elements in y-direction
double h() const; // edge size
// type
ElementType getElementType() const;
// mesh
xt::xtensor<double, 2> coor() const; // nodal positions [nnode, ndim]
xt::xtensor<size_t, 2> conn() const; // connectivity [nelem, nne]
// element sets
xt::xtensor<size_t, 1> elementsMiddleLayer() const; // elements in the middle (fine) layer
// boundary nodes: edges
xt::xtensor<size_t, 1> nodesBottomEdge() const;
xt::xtensor<size_t, 1> nodesTopEdge() const;
xt::xtensor<size_t, 1> nodesLeftEdge() const;
xt::xtensor<size_t, 1> nodesRightEdge() const;
// boundary nodes: edges, without corners
xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
xt::xtensor<size_t, 1> nodesRightOpenEdge() const;
// boundary nodes: corners (including aliases)
size_t nodesBottomLeftCorner() const;
size_t nodesBottomRightCorner() const;
size_t nodesTopLeftCorner() const;
size_t nodesTopRightCorner() const;
size_t nodesLeftBottomCorner() const;
size_t nodesLeftTopCorner() const;
size_t nodesRightBottomCorner() const;
size_t nodesRightTopCorner() const;
// DOF-numbers for each component of each node (sequential)
xt::xtensor<size_t, 2> dofs() const;
// DOF-numbers for the case that the periodicity if fully eliminated
xt::xtensor<size_t, 2> dofsPeriodic() const;
// periodic node pairs [:,2]: (independent, dependent)
xt::xtensor<size_t, 2> nodesPeriodic() const;
// front-bottom-left node, used as reference for periodicity
size_t nodesOrigin() const;
- // mapping to 'roll' periodically,
- // return element mapping, to get the new connectivity use "comm[mesh.roll(), :]"
- xt::xtensor<size_t, 1> roll(size_t n, size_t axis = 0);
+ // mapping to 'roll' periodically in the x-direction,
+ // returns element mapping, such that: new_elemvar = elemvar[elem_map]
+ xt::xtensor<size_t, 1> roll(size_t n);
private:
double m_h; // elementary element edge-size (in all directions)
double m_Lx; // mesh size in "x"
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
static const size_t m_nne = 4; // number of nodes-per-element
static const size_t m_ndim = 2; // number of dimensions
xt::xtensor<size_t, 1> m_nelx; // number of elements in "x" (*)
xt::xtensor<size_t, 1> m_nnd; // total number of nodes in the main node layer (**)
xt::xtensor<size_t, 1> m_nhx; // element size in x-direction (*)
xt::xtensor<size_t, 1> m_nhy; // element size in y-direction (*)
xt::xtensor<int, 1> m_refine; // refine direction (-1:no refine, 0:"x" (*)
xt::xtensor<size_t, 1> m_startElem; // start element (*)
xt::xtensor<size_t, 1> m_startNode; // start node (**)
// (*) per element layer in "y"
// (**) per node layer in "y"
void init(size_t nelx, size_t nely, double h, size_t nfine = 1);
void map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular;
};
namespace Map {
// Return "FineLayer"-class responsible for generating a connectivity
// Throws if conversion is not possible
GooseFEM::Mesh::Quad4::FineLayer FineLayer(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn);
class RefineRegular {
public:
// constructor
RefineRegular() = default;
RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny);
// return the one of the two meshes
GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const;
GooseFEM::Mesh::Quad4::Regular getFineMesh() const;
// elements of the Fine mesh per element of the Coarse mesh
xt::xtensor<size_t, 2> getMap() const;
// map field
xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 1>& data) const; // scalar per el
xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
xt::xtensor<double, 4> mapToCoarse(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
// map field
xt::xtensor<double, 1> mapToFine(const xt::xtensor<double, 1>& data) const; // scalar per el
xt::xtensor<double, 2> mapToFine(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
xt::xtensor<double, 4> mapToFine(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
private:
// the meshes
GooseFEM::Mesh::Quad4::Regular m_coarse;
GooseFEM::Mesh::Quad4::Regular m_fine;
// mapping
xt::xtensor<size_t, 1> m_fine2coarse;
xt::xtensor<size_t, 1> m_fine2coarse_index;
xt::xtensor<size_t, 2> m_coarse2fine;
};
class FineLayer2Regular {
public:
// constructor
FineLayer2Regular() = default;
FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh);
// return either of the meshes
GooseFEM::Mesh::Quad4::Regular getRegularMesh() const;
GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const;
// elements of the Regular mesh per element of the FineLayer mesh
// and the fraction by which the overlap is
std::vector<std::vector<size_t>> getMap() const;
std::vector<std::vector<double>> getMapFraction() const;
// map field
xt::xtensor<double, 1> mapToRegular(const xt::xtensor<double, 1>& data) const; // scalar per el
xt::xtensor<double, 2> mapToRegular(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
xt::xtensor<double, 4> mapToRegular(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
private:
// the "FineLayer" mesh to map
GooseFEM::Mesh::Quad4::FineLayer m_finelayer;
// the new "Regular" mesh to which to map
GooseFEM::Mesh::Quad4::Regular m_regular;
// mapping
std::vector<std::vector<size_t>> m_elem_regular;
std::vector<std::vector<double>> m_frac_regular;
};
} // namespace Map
} // namespace Quad4
} // namespace Mesh
} // namespace GooseFEM
#include "MeshQuad4.hpp"
#endif
diff --git a/include/GooseFEM/MeshQuad4.hpp b/include/GooseFEM/MeshQuad4.hpp
index b6d151b..77557a5 100644
--- a/include/GooseFEM/MeshQuad4.hpp
+++ b/include/GooseFEM/MeshQuad4.hpp
@@ -1,1411 +1,1410 @@
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_MESHQUAD4_HPP
#define GOOSEFEM_MESHQUAD4_HPP
#include "MeshQuad4.h"
namespace GooseFEM {
namespace Mesh {
namespace Quad4 {
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;
}
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 size_t Regular::nelx() const
{
return m_nelx;
}
inline size_t Regular::nely() const
{
return m_nely;
}
inline double Regular::h() const
{
return m_h;
}
inline ElementType Regular::getElementType() const
{
return ElementType::Quad4;
}
inline xt::xtensor<double, 2> Regular::coor() const
{
xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
xt::xtensor<double, 1> x =
xt::linspace<double>(0.0, m_h * static_cast<double>(m_nelx), m_nelx + 1);
xt::xtensor<double, 1> y =
xt::linspace<double>(0.0, m_h * static_cast<double>(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) {
ret(inode, 0) = x(ix);
ret(inode, 1) = y(iy);
++inode;
}
}
return ret;
}
inline xt::xtensor<size_t, 2> Regular::conn() const
{
xt::xtensor<size_t, 2> ret = xt::empty<size_t>({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) {
ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix);
ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1);
ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix);
ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1);
++ielem;
}
}
return ret;
}
inline xt::xtensor<size_t, 1> Regular::nodesBottomEdge() const
{
return xt::arange<size_t>(m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesTopEdge() const
{
return xt::arange<size_t>(m_nelx + 1) + m_nely * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesLeftEdge() const
{
return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesRightEdge() const
{
return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1) + m_nelx;
}
inline xt::xtensor<size_t, 1> Regular::nodesBottomOpenEdge() const
{
return xt::arange<size_t>(1, m_nelx);
}
inline xt::xtensor<size_t, 1> Regular::nodesTopOpenEdge() const
{
return xt::arange<size_t>(1, m_nelx) + m_nely * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesLeftOpenEdge() const
{
return xt::arange<size_t>(1, m_nely) * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesRightOpenEdge() const
{
return xt::arange<size_t>(1, m_nely) * (m_nelx + 1) + m_nelx;
}
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<size_t, 2> Regular::nodesPeriodic() const
{
xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);
ret(0, 0) = nodesBottomLeftCorner();
ret(0, 1) = nodesBottomRightCorner();
ret(1, 0) = nodesBottomLeftCorner();
ret(1, 1) = nodesTopRightCorner();
ret(2, 0) = nodesBottomLeftCorner();
ret(2, 1) = nodesTopLeftCorner();
size_t i = 3;
xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
i += bot.size();
xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
return ret;
}
inline size_t Regular::nodesOrigin() const
{
return nodesBottomLeftCorner();
}
inline xt::xtensor<size_t, 2> Regular::dofs() const
{
return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
}
inline xt::xtensor<size_t, 2> Regular::dofsPeriodic() const
{
xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
for (size_t j = 0; j < m_ndim; ++j) {
xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
}
return GooseFEM::Mesh::renumber(ret);
}
inline xt::xtensor<size_t, 2> Regular::elementMatrix() const
{
return xt::arange<size_t>(m_nelem).reshape({m_nely, m_nelx});
}
inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine)
{
this->init(nelx, nely, h, nfine);
}
inline FineLayer::FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
{
this->map(coor, conn);
}
inline void FineLayer::init(size_t nelx, size_t nely, double h, size_t nfine)
{
GOOSEFEM_ASSERT(nelx >= 1ul);
GOOSEFEM_ASSERT(nely >= 1ul);
m_h = h;
m_Lx = m_h * static_cast<double>(nelx);
// compute element size in y-direction (use symmetry, compute upper half)
// temporary variables
size_t nmin, ntot;
xt::xtensor<size_t, 1> nhx = xt::ones<size_t>({nely});
xt::xtensor<size_t, 1> nhy = xt::ones<size_t>({nely});
xt::xtensor<int, 1> refine = -1 * xt::ones<int>({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;
}
// 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 || ntot >= nmin) {
nely = iy;
break;
}
// rules (1,2) satisfied: coarsen in x-direction
if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) {
refine(iy) = 0;
nhy(iy) *= 2;
auto vnhy = xt::view(nhy, xt::range(iy + 1, _));
auto vnhx = xt::view(nhx, xt::range(iy, _));
vnhy *= 3;
vnhx *= 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 || ntot >= nmin) {
nely = iy;
break;
}
}
// symmetrize, compute full information
// allocate mesh constructor parameters
m_nhx = xt::empty<size_t>({nely * 2 - 1});
m_nhy = xt::empty<size_t>({nely * 2 - 1});
m_refine = xt::empty<int>({nely * 2 - 1});
m_nelx = xt::empty<size_t>({nely * 2 - 1});
m_nnd = xt::empty<size_t>({nely * 2});
m_startElem = xt::empty<size_t>({nely * 2 - 1});
m_startNode = xt::empty<size_t>({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_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_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);
}
// 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;
}
// 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);
}
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 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);
}
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
m_nnode += m_nelx(nely - 1) + 1;
}
inline size_t FineLayer::nelem() const
{
return m_nelem;
}
inline size_t FineLayer::nnode() const
{
return m_nnode;
}
inline size_t FineLayer::nne() const
{
return m_nne;
}
inline size_t FineLayer::ndim() const
{
return m_ndim;
}
inline size_t FineLayer::nelx() const
{
return xt::amax(m_nelx)();
}
inline size_t FineLayer::nely() const
{
return xt::sum(m_nhy)();
}
inline double FineLayer::h() const
{
return m_h;
}
inline ElementType FineLayer::getElementType() const
{
return ElementType::Quad4;
}
inline xt::xtensor<double, 2> FineLayer::coor() const
{
// allocate output
xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
// current node, number of element layers
size_t inode = 0;
size_t nely = static_cast<size_t>(m_nhy.size());
// y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
// - allocate
xt::xtensor<double, 1> y = xt::empty<double>({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
xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);
// add nodes of the bottom layer of this element
for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
ret(inode, 0) = x(ix);
ret(inode, 1) = y(iy);
++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<double>(m_nhx(iy) / 3);
double dy = m_h * static_cast<double>(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) {
ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
ret(inode, 1) = y(iy) + dy;
++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
xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);
// 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<double>(m_nhx(iy) / 3);
double dy = m_h * static_cast<double>(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) {
ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
ret(inode, 1) = y(iy) + dy;
++inode;
}
}
}
// add nodes of the top layer of this element
for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
ret(inode, 0) = x(ix);
ret(inode, 1) = y(iy + 1);
++inode;
}
}
return ret;
}
inline xt::xtensor<size_t, 2> FineLayer::conn() const
{
// allocate output
xt::xtensor<size_t, 2> ret = xt::empty<size_t>({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<size_t>(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 ix = 0; ix < m_nelx(iy); ++ix) {
ret(ielem, 0) = bot + (ix);
ret(ielem, 1) = bot + (ix + 1);
ret(ielem, 2) = top + (ix + 1);
ret(ielem, 3) = top + (ix);
ielem++;
}
}
// - define connectivity: refinement along the x-direction (below the middle layer)
else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
// -- bottom element
ret(ielem, 0) = bot + (ix);
ret(ielem, 1) = bot + (ix + 1);
ret(ielem, 2) = mid + (2 * ix + 1);
ret(ielem, 3) = mid + (2 * ix);
ielem++;
// -- top-right element
ret(ielem, 0) = bot + (ix + 1);
ret(ielem, 1) = top + (3 * ix + 3);
ret(ielem, 2) = top + (3 * ix + 2);
ret(ielem, 3) = mid + (2 * ix + 1);
ielem++;
// -- top-center element
ret(ielem, 0) = mid + (2 * ix);
ret(ielem, 1) = mid + (2 * ix + 1);
ret(ielem, 2) = top + (3 * ix + 2);
ret(ielem, 3) = top + (3 * ix + 1);
ielem++;
// -- top-left element
ret(ielem, 0) = bot + (ix);
ret(ielem, 1) = mid + (2 * ix);
ret(ielem, 2) = top + (3 * ix + 1);
ret(ielem, 3) = top + (3 * ix);
ielem++;
}
}
// - define connectivity: coarsening along the x-direction (above the middle layer)
else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) {
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
// -- lower-left element
ret(ielem, 0) = bot + (3 * ix);
ret(ielem, 1) = bot + (3 * ix + 1);
ret(ielem, 2) = mid + (2 * ix);
ret(ielem, 3) = top + (ix);
ielem++;
// -- lower-center element
ret(ielem, 0) = bot + (3 * ix + 1);
ret(ielem, 1) = bot + (3 * ix + 2);
ret(ielem, 2) = mid + (2 * ix + 1);
ret(ielem, 3) = mid + (2 * ix);
ielem++;
// -- lower-right element
ret(ielem, 0) = bot + (3 * ix + 2);
ret(ielem, 1) = bot + (3 * ix + 3);
ret(ielem, 2) = top + (ix + 1);
ret(ielem, 3) = mid + (2 * ix + 1);
ielem++;
// -- upper element
ret(ielem, 0) = mid + (2 * ix);
ret(ielem, 1) = mid + (2 * ix + 1);
ret(ielem, 2) = top + (ix + 1);
ret(ielem, 3) = top + (ix);
ielem++;
}
}
}
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::elementsMiddleLayer() const
{
size_t nely = m_nhy.size();
size_t iy = (nely - 1) / 2;
return m_startElem(iy) + xt::arange<size_t>(m_nelx(iy));
}
inline xt::xtensor<size_t, 1> FineLayer::nodesBottomEdge() const
{
return m_startNode(0) + xt::arange<size_t>(m_nelx(0) + 1);
}
inline xt::xtensor<size_t, 1> FineLayer::nodesTopEdge() const
{
size_t nely = m_nhy.size();
return m_startNode(nely) + xt::arange<size_t>(m_nelx(nely - 1) + 1);
}
inline xt::xtensor<size_t, 1> FineLayer::nodesLeftEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j));
xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1));
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::nodesRightEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j)) =
xt::view(m_startNode, xt::range(i, j)) + xt::view(m_nelx, xt::range(i, j));
xt::view(ret, xt::range(k + 1, l + 1)) =
xt::view(m_startNode, xt::range(k + 1, l + 1)) + xt::view(m_nelx, xt::range(k, l));
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::nodesBottomOpenEdge() const
{
return m_startNode(0) + xt::arange<size_t>(1, m_nelx(0));
}
inline xt::xtensor<size_t, 1> FineLayer::nodesTopOpenEdge() const
{
size_t nely = m_nhy.size();
return m_startNode(nely) + xt::arange<size_t>(1, m_nelx(nely - 1));
}
inline xt::xtensor<size_t, 1> FineLayer::nodesLeftOpenEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j));
xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l));
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::nodesRightOpenEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j - 1)) =
xt::view(m_startNode, xt::range(i + 1, j)) + xt::view(m_nelx, xt::range(i + 1, j));
xt::view(ret, xt::range(k, l - 1)) =
xt::view(m_startNode, xt::range(k + 1, l)) + xt::view(m_nelx, xt::range(k, l - 1));
return ret;
}
inline size_t FineLayer::nodesBottomLeftCorner() const
{
return m_startNode(0);
}
inline size_t FineLayer::nodesBottomRightCorner() const
{
return m_startNode(0) + m_nelx(0);
}
inline size_t FineLayer::nodesTopLeftCorner() const
{
size_t nely = m_nhy.size();
return m_startNode(nely);
}
inline size_t FineLayer::nodesTopRightCorner() const
{
size_t nely = m_nhy.size();
return m_startNode(nely) + m_nelx(nely - 1);
}
inline size_t FineLayer::nodesLeftBottomCorner() const
{
return nodesBottomLeftCorner();
}
inline size_t FineLayer::nodesRightBottomCorner() const
{
return nodesBottomRightCorner();
}
inline size_t FineLayer::nodesLeftTopCorner() const
{
return nodesTopLeftCorner();
}
inline size_t FineLayer::nodesRightTopCorner() const
{
return nodesTopRightCorner();
}
inline xt::xtensor<size_t, 2> FineLayer::nodesPeriodic() const
{
xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);
ret(0, 0) = nodesBottomLeftCorner();
ret(0, 1) = nodesBottomRightCorner();
ret(1, 0) = nodesBottomLeftCorner();
ret(1, 1) = nodesTopRightCorner();
ret(2, 0) = nodesBottomLeftCorner();
ret(2, 1) = nodesTopLeftCorner();
size_t i = 3;
xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
i += bot.size();
xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
return ret;
}
inline size_t FineLayer::nodesOrigin() const
{
return nodesBottomLeftCorner();
}
inline xt::xtensor<size_t, 2> FineLayer::dofs() const
{
return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
}
inline xt::xtensor<size_t, 2> FineLayer::dofsPeriodic() const
{
xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
for (size_t j = 0; j < m_ndim; ++j) {
xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
}
return GooseFEM::Mesh::renumber(ret);
}
-inline xt::xtensor<size_t, 1> FineLayer::roll(size_t n, size_t axis)
+inline xt::xtensor<size_t, 1> FineLayer::roll(size_t n)
{
- GOOSEFEM_ASSERT(axis == 0);
auto conn = this->conn();
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelem});
// loop over all element layers
for (size_t iy = 0; iy < nely; ++iy) {
// no refinement
size_t shift = n * (m_nelx(iy) / m_nelx(0));
size_t nel = m_nelx(iy);
// refinement
if (m_refine(iy) != -1) {
shift = n * (m_nelx(iy) / m_nelx(0)) * 4;
nel = m_nelx(iy) * 4;
}
// element numbers of the layer, and roll them
auto e = m_startElem(iy) + xt::arange<size_t>(nel);
xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift);
}
return ret;
}
inline void FineLayer::map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
{
GOOSEFEM_ASSERT(coor.shape(1) == 2);
GOOSEFEM_ASSERT(conn.shape(1) == 4);
GOOSEFEM_ASSERT(conn.shape(0) > 0);
GOOSEFEM_ASSERT(coor.shape(0) >= 4);
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
auto n0 = xt::view(conn, xt::all(), 0);
auto n1 = xt::view(conn, xt::all(), 1);
auto n2 = xt::view(conn, xt::all(), 2);
auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0);
auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1);
auto hx = xt::amin(xt::where(dx > 0.0, dx, std::numeric_limits<double>::max()))();
auto hy = xt::amin(xt::where(dy > 0.0, dx, std::numeric_limits<double>::max()))();
GOOSEFEM_CHECK(xt::allclose(hx, hy));
if (conn.shape(0) == 1) {
this->init(1, 1, hx);
GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
return;
}
size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2;
size_t eleft = emid;
size_t eright = emid;
while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) {
eleft--;
}
while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) {
eright++;
}
GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0));
size_t nelx = eright - eleft + 1;
size_t nely = static_cast<size_t>((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx);
this->init(nelx, nely, hx);
GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
GOOSEFEM_CHECK(xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange<size_t>(nelx))));
}
namespace Map {
inline RefineRegular::RefineRegular(
const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny)
: m_coarse(mesh)
{
m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h());
xt::xtensor<size_t, 2> elmat_coarse = m_coarse.elementMatrix();
xt::xtensor<size_t, 2> elmat_fine = m_fine.elementMatrix();
m_coarse2fine = xt::empty<size_t>({m_coarse.nelem(), nx * ny});
for (size_t i = 0; i < elmat_coarse.shape(0); ++i) {
for (size_t j = 0; j < elmat_coarse.shape(1); ++j) {
xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view(
elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx)));
}
}
}
inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getCoarseMesh() const
{
return m_coarse;
}
inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getFineMesh() const
{
return m_fine;
}
inline xt::xtensor<size_t, 2> RefineRegular::getMap() const
{
return m_coarse2fine;
}
inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 1>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
size_t m = m_coarse2fine.shape(0);
size_t n = m_coarse2fine.shape(1);
xt::xtensor<double, 2> ret = xt::empty<double>({m, n});
for (size_t i = 0; i < m; ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, i) = xt::view(data, xt::keep(e));
}
return ret;
}
inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 2>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
size_t m = m_coarse2fine.shape(0);
size_t n = m_coarse2fine.shape(1);
size_t N = data.shape(1);
xt::xtensor<double, 2> ret = xt::empty<double>({m, n * N});
for (size_t i = 0; i < m; ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
for (size_t q = 0; q < data.shape(1); ++q) {
xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
}
}
return ret;
}
inline xt::xtensor<double, 4> RefineRegular::mapToCoarse(const xt::xtensor<double, 4>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
size_t m = m_coarse2fine.shape(0);
size_t n = m_coarse2fine.shape(1);
size_t N = data.shape(1);
xt::xtensor<double, 4> ret = xt::empty<double>({m, n * N, data.shape(2), data.shape(3)});
for (size_t i = 0; i < m; ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
for (size_t q = 0; q < data.shape(1); ++q) {
xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
}
}
return ret;
}
inline xt::xtensor<double, 1> RefineRegular::mapToFine(const xt::xtensor<double, 1>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
xt::xtensor<double, 1> ret = xt::empty<double>({m_coarse2fine.size()});
for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, xt::keep(e)) = data(i);
}
return ret;
}
inline xt::xtensor<double, 2> RefineRegular::mapToFine(const xt::xtensor<double, 2>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
xt::xtensor<double, 2> ret = xt::empty<double>({m_coarse2fine.size(), data.shape(1)});
for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, xt::keep(e)) = xt::view(data, i);
}
return ret;
}
inline xt::xtensor<double, 4> RefineRegular::mapToFine(const xt::xtensor<double, 4>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
xt::xtensor<double, 4> ret =
xt::empty<double>({m_coarse2fine.size(), data.shape(1), data.shape(2), data.shape(3)});
for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, xt::keep(e)) = xt::view(data, i);
}
return ret;
}
inline FineLayer2Regular::FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh)
: m_finelayer(mesh)
{
// ------------
// Regular-mesh
// ------------
m_regular = GooseFEM::Mesh::Quad4::Regular(
xt::amax(m_finelayer.m_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h);
// -------
// mapping
// -------
// allocate mapping
m_elem_regular.resize(m_finelayer.m_nelem);
m_frac_regular.resize(m_finelayer.m_nelem);
// alias
xt::xtensor<size_t, 1> nhx = m_finelayer.m_nhx;
xt::xtensor<size_t, 1> nhy = m_finelayer.m_nhy;
xt::xtensor<size_t, 1> nelx = m_finelayer.m_nelx;
xt::xtensor<size_t, 1> start = m_finelayer.m_startElem;
// 'matrix' of element numbers of the Regular-mesh
xt::xtensor<size_t, 2> elementMatrix = m_regular.elementMatrix();
// cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh
xt::xtensor<size_t, 1> cum_nhy =
xt::concatenate(xt::xtuple(xt::zeros<size_t>({1}), xt::cumsum(nhy)));
// number of element layers in y-direction of the FineLayer-mesh
size_t nely = nhy.size();
// loop over layers of the FineLayer-mesh
for (size_t iy = 0; iy < nely; ++iy) {
// element numbers of the Regular-mesh along this layer of the FineLayer-mesh
auto el_new = xt::view(elementMatrix, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all());
// no coarsening/refinement
// ------------------------
if (m_finelayer.m_refine(iy) == -1) {
// element numbers of the FineLayer-mesh along this layer
xt::xtensor<size_t, 1> el_old = start(iy) + xt::arange<size_t>(nelx(iy));
// loop along this layer of the FineLayer-mesh
for (size_t ix = 0; ix < nelx(iy); ++ix) {
// get the element numbers of the Regular-mesh for this element of the
// FineLayer-mesh
auto block =
xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
// write to mapping
for (auto& i : block) {
m_elem_regular[el_old(ix)].push_back(i);
m_frac_regular[el_old(ix)].push_back(1.0);
}
}
}
// refinement along the x-direction (below the middle layer)
// ---------------------------------------------------------
else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
// element numbers of the FineLayer-mesh along this layer
// rows: coarse block, columns element numbers per block
xt::xtensor<size_t, 2> el_old =
start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
// loop along this layer of the FineLayer-mesh
for (size_t ix = 0; ix < nelx(iy); ++ix) {
// get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
auto block =
xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
// bottom: wide-to-narrow
{
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(block, j, xt::range(j, nhx(iy) - j));
m_elem_regular[el_old(ix, 0)].push_back(e(0));
m_frac_regular[el_old(ix, 0)].push_back(0.5);
for (size_t k = 1; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 0)].push_back(e(k));
m_frac_regular[el_old(ix, 0)].push_back(1.0);
}
m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 0)].push_back(0.5);
}
}
// top: regular small element
{
auto e = xt::view(
block,
xt::range(nhy(iy) / 2, nhy(iy)),
xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 2)].push_back(i);
m_frac_regular[el_old(ix, 2)].push_back(1.0);
}
}
// left
{
// left-bottom: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(block, j, xt::range(0, j + 1));
for (size_t k = 0; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 3)].push_back(e(k));
m_frac_regular[el_old(ix, 3)].push_back(1.0);
}
m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 3)].push_back(0.5);
}
// left-top: regular
{
auto e = xt::view(
block,
xt::range(nhy(iy) / 2, nhy(iy)),
xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 3)].push_back(i);
m_frac_regular[el_old(ix, 3)].push_back(1.0);
}
}
}
// right
{
// right-bottom: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy)));
m_elem_regular[el_old(ix, 1)].push_back(e(0));
m_frac_regular[el_old(ix, 1)].push_back(0.5);
for (size_t k = 1; k < e.size(); ++k) {
m_elem_regular[el_old(ix, 1)].push_back(e(k));
m_frac_regular[el_old(ix, 1)].push_back(1.0);
}
}
// right-top: regular
{
auto e = xt::view(
block,
xt::range(nhy(iy) / 2, nhy(iy)),
xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 1)].push_back(i);
m_frac_regular[el_old(ix, 1)].push_back(1.0);
}
}
}
}
}
// coarsening along the x-direction (above the middle layer)
else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) {
// element numbers of the FineLayer-mesh along this layer
// rows: coarse block, columns element numbers per block
xt::xtensor<size_t, 2> el_old =
start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
// loop along this layer of the FineLayer-mesh
for (size_t ix = 0; ix < nelx(iy); ++ix) {
// get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
auto block =
xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
// top: narrow-to-wide
{
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(
block,
nhy(iy) / 2 + j,
xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1));
m_elem_regular[el_old(ix, 3)].push_back(e(0));
m_frac_regular[el_old(ix, 3)].push_back(0.5);
for (size_t k = 1; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 3)].push_back(e(k));
m_frac_regular[el_old(ix, 3)].push_back(1.0);
}
m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 3)].push_back(0.5);
}
}
// bottom: regular small element
{
auto e = xt::view(
block,
xt::range(0, nhy(iy) / 2),
xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 1)].push_back(i);
m_frac_regular[el_old(ix, 1)].push_back(1.0);
}
}
// left
{
// left-bottom: regular
{
auto e = xt::view(
block,
xt::range(0, nhy(iy) / 2),
xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 0)].push_back(i);
m_frac_regular[el_old(ix, 0)].push_back(1.0);
}
}
// left-top: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e =
xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j));
for (size_t k = 0; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 0)].push_back(e(k));
m_frac_regular[el_old(ix, 0)].push_back(1.0);
}
m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 0)].push_back(0.5);
}
}
// right
{
// right-bottom: regular
{
auto e = xt::view(
block,
xt::range(0, nhy(iy) / 2),
xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 2)].push_back(i);
m_frac_regular[el_old(ix, 2)].push_back(1.0);
}
}
// right-top: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(
block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy)));
m_elem_regular[el_old(ix, 2)].push_back(e(0));
m_frac_regular[el_old(ix, 2)].push_back(0.5);
for (size_t k = 1; k < e.size(); ++k) {
m_elem_regular[el_old(ix, 2)].push_back(e(k));
m_frac_regular[el_old(ix, 2)].push_back(1.0);
}
}
}
}
}
}
}
inline GooseFEM::Mesh::Quad4::Regular FineLayer2Regular::getRegularMesh() const
{
return m_regular;
}
inline GooseFEM::Mesh::Quad4::FineLayer FineLayer2Regular::getFineLayerMesh() const
{
return m_finelayer;
}
inline std::vector<std::vector<size_t>> FineLayer2Regular::getMap() const
{
return m_elem_regular;
}
inline std::vector<std::vector<double>> FineLayer2Regular::getMapFraction() const
{
return m_frac_regular;
}
inline xt::xtensor<double, 1>
FineLayer2Regular::mapToRegular(const xt::xtensor<double, 1>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
xt::xtensor<double, 1> ret = xt::zeros<double>({m_regular.nelem()});
for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
ret(m_elem_regular[e][i]) += m_frac_regular[e][i] * data(e);
}
}
return ret;
}
inline xt::xtensor<double, 2>
FineLayer2Regular::mapToRegular(const xt::xtensor<double, 2>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
xt::xtensor<double, 2> ret = xt::zeros<double>({m_regular.nelem(), data.shape(1)});
for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
}
}
return ret;
}
inline xt::xtensor<double, 4>
FineLayer2Regular::mapToRegular(const xt::xtensor<double, 4>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
xt::xtensor<double, 4> ret =
xt::zeros<double>({m_regular.nelem(), data.shape(1), data.shape(2), data.shape(3)});
for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
}
}
return ret;
}
} // namespace Map
} // namespace Quad4
} // namespace Mesh
} // namespace GooseFEM
#endif
diff --git a/python/Mesh.hpp b/python/Mesh.hpp
index 5fed2d7..93f69d0 100644
--- a/python/Mesh.hpp
+++ b/python/Mesh.hpp
@@ -1,114 +1,157 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#include <GooseFEM/GooseFEM.h>
#include <pybind11/pybind11.h>
#include <pyxtensor/pyxtensor.hpp>
namespace py = pybind11;
void init_Mesh(py::module& m)
{
py::enum_<GooseFEM::Mesh::ElementType>(m, "ElementType", "ElementType")
.value("Tri3", GooseFEM::Mesh::ElementType::Tri3)
.value("Quad4", GooseFEM::Mesh::ElementType::Quad4)
.value("Hex8", GooseFEM::Mesh::ElementType::Hex8)
.export_values();
py::class_<GooseFEM::Mesh::Renumber>(m, "Renumber")
.def(
py::init<const xt::xarray<size_t>&>(),
"Class to renumber to the lowest possible indices. Use ``Renumber(...).get()`` to get "
"the renumbered result",
py::arg("dofs"))
.def("get", &GooseFEM::Mesh::Renumber::get, "Get result of renumbering")
.def(
"index",
&GooseFEM::Mesh::Renumber::index,
"Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
.def(
"__repr__", [](const GooseFEM::Mesh::Renumber&) { return "<GooseFEM.Mesh.Renumber>"; });
py::class_<GooseFEM::Mesh::Reorder>(m, "Reorder")
.def(py::init([](xt::xtensor<size_t, 1>& a) { return new GooseFEM::Mesh::Reorder({a}); }))
.def(py::init([](xt::xtensor<size_t, 1>& a, xt::xtensor<size_t, 1>& b) {
return new GooseFEM::Mesh::Reorder({a, b});
}))
.def(py::init(
[](xt::xtensor<size_t, 1>& a, xt::xtensor<size_t, 1>& b, xt::xtensor<size_t, 1>& c) {
return new GooseFEM::Mesh::Reorder({a, b, c});
}))
.def(py::init([](xt::xtensor<size_t, 1>& a,
xt::xtensor<size_t, 1>& b,
xt::xtensor<size_t, 1>& c,
xt::xtensor<size_t, 1>& d) {
return new GooseFEM::Mesh::Reorder({a, b, c, d});
}))
.def("get", &GooseFEM::Mesh::Reorder::get, "Reorder matrix (e.g. ``dofs``)")
.def(
"index",
&GooseFEM::Mesh::Reorder::index,
"Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
.def("__repr__", [](const GooseFEM::Mesh::Reorder&) { return "<GooseFEM.Mesh.Reorder>"; });
m.def(
"dofs",
&GooseFEM::Mesh::dofs,
"List with DOF-numbers (in sequential order)",
py::arg("nnode"),
py::arg("ndim"));
m.def(
"renumber",
&GooseFEM::Mesh::renumber,
"Renumber to lowest possible indices. Use ``GooseFEM.Mesh.Renumber`` for advanced "
"functionality",
py::arg("dofs"));
m.def(
"coordination",
&GooseFEM::Mesh::coordination,
"Coordination number of each node (number of elements connected to each node)",
py::arg("conn"));
m.def(
"elem2node",
&GooseFEM::Mesh::elem2node,
"Element-numbers connected to each node",
py::arg("conn"),
py::arg("sorted") = true);
m.def(
"edgesize",
py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(
&GooseFEM::Mesh::edgesize),
"Get the edge size of all elements",
py::arg("coor"),
py::arg("conn"));
m.def(
"edgesize",
py::overload_cast<
const xt::xtensor<double, 2>&,
const xt::xtensor<size_t, 2>&,
GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize),
"Get the edge size of all elements",
py::arg("coor"),
py::arg("conn"),
py::arg("type"));
+
+ m.def(
+ "centers",
+ py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(
+ &GooseFEM::Mesh::centers),
+ "Coordinates of the center of each element",
+ py::arg("coor"),
+ py::arg("conn"));
+
+ m.def(
+ "centers",
+ py::overload_cast<
+ const xt::xtensor<double, 2>&,
+ const xt::xtensor<size_t, 2>&,
+ GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::centers),
+ "Coordinates of the center of each element",
+ py::arg("coor"),
+ py::arg("conn"),
+ py::arg("type"));
+
+ m.def(
+ "elemmap2nodemap",
+ py::overload_cast<
+ const xt::xtensor<size_t, 1>&,
+ const xt::xtensor<double, 2>&,
+ const xt::xtensor<size_t, 2>&>(&GooseFEM::Mesh::elemmap2nodemap),
+ "Convert an element-map to a node-map",
+ py::arg("elem_map"),
+ py::arg("coor"),
+ py::arg("conn"));
+
+ m.def(
+ "elemmap2nodemap",
+ py::overload_cast<
+ const xt::xtensor<size_t, 1>&,
+ const xt::xtensor<double, 2>&,
+ const xt::xtensor<size_t, 2>&,
+ GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::elemmap2nodemap),
+ "Convert an element-map to a node-map",
+ py::arg("elem_map"),
+ py::arg("coor"),
+ py::arg("conn"),
+ py::arg("type"));
}
diff --git a/python/MeshQuad4.hpp b/python/MeshQuad4.hpp
index 8cdb385..26a0985 100644
--- a/python/MeshQuad4.hpp
+++ b/python/MeshQuad4.hpp
@@ -1,261 +1,263 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#include <GooseFEM/GooseFEM.h>
#include <pybind11/pybind11.h>
#include <pyxtensor/pyxtensor.hpp>
namespace py = pybind11;
void init_MeshQuad4(py::module& m)
{
py::class_<GooseFEM::Mesh::Quad4::Regular>(m, "Regular")
.def(
py::init<size_t, size_t, double>(),
"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("nelx", &GooseFEM::Mesh::Quad4::Regular::nelx)
.def("nely", &GooseFEM::Mesh::Quad4::Regular::nely)
.def("h", &GooseFEM::Mesh::Quad4::Regular::h)
.def("getElementType", &GooseFEM::Mesh::Quad4::Regular::getElementType)
.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("dofs", &GooseFEM::Mesh::Quad4::Regular::dofs)
.def("nodesPeriodic", &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic)
.def("nodesOrigin", &GooseFEM::Mesh::Quad4::Regular::nodesOrigin)
.def("dofsPeriodic", &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic)
.def("elementMatrix", &GooseFEM::Mesh::Quad4::Regular::elementMatrix)
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) {
return "<GooseFEM.Mesh.Quad4.Regular>";
});
py::class_<GooseFEM::Mesh::Quad4::FineLayer>(m, "FineLayer")
.def(
py::init<size_t, size_t, double, size_t>(),
"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)
.def(
py::init<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(),
"Map connectivity to generating FineLayer-object.",
py::arg("coor"),
py::arg("conn"))
.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("nelx", &GooseFEM::Mesh::Quad4::FineLayer::nelx)
.def("nely", &GooseFEM::Mesh::Quad4::FineLayer::nely)
.def("h", &GooseFEM::Mesh::Quad4::FineLayer::h)
.def("getElementType", &GooseFEM::Mesh::Quad4::FineLayer::getElementType)
.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("dofs", &GooseFEM::Mesh::Quad4::FineLayer::dofs)
.def("nodesPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic)
.def("nodesOrigin", &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin)
.def("dofsPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic)
+ .def("roll", &GooseFEM::Mesh::Quad4::FineLayer::roll)
+
.def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer&) {
return "<GooseFEM.Mesh.Quad4.FineLayer>";
});
}
void init_MeshQuad4Map(py::module& m)
{
py::class_<GooseFEM::Mesh::Quad4::Map::RefineRegular>(m, "RefineRegular")
.def(
py::init<const GooseFEM::Mesh::Quad4::Regular&, size_t, size_t>(),
"Refine a regular mesh",
py::arg("mesh"),
py::arg("nx"),
py::arg("ny"))
.def("getCoarseMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getCoarseMesh)
.def("getFineMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getFineMesh)
.def("getMap", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getMap)
.def(
"mapToCoarse",
py::overload_cast<const xt::xtensor<double, 1>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
.def(
"mapToCoarse",
py::overload_cast<const xt::xtensor<double, 2>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
.def(
"mapToCoarse",
py::overload_cast<const xt::xtensor<double, 4>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
.def(
"mapToFine",
py::overload_cast<const xt::xtensor<double, 1>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
.def(
"mapToFine",
py::overload_cast<const xt::xtensor<double, 2>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
.def(
"mapToFine",
py::overload_cast<const xt::xtensor<double, 4>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::RefineRegular&) {
return "<GooseFEM.Mesh.Quad4.Map.RefineRegular>";
});
py::class_<GooseFEM::Mesh::Quad4::Map::FineLayer2Regular>(m, "FineLayer2Regular")
.def(
py::init<const GooseFEM::Mesh::Quad4::FineLayer&>(),
"Map a FineLayer mesh to a Regular mesh",
py::arg("mesh"))
.def("getRegularMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getRegularMesh)
.def("getFineLayerMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getFineLayerMesh)
.def("getMap", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMap)
.def("getMapFraction", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMapFraction)
.def(
"mapToRegular",
py::overload_cast<const xt::xtensor<double, 1>&>(
&GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
.def(
"mapToRegular",
py::overload_cast<const xt::xtensor<double, 2>&>(
&GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
.def(
"mapToRegular",
py::overload_cast<const xt::xtensor<double, 4>&>(
&GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular&) {
return "<GooseFEM.Mesh.Quad4.Map.FineLayer2Regular>";
});
}
diff --git a/test/basic/Mesh.cpp b/test/basic/Mesh.cpp
index b6f126b..b51fe14 100644
--- a/test/basic/Mesh.cpp
+++ b/test/basic/Mesh.cpp
@@ -1,62 +1,278 @@
#include <catch2/catch.hpp>
#include <xtensor/xrandom.hpp>
#include <xtensor/xmath.hpp>
#include <GooseFEM/GooseFEM.h>
#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12));
TEST_CASE("GooseFEM::Mesh", "Mesh.h")
{
SECTION("edgesize")
{
GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 10.0);
auto s = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn());
auto t = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn(), mesh.getElementType());
REQUIRE(xt::allclose(s, 10.0));
REQUIRE(xt::allclose(t, 10.0));
}
SECTION("coordination")
{
GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
auto N = GooseFEM::Mesh::coordination(mesh.conn());
xt::xtensor<size_t, 1> ret = {1, 2, 2, 1, 2, 4, 4, 2, 2, 4, 4, 2, 1, 2, 2, 1};
REQUIRE(xt::all(xt::equal(N, ret)));
}
SECTION("centers")
{
GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 2.0);
xt::xtensor<double, 2> c = {
{1.0, 1.0},
{3.0, 1.0},
{1.0, 3.0},
{3.0, 3.0}};
REQUIRE(xt::allclose(GooseFEM::Mesh::centers(mesh.coor(), mesh.conn()), c));
}
SECTION("elem2node")
{
GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
auto tonode = GooseFEM::Mesh::elem2node(mesh.conn());
REQUIRE(tonode.size() == 16);
REQUIRE(tonode[0] == std::vector<size_t>{0});
REQUIRE(tonode[1] == std::vector<size_t>{0, 1});
REQUIRE(tonode[2] == std::vector<size_t>{1, 2});
REQUIRE(tonode[3] == std::vector<size_t>{2});
REQUIRE(tonode[4] == std::vector<size_t>{0, 3});
REQUIRE(tonode[5] == std::vector<size_t>{0, 1, 3, 4});
REQUIRE(tonode[6] == std::vector<size_t>{1, 2, 4, 5});
REQUIRE(tonode[7] == std::vector<size_t>{2, 5});
REQUIRE(tonode[8] == std::vector<size_t>{3, 6});
REQUIRE(tonode[9] == std::vector<size_t>{3, 4, 6, 7});
REQUIRE(tonode[10] == std::vector<size_t>{4, 5, 7, 8});
REQUIRE(tonode[11] == std::vector<size_t>{5, 8});
REQUIRE(tonode[12] == std::vector<size_t>{6});
REQUIRE(tonode[13] == std::vector<size_t>{6, 7});
REQUIRE(tonode[14] == std::vector<size_t>{7, 8});
REQUIRE(tonode[15] == std::vector<size_t>{8});
}
+
+ SECTION("elemmap2nodemap")
+ {
+ GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
+
+ xt::xtensor<size_t, 1> elmap0 = {
+ 0, 1, 2,
+ 3, 4, 5,
+ 6, 7, 8
+ };
+ xt::xtensor<size_t, 1> elmap1 = {
+ 2, 0, 1,
+ 5, 3, 4,
+ 8, 6, 7
+ };
+ xt::xtensor<size_t, 1> elmap2 = {
+ 1, 2, 0,
+ 4, 5, 3,
+ 7, 8, 6
+ };
+
+ xt::xtensor<size_t, 1> nodemap0 = {
+ 0, 1, 2, 3,
+ 4, 5, 6, 7,
+ 8, 9, 10, 11,
+ 12, 13, 14, 15
+ };
+ xt::xtensor<size_t, 1> nodemap1 = {
+ 2, 0, 1, 2,
+ 6, 4, 5, 6,
+ 10, 8, 9, 10,
+ 14, 15, 13, 14
+ };
+ xt::xtensor<size_t, 1> nodemap2 = {
+ 1, 2, 0, 1,
+ 5, 6, 4, 5,
+ 9, 10, 8, 9,
+ 13, 14, 15, 13
+ };
+
+ REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap0, mesh.coor(), mesh.conn()), nodemap0)));
+ REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap1, mesh.coor(), mesh.conn()), nodemap1)));
+ REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap2, mesh.coor(), mesh.conn()), nodemap2)));
+ }
+
+ SECTION("elemmap2nodemap - example 1")
+ {
+ GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3);
+
+ xt::xtensor<int, 1> elemval = {
+ 1, 0, 0,
+ 1, 0, 0,
+ 1, 0, 0
+ };
+
+ xt::xtensor<int, 1> elemval_r1 = {
+ 0, 1, 0,
+ 0, 1, 0,
+ 0, 1, 0
+ };
+
+ xt::xtensor<int, 1> elemval_r2 = {
+ 0, 0, 1,
+ 0, 0, 1,
+ 0, 0, 1
+ };
+
+ xt::xtensor<int, 1> nodeval = {
+ 1, 0, 0, 1,
+ 1, 0, 0, 1,
+ 1, 0, 0, 1,
+ 1, 0, 0, 1
+ };
+
+ xt::xtensor<int, 1> nodeval_r1 = {
+ 0, 1, 0, 0,
+ 0, 1, 0, 0,
+ 0, 1, 0, 0,
+ 0, 1, 0, 0
+ };
+
+ xt::xtensor<int, 1> nodeval_r2 = {
+ 0, 0, 1, 0,
+ 0, 0, 1, 0,
+ 0, 0, 1, 0,
+ 0, 0, 1, 0
+ };
+
+ {
+ auto elemmap = mesh.roll(0);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(1);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(2);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(3);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(4);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(5);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+ }
+
+ SECTION("elemmap2nodemap - example 2")
+ {
+ GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3);
+
+ xt::xtensor<int, 1> elemval = {
+ 1, 0, 0,
+ 0, 1, 0,
+ 0, 0, 1
+ };
+
+ xt::xtensor<int, 1> elemval_r1 = {
+ 0, 1, 0,
+ 0, 0, 1,
+ 1, 0, 0
+ };
+
+ xt::xtensor<int, 1> elemval_r2 = {
+ 0, 0, 1,
+ 1, 0, 0,
+ 0, 1, 0
+ };
+
+ xt::xtensor<int, 1> nodeval = {
+ 1, 0, 0, 1,
+ 0, 1, 0, 0,
+ 0, 0, 1, 0,
+ 1, 0, 0, 1
+ };
+
+ xt::xtensor<int, 1> nodeval_r1 = {
+ 0, 1, 0, 0,
+ 0, 0, 1, 0,
+ 1, 0, 0, 1,
+ 0, 1, 0, 0
+ };
+
+ xt::xtensor<int, 1> nodeval_r2 = {
+ 0, 0, 1, 0,
+ 1, 0, 0, 1,
+ 0, 1, 0, 0,
+ 0, 0, 1, 0
+ };
+
+ {
+ auto elemmap = mesh.roll(0);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(1);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(2);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(3);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(4);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+
+ {
+ auto elemmap = mesh.roll(5);
+ auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
+ REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
+ }
+ }
}