diff --git a/include/GooseFEM/Allocate.h b/include/GooseFEM/Allocate.h index 5095de3..8c86809 100644 --- a/include/GooseFEM/Allocate.h +++ b/include/GooseFEM/Allocate.h @@ -1,247 +1,247 @@ /** * Common allocation methods. * * @file * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ALLOCATE_H #define GOOSEFEM_ALLOCATE_H #include "config.h" namespace GooseFEM { template inline void asTensor(const T& arg, R& ret); namespace detail { /** * Check that two shapes partly overlap. If `s` has more dimensions that `t` the excess dimensions * of `s` are ignored and the first `t.size()` dimensions are checked for equality. */ template inline bool has_shape_begin(const T& t, const S& s) { return s.dimension() >= t.dimension() && std::equal(t.shape().cbegin(), t.shape().cend(), s.shape().begin()); } /** * Static identification of an std::array */ template struct is_std_array : std::false_type { }; template struct is_std_array> : std::true_type { }; /** * Helper for std_array_size */ template auto std_array_size_impl(const std::array&) -> std::integral_constant; /** * Get the size of an std:array (T::size is not static) */ template using std_array_size = decltype(std_array_size_impl(std::declval())); /** * Return as std::array. */ template std::array to_std_array(const I (&shape)[L]) { std::array r; std::copy(&shape[0], &shape[0] + L, r.begin()); return r; } /** * asTensor for array_type::array. */ template struct asTensor_write { static void impl(const T& arg, R& ret) { GOOSEFEM_ASSERT(arg.dimension() <= ret.dimension()); GOOSEFEM_ASSERT(detail::has_shape_begin(arg, ret)); using strides_type = typename T::strides_type::value_type; std::vector ret_strides(ret.dimension()); std::copy(arg.strides().begin(), arg.strides().end(), ret_strides.begin()); std::fill(ret_strides.begin() + arg.dimension(), ret_strides.end(), 0); ret = xt::strided_view( arg, ret.shape(), std::move(ret_strides), 0ul, xt::layout_type::dynamic); } }; /** * asTensor for xt::tensor. */ template struct asTensor_write< T, R, typename std::enable_if_t::value && xt::has_fixed_rank_t::value>> { static void impl(const T& arg, R& ret) { static_assert(T::rank <= R::rank, "Return must be fixed rank too"); GOOSEFEM_ASSERT(detail::has_shape_begin(arg, ret)); using strides_type = typename T::strides_type::value_type; std::array ret_strides; std::copy(arg.strides().begin(), arg.strides().end(), ret_strides.begin()); std::fill(ret_strides.begin() + T::rank, ret_strides.end(), 0); ret = xt::strided_view( arg, ret.shape(), std::move(ret_strides), 0ul, xt::layout_type::dynamic); } }; /** * AsTensor for array_type::array. */ template struct asTensor_allocate { static auto impl(const T& arg, const S& shape) { using value_type = typename T::value_type; size_t dim = arg.dimension(); size_t rank = shape.size(); std::vector ret_shape(dim + rank); std::copy(arg.shape().begin(), arg.shape().end(), ret_shape.begin()); std::copy(shape.begin(), shape.end(), ret_shape.begin() + dim); xt::xarray ret(ret_shape); GooseFEM::asTensor(arg, ret); return ret; } }; /** * AsTensor for xt::tensor. */ template struct asTensor_allocate::value>> { static auto impl(const T& arg, const S& shape) { using value_type = typename T::value_type; static constexpr size_t dim = T::rank; static constexpr size_t rank = std_array_size::value; std::array ret_shape; std::copy(arg.shape().begin(), arg.shape().end(), ret_shape.begin()); std::copy(shape.begin(), shape.end(), ret_shape.begin() + dim); array_type::tensor ret = xt::empty(ret_shape); detail::asTensor_write, decltype(ret)>::impl(arg, ret); return ret; } }; } // namespace detail /** - * "Broadcast" a scalar stored in an array (e.g. ``[r, s]``) to the same scalar of all - * tensor components of a tensor of certain rank (e.g. for rank 2: ``[r, s, i, j]``). + * "Broadcast" a scalar stored in an array (e.g. `[r, s]`) to the same scalar of all + * tensor components of a tensor of certain rank (e.g. for rank 2: `[r, s, i, j]`). * * @param arg An array with scalars. * @param ret Corresponding array with tensors. */ template inline void asTensor(const T& arg, R& ret) { detail::asTensor_write, std::decay_t>::impl(arg, ret); } /** - * "Broadcast" a scalar stored in an array (e.g. ``[r, s]``) to the same scalar of all - * tensor components of a tensor of certain rank (e.g. for rank 2: ``[r, s, i, j]``). + * "Broadcast" a scalar stored in an array (e.g. `[r, s]`) to the same scalar of all + * tensor components of a tensor of certain rank (e.g. for rank 2: `[r, s, i, j]`). * * @param arg An array with scalars. - * @param shape The shape of the added tensor dimensions (e.g.: ``[i, j]``). + * @param shape The shape of the added tensor dimensions (e.g.: `[i, j]`). * @return Corresponding array with tensors. */ template inline auto AsTensor(const T& arg, const S& shape) { return detail::asTensor_allocate, std::decay_t>::impl(arg, shape); } /** * @copydoc AsTensor(const T& arg, const S& shape) */ template inline auto AsTensor(const T& arg, const I (&shape)[L]) { auto s = detail::to_std_array(shape); return detail::asTensor_allocate, decltype(s)>::impl(arg, s); } /** - * "Broadcast" a scalar stored in an array (e.g. ``[r, s]``) to the same scalar of all - * tensor components of a tensor of certain rank (e.g. for rank 2: ``[r, s, n, n]``). + * "Broadcast" a scalar stored in an array (e.g. `[r, s]`) to the same scalar of all + * tensor components of a tensor of certain rank (e.g. for rank 2: `[r, s, n, n]`). * * @tparam rank Number of tensor dimensions (number of dimensions to add to the input). * @param arg An array with scalars. * @param n The shape along each of the added dimensions. * @return Corresponding array with tensors. */ template inline auto AsTensor(const T& arg, size_t n) { std::array shape; std::fill(shape.begin(), shape.end(), n); return detail::asTensor_allocate, decltype(shape)>::impl(arg, shape); } /** - * "Broadcast" a scalar stored in an array (e.g. ``[r, s]``) to the same scalar of all - * tensor components of a tensor of certain rank (e.g. for rank 2: ``[r, s, n, n]``). + * "Broadcast" a scalar stored in an array (e.g. `[r, s]`) to the same scalar of all + * tensor components of a tensor of certain rank (e.g. for rank 2: `[r, s, n, n]`). * * @param rank Number of tensor dimensions (number of dimensions to add to the input). * @param arg An array with scalars. * @param n The shape along each of the added dimensions. * @return Corresponding array with tensors. */ template inline auto AsTensor(size_t rank, const T& arg, size_t n) { std::vector shape(rank); std::fill(shape.begin(), shape.end(), n); return detail::asTensor_allocate, decltype(shape)>::impl(arg, shape); } /** - * Zero-pad columns to a matrix until is that shape ``[m, 3]``. + * Zero-pad columns to a matrix until is that shape `[m, 3]`. * * @param arg A "nodevec" (``arg.shape(1) <= 3``). * @return Corresponding "nodevec" in 3-d (``ret.shape(1) == 3``) */ template inline T as3d(const T& arg) { GOOSEFEM_ASSERT(arg.dimension() == 2); GOOSEFEM_ASSERT(arg.shape(1) > 0 && arg.shape(1) < 4); if (arg.shape(1) == 3ul) { return arg; } T ret = xt::zeros(std::array{arg.shape(0), 3ul}); if (arg.shape(1) == 2ul) { xt::view(ret, xt::all(), xt::keep(0, 1)) = arg; } if (arg.shape(1) == 1ul) { xt::view(ret, xt::all(), xt::keep(0)) = arg; } return ret; } } // namespace GooseFEM #endif diff --git a/include/GooseFEM/Element.h b/include/GooseFEM/Element.h index 519db26..97479d4 100644 --- a/include/GooseFEM/Element.h +++ b/include/GooseFEM/Element.h @@ -1,1229 +1,1229 @@ /** * Convenience methods for integration point data. * * @file Element.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENT_H #define GOOSEFEM_ELEMENT_H #include "Allocate.h" #include "config.h" #include "detail.h" namespace GooseFEM { /** * Element quadrature and interpolation. */ namespace Element { /** - * Convert nodal vector with ("nodevec", shape:``[nnode, ndim]``) to nodal vector stored per element - * ("elemvec", shape: ``[nelem, nne, ndim]``). + * Convert nodal vector with ("nodevec", shape:`[nnode, ndim]`) to nodal vector stored per element + * ("elemvec", shape: `[nelem, nne, ndim]`). * * @param conn Connectivity. * @param nodevec "nodevec". * @return "elemvec". */ inline array_type::tensor asElementVector( const array_type::tensor& conn, const array_type::tensor& nodevec) { size_t nelem = conn.shape(0); size_t nne = conn.shape(1); size_t ndim = nodevec.shape(1); array_type::tensor elemvec = xt::empty({nelem, nne, ndim}); #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; } /** - * Assemble nodal vector stored per element ("elemvec", shape ``[nelem, nne, ndim]``) to nodal + * Assemble nodal vector stored per element ("elemvec", shape `[nelem, nne, ndim]`) to nodal * vector - * ("nodevec", shape ``[nnode, ndim]``). + * ("nodevec", shape `[nnode, ndim]`). * * @param conn Connectivity. * @param elemvec "elemvec". * @return "nodevec". */ inline array_type::tensor assembleNodeVector( const array_type::tensor& conn, const array_type::tensor& elemvec) { size_t nelem = conn.shape(0); size_t nne = conn.shape(1); size_t ndim = elemvec.shape(2); size_t nnode = xt::amax(conn)() + 1; GOOSEFEM_ASSERT(elemvec.shape(0) == nelem); GOOSEFEM_ASSERT(elemvec.shape(1) == nne); array_type::tensor nodevec = xt::zeros({nnode, ndim}); 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; } /** * Check that DOFs leave no holes. * * @param dofs DOFs ("nodevec") - * @return ``true`` if there are no holds. + * @return `true` if there are no holds. */ template inline bool isSequential(const E& dofs) { size_t ndof = xt::amax(dofs)() + 1; array_type::tensor exists = xt::zeros({ndof}); for (auto& i : dofs) { exists[i]++; } for (auto& i : dofs) { if (exists[i] == 0) { return false; } } return true; } /** - * Check that all of the matrices stored per elemmat (shape: ``[nelem, nne * ndim, nne * ndim]``) + * Check that all of the matrices stored per elemmat (shape: `[nelem, nne * ndim, nne * ndim]`) * are diagonal. * * @param elemmat Element-vectors ("elemmat") - * @return ``true`` if all element matrices are diagonal. + * @return `true` if all element matrices are diagonal. */ bool isDiagonal(const array_type::tensor& elemmat) { GOOSEFEM_ASSERT(elemmat.shape(1) == elemmat.shape(2)); size_t nelem = elemmat.shape(0); size_t N = elemmat.shape(1); double eps = std::numeric_limits::epsilon(); #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; } } } } } return true; } /** * CRTP base class for quadrature. */ template class QuadratureBase { public: /** * Underlying type. */ using derived_type = D; /** * Number of elements. * * @return Scalar. */ auto nelem() const { return derived_cast().m_nelem; } /** * Number of nodes per element. * * @return Scalar. */ auto nne() const { return D::s_nne; } /** * Number of dimensions for node vectors. * * @return Scalar. */ auto ndim() const { return D::s_ndim; } /** * Number of dimensions for integration point tensors. * * @return Scalar. */ auto tdim() const { return D::s_tdim; } /** * Number of integration points. * * @return Scalar. */ auto nip() const { return derived_cast().m_nip; } /** * Convert "qscalar" to "qtensor" of certain rank. * Fully allocated output passed as reference, use AsTensor to allocate and return data. * * @param arg A "qscalar". * @param ret A "qtensor". */ template void asTensor(const T& arg, R& ret) const { GOOSEFEM_ASSERT(xt::has_shape(arg, this->shape_qscalar())); GooseFEM::asTensor(arg, ret); } /** * Convert "qscalar" to "qtensor" of certain rank. * * @param arg A "qscalar". * @return "qtensor". */ template auto AsTensor(const T& arg) const { return GooseFEM::AsTensor(arg, derived_cast().m_tdim); } /** * Convert "qscalar" to "qtensor" of certain rank. * * @param rank Tensor rank. * @param arg A "qscalar". * @return "qtensor". */ template auto AsTensor(size_t rank, const T& arg) const { return GooseFEM::AsTensor(rank, arg, derived_cast().m_tdim); } /** * Get the shape of an "elemvec". * * @returns [#nelem, #nne, #ndim]. */ auto shape_elemvec() const -> std::array { return std::array{derived_cast().m_nelem, D::s_nne, D::s_ndim}; } /** * Get the shape of an "elemvec". * * @param arg The vector dimension. * @returns [#nelem, #nne, tdim]. */ auto shape_elemvec(size_t arg) const -> std::array { return std::array{derived_cast().m_nelem, D::s_nne, arg}; } /** * Get the shape of an "elemmat". * * @returns [#nelem, #nne * #ndim, #nne * #ndim]. */ auto shape_elemmat() const -> std::array { return std::array{ derived_cast().m_nelem, D::s_nne * D::s_ndim, D::s_nne * D::s_ndim}; } /** * Get the shape of a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * * @tparam rank Rank of the tensor. Output is fixed-size: `std::array`. * @returns [#nelem, #nip, #tdim, ...]. */ template auto shape_qtensor() const -> std::array { std::array shape; shape[0] = derived_cast().m_nelem; shape[1] = derived_cast().m_nip; std::fill(shape.begin() + 2, shape.end(), derived_cast().m_tdim); return shape; } /** * Get the shape of a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * * @param rank The tensor rank. * @returns [#nelem, #nip, #tdim, ...]. */ auto shape_qtensor(size_t rank) const -> std::vector { std::vector shape(2 + rank); shape[0] = derived_cast().m_nelem; shape[1] = derived_cast().m_nip; std::fill(shape.begin() + 2, shape.end(), derived_cast().m_tdim); return shape; } /** * Get the shape of a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * * @tparam rank Rank of the tensor. Output is fixed-size: `std::array`. * @param rank The tensor rank. Effectively useless: to distinguish from the dynamic-sized. * @param arg The tensor dimension. * @returns [#nelem, #nip, tdim, ...]. */ template auto shape_qtensor(size_t rank, size_t arg) const -> std::array { GOOSEFEM_ASSERT(trank == rank); std::array shape; shape[0] = derived_cast().m_nelem; shape[1] = derived_cast().m_nip; std::fill(shape.begin() + 2, shape.end(), arg); return shape; } /** * Get the shape of a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * * @param rank The tensor rank. * @param arg The tensor dimension. * @returns [#nelem, #nip, tdim, ...]. */ auto shape_qtensor(size_t rank, size_t arg) const -> std::vector { std::vector shape(2 + rank); shape[0] = derived_cast().m_nelem; shape[1] = derived_cast().m_nip; std::fill(shape.begin() + 2, shape.end(), arg); return shape; } /** * Get the shape of a "qscalar" (a "qtensor" of rank 0) * @returns [#nelem, #nip]. */ auto shape_qscalar() const -> std::array { return std::array{derived_cast().m_nelem, derived_cast().m_nip}; } /** * Get the shape of a "qvector" (a "qtensor" of rank 1) * @returns [#nelem, #nip, #tdim]. */ auto shape_qvector() const -> std::array { return std::array{derived_cast().m_nelem, derived_cast().m_nip, D::s_tdim}; } /** * Get the shape of a "qvector" (a "qtensor" of rank 1) * @param arg Tensor dimension. * @returns [#nelem, #nip, #tdim]. */ auto shape_qvector(size_t arg) const -> std::array { return std::array{derived_cast().m_nelem, derived_cast().m_nip, arg}; } /** * Get an allocated `array_type::tensor` to store a "elemvec". * Note: the container is not (zero-)initialised. * * @tparam R value-type of the array, e.g. `double`. * @returns `xt::xarray` container of the correct shape. */ template auto allocate_elemvec() const { return array_type::tensor::from_shape(this->shape_elemvec()); } /** * Get an allocated and initialised `xt::xarray` to store a "elemvec". * * @tparam R value-type of the array, e.g. `double`. * @param val The value to which to initialise all items. * @returns `array_type::tensor` container of the correct shape. */ template auto allocate_elemvec(R val) const { auto ret = array_type::tensor::from_shape(this->shape_elemvec()); ret.fill(val); return ret; } /** * Get an allocated `array_type::tensor` to store a "elemmat". * Note: the container is not (zero-)initialised. * * @tparam R value-type of the array, e.g. `double`. * @returns `xt::xarray` container of the correct shape. */ template auto allocate_elemmat() const { return array_type::tensor::from_shape(this->shape_elemmat()); } /** * Get an allocated and initialised `xt::xarray` to store a "elemmat". * * @tparam R value-type of the array, e.g. `double`. * @param val The value to which to initialise all items. * @returns `array_type::tensor` container of the correct shape. */ template auto allocate_elemmat(R val) const { auto ret = array_type::tensor::from_shape(this->shape_elemmat()); ret.fill(val); return ret; } /** * Get an allocated `array_type::tensor` to store a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * Default: rank = 0, a.k.a. scalar. * Note: the container is not (zero-)initialised. * * @tparam R value-type of the array, e.g. `double`. * @returns [#nelem, #nip]. */ template auto allocate_qtensor() const { return array_type::tensor::from_shape(this->shape_qtensor()); } /** * Get an allocated and initialised `array_type::tensor` to store a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * Default: rank = 0, a.k.a. scalar. * * @tparam R value-type of the array, e.g. `double`. * @param val The value to which to initialise all items. * @returns `array_type::tensor` container of the correct shape (and rank). */ template auto allocate_qtensor(R val) const { auto ret = array_type::tensor::from_shape(this->shape_qtensor()); ret.fill(val); return ret; } /** * Get an allocated `xt::xarray` to store a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * Note: the container is not (zero-)initialised. * * @tparam R value-type of the array, e.g. `double`. * @param rank The tensor rank. * @returns `xt::xarray` container of the correct shape. */ template auto allocate_qtensor(size_t rank) const { return xt::xarray::from_shape(this->shape_qtensor(rank)); } /** * Get an allocated and initialised `xt::xarray` to store a "qtensor" of a certain rank * (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). * * @tparam R value-type of the array, e.g. `double`. * @param rank The tensor rank. * @param val The value to which to initialise all items. * @returns `array_type::tensor` container of the correct shape (and rank). */ template auto allocate_qtensor(size_t rank, R val) const { auto ret = xt::xarray::from_shape(this->shape_qtensor(rank)); ret.fill(val); return ret; } /** * Get an allocated `array_type::tensor` to store a "qscalar" (a "qtensor" of rank 0). * Note: the container is not (zero-)initialised. * * @tparam R value-type of the array, e.g. `double`. * @returns `xt::xarray` container of the correct shape. */ template auto allocate_qscalar() const { return this->allocate_qtensor<0, R>(); } /** * Get an allocated and initialised `xt::xarray` to store a "qscalar" (a "qtensor" of rank 0). * * @tparam R value-type of the array, e.g. `double`. * @param val The value to which to initialise all items. * @returns `array_type::tensor` container of the correct shape (and rank). */ template auto allocate_qscalar(R val) const { return this->allocate_qtensor<0, R>(val); } private: auto derived_cast() -> derived_type& { return *static_cast(this); } auto derived_cast() const -> const derived_type& { return *static_cast(this); } }; /** * CRTP base class for interpolation and quadrature for a generic element in Cartesian coordinates. * * Naming convention: - * - ``elemmat``: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] - * - ``elemvec``: nodal vectors stored per element, [#nelem, #nne, #ndim] - * - ``qtensor``: integration point tensor, [#nelem, #nip, #tdim, #tdim] - * - ``qscalar``: integration point scalar, [#nelem, #nip] + * - `elemmat`: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] + * - `elemvec`: nodal vectors stored per element, [#nelem, #nne, #ndim] + * - `qtensor`: integration point tensor, [#nelem, #nip, #tdim, #tdim] + * - `qscalar`: integration point scalar, [#nelem, #nip] */ template class QuadratureBaseCartesian : public QuadratureBase { public: /** * Underlying type. */ using derived_type = D; /** * Update the nodal positions. * This recomputes: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * Under the small deformations assumption this function should not be called. * - * @param x nodal coordinates (``elemvec``). Shape should match the earlier definition. + * @param x nodal coordinates (`elemvec`). Shape should match the earlier definition. */ template void update_x(const T& x) { GOOSEFEM_ASSERT(xt::has_shape(x, derived_cast().m_x.shape())); xt::noalias(derived_cast().m_x) = x; derived_cast().compute_dN_impl(); } /** * Shape function gradients (in global coordinates). - * @return ``gradN`` stored per element, per integration point [#nelem, #nip, #nne, #ndim]. + * @return `gradN` stored per element, per integration point [#nelem, #nip, #nne, #ndim]. */ auto GradN() const -> const array_type::tensor& { return derived_cast().m_dNx; } /** * Integration volume. * @return volume stored per element, per integration point [#nelem, #nip]. */ auto dV() const -> const array_type::tensor& { return derived_cast().m_vol; } /** * Interpolate element vector and evaluate at each quadrature point. * * \f$ \vec{u}(\vec{x}_q) = N_i^e(\vec{x}) \vec{u}_i^e \f$ * * @param elemvec nodal vector stored per element [#nelem, #nne, #ndim]. * @return qvector [#nelem, #nip, #ndim]. */ template auto InterpQuad_vector(const T& elemvec) const -> array_type::tensor { size_t n = elemvec.shape(2); auto qvector = array_type::tensor::from_shape(this->shape_qvector(n)); derived_cast().interpQuad_vector_impl(elemvec, qvector); return qvector; } /** * Same as InterpQuad_vector(), but writing to preallocated return. * * @param elemvec nodal vector stored per element [#nelem, #nne, #ndim]. * @param qvector [#nelem, #nip, #ndim]. */ template void interpQuad_vector(const T& elemvec, R& qvector) const { derived_cast().interpQuad_vector_impl(elemvec, qvector); } /** * Element-by-element: dyadic product of the shape function gradients and a nodal vector. * Typical input: nodal displacements. Typical output: quadrature point strains. * Within one element:: * * for e in range(nelem): * for q in range(nip): * for m in range(nne): * qtensor(e, q, i, j) += dNdx(e, q, m, i) * elemvec(e, m, j) * * Note that the functions and their gradients are precomputed upon construction, * or updated when calling update_x(). * * @param elemvec [#nelem, #nne, #ndim] * @return qtensor [#nelem, #nip, #tdim, #tdim] */ template auto GradN_vector(const T& elemvec) const -> array_type::tensor { auto qtensor = array_type::tensor::from_shape(this->template shape_qtensor<2>()); derived_cast().gradN_vector_impl(elemvec, qtensor); return qtensor; } /** * Same as GradN_vector(), but writing to preallocated return. * * @param elemvec [#nelem, #nne, #ndim] * @param qtensor overwritten [#nelem, #nip, #tdim, #tdim] */ template void gradN_vector(const T& elemvec, R& qtensor) const { derived_cast().gradN_vector_impl(elemvec, qtensor); } /** * The transposed output of GradN_vector(). * Within one element:: * * for e in range(nelem): * for q in range(nip): * for m in range(nne): * qtensor(e, q, j, i) += dNdx(e, q, m, i) * elemvec(e, m, j) * * @param elemvec [#nelem, #nne, #ndim] * @return qtensor [#nelem, #nip, #tdim, #tdim] */ template auto GradN_vector_T(const T& elemvec) const -> array_type::tensor { auto qtensor = array_type::tensor::from_shape(this->template shape_qtensor<2>()); derived_cast().gradN_vector_T_impl(elemvec, qtensor); return qtensor; } /** * Same as GradN_vector_T(), but writing to preallocated return. * * @param elemvec [#nelem, #nne, #ndim] * @param qtensor overwritten [#nelem, #nip, #tdim, #tdim] */ template void gradN_vector_T(const T& elemvec, R& qtensor) const { derived_cast().gradN_vector_T_impl(elemvec, qtensor); } /** * The symmetric output of GradN_vector(). * Without one element:: * * for e in range(nelem): * for q in range(nip): * for m in range(nne): * qtensor(e, q, i, j) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j) * qtensor(e, q, j, i) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j) * * @param elemvec [#nelem, #nne, #ndim] * @return qtensor [#nelem, #nip, #tdim, #tdim] */ template auto SymGradN_vector(const T& elemvec) const -> array_type::tensor { auto qtensor = array_type::tensor::from_shape(this->template shape_qtensor<2>()); derived_cast().symGradN_vector_impl(elemvec, qtensor); return qtensor; } /** * Same as SymGradN_vector(), but writing to preallocated return. * * @param elemvec [#nelem, #nne, #ndim] * @param qtensor overwritten [#nelem, #nip, #tdim, #tdim] */ template void symGradN_vector(const T& elemvec, R& qtensor) const { derived_cast().symGradN_vector_impl(elemvec, qtensor); } /** * Element-by-element: integral of a continuous vector-field. * * \f$ \vec{f}_i^e = \int N_i^e(\vec{x}) \vec{f}(\vec{x}) d\Omega_e \f$ * * which is integration numerically as follows * * \f$ \vec{f}_i^e = \sum\limits_q N_i^e(\vec{x}_q) \vec{f}(\vec{x}_q) \f$ * * @param qvector [#nelem, #nip. #ndim] * @return elemvec [#nelem, #nne. #ndim] */ template auto Int_N_vector_dV(const T& qvector) const -> array_type::tensor { size_t n = qvector.shape(2); auto elemvec = array_type::tensor::from_shape(this->shape_elemvec(n)); derived_cast().int_N_vector_dV_impl(qvector, elemvec); return elemvec; } /** * Same as Int_N_vector_dV(), but writing to preallocated return. * * @param qvector [#nelem, #nip. #ndim] * @param elemvec overwritten [#nelem, #nne. #ndim] */ template void int_N_vector_dV(const T& qvector, R& elemvec) const { derived_cast().int_N_vector_dV_impl(qvector, elemvec); } /** * Element-by-element: integral of the scalar product of the shape function with a scalar. * Within one one element:: * * for e in range(nelem): * for q in range(nip): * for m in range(nne): * for n in range(nne): * elemmat(e, m * ndim + i, n * ndim + i) += * N(e, q, m) * qscalar(e, q) * N(e, q, n) * dV(e, q) * - * with ``i`` a tensor dimension. + * with `i` a tensor dimension. * Note that the functions and their gradients are precomputed upon construction, * or updated when calling update_x(). * * @param qscalar [#nelem, #nip] * @return elemmat [#nelem, #nne * #ndim, #nne * #ndim] */ template auto Int_N_scalar_NT_dV(const T& qscalar) const -> array_type::tensor { auto elemmat = array_type::tensor::from_shape(this->shape_elemmat()); derived_cast().int_N_scalar_NT_dV_impl(qscalar, elemmat); return elemmat; } /** * Same as Int_N_scalar_NT_dV(), but writing to preallocated return. * * @param qscalar [#nelem, #nip] * @param elemmat overwritten [#nelem, #nne * #ndim, #nne * #ndim] */ template void int_N_scalar_NT_dV(const T& qscalar, R& elemmat) const { derived_cast().int_N_scalar_NT_dV_impl(qscalar, elemmat); } /** * Element-by-element: integral of the dot product of the shape function gradients with * a second order tensor. Typical input: stress. Typical output: nodal force. * Within one one element:: * * for e in range(nelem): * for q in range(nip): * for m in range(nne): * elemvec(e, m, j) += dNdx(e, q, m, i) * qtensor(e, q, i, j) * dV(e, q) * - * with ``i`` and ``j`` tensor dimensions. + * with `i` and `j` tensor dimensions. * Note that the functions and their gradients are precomputed upon construction, * or updated when calling update_x(). * * @param qtensor [#nelem, #nip, #ndim, #ndim] * @return elemvec [#nelem, #nne. #ndim] */ template auto Int_gradN_dot_tensor2_dV(const T& qtensor) const -> array_type::tensor { auto elemvec = array_type::tensor::from_shape(this->shape_elemvec()); derived_cast().int_gradN_dot_tensor2_dV_impl(qtensor, elemvec); return elemvec; } /** * Same as Int_gradN_dot_tensor2_dV(), but writing to preallocated return. * * @param qtensor [#nelem, #nip, #ndim, #ndim] * @param elemvec overwritten [#nelem, #nne. #ndim] */ template void int_gradN_dot_tensor2_dV(const T& qtensor, R& elemvec) const { derived_cast().int_gradN_dot_tensor2_dV_impl(qtensor, elemvec); } // Integral of the dot product // elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV /** * Element-by-element: integral of the dot products of the shape function gradients with * a fourth order tensor. Typical input: stiffness tensor. Typical output: stiffness matrix. * Within one one element:: * * for e in range(nelem): * for q in range(nip): * for m in range(nne): * for n in range(nne): * elemmat(e, m * ndim + j, n * ndim + k) += * dNdx(e,q,m,i) * qtensor(e,q,i,j,k,l) * dNdx(e,q,n,l) * dV(e,q) * - * with ``i``, ``j``, ``k``, and ``l`` tensor dimensions. + * with `i`, `j`, `k`, and `l` tensor dimensions. * Note that the functions and their gradients are precomputed upon construction, * or updated when calling update_x(). * * @param qtensor [#nelem, #nip, #ndim, #ndim, #ndim, #ndim] * @return elemmat [#nelem, #nne * #ndim, #nne * #ndim] */ template auto Int_gradN_dot_tensor4_dot_gradNT_dV(const T& qtensor) const -> array_type::tensor { auto elemmat = array_type::tensor::from_shape(this->shape_elemmat()); derived_cast().int_gradN_dot_tensor4_dot_gradNT_dV_impl(qtensor, elemmat); return elemmat; } /** * Same as Int_gradN_dot_tensor4_dot_gradNT_dV(), but writing to preallocated return. * * @param qtensor [#nelem, #nip, #ndim, #ndim, #ndim, #ndim] * @param elemmat overwritten [#nelem, #nne * #ndim, #nne * #ndim] */ template void int_gradN_dot_tensor4_dot_gradNT_dV(const T& qtensor, R& elemmat) const { derived_cast().int_gradN_dot_tensor4_dot_gradNT_dV_impl(qtensor, elemmat); } protected: /** * Update the shape function gradients (called when the nodal positions are updated). */ void compute_dN() { derived_cast().compute_dN_impl(); } private: auto derived_cast() -> derived_type& { return *static_cast(this); } auto derived_cast() const -> const derived_type& { return *static_cast(this); } friend class QuadratureBase; template void interpQuad_vector_impl(const T& elemvec, R& qvector) const { size_t n = elemvec.shape(2); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec(n))); GOOSEFEM_ASSERT(xt::has_shape(qvector, this->shape_qvector(n))); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& N = derived_cast().m_N; qvector.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto fq = &elemvec(e, 0, 0); for (size_t q = 0; q < nip; ++q) { auto Nq = &N(q, 0); auto tq = &qvector(e, q, 0); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < n; ++i) { tq[i] += Nq[m] * fq[m * n + i]; } } } } } template void gradN_vector_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->template shape_qtensor<2>())); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& dNx = derived_cast().m_dNx; qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto ue = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape()); auto graduq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { graduq(i, j) += dNxq(m, i) * ue(m, j); } } } } } } template void gradN_vector_T_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->template shape_qtensor<2>())); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& dNx = derived_cast().m_dNx; qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto ue = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape()); auto graduq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { graduq(j, i) += dNxq(m, i) * ue(m, j); } } } } } } template void symGradN_vector_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->template shape_qtensor<2>())); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& dNx = derived_cast().m_dNx; qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto ue = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape()); auto epsq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { epsq(i, j) += 0.5 * dNxq(m, i) * ue(m, j); epsq(j, i) += 0.5 * dNxq(m, i) * ue(m, j); } } } } } } template void int_N_vector_dV_impl(const T& qvector, R& elemvec) const { size_t n = qvector.shape(2); GOOSEFEM_ASSERT(xt::has_shape(qvector, this->shape_qvector(n))); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec(n))); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& N = derived_cast().m_N; auto& vol = derived_cast().m_vol; elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto f = &elemvec(e, 0, 0); for (size_t q = 0; q < nip; ++q) { auto Ne = &N(q, 0); auto tq = &qvector(e, q, 0); auto& volq = vol(e, q); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < n; ++i) { f[m * n + i] += Ne[m] * tq[i] * volq; } } } } } template void int_N_scalar_NT_dV_impl(const T& qscalar, R& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, this->shape_qscalar())); GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat())); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& N = derived_cast().m_N; auto& vol = derived_cast().m_vol; elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto Me = xt::adapt( &elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto Ne = xt::adapt(&N(q, 0), xt::xshape()); auto& volq = vol(e, q); auto& rho = qscalar(e, q); // M(m * D::s_ndim + i, n * D::s_ndim + i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < D::s_nne; ++m) { for (size_t n = 0; n < D::s_nne; ++n) { for (size_t i = 0; i < D::s_ndim; ++i) { Me(m * D::s_ndim + i, n * D::s_ndim + i) += Ne(m) * rho * Ne(n) * volq; } } } } } } template void int_gradN_dot_tensor2_dV_impl(const T& qtensor, R& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->template shape_qtensor<2>())); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& dNx = derived_cast().m_m_dNx; auto& vol = derived_cast().m_vol; elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto fe = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape()); auto sigq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& volq = vol(e, q); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { fe(m, j) += dNxq(m, i) * sigq(i, j) * volq; } } } } } } template void int_gradN_dot_tensor4_dot_gradNT_dV_impl(const T& qtensor, R& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->template shape_qtensor<4>())); GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat())); auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& dNx = derived_cast().m_dNx; auto& vol = derived_cast().m_vol; elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < nelem; ++e) { auto K = xt::adapt( &elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape()); auto Cq = xt::adapt( &qtensor(e, q, 0, 0, 0, 0), xt::xshape()); auto& volq = vol(e, q); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t n = 0; n < D::s_nne; ++n) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { for (size_t k = 0; k < D::s_ndim; ++k) { for (size_t l = 0; l < D::s_ndim; ++l) { K(m * D::s_ndim + j, n * D::s_ndim + k) += dNxq(m, i) * Cq(i, j, k, l) * dNxq(n, l) * volq; } } } } } } } } } void compute_dN_impl() { auto nelem = derived_cast().m_nelem; auto nip = derived_cast().m_nip; auto& vol = derived_cast().m_vol; auto& w = derived_cast().m_w; auto& dNxi = derived_cast().m_dNxi; auto& dNx = derived_cast().m_dNx; auto& x = derived_cast().m_x; dNx.fill(0.0); #pragma omp parallel { auto J = array_type::tensor::from_shape({D::s_ndim, D::s_ndim}); auto Jinv = array_type::tensor::from_shape({D::s_ndim, D::s_ndim}); #pragma omp for for (size_t e = 0; e < nelem; ++e) { auto xe = xt::adapt(&x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < nip; ++q) { auto dNxiq = xt::adapt(&dNxi(q, 0, 0), xt::xshape()); auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape()); J.fill(0.0); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { J(i, j) += dNxiq(m, i) * xe(m, j); } } } double Jdet = detail::tensor::inv(J, Jinv); for (size_t m = 0; m < D::s_nne; ++m) { for (size_t i = 0; i < D::s_ndim; ++i) { for (size_t j = 0; j < D::s_ndim; ++j) { dNxq(m, i) += Jinv(i, j) * dNxiq(m, i); } } } vol(e, q) = w(q) * Jdet; } } } } }; } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/ElementHex8.h b/include/GooseFEM/ElementHex8.h index e570101..a50662c 100644 --- a/include/GooseFEM/ElementHex8.h +++ b/include/GooseFEM/ElementHex8.h @@ -1,404 +1,404 @@ /** * Quadrature for 8-noded hexahedral element in 3d (GooseFEM::Mesh::ElementType::Hex8), * in a Cartesian coordinate system. * * @file ElementHex8.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENTHEX8_H #define GOOSEFEM_ELEMENTHEX8_H #include "config.h" namespace GooseFEM { namespace Element { /** * 8-noded hexahedral element in 3d (GooseFEM::Mesh::ElementType::Hex8). */ namespace Hex8 { /** * gauss quadrature: quadrature points such that integration is exact for these bi-linear elements:: */ namespace Gauss { /** * Number of integration points: * * nip = nne = 8 * * @return unsigned int */ inline size_t nip() { return 8; } /** * Integration point coordinates (local coordinates). * * @return Coordinates [#nip, ndim], with `ndim = 3`. */ inline array_type::tensor xi() { size_t nip = 8; size_t ndim = 3; array_type::tensor xi = xt::empty({nip, ndim}); xi(0, 0) = -1.0 / std::sqrt(3.0); xi(0, 1) = -1.0 / std::sqrt(3.0); xi(0, 2) = -1.0 / std::sqrt(3.0); xi(1, 0) = +1.0 / std::sqrt(3.0); xi(1, 1) = -1.0 / std::sqrt(3.0); xi(1, 2) = -1.0 / std::sqrt(3.0); xi(2, 0) = +1.0 / std::sqrt(3.0); xi(2, 1) = +1.0 / std::sqrt(3.0); xi(2, 2) = -1.0 / std::sqrt(3.0); xi(3, 0) = -1.0 / std::sqrt(3.0); xi(3, 1) = +1.0 / std::sqrt(3.0); xi(3, 2) = -1.0 / std::sqrt(3.0); xi(4, 0) = -1.0 / std::sqrt(3.0); xi(4, 1) = -1.0 / std::sqrt(3.0); xi(4, 2) = +1.0 / std::sqrt(3.0); xi(5, 0) = +1.0 / std::sqrt(3.0); xi(5, 1) = -1.0 / std::sqrt(3.0); xi(5, 2) = +1.0 / std::sqrt(3.0); xi(6, 0) = +1.0 / std::sqrt(3.0); xi(6, 1) = +1.0 / std::sqrt(3.0); xi(6, 2) = +1.0 / std::sqrt(3.0); xi(7, 0) = -1.0 / std::sqrt(3.0); xi(7, 1) = +1.0 / std::sqrt(3.0); xi(7, 2) = +1.0 / std::sqrt(3.0); return xi; } /** * Integration point weights. * * @return Coordinates [#nip]. */ inline array_type::tensor w() { size_t nip = 8; array_type::tensor w = xt::empty({nip}); w(0) = 1.0; w(1) = 1.0; w(2) = 1.0; w(3) = 1.0; w(4) = 1.0; w(5) = 1.0; w(6) = 1.0; w(7) = 1.0; return w; } } // namespace Gauss /** * nodal quadrature: quadrature points coincide with the nodes. * The order is the same as in the connectivity. */ namespace Nodal { /** * Number of integration points: * * nip = nne = 8 * * @return unsigned int */ inline size_t nip() { return 8; } /** * Integration point coordinates (local coordinates). * - * @return Coordinates [#nip, `ndim`], with ``ndim = 3``. + * @return Coordinates [#nip, `ndim`], with `ndim = 3`. */ inline array_type::tensor xi() { size_t nip = 8; size_t ndim = 3; array_type::tensor xi = xt::empty({nip, ndim}); xi(0, 0) = -1.0; xi(0, 1) = -1.0; xi(0, 2) = -1.0; xi(1, 0) = +1.0; xi(1, 1) = -1.0; xi(1, 2) = -1.0; xi(2, 0) = +1.0; xi(2, 1) = +1.0; xi(2, 2) = -1.0; xi(3, 0) = -1.0; xi(3, 1) = +1.0; xi(3, 2) = -1.0; xi(4, 0) = -1.0; xi(4, 1) = -1.0; xi(4, 2) = +1.0; xi(5, 0) = +1.0; xi(5, 1) = -1.0; xi(5, 2) = +1.0; xi(6, 0) = +1.0; xi(6, 1) = +1.0; xi(6, 2) = +1.0; xi(7, 0) = -1.0; xi(7, 1) = +1.0; xi(7, 2) = +1.0; return xi; } /** * Integration point weights. * * @return Coordinates [#nip]. */ inline array_type::tensor w() { size_t nip = 8; array_type::tensor w = xt::empty({nip}); w(0) = 1.0; w(1) = 1.0; w(2) = 1.0; w(3) = 1.0; w(4) = 1.0; w(5) = 1.0; w(6) = 1.0; w(7) = 1.0; return w; } } // namespace Nodal /** * Interpolation and quadrature. * * Fixed dimensions: * - `ndim = 3`: number of dimensions. * - `nne = 8`: number of nodes per element. * * Naming convention: * - `elemmat`: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] * - `elemvec`: nodal vectors stored per element, [#nelem, #nne, #ndim] * - `qtensor`: integration point tensor, [#nelem, #nip, #ndim, #ndim] * - `qscalar`: integration point scalar, [#nelem, #nip] */ class Quadrature : public QuadratureBaseCartesian { public: Quadrature() = default; /** * Constructor: use default Gauss integration. * The following is pre-computed during construction: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * They can be reused without any cost. * They only have to be recomputed when the nodal position changes * (note that they are assumed to be constant under a small-strain assumption). * In that case use update_x() to update the nodal positions and * to recompute the above listed quantities. * - * @param x nodal coordinates (``elemvec``). + * @param x nodal coordinates (`elemvec`). */ template Quadrature(const T& x) : Quadrature(x, Gauss::xi(), Gauss::w()) { } /** * Constructor with custom integration. * The following is pre-computed during construction: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * They can be reused without any cost. * They only have to be recomputed when the nodal position changes * (note that they are assumed to be constant under a small-strain assumption). * In that case use update_x() to update the nodal positions and * to recompute the above listed quantities. * - * @param x nodal coordinates (``elemvec``). + * @param x nodal coordinates (`elemvec`). * @param xi Integration point coordinates (local coordinates) [#nip]. * @param w Integration point weights [#nip]. */ template Quadrature(const T& x, const X& xi, const W& w) { m_x = x; m_w = w; m_xi = xi; m_nip = w.size(); m_nelem = m_x.shape(0); m_N = xt::empty({m_nip, s_nne}); m_dNxi = xt::empty({m_nip, s_nne, s_ndim}); for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2)); m_N(q, 1) = 0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2)); m_N(q, 2) = 0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2)); m_N(q, 3) = 0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2)); m_N(q, 4) = 0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2)); m_N(q, 5) = 0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2)); m_N(q, 6) = 0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2)); m_N(q, 7) = 0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.125 * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2)); m_dNxi(q, 1, 0) = +0.125 * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2)); m_dNxi(q, 2, 0) = +0.125 * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2)); m_dNxi(q, 3, 0) = -0.125 * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2)); m_dNxi(q, 4, 0) = -0.125 * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2)); m_dNxi(q, 5, 0) = +0.125 * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2)); m_dNxi(q, 6, 0) = +0.125 * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2)); m_dNxi(q, 7, 0) = -0.125 * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 2)); m_dNxi(q, 1, 1) = -0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 2)); m_dNxi(q, 2, 1) = +0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 2)); m_dNxi(q, 3, 1) = +0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 2)); m_dNxi(q, 4, 1) = -0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 2)); m_dNxi(q, 5, 1) = -0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 2)); m_dNxi(q, 6, 1) = +0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 2)); m_dNxi(q, 7, 1) = +0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 2)); // - dN / dxi_2 m_dNxi(q, 0, 2) = -0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)); m_dNxi(q, 1, 2) = -0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)); m_dNxi(q, 2, 2) = -0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)); m_dNxi(q, 3, 2) = -0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)); m_dNxi(q, 4, 2) = +0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)); m_dNxi(q, 5, 2) = +0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)); m_dNxi(q, 6, 2) = +0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)); m_dNxi(q, 7, 2) = +0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)); } GOOSEFEM_ASSERT(m_x.shape(1) == s_nne); GOOSEFEM_ASSERT(m_x.shape(2) == s_ndim); GOOSEFEM_ASSERT(xt::has_shape(m_xi, {m_nip, s_ndim})); GOOSEFEM_ASSERT(xt::has_shape(m_w, {m_nip})); GOOSEFEM_ASSERT(xt::has_shape(m_N, {m_nip, s_nne})); GOOSEFEM_ASSERT(xt::has_shape(m_dNxi, {m_nip, s_nne, s_ndim})); m_dNx = xt::empty({m_nelem, m_nip, s_nne, s_ndim}); m_vol = xt::empty(this->shape_qscalar()); this->compute_dN(); } private: friend QuadratureBase; friend QuadratureBaseCartesian; template void int_N_scalar_NT_dV_impl(const T& qscalar, R& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, this->shape_qscalar())); GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat())); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m * ndim + i, n * ndim + i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < s_nne; ++m) { for (size_t n = 0; n < s_nne; ++n) { M(m * s_ndim + 0, n * s_ndim + 0) += N(m) * rho * N(n) * vol; M(m * s_ndim + 1, n * s_ndim + 1) += N(m) * rho * N(n) * vol; M(m * s_ndim + 2, n * s_ndim + 2) += N(m) * rho * N(n) * vol; } } } } } template void int_gradN_dot_tensor2_dV_impl(const T& qtensor, R& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& v = m_vol(e, q); for (size_t m = 0; m < s_nne; ++m) { f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0) + dNx(m, 2) * sig(2, 0)) * v; f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1) + dNx(m, 2) * sig(2, 1)) * v; f(m, 2) += (dNx(m, 0) * sig(0, 2) + dNx(m, 1) * sig(1, 2) + dNx(m, 2) * sig(2, 2)) * v; } } } } constexpr static size_t s_nne = 8; ///< Number of nodes per element. constexpr static size_t s_ndim = 3; ///< Number of dimensions for nodal vectors. constexpr static size_t s_tdim = 3; ///< Number of dimensions for tensors. size_t m_tdim = 3; ///< Dynamic alias of s_tdim (remove in C++17) size_t m_nelem; ///< Number of elements. size_t m_nip; ///< Number of integration points per element. array_type::tensor m_x; ///< nodal positions stored per element [#nelem, #nne, #ndim] array_type::tensor m_w; ///< weight of each integration point [nip] array_type::tensor m_xi; ///< local coordinate per integration point [#nip, #ndim] array_type::tensor m_N; ///< shape functions [#nip, #nne] array_type::tensor m_dNxi; ///< local shape func grad [#nip, #nne, #ndim] array_type::tensor m_dNx; ///< global shape func grad [#nelem, #nip, #nne, #ndim] array_type::tensor m_vol; ///< integration point volume [#nelem, #nip] }; } // namespace Hex8 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/ElementQuad4Axisymmetric.h b/include/GooseFEM/ElementQuad4Axisymmetric.h index f9cd9a7..4880005 100644 --- a/include/GooseFEM/ElementQuad4Axisymmetric.h +++ b/include/GooseFEM/ElementQuad4Axisymmetric.h @@ -1,465 +1,465 @@ /** * Quadrature for 4-noded quadrilateral element in 2d (GooseFEM::Mesh::ElementType::Quad4), * in an axisymmetric coordinated system. * * @file ElementQuad4Axisymmetric.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENTQUAD4AXISYMMETRIC_H #define GOOSEFEM_ELEMENTQUAD4AXISYMMETRIC_H #include "config.h" namespace GooseFEM { namespace Element { namespace Quad4 { /** * Interpolation and quadrature. * * Fixed dimensions: - * - ``ndim = 2``: number of dimensions. - * - ``tdim = 3``: number of dimensions or tensor. - * - ``nne = 4``: number of nodes per element. + * - `ndim = 2`: number of dimensions. + * - `tdim = 3`: number of dimensions or tensor. + * - `nne = 4`: number of nodes per element. * * Naming convention: - * - ``elemmat``: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] - * - ``elemvec``: nodal vectors stored per element, [#nelem, #nne, #ndim] - * - ``qtensor``: integration point tensor, [#nelem, #nip, #tdim, #tdim] - * - ``qscalar``: integration point scalar, [#nelem, #nip] + * - `elemmat`: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] + * - `elemvec`: nodal vectors stored per element, [#nelem, #nne, #ndim] + * - `qtensor`: integration point tensor, [#nelem, #nip, #tdim, #tdim] + * - `qscalar`: integration point scalar, [#nelem, #nip] */ class QuadratureAxisymmetric : public QuadratureBaseCartesian { public: QuadratureAxisymmetric() = default; /** * Constructor: use default Gauss integration. * The following is pre-computed during construction: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * They can be reused without any cost. * They only have to be recomputed when the nodal position changes * (note that they are assumed to be constant under a small-strain assumption). * In that case use update_x() to update the nodal positions and * to recompute the above listed quantities. * - * @param x nodal coordinates (``elemvec``). + * @param x nodal coordinates (`elemvec`). */ template QuadratureAxisymmetric(const T& x) : QuadratureAxisymmetric(x, Gauss::xi(), Gauss::w()) { } /** * Constructor with custom integration. * The following is pre-computed during construction: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * They can be reused without any cost. * They only have to be recomputed when the nodal position changes * (note that they are assumed to be constant under a small-strain assumption). * In that case use update_x() to update the nodal positions and * to recompute the above listed quantities. * - * @param x nodal coordinates (``elemvec``). + * @param x nodal coordinates (`elemvec`). * @param xi Integration point coordinates (local coordinates) [#nip]. * @param w Integration point weights [#nip]. */ template QuadratureAxisymmetric(const T& x, const X& xi, const W& w) { m_x = x; m_w = w; m_xi = xi; m_nip = w.size(); m_nelem = m_x.shape(0); GOOSEFEM_ASSERT(m_x.shape(1) == s_nne); GOOSEFEM_ASSERT(m_x.shape(2) == s_ndim); GOOSEFEM_ASSERT(xt::has_shape(m_xi, {m_nip, s_ndim})); GOOSEFEM_ASSERT(xt::has_shape(m_w, {m_nip})); m_N = xt::empty({m_nip, s_nne}); m_dNxi = xt::empty({m_nip, s_nne, s_ndim}); m_B = xt::empty({m_nelem, m_nip, s_nne, s_tdim, s_tdim, s_tdim}); m_vol = xt::empty(this->shape_qscalar()); for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 1) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 2) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_N(q, 3) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 1, 0) = +0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 2, 0) = +0.25 * (1.0 + m_xi(q, 1)); m_dNxi(q, 3, 0) = -0.25 * (1.0 + m_xi(q, 1)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.25 * (1.0 - m_xi(q, 0)); m_dNxi(q, 1, 1) = -0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 2, 1) = +0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 3, 1) = +0.25 * (1.0 - m_xi(q, 0)); } this->compute_dN_impl(); } /** * Get the B-matrix (shape function gradients) (in global coordinates). * Note that the functions and their gradients are precomputed upon construction, * or updated when calling update_x(). * - * @return ``B`` matrix stored per element, per integration point [#nelem, #nne, #tdim, #tdim, + * @return `B` matrix stored per element, per integration point [#nelem, #nne, #tdim, #tdim, * #tdim] */ const array_type::tensor& B() const { return m_B; } private: friend QuadratureBase; friend QuadratureBaseCartesian; // qtensor(e, q, i, j) += B(e, q, m, i, j, k) * elemvec(e, q, m, k) template void gradN_vector_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += B(m,i,j,k) * u(m,perm(k)) // (where perm(0) = 1, perm(2) = 0) gradu(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); gradu(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); gradu(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); gradu(0, 2) = B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0); gradu(2, 0) = B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1); } } } template void gradN_vector_T_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += B(m,i,j,k) * u(m,perm(k)) // (where perm(0) = 1, perm(2) = 0) gradu(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); gradu(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); gradu(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); gradu(2, 0) = B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0); gradu(0, 2) = B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1); } } } template void symGradN_vector_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += B(m,i,j,k) * u(m,perm(k)) // eps(j,i) = 0.5 * (gradu(i,j) + gradu(j,i)) // (where perm(0) = 1, perm(2) = 0) eps(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); eps(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); eps(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); eps(2, 0) = 0.5 * (B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0) + B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1)); eps(0, 2) = eps(2, 0); } } } // elemmat(e, q, m * ndim + i, n * ndim + i) += // N(e, q, m) * qscalar(e, q) * N(e, q, n) * dV(e, q) template void int_N_scalar_NT_dV_impl(const T& qscalar, R& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, this->shape_qscalar())); GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat())); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < s_nne; ++m) { for (size_t n = 0; n < s_nne; ++n) { M(m * s_ndim + 0, n * s_ndim + 0) += N(m) * rho * N(n) * vol; M(m * s_ndim + 1, n * s_ndim + 1) += N(m) * rho * N(n) * vol; } } } } } // fm = ( Bm^T : qtensor ) dV template void int_gradN_dot_tensor2_dV_impl(const T& qtensor, R& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); // f(m,i) += B(m,i,j,perm(k)) * sig(i,j) * dV // (where perm(0) = 1, perm(2) = 0) for (size_t m = 0; m < s_nne; ++m) { f(m, 0) += vol * (B(m, 2, 2, 2) * sig(2, 2) + B(m, 0, 2, 2) * sig(0, 2)); f(m, 1) += vol * (B(m, 0, 0, 0) * sig(0, 0) + B(m, 1, 1, 0) * sig(1, 1) + B(m, 2, 0, 0) * sig(2, 0)); } } } } // Kmn = ( Bm^T : qtensor : Bn ) dV template void int_gradN_dot_tensor4_dot_gradNT_dV_impl(const T& qtensor, R& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<4>())); GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat())); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto C = xt::adapt( &qtensor(e, q, 0, 0, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); // K(m*s_ndim+perm(c), n*s_ndim+perm(f)) = B(m,a,b,c) * C(a,b,d,e) * B(n,e,d,f) * // vol; (where perm(0) = 1, perm(2) = 0) for (size_t m = 0; m < s_nne; ++m) { for (size_t n = 0; n < s_nne; ++n) { K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 0, 0) * B(n, 0, 0, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 1, 1) * B(n, 1, 1, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 0) += B(m, 0, 0, 0) * C(0, 0, 2, 2) * B(n, 2, 2, 2) * vol; K(m * s_ndim + 1, n * s_ndim + 0) += B(m, 0, 0, 0) * C(0, 0, 2, 0) * B(n, 0, 2, 2) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 0, 2) * B(n, 2, 0, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 0, 0) * B(n, 0, 0, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 1, 1) * B(n, 1, 1, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 0) += B(m, 1, 1, 0) * C(1, 1, 2, 2) * B(n, 2, 2, 2) * vol; K(m * s_ndim + 1, n * s_ndim + 0) += B(m, 1, 1, 0) * C(1, 1, 2, 0) * B(n, 0, 2, 2) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 0, 2) * B(n, 2, 0, 0) * vol; K(m * s_ndim + 0, n * s_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 0, 0) * B(n, 0, 0, 0) * vol; K(m * s_ndim + 0, n * s_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 1, 1) * B(n, 1, 1, 0) * vol; K(m * s_ndim + 0, n * s_ndim + 0) += B(m, 2, 2, 2) * C(2, 2, 2, 2) * B(n, 2, 2, 2) * vol; K(m * s_ndim + 0, n * s_ndim + 0) += B(m, 2, 2, 2) * C(2, 2, 2, 0) * B(n, 0, 2, 2) * vol; K(m * s_ndim + 0, n * s_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 0, 2) * B(n, 2, 0, 0) * vol; K(m * s_ndim + 0, n * s_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 0, 0) * B(n, 0, 0, 0) * vol; K(m * s_ndim + 0, n * s_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 1, 1) * B(n, 1, 1, 0) * vol; K(m * s_ndim + 0, n * s_ndim + 0) += B(m, 0, 2, 2) * C(0, 2, 2, 2) * B(n, 2, 2, 2) * vol; K(m * s_ndim + 0, n * s_ndim + 0) += B(m, 0, 2, 2) * C(0, 2, 2, 0) * B(n, 0, 2, 2) * vol; K(m * s_ndim + 0, n * s_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 0, 2) * B(n, 2, 0, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 0, 0) * B(n, 0, 0, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 1, 1) * B(n, 1, 1, 0) * vol; K(m * s_ndim + 1, n * s_ndim + 0) += B(m, 2, 0, 0) * C(2, 0, 2, 2) * B(n, 2, 2, 2) * vol; K(m * s_ndim + 1, n * s_ndim + 0) += B(m, 2, 0, 0) * C(2, 0, 2, 0) * B(n, 0, 2, 2) * vol; K(m * s_ndim + 1, n * s_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 0, 2) * B(n, 2, 0, 0) * vol; } } } } } void compute_dN_impl() { // most components remain zero, and are not written m_B.fill(0.0); #pragma omp parallel { array_type::tensor J = xt::empty({2, 2}); array_type::tensor Jinv = xt::empty({2, 2}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto B = xt::adapt( &m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto N = xt::adapt(&m_N(q, 0), xt::xshape()); // J(i,j) += dNxi(m,i) * x(m,j); J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) + dNxi(3, 0) * x(3, 0); J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) + dNxi(3, 0) * x(3, 1); J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) + dNxi(3, 1) * x(3, 0); J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) + dNxi(3, 1) * x(3, 1); double Jdet = detail::tensor<2>::inv(J, Jinv); // radius for computation of volume double rq = N(0) * x(0, 1) + N(1) * x(1, 1) + N(2) * x(2, 1) + N(3) * x(3, 1); // dNx(m,i) += Jinv(i,j) * dNxi(m,j) for (size_t m = 0; m < s_nne; ++m) { // B(m, r, r, r) = dNdx(m,1) B(m, 0, 0, 0) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); // B(m, r, z, z) = dNdx(m,1) B(m, 0, 2, 2) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); // B(m, t, t, r) B(m, 1, 1, 0) = 1.0 / rq * N(m); // B(m, z, r, r) = dNdx(m,0) B(m, 2, 0, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); // B(m, z, z, z) = dNdx(m,0) B(m, 2, 2, 2) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); } m_vol(e, q) = m_w(q) * Jdet * 2.0 * M_PI * rq; } } } } constexpr static size_t s_nne = 4; ///< Number of nodes per element. constexpr static size_t s_ndim = 2; ///< Number of dimensions for nodal vectors. constexpr static size_t s_tdim = 3; ///< Number of dimensions for tensors. size_t m_tdim = 3; ///< Dynamic alias of s_tdim (remove in C++17) size_t m_nelem; ///< Number of elements. size_t m_nip; ///< Number of integration points per element. array_type::tensor m_x; ///< nodal positions stored per element [#nelem, #nne, #ndim] array_type::tensor m_w; ///< weight per integration point [nip] array_type::tensor m_xi; ///< local coordinate per integration point [#nip, #ndim] array_type::tensor m_N; ///< shape functions [#nip, #nne] array_type::tensor m_dNxi; ///< local shape func grad [#nip, #nne, #ndim] array_type::tensor m_vol; ///< integration point volume [#nelem, #nip] array_type::tensor m_B; ///< B-matrix [#nelem, #nne, #tdim, #tdim, #tdim] }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/ElementQuad4Planar.h b/include/GooseFEM/ElementQuad4Planar.h index f1a4d07..cf278ee 100644 --- a/include/GooseFEM/ElementQuad4Planar.h +++ b/include/GooseFEM/ElementQuad4Planar.h @@ -1,341 +1,341 @@ /** * Quadrature for 4-noded quadrilateral element in 2d (GooseFEM::Mesh::ElementType::Quad4), * in a Cartesian coordinate system. * The different with ElementQuad4.h is that here the tensors live in 3d and are assumed plane * strain. * * @file ElementQuad4Planar.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENTQUAD4PLANAR_H #define GOOSEFEM_ELEMENTQUAD4PLANAR_H #include "config.h" #include "detail.h" namespace GooseFEM { namespace Element { namespace Quad4 { /** * Interpolation and quadrature. * Similar to Element::Quad4::Quadrature with the only different that quadrature point tensors * are 3d ("plane strain") while the mesh is 2d. * * Fixed dimensions: - * - ``ndim = 2``: number of dimensions. - * - ``tdim = 3``: number of dimensions or tensor. - * - ``nne = 4``: number of nodes per element. + * - `ndim = 2`: number of dimensions. + * - `tdim = 3`: number of dimensions or tensor. + * - `nne = 4`: number of nodes per element. * * Naming convention: - * - ``elemmat``: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] - * - ``elemvec``: nodal vectors stored per element, [#nelem, #nne, #ndim] - * - ``qtensor``: integration point tensor, [#nelem, #nip, #tdim, #tdim] - * - ``qscalar``: integration point scalar, [#nelem, #nip] + * - `elemmat`: matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim] + * - `elemvec`: nodal vectors stored per element, [#nelem, #nne, #ndim] + * - `qtensor`: integration point tensor, [#nelem, #nip, #tdim, #tdim] + * - `qscalar`: integration point scalar, [#nelem, #nip] */ class QuadraturePlanar : public QuadratureBaseCartesian { public: QuadraturePlanar() = default; /** * Constructor: use default Gauss integration. * The following is pre-computed during construction: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * They can be reused without any cost. * They only have to be recomputed when the nodal position changes * (note that they are assumed to be constant under a small-strain assumption). * In that case use update_x() to update the nodal positions and * to recompute the above listed quantities. * - * @param x nodal coordinates (``elemvec``). + * @param x nodal coordinates (`elemvec`). * @param thick out-of-plane thickness (incorporated in the element volume). */ template QuadraturePlanar(const T& x, double thick = 1.0) : QuadraturePlanar(x, Gauss::xi(), Gauss::w(), thick) { } /** * Constructor with custom integration. * The following is pre-computed during construction: * - the shape functions, * - the shape function gradients (in local and global) coordinates, * - the integration points volumes. * They can be reused without any cost. * They only have to be recomputed when the nodal position changes * (note that they are assumed to be constant under a small-strain assumption). * In that case use update_x() to update the nodal positions and * to recompute the above listed quantities. * - * @param x nodal coordinates (``elemvec``). + * @param x nodal coordinates (`elemvec`). * @param xi Integration point coordinates (local coordinates) [#nip]. * @param w Integration point weights [#nip]. * @param thick out-of-plane thickness (incorporated in the element volume). */ template QuadraturePlanar(const T& x, const X& xi, const W& w, double thick = 1.0) { m_x = x; m_w = w; m_xi = xi; m_nip = w.size(); m_nelem = m_x.shape(0); m_N = xt::empty({m_nip, s_nne}); m_dNxi = xt::empty({m_nip, s_nne, s_ndim}); m_thick = thick; for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.25 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)); m_N(q, 1) = 0.25 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)); m_N(q, 2) = 0.25 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)); m_N(q, 3) = 0.25 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.25 * (1.0 - xi(q, 1)); m_dNxi(q, 1, 0) = +0.25 * (1.0 - xi(q, 1)); m_dNxi(q, 2, 0) = +0.25 * (1.0 + xi(q, 1)); m_dNxi(q, 3, 0) = -0.25 * (1.0 + xi(q, 1)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.25 * (1.0 - xi(q, 0)); m_dNxi(q, 1, 1) = -0.25 * (1.0 + xi(q, 0)); m_dNxi(q, 2, 1) = +0.25 * (1.0 + xi(q, 0)); m_dNxi(q, 3, 1) = +0.25 * (1.0 - xi(q, 0)); } GOOSEFEM_ASSERT(m_x.shape(1) == s_nne); GOOSEFEM_ASSERT(m_x.shape(2) == s_ndim); GOOSEFEM_ASSERT(xt::has_shape(m_xi, {m_nip, s_ndim})); GOOSEFEM_ASSERT(xt::has_shape(m_w, {m_nip})); GOOSEFEM_ASSERT(xt::has_shape(m_N, {m_nip, s_nne})); GOOSEFEM_ASSERT(xt::has_shape(m_dNxi, {m_nip, s_nne, s_ndim})); m_dNx = xt::empty({m_nelem, m_nip, s_nne, s_ndim}); m_vol = xt::empty(this->shape_qscalar()); this->compute_dN_impl(); } private: friend QuadratureBase; friend QuadratureBaseCartesian; template void gradN_vector_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += dNx(m,i) * u(m,j) gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); gradu(0, 1) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1); gradu(1, 0) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0); gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); } } } template void gradN_vector_T_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += dNx(m,i) * u(m,j) gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); gradu(1, 0) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1); gradu(0, 1) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0); gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); } } } template void symGradN_vector_impl(const T& elemvec, R& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += dNx(m,i) * u(m,j) // eps(j,i) = 0.5 * (gradu(i,j) + gradu(j,i)) eps(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); eps(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); eps(0, 1) = 0.5 * (dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1) + dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0)); eps(1, 0) = eps(0, 1); } } } template void int_N_scalar_NT_dV_impl(const T& qscalar, R& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, this->shape_qscalar())); GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat())); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < s_nne; ++m) { for (size_t n = 0; n < s_nne; ++n) { M(m * s_ndim + 0, n * s_ndim + 0) += N(m) * rho * N(n) * vol; M(m * s_ndim + 1, n * s_ndim + 1) += N(m) * rho * N(n) * vol; } } } } } template void int_gradN_dot_tensor2_dV_impl(const T& qtensor, R& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>())); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); for (size_t m = 0; m < s_nne; ++m) { f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0)) * vol; f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1)) * vol; } } } } void compute_dN_impl() { #pragma omp parallel { array_type::tensor J = xt::empty({2, 2}); array_type::tensor Jinv = xt::empty({2, 2}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); // J(i,j) += dNxi(m,i) * x(m,j); J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) + dNxi(3, 0) * x(3, 0); J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) + dNxi(3, 0) * x(3, 1); J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) + dNxi(3, 1) * x(3, 0); J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) + dNxi(3, 1) * x(3, 1); double Jdet = detail::tensor<2>::inv(J, Jinv); // dNx(m,i) += Jinv(i,j) * dNxi(m,j); for (size_t m = 0; m < s_nne; ++m) { dNx(m, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); dNx(m, 1) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); } m_vol(e, q) = m_w(q) * Jdet * m_thick; } } } } constexpr static size_t s_nne = 4; ///< Number of nodes per element. constexpr static size_t s_ndim = 2; ///< Number of dimensions for nodal vectors. constexpr static size_t s_tdim = 3; ///< Dynamic alias of s_tdim (remove in C++17) size_t m_tdim = 3; ///< Number of dimensions for tensors. size_t m_nelem; ///< Number of elements. size_t m_nip; ///< Number of integration points per element. array_type::tensor m_x; ///< nodal positions stored per element [#nelem, #nne, #ndim] array_type::tensor m_w; ///< weight of each integration point [nip] array_type::tensor m_xi; ///< local coordinate per integration point [#nip, #ndim] array_type::tensor m_N; ///< shape functions [#nip, #nne] array_type::tensor m_dNxi; ///< local shape func grad [#nip, #nne, #ndim] array_type::tensor m_dNx; ///< global shape func grad [#nelem, #nip, #nne, #ndim] array_type::tensor m_vol; ///< integration point volume [#nelem, #nip] double m_thick; ///< out-of-plane thickness }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h index 69931f4..a7b330d 100644 --- a/include/GooseFEM/Mesh.h +++ b/include/GooseFEM/Mesh.h @@ -1,2862 +1,2862 @@ /** * Generic mesh operations. * * @file Mesh.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESH_H #define GOOSEFEM_MESH_H #include "ElementQuad4.h" #include "MatrixDiagonal.h" #include "Vector.h" #include "assertions.h" #include "config.h" namespace GooseFEM { /** * Generic mesh operations, and simple mesh definitions. */ namespace Mesh { template inline std::vector> nodaltyings(const D& dofs); /** * Enumerator for element-types */ enum class ElementType { Unknown, ///< Unknown element-type Quad4, ///< Quadrilateral: 4-noded element in 2-d Hex8, ///< Hexahedron: 8-noded element in 3-d Tri3 ///< Triangle: 3-noded element in 2-d }; /** * Extract the element type based on the connectivity. * * @param coor Nodal coordinates [nnode, ndim]. * @param conn Connectivity [nelem, nne]. * @return ElementType(). */ template 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"); } 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 /** * List with DOF-numbers in sequential order. * The output is a sequential list of DOF-numbers for each vector-component of each node. * For example for 3 nodes in 2 dimensions the output is * * \f$ \begin{bmatrix} 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \f$ * * @param nnode Number of nodes. * @param ndim Number of dimensions. * @return DOF-numbers. */ inline array_type::tensor dofs(size_t nnode, size_t ndim) { return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } /** * Renumber indices to lowest possible index. For example: * * \f$ \begin{bmatrix} 0 & 1 \\ 5 & 4 \end{bmatrix} \f$ * * is renumbered to * * \f$ \begin{bmatrix} 0 & 1 \\ 3 & 2 \end{bmatrix} \f$ * * Or, in pseudo-code, the result of this function is that: * * dofs = renumber(dofs) * sort(unique(dofs[:])) == range(max(dofs+1)) * * \note One can use the wrapper function renumber(). This class gives more advanced features. */ class Renumber { public: Renumber() = default; /** * @param dofs DOF-numbers. */ template Renumber(const T& dofs) { size_t n = xt::amax(dofs)() + 1; size_t i = 0; array_type::tensor unique = xt::unique(dofs); m_renum = xt::empty({n}); for (auto& j : unique) { m_renum(j) = i; ++i; } } /** * Apply renumbering to other set. * * @param list List of (DOF-)numbers. * @return Renumbered list of (DOF-)numbers. */ template T apply(const T& list) const { return detail::renum(list, m_renum); } /** * Get the list needed to renumber, e.g.: * * dofs_renumbered(i, j) = index(dofs(i, j)) * * @return Renumber-index. */ const array_type::tensor& index() const { return m_renum; } private: array_type::tensor m_renum; }; /** * Renumber to lowest possible index (see GooseFEM::Mesh::Renumber). * * @param dofs DOF-numbers [nnode, ndim]. * @return Renumbered DOF-numbers. */ template inline T renumber(const T& dofs) { return Renumber(dofs).apply(dofs); } /** * CRTP base class for regular meshes. */ template class RegularBase { public: /** * Underlying type. */ using derived_type = D; /** * Number of elements. * @return unsigned int */ auto nelem() const { return derived_cast().m_nelem; } /** * Number of nodes. * @return unsigned int */ auto nnode() const { return derived_cast().m_nnode; } /** * Number of nodes-per-element == 4. * @return unsigned int */ auto nne() const { return derived_cast().m_nne; } /** * Number of dimensions == 2. * @return unsigned int */ auto ndim() const { return derived_cast().m_ndim; } /** * Number of elements in x-direction == width of the mesh in units of #h. * @return unsigned int */ auto nelx() const { return derived_cast().nelx_impl(); } /** * Number of elements in y-direction == height of the mesh, in units of #h, * @return unsigned int */ auto nely() const { return derived_cast().nely_impl(); } /** * Linear edge size of one 'block'. * @return double */ auto h() const { return derived_cast().m_h; } /** * The ElementType(). * @return element type */ auto getElementType() const { return derived_cast().getElementType_impl(); } /** * Nodal coordinates [#nnode, #ndim]. * @return coordinates per node */ auto coor() const { return derived_cast().coor_impl(); } /** * Connectivity [#nelem, #nne]. * @return nodes per element */ auto conn() const { return derived_cast().conn_impl(); } /** * DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. * @return DOFs per node */ auto dofs() const { return GooseFEM::Mesh::dofs(this->nnode(), this->ndim()); } /** * DOF-numbers for the case that the periodicity if fully eliminated. * @return DOF numbers for each node [#nnode, #ndim]. */ auto dofsPeriodic() const { array_type::tensor ret = this->dofs(); array_type::tensor nodePer = this->nodesPeriodic(); array_type::tensor independent = xt::view(nodePer, xt::all(), 0); array_type::tensor 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); } /** * Periodic node pairs, in two columns: (independent, dependent). * @return [ntyings, #ndim]. */ auto nodesPeriodic() const { return derived_cast().nodesPeriodic_impl(); } /** * Reference node to use for periodicity, because all corners are tied to it. * @return Node number. */ auto nodesOrigin() const { return derived_cast().nodesOrigin_impl(); } private: auto derived_cast() -> derived_type& { return *static_cast(this); } auto derived_cast() const -> const derived_type& { return *static_cast(this); } }; /** * CRTP base class for regular meshes in 2d. */ template class RegularBase2d : public RegularBase { public: /** * Underlying type. */ using derived_type = D; /** * Nodes along the bottom edge (y = 0), in order of increasing x. * @return List of node numbers. */ auto nodesBottomEdge() const { return derived_cast().nodesBottomEdge_impl(); } /** * Nodes along the top edge (y = #nely * #h), in order of increasing x. * @return List of node numbers. */ auto nodesTopEdge() const { return derived_cast().nodesTopEdge_impl(); } /** * Nodes along the left edge (x = 0), in order of increasing y. * @return List of node numbers. */ auto nodesLeftEdge() const { return derived_cast().nodesLeftEdge_impl(); } /** * Nodes along the right edge (x = #nelx * #h), in order of increasing y. * @return List of node numbers. */ auto nodesRightEdge() const { return derived_cast().nodesRightEdge_impl(); } /** * Nodes along the bottom edge (y = 0), without the corners (at x = 0 and x = #nelx * #h). * Same as: nodesBottomEdge()[1: -1]. * @return List of node numbers. */ auto nodesBottomOpenEdge() const { return derived_cast().nodesBottomOpenEdge_impl(); } /** * Nodes along the top edge (y = #nely * #h), without the corners (at x = 0 and x = #nelx * #h). * Same as: nodesTopEdge()[1: -1]. * @return List of node numbers. */ auto nodesTopOpenEdge() const { return derived_cast().nodesTopOpenEdge_impl(); } /** * Nodes along the left edge (x = 0), without the corners (at y = 0 and y = #nely * #h). * Same as: nodesLeftEdge()[1: -1]. * @return List of node numbers. */ auto nodesLeftOpenEdge() const { return derived_cast().nodesLeftOpenEdge_impl(); } /** * Nodes along the right edge (x = #nelx * #h), without the corners (at y = 0 and y = #nely * * #h). Same as: nodesRightEdge()[1: -1]. * @return List of node numbers. */ auto nodesRightOpenEdge() const { return derived_cast().nodesRightOpenEdge_impl(); } /** * The bottom-left corner node (at x = 0, y = 0). * Same as nodesBottomEdge()[0] and nodesLeftEdge()[0]. * @return Node number. */ auto nodesBottomLeftCorner() const { return derived_cast().nodesBottomLeftCorner_impl(); } /** * The bottom-right corner node (at x = #nelx * #h, y = 0). * Same as nodesBottomEdge()[-1] and nodesRightEdge()[0]. * @return Node number. */ auto nodesBottomRightCorner() const { return derived_cast().nodesBottomRightCorner_impl(); } /** * The top-left corner node (at x = 0, y = #nely * #h). * Same as nodesTopEdge()[0] and nodesRightEdge()[-1]. * @return Node number. */ auto nodesTopLeftCorner() const { return derived_cast().nodesTopLeftCorner_impl(); } /** * The top-right corner node (at x = #nelx * #h, y = #nely * #h). * Same as nodesTopEdge()[-1] and nodesRightEdge()[-1]. * @return Node number. */ auto nodesTopRightCorner() const { return derived_cast().nodesTopRightCorner_impl(); } /** * Alias of nodesBottomLeftCorner(). * @return Node number. */ auto nodesLeftBottomCorner() const { return derived_cast().nodesBottomLeftCorner_impl(); } /** * Alias of nodesTopLeftCorner(). * @return Node number. */ auto nodesLeftTopCorner() const { return derived_cast().nodesTopLeftCorner_impl(); } /** * Alias of nodesBottomRightCorner(). * @return Node number. */ auto nodesRightBottomCorner() const { return derived_cast().nodesBottomRightCorner_impl(); } /** * Alias of nodesTopRightCorner(). * @return Node number. */ auto nodesRightTopCorner() const { return derived_cast().nodesTopRightCorner_impl(); } private: auto derived_cast() -> derived_type& { return *static_cast(this); } auto derived_cast() const -> const derived_type& { return *static_cast(this); } friend class RegularBase; array_type::tensor nodesPeriodic_impl() const { array_type::tensor bot = derived_cast().nodesBottomOpenEdge_impl(); array_type::tensor top = derived_cast().nodesTopOpenEdge_impl(); array_type::tensor lft = derived_cast().nodesLeftOpenEdge_impl(); array_type::tensor rgt = derived_cast().nodesRightOpenEdge_impl(); std::array shape = {bot.size() + lft.size() + size_t(3), size_t(2)}; auto ret = array_type::tensor::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; } auto nodesOrigin_impl() const { return derived_cast().nodesBottomLeftCorner_impl(); } }; /** * CRTP base class for regular meshes in 3d. */ template class RegularBase3d : public RegularBase { public: /** * Underlying type. */ using derived_type = D; /** * Number of elements in y-direction == height of the mesh, in units of #h, * @return unsigned int */ auto nelz() const { return derived_cast().nelz_impl(); } /** * Nodes along the bottom face (y = 0). * @return List of node numbers. */ auto nodesBottom() const { return derived_cast().nodesBottom_impl(); } /** * Nodes along the top face (y = #nely * #h). * @return List of node numbers. */ auto nodesTop() const { return derived_cast().nodesTop_impl(); } /** * Nodes along the left face (x = 0). * @return List of node numbers. */ auto nodesLeft() const { return derived_cast().nodesLeft_impl(); } /** * Nodes along the right face (x = #nelx * #h). * @return List of node numbers. */ auto nodesRight() const { return derived_cast().nodesRight_impl(); } /** * Nodes along the front face (z = 0). * @return List of node numbers. */ auto nodesFront() const { return derived_cast().nodesFront_impl(); } /** * Nodes along the back face (z = #nelz * #h). * @return List of node numbers. */ auto nodesBack() const { return derived_cast().nodesBack_impl(); } /** * Nodes along the edge at the intersection of the front and bottom faces * (z = 0 and y = 0). * @return List of node numbers. */ auto nodesFrontBottomEdge() const { return derived_cast().nodesFrontBottomEdge_impl(); } /** * Nodes along the edge at the intersection of the front and top faces * (z = 0 and y = #nely * #h). * @return List of node numbers. */ auto nodesFrontTopEdge() const { return derived_cast().nodesFrontTopEdge_impl(); } /** * Nodes along the edge at the intersection of the front and left faces * (z = 0 and x = 0). * @return List of node numbers. */ auto nodesFrontLeftEdge() const { return derived_cast().nodesFrontLeftEdge_impl(); } /** * Nodes along the edge at the intersection of the front and right faces * (z = 0 and x = #nelx * #h). * @return List of node numbers. */ auto nodesFrontRightEdge() const { return derived_cast().nodesFrontRightEdge_impl(); } /** * Nodes along the edge at the intersection of the back and bottom faces * (z = #nelz * #h and y = #nely * #h). * @return List of node numbers. */ auto nodesBackBottomEdge() const { return derived_cast().nodesBackBottomEdge_impl(); } /** * Nodes along the edge at the intersection of the back and top faces * (z = #nelz * #h and x = 0). * @return List of node numbers. */ auto nodesBackTopEdge() const { return derived_cast().nodesBackTopEdge_impl(); } /** * Nodes along the edge at the intersection of the back and left faces * (z = #nelz * #h and x = #nelx * #h). * @return List of node numbers. */ auto nodesBackLeftEdge() const { return derived_cast().nodesBackLeftEdge_impl(); } /** * Nodes along the edge at the intersection of the back and right faces * (? = #nelz * #h and ? = ?). * @return List of node numbers. */ auto nodesBackRightEdge() const { return derived_cast().nodesBackRightEdge_impl(); } /** * Nodes along the edge at the intersection of the bottom and left faces * (y = 0 and x = 0). * @return List of node numbers. */ auto nodesBottomLeftEdge() const { return derived_cast().nodesBottomLeftEdge_impl(); } /** * Nodes along the edge at the intersection of the bottom and right faces * (y = 0 and x = #nelx * #h). * @return List of node numbers. */ auto nodesBottomRightEdge() const { return derived_cast().nodesBottomRightEdge_impl(); } /** * Nodes along the edge at the intersection of the top and left faces * (y = 0 and x = #nelx * #h). * @return List of node numbers. */ auto nodesTopLeftEdge() const { return derived_cast().nodesTopLeftEdge_impl(); } /** * Nodes along the edge at the intersection of the top and right faces * (y = #nely * #h and x = #nelx * #h). * @return List of node numbers. */ auto nodesTopRightEdge() const { return derived_cast().nodesTopRightEdge_impl(); } /** * Alias of nodesFrontBottomEdge() * @return List of node numbers. */ auto nodesBottomFrontEdge() const { return derived_cast().nodesFrontBottomEdge_impl(); } /** * Alias of nodesBackBottomEdge() * @return List of node numbers. */ auto nodesBottomBackEdge() const { return derived_cast().nodesBackBottomEdge_impl(); } /** * Alias of nodesFrontTopEdge() * @return List of node numbers. */ auto nodesTopFrontEdge() const { return derived_cast().nodesFrontTopEdge_impl(); } /** * Alias of nodesBackTopEdge() * @return List of node numbers. */ auto nodesTopBackEdge() const { return derived_cast().nodesBackTopEdge_impl(); } /** * Alias of nodesBottomLeftEdge() * @return List of node numbers. */ auto nodesLeftBottomEdge() const { return derived_cast().nodesBottomLeftEdge_impl(); } /** * Alias of nodesFrontLeftEdge() * @return List of node numbers. */ auto nodesLeftFrontEdge() const { return derived_cast().nodesFrontLeftEdge_impl(); } /** * Alias of nodesBackLeftEdge() * @return List of node numbers. */ auto nodesLeftBackEdge() const { return derived_cast().nodesBackLeftEdge_impl(); } /** * Alias of nodesTopLeftEdge() * @return List of node numbers. */ auto nodesLeftTopEdge() const { return derived_cast().nodesTopLeftEdge_impl(); } /** * Alias of nodesBottomRightEdge() * @return List of node numbers. */ auto nodesRightBottomEdge() const { return derived_cast().nodesBottomRightEdge_impl(); } /** * Alias of nodesTopRightEdge() * @return List of node numbers. */ auto nodesRightTopEdge() const { return derived_cast().nodesTopRightEdge_impl(); } /** * Alias of nodesFrontRightEdge() * @return List of node numbers. */ auto nodesRightFrontEdge() const { return derived_cast().nodesFrontRightEdge_impl(); } /** * Alias of nodesBackRightEdge() * @return List of node numbers. */ auto nodesRightBackEdge() const { return derived_cast().nodesBackRightEdge_impl(); } /** * Nodes along the front face excluding edges. * Same as different between nodesFront() and * [nodesFrontBottomEdge(), nodesFrontTopEdge(), nodesFrontLeftEdge(), nodesFrontRightEdge()] * @return list of node numbers. */ auto nodesFrontFace() const { return derived_cast().nodesFrontFace_impl(); } /** * Nodes along the back face excluding edges. * Same as different between nodesBack() and * [nodesBackBottomEdge(), nodesBackTopEdge(), nodesBackLeftEdge(), nodesBackRightEdge()] * @return list of node numbers. */ auto nodesBackFace() const { return derived_cast().nodesBackFace_impl(); } /** * Nodes along the left face excluding edges. * Same as different between nodesLeft() and * [nodesFrontLeftEdge(), nodesBackLeftEdge(), nodesBottomLeftEdge(), nodesTopLeftEdge()] * @return list of node numbers. */ auto nodesLeftFace() const { return derived_cast().nodesLeftFace_impl(); } /** * Nodes along the right face excluding edges. * Same as different between nodesRight() and * [nodesFrontRightEdge(), nodesBackRightEdge(), nodesBottomRightEdge(), nodesTopRightEdge()] * @return list of node numbers. */ auto nodesRightFace() const { return derived_cast().nodesRightFace_impl(); } /** * Nodes along the bottom face excluding edges. * Same as different between nodesBottom() and * [nodesBackBottomEdge(), nodesBackTopEdge(), nodesBackLeftEdge(), nodesBackRightEdge()] * @return list of node numbers. */ auto nodesBottomFace() const { return derived_cast().nodesBottomFace_impl(); } /** * Nodes along the top face excluding edges. * Same as different between nodesTop() and * [nodesFrontBottomEdge(), nodesFrontTopEdge(), nodesFrontLeftEdge(), nodesFrontRightEdge()] * @return list of node numbers. */ auto nodesTopFace() const { return derived_cast().nodesTopFace_impl(); } /** * Same as nodesFrontBottomEdge() but without corners. * @return List of node numbers. */ auto nodesFrontBottomOpenEdge() const { return derived_cast().nodesFrontBottomOpenEdge_impl(); } /** * Same as nodesFrontTopEdge() but without corners. * @return List of node numbers. */ auto nodesFrontTopOpenEdge() const { return derived_cast().nodesFrontTopOpenEdge_impl(); } /** * Same as nodesFrontLeftEdge() but without corners. * @return List of node numbers. */ auto nodesFrontLeftOpenEdge() const { return derived_cast().nodesFrontLeftOpenEdge_impl(); } /** * Same as nodesFrontRightEdge() but without corners. * @return List of node numbers. */ auto nodesFrontRightOpenEdge() const { return derived_cast().nodesFrontRightOpenEdge_impl(); } /** * Same as nodesBackBottomEdge() but without corners. * @return List of node numbers. */ auto nodesBackBottomOpenEdge() const { return derived_cast().nodesBackBottomOpenEdge_impl(); } /** * Same as nodesBackTopEdge() but without corners. * @return List of node numbers. */ auto nodesBackTopOpenEdge() const { return derived_cast().nodesBackTopOpenEdge_impl(); } /** * Same as nodesBackLeftEdge() but without corners. * @return List of node numbers. */ auto nodesBackLeftOpenEdge() const { return derived_cast().nodesBackLeftOpenEdge_impl(); } /** * Same as nodesBackRightEdge() but without corners. * @return List of node numbers. */ auto nodesBackRightOpenEdge() const { return derived_cast().nodesBackRightOpenEdge_impl(); } /** * Same as nodesBottomLeftEdge() but without corners. * @return List of node numbers. */ auto nodesBottomLeftOpenEdge() const { return derived_cast().nodesBottomLeftOpenEdge_impl(); } /** * Same as nodesBottomRightEdge() but without corners. * @return List of node numbers. */ auto nodesBottomRightOpenEdge() const { return derived_cast().nodesBottomRightOpenEdge_impl(); } /** * Same as nodesTopLeftEdge() but without corners. * @return List of node numbers. */ auto nodesTopLeftOpenEdge() const { return derived_cast().nodesTopLeftOpenEdge_impl(); } /** * Same as nodesTopRightEdge() but without corners. * @return List of node numbers. */ auto nodesTopRightOpenEdge() const { return derived_cast().nodesTopRightOpenEdge_impl(); } /** * Alias of nodesFrontBottomOpenEdge(). * @return List of node numbers. */ auto nodesBottomFrontOpenEdge() const { return derived_cast().nodesFrontBottomOpenEdge_impl(); } /** * Alias of nodesBackBottomOpenEdge(). * @return List of node numbers. */ auto nodesBottomBackOpenEdge() const { return derived_cast().nodesBackBottomOpenEdge_impl(); } /** * Alias of nodesFrontTopOpenEdge(). * @return List of node numbers. */ auto nodesTopFrontOpenEdge() const { return derived_cast().nodesFrontTopOpenEdge_impl(); } /** * Alias of nodesBackTopOpenEdge(). * @return List of node numbers. */ auto nodesTopBackOpenEdge() const { return derived_cast().nodesBackTopOpenEdge_impl(); } /** * Alias of nodesBottomLeftOpenEdge(). * @return List of node numbers. */ auto nodesLeftBottomOpenEdge() const { return derived_cast().nodesBottomLeftOpenEdge_impl(); } /** * Alias of nodesFrontLeftOpenEdge(). * @return List of node numbers. */ auto nodesLeftFrontOpenEdge() const { return derived_cast().nodesFrontLeftOpenEdge_impl(); } /** * Alias of nodesBackLeftOpenEdge(). * @return List of node numbers. */ auto nodesLeftBackOpenEdge() const { return derived_cast().nodesBackLeftOpenEdge_impl(); } /** * Alias of nodesTopLeftOpenEdge(). * @return List of node numbers. */ auto nodesLeftTopOpenEdge() const { return derived_cast().nodesTopLeftOpenEdge_impl(); } /** * Alias of nodesBottomRightOpenEdge(). * @return List of node numbers. */ auto nodesRightBottomOpenEdge() const { return derived_cast().nodesBottomRightOpenEdge_impl(); } /** * Alias of nodesTopRightOpenEdge(). * @return List of node numbers. */ auto nodesRightTopOpenEdge() const { return derived_cast().nodesTopRightOpenEdge_impl(); } /** * Alias of nodesFrontRightOpenEdge(). * @return List of node numbers. */ auto nodesRightFrontOpenEdge() const { return derived_cast().nodesFrontRightOpenEdge_impl(); } /** * Alias of nodesBackRightOpenEdge(). * @return List of node numbers. */ auto nodesRightBackOpenEdge() const { return derived_cast().nodesBackRightOpenEdge_impl(); } /** * Front-Bottom-Left corner node. * @return Node number. */ auto nodesFrontBottomLeftCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } /** * Front-Bottom-Right corner node. * @return Node number. */ auto nodesFrontBottomRightCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } /** * Front-Top-Left corner node. * @return Node number. */ auto nodesFrontTopLeftCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } /** * Front-Top-Right corner node. * @return Node number. */ auto nodesFrontTopRightCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } /** * Back-Bottom-Left corner node. * @return Node number. */ auto nodesBackBottomLeftCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } /** * Back-Bottom-Right corner node. * @return Node number. */ auto nodesBackBottomRightCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } /** * Back-Top-Left corner node. * @return Node number. */ auto nodesBackTopLeftCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } /** * Back-Top-Right corner node. * @return Node number. */ auto nodesBackTopRightCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } /** * Alias of nodesFrontBottomLeftCorner(). * @return Node number. */ auto nodesFrontLeftBottomCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } /** * Alias of nodesFrontBottomLeftCorner(). * @return Node number. */ auto nodesBottomFrontLeftCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } /** * Alias of nodesFrontBottomLeftCorner(). * @return Node number. */ auto nodesBottomLeftFrontCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } /** * Alias of nodesFrontBottomLeftCorner(). * @return Node number. */ auto nodesLeftFrontBottomCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } /** * Alias of nodesFrontBottomLeftCorner(). * @return Node number. */ auto nodesLeftBottomFrontCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } /** * Alias of nodesFrontBottomRightCorner(). * @return Node number. */ auto nodesFrontRightBottomCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } /** * Alias of nodesFrontBottomRightCorner(). * @return Node number. */ auto nodesBottomFrontRightCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } /** * Alias of nodesFrontBottomRightCorner(). * @return Node number. */ auto nodesBottomRightFrontCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } /** * Alias of nodesFrontBottomRightCorner(). * @return Node number. */ auto nodesRightFrontBottomCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } /** * Alias of nodesFrontBottomRightCorner(). * @return Node number. */ auto nodesRightBottomFrontCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } /** * Alias of nodesFrontTopLeftCorner(). * @return Node number. */ auto nodesFrontLeftTopCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } /** * Alias of nodesFrontTopLeftCorner(). * @return Node number. */ auto nodesTopFrontLeftCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } /** * Alias of nodesFrontTopLeftCorner(). * @return Node number. */ auto nodesTopLeftFrontCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } /** * Alias of nodesFrontTopLeftCorner(). * @return Node number. */ auto nodesLeftFrontTopCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } /** * Alias of nodesFrontTopLeftCorner(). * @return Node number. */ auto nodesLeftTopFrontCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } /** * Alias of nodesFrontTopRightCorner(). * @return Node number. */ auto nodesFrontRightTopCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } /** * Alias of nodesFrontTopRightCorner(). * @return Node number. */ auto nodesTopFrontRightCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } /** * Alias of nodesFrontTopRightCorner(). * @return Node number. */ auto nodesTopRightFrontCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } /** * Alias of nodesFrontTopRightCorner(). * @return Node number. */ auto nodesRightFrontTopCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } /** * Alias of nodesFrontTopRightCorner(). * @return Node number. */ auto nodesRightTopFrontCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } /** * Alias of nodesBackBottomLeftCorner(). * @return Node number. */ auto nodesBackLeftBottomCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } /** * Alias of nodesBackBottomLeftCorner(). * @return Node number. */ auto nodesBottomBackLeftCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } /** * Alias of nodesBackBottomLeftCorner(). * @return Node number. */ auto nodesBottomLeftBackCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } /** * Alias of nodesBackBottomLeftCorner(). * @return Node number. */ auto nodesLeftBackBottomCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } /** * Alias of nodesBackBottomLeftCorner(). * @return Node number. */ auto nodesLeftBottomBackCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } /** * Alias of nodesBackBottomRightCorner(). * @return Node number. */ auto nodesBackRightBottomCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } /** * Alias of nodesBackBottomRightCorner(). * @return Node number. */ auto nodesBottomBackRightCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } /** * Alias of nodesBackBottomRightCorner(). * @return Node number. */ auto nodesBottomRightBackCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } /** * Alias of nodesBackBottomRightCorner(). * @return Node number. */ auto nodesRightBackBottomCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } /** * Alias of nodesBackBottomRightCorner(). * @return Node number. */ auto nodesRightBottomBackCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } /** * Alias of nodesBackTopLeftCorner(). * @return Node number. */ auto nodesBackLeftTopCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } /** * Alias of nodesBackTopLeftCorner(). * @return Node number. */ auto nodesTopBackLeftCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } /** * Alias of nodesBackTopLeftCorner(). * @return Node number. */ auto nodesTopLeftBackCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } /** * Alias of nodesBackTopLeftCorner(). * @return Node number. */ auto nodesLeftBackTopCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } /** * Alias of nodesBackTopLeftCorner(). * @return Node number. */ auto nodesLeftTopBackCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } /** * Alias of nodesBackTopRightCorner(). * @return Node number. */ auto nodesBackRightTopCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } /** * Alias of nodesBackTopRightCorner(). * @return Node number. */ auto nodesTopBackRightCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } /** * Alias of nodesBackTopRightCorner(). * @return Node number. */ auto nodesTopRightBackCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } /** * Alias of nodesBackTopRightCorner(). * @return Node number. */ auto nodesRightBackTopCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } /** * Alias of nodesBackTopRightCorner(). * @return Node number. */ auto nodesRightTopBackCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } private: auto derived_cast() -> derived_type& { return *static_cast(this); } auto derived_cast() const -> const derived_type& { return *static_cast(this); } friend class RegularBase; array_type::tensor nodesPeriodic_impl() const { array_type::tensor fro = derived_cast().nodesFrontFace_impl(); array_type::tensor bck = derived_cast().nodesBackFace_impl(); array_type::tensor lft = derived_cast().nodesLeftFace_impl(); array_type::tensor rgt = derived_cast().nodesRightFace_impl(); array_type::tensor bot = derived_cast().nodesBottomFace_impl(); array_type::tensor top = derived_cast().nodesTopFace_impl(); array_type::tensor froBot = derived_cast().nodesFrontBottomOpenEdge_impl(); array_type::tensor froTop = derived_cast().nodesFrontTopOpenEdge_impl(); array_type::tensor froLft = derived_cast().nodesFrontLeftOpenEdge_impl(); array_type::tensor froRgt = derived_cast().nodesFrontRightOpenEdge_impl(); array_type::tensor bckBot = derived_cast().nodesBackBottomOpenEdge_impl(); array_type::tensor bckTop = derived_cast().nodesBackTopOpenEdge_impl(); array_type::tensor bckLft = derived_cast().nodesBackLeftOpenEdge_impl(); array_type::tensor bckRgt = derived_cast().nodesBackRightOpenEdge_impl(); array_type::tensor botLft = derived_cast().nodesBottomLeftOpenEdge_impl(); array_type::tensor botRgt = derived_cast().nodesBottomRightOpenEdge_impl(); array_type::tensor topLft = derived_cast().nodesTopLeftOpenEdge_impl(); array_type::tensor 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 ret = xt::empty({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; } auto nodesOrigin_impl() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } }; /** * Find overlapping nodes. The output has the following structure: * * [ * [nodes_from_mesh_a], * [nodes_from_mesh_b], * ] * * @param coor_a Nodal coordinates of mesh "a" [nnode, ndim]. * @param coor_b Nodal coordinates of mesh "b" [nnode, ndim]. * @param rtol Relative tolerance for position match. * @param atol Absolute tolerance for position match. * @return Overlapping nodes. */ template inline array_type::tensor overlapping(const S& coor_a, const T& coor_b, double rtol = 1e-5, double atol = 1e-8) { GOOSEFEM_ASSERT(coor_a.dimension() == 2); GOOSEFEM_ASSERT(coor_b.dimension() == 2); 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); } } array_type::tensor 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; } /** * Stitch two mesh objects, specifying overlapping nodes by hand. */ class ManualStitch { public: ManualStitch() = default; /** * @param coor_a Nodal coordinates of mesh "a" [nnode, ndim]. * @param conn_a Connectivity of mesh "a" [nelem, nne]. * @param overlapping_nodes_a Node-numbers of mesh "a" that overlap with mesh "b" [n]. * @param coor_b Nodal coordinates of mesh "b" [nnode, ndim]. * @param conn_b Connectivity of mesh "b" [nelem, nne]. * @param overlapping_nodes_b Node-numbers of mesh "b" that overlap with mesh "a" [n]. - * @param check_position If ``true`` the nodes are checked for position overlap. + * @param check_position If `true` the nodes are checked for position overlap. * @param rtol Relative tolerance for check on position overlap. * @param atol Absolute tolerance for check on position overlap. */ template 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 = true, double rtol = 1e-5, double atol = 1e-8) { 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 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()); } /** * Number of sub meshes == 2. * @return unsigned int */ size_t nmesh() const { return 2; } /** * Number of elements. * @return unsigned int */ size_t nelem() const { return m_conn.shape(0); } /** * Number of nodes. * @return unsigned int */ size_t nnode() const { return m_coor.shape(0); } /** * Number of nodes-per-element. * @return unsigned int */ size_t nne() const { return m_conn.shape(1); } /** * Number of dimensions. * @return unsigned int */ size_t ndim() const { return m_coor.shape(1); } /** * Nodal coordinates [#nnode, #ndim]. * @return coordinates per node */ const array_type::tensor& coor() const { return m_coor; } /** * Connectivity [#nelem, #nne]. * @return nodes per element */ const array_type::tensor& conn() const { return m_conn; } /** * DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. * @return DOFs per node */ array_type::tensor dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } /** * Node-map per sub-mesh. * @return nodes per mesh */ std::vector> nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } /** * Element-map per sub-mesh. * @return elements per mesh */ std::vector> elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } /** * @param mesh_index Index of the mesh ("a" = 1, "b" = 1). * @return Node-map for a given mesh. */ array_type::tensor nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nnd_a); } return m_map_b; } /** * @param mesh_index Index of the mesh ("a" = 1, "b" = 1). * @return Element-map for a given mesh. */ array_type::tensor elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nel_a); } return xt::arange(m_nel_b) + m_nel_a; } /** * Convert set of node numbers for an original mesh to the stitched mesh. * * @param set List of node numbers. * @param mesh_index Index of the mesh ("a" = 1, "b" = 1). * @return List of node numbers for the stitched mesh. */ template T 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); } /** * Convert set of element numbers for an original mesh to the stitched mesh. * * @param set List of element numbers. * @param mesh_index Index of the mesh ("a" = 1, "b" = 1). * @return List of element numbers for the stitched mesh. */ template T 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; } private: array_type::tensor m_coor; array_type::tensor m_conn; array_type::tensor m_map_b; size_t m_nnd_a; size_t m_nel_a; size_t m_nel_b; }; /** * Stitch mesh objects, automatically searching for overlapping nodes. */ class Stitch { public: /** * @param rtol Relative tolerance for position match. * @param atol Absolute tolerance for position match. */ Stitch(double rtol = 1e-5, double atol = 1e-8) { m_rtol = rtol; m_atol = atol; } /** * Add mesh to be stitched. * * @param coor Nodal coordinates [nnode, ndim]. * @param conn Connectivity [nelem, nne]. */ template void 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(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]); } /** * Number of sub meshes. * @return unsigned int */ size_t nmesh() const { return m_map.size(); } /** * Number of elements. * @return unsigned int */ size_t nelem() const { return m_conn.shape(0); } /** * Number of nodes. * @return unsigned int */ size_t nnode() const { return m_coor.shape(0); } /** * Number of nodes-per-element. * @return unsigned int */ size_t nne() const { return m_conn.shape(1); } /** * Number of dimensions. * @return unsigned int */ size_t ndim() const { return m_coor.shape(1); } /** * Nodal coordinates [#nnode, #ndim]. * @return coordinates per node */ const array_type::tensor& coor() const { return m_coor; } /** * Connectivity [#nelem, #nne]. * @return nodes per element */ const array_type::tensor& conn() const { return m_conn; } /** * DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. * @return DOFs per node */ array_type::tensor dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } /** * Node-map per sub-mesh. * @return nodes per mesh */ std::vector> nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } /** * Element-map per sub-mesh. * @return elements per mesh */ std::vector> elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } /** * The node numbers in the stitched mesh that are coming from a specific sub-mesh. * * @param mesh_index Index of the sub-mesh. * @return List of node numbers. */ array_type::tensor nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return m_map[mesh_index]; } /** * The element numbers in the stitched mesh that are coming from a specific sub-mesh. * * @param mesh_index Index of the sub-mesh. * @return List of element numbers. */ array_type::tensor elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return xt::arange(m_nel[mesh_index]) + m_el_offset[mesh_index]; } /** * Convert set of node-numbers for a sub-mesh to the stitched mesh. * * @param set List of node numbers. * @param mesh_index Index of the sub-mesh. * @return List of node numbers for the stitched mesh. */ template T 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]); } /** * Convert set of element-numbers for a sub-mesh to the stitched mesh. * * @param set List of element numbers. * @param mesh_index Index of the sub-mesh. * @return List of element numbers for the stitched mesh. */ template T 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]; } /** * Combine set of node numbers for an original to the final mesh (removes duplicates). * * @param set List of node numbers per mesh. * @return List of node numbers for the stitched mesh. */ template T 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(); } array_type::tensor 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); } * / template T nodeset(std::initializer_list set) const { return this->nodeset(std::vector(set)); } /** * Combine set of element numbers for an original to the final mesh. * * @param set List of element numbers per mesh. * @return List of element numbers for the stitched mesh. */ template T 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(); } array_type::tensor 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; } * / template T elemset(std::initializer_list set) const { return this->elemset(std::vector(set)); } protected: array_type::tensor m_coor; ///< Nodal coordinates [#nnode, #ndim] array_type::tensor m_conn; ///< Connectivity [#nelem, #nne] std::vector> m_map; ///< See nodemap(size_t) std::vector m_nel; ///< Number of elements per sub-mesh. std::vector m_el_offset; ///< First element of every sub-mesh. double m_rtol; ///< Relative tolerance to find overlapping nodes. double m_atol; ///< Absolute tolerance to find overlapping nodes. }; /** * Vertically stack meshes. */ class Vstack : public Stitch { public: /** * @param check_overlap Check if nodes are overlapping when adding a mesh. * @param rtol Relative tolerance for position match. * @param atol Absolute tolerance for position match. */ Vstack(bool check_overlap = true, double rtol = 1e-5, double atol = 1e-8) { m_check_overlap = check_overlap; m_rtol = rtol; m_atol = atol; } /** * Add a mesh to the top of the current stack. * Each time the current `nodes_bot` are stitched with the then highest `nodes_top`. * * @param coor Nodal coordinates [nnode, ndim]. * @param conn Connectivity [nelem, nne]. * @param nodes_bot Nodes along the bottom edge [n]. * @param nodes_top Nodes along the top edge [n]. */ template void 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(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]); } private: std::vector> m_nodes_bot; ///< Bottom nodes of each mesh (renumbered). std::vector> m_nodes_top; ///< Top nodes of each mesh (renumbered). bool m_check_overlap; ///< Check if nodes are overlapping when adding a mesh. }; /** * Reorder to lowest possible index, in specific order. * * For example for ``Reorder({iiu, iip})`` after reordering: * * iiu = xt::range(nnu); * iip = xt::range(nnp) + nnu; */ class Reorder { public: Reorder() = default; /** * @param args List of (DOF-)numbers. */ template 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(is_unique(arg)); } #endif m_renum = xt::empty({n}); for (auto& arg : args) { for (auto& j : arg) { m_renum(j) = i; ++i; } } } /** * Apply reordering to other set. * * @param list List of (DOF-)numbers. * @return Reordered list of (DOF-)numbers. */ template T 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; } /** * Get the list needed to reorder, e.g.: * * dofs_reordered(i, j) = index(dofs(i, j)) * * @return Reorder-index. */ const array_type::tensor& index() const { return m_renum; } private: array_type::tensor m_renum; }; /** * Number of elements connected to each node. * * @param conn Connectivity [nelem, nne]. * @return Coordination per node. */ template inline array_type::tensor coordination(const E& conn) { GOOSEFEM_ASSERT(conn.dimension() == 2); size_t nnode = xt::amax(conn)() + 1; array_type::tensor N = xt::zeros({nnode}); for (auto it = conn.begin(); it != conn.end(); ++it) { N(*it) += 1; } return N; } /** * Elements connected to each node. * * @param conn Connectivity [nelem, nne]. - * @param sorted If ``true`` the list of elements for each node is sorted. + * @param sorted If `true` the list of elements for each node is sorted. * @return Elements per node [nnode, ...]. */ template inline std::vector> elem2node(const E& conn, bool sorted = true) { 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; } /** * @copydoc elem2node(const E&, bool) * * @param dofs DOFs per node, allowing accounting for periodicity [nnode, ndim]. */ template inline std::vector> elem2node(const E& conn, const D& dofs, bool sorted = true) { 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; } /** * Nodes connected to each DOF. * * @param dofs DOFs per node [nnode, ndim]. - * @param sorted If ``true`` the list of nodes for each DOF is sorted. + * @param sorted If `true` the list of nodes for each DOF is sorted. * @return Nodes per DOF [ndof, ...]. */ template inline std::vector> node2dof(const D& dofs, bool sorted = true) { return elem2node(dofs, sorted); } /** * Return size of each element edge. * * @param coor Nodal coordinates. * @param conn Connectivity. * @param type ElementType. * @return Edge-sizes per element. */ template inline array_type::tensor 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 n0 = xt::view(conn, xt::all(), 0); array_type::tensor n1 = xt::view(conn, xt::all(), 1); array_type::tensor n2 = xt::view(conn, xt::all(), 2); array_type::tensor n3 = xt::view(conn, xt::all(), 3); array_type::tensor x0 = xt::view(coor, xt::keep(n0), 0); array_type::tensor x1 = xt::view(coor, xt::keep(n1), 0); array_type::tensor x2 = xt::view(coor, xt::keep(n2), 0); array_type::tensor x3 = xt::view(coor, xt::keep(n3), 0); array_type::tensor y0 = xt::view(coor, xt::keep(n0), 1); array_type::tensor y1 = xt::view(coor, xt::keep(n1), 1); array_type::tensor y2 = xt::view(coor, xt::keep(n2), 1); array_type::tensor y3 = xt::view(coor, xt::keep(n3), 1); array_type::tensor 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"); } /** * Return size of each element edge. * The element-type is automatically determined, see defaultElementType(). * * @param coor Nodal coordinates. * @param conn Connectivity. * @return Edge-sizes per element. */ template inline array_type::tensor edgesize(const C& coor, const E& conn) { return edgesize(coor, conn, defaultElementType(coor, conn)); } /** * Coordinates of the center of each element. * * @param coor Nodal coordinates. * @param conn Connectivity. * @param type ElementType. * @return Center of each element. */ template inline array_type::tensor 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 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"); } /** * Coordinates of the center of each element. * The element-type is automatically determined, see defaultElementType(). * * @param coor Nodal coordinates. * @param conn Connectivity. * @return Center of each element. */ template inline array_type::tensor centers(const C& coor, const E& conn) { return centers(coor, conn, defaultElementType(coor, conn)); } /** * Convert an element-map to a node-map. * - * @param elem_map Element-map such that ``new_elvar = elvar[elem_map]``. + * @param elem_map Element-map such that `new_elvar = elvar[elem_map]`. * @param coor Nodal coordinates. * @param conn Connectivity. * @param type ElementType. - * @return Node-map such that ``new_nodevar = nodevar[node_map]`` + * @return Node-map such that `new_nodevar = nodevar[node_map]` */ template inline array_type::tensor 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 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) { array_type::tensor 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"); } /** * Convert an element-map to a node-map. * The element-type is automatically determined, see defaultElementType(). * - * @param elem_map Element-map such that ``new_elvar = elvar[elem_map]``. + * @param elem_map Element-map such that `new_elvar = elvar[elem_map]`. * @param coor Nodal coordinates. * @param conn Connectivity. - * @return Node-map such that ``new_nodevar = nodevar[node_map]`` + * @return Node-map such that `new_nodevar = nodevar[node_map]` */ template inline array_type::tensor elemmap2nodemap(const T& elem_map, const C& coor, const E& conn) { return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn)); } /** * Compute the mass of each node in the mesh. * If nodes are not part of the connectivity the mass is set to zero, * such that the center of gravity is simply:: * * average(coor, GooseFEM.Mesh.nodal_mass(coor, conn, type), axis=0); * * @tparam C e.g. `array_type::tensor` * @tparam E e.g. `array_type::tensor` * @param coor Nodal coordinates `[nnode, ndim]`. * @param conn Connectivity `[nelem, nne]`. * @param type ElementType. * @return Center of gravity `[ndim]`. */ template inline array_type::tensor nodal_mass(const C& coor, const E& conn, ElementType type) { auto dof = dofs(coor.shape(0), coor.shape(1)); GooseFEM::MatrixDiagonal M(conn, dof); GooseFEM::Vector vector(conn, dof); array_type::tensor rho = xt::ones(conn.shape()); if (type == ElementType::Quad4) { GooseFEM::Element::Quad4::Quadrature quad( vector.AsElement(coor), GooseFEM::Element::Quad4::Nodal::xi(), GooseFEM::Element::Quad4::Nodal::w()); M.assemble(quad.Int_N_scalar_NT_dV(rho)); } else { throw std::runtime_error("Element-type not implemented"); } return vector.AsNode(M.data()); } /** * Compute the mass of each node in the mesh. * If nodes are not part of the connectivity the mass is set to zero, * such that the center of gravity is simply:: * * average(coor, GooseFEM.Mesh.nodal_mass(coor, conn), axis=0); * * @tparam C e.g. `array_type::tensor` * @tparam E e.g. `array_type::tensor` * @param coor Nodal coordinates `[nnode, ndim]`. * @param conn Connectivity `[nelem, nne]`. * @return Center of gravity `[ndim]`. */ template inline array_type::tensor nodal_mass(const C& coor, const E& conn) { return nodal_mass(coor, conn, defaultElementType(coor, conn)); } namespace detail { // todo: remove after upstream fix template array_type::tensor average_axis_0(const T& data, const T& weights) { GOOSEFEM_ASSERT(data.dimension() == 2); GOOSEFEM_ASSERT(xt::has_shape(data, weights.shape())); array_type::tensor ret = xt::zeros({weights.shape(1)}); for (size_t j = 0; j < weights.shape(1); ++j) { double norm = 0.0; for (size_t i = 0; i < weights.shape(0); ++i) { ret(j) += data(i, j) * weights(i, j); norm += weights(i, j); } ret(j) /= norm; } return ret; } } // namespace detail /** * Compute the center of gravity of a mesh. * * @tparam C e.g. `array_type::tensor` * @tparam E e.g. `array_type::tensor` * @param coor Nodal coordinates `[nnode, ndim]`. * @param conn Connectivity `[nelem, nne]`. * @param type ElementType. * @return Center of gravity `[ndim]`. */ template inline array_type::tensor center_of_gravity(const C& coor, const E& conn, ElementType type) { // todo: remove after upstream fix return detail::average_axis_0(coor, nodal_mass(coor, conn, type)); // return xt::average(coor, nodal_mass(coor, conn, type), 0); } /** * Compute the center of gravity of a mesh. * * @tparam C e.g. `array_type::tensor` * @tparam E e.g. `array_type::tensor` * @param coor Nodal coordinates `[nnode, ndim]`. * @param conn Connectivity `[nelem, nne]`. * @return Center of gravity `[ndim]`. */ template inline array_type::tensor center_of_gravity(const C& coor, const E& conn) { // todo: remove after upstream fix return detail::average_axis_0(coor, nodal_mass(coor, conn, defaultElementType(coor, conn))); // return xt::average(coor, nodal_mass(coor, conn, defaultElementType(coor, conn)), 0); } /** * List nodal tyings based on DOF-numbers per node. * * @param dofs DOFs per node [nnode, ndim]. * @return Nodes to which the nodes is connected (sorted) [nnode, ...] */ template inline std::vector> nodaltyings(const D& dofs) { size_t nnode = dofs.shape(0); size_t ndim = dofs.shape(1); auto nodemap = node2dof(dofs); std::vector> 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; } } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h index 25397ab..c94bbf1 100644 --- a/include/GooseFEM/MeshQuad4.h +++ b/include/GooseFEM/MeshQuad4.h @@ -1,1774 +1,1774 @@ /** * Generate simple meshes of 4-noded quadrilateral elements in 2d * (GooseFEM::Mesh::ElementType::Quad4). * * @file MeshQuad4.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHQUAD4_H #define GOOSEFEM_MESHQUAD4_H #include "Mesh.h" #include "config.h" namespace GooseFEM { namespace Mesh { /** * Simple meshes of 4-noded quadrilateral elements in 2d (ElementType::Quad4). */ namespace Quad4 { // pre-allocation namespace Map { class FineLayer2Regular; } /** * Regular mesh: equi-sized elements. */ class Regular : public RegularBase2d { public: Regular() = default; /** * Constructor. * * @param nelx Number of elements in horizontal (x) direction. * @param nely Number of elements in vertical (y) direction. * @param h Edge size (width == height). */ Regular(size_t nelx, size_t nely, double h = 1.0) { m_h = h; m_nelx = nelx; m_nely = nely; m_ndim = 2; m_nne = 4; GOOSEFEM_ASSERT(m_nelx >= 1); GOOSEFEM_ASSERT(m_nely >= 1); m_nnode = (m_nelx + 1) * (m_nely + 1); m_nelem = m_nelx * m_nely; } /** * Element numbers as 'matrix'. * * @return [#nely, #nelx]. */ array_type::tensor elementgrid() const { return xt::arange(m_nelem).reshape({m_nely, m_nelx}); } private: friend class RegularBase; friend class RegularBase2d; size_t nelx_impl() const { return m_nelx; } size_t nely_impl() const { return m_nely; } ElementType getElementType_impl() const { return ElementType::Quad4; } array_type::tensor coor_impl() const { array_type::tensor ret = xt::empty({m_nnode, m_ndim}); array_type::tensor x = xt::linspace(0.0, m_h * static_cast(m_nelx), m_nelx + 1); array_type::tensor 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; } array_type::tensor conn_impl() const { array_type::tensor 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; } array_type::tensor nodesBottomEdge_impl() const { return xt::arange(m_nelx + 1); } array_type::tensor nodesTopEdge_impl() const { return xt::arange(m_nelx + 1) + m_nely * (m_nelx + 1); } array_type::tensor nodesLeftEdge_impl() const { return xt::arange(m_nely + 1) * (m_nelx + 1); } array_type::tensor nodesRightEdge_impl() const { return xt::arange(m_nely + 1) * (m_nelx + 1) + m_nelx; } array_type::tensor nodesBottomOpenEdge_impl() const { return xt::arange(1, m_nelx); } array_type::tensor nodesTopOpenEdge_impl() const { return xt::arange(1, m_nelx) + m_nely * (m_nelx + 1); } array_type::tensor nodesLeftOpenEdge_impl() const { return xt::arange(1, m_nely) * (m_nelx + 1); } array_type::tensor nodesRightOpenEdge_impl() const { return xt::arange(1, m_nely) * (m_nelx + 1) + m_nelx; } size_t nodesBottomLeftCorner_impl() const { return 0; } size_t nodesBottomRightCorner_impl() const { return m_nelx; } size_t nodesTopLeftCorner_impl() const { return m_nely * (m_nelx + 1); } size_t nodesTopRightCorner_impl() const { return m_nely * (m_nelx + 1) + m_nelx; } double m_h; ///< See h() size_t m_nelx; ///< See nelx() size_t m_nely; ///< See nely() size_t m_nelem; ///< See nelem() size_t m_nnode; ///< See nnode() size_t m_nne; ///< See nne() size_t m_ndim; ///< See ndim() }; /** * Mesh with fine middle layer, and coarser elements towards the top and bottom. */ class FineLayer : public RegularBase2d { public: FineLayer() = default; /** * Constructor. * * @param nelx Number of elements (along the middle layer) in horizontal (x) direction. * @param nely Approximate equivalent number of elements in vertical (y) direction. * @param h Edge size (width == height) of elements along the weak layer. * * @param nfine * Extra number of fine layers around the middle layer. * By default the element size is kept smaller than the distance to the middle layer. */ FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1) { this->init(nelx, nely, h, nfine); } /** * Reconstruct class for given coordinates / connectivity. * * @tparam C e.g. `array_type::tensor` * @tparam E e.g. `array_type::tensor` - * @param coor Nodal coordinates ``[nnode, ndim]`` with ``ndim == 2``. - * @param conn Connectivity ``[nne, nne]`` with ``nne == 4``. + * @param coor Nodal coordinates `[nnode, ndim]` with `ndim == 2`. + * @param conn Connectivity `[nne, nne]` with `nne == 4`. */ template ::value, bool> = true> FineLayer(const C& coor, const E& conn) { this->init_by_mapping(coor, conn); } /** * Edge size in x-direction of a block, in units of #h, per row of blocks. * Note that a block is equal to an element except in refinement layers * where it contains three elements. * * @return List of size equal to the number of rows of blocks. */ const array_type::tensor& elemrow_nhx() const { return m_nhx; } /** * Edge size in y-direction of a block, in units of #h, per row of blocks. * Note that a block is equal to an element except in refinement layers * where it contains three elements. * * @return List of size equal to the number of rows of blocks. */ const array_type::tensor& elemrow_nhy() const { return m_nhy; } /** * Per row of blocks: * `-1`: normal layer * `0`: transition layer to match coarse and finer element on the previous/next row. * * @return List of size equal to the number of rows of blocks. */ const array_type::tensor& elemrow_type() const { return m_refine; } /** * Number of elements per row of blocks. * Note that a block is equal to an element except in refinement layers * where it contains three elements. * * @return List of size equal to the number of rows of blocks. */ const array_type::tensor& elemrow_nelem() const { return m_layer_nelx; } /** * Elements in the middle (fine) layer. * * @return List of element numbers. */ array_type::tensor elementsMiddleLayer() const { size_t nely = m_nhy.size(); size_t iy = (nely - 1) / 2; return m_startElem(iy) + xt::arange(m_layer_nelx(iy)); } /** * Elements along a layer. * * @return List of element numbers. */ array_type::tensor elementsLayer(size_t layer) const { GOOSEFEM_ASSERT(layer < m_layer_nelx.size()); size_t n = m_layer_nelx(layer); if (m_refine(layer) != -1) { n *= 4; } return m_startElem(layer) + xt::arange(n); } /** * Select region of elements from 'matrix' of element numbers. * * @return List of element numbers. */ array_type::tensor 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]) { array_type::tensor 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; } array_type::tensor 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; } array_type::tensor 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; } /** * Select region of elements from 'matrix' of element numbers around an element: - * square box with edge-size ``(2 * size + 1) * h``, around ``element``. + * square box with edge-size `(2 * size + 1) * h`, around `element`. * * @param e The element around which to select elements. * @param size Edge size of the square box encapsulating the selected element. * @param periodic Assume the mesh periodic. * @return List of elements. */ array_type::tensor elementgrid_around_ravel(size_t e, size_t size, bool periodic = true) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_layer_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)); } /** * Select region of elements from 'matrix' of element numbers around an element: - * left/right from ``element`` (on the same layer). + * left/right from `element` (on the same layer). * * @param e The element around which to select elements. * @param left Number of elements to select to the left. * @param right Number of elements to select to the right. * @param periodic Assume the mesh periodic. * @return List of elements. */ // - array_type::tensor elementgrid_leftright(size_t e, size_t left, size_t right, bool periodic = true) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_layer_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)); } /** * Mapping to 'roll' periodically in the x-direction, * * @return element mapping, such that: new_elemvar = elemvar[elem_map] */ array_type::tensor roll(size_t n) { auto conn = this->conn(); size_t nely = static_cast(m_nhy.size()); array_type::tensor 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_layer_nelx(iy) / m_layer_nelx(0)); size_t nel = m_layer_nelx(iy); // refinement if (m_refine(iy) != -1) { shift = n * (m_layer_nelx(iy) / m_layer_nelx(0)) * 4; nel = m_layer_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; } private: friend class RegularBase; friend class RegularBase2d; friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular; size_t nelx_impl() const { return xt::amax(m_layer_nelx)(); } size_t nely_impl() const { return xt::sum(m_nhy)(); } ElementType getElementType_impl() const { return ElementType::Quad4; } array_type::tensor coor_impl() const { // allocate output array_type::tensor 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 array_type::tensor 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 array_type::tensor x = xt::linspace(0.0, m_Lx, m_layer_nelx(iy) + 1); // add nodes of the bottom layer of this element for (size_t ix = 0; ix < m_layer_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_layer_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 array_type::tensor x = xt::linspace(0.0, m_Lx, m_layer_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_layer_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_layer_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy + 1); ++inode; } } return ret; } array_type::tensor conn_impl() const { // allocate output array_type::tensor 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_layer_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_layer_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_layer_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; } array_type::tensor nodesBottomEdge_impl() const { return m_startNode(0) + xt::arange(m_layer_nelx(0) + 1); } array_type::tensor nodesTopEdge_impl() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(m_layer_nelx(nely - 1) + 1); } array_type::tensor nodesLeftEdge_impl() const { size_t nely = m_nhy.size(); array_type::tensor 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; } array_type::tensor nodesRightEdge_impl() const { size_t nely = m_nhy.size(); array_type::tensor 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_layer_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_layer_nelx, xt::range(k, l)); return ret; } array_type::tensor nodesBottomOpenEdge_impl() const { return m_startNode(0) + xt::arange(1, m_layer_nelx(0)); } array_type::tensor nodesTopOpenEdge_impl() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(1, m_layer_nelx(nely - 1)); } array_type::tensor nodesLeftOpenEdge_impl() const { size_t nely = m_nhy.size(); array_type::tensor 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; } array_type::tensor nodesRightOpenEdge_impl() const { size_t nely = m_nhy.size(); array_type::tensor 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_layer_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_layer_nelx, xt::range(k, l - 1)); return ret; } size_t nodesBottomLeftCorner_impl() const { return m_startNode(0); } size_t nodesBottomRightCorner_impl() const { return m_startNode(0) + m_layer_nelx(0); } size_t nodesTopLeftCorner_impl() const { size_t nely = m_nhy.size(); return m_startNode(nely); } size_t nodesTopRightCorner_impl() const { size_t nely = m_nhy.size(); return m_startNode(nely) + m_layer_nelx(nely - 1); } double m_h; ///< See h() size_t m_nelem; ///< See nelem() size_t m_nnode; ///< See nnode() size_t m_nne; ///< See nne() size_t m_ndim; ///< See ndim() double m_Lx; ///< Mesh size in x-direction. array_type::tensor m_layer_nelx; ///< See elemrow_nelem(). array_type::tensor m_nhx; ///< See elemrow_nhx(). array_type::tensor m_nhy; ///< See elemrow_nhy(). array_type::tensor m_nnd; ///< num. nodes in main node layer (per node layer in "y") array_type::tensor m_refine; ///< See elemrow_type(). array_type::tensor m_startElem; ///< start element (per element layer in "y") array_type::tensor m_startNode; ///< start node (per node layer in "y") /** * @copydoc FineLayer::FineLayer(size_t, size_t, double, size_t) */ void init(size_t nelx, size_t nely, double h, size_t nfine = 1) { GOOSEFEM_ASSERT(nelx >= 1ul); GOOSEFEM_ASSERT(nely >= 1ul); m_h = h; m_ndim = 2; m_nne = 4; m_Lx = m_h * static_cast(nelx); // compute element size in y-direction (use symmetry, compute upper half) // temporary variables size_t nmin, ntot; array_type::tensor nhx = xt::ones({nely}); array_type::tensor nhy = xt::ones({nely}); array_type::tensor 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_layer_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_layer_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_layer_nelx(iy) + 1; } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { m_nnd(iy + 1) = m_layer_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_layer_nelx(i) + 1); } else { m_nnode += (m_layer_nelx(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_layer_nelx(i)); } else { m_nelem += (m_layer_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_layer_nelx(i) + 1); } else { m_nnode += (m_layer_nelx(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_layer_nelx(i)); } else { m_nelem += (m_layer_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_layer_nelx(nely - 1) + 1; } /** * @copydoc FineLayer::FineLayer(const C&, const E&) */ template void init_by_mapping(const C& coor, const E& conn) { GOOSEFEM_ASSERT(coor.dimension() == 2); GOOSEFEM_ASSERT(conn.dimension() == 2); 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)))); } }; /** * Mesh mappings. */ namespace Map { /** * Refine a Regular mesh: subdivide elements in several smaller elements. */ class RefineRegular { public: RefineRegular() = default; /** * Constructor. * * @param mesh the coarse mesh. * @param nx for each coarse element: number of fine elements in x-direction. * @param ny for each coarse element: number of fine elements in y-direction. */ RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny) : m_coarse(mesh), m_nx(nx), m_ny(ny) { m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h()); array_type::tensor elmat_coarse = m_coarse.elementgrid(); array_type::tensor 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))); } } } /** * For each coarse element: number of fine elements in x-direction. * * @return unsigned int (same as used in constructor) */ size_t nx() const { return m_nx; } /** * For each coarse element: number of fine elements in y-direction. * * @return unsigned int (same as used in constructor) */ size_t ny() const { return m_ny; } /** * Obtain the coarse mesh (copy of the mesh passed to the constructor). * @return mesh */ GooseFEM::Mesh::Quad4::Regular coarseMesh() const { return m_coarse; } /** * Obtain the fine mesh. * @return mesh */ GooseFEM::Mesh::Quad4::Regular fineMesh() const { return m_fine; } /** * Get element-mapping: elements of the fine mesh per element of the coarse mesh. * @return [nelem_coarse, nx() * ny()] */ const array_type::tensor& map() const { return m_coarse2fine; } /** * Obtain the coarse mesh (copy of the mesh passed to the constructor). * @return mesh */ [[deprecated]] GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const { return m_coarse; } /** * Obtain the fine mesh. * @return mesh */ [[deprecated]] GooseFEM::Mesh::Quad4::Regular getFineMesh() const { return m_fine; } /** * Get element-mapping: elements of the fine mesh per element of the coarse mesh. * @return [nelem_coarse, nx() * ny()] */ [[deprecated]] const array_type::tensor& getMap() const { return m_coarse2fine; } /** * Compute the mean of the quantity define on the fine mesh when mapped on the coarse mesh. * - * @tparam T type of the data (e.g. ``double``). + * @tparam T type of the data (e.g. `double`). * @tparam rank rank of the data. * @param data the data [nelem_fine, ...] * @return the average data of the coarse mesh [nelem_coarse, ...] */ template array_type::tensor meanToCoarse(const array_type::tensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_coarse2fine.shape(0); array_type::tensor ret = xt::empty(shape); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); auto d = xt::view(data, xt::keep(e)); xt::view(ret, i) = xt::mean(d, 0); } return ret; } /** * Compute the average of the quantity define on the fine mesh when mapped on the coarse mesh. * - * @tparam T type of the data (e.g. ``double``). + * @tparam T type of the data (e.g. `double`). * @tparam rank rank of the data. - * @tparam S type of the weights (e.g. ``double``). + * @tparam S type of the weights (e.g. `double`). * @param data the data [nelem_fine, ...] * @param weights the weights [nelem_fine, ...] * @return the average data of the coarse mesh [nelem_coarse, ...] */ template array_type::tensor averageToCoarse( const array_type::tensor& data, const array_type::tensor& weights) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_coarse2fine.shape(0); array_type::tensor ret = xt::empty(shape); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); array_type::tensor d = xt::view(data, xt::keep(e)); array_type::tensor w = xt::view(weights, xt::keep(e)); xt::view(ret, i) = xt::average(d, w, {0}); } return ret; } /** * Map element quantities to the fine mesh. * The mapping is a bit simplistic: no interpolation is involved. * The mapping is such that:: * * ret[e_fine, ...] <- data[e_coarse, ...] * - * @tparam T type of the data (e.g. ``double``). + * @tparam T type of the data (e.g. `double`). * @tparam rank rank of the data. * @param data the data. * @return mapped data. */ template array_type::tensor mapToFine(const array_type::tensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_coarse2fine.size(); array_type::tensor ret = xt::empty(shape); for (size_t e = 0; e < m_coarse2fine.shape(0); ++e) { for (size_t i = 0; i < m_coarse2fine.shape(1); ++i) { xt::view(ret, m_coarse2fine(e, i)) = xt::view(data, e); } } return ret; } private: GooseFEM::Mesh::Quad4::Regular m_coarse; ///< the coarse mesh GooseFEM::Mesh::Quad4::Regular m_fine; ///< the fine mesh size_t m_nx; ///< see nx() size_t m_ny; ///< see ny() array_type::tensor m_coarse2fine; ///< see getMap() }; /** * Map a FineLayer mesh to a Regular mesh. * The element size of the Regular corresponds to the smallest elements of the FineLayer mesh * (along the middle layer). */ class FineLayer2Regular { public: FineLayer2Regular() = default; /** * Constructors. * * @param mesh The FineLayer mesh. */ FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh) : m_finelayer(mesh) { // ------------ // Regular-mesh // ------------ m_regular = GooseFEM::Mesh::Quad4::Regular( xt::amax(m_finelayer.m_layer_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 array_type::tensor nhx = m_finelayer.m_nhx; array_type::tensor nhy = m_finelayer.m_nhy; array_type::tensor nelx = m_finelayer.m_layer_nelx; array_type::tensor start = m_finelayer.m_startElem; // 'matrix' of element numbers of the Regular-mesh array_type::tensor elementgrid = m_regular.elementgrid(); // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh array_type::tensor 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 array_type::tensor 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 array_type::tensor 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 array_type::tensor 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); } } } } } } } /** * Obtain the Regular mesh. * * @return mesh. */ GooseFEM::Mesh::Quad4::Regular regularMesh() const { return m_regular; } /** * Obtain the FineLayer mesh (copy of the mesh passed to the constructor). * * @return mesh. */ GooseFEM::Mesh::Quad4::FineLayer fineLayerMesh() const { return m_finelayer; } // elements of the Regular mesh per element of the FineLayer mesh // and the fraction by which the overlap is /** * Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh. * The number of Regular elements varies between elements of the FineLayer mesh. * * @return [nelem_finelayer, ?] */ std::vector> map() const { return m_elem_regular; } /** * To overlap fraction for each item in the mapping in map(). * * @return [nelem_finelayer, ?] */ std::vector> mapFraction() const { return m_frac_regular; } /** * Obtain the Regular mesh. * * @return mesh. */ [[deprecated]] GooseFEM::Mesh::Quad4::Regular getRegularMesh() const { return m_regular; } /** * Obtain the FineLayer mesh (copy of the mesh passed to the constructor). * * @return mesh. */ [[deprecated]] GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const { return m_finelayer; } // elements of the Regular mesh per element of the FineLayer mesh // and the fraction by which the overlap is /** * Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh. * The number of Regular elements varies between elements of the FineLayer mesh. * * @return [nelem_finelayer, ?] */ [[deprecated]] std::vector> getMap() const { return m_elem_regular; } /** * To overlap fraction for each item in the mapping in getMap(). * * @return [nelem_finelayer, ?] */ [[deprecated]] std::vector> getMapFraction() const { return m_frac_regular; } /** * Map element quantities to Regular. * The mapping is a bit simplistic: no interpolation is involved, the function just * accounts the fraction of overlap between the FineLayer element and the Regular element. * The mapping is such that:: * * ret[e_regular, ...] <- arg[e_finelayer, ...] * - * @tparam T type of the data (e.g. ``double``). + * @tparam T type of the data (e.g. `double`). * @tparam rank rank of the data. * @param data data. * @return mapped data. */ template array_type::tensor mapToRegular(const array_type::tensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_regular.nelem(); array_type::tensor ret = xt::zeros(shape); 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; } private: GooseFEM::Mesh::Quad4::FineLayer m_finelayer; ///< the FineLayer mesh to map GooseFEM::Mesh::Quad4::Regular m_regular; ///< the new Regular mesh to which to map std::vector> m_elem_regular; ///< see getMap() std::vector> m_frac_regular; ///< see getMapFraction() }; } // namespace Map } // namespace Quad4 } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/VectorPartitioned.h b/include/GooseFEM/VectorPartitioned.h index ac2937a..e73e115 100644 --- a/include/GooseFEM/VectorPartitioned.h +++ b/include/GooseFEM/VectorPartitioned.h @@ -1,648 +1,648 @@ /** * Methods to switch between storage types based on a mesh and DOFs that are partitioned in: * - unknown DOFs * - prescribed DOFs * * @file VectorPartitioned.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_VECTORPARTITIONED_H #define GOOSEFEM_VECTORPARTITIONED_H #include "Mesh.h" #include "Vector.h" #include "assertions.h" #include "config.h" namespace GooseFEM { /** * Class to switch between storage types, * based on a mesh and DOFs that are partitioned in: * * - unknown DOFs (iiu()), indicated with "u". * - prescribed DOFs (iip()), indicated with "p". * * To this end some internal re-ordering of the DOFs has to be done, as follows: * * iiu() -> arange(nnu()) * iip() -> nnu() + arange(nnp()) * * which is relevant only if you interact using partitioned DOF-lists ("dofval_u" or "dofval_p"). * * The "dofval", "nodevec", and "elemvec" are all stored in the 'normal' order. * * For reference: * * - "dofval": DOF values [#ndof]. * - "dofval_u": unknown DOF values, `== dofval[iiu()]`, [#nnu]. * - "dofval_p": prescribed DOF values, `== dofval[iip()]`, [#nnp]. * - "nodevec": nodal vectors [#nnode, #ndim]. * - "elemvec": nodal vectors stored per element [#nelem, #nne, #ndim]. */ class VectorPartitioned : public Vector { protected: array_type::tensor m_iiu; ///< See iiu() array_type::tensor m_iip; ///< See iip() size_t m_nnu; ///< See #nnu size_t m_nnp; ///< See #nnp /** * Renumbered DOFs per node, such that * * iiu = arange(nnu) * iip = nnu + arange(nnp) * * making is much simpler to slice. */ array_type::tensor m_part; public: VectorPartitioned() = default; /** * Constructor. * * @param conn connectivity [#nelem, #nne]. * @param dofs DOFs per node [#nnode, #ndim]. * @param iip prescribed DOFs [#nnp]. */ VectorPartitioned( const array_type::tensor& conn, const array_type::tensor& dofs, const array_type::tensor& iip) : Vector(conn, dofs), m_iip(iip) { GOOSEFEM_ASSERT(is_unique(iip)); m_iiu = xt::setdiff1d(m_dofs, m_iip); m_nnp = m_iip.size(); m_nnu = m_iiu.size(); m_part = Mesh::Reorder({m_iiu, m_iip}).apply(m_dofs); GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)()); } /** * @return Number of unknown DOFs. */ size_t nnu() const { return m_nnu; } /** * @return Number of prescribed DOFs. */ size_t nnp() const { return m_nnp; } /** * @return Unknown DOFs [#nnu]. */ const array_type::tensor& iiu() const { return m_iiu; } /** * @return Prescribed DOFs [#nnp]. */ const array_type::tensor& iip() const { return m_iip; } /** * Per DOF (see Vector::dofs()) list if unknown ("u"). * * @return Boolean "nodevec". */ array_type::tensor dofs_is_u() const { array_type::tensor ret = xt::zeros(this->shape_nodevec()); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { ret(m, i) = true; } } } return ret; } /** * Per DOF (see Vector::dofs()) list if prescribed ("p"). * * @return Boolean "nodevec". */ array_type::tensor dofs_is_p() const { array_type::tensor ret = xt::zeros(this->shape_nodevec()); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) >= m_nnu) { ret(m, i) = true; } } } return ret; } /** * Copy unknown DOFs from "nodevec" to another "nodvec": * * nodevec_dest[vector.dofs_is_u()] = nodevec_src * - * the other DOFs are taken from ``nodevec_dest``: + * the other DOFs are taken from `nodevec_dest`: * * nodevec_dest[vector.dofs_is_p()] = nodevec_dest * * @param nodevec_src input [#nnode, #ndim] * @param nodevec_dest input [#nnode, #ndim] * @return nodevec output [#nnode, #ndim] */ array_type::tensor Copy_u( const array_type::tensor& nodevec_src, const array_type::tensor& nodevec_dest) const { array_type::tensor ret = nodevec_dest; this->copy_u(nodevec_src, ret); return ret; } /** * Copy unknown DOFs from "nodevec" to another "nodvec": * * nodevec_dest[vector.dofs_is_u()] = nodevec_src * - * the other DOFs are taken from ``nodevec_dest``: + * the other DOFs are taken from `nodevec_dest`: * * nodevec_dest[vector.dofs_is_p()] = nodevec_dest * * @param nodevec_src input [#nnode, #ndim] * @param nodevec_dest input/output [#nnode, #ndim] */ void copy_u( const array_type::tensor& nodevec_src, array_type::tensor& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { nodevec_dest(m, i) = nodevec_src(m, i); } } } } /** * Copy prescribed DOFs from "nodevec" to another "nodvec": * * nodevec_dest[vector.dofs_is_p()] = nodevec_src * - * the other DOFs are taken from ``nodevec_dest``: + * the other DOFs are taken from `nodevec_dest`: * * nodevec_dest[vector.dofs_is_u()] = nodevec_dest * * @param nodevec_src input [#nnode, #ndim] * @param nodevec_dest input [#nnode, #ndim] * @return nodevec output [#nnode, #ndim] */ array_type::tensor Copy_p( const array_type::tensor& nodevec_src, const array_type::tensor& nodevec_dest) const { array_type::tensor ret = nodevec_dest; this->copy_p(nodevec_src, ret); return ret; } /** * Copy prescribed DOFs from "nodevec" to another "nodvec": * * nodevec_dest[vector.dofs_is_p()] = nodevec_src * - * the other DOFs are taken from ``nodevec_dest``: + * the other DOFs are taken from `nodevec_dest`: * * nodevec_dest[vector.dofs_is_u()] = nodevec_dest * * @param nodevec_src input [#nnode, #ndim] * @param nodevec_dest input/output [#nnode, #ndim] */ void copy_p( const array_type::tensor& nodevec_src, array_type::tensor& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) >= m_nnu) { nodevec_dest(m, i) = nodevec_src(m, i); } } } } /** * Combine unknown and prescribed "dofval" into a single "dofval" list. * * @param dofval_u input [#nnu] * @param dofval_p input [#nnp] * @return dofval output [#ndof] */ array_type::tensor DofsFromParitioned( const array_type::tensor& dofval_u, const array_type::tensor& dofval_p) const { array_type::tensor dofval = xt::empty({m_ndof}); this->dofsFromParitioned(dofval_u, dofval_p, dofval); return dofval; } /** * Combine unknown and prescribed "dofval" into a single "dofval" list. * * @param dofval_u input [#nnu] * @param dofval_p input [#nnp] * @param dofval output [#ndof] */ void dofsFromParitioned( const array_type::tensor& dofval_u, const array_type::tensor& dofval_p, array_type::tensor& dofval) const { GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); #pragma omp parallel for for (size_t d = 0; d < m_nnu; ++d) { dofval(m_iiu(d)) = dofval_u(d); } #pragma omp parallel for for (size_t d = 0; d < m_nnp; ++d) { dofval(m_iip(d)) = dofval_p(d); } } /** * Combine unknown and prescribed "dofval" into a single "dofval" list * and directly convert to "nodeval" without a temporary * (overwrite entries that occur more than once). * * @param dofval_u input [#nnu] * @param dofval_p input [#nnp] * @return nodevec output [#nnode, #ndim] */ array_type::tensor NodeFromPartitioned( const array_type::tensor& dofval_u, const array_type::tensor& dofval_p) const { array_type::tensor nodevec = xt::empty({m_nnode, m_ndim}); this->nodeFromPartitioned(dofval_u, dofval_p, nodevec); return nodevec; } /** * Combine unknown and prescribed "dofval" into a single "dofval" list * and directly convert to "nodeval" without a temporary * (overwrite entries that occur more than once). * * @param dofval_u input [#nnu] * @param dofval_p input [#nnp] * @param nodevec output [#nnode, #ndim] */ void nodeFromPartitioned( const array_type::tensor& dofval_u, const array_type::tensor& dofval_p, array_type::tensor& nodevec) const { GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { nodevec(m, i) = dofval_u(m_part(m, i)); } else { nodevec(m, i) = dofval_p(m_part(m, i) - m_nnu); } } } } /** * Combine unknown and prescribed "dofval" into a single "dofval" list * and directly convert to "elemvec" without a temporary * (overwrite entries that occur more than once). * * @param dofval_u input [#nnu] * @param dofval_p input [#nnp] * @return elemvec output [#nelem, #nne, #ndim] */ array_type::tensor ElementFromPartitioned( const array_type::tensor& dofval_u, const array_type::tensor& dofval_p) const { array_type::tensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); this->elementFromPartitioned(dofval_u, dofval_p, elemvec); return elemvec; } /** * Combine unknown and prescribed "dofval" into a single "dofval" list * and directly convert to "elemvec" without a temporary * (overwrite entries that occur more than once). * * @param dofval_u input [#nnu] * @param dofval_p input [#nnp] * @param elemvec output [#nelem, #nne, #ndim] */ void elementFromPartitioned( const array_type::tensor& dofval_u, const array_type::tensor& dofval_p, array_type::tensor& elemvec) const { GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m_conn(e, m), i) < m_nnu) { elemvec(e, m, i) = dofval_u(m_part(m_conn(e, m), i)); } else { elemvec(e, m, i) = dofval_p(m_part(m_conn(e, m), i) - m_nnu); } } } } } /** * Extract the unknown "dofval": * * dofval[iiu()] * * @param dofval input [#ndof] * @return dofval_u input [#nnu] */ array_type::tensor AsDofs_u(const array_type::tensor& dofval) const { array_type::tensor dofval_u = xt::empty({m_nnu}); this->asDofs_u(dofval, dofval_u); return dofval_u; } /** * Extract the unknown "dofval": * * dofval[iiu()] * * @param dofval input [#ndof] * @param dofval_u input [#nnu] */ void asDofs_u( const array_type::tensor& dofval, array_type::tensor& dofval_u) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); #pragma omp parallel for for (size_t d = 0; d < m_nnu; ++d) { dofval_u(d) = dofval(m_iiu(d)); } } /** * Convert "nodevec" to "dofval" (overwrite entries that occur more than once) * and extract the unknown "dofval" without a temporary. * * @param nodevec input [#nnode, #ndim] * @return dofval_u input [#nnu] */ array_type::tensor AsDofs_u(const array_type::tensor& nodevec) const { array_type::tensor dofval_u = xt::empty({m_nnu}); this->asDofs_u(nodevec, dofval_u); return dofval_u; } /** * Convert "nodevec" to "dofval" (overwrite entries that occur more than once) * and extract the unknown "dofval" without a temporary. * * @param nodevec input [#nnode, #ndim] * @param dofval_u input [#nnu] */ void asDofs_u( const array_type::tensor& nodevec, array_type::tensor& dofval_u) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); dofval_u.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { dofval_u(m_part(m, i)) = nodevec(m, i); } } } } /** * Convert "elemvec" to "dofval" (overwrite entries that occur more than once) * and extract the unknown "dofval" without a temporary. * * @param elemvec input [#nelem, #nne, #ndim] * @return dofval_u input [#nnu] */ array_type::tensor AsDofs_u(const array_type::tensor& elemvec) const { array_type::tensor dofval_u = xt::empty({m_nnu}); this->asDofs_u(elemvec, dofval_u); return dofval_u; } /** * Convert "elemvec" to "dofval" (overwrite entries that occur more than once) * and extract the unknown "dofval" without a temporary. * * @param elemvec input [#nelem, #nne, #ndim] * @param dofval_u input [#nnu] */ void asDofs_u( const array_type::tensor& elemvec, array_type::tensor& dofval_u) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); dofval_u.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m_conn(e, m), i) < m_nnu) { dofval_u(m_part(m_conn(e, m), i)) = elemvec(e, m, i); } } } } } /** * Extract the prescribed "dofval": * * dofval[iip()] * * @param dofval input [#ndof] * @return dofval_p input [#nnp] */ array_type::tensor AsDofs_p(const array_type::tensor& dofval) const { array_type::tensor dofval_p = xt::empty({m_nnp}); this->asDofs_p(dofval, dofval_p); return dofval_p; } /** * Extract the prescribed "dofval": * * dofval[iip()] * * @param dofval input [#ndof] * @param dofval_p input [#nnp] */ void asDofs_p( const array_type::tensor& dofval, array_type::tensor& dofval_p) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); #pragma omp parallel for for (size_t d = 0; d < m_nnp; ++d) { dofval_p(d) = dofval(m_iip(d)); } } /** * Convert "nodevec" to "dofval" (overwrite entries that occur more than once) * and extract the prescribed "dofval" without a temporary. * * @param nodevec input [#nnode, #ndim] * @return dofval_p input [#nnp] */ array_type::tensor AsDofs_p(const array_type::tensor& nodevec) const { array_type::tensor dofval_p = xt::empty({m_nnp}); this->asDofs_p(nodevec, dofval_p); return dofval_p; } /** * Convert "nodevec" to "dofval" (overwrite entries that occur more than once) * and extract the prescribed "dofval" without a temporary. * * @param nodevec input [#nnode, #ndim] * @param dofval_p input [#nnp] */ void asDofs_p( const array_type::tensor& nodevec, array_type::tensor& dofval_p) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); dofval_p.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) >= m_nnu) { dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i); } } } } /** * Convert "elemvec" to "dofval" (overwrite entries that occur more than once) * and extract the prescribed "dofval" without a temporary. * * @param elemvec input [#nelem, #nne, #ndim] * @return dofval_p input [#nnp] */ array_type::tensor AsDofs_p(const array_type::tensor& elemvec) const { array_type::tensor dofval_p = xt::empty({m_nnp}); this->asDofs_p(elemvec, dofval_p); return dofval_p; } /** * Convert "elemvec" to "dofval" (overwrite entries that occur more than once) * and extract the prescribed "dofval" without a temporary. * * @param elemvec input [#nelem, #nne, #ndim] * @param dofval_p input [#nnp] */ void asDofs_p( const array_type::tensor& elemvec, array_type::tensor& dofval_p) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); dofval_p.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m_conn(e, m), i) >= m_nnu) { dofval_p(m_part(m_conn(e, m), i) - m_nnu) = elemvec(e, m, i); } } } } } }; } // namespace GooseFEM #endif diff --git a/include/GooseFEM/version.h b/include/GooseFEM/version.h index 94902b0..769eabb 100644 --- a/include/GooseFEM/version.h +++ b/include/GooseFEM/version.h @@ -1,95 +1,95 @@ /** * Version information. * * @file version.h * @copyright Copyright 2017. Tom de Geus. All rights reserved. * @license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_VERSION_H #define GOOSEFEM_VERSION_H #include "config.h" /** * Current version. * * Either: * * - Configure using CMake at install time. Internally uses:: * * python -c "from setuptools_scm import get_version; print(get_version())" * * - Define externally using:: * * MYVERSION=`python -c "from setuptools_scm import get_version; print(get_version())"` * -DGOOSEFEM_VERSION="$MYVERSION" * - * From the root of this project. This is what ``setup.py`` does. + * From the root of this project. This is what `setup.py` does. * - * Note that both ``CMakeLists.txt`` and ``setup.py`` will construct the version using - * ``setuptools_scm``. Tip: use the environment variable ``SETUPTOOLS_SCM_PRETEND_VERSION`` to + * Note that both `CMakeLists.txt` and `setup.py` will construct the version using + * `setuptools_scm`. Tip: use the environment variable `SETUPTOOLS_SCM_PRETEND_VERSION` to * overwrite the automatic version. */ #ifndef GOOSEFEM_VERSION #define GOOSEFEM_VERSION "@PROJECT_VERSION@" #endif namespace GooseFEM { namespace detail { inline std::string unquote(const std::string& arg) { std::string ret = arg; ret.erase(std::remove(ret.begin(), ret.end(), '\"'), ret.end()); return ret; } } // namespace detail /** * Return version string, e.g. `"0.8.0"` * @return String. */ inline std::string version() { return detail::unquote(std::string(QUOTE(GOOSEFEM_VERSION))); } /** * Return versions of this library and of all of its dependencies. * The output is a list of strings, e.g.:: * * "goosefem=0.7.0", * "xtensor=0.20.1" * ... * * @return List of strings. */ inline std::vector version_dependencies() { std::vector ret; ret.push_back("goosefem=" + version()); ret.push_back( "xtensor=" + detail::unquote(std::string(QUOTE(XTENSOR_VERSION_MAJOR))) + "." + detail::unquote(std::string(QUOTE(XTENSOR_VERSION_MINOR))) + "." + detail::unquote(std::string(QUOTE(XTENSOR_VERSION_PATCH)))); #if defined(GOOSEFEM_EIGEN) || defined(EIGEN_WORLD_VERSION) ret.push_back( "eigen=" + detail::unquote(std::string(QUOTE(EIGEN_WORLD_VERSION))) + "." + detail::unquote(std::string(QUOTE(EIGEN_MAJOR_VERSION))) + "." + detail::unquote(std::string(QUOTE(EIGEN_MINOR_VERSION)))); #endif return ret; } } // namespace GooseFEM #endif