Page MenuHomec4science
Files F64486321

Mesh.hpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 27, 05:45

Mesh.hpp
View Options

/**
Implementation of Mesh.h
\file Mesh.hpp
\copyright Copyright 2017. Tom de Geus. All rights reserved.
\license This project is released under the GNU Public License (GPLv3).
*/
#ifndef GOOSEFEM_MESH_HPP
#define GOOSEFEM_MESH_HPP
#include "Mesh.h"
#include "assertions.h"
namespace GooseFEM {
namespace Mesh {
template <class S, class T>
inline ElementType defaultElementType(const S& coor, const T& conn)
{
GOOSEFEM_ASSERT(coor.dimension() == 2);
GOOSEFEM_ASSERT(conn.dimension() == 2);
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");
}
template <class D>
inline auto RegularBase<D>::nelem() const
{
return derived_cast().m_nelem;
}
template <class D>
inline auto RegularBase<D>::nnode() const
{
return derived_cast().m_nnode;
}
template <class D>
inline auto RegularBase<D>::nne() const
{
return derived_cast().m_nne;
}
template <class D>
inline auto RegularBase<D>::ndim() const
{
return derived_cast().m_ndim;
}
template <class D>
inline auto RegularBase<D>::nelx() const
{
return derived_cast().nelx_impl();
}
template <class D>
inline auto RegularBase<D>::nely() const
{
return derived_cast().nely_impl();
}
template <class D>
inline auto RegularBase<D>::h() const
{
return derived_cast().m_h;
}
template <class D>
inline auto RegularBase<D>::getElementType() const
{
return derived_cast().getElementType_impl();
}
template <class D>
inline auto RegularBase<D>::coor() const
{
return derived_cast().coor_impl();
}
template <class D>
inline auto RegularBase<D>::conn() const
{
return derived_cast().conn_impl();
}
template <class D>
inline auto RegularBase<D>::dofs() const
{
return GooseFEM::Mesh::dofs(this->nnode(), this->ndim());
}
template <class D>
inline auto RegularBase<D>::dofsPeriodic() const
{
array_type::tensor<size_t, 2> ret = this->dofs();
array_type::tensor<size_t, 2> nodePer = this->nodesPeriodic();
array_type::tensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
array_type::tensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
for (size_t j = 0; j < this->ndim(); ++j) {
xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
}
return GooseFEM::Mesh::renumber(ret);
}
template <class D>
inline auto RegularBase<D>::derived_cast() -> derived_type&
{
return *static_cast<derived_type*>(this);
}
template <class D>
inline auto RegularBase<D>::derived_cast() const -> const derived_type&
{
return *static_cast<const derived_type*>(this);
}
template <class D>
inline auto RegularBase<D>::nodesPeriodic() const
{
return derived_cast().nodesPeriodic_impl();
}
template <class D>
inline auto RegularBase<D>::nodesOrigin() const
{
return derived_cast().nodesOrigin_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesBottomEdge() const
{
return derived_cast().nodesBottomEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesTopEdge() const
{
return derived_cast().nodesTopEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesLeftEdge() const
{
return derived_cast().nodesLeftEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesRightEdge() const
{
return derived_cast().nodesRightEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesBottomOpenEdge() const
{
return derived_cast().nodesBottomOpenEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesTopOpenEdge() const
{
return derived_cast().nodesTopOpenEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesLeftOpenEdge() const
{
return derived_cast().nodesLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesRightOpenEdge() const
{
return derived_cast().nodesRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesBottomLeftCorner() const
{
return derived_cast().nodesBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesTopLeftCorner() const
{
return derived_cast().nodesTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesBottomRightCorner() const
{
return derived_cast().nodesBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesTopRightCorner() const
{
return derived_cast().nodesTopRightCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesLeftBottomCorner() const
{
return derived_cast().nodesBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesLeftTopCorner() const
{
return derived_cast().nodesTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesRightBottomCorner() const
{
return derived_cast().nodesBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::nodesRightTopCorner() const
{
return derived_cast().nodesTopRightCorner_impl();
}
template <class D>
inline array_type::tensor<size_t, 2> RegularBase2d<D>::nodesPeriodic_impl() const
{
array_type::tensor<size_t, 1> bot = derived_cast().nodesBottomOpenEdge_impl();
array_type::tensor<size_t, 1> top = derived_cast().nodesTopOpenEdge_impl();
array_type::tensor<size_t, 1> lft = derived_cast().nodesLeftOpenEdge_impl();
array_type::tensor<size_t, 1> rgt = derived_cast().nodesRightOpenEdge_impl();
std::array<size_t, 2> shape = {bot.size() + lft.size() + size_t(3), size_t(2)};
auto ret = array_type::tensor<size_t, 2>::from_shape(shape);
ret(0, 0) = derived_cast().nodesBottomLeftCorner_impl();
ret(0, 1) = derived_cast().nodesBottomRightCorner_impl();
ret(1, 0) = derived_cast().nodesBottomLeftCorner_impl();
ret(1, 1) = derived_cast().nodesTopRightCorner_impl();
ret(2, 0) = derived_cast().nodesBottomLeftCorner_impl();
ret(2, 1) = derived_cast().nodesTopLeftCorner_impl();
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;
}
template <class D>
inline auto RegularBase2d<D>::nodesOrigin_impl() const
{
return derived_cast().nodesBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase2d<D>::derived_cast() -> derived_type&
{
return *static_cast<derived_type*>(this);
}
template <class D>
inline auto RegularBase2d<D>::derived_cast() const -> const derived_type&
{
return *static_cast<const derived_type*>(this);
}
template <class D>
inline auto RegularBase3d<D>::nelz() const
{
return derived_cast().nelz_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottom() const
{
return derived_cast().nodesBottom_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTop() const
{
return derived_cast().nodesTop_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeft() const
{
return derived_cast().nodesLeft_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRight() const
{
return derived_cast().nodesRight_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFront() const
{
return derived_cast().nodesFront_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBack() const
{
return derived_cast().nodesBack_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontBottomEdge() const
{
return derived_cast().nodesFrontBottomEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontTopEdge() const
{
return derived_cast().nodesFrontTopEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontLeftEdge() const
{
return derived_cast().nodesFrontLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontRightEdge() const
{
return derived_cast().nodesFrontRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackBottomEdge() const
{
return derived_cast().nodesBackBottomEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackTopEdge() const
{
return derived_cast().nodesBackTopEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackLeftEdge() const
{
return derived_cast().nodesBackLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackRightEdge() const
{
return derived_cast().nodesBackRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomLeftEdge() const
{
return derived_cast().nodesBottomLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomRightEdge() const
{
return derived_cast().nodesBottomRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopLeftEdge() const
{
return derived_cast().nodesTopLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopRightEdge() const
{
return derived_cast().nodesTopRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomFrontEdge() const
{
return derived_cast().nodesFrontBottomEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomBackEdge() const
{
return derived_cast().nodesBackBottomEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopFrontEdge() const
{
return derived_cast().nodesFrontTopEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopBackEdge() const
{
return derived_cast().nodesBackTopEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBottomEdge() const
{
return derived_cast().nodesBottomLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftFrontEdge() const
{
return derived_cast().nodesFrontLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBackEdge() const
{
return derived_cast().nodesBackLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftTopEdge() const
{
return derived_cast().nodesTopLeftEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBottomEdge() const
{
return derived_cast().nodesBottomRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightTopEdge() const
{
return derived_cast().nodesTopRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightFrontEdge() const
{
return derived_cast().nodesFrontRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBackEdge() const
{
return derived_cast().nodesBackRightEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontFace() const
{
return derived_cast().nodesFrontFace_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackFace() const
{
return derived_cast().nodesBackFace_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftFace() const
{
return derived_cast().nodesLeftFace_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightFace() const
{
return derived_cast().nodesRightFace_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomFace() const
{
return derived_cast().nodesBottomFace_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopFace() const
{
return derived_cast().nodesTopFace_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontBottomOpenEdge() const
{
return derived_cast().nodesFrontBottomOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontTopOpenEdge() const
{
return derived_cast().nodesFrontTopOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontLeftOpenEdge() const
{
return derived_cast().nodesFrontLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontRightOpenEdge() const
{
return derived_cast().nodesFrontRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackBottomOpenEdge() const
{
return derived_cast().nodesBackBottomOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackTopOpenEdge() const
{
return derived_cast().nodesBackTopOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackLeftOpenEdge() const
{
return derived_cast().nodesBackLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackRightOpenEdge() const
{
return derived_cast().nodesBackRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomLeftOpenEdge() const
{
return derived_cast().nodesBottomLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomRightOpenEdge() const
{
return derived_cast().nodesBottomRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopLeftOpenEdge() const
{
return derived_cast().nodesTopLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopRightOpenEdge() const
{
return derived_cast().nodesTopRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomFrontOpenEdge() const
{
return derived_cast().nodesFrontBottomOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomBackOpenEdge() const
{
return derived_cast().nodesBackBottomOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopFrontOpenEdge() const
{
return derived_cast().nodesFrontTopOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopBackOpenEdge() const
{
return derived_cast().nodesBackTopOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBottomOpenEdge() const
{
return derived_cast().nodesBottomLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftFrontOpenEdge() const
{
return derived_cast().nodesFrontLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBackOpenEdge() const
{
return derived_cast().nodesBackLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftTopOpenEdge() const
{
return derived_cast().nodesTopLeftOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBottomOpenEdge() const
{
return derived_cast().nodesBottomRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightTopOpenEdge() const
{
return derived_cast().nodesTopRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightFrontOpenEdge() const
{
return derived_cast().nodesFrontRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBackOpenEdge() const
{
return derived_cast().nodesBackRightOpenEdge_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontBottomLeftCorner() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontBottomRightCorner() const
{
return derived_cast().nodesFrontBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontTopLeftCorner() const
{
return derived_cast().nodesFrontTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontTopRightCorner() const
{
return derived_cast().nodesFrontTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackBottomLeftCorner() const
{
return derived_cast().nodesBackBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackBottomRightCorner() const
{
return derived_cast().nodesBackBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackTopLeftCorner() const
{
return derived_cast().nodesBackTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackTopRightCorner() const
{
return derived_cast().nodesBackTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontLeftBottomCorner() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomFrontLeftCorner() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomLeftFrontCorner() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftFrontBottomCorner() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBottomFrontCorner() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontRightBottomCorner() const
{
return derived_cast().nodesFrontBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomFrontRightCorner() const
{
return derived_cast().nodesFrontBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomRightFrontCorner() const
{
return derived_cast().nodesFrontBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightFrontBottomCorner() const
{
return derived_cast().nodesFrontBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBottomFrontCorner() const
{
return derived_cast().nodesFrontBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontLeftTopCorner() const
{
return derived_cast().nodesFrontTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopFrontLeftCorner() const
{
return derived_cast().nodesFrontTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopLeftFrontCorner() const
{
return derived_cast().nodesFrontTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftFrontTopCorner() const
{
return derived_cast().nodesFrontTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftTopFrontCorner() const
{
return derived_cast().nodesFrontTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesFrontRightTopCorner() const
{
return derived_cast().nodesFrontTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopFrontRightCorner() const
{
return derived_cast().nodesFrontTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopRightFrontCorner() const
{
return derived_cast().nodesFrontTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightFrontTopCorner() const
{
return derived_cast().nodesFrontTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightTopFrontCorner() const
{
return derived_cast().nodesFrontTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackLeftBottomCorner() const
{
return derived_cast().nodesBackBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomBackLeftCorner() const
{
return derived_cast().nodesBackBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomLeftBackCorner() const
{
return derived_cast().nodesBackBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBackBottomCorner() const
{
return derived_cast().nodesBackBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBottomBackCorner() const
{
return derived_cast().nodesBackBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackRightBottomCorner() const
{
return derived_cast().nodesBackBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomBackRightCorner() const
{
return derived_cast().nodesBackBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBottomRightBackCorner() const
{
return derived_cast().nodesBackBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBackBottomCorner() const
{
return derived_cast().nodesBackBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBottomBackCorner() const
{
return derived_cast().nodesBackBottomRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackLeftTopCorner() const
{
return derived_cast().nodesBackTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopBackLeftCorner() const
{
return derived_cast().nodesBackTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopLeftBackCorner() const
{
return derived_cast().nodesBackTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftBackTopCorner() const
{
return derived_cast().nodesBackTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesLeftTopBackCorner() const
{
return derived_cast().nodesBackTopLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesBackRightTopCorner() const
{
return derived_cast().nodesBackTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopBackRightCorner() const
{
return derived_cast().nodesBackTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesTopRightBackCorner() const
{
return derived_cast().nodesBackTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightBackTopCorner() const
{
return derived_cast().nodesBackTopRightCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::nodesRightTopBackCorner() const
{
return derived_cast().nodesBackTopRightCorner_impl();
}
template <class D>
inline array_type::tensor<size_t, 2> RegularBase3d<D>::nodesPeriodic_impl() const
{
array_type::tensor<size_t, 1> fro = derived_cast().nodesFrontFace_impl();
array_type::tensor<size_t, 1> bck = derived_cast().nodesBackFace_impl();
array_type::tensor<size_t, 1> lft = derived_cast().nodesLeftFace_impl();
array_type::tensor<size_t, 1> rgt = derived_cast().nodesRightFace_impl();
array_type::tensor<size_t, 1> bot = derived_cast().nodesBottomFace_impl();
array_type::tensor<size_t, 1> top = derived_cast().nodesTopFace_impl();
array_type::tensor<size_t, 1> froBot = derived_cast().nodesFrontBottomOpenEdge_impl();
array_type::tensor<size_t, 1> froTop = derived_cast().nodesFrontTopOpenEdge_impl();
array_type::tensor<size_t, 1> froLft = derived_cast().nodesFrontLeftOpenEdge_impl();
array_type::tensor<size_t, 1> froRgt = derived_cast().nodesFrontRightOpenEdge_impl();
array_type::tensor<size_t, 1> bckBot = derived_cast().nodesBackBottomOpenEdge_impl();
array_type::tensor<size_t, 1> bckTop = derived_cast().nodesBackTopOpenEdge_impl();
array_type::tensor<size_t, 1> bckLft = derived_cast().nodesBackLeftOpenEdge_impl();
array_type::tensor<size_t, 1> bckRgt = derived_cast().nodesBackRightOpenEdge_impl();
array_type::tensor<size_t, 1> botLft = derived_cast().nodesBottomLeftOpenEdge_impl();
array_type::tensor<size_t, 1> botRgt = derived_cast().nodesBottomRightOpenEdge_impl();
array_type::tensor<size_t, 1> topLft = derived_cast().nodesTopLeftOpenEdge_impl();
array_type::tensor<size_t, 1> topRgt = derived_cast().nodesTopRightOpenEdge_impl();
size_t tface = fro.size() + lft.size() + bot.size();
size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size();
size_t tnode = 7;
array_type::tensor<size_t, 2> ret = xt::empty<size_t>({tface + tedge + tnode, std::size_t(2)});
size_t i = 0;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesFrontBottomRightCorner_impl();
++i;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesBackBottomRightCorner_impl();
++i;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesBackBottomLeftCorner_impl();
++i;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesFrontTopLeftCorner_impl();
++i;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesFrontTopRightCorner_impl();
++i;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesBackTopRightCorner_impl();
++i;
ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
ret(i, 1) = derived_cast().nodesBackTopLeftCorner_impl();
++i;
for (size_t j = 0; j < froBot.size(); ++j) {
ret(i, 0) = froBot(j);
ret(i, 1) = bckBot(j);
++i;
}
for (size_t j = 0; j < froBot.size(); ++j) {
ret(i, 0) = froBot(j);
ret(i, 1) = bckTop(j);
++i;
}
for (size_t j = 0; j < froBot.size(); ++j) {
ret(i, 0) = froBot(j);
ret(i, 1) = froTop(j);
++i;
}
for (size_t j = 0; j < botLft.size(); ++j) {
ret(i, 0) = botLft(j);
ret(i, 1) = botRgt(j);
++i;
}
for (size_t j = 0; j < botLft.size(); ++j) {
ret(i, 0) = botLft(j);
ret(i, 1) = topRgt(j);
++i;
}
for (size_t j = 0; j < botLft.size(); ++j) {
ret(i, 0) = botLft(j);
ret(i, 1) = topLft(j);
++i;
}
for (size_t j = 0; j < froLft.size(); ++j) {
ret(i, 0) = froLft(j);
ret(i, 1) = froRgt(j);
++i;
}
for (size_t j = 0; j < froLft.size(); ++j) {
ret(i, 0) = froLft(j);
ret(i, 1) = bckRgt(j);
++i;
}
for (size_t j = 0; j < froLft.size(); ++j) {
ret(i, 0) = froLft(j);
ret(i, 1) = bckLft(j);
++i;
}
for (size_t j = 0; j < fro.size(); ++j) {
ret(i, 0) = fro(j);
ret(i, 1) = bck(j);
++i;
}
for (size_t j = 0; j < lft.size(); ++j) {
ret(i, 0) = lft(j);
ret(i, 1) = rgt(j);
++i;
}
for (size_t j = 0; j < bot.size(); ++j) {
ret(i, 0) = bot(j);
ret(i, 1) = top(j);
++i;
}
return ret;
}
template <class D>
inline auto RegularBase3d<D>::nodesOrigin_impl() const
{
return derived_cast().nodesFrontBottomLeftCorner_impl();
}
template <class D>
inline auto RegularBase3d<D>::derived_cast() -> derived_type&
{
return *static_cast<derived_type*>(this);
}
template <class D>
inline auto RegularBase3d<D>::derived_cast() const -> const derived_type&
{
return *static_cast<const derived_type*>(this);
}
namespace detail {
template <class T, class R>
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
template <class S, class T>
inline array_type::tensor<size_t, 2>
overlapping(const S& coor_a, const T& coor_b, double rtol, double atol)
{
GOOSEFEM_ASSERT(coor_a.dimension() == 2);
GOOSEFEM_ASSERT(coor_b.dimension() == 2);
GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
std::vector<size_t> ret_a;
std::vector<size_t> 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);
}
}
array_type::tensor<size_t, 2> ret = xt::empty<size_t>({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;
}
template <class CA, class EA, class NA, class CB, class EB, class NB>
inline ManualStitch::ManualStitch(
const CA& coor_a,
const EA& conn_a,
const NA& overlapping_nodes_a,
const CB& coor_b,
const EB& conn_b,
const NB& overlapping_nodes_b,
bool check_position,
double rtol,
double atol)
{
UNUSED(rtol);
UNUSED(atol);
GOOSEFEM_ASSERT(coor_a.dimension() == 2);
GOOSEFEM_ASSERT(conn_a.dimension() == 2);
GOOSEFEM_ASSERT(overlapping_nodes_a.dimension() == 1);
GOOSEFEM_ASSERT(coor_b.dimension() == 2);
GOOSEFEM_ASSERT(conn_b.dimension() == 2);
GOOSEFEM_ASSERT(overlapping_nodes_b.dimension() == 1);
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_CHECK(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;
array_type::tensor<size_t, 1> keep_b =
xt::setdiff1d(xt::arange<size_t>(nndb), overlapping_nodes_b);
m_map_b = xt::empty<size_t>({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<size_t>(keep_b.size()) + nnda;
m_conn = xt::empty<size_t>({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<size_t>({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 array_type::tensor<double, 2> ManualStitch::coor() const
{
return m_coor;
}
inline array_type::tensor<size_t, 2> ManualStitch::conn() const
{
return m_conn;
}
inline size_t ManualStitch::nmesh() const
{
return 2;
}
inline size_t ManualStitch::nelem() const
{
return m_conn.shape(0);
}
inline size_t ManualStitch::nnode() const
{
return m_coor.shape(0);
}
inline size_t ManualStitch::nne() const
{
return m_conn.shape(1);
}
inline size_t ManualStitch::ndim() const
{
return m_coor.shape(1);
}
inline array_type::tensor<size_t, 2> ManualStitch::dofs() const
{
size_t nnode = this->nnode();
size_t ndim = this->ndim();
return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
}
inline std::vector<array_type::tensor<size_t, 1>> ManualStitch::nodemap() const
{
std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
for (size_t i = 0; i < this->nmesh(); ++i) {
ret[i] = this->nodemap(i);
}
return ret;
}
inline std::vector<array_type::tensor<size_t, 1>> ManualStitch::elemmap() const
{
std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
for (size_t i = 0; i < this->nmesh(); ++i) {
ret[i] = this->elemmap(i);
}
return ret;
}
inline array_type::tensor<size_t, 1> ManualStitch::nodemap(size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index <= 1);
if (mesh_index == 0) {
return xt::arange<size_t>(m_nnd_a);
}
return m_map_b;
}
inline array_type::tensor<size_t, 1> ManualStitch::elemmap(size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index <= 1);
if (mesh_index == 0) {
return xt::arange<size_t>(m_nel_a);
}
return xt::arange<size_t>(m_nel_b) + m_nel_a;
}
template <class T>
inline T ManualStitch::nodeset(const T& set, size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index <= 1);
if (mesh_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);
}
template <class T>
inline T ManualStitch::elemset(const T& set, size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index <= 1);
if (mesh_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;
}
template <class C, class E>
inline void Stitch::push_back(const C& coor, const E& conn)
{
GOOSEFEM_ASSERT(coor.dimension() == 2);
GOOSEFEM_ASSERT(conn.dimension() == 2);
if (m_map.size() == 0) {
m_coor = coor;
m_conn = conn;
m_map.push_back(xt::eval(xt::arange<size_t>(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 stitch(
m_coor,
m_conn,
xt::eval(xt::view(overlap, 0, xt::all())),
coor,
conn,
xt::eval(xt::view(overlap, 1, xt::all())),
false);
m_coor = stitch.coor();
m_conn = stitch.conn();
m_map.push_back(stitch.nodemap(1));
m_nel.push_back(conn.shape(0));
m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
}
inline size_t Stitch::nmesh() const
{
return m_map.size();
}
inline array_type::tensor<double, 2> Stitch::coor() const
{
return m_coor;
}
inline array_type::tensor<size_t, 2> Stitch::conn() const
{
return m_conn;
}
inline size_t Stitch::nelem() const
{
return m_conn.shape(0);
}
inline size_t Stitch::nnode() const
{
return m_coor.shape(0);
}
inline size_t Stitch::nne() const
{
return m_conn.shape(1);
}
inline size_t Stitch::ndim() const
{
return m_coor.shape(1);
}
inline array_type::tensor<size_t, 2> Stitch::dofs() const
{
size_t nnode = this->nnode();
size_t ndim = this->ndim();
return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
}
inline std::vector<array_type::tensor<size_t, 1>> Stitch::nodemap() const
{
std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
for (size_t i = 0; i < this->nmesh(); ++i) {
ret[i] = this->nodemap(i);
}
return ret;
}
inline std::vector<array_type::tensor<size_t, 1>> Stitch::elemmap() const
{
std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
for (size_t i = 0; i < this->nmesh(); ++i) {
ret[i] = this->elemmap(i);
}
return ret;
}
inline array_type::tensor<size_t, 1> Stitch::nodemap(size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index < m_map.size());
return m_map[mesh_index];
}
inline array_type::tensor<size_t, 1> Stitch::elemmap(size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index < m_map.size());
return xt::arange<size_t>(m_nel[mesh_index]) + m_el_offset[mesh_index];
}
template <class T>
inline T Stitch::nodeset(const T& set, size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index < m_map.size());
GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size());
return detail::renum(set, m_map[mesh_index]);
}
template <class T>
inline T Stitch::elemset(const T& set, size_t mesh_index) const
{
GOOSEFEM_ASSERT(mesh_index < m_map.size());
GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]);
return set + m_el_offset[mesh_index];
}
template <class T>
inline T Stitch::nodeset(const std::vector<T>& 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();
}
array_type::tensor<size_t, 1> ret = xt::empty<size_t>({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);
}
template <class T>
inline T Stitch::elemset(const std::vector<T>& 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();
}
array_type::tensor<size_t, 1> ret = xt::empty<size_t>({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;
}
template <class T>
inline T Stitch::nodeset(std::initializer_list<T> set) const
{
return this->nodeset(std::vector<T>(set));
}
template <class T>
inline T Stitch::elemset(std::initializer_list<T> set) const
{
return this->elemset(std::vector<T>(set));
}
inline Vstack::Vstack(bool check_overlap, double rtol, double atol)
{
m_check_overlap = check_overlap;
m_rtol = rtol;
m_atol = atol;
}
template <class C, class E, class N>
inline void Vstack::push_back(const C& coor, const E& conn, const N& nodes_bot, const N& nodes_top)
{
if (m_map.size() == 0) {
m_coor = coor;
m_conn = conn;
m_map.push_back(xt::eval(xt::arange<size_t>(coor.shape(0))));
m_nel.push_back(conn.shape(0));
m_el_offset.push_back(0);
m_nodes_bot.push_back(nodes_bot);
m_nodes_top.push_back(nodes_top);
return;
}
GOOSEFEM_ASSERT(nodes_bot.size() == m_nodes_top.back().size());
size_t index = m_map.size();
double shift = xt::amax(xt::view(m_coor, xt::all(), 1))();
auto x = coor;
xt::view(x, xt::all(), 1) += shift;
ManualStitch stitch(
m_coor, m_conn, m_nodes_top.back(), x, conn, nodes_bot, m_check_overlap, m_rtol, m_atol);
m_nodes_bot.push_back(stitch.nodeset(nodes_bot, 1));
m_nodes_top.push_back(stitch.nodeset(nodes_top, 1));
m_coor = stitch.coor();
m_conn = stitch.conn();
m_map.push_back(stitch.nodemap(1));
m_nel.push_back(conn.shape(0));
m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
}
template <class T>
inline Renumber::Renumber(const T& dofs)
{
size_t n = xt::amax(dofs)() + 1;
size_t i = 0;
array_type::tensor<size_t, 1> unique = xt::unique(dofs);
m_renum = xt::empty<size_t>({n});
for (auto& j : unique) {
m_renum(j) = i;
++i;
}
}
template <class T>
inline T Renumber::apply(const T& list) const
{
return detail::renum(list, m_renum);
}
inline array_type::tensor<size_t, 1> Renumber::index() const
{
return m_renum;
}
template <class T>
inline T renumber(const T& dofs)
{
return Renumber(dofs).apply(dofs);
}
template <class T>
inline Reorder::Reorder(const std::initializer_list<T> 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(is_unique(arg));
}
#endif
m_renum = xt::empty<size_t>({n});
for (auto& arg : args) {
for (auto& j : arg) {
m_renum(j) = i;
++i;
}
}
}
template <class T>
inline 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 array_type::tensor<size_t, 1> Reorder::index() const
{
return m_renum;
}
inline array_type::tensor<size_t, 2> dofs(size_t nnode, size_t ndim)
{
return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
}
template <class D>
inline std::vector<std::vector<size_t>> nodaltyings(const D& dofs)
{
size_t nnode = dofs.shape(0);
size_t ndim = dofs.shape(1);
auto nodemap = node2dof(dofs);
std::vector<std::vector<size_t>> ret(nnode);
for (size_t m = 0; m < nnode; ++m) {
auto r = nodemap[dofs(m, 0)];
std::sort(r.begin(), r.end());
ret[m] = r;
#ifdef GOOSEFEM_ENABLE_ASSERT
for (size_t i = 1; i < ndim; ++i) {
auto t = nodemap[dofs(m, i)];
std::sort(t.begin(), t.end());
GOOSEFEM_ASSERT(std::equal(r.begin(), r.end(), t.begin()));
}
#endif
}
return ret;
}
template <class E>
inline array_type::tensor<size_t, 1> coordination(const E& conn)
{
GOOSEFEM_ASSERT(conn.dimension() == 2);
size_t nnode = xt::amax(conn)() + 1;
array_type::tensor<size_t, 1> N = xt::zeros<size_t>({nnode});
for (auto it = conn.begin(); it != conn.end(); ++it) {
N(*it) += 1;
}
return N;
}
template <class D>
inline std::vector<std::vector<size_t>> node2dof(const D& dofs, bool sorted)
{
return elem2node(dofs, sorted);
}
template <class E>
inline std::vector<std::vector<size_t>> elem2node(const E& conn, bool sorted)
{
auto N = coordination(conn);
auto nnode = N.size();
std::vector<std::vector<size_t>> ret(nnode);
for (size_t i = 0; i < nnode; ++i) {
ret[i].reserve(N(i));
}
for (size_t e = 0; e < conn.shape(0); ++e) {
for (size_t m = 0; m < conn.shape(1); ++m) {
ret[conn(e, m)].push_back(e);
}
}
if (sorted) {
for (auto& row : ret) {
std::sort(row.begin(), row.end());
}
}
return ret;
}
template <class E, class D>
inline std::vector<std::vector<size_t>> elem2node(const E& conn, const D& dofs, bool sorted)
{
size_t nnode = dofs.shape(0);
auto ret = elem2node(conn, sorted);
auto nties = nodaltyings(dofs);
for (size_t m = 0; m < nnode; ++m) {
if (nties[m].size() <= 1) {
continue;
}
if (m > nties[m][0]) {
continue;
}
size_t n = ret[m].size();
for (size_t j = 1; j < nties[m].size(); ++j) {
n += ret[nties[m][j]].size();
}
ret[m].reserve(n);
for (size_t j = 1; j < nties[m].size(); ++j) {
ret[m].insert(ret[m].end(), ret[nties[m][j]].begin(), ret[nties[m][j]].end());
}
if (sorted) {
std::sort(ret[m].begin(), ret[m].end());
}
for (size_t j = 1; j < nties[m].size(); ++j) {
ret[nties[m][j]] = ret[m];
}
}
return ret;
}
template <class C, class E>
inline array_type::tensor<double, 2> edgesize(const C& coor, const E& conn, ElementType type)
{
GOOSEFEM_ASSERT(coor.dimension() == 2);
GOOSEFEM_ASSERT(conn.dimension() == 2);
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
if (type == ElementType::Quad4) {
GOOSEFEM_ASSERT(coor.shape(1) == 2ul);
GOOSEFEM_ASSERT(conn.shape(1) == 4ul);
array_type::tensor<size_t, 1> n0 = xt::view(conn, xt::all(), 0);
array_type::tensor<size_t, 1> n1 = xt::view(conn, xt::all(), 1);
array_type::tensor<size_t, 1> n2 = xt::view(conn, xt::all(), 2);
array_type::tensor<size_t, 1> n3 = xt::view(conn, xt::all(), 3);
array_type::tensor<double, 1> x0 = xt::view(coor, xt::keep(n0), 0);
array_type::tensor<double, 1> x1 = xt::view(coor, xt::keep(n1), 0);
array_type::tensor<double, 1> x2 = xt::view(coor, xt::keep(n2), 0);
array_type::tensor<double, 1> x3 = xt::view(coor, xt::keep(n3), 0);
array_type::tensor<double, 1> y0 = xt::view(coor, xt::keep(n0), 1);
array_type::tensor<double, 1> y1 = xt::view(coor, xt::keep(n1), 1);
array_type::tensor<double, 1> y2 = xt::view(coor, xt::keep(n2), 1);
array_type::tensor<double, 1> y3 = xt::view(coor, xt::keep(n3), 1);
array_type::tensor<double, 2> ret = xt::empty<double>(conn.shape());
xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0));
xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0));
xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0));
xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0));
return ret;
}
throw std::runtime_error("Element-type not implemented");
}
template <class C, class E>
inline array_type::tensor<double, 2> edgesize(const C& coor, const E& conn)
{
return edgesize(coor, conn, defaultElementType(coor, conn));
}
template <class C, class E>
inline array_type::tensor<double, 2> centers(const C& coor, const E& conn, ElementType type)
{
GOOSEFEM_ASSERT(coor.dimension() == 2);
GOOSEFEM_ASSERT(conn.dimension() == 2);
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
array_type::tensor<double, 2> ret = xt::zeros<double>({conn.shape(0), coor.shape(1)});
if (type == ElementType::Quad4) {
GOOSEFEM_ASSERT(coor.shape(1) == 2);
GOOSEFEM_ASSERT(conn.shape(1) == 4);
for (size_t i = 0; i < 4; ++i) {
auto n = xt::view(conn, xt::all(), i);
ret += xt::view(coor, xt::keep(n), xt::all());
}
ret /= 4.0;
return ret;
}
throw std::runtime_error("Element-type not implemented");
}
template <class C, class E>
inline array_type::tensor<double, 2> centers(const C& coor, const E& conn)
{
return centers(coor, conn, defaultElementType(coor, conn));
}
template <class T, class C, class E>
inline array_type::tensor<size_t, 1>
elemmap2nodemap(const T& elem_map, const C& coor, const E& conn, ElementType type)
{
GOOSEFEM_ASSERT(elem_map.dimension() == 1);
GOOSEFEM_ASSERT(coor.dimension() == 2);
GOOSEFEM_ASSERT(conn.dimension() == 2);
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0));
size_t N = coor.shape(0);
array_type::tensor<size_t, 1> ret = N * xt::ones<size_t>({N});
if (type == ElementType::Quad4) {
GOOSEFEM_ASSERT(coor.shape(1) == 2);
GOOSEFEM_ASSERT(conn.shape(1) == 4);
for (size_t i = 0; i < 4; ++i) {
array_type::tensor<size_t, 1> t = N * xt::ones<size_t>({N});
auto old_nd = xt::view(conn, xt::all(), i);
auto new_nd = xt::view(conn, xt::keep(elem_map), i);
xt::view(t, xt::keep(old_nd)) = new_nd;
ret = xt::where(xt::equal(ret, N), t, ret);
}
return ret;
}
throw std::runtime_error("Element-type not implemented");
}
template <class T, class C, class E>
inline array_type::tensor<size_t, 1>
elemmap2nodemap(const T& elem_map, const C& coor, const E& conn)
{
return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn));
}
} // namespace Mesh
} // namespace GooseFEM
#endif

Event Timeline

c4science · Help