diff --git a/docs/details/figures/MeshQuad4/FineLayer/stitch.py b/docs/details/figures/MeshQuad4/FineLayer/stitch.py index f3297f2..a007d04 100644 --- a/docs/details/figures/MeshQuad4/FineLayer/stitch.py +++ b/docs/details/figures/MeshQuad4/FineLayer/stitch.py @@ -1,44 +1,41 @@ import GooseFEM import matplotlib.pyplot as plt import GooseMPL as gplt import numpy as np plt.style.use(['goose']) fig, ax = plt.subplots() mesh = GooseFEM.Mesh.Quad4.FineLayer(6 * 9, 51) -coor_a = mesh.coor() -conn_a = mesh.conn() +coor_0 = mesh.coor() +coor_1 = mesh.coor() +coor_2 = mesh.coor() -coor_b = mesh.coor() -conn_b = mesh.conn() +H = np.max(coor_0[:, 1]) -coor_c = mesh.coor() -conn_c = mesh.conn() - -coor_b[:, 1] += np.max(coor_a[:, 1]) -coor_c[:, 1] += np.max(coor_a[:, 1]) * 2 +coor_1[:, 1] += H +coor_2[:, 1] += H * 2 stitch = GooseFEM.Mesh.Stitch() -stitch.push_back(coor_a, conn_a) -stitch.push_back(coor_b, conn_b) -stitch.push_back(coor_c, conn_c) +stitch.push_back(coor_0, mesh.conn()) +stitch.push_back(coor_1, mesh.conn()) +stitch.push_back(coor_2, mesh.conn()) coor = stitch.coor() conn = stitch.conn() cindex = np.zeros(conn.shape[0]) cindex[stitch.elemset(np.arange(mesh.nelem()), 0)] = 1 cindex[stitch.elemset(np.arange(mesh.nelem()), 1)] = 2 cindex[stitch.elemset(np.arange(mesh.nelem()), 2)] = 3 im = gplt.patch(coor=coor, conn=conn, cindex=cindex, cmap='jet', axis=ax) ax.set_aspect('equal') ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show() diff --git a/docs/details/figures/MeshQuad4/FineLayer/stitch_hybrid.py b/docs/details/figures/MeshQuad4/FineLayer/stitch_hybrid.py new file mode 100644 index 0000000..95a26e2 --- /dev/null +++ b/docs/details/figures/MeshQuad4/FineLayer/stitch_hybrid.py @@ -0,0 +1,59 @@ + +import GooseFEM +import matplotlib.pyplot as plt +import GooseMPL as gplt +import numpy as np + +plt.style.use(['goose']) + +fig, ax = plt.subplots() + +finelayer = GooseFEM.Mesh.Quad4.FineLayer(6 * 9, 51) + +coor_l0 = finelayer.coor() +coor_l1 = finelayer.coor() +coor_l2 = finelayer.coor() + +Hl = np.max(coor_l0[:, 1]) + +coor_l1[:, 1] += Hl +coor_l2[:, 1] += Hl * 2 + +h = finelayer.elemrow_nhx()[0] +nx = finelayer.elemrow_nelem()[0] +ny = np.ceil(((finelayer.nely() - 1) / 2) / h) + +regular = GooseFEM.Mesh.Quad4.Regular(int(nx), int(ny), float(h)) + +coor_r0 = regular.coor() +coor_r1 = regular.coor() + +Hr = np.max(coor_r0[:, 1]) + +coor_r0[:, 1] -= Hr +coor_r1[:, 1] += 3 * Hl + +stitch = GooseFEM.Mesh.Stitch() +stitch.push_back(coor_r0, regular.conn()) +stitch.push_back(coor_l0, finelayer.conn()) +stitch.push_back(coor_l1, finelayer.conn()) +stitch.push_back(coor_l2, finelayer.conn()) +stitch.push_back(coor_r1, regular.conn()) + +coor = stitch.coor() +conn = stitch.conn() + +cindex = np.zeros(conn.shape[0]) +cindex[stitch.elemset(np.arange(regular.nelem()), 0)] = 1 +cindex[stitch.elemset(np.arange(finelayer.nelem()), 1)] = 2 +cindex[stitch.elemset(np.arange(finelayer.nelem()), 2)] = 3 +cindex[stitch.elemset(np.arange(finelayer.nelem()), 3)] = 4 +cindex[stitch.elemset(np.arange(regular.nelem()), 4)] = 5 + +im = gplt.patch(coor=coor, conn=conn, cindex=cindex, cmap='jet', axis=ax) + +ax.set_aspect('equal') +ax.get_xaxis().set_visible(False) +ax.get_yaxis().set_visible(False) + +plt.show() diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp index cc97568..f2aef98 100644 --- a/include/GooseFEM/Mesh.hpp +++ b/include/GooseFEM/Mesh.hpp @@ -1,542 +1,552 @@ /* (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& coor, const xt::xtensor& 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"); } namespace detail { template inline T renum(const T& arg, const R& mapping) { T ret = T::from_shape(arg.shape()); auto jt = ret.begin(); for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) { *jt = mapping(*it); } return ret; } } // namespace detail inline ManualStitch::ManualStitch( const xt::xtensor& coor_a, const xt::xtensor& conn_a, const xt::xtensor& overlapping_nodes_a, const xt::xtensor& coor_b, const xt::xtensor& conn_b, const xt::xtensor& overlapping_nodes_b, bool check_position, double rtol, double atol) { UNUSED(rtol); UNUSED(atol); GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape())); GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1)); if (check_position) { GOOSEFEM_ASSERT(xt::allclose( xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()), xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()), rtol, atol)); } size_t nnda = coor_a.shape(0); size_t nndb = coor_b.shape(0); size_t ndim = coor_a.shape(1); size_t nelim = overlapping_nodes_a.size(); size_t nela = conn_a.shape(0); size_t nelb = conn_b.shape(0); size_t nne = conn_a.shape(1); m_nel_a = nela; m_nel_b = nelb; m_nnd_a = nnda; xt::xtensor keep_b = xt::setdiff1d(xt::arange(nndb), overlapping_nodes_b); m_map_b = xt::empty({nndb}); xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a; xt::view(m_map_b, xt::keep(keep_b)) = xt::arange(keep_b.size()) + nnda; m_conn = xt::empty({nela + nelb, nne}); xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a; xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b); m_coor = xt::empty({nnda + nndb - nelim, ndim}); xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a; xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) = xt::view(coor_b, xt::keep(keep_b), xt::all()); } inline xt::xtensor ManualStitch::coor() const { return m_coor; } inline xt::xtensor ManualStitch::conn() const { return m_conn; } inline xt::xtensor ManualStitch::nodemap(size_t index) const { GOOSEFEM_ASSERT(index <= 1); if (index == 0) { return xt::arange(m_nnd_a); } return m_map_b; } inline xt::xtensor ManualStitch::elemmap(size_t index) const { GOOSEFEM_ASSERT(index <= 1); if (index == 0) { return xt::arange(m_nel_a); } return xt::arange(m_nel_b) + m_nel_a; } inline xt::xtensor ManualStitch::nodeset(const xt::xtensor& set, size_t index) const { GOOSEFEM_ASSERT(index <= 1); if (index == 0) { + GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a); return set; } + GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size()); return detail::renum(set, m_map_b); } inline xt::xtensor ManualStitch::elemset(const xt::xtensor& set, size_t index) const { GOOSEFEM_ASSERT(index <= 1); if (index == 0) { + GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a); return set; } + GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b); return set + m_nel_a; } inline Stitch::Stitch(double rtol, double atol) { m_rtol = rtol; m_atol = atol; } inline void Stitch::push_back(const xt::xtensor& coor, const xt::xtensor& conn) { if (m_map.size() == 0) { m_coor = coor; m_conn = conn; m_map.push_back(xt::eval(xt::arange(coor.shape(0)))); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(0); return; } auto overlap = overlapping(m_coor, coor, m_rtol, m_atol); size_t index = m_map.size(); ManualStitch stich( m_coor, m_conn, xt::view(overlap, 0, xt::all()), coor, conn, xt::view(overlap, 1, xt::all()), false); m_coor = stich.coor(); m_conn = stich.conn(); m_map.push_back(stich.nodemap(1)); m_nel.push_back(conn.shape(0)); - m_el_offset.push_back(m_el_offset[index - 1] + conn.shape(0)); + m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]); } inline xt::xtensor Stitch::coor() const { return m_coor; } inline xt::xtensor Stitch::conn() const { return m_conn; } inline xt::xtensor Stitch::nodemap(size_t index) const { GOOSEFEM_ASSERT(index < m_map.size()); return m_map[index]; } inline xt::xtensor Stitch::elemmap(size_t index) const { GOOSEFEM_ASSERT(index < m_map.size()); return xt::arange(m_nel[index]) + m_el_offset[index]; } inline xt::xtensor Stitch::nodeset(const xt::xtensor& set, size_t index) const { GOOSEFEM_ASSERT(index < m_map.size()); + GOOSEFEM_ASSERT(xt::amax(set)() < m_map[index].size()); return detail::renum(set, m_map[index]); } inline xt::xtensor Stitch::elemset(const xt::xtensor& set, size_t index) const { GOOSEFEM_ASSERT(index < m_map.size()); + GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[index]); return set + m_el_offset[index]; } inline xt::xtensor Stitch::nodeset(const std::vector>& set) const { + GOOSEFEM_ASSERT(set.size() == m_map.size()); + size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->nodeset(set[i], i); n += set[i].size(); } return xt::unique(ret); } inline xt::xtensor Stitch::elemset(const std::vector>& set) const { + GOOSEFEM_ASSERT(set.size() == m_map.size()); + size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i); n += set[i].size(); } return ret; } inline Renumber::Renumber(const xt::xarray& dofs) { size_t n = xt::amax(dofs)() + 1; size_t i = 0; xt::xtensor unique = xt::unique(dofs); m_renum = xt::empty({n}); for (auto& j : unique) { m_renum(j) = i; ++i; } } // ret(i,j) = renum(list(i,j)) template T Renumber::apply(const T& list) const { return detail::renum(list, m_renum); } inline xt::xtensor Renumber::get(const xt::xtensor& dofs) const { return this->apply(dofs); } inline xt::xtensor Renumber::index() const { return m_renum; } inline Reorder::Reorder(const std::initializer_list> 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({n}); for (auto& arg : args) { for (auto& j : arg) { m_renum(j) = i; ++i; } } } inline xt::xtensor Reorder::get(const xt::xtensor& dofs) const { return this->apply(dofs); } inline xt::xtensor Reorder::index() const { return m_renum; } // apply renumbering, e.g. for a matrix: // // ret(i,j) = renum(list(i,j)) template 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 renumber(const xt::xtensor& dofs) { return Renumber(dofs).get(dofs); } inline xt::xtensor dofs(size_t nnode, size_t ndim) { return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline xt::xtensor coordination(const xt::xtensor& conn) { size_t nnode = xt::amax(conn)() + 1; xt::xtensor N = xt::zeros({nnode}); for (auto it = conn.begin(); it != conn.end(); ++it) { N(*it) += 1; } return N; } inline std::vector> elem2node(const xt::xtensor& conn, bool sorted) { auto N = coordination(conn); auto nnode = N.size(); std::vector> 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 edgesize( const xt::xtensor& coor, const xt::xtensor& 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 n0 = xt::view(conn, xt::all(), 0); xt::xtensor n1 = xt::view(conn, xt::all(), 1); xt::xtensor n2 = xt::view(conn, xt::all(), 2); xt::xtensor n3 = xt::view(conn, xt::all(), 3); xt::xtensor x0 = xt::view(coor, xt::keep(n0), 0); xt::xtensor x1 = xt::view(coor, xt::keep(n1), 0); xt::xtensor x2 = xt::view(coor, xt::keep(n2), 0); xt::xtensor x3 = xt::view(coor, xt::keep(n3), 0); xt::xtensor y0 = xt::view(coor, xt::keep(n0), 1); xt::xtensor y1 = xt::view(coor, xt::keep(n1), 1); xt::xtensor y2 = xt::view(coor, xt::keep(n2), 1); xt::xtensor y3 = xt::view(coor, xt::keep(n3), 1); xt::xtensor ret = xt::empty(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 edgesize( const xt::xtensor& coor, const xt::xtensor& conn) { return edgesize(coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); xt::xtensor ret = xt::zeros({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 centers( const xt::xtensor& coor, const xt::xtensor& conn) { return centers(coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& 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 ret = N * xt::ones({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 t = N * xt::ones({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 elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn) { return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor overlapping( const xt::xtensor& coor_a, const xt::xtensor& coor_b, double rtol, double atol) { GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); std::vector ret_a; std::vector ret_b; for (size_t i = 0; i < coor_a.shape(0); ++i) { auto idx = xt::flatten_indices(xt::argwhere(xt::prod(xt::isclose( coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1))); for (auto& j : idx) { ret_a.push_back(i); ret_b.push_back(j); } } xt::xtensor ret = xt::empty({size_t(2), ret_a.size()}); for (size_t i = 0; i < ret_a.size(); ++i) { ret(0, i) = ret_a[i]; ret(1, i) = ret_b[i]; } return ret; } } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h index ccb0a77..d013098 100644 --- a/include/GooseFEM/MeshQuad4.h +++ b/include/GooseFEM/MeshQuad4.h @@ -1,286 +1,292 @@ /* (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 { // pre-allocation namespace Map { class FineLayer2Regular; } // Regular mesh: equi-sized elements 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 coor() const; // nodal positions [nnode, ndim] xt::xtensor conn() const; // connectivity [nelem, nne] // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor 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 dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; // element numbers as matrix xt::xtensor elementgrid() 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 }; // Mesh with fine middle layer, and coarser elements towards the top and bottom 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& coor, const xt::xtensor& 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 + // edge size, per row of elements (in units of "h") + xt::xtensor elemrow_nhx() const; + xt::xtensor elemrow_nhy() const; + xt::xtensor elemrow_nelem() const; + // type ElementType getElementType() const; // mesh xt::xtensor coor() const; // nodal positions [nnode, ndim] xt::xtensor conn() const; // connectivity [nelem, nne] // elements in the middle (fine) layer xt::xtensor elementsMiddleLayer() const; // select region of elements from 'matrix' of element numbers xt::xtensor elementgrid_ravel( - std::vector rows_range, - std::vector cols_range) const; + std::vector rows_start_stop, + std::vector cols_start_stop) const; // select region of elements from 'matrix' of element numbers around an element + // - square box with edge-size (2 * "size" + 1) * "h", around "element" xt::xtensor elementgrid_around_ravel( size_t element, size_t size, bool periodic = true); - + // - left/right from "element" (on the same layer) xt::xtensor elementgrid_leftright( size_t element, size_t left, size_t right, bool periodic = true); // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor 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 dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; // mapping to 'roll' periodically in the x-direction, // returns element mapping, such that: new_elemvar = elemvar[elem_map] xt::xtensor 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 m_nelx; // number of elements in "x" (*) xt::xtensor m_nnd; // total number of nodes in the main node layer (**) xt::xtensor m_nhx; // element size in x-direction (*) xt::xtensor m_nhy; // element size in y-direction (*) xt::xtensor m_refine; // refine direction (-1:no refine, 0:"x" (*) xt::xtensor m_startElem; // start element (*) xt::xtensor 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& coor, const xt::xtensor& conn); friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular; }; // Mesh mappings namespace Map { // Return "FineLayer"-class responsible for generating a connectivity // Throws if conversion is not possible GooseFEM::Mesh::Quad4::FineLayer FineLayer( const xt::xtensor& coor, const xt::xtensor& conn); // Refine a regular mesh: sub-divide elements in several smaller elements class RefineRegular { public: // Constructors RefineRegular() = default; RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny); // return the coarse or the fine mesh objects 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 getMap() const; // map field xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToCoarse(const xt::xtensor& data) const; // tensor per intpnt // map field xt::xtensor mapToFine(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToFine(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToFine(const xt::xtensor& data) const; // tensor per intpnt private: // the meshes GooseFEM::Mesh::Quad4::Regular m_coarse; GooseFEM::Mesh::Quad4::Regular m_fine; // mapping xt::xtensor m_fine2coarse; xt::xtensor m_fine2coarse_index; xt::xtensor 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> getMap() const; std::vector> getMapFraction() const; // map field xt::xtensor mapToRegular(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToRegular(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToRegular(const xt::xtensor& 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> m_elem_regular; std::vector> 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 949fd8a..7d605d1 100644 --- a/include/GooseFEM/MeshQuad4.hpp +++ b/include/GooseFEM/MeshQuad4.hpp @@ -1,1621 +1,1636 @@ /* (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 Regular::coor() const { xt::xtensor ret = xt::empty({m_nnode, m_ndim}); xt::xtensor x = xt::linspace(0.0, m_h * static_cast(m_nelx), m_nelx + 1); xt::xtensor y = xt::linspace(0.0, m_h * static_cast(m_nely), m_nely + 1); size_t inode = 0; for (size_t iy = 0; iy < m_nely + 1; ++iy) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ++inode; } } return ret; } inline xt::xtensor Regular::conn() const { xt::xtensor ret = xt::empty({m_nelem, m_nne}); size_t ielem = 0; for (size_t iy = 0; iy < m_nely; ++iy) { for (size_t ix = 0; ix < m_nelx; ++ix) { 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 Regular::nodesBottomEdge() const { return xt::arange(m_nelx + 1); } inline xt::xtensor Regular::nodesTopEdge() const { return xt::arange(m_nelx + 1) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1) + m_nelx; } inline xt::xtensor Regular::nodesBottomOpenEdge() const { return xt::arange(1, m_nelx); } inline xt::xtensor Regular::nodesTopOpenEdge() const { return xt::arange(1, m_nelx) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftOpenEdge() const { return xt::arange(1, m_nely) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightOpenEdge() const { return xt::arange(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 Regular::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(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 Regular::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor Regular::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor 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 Regular::elementgrid() const { return xt::arange(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& coor, const xt::xtensor& 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(nelx); // compute element size in y-direction (use symmetry, compute upper half) // temporary variables size_t nmin, ntot; xt::xtensor nhx = xt::ones({nely}); xt::xtensor nhy = xt::ones({nely}); xt::xtensor refine = -1 * xt::ones({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({nely * 2 - 1}); m_nhy = xt::empty({nely * 2 - 1}); m_refine = xt::empty({nely * 2 - 1}); m_nelx = xt::empty({nely * 2 - 1}); m_nnd = xt::empty({nely * 2}); m_startElem = xt::empty({nely * 2 - 1}); m_startNode = xt::empty({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 xt::xtensor FineLayer::elemrow_nhx() const +{ + return m_nhx; +} + +inline xt::xtensor FineLayer::elemrow_nhy() const +{ + return m_nhy; +} + +inline xt::xtensor FineLayer::elemrow_nelem() const +{ + return m_nelx; +} + inline ElementType FineLayer::getElementType() const { return ElementType::Quad4; } inline xt::xtensor FineLayer::coor() const { // allocate output xt::xtensor ret = xt::empty({m_nnode, m_ndim}); // current node, number of element layers size_t inode = 0; size_t nely = static_cast(m_nhy.size()); // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) // - allocate xt::xtensor y = xt::empty({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 x = xt::linspace(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(m_nhx(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t ix = 0; ix < m_nelx(iy); ++ix) { for (size_t j = 0; j < 2; ++j) { ret(inode, 0) = x(ix) + dx * static_cast(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 x = xt::linspace(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(m_nhx(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t ix = 0; ix < m_nelx(iy); ++ix) { for (size_t j = 0; j < 2; ++j) { ret(inode, 0) = x(ix) + dx * static_cast(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 FineLayer::conn() const { // allocate output xt::xtensor ret = xt::empty({m_nelem, m_nne}); // current element, number of element layers, starting nodes of each node layer size_t ielem = 0; size_t nely = static_cast(m_nhy.size()); size_t bot, mid, top; // loop over all element layers for (size_t iy = 0; iy < nely; ++iy) { // - get: starting nodes of bottom(, middle) and top layer bot = m_startNode(iy); mid = m_startNode(iy) + m_nnd(iy); top = m_startNode(iy + 1); // - define connectivity: no coarsening/refinement if (m_refine(iy) == -1) { for (size_t 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 FineLayer::elementsMiddleLayer() const { size_t nely = m_nhy.size(); size_t iy = (nely - 1) / 2; return m_startElem(iy) + xt::arange(m_nelx(iy)); } inline xt::xtensor FineLayer::elementgrid_ravel( std::vector start_stop_rows, std::vector start_stop_cols) const { GOOSEFEM_ASSERT(start_stop_rows.size() == 0 || start_stop_rows.size() == 2); GOOSEFEM_ASSERT(start_stop_cols.size() == 0 || start_stop_cols.size() == 2); std::array rows; std::array cols; if (start_stop_rows.size() == 2) { std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin()); GOOSEFEM_ASSERT(rows[0] <= this->nely()); GOOSEFEM_ASSERT(rows[1] <= this->nely()); } else { rows[0] = 0; rows[1] = this->nely(); } if (start_stop_cols.size() == 2) { std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin()); GOOSEFEM_ASSERT(cols[0] <= this->nelx()); GOOSEFEM_ASSERT(cols[1] <= this->nelx()); } else { cols[0] = 0; cols[1] = this->nelx(); } if (rows[0] == rows[1] || cols[0] == cols[1]) { xt::xtensor ret = xt::empty({0}); return ret; } // Compute dimensions auto H = xt::cumsum(m_nhy); size_t yl = 0; if (rows[0] > 0) { yl = xt::argmax(H > rows[0])(); } size_t yu = xt::argmax(H >= rows[1])(); size_t hx = std::max(m_nhx(yl), m_nhx(yu)); size_t xl = (cols[0] - cols[0] % hx) / hx; size_t xu = (cols[1] - cols[1] % hx) / hx; // Allocate output size_t N = 0; for (size_t i = yl; i <= yu; ++i) { // no refinement size_t n = (xu - xl) * hx / m_nhx(i); // refinement if (m_refine(i) != -1) { n *= 4; } N += n; } xt::xtensor ret = xt::empty({N}); // Write output N = 0; for (size_t i = yl; i <= yu; ++i) { // no refinement size_t n = (xu - xl) * hx / m_nhx(i); size_t h = hx; // refinement if (m_refine(i) != -1) { n *= 4; h *= 4; } xt::xtensor e = m_startElem(i) + xl * h / m_nhx(i) + xt::arange(n); xt::view(ret, xt::range(N, N + n)) = e; N += n; } return ret; } inline xt::xtensor FineLayer::elementgrid_around_ravel( size_t e, size_t size, bool periodic) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_nelx(iy); GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2); auto H = xt::cumsum(m_nhy); if (2 * size >= H(H.size() - 1)) { return xt::arange(this->nelem()); } size_t hy = H(iy); size_t l = xt::argmax(H > (hy - size - 1))(); size_t u = xt::argmax(H >= (hy + size))(); size_t lh = 0; if (l > 0) { lh = H(l - 1); } size_t uh = H(u); size_t step = xt::amax(m_nhx)(); size_t relx = (e - m_startElem(iy)) % step; size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx; size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step; size_t dx = m_nhx(u); size_t xl = 0; size_t xu = nel; if (mid >= size) { xl = mid - size; } if (mid + size < nel) { xu = mid + size + 1; } xl = xl - xl % dx; xu = xu - xu % dx; if (mid - xl < size) { if (xl < dx) { xl = 0; } else { xl -= dx; } } if (xu - mid < size) { if (xu > nel - dx) { xu = nel; } else { xu += dx; } } auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu}); auto map = this->roll(nroll); return xt::view(map, xt::keep(ret)); } inline xt::xtensor FineLayer::elementgrid_leftright( size_t e, size_t left, size_t right, bool periodic) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_nelx(iy); GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2); size_t step = xt::amax(m_nhx)(); size_t relx = (e - m_startElem(iy)) % step; size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx; size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step; size_t dx = m_nhx(iy); size_t xl = 0; size_t xu = nel; if (mid >= left) { xl = mid - left; } if (mid + right < nel) { xu = mid + right + 1; } xl = xl - xl % dx; xu = xu - xu % dx; if (mid - xl < left) { if (xl < dx) { xl = 0; } else { xl -= dx; } } if (xu - mid < right) { if (xu > nel - dx) { xu = nel; } else { xu += dx; } } auto H = xt::cumsum(m_nhy); size_t yl = 0; if (iy > 0) { yl = H(iy - 1); } auto ret = this->elementgrid_ravel({yl, H(iy)}, {xl, xu}); auto map = this->roll(nroll); return xt::view(map, xt::keep(ret)); } inline xt::xtensor FineLayer::nodesBottomEdge() const { return m_startNode(0) + xt::arange(m_nelx(0) + 1); } inline xt::xtensor FineLayer::nodesTopEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(m_nelx(nely - 1) + 1); } inline xt::xtensor FineLayer::nodesLeftEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesRightEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesBottomOpenEdge() const { return m_startNode(0) + xt::arange(1, m_nelx(0)); } inline xt::xtensor FineLayer::nodesTopOpenEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(1, m_nelx(nely - 1)); } inline xt::xtensor FineLayer::nodesLeftOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesRightOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(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 FineLayer::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor FineLayer::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor 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 FineLayer::roll(size_t n) { auto conn = this->conn(); size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({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(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& coor, const xt::xtensor& 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)); if (conn.shape(0) == 1) { this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0)); GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn))); GOOSEFEM_CHECK(xt::allclose(this->coor(), coor)); return; } // Identify the middle layer 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)); // Get element sizes along the middle layer auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0); auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1); auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 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(dx)(); auto hy = xt::amin(dy)(); GOOSEFEM_CHECK(xt::allclose(hx, hy)); GOOSEFEM_CHECK(xt::allclose(dx, hx)); GOOSEFEM_CHECK(xt::allclose(dy, hy)); // Extract shape and initialise size_t nelx = eright - eleft + 1; size_t nely = static_cast((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(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 elmat_coarse = m_coarse.elementgrid(); xt::xtensor elmat_fine = m_fine.elementgrid(); m_coarse2fine = xt::empty({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 RefineRegular::getMap() const { return m_coarse2fine; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& 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 ret = xt::empty({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 RefineRegular::mapToCoarse(const xt::xtensor& 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 ret = xt::empty({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 RefineRegular::mapToCoarse(const xt::xtensor& 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 ret = xt::empty({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 RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({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 RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({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 RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({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 nhx = m_finelayer.m_nhx; xt::xtensor nhy = m_finelayer.m_nhy; xt::xtensor nelx = m_finelayer.m_nelx; xt::xtensor start = m_finelayer.m_startElem; // 'matrix' of element numbers of the Regular-mesh xt::xtensor elementgrid = m_regular.elementgrid(); // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh xt::xtensor cum_nhy = xt::concatenate(xt::xtuple(xt::zeros({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(elementgrid, 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 el_old = start(iy) + xt::arange(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 el_old = start(iy) + xt::arange(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 el_old = start(iy) + xt::arange(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> FineLayer2Regular::getMap() const { return m_elem_regular; } inline std::vector> FineLayer2Regular::getMapFraction() const { return m_frac_regular; } inline xt::xtensor FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); xt::xtensor ret = xt::zeros({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 FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); xt::xtensor ret = xt::zeros({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 FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); xt::xtensor ret = xt::zeros({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 02cdb54..ea429a4 100644 --- a/python/MeshQuad4.hpp +++ b/python/MeshQuad4.hpp @@ -1,255 +1,228 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include namespace py = pybind11; void init_MeshQuad4(py::module& m) { py::class_(m, "Regular") .def( py::init(), "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge " "size 'h'", py::arg("nx"), py::arg("ny"), py::arg("h") = 1.) .def("coor", &GooseFEM::Mesh::Quad4::Regular::coor) - .def("conn", &GooseFEM::Mesh::Quad4::Regular::conn) - .def("nelem", &GooseFEM::Mesh::Quad4::Regular::nelem) - .def("nnode", &GooseFEM::Mesh::Quad4::Regular::nnode) - .def("nne", &GooseFEM::Mesh::Quad4::Regular::nne) - .def("ndim", &GooseFEM::Mesh::Quad4::Regular::ndim) - .def("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("elementgrid", &GooseFEM::Mesh::Quad4::Regular::elementgrid) .def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) { return ""; }); py::class_(m, "FineLayer") .def( py::init(), "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical " "direction", py::arg("nx"), py::arg("ny"), py::arg("h") = 1., py::arg("nfine") = 1) .def( py::init&, const xt::xtensor&>(), "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("elemrow_nhx", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhx) + .def("elemrow_nhy", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhy) + .def("elemrow_nelem", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nelem) .def("getElementType", &GooseFEM::Mesh::Quad4::FineLayer::getElementType) .def("elementsMiddleLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer) .def( "elementgrid_ravel", &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_ravel, py::arg("rows_range"), py::arg("cols_range")) .def( "elementgrid_around_ravel", &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_around_ravel, py::arg("element"), py::arg("size"), py::arg("periodic") = true) .def( "elementgrid_leftright", &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_leftright, py::arg("element"), py::arg("left"), py::arg("right"), py::arg("periodic") = true) .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 ""; }); } void init_MeshQuad4Map(py::module& m) { py::class_(m, "RefineRegular") .def( py::init(), "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&>( &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_)) .def( "mapToCoarse", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_)) .def( "mapToCoarse", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_)) .def( "mapToFine", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_)) .def( "mapToFine", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_)) .def( "mapToFine", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_)) .def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::RefineRegular&) { return ""; }); py::class_(m, "FineLayer2Regular") .def( py::init(), "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&>( &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_)) .def( "mapToRegular", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_)) .def( "mapToRegular", py::overload_cast&>( &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_)) .def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular&) { return ""; }); }