9#ifndef GOOSEFEM_ELEMENT_H
10#define GOOSEFEM_ELEMENT_H
35 size_t nelem = conn.shape(0);
36 size_t nne = conn.shape(1);
37 size_t ndim = nodevec.shape(1);
41#pragma omp parallel for
42 for (
size_t e = 0; e < nelem; ++e) {
43 for (
size_t m = 0; m < nne; ++m) {
44 for (
size_t i = 0; i < ndim; ++i) {
45 elemvec(e, m, i) = nodevec(conn(e, m), i);
66 size_t nelem = conn.shape(0);
67 size_t nne = conn.shape(1);
68 size_t ndim = elemvec.shape(2);
69 size_t nnode = xt::amax(conn)() + 1;
76 for (
size_t e = 0; e < nelem; ++e) {
77 for (
size_t m = 0; m < nne; ++m) {
78 for (
size_t i = 0; i < ndim; ++i) {
79 nodevec(conn(e, m), i) += elemvec(e, m, i);
96 size_t ndof = xt::amax(dofs)() + 1;
100 for (
auto& i : dofs) {
104 for (
auto& i : dofs) {
105 if (exists[i] == 0) {
124 size_t nelem = elemmat.shape(0);
125 size_t N = elemmat.shape(1);
127 double eps = std::numeric_limits<double>::epsilon();
129#pragma omp parallel for
130 for (
size_t e = 0; e < nelem; ++e) {
131 for (
size_t i = 0; i < N; ++i) {
132 for (
size_t j = 0; j < N; ++j) {
134 if (std::abs(elemmat(e, i, j)) > eps) {
163 return derived_cast().m_nelem;
203 return derived_cast().m_nip;
213 template <
class T,
class R>
226 template <
size_t rank,
class T>
229 return GooseFEM::AsTensor<rank>(arg, derived_cast().m_tdim);
252 return std::array<size_t, 3>{derived_cast().m_nelem, D::s_nne, D::s_ndim};
263 return std::array<size_t, 3>{derived_cast().m_nelem, D::s_nne, arg};
273 return std::array<size_t, 3>{
274 derived_cast().m_nelem, D::s_nne * D::s_ndim, D::s_nne * D::s_ndim};
284 template <
size_t rank = 0>
287 std::array<size_t, 2 + rank> shape;
288 shape[0] = derived_cast().m_nelem;
289 shape[1] = derived_cast().m_nip;
290 std::fill(shape.begin() + 2, shape.end(), derived_cast().m_tdim);
303 std::vector<size_t> shape(2 + rank);
304 shape[0] = derived_cast().m_nelem;
305 shape[1] = derived_cast().m_nip;
306 std::fill(shape.begin() + 2, shape.end(), derived_cast().m_tdim);
319 template <
size_t trank>
320 auto shape_qtensor(
size_t rank,
size_t arg)
const -> std::array<size_t, 2 + trank>
323 std::array<size_t, 2 + trank> shape;
324 shape[0] = derived_cast().m_nelem;
325 shape[1] = derived_cast().m_nip;
326 std::fill(shape.begin() + 2, shape.end(), arg);
340 std::vector<size_t> shape(2 + rank);
341 shape[0] = derived_cast().m_nelem;
342 shape[1] = derived_cast().m_nip;
343 std::fill(shape.begin() + 2, shape.end(), arg);
353 return std::array<size_t, 2>{derived_cast().m_nelem, derived_cast().m_nip};
362 return std::array<size_t, 3>{derived_cast().m_nelem, derived_cast().m_nip, D::s_tdim};
372 return std::array<size_t, 3>{derived_cast().m_nelem, derived_cast().m_nip, arg};
440 template <
size_t rank = 0,
class R>
455 template <
size_t rank = 0,
class R>
490 auto ret = xt::xarray<R>::from_shape(this->
shape_qtensor(rank));
505 return this->allocate_qtensor<0, R>();
518 return this->allocate_qtensor<0, R>(val);
564 xt::noalias(derived_cast().m_x) = x;
565 derived_cast().compute_dN_impl();
572 auto GradN() const -> const array_type::tensor<
double, 4>&
574 return derived_cast().m_dNx;
581 auto dV() const -> const array_type::tensor<
double, 2>&
583 return derived_cast().m_vol;
597 size_t n = elemvec.shape(2);
599 derived_cast().interpQuad_vector_impl(elemvec, qvector);
609 template <
class T,
class R>
612 derived_cast().interpQuad_vector_impl(elemvec, qvector);
635 derived_cast().gradN_vector_impl(elemvec, qtensor);
645 template <
class T,
class R>
648 derived_cast().gradN_vector_impl(elemvec, qtensor);
667 derived_cast().gradN_vector_T_impl(elemvec, qtensor);
677 template <
class T,
class R>
680 derived_cast().gradN_vector_T_impl(elemvec, qtensor);
700 derived_cast().symGradN_vector_impl(elemvec, qtensor);
710 template <
class T,
class R>
713 derived_cast().symGradN_vector_impl(elemvec, qtensor);
731 size_t n = qvector.shape(2);
733 derived_cast().int_N_vector_dV_impl(qvector, elemvec);
743 template <
class T,
class R>
746 derived_cast().int_N_vector_dV_impl(qvector, elemvec);
771 derived_cast().int_N_scalar_NT_dV_impl(qscalar, elemmat);
781 template <
class T,
class R>
784 derived_cast().int_N_scalar_NT_dV_impl(qscalar, elemmat);
808 derived_cast().int_gradN_dot_tensor2_dV_impl(qtensor, elemvec);
818 template <
class T,
class R>
821 derived_cast().int_gradN_dot_tensor2_dV_impl(qtensor, elemvec);
851 derived_cast().int_gradN_dot_tensor4_dot_gradNT_dV_impl(qtensor, elemmat);
861 template <
class T,
class R>
864 derived_cast().int_gradN_dot_tensor4_dot_gradNT_dV_impl(qtensor, elemmat);
873 derived_cast().compute_dN_impl();
887 friend class QuadratureBase<D>;
889 template <
class T,
class R>
890 void interpQuad_vector_impl(
const T& elemvec, R& qvector)
const
892 size_t n = elemvec.shape(2);
896 auto nelem = derived_cast().m_nelem;
897 auto nip = derived_cast().m_nip;
898 auto& N = derived_cast().m_N;
902#pragma omp parallel for
903 for (
size_t e = 0; e <
nelem; ++e) {
905 auto fq = &elemvec(e, 0, 0);
907 for (
size_t q = 0; q <
nip; ++q) {
910 auto tq = &qvector(e, q, 0);
912 for (
size_t m = 0; m < D::s_nne; ++m) {
913 for (
size_t i = 0; i < n; ++i) {
914 tq[i] += Nq[m] * fq[m * n + i];
921 template <
class T,
class R>
922 void gradN_vector_impl(
const T& elemvec, R& qtensor)
const
925 GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->
template shape_qtensor<2>()));
927 auto nelem = derived_cast().m_nelem;
928 auto nip = derived_cast().m_nip;
929 auto& dNx = derived_cast().m_dNx;
933#pragma omp parallel for
934 for (
size_t e = 0; e <
nelem; ++e) {
936 auto ue = xt::adapt(&elemvec(e, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
938 for (
size_t q = 0; q <
nip; ++q) {
940 auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
941 auto graduq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<D::s_tdim, D::s_tdim>());
943 for (
size_t m = 0; m < D::s_nne; ++m) {
944 for (
size_t i = 0; i < D::s_ndim; ++i) {
945 for (
size_t j = 0; j < D::s_ndim; ++j) {
946 graduq(i, j) += dNxq(m, i) * ue(m, j);
954 template <
class T,
class R>
955 void gradN_vector_T_impl(
const T& elemvec, R& qtensor)
const
958 GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->
template shape_qtensor<2>()));
960 auto nelem = derived_cast().m_nelem;
961 auto nip = derived_cast().m_nip;
962 auto& dNx = derived_cast().m_dNx;
966#pragma omp parallel for
967 for (
size_t e = 0; e <
nelem; ++e) {
969 auto ue = xt::adapt(&elemvec(e, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
971 for (
size_t q = 0; q <
nip; ++q) {
973 auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
974 auto graduq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<D::s_tdim, D::s_tdim>());
976 for (
size_t m = 0; m < D::s_nne; ++m) {
977 for (
size_t i = 0; i < D::s_ndim; ++i) {
978 for (
size_t j = 0; j < D::s_ndim; ++j) {
979 graduq(j, i) += dNxq(m, i) * ue(m, j);
987 template <
class T,
class R>
988 void symGradN_vector_impl(
const T& elemvec, R& qtensor)
const
991 GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->
template shape_qtensor<2>()));
993 auto nelem = derived_cast().m_nelem;
994 auto nip = derived_cast().m_nip;
995 auto& dNx = derived_cast().m_dNx;
999#pragma omp parallel for
1000 for (
size_t e = 0; e <
nelem; ++e) {
1002 auto ue = xt::adapt(&elemvec(e, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1004 for (
size_t q = 0; q <
nip; ++q) {
1006 auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1007 auto epsq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<D::s_tdim, D::s_tdim>());
1009 for (
size_t m = 0; m < D::s_nne; ++m) {
1010 for (
size_t i = 0; i < D::s_ndim; ++i) {
1011 for (
size_t j = 0; j < D::s_ndim; ++j) {
1012 epsq(i, j) += 0.5 * dNxq(m, i) * ue(m, j);
1013 epsq(j, i) += 0.5 * dNxq(m, i) * ue(m, j);
1021 template <
class T,
class R>
1022 void int_N_vector_dV_impl(
const T& qvector, R& elemvec)
const
1024 size_t n = qvector.shape(2);
1028 auto nelem = derived_cast().m_nelem;
1029 auto nip = derived_cast().m_nip;
1030 auto& N = derived_cast().m_N;
1031 auto& vol = derived_cast().m_vol;
1035#pragma omp parallel for
1036 for (
size_t e = 0; e <
nelem; ++e) {
1038 auto f = &elemvec(e, 0, 0);
1040 for (
size_t q = 0; q <
nip; ++q) {
1043 auto tq = &qvector(e, q, 0);
1044 auto& volq = vol(e, q);
1046 for (
size_t m = 0; m < D::s_nne; ++m) {
1047 for (
size_t i = 0; i < n; ++i) {
1048 f[m * n + i] += Ne[m] * tq[i] * volq;
1055 template <
class T,
class R>
1056 void int_N_scalar_NT_dV_impl(
const T& qscalar, R& elemmat)
const
1061 auto nelem = derived_cast().m_nelem;
1062 auto nip = derived_cast().m_nip;
1063 auto& N = derived_cast().m_N;
1064 auto& vol = derived_cast().m_vol;
1068#pragma omp parallel for
1069 for (
size_t e = 0; e <
nelem; ++e) {
1071 auto Me = xt::adapt(
1072 &elemmat(e, 0, 0), xt::xshape<D::s_nne * D::s_ndim, D::s_nne * D::s_ndim>());
1074 for (
size_t q = 0; q <
nip; ++q) {
1076 auto Ne = xt::adapt(&N(q, 0), xt::xshape<D::s_nne>());
1077 auto& volq = vol(e, q);
1078 auto& rho = qscalar(e, q);
1081 for (
size_t m = 0; m < D::s_nne; ++m) {
1082 for (
size_t n = 0; n < D::s_nne; ++n) {
1083 for (
size_t i = 0; i < D::s_ndim; ++i) {
1084 Me(m * D::s_ndim + i, n * D::s_ndim + i) += Ne(m) * rho * Ne(n) * volq;
1092 template <
class T,
class R>
1093 void int_gradN_dot_tensor2_dV_impl(
const T& qtensor, R& elemvec)
const
1095 GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->
template shape_qtensor<2>()));
1098 auto nelem = derived_cast().m_nelem;
1099 auto nip = derived_cast().m_nip;
1100 auto& dNx = derived_cast().m_m_dNx;
1101 auto& vol = derived_cast().m_vol;
1105#pragma omp parallel for
1106 for (
size_t e = 0; e <
nelem; ++e) {
1108 auto fe = xt::adapt(&elemvec(e, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1110 for (
size_t q = 0; q <
nip; ++q) {
1112 auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1113 auto sigq = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<D::s_tdim, D::s_tdim>());
1114 auto& volq = vol(e, q);
1116 for (
size_t m = 0; m < D::s_nne; ++m) {
1117 for (
size_t i = 0; i < D::s_ndim; ++i) {
1118 for (
size_t j = 0; j < D::s_ndim; ++j) {
1119 fe(m, j) += dNxq(m, i) * sigq(i, j) * volq;
1127 template <
class T,
class R>
1128 void int_gradN_dot_tensor4_dot_gradNT_dV_impl(
const T& qtensor, R& elemmat)
const
1130 GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->
template shape_qtensor<4>()));
1133 auto nelem = derived_cast().m_nelem;
1134 auto nip = derived_cast().m_nip;
1135 auto& dNx = derived_cast().m_dNx;
1136 auto& vol = derived_cast().m_vol;
1140#pragma omp parallel for
1141 for (
size_t e = 0; e <
nelem; ++e) {
1144 &elemmat(e, 0, 0), xt::xshape<D::s_nne * D::s_ndim, D::s_nne * D::s_ndim>());
1146 for (
size_t q = 0; q <
nip; ++q) {
1148 auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1149 auto Cq = xt::adapt(
1150 &qtensor(e, q, 0, 0, 0, 0),
1151 xt::xshape<D::s_tdim, D::s_tdim, D::s_tdim, D::s_tdim>());
1152 auto& volq = vol(e, q);
1154 for (
size_t m = 0; m < D::s_nne; ++m) {
1155 for (
size_t n = 0; n < D::s_nne; ++n) {
1156 for (
size_t i = 0; i < D::s_ndim; ++i) {
1157 for (
size_t j = 0; j < D::s_ndim; ++j) {
1158 for (
size_t k = 0; k < D::s_ndim; ++k) {
1159 for (
size_t l = 0; l < D::s_ndim; ++l) {
1160 K(m * D::s_ndim + j, n * D::s_ndim + k) +=
1161 dNxq(m, i) * Cq(i, j, k, l) * dNxq(n, l) * volq;
1172 void compute_dN_impl()
1174 auto nelem = derived_cast().m_nelem;
1175 auto nip = derived_cast().m_nip;
1176 auto& vol = derived_cast().m_vol;
1177 auto&
w = derived_cast().m_w;
1178 auto& dNxi = derived_cast().m_dNxi;
1179 auto& dNx = derived_cast().m_dNx;
1180 auto& x = derived_cast().m_x;
1186 auto J = array_type::tensor<double, 2>::from_shape({D::s_ndim, D::s_ndim});
1187 auto Jinv = array_type::tensor<double, 2>::from_shape({D::s_ndim, D::s_ndim});
1190 for (
size_t e = 0; e <
nelem; ++e) {
1192 auto xe = xt::adapt(&x(e, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1194 for (
size_t q = 0; q <
nip; ++q) {
1196 auto dNxiq = xt::adapt(&dNxi(q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1197 auto dNxq = xt::adapt(&dNx(e, q, 0, 0), xt::xshape<D::s_nne, D::s_ndim>());
1201 for (
size_t m = 0; m < D::s_nne; ++m) {
1202 for (
size_t i = 0; i < D::s_ndim; ++i) {
1203 for (
size_t j = 0; j < D::s_ndim; ++j) {
1204 J(i, j) += dNxiq(m, i) * xe(m, j);
1209 double Jdet = detail::tensor<D::s_ndim>::inv(J, Jinv);
1211 for (
size_t m = 0; m < D::s_nne; ++m) {
1212 for (
size_t i = 0; i < D::s_ndim; ++i) {
1213 for (
size_t j = 0; j < D::s_ndim; ++j) {
1214 dNxq(m, i) += Jinv(i, j) * dNxiq(m, i);
1219 vol(e, q) =
w(q) * Jdet;
Common allocation methods.
CRTP base class for interpolation and quadrature for a generic element in Cartesian coordinates.
auto Int_gradN_dot_tensor2_dV(const T &qtensor) const -> array_type::tensor< double, 3 >
Element-by-element: integral of the dot product of the shape function gradients with a second order t...
auto Int_N_scalar_NT_dV(const T &qscalar) const -> array_type::tensor< double, 3 >
Element-by-element: integral of the scalar product of the shape function with a scalar.
void gradN_vector(const T &elemvec, R &qtensor) const
Same as GradN_vector(), but writing to preallocated return.
auto SymGradN_vector(const T &elemvec) const -> array_type::tensor< double, 4 >
The symmetric output of GradN_vector().
void symGradN_vector(const T &elemvec, R &qtensor) const
Same as SymGradN_vector(), but writing to preallocated return.
void int_gradN_dot_tensor2_dV(const T &qtensor, R &elemvec) const
Same as Int_gradN_dot_tensor2_dV(), but writing to preallocated return.
auto InterpQuad_vector(const T &elemvec) const -> array_type::tensor< double, 3 >
Interpolate element vector and evaluate at each quadrature point.
void int_gradN_dot_tensor4_dot_gradNT_dV(const T &qtensor, R &elemmat) const
Same as Int_gradN_dot_tensor4_dot_gradNT_dV(), but writing to preallocated return.
auto dV() const -> const array_type::tensor< double, 2 > &
Integration volume.
void gradN_vector_T(const T &elemvec, R &qtensor) const
Same as GradN_vector_T(), but writing to preallocated return.
auto GradN_vector_T(const T &elemvec) const -> array_type::tensor< double, 4 >
The transposed output of GradN_vector().
D derived_type
Underlying type.
void int_N_vector_dV(const T &qvector, R &elemvec) const
Same as Int_N_vector_dV(), but writing to preallocated return.
auto Int_gradN_dot_tensor4_dot_gradNT_dV(const T &qtensor) const -> array_type::tensor< double, 3 >
Element-by-element: integral of the dot products of the shape function gradients with a fourth order ...
auto GradN_vector(const T &elemvec) const -> array_type::tensor< double, 4 >
Element-by-element: dyadic product of the shape function gradients and a nodal vector.
void interpQuad_vector(const T &elemvec, R &qvector) const
Same as InterpQuad_vector(), but writing to preallocated return.
void compute_dN()
Update the shape function gradients (called when the nodal positions are updated).
void int_N_scalar_NT_dV(const T &qscalar, R &elemmat) const
Same as Int_N_scalar_NT_dV(), but writing to preallocated return.
void update_x(const T &x)
Update the nodal positions.
auto GradN() const -> const array_type::tensor< double, 4 > &
Shape function gradients (in global coordinates).
auto Int_N_vector_dV(const T &qvector) const -> array_type::tensor< double, 3 >
Element-by-element: integral of a continuous vector-field.
CRTP base class for quadrature.
auto shape_qscalar() const -> std::array< size_t, 2 >
Get the shape of a "qscalar" (a "qtensor" of rank 0)
auto shape_qtensor(size_t rank, size_t arg) const -> std::vector< size_t >
Get the shape of a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor,...
auto shape_elemvec(size_t arg) const -> std::array< size_t, 3 >
Get the shape of an "elemvec".
auto shape_qvector() const -> std::array< size_t, 3 >
Get the shape of a "qvector" (a "qtensor" of rank 1)
auto nip() const
Number of integration points.
auto allocate_elemmat() const
Get an allocated array_type::tensor to store a "elemmat".
auto nelem() const
Number of elements.
auto nne() const
Number of nodes per element.
void asTensor(const T &arg, R &ret) const
Convert "qscalar" to "qtensor" of certain rank.
auto shape_qtensor() const -> std::array< size_t, 2+rank >
Get the shape of a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor,...
auto allocate_elemvec() const
Get an allocated array_type::tensor to store a "elemvec".
auto shape_elemmat() const -> std::array< size_t, 3 >
Get the shape of an "elemmat".
auto allocate_elemvec(R val) const
Get an allocated and initialised xt::xarray to store a "elemvec".
auto shape_qvector(size_t arg) const -> std::array< size_t, 3 >
Get the shape of a "qvector" (a "qtensor" of rank 1)
auto allocate_elemmat(R val) const
Get an allocated and initialised xt::xarray to store a "elemmat".
auto shape_elemvec() const -> std::array< size_t, 3 >
Get the shape of an "elemvec".
auto allocate_qscalar(R val) const
Get an allocated and initialised xt::xarray to store a "qscalar" (a "qtensor" of rank 0).
auto AsTensor(const T &arg) const
Convert "qscalar" to "qtensor" of certain rank.
auto tdim() const
Number of dimensions for integration point tensors.
auto allocate_qtensor(size_t rank) const
Get an allocated xt::xarray to store a "qtensor" of a certain rank (0 = scalar, 1,...
auto shape_qtensor(size_t rank) const -> std::vector< size_t >
Get the shape of a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor,...
auto allocate_qtensor() const
Get an allocated array_type::tensor to store a "qtensor" of a certain rank (0 = scalar,...
auto allocate_qtensor(size_t rank, R val) const
Get an allocated and initialised xt::xarray to store a "qtensor" of a certain rank (0 = scalar,...
auto shape_qtensor(size_t rank, size_t arg) const -> std::array< size_t, 2+trank >
Get the shape of a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor,...
auto ndim() const
Number of dimensions for node vectors.
auto allocate_qscalar() const
Get an allocated array_type::tensor to store a "qscalar" (a "qtensor" of rank 0).
auto AsTensor(size_t rank, const T &arg) const
Convert "qscalar" to "qtensor" of certain rank.
auto allocate_qtensor(R val) const
Get an allocated and initialised array_type::tensor to store a "qtensor" of a certain rank (0 = scala...
D derived_type
Underlying type.
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
array_type::tensor< double, 1 > w()
Integration point weights.
bool isDiagonal(const array_type::tensor< double, 3 > &elemmat)
Check that all of the matrices stored per elemmat (shape: [nelem, nne * ndim, nne * ndim]) are diagon...
array_type::tensor< double, 3 > asElementVector(const array_type::tensor< size_t, 2 > &conn, const array_type::tensor< double, 2 > &nodevec)
Convert nodal vector with ("nodevec", shape:[nnode, ndim]) to nodal vector stored per element ("elemv...
bool isSequential(const E &dofs)
Check that DOFs leave no holes.
array_type::tensor< double, 2 > assembleNodeVector(const array_type::tensor< size_t, 2 > &conn, const array_type::tensor< double, 3 > &elemvec)
Assemble nodal vector stored per element ("elemvec", shape [nelem, nne, ndim]) to nodal vector ("node...
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Toolbox to perform finite element computations.
auto AsTensor(const T &arg, const S &shape)
"Broadcast" a scalar stored in an array (e.g.
void asTensor(const T &arg, R &ret)
"Broadcast" a scalar stored in an array (e.g.