Page MenuHomec4science
Files F61330485

Element.hpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 00:09

Element.hpp
View Options

/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_ELEMENT_CPP
#define GOOSEFEM_ELEMENT_CPP
// -------------------------------------------------------------------------------------------------
#include "Element.h"
// =================================================================================================
namespace GooseFEM {
namespace Element {
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,3> asElementVector(
const xt::xtensor<size_t,2> &conn, const xt::xtensor<double,2> &nodevec)
{
// extract dimensions
size_t nelem = conn .shape()[0]; // number of elements
size_t nne = conn .shape()[1]; // number of nodes per element
size_t ndim = nodevec.shape()[1]; // number of dimensions
// allocate output: nodal vectors stored per element
xt::xtensor<double,3> elemvec = xt::empty<double>({nelem, nne, ndim});
// read from nodal vectors
#pragma omp parallel for
for ( size_t e = 0 ; e < nelem ; ++e )
for ( size_t m = 0 ; m < nne ; ++m )
for ( size_t i = 0 ; i < ndim ; ++i )
elemvec(e,m,i) = nodevec(conn(e,m),i);
return elemvec;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,2> assembleNodeVector(
const xt::xtensor<size_t,2> &conn, const xt::xtensor<double,3> &elemvec)
{
// extract dimensions
size_t nelem = conn.shape()[0]; // number of elements
size_t nne = conn.shape()[1]; // number of nodes per element
size_t ndim = elemvec.shape()[2]; // number of dimensions
size_t nnode = xt::amax(conn)[0]+1; // number of nodes
// check input
assert( elemvec.shape()[0] == nelem );
assert( elemvec.shape()[1] == nne );
// zero-initialize output: nodal vectors
xt::xtensor<double,2> nodevec = xt::zeros<double>({nnode, ndim});
// assemble from nodal vectors stored per element
for ( size_t e = 0 ; e < nelem ; ++e )
for ( size_t m = 0 ; m < nne ; ++m )
for ( size_t i = 0 ; i < ndim ; ++i )
nodevec(conn(e,m),i) += elemvec(e,m,i);
return nodevec;
}
// -------------------------------------------------------------------------------------------------
template<class E>
inline bool isSequential(const E &dofs)
{
// number of DOFs
size_t ndof = xt::amax(dofs)[0] + 1;
// list to check if all DOFs are present
// - allocate
xt::xtensor<int,1> exists = xt::zeros<int>({ndof});
// - fill
for ( auto &i : dofs ) exists[i]++;
// check
for ( auto &i : dofs )
if ( exists[i] == 0 )
return false;
// checks passed successfully
return true;
}
// -------------------------------------------------------------------------------------------------
inline bool isDiagonal(const xt::xtensor<double,3> &elemmat)
{
// check input
assert( elemmat.shape()[1] == elemmat.shape()[2] );
// get dimensions
size_t nelem = elemmat.shape()[0];
size_t N = elemmat.shape()[1];
// get numerical precision
double eps = std::numeric_limits<double>::epsilon();
// loop over all entries
#pragma omp parallel for
for ( size_t e = 0 ; e < nelem ; ++e )
for ( size_t i = 0 ; i < N ; ++i )
for ( size_t j = 0 ; j < N ; ++j )
if ( i != j )
if ( std::abs(elemmat(e,i,j)) > eps )
return false;
// all checks passed successful
return true;
}
// -------------------------------------------------------------------------------------------------
}} // namespace ...
// =================================================================================================
#endif

Event Timeline

c4science · Help