diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 42721ae..19d0aed 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,72 +1,76 @@
name: CI
on:
workflow_dispatch:
pull_request:
push:
branches:
- master
jobs:
default-shell:
strategy:
fail-fast: false
matrix:
runs-on: [ubuntu-latest, macos-latest, windows-latest]
defaults:
run:
shell: bash -l {0}
name: "${{ matrix.runs-on }} • x64 ${{ matrix.args }}"
runs-on: ${{ matrix.runs-on }}
steps:
- name: Basic GitHub action setup
uses: actions/checkout@v2
- name: Set conda environment "test"
uses: conda-incubator/setup-miniconda@v2
with:
mamba-version: "*"
channels: conda-forge,defaults
channel-priority: true
environment-file: environment.yaml
activate-environment: test
auto-activate-base: false
- name: Extending conda environment for testing
run: conda env update --file environment_test.yaml
- name: Extending conda environment to run examples
if: runner.os == 'Linux'
run: conda env update --file environment_examples.yaml
- name: Configure using CMake (with examples)
if: runner.os == 'Linux'
run: cmake . -DBUILD_TESTS=ON -DBUILD_EXAMPLES=ON
- name: Configure using CMake (without examples)
if: runner.os != 'Linux'
run: cmake . -DBUILD_TESTS=ON -DBUILD_EXAMPLES=OFF
- name: Build C++ tests & examples
run: cmake --build .
- name: Run C++ tests & examples
run: cmake --build . --target "RUN_TESTS_AND_EXAMPLES"
- name: Build and install Python module
run: |
python setup.py build
python setup.py install
+ - name: Python tests
+ run: |
+ python test/basic-python/MeshQuad4.py
+
- name: Run Python examples
if: runner.os == 'Linux'
run: |
python ./docs/examples/statics/FixedDisplacements_LinearElastic/example.py --no-plot
python ./docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py --no-plot
diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h
index 1aafc67..000e9d7 100644
--- a/include/GooseFEM/MeshQuad4.h
+++ b/include/GooseFEM/MeshQuad4.h
@@ -1,265 +1,254 @@
/*
(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 Quad4
-} // namespace Mesh
-} // namespace GooseFEM
+ class FineLayer2Regular;
+} // namespace Map
-namespace GooseFEM {
-namespace Mesh {
-namespace Quad4 {
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
};
-} // namespace Quad4
-} // namespace Mesh
-} // namespace GooseFEM
-namespace GooseFEM {
-namespace Mesh {
-namespace Quad4 {
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;
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 Quad4
-} // namespace Mesh
-} // namespace GooseFEM
-namespace GooseFEM {
-namespace Mesh {
-namespace Quad4 {
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;
};
-} // namespace Map
-} // namespace Quad4
-} // namespace Mesh
-} // namespace GooseFEM
-namespace GooseFEM {
-namespace Mesh {
-namespace Quad4 {
-namespace Map {
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 7f1758c..0fc8465 100644
--- a/include/GooseFEM/MeshQuad4.hpp
+++ b/include/GooseFEM/MeshQuad4.hpp
@@ -1,1325 +1,1383 @@
/*
(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) : m_h(h)
+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 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/MeshQuad4.hpp b/python/MeshQuad4.hpp
index 2a038c1..8cdb385 100644
--- a/python/MeshQuad4.hpp
+++ b/python/MeshQuad4.hpp
@@ -1,255 +1,261 @@
/* =================================================================================================
(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("__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-python/MeshQuad4.py b/test/basic-python/MeshQuad4.py
new file mode 100644
index 0000000..f650392
--- /dev/null
+++ b/test/basic-python/MeshQuad4.py
@@ -0,0 +1,16 @@
+import unittest
+import numpy as np
+import GooseFEM as gf
+
+class Test_MeshQuad4(unittest.TestCase):
+
+ def test_FineLayer_replica(self):
+
+ mesh = gf.Mesh.Quad4.FineLayer(27, 27)
+ replica = gf.Mesh.Quad4.FineLayer(mesh.coor(), mesh.conn())
+ self.assertTrue(np.allclose(mesh.coor(), replica.coor()))
+ self.assertTrue(np.all(np.equal(mesh.conn(), replica.conn())))
+
+if __name__ == '__main__':
+
+ unittest.main()
diff --git a/test/basic/MeshQuad4.cpp b/test/basic/MeshQuad4.cpp
index 7de0cf1..e1e89ea 100644
--- a/test/basic/MeshQuad4.cpp
+++ b/test/basic/MeshQuad4.cpp
@@ -1,372 +1,396 @@
#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::MeshQuad4", "MeshQuad4.h")
{
SECTION("Regular")
{
xt::xtensor<double, 2> coor_ = {{0., 0.}, {1., 0.}, {2., 0.}, {3., 0.}, {4., 0.}, {5., 0.},
{0., 1.}, {1., 1.}, {2., 1.}, {3., 1.}, {4., 1.}, {5., 1.},
{0., 2.}, {1., 2.}, {2., 2.}, {3., 2.}, {4., 2.}, {5., 2.},
{0., 3.}, {1., 3.}, {2., 3.}, {3., 3.}, {4., 3.}, {5., 3.}};
xt::xtensor<double, 2> conn_ = {{0, 1, 7, 6},
{1, 2, 8, 7},
{2, 3, 9, 8},
{3, 4, 10, 9},
{4, 5, 11, 10},
{6, 7, 13, 12},
{7, 8, 14, 13},
{8, 9, 15, 14},
{9, 10, 16, 15},
{10, 11, 17, 16},
{12, 13, 19, 18},
{13, 14, 20, 19},
{14, 15, 21, 20},
{15, 16, 22, 21},
{16, 17, 23, 22}};
size_t nelem_ = 15;
size_t nnode_ = 24;
size_t nne_ = 4;
size_t ndim_ = 2;
size_t nnodePeriodic_ = 15;
xt::xtensor<size_t, 1> nodesBottomEdge_ = {0, 1, 2, 3, 4, 5};
xt::xtensor<size_t, 1> nodesTopEdge_ = {18, 19, 20, 21, 22, 23};
xt::xtensor<size_t, 1> nodesLeftEdge_ = {0, 6, 12, 18};
xt::xtensor<size_t, 1> nodesRightEdge_ = {5, 11, 17, 23};
xt::xtensor<size_t, 1> nodesBottomOpenEdge_ = {1, 2, 3, 4};
xt::xtensor<size_t, 1> nodesTopOpenEdge_ = {19, 20, 21, 22};
xt::xtensor<size_t, 1> nodesLeftOpenEdge_ = {6, 12};
xt::xtensor<size_t, 1> nodesRightOpenEdge_ = {11, 17};
size_t nodesBottomLeftCorner_ = 0;
size_t nodesBottomRightCorner_ = 5;
size_t nodesTopLeftCorner_ = 18;
size_t nodesTopRightCorner_ = 23;
size_t nodesLeftBottomCorner_ = 0;
size_t nodesLeftTopCorner_ = 18;
size_t nodesRightBottomCorner_ = 5;
size_t nodesRightTopCorner_ = 23;
xt::xtensor<size_t, 2> dofs_ = {{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11},
{12, 13}, {14, 15}, {16, 17}, {18, 19}, {20, 21}, {22, 23},
{24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35},
{36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}};
xt::xtensor<size_t, 2> nodesPeriodic_ = {
{0, 5}, {0, 23}, {0, 18}, {1, 19}, {2, 20}, {3, 21}, {4, 22}, {6, 11}, {12, 17}};
size_t nodesOrigin_ = 0;
xt::xtensor<size_t, 2> dofsPeriodic_ = {
{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}, {10, 11}, {12, 13},
{14, 15}, {16, 17}, {18, 19}, {10, 11}, {20, 21}, {22, 23}, {24, 25}, {26, 27},
{28, 29}, {20, 21}, {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}};
GooseFEM::Mesh::Quad4::Regular mesh(5, 3);
xt::xtensor<double, 2> coor = mesh.coor();
xt::xtensor<size_t, 2> conn = mesh.conn();
size_t nelem = mesh.nelem();
size_t nnode = mesh.nnode();
size_t nne = mesh.nne();
size_t ndim = mesh.ndim();
size_t nnodePeriodic = mesh.nnode() - mesh.nodesPeriodic().shape(0);
xt::xtensor<size_t, 1> nodesBottomEdge = mesh.nodesBottomEdge();
xt::xtensor<size_t, 1> nodesTopEdge = mesh.nodesTopEdge();
xt::xtensor<size_t, 1> nodesLeftEdge = mesh.nodesLeftEdge();
xt::xtensor<size_t, 1> nodesRightEdge = mesh.nodesRightEdge();
xt::xtensor<size_t, 1> nodesBottomOpenEdge = mesh.nodesBottomOpenEdge();
xt::xtensor<size_t, 1> nodesTopOpenEdge = mesh.nodesTopOpenEdge();
xt::xtensor<size_t, 1> nodesLeftOpenEdge = mesh.nodesLeftOpenEdge();
xt::xtensor<size_t, 1> nodesRightOpenEdge = mesh.nodesRightOpenEdge();
size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner();
size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner();
size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner();
size_t nodesTopRightCorner = mesh.nodesTopRightCorner();
size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner();
size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner();
size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner();
size_t nodesRightTopCorner = mesh.nodesRightTopCorner();
xt::xtensor<size_t, 2> dofs = mesh.dofs();
xt::xtensor<size_t, 2> nodesPeriodic = mesh.nodesPeriodic();
size_t nodesOrigin = mesh.nodesOrigin();
xt::xtensor<size_t, 2> dofsPeriodic = mesh.dofsPeriodic();
REQUIRE(nelem == nelem_);
REQUIRE(nnode == nnode_);
REQUIRE(nne == nne_);
REQUIRE(ndim == ndim_);
REQUIRE(nnodePeriodic == nnodePeriodic_);
REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_);
REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_);
REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_);
REQUIRE(nodesTopRightCorner == nodesTopRightCorner_);
REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_);
REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_);
REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_);
REQUIRE(nodesRightTopCorner == nodesRightTopCorner_);
REQUIRE(nodesOrigin == nodesOrigin_);
REQUIRE(xt::allclose(coor, coor_));
REQUIRE(xt::all(xt::equal(conn, conn_)));
REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_)));
REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_)));
REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_)));
}
SECTION("FineLayer")
{
xt::xtensor<double, 2> coor_ = {
{0., 0.}, {3., 0.}, {6., 0.}, {9., 0.}, {0., 3.}, {3., 3.}, {6., 3.}, {9., 3.},
{0., 6.}, {3., 6.}, {6., 6.}, {9., 6.}, {1., 7.}, {2., 7.}, {4., 7.}, {5., 7.},
{7., 7.}, {8., 7.}, {0., 8.}, {1., 8.}, {2., 8.}, {3., 8.}, {4., 8.}, {5., 8.},
{6., 8.}, {7., 8.}, {8., 8.}, {9., 8.}, {0., 9.}, {1., 9.}, {2., 9.}, {3., 9.},
{4., 9.}, {5., 9.}, {6., 9.}, {7., 9.}, {8., 9.}, {9., 9.}, {0., 10.}, {1., 10.},
{2., 10.}, {3., 10.}, {4., 10.}, {5., 10.}, {6., 10.}, {7., 10.}, {8., 10.}, {9., 10.},
{0., 11.}, {1., 11.}, {2., 11.}, {3., 11.}, {4., 11.}, {5., 11.}, {6., 11.}, {7., 11.},
{8., 11.}, {9., 11.}, {0., 12.}, {1., 12.}, {2., 12.}, {3., 12.}, {4., 12.}, {5., 12.},
{6., 12.}, {7., 12.}, {8., 12.}, {9., 12.}, {0., 13.}, {1., 13.}, {2., 13.}, {3., 13.},
{4., 13.}, {5., 13.}, {6., 13.}, {7., 13.}, {8., 13.}, {9., 13.}, {1., 14.}, {2., 14.},
{4., 14.}, {5., 14.}, {7., 14.}, {8., 14.}, {0., 15.}, {3., 15.}, {6., 15.}, {9., 15.},
{0., 18.}, {3., 18.}, {6., 18.}, {9., 18.}, {0., 21.}, {3., 21.}, {6., 21.}, {9., 21.}};
xt::xtensor<double, 2> conn_ = {
{0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 9, 8},
{5, 6, 10, 9}, {6, 7, 11, 10}, {8, 9, 13, 12}, {9, 21, 20, 13},
{12, 13, 20, 19}, {8, 12, 19, 18}, {9, 10, 15, 14}, {10, 24, 23, 15},
{14, 15, 23, 22}, {9, 14, 22, 21}, {10, 11, 17, 16}, {11, 27, 26, 17},
{16, 17, 26, 25}, {10, 16, 25, 24}, {18, 19, 29, 28}, {19, 20, 30, 29},
{20, 21, 31, 30}, {21, 22, 32, 31}, {22, 23, 33, 32}, {23, 24, 34, 33},
{24, 25, 35, 34}, {25, 26, 36, 35}, {26, 27, 37, 36}, {28, 29, 39, 38},
{29, 30, 40, 39}, {30, 31, 41, 40}, {31, 32, 42, 41}, {32, 33, 43, 42},
{33, 34, 44, 43}, {34, 35, 45, 44}, {35, 36, 46, 45}, {36, 37, 47, 46},
{38, 39, 49, 48}, {39, 40, 50, 49}, {40, 41, 51, 50}, {41, 42, 52, 51},
{42, 43, 53, 52}, {43, 44, 54, 53}, {44, 45, 55, 54}, {45, 46, 56, 55},
{46, 47, 57, 56}, {48, 49, 59, 58}, {49, 50, 60, 59}, {50, 51, 61, 60},
{51, 52, 62, 61}, {52, 53, 63, 62}, {53, 54, 64, 63}, {54, 55, 65, 64},
{55, 56, 66, 65}, {56, 57, 67, 66}, {58, 59, 69, 68}, {59, 60, 70, 69},
{60, 61, 71, 70}, {61, 62, 72, 71}, {62, 63, 73, 72}, {63, 64, 74, 73},
{64, 65, 75, 74}, {65, 66, 76, 75}, {66, 67, 77, 76}, {68, 69, 78, 84},
{69, 70, 79, 78}, {70, 71, 85, 79}, {78, 79, 85, 84}, {71, 72, 80, 85},
{72, 73, 81, 80}, {73, 74, 86, 81}, {80, 81, 86, 85}, {74, 75, 82, 86},
{75, 76, 83, 82}, {76, 77, 87, 83}, {82, 83, 87, 86}, {84, 85, 89, 88},
{85, 86, 90, 89}, {86, 87, 91, 90}, {88, 89, 93, 92}, {89, 90, 94, 93},
{90, 91, 95, 94}};
size_t nelem_ = 81;
size_t nnode_ = 96;
size_t nne_ = 4;
size_t ndim_ = 2;
size_t shape_x_ = 9;
size_t shape_y_ = 21;
xt::xtensor<size_t, 1> elementsMiddleLayer_ = {36, 37, 38, 39, 40, 41, 42, 43, 44};
xt::xtensor<size_t, 1> nodesBottomEdge_ = {0, 1, 2, 3};
xt::xtensor<size_t, 1> nodesTopEdge_ = {92, 93, 94, 95};
xt::xtensor<size_t, 1> nodesLeftEdge_ = {0, 4, 8, 18, 28, 38, 48, 58, 68, 84, 88, 92};
xt::xtensor<size_t, 1> nodesRightEdge_ = {3, 7, 11, 27, 37, 47, 57, 67, 77, 87, 91, 95};
xt::xtensor<size_t, 1> nodesBottomOpenEdge_ = {1, 2};
xt::xtensor<size_t, 1> nodesTopOpenEdge_ = {93, 94};
xt::xtensor<size_t, 1> nodesLeftOpenEdge_ = {4, 8, 18, 28, 38, 48, 58, 68, 84, 88};
xt::xtensor<size_t, 1> nodesRightOpenEdge_ = {7, 11, 27, 37, 47, 57, 67, 77, 87, 91};
size_t nodesBottomLeftCorner_ = 0;
size_t nodesBottomRightCorner_ = 3;
size_t nodesTopLeftCorner_ = 92;
size_t nodesTopRightCorner_ = 95;
size_t nodesLeftBottomCorner_ = 0;
size_t nodesLeftTopCorner_ = 92;
size_t nodesRightBottomCorner_ = 3;
size_t nodesRightTopCorner_ = 95;
xt::xtensor<size_t, 2> dofs_ = {
{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11}, {12, 13},
{14, 15}, {16, 17}, {18, 19}, {20, 21}, {22, 23}, {24, 25}, {26, 27},
{28, 29}, {30, 31}, {32, 33}, {34, 35}, {36, 37}, {38, 39}, {40, 41},
{42, 43}, {44, 45}, {46, 47}, {48, 49}, {50, 51}, {52, 53}, {54, 55},
{56, 57}, {58, 59}, {60, 61}, {62, 63}, {64, 65}, {66, 67}, {68, 69},
{70, 71}, {72, 73}, {74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83},
{84, 85}, {86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97},
{98, 99}, {100, 101}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111},
{112, 113}, {114, 115}, {116, 117}, {118, 119}, {120, 121}, {122, 123}, {124, 125},
{126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137}, {138, 139},
{140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149}, {150, 151}, {152, 153},
{154, 155}, {156, 157}, {158, 159}, {160, 161}, {162, 163}, {164, 165}, {166, 167},
{168, 169}, {170, 171}, {172, 173}, {174, 175}, {176, 177}, {178, 179}, {180, 181},
{182, 183}, {184, 185}, {186, 187}, {188, 189}, {190, 191}};
xt::xtensor<size_t, 2> nodesPeriodic_ = {{0, 3},
{0, 95},
{0, 92},
{1, 93},
{2, 94},
{4, 7},
{8, 11},
{18, 27},
{28, 37},
{38, 47},
{48, 57},
{58, 67},
{68, 77},
{84, 87},
{88, 91}};
size_t nodesOrigin_ = 0;
xt::xtensor<size_t, 2> dofsPeriodic_ = {
{0, 1}, {2, 3}, {4, 5}, {0, 1}, {6, 7}, {8, 9}, {10, 11},
{6, 7}, {12, 13}, {14, 15}, {16, 17}, {12, 13}, {18, 19}, {20, 21},
{22, 23}, {24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35},
{36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}, {30, 31},
{48, 49}, {50, 51}, {52, 53}, {54, 55}, {56, 57}, {58, 59}, {60, 61},
{62, 63}, {64, 65}, {48, 49}, {66, 67}, {68, 69}, {70, 71}, {72, 73},
{74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83}, {66, 67}, {84, 85},
{86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97}, {98, 99},
{100, 101}, {84, 85}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111},
{112, 113}, {114, 115}, {116, 117}, {118, 119}, {102, 103}, {120, 121}, {122, 123},
{124, 125}, {126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137},
{120, 121}, {138, 139}, {140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149},
{150, 151}, {152, 153}, {154, 155}, {150, 151}, {156, 157}, {158, 159}, {160, 161},
{156, 157}, {0, 1}, {2, 3}, {4, 5}, {0, 1}};
GooseFEM::Mesh::Quad4::FineLayer mesh(9, 18);
xt::xtensor<double, 2> coor = mesh.coor();
xt::xtensor<size_t, 2> conn = mesh.conn();
size_t nelem = mesh.nelem();
size_t nnode = mesh.nnode();
size_t nne = mesh.nne();
size_t ndim = mesh.ndim();
size_t shape_x = mesh.nelx();
size_t shape_y = mesh.nely();
xt::xtensor<size_t, 1> elementsMiddleLayer = mesh.elementsMiddleLayer();
xt::xtensor<size_t, 1> nodesBottomEdge = mesh.nodesBottomEdge();
xt::xtensor<size_t, 1> nodesTopEdge = mesh.nodesTopEdge();
xt::xtensor<size_t, 1> nodesLeftEdge = mesh.nodesLeftEdge();
xt::xtensor<size_t, 1> nodesRightEdge = mesh.nodesRightEdge();
xt::xtensor<size_t, 1> nodesBottomOpenEdge = mesh.nodesBottomOpenEdge();
xt::xtensor<size_t, 1> nodesTopOpenEdge = mesh.nodesTopOpenEdge();
xt::xtensor<size_t, 1> nodesLeftOpenEdge = mesh.nodesLeftOpenEdge();
xt::xtensor<size_t, 1> nodesRightOpenEdge = mesh.nodesRightOpenEdge();
size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner();
size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner();
size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner();
size_t nodesTopRightCorner = mesh.nodesTopRightCorner();
size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner();
size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner();
size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner();
size_t nodesRightTopCorner = mesh.nodesRightTopCorner();
xt::xtensor<size_t, 2> dofs = mesh.dofs();
xt::xtensor<size_t, 2> nodesPeriodic = mesh.nodesPeriodic();
size_t nodesOrigin = mesh.nodesOrigin();
xt::xtensor<size_t, 2> dofsPeriodic = mesh.dofsPeriodic();
REQUIRE(nelem == nelem_);
REQUIRE(nnode == nnode_);
REQUIRE(nne == nne_);
REQUIRE(ndim == ndim_);
REQUIRE(shape_x == shape_x_);
REQUIRE(shape_y == shape_y_);
REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_);
REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_);
REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_);
REQUIRE(nodesTopRightCorner == nodesTopRightCorner_);
REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_);
REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_);
REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_);
REQUIRE(nodesRightTopCorner == nodesRightTopCorner_);
REQUIRE(nodesOrigin == nodesOrigin_);
REQUIRE(xt::allclose(coor, coor_));
REQUIRE(xt::all(xt::equal(elementsMiddleLayer, elementsMiddleLayer_)));
REQUIRE(xt::all(xt::equal(conn, conn_)));
REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_)));
REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_)));
REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_)));
}
+ SECTION("FineLayer - replica - trivial")
+ {
+ GooseFEM::Mesh::Quad4::FineLayer mesh(1, 1);
+ GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
+ REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
+ }
+
+ SECTION("FineLayer - replica - equi-sized")
+ {
+ GooseFEM::Mesh::Quad4::FineLayer mesh(4, 4);
+ GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
+ REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
+ }
+
+ SECTION("FineLayer - replica")
+ {
+ GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27);
+ GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
+ REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
+ REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
+ }
+
SECTION("Map::RefineRegular")
{
GooseFEM::Mesh::Quad4::Regular mesh(5, 4);
GooseFEM::Mesh::Quad4::Map::RefineRegular refine(mesh, 5, 3);
xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nelem()});
auto a_ = refine.mapToCoarse(refine.mapToFine(a));
REQUIRE(xt::allclose(a, xt::mean(a_, {1})));
xt::xtensor<double, 2> b =
xt::random::rand<double>(std::array<size_t, 2>{mesh.nelem(), 4ul});
auto b_ = refine.mapToCoarse(refine.mapToFine(b));
REQUIRE(xt::allclose(xt::mean(b, {1}), xt::mean(b_, {1})));
xt::xtensor<double, 4> c =
xt::random::rand<double>(std::array<size_t, 4>{mesh.nelem(), 4ul, 3ul, 3ul});
auto c_ = refine.mapToCoarse(refine.mapToFine(c));
REQUIRE(xt::allclose(xt::mean(c, {1}), xt::mean(c_, {1})));
}
SECTION("Map::FineLayer2Regular")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);
GooseFEM::Mesh::Quad4::Map::FineLayer2Regular map(mesh);
xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nelem()});
auto a_ = map.mapToRegular(a);
REQUIRE(xt::allclose(a, a_));
xt::xtensor<double, 2> b =
xt::random::rand<double>(std::array<size_t, 2>{mesh.nelem(), 4ul});
auto b_ = map.mapToRegular(b);
REQUIRE(xt::allclose(b, b_));
xt::xtensor<double, 4> c =
xt::random::rand<double>(std::array<size_t, 4>{mesh.nelem(), 4ul, 3ul, 3ul});
auto c_ = map.mapToRegular(c);
REQUIRE(xt::allclose(c, c_));
}
}