12#ifndef GOOSEFEM_ELEMENTQUAD4PLANAR_H
13#define GOOSEFEM_ELEMENTQUAD4PLANAR_H
80 template <
class T,
class X,
class W>
87 m_nelem = m_x.shape(0);
88 m_N = xt::empty<double>({m_nip, s_nne});
89 m_dNxi = xt::empty<double>({m_nip, s_nne, s_ndim});
92 for (
size_t q = 0; q < m_nip; ++q) {
93 m_N(q, 0) = 0.25 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1));
94 m_N(q, 1) = 0.25 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1));
95 m_N(q, 2) = 0.25 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1));
96 m_N(q, 3) = 0.25 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1));
99 for (
size_t q = 0; q < m_nip; ++q) {
101 m_dNxi(q, 0, 0) = -0.25 * (1.0 - xi(q, 1));
102 m_dNxi(q, 1, 0) = +0.25 * (1.0 - xi(q, 1));
103 m_dNxi(q, 2, 0) = +0.25 * (1.0 + xi(q, 1));
104 m_dNxi(q, 3, 0) = -0.25 * (1.0 + xi(q, 1));
106 m_dNxi(q, 0, 1) = -0.25 * (1.0 - xi(q, 0));
107 m_dNxi(q, 1, 1) = -0.25 * (1.0 + xi(q, 0));
108 m_dNxi(q, 2, 1) = +0.25 * (1.0 + xi(q, 0));
109 m_dNxi(q, 3, 1) = +0.25 * (1.0 - xi(q, 0));
119 m_dNx = xt::empty<double>({m_nelem, m_nip, s_nne, s_ndim});
122 this->compute_dN_impl();
129 template <
class T,
class R>
130 void gradN_vector_impl(
const T& elemvec, R& qtensor)
const
137#pragma omp parallel for
138 for (
size_t e = 0; e < m_nelem; ++e) {
140 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
142 for (
size_t q = 0; q < m_nip; ++q) {
144 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
145 auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_tdim, s_tdim>());
148 gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) +
150 gradu(0, 1) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) +
152 gradu(1, 0) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) +
154 gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) +
160 template <
class T,
class R>
161 void gradN_vector_T_impl(
const T& elemvec, R& qtensor)
const
168#pragma omp parallel for
169 for (
size_t e = 0; e < m_nelem; ++e) {
171 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
173 for (
size_t q = 0; q < m_nip; ++q) {
175 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
176 auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_tdim, s_tdim>());
179 gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) +
181 gradu(1, 0) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) +
183 gradu(0, 1) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) +
185 gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) +
191 template <
class T,
class R>
192 void symGradN_vector_impl(
const T& elemvec, R& qtensor)
const
199#pragma omp parallel for
200 for (
size_t e = 0; e < m_nelem; ++e) {
202 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
204 for (
size_t q = 0; q < m_nip; ++q) {
206 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
207 auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_tdim, s_tdim>());
211 eps(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) +
213 eps(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) +
215 eps(0, 1) = 0.5 * (dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) +
216 dNx(3, 0) * u(3, 1) + dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) +
217 dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0));
218 eps(1, 0) = eps(0, 1);
223 template <
class T,
class R>
224 void int_N_scalar_NT_dV_impl(
const T& qscalar, R& elemmat)
const
231#pragma omp parallel for
232 for (
size_t e = 0; e < m_nelem; ++e) {
234 auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape<s_nne * s_ndim, s_nne * s_ndim>());
236 for (
size_t q = 0; q < m_nip; ++q) {
238 auto N = xt::adapt(&m_N(q, 0), xt::xshape<s_nne>());
239 auto& vol = m_vol(e, q);
240 auto& rho = qscalar(e, q);
243 for (
size_t m = 0; m < s_nne; ++m) {
244 for (
size_t n = 0; n < s_nne; ++n) {
245 M(m * s_ndim + 0, n * s_ndim + 0) += N(m) * rho * N(n) * vol;
246 M(m * s_ndim + 1, n * s_ndim + 1) += N(m) * rho * N(n) * vol;
253 template <
class T,
class R>
254 void int_gradN_dot_tensor2_dV_impl(
const T& qtensor, R& elemvec)
const
261#pragma omp parallel for
262 for (
size_t e = 0; e < m_nelem; ++e) {
264 auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
266 for (
size_t q = 0; q < m_nip; ++q) {
268 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
269 auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_tdim, s_tdim>());
270 auto& vol = m_vol(e, q);
272 for (
size_t m = 0; m < s_nne; ++m) {
273 f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0)) * vol;
274 f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1)) * vol;
280 void compute_dN_impl()
284 array_type::tensor<double, 2> J = xt::empty<double>({2, 2});
285 array_type::tensor<double, 2> Jinv = xt::empty<double>({2, 2});
288 for (
size_t e = 0; e < m_nelem; ++e) {
290 auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape<s_nne, s_ndim>());
292 for (
size_t q = 0; q < m_nip; ++q) {
294 auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape<s_nne, s_ndim>());
295 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
298 J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) +
299 dNxi(3, 0) * x(3, 0);
300 J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) +
301 dNxi(3, 0) * x(3, 1);
302 J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) +
303 dNxi(3, 1) * x(3, 0);
304 J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) +
305 dNxi(3, 1) * x(3, 1);
307 double Jdet = detail::tensor<2>::inv(J, Jinv);
310 for (
size_t m = 0; m < s_nne; ++m) {
311 dNx(m, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1);
312 dNx(m, 1) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1);
315 m_vol(e, q) = m_w(q) * Jdet * m_thick;
321 constexpr static size_t s_nne = 4;
322 constexpr static size_t s_ndim = 2;
323 constexpr static size_t s_tdim = 3;
327 array_type::tensor<double, 3> m_x;
328 array_type::tensor<double, 1> m_w;
329 array_type::tensor<double, 2> m_xi;
330 array_type::tensor<double, 2> m_N;
331 array_type::tensor<double, 3> m_dNxi;
332 array_type::tensor<double, 4> m_dNx;
333 array_type::tensor<double, 2> m_vol;
Interpolation and quadrature.
QuadraturePlanar(const T &x, const X &xi, const W &w, double thick=1.0)
Constructor with custom integration.
QuadraturePlanar(const T &x, double thick=1.0)
Constructor: use default Gauss integration.
CRTP base class for interpolation and quadrature for a generic element in Cartesian coordinates.
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_elemmat() const -> std::array< size_t, 3 >
Get the shape of an "elemmat".
auto shape_elemvec() const -> std::array< size_t, 3 >
Get the shape of an "elemvec".
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
@ Quad4
Quadrilateral: 4-noded element in 2-d.
Toolbox to perform finite element computations.