10#ifndef GOOSEFEM_ELEMENTQUAD4_H
11#define GOOSEFEM_ELEMENTQUAD4_H
62 xi(0, 0) = -1.0 / std::sqrt(3.0);
63 xi(0, 1) = -1.0 / std::sqrt(3.0);
64 xi(1, 0) = +1.0 / std::sqrt(3.0);
65 xi(1, 1) = -1.0 / std::sqrt(3.0);
66 xi(2, 0) = +1.0 / std::sqrt(3.0);
67 xi(2, 1) = +1.0 / std::sqrt(3.0);
68 xi(3, 0) = -1.0 / std::sqrt(3.0);
69 xi(3, 1) = +1.0 / std::sqrt(3.0);
276 template <
class T,
class X,
class W>
283 m_nelem = m_x.shape(0);
284 m_N = xt::empty<double>({m_nip, s_nne});
285 m_dNxi = xt::empty<double>({m_nip, s_nne, s_ndim});
287 for (
size_t q = 0; q < m_nip; ++q) {
288 m_N(q, 0) = 0.25 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1));
289 m_N(q, 1) = 0.25 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1));
290 m_N(q, 2) = 0.25 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1));
291 m_N(q, 3) = 0.25 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1));
294 for (
size_t q = 0; q < m_nip; ++q) {
296 m_dNxi(q, 0, 0) = -0.25 * (1.0 - xi(q, 1));
297 m_dNxi(q, 1, 0) = +0.25 * (1.0 - xi(q, 1));
298 m_dNxi(q, 2, 0) = +0.25 * (1.0 + xi(q, 1));
299 m_dNxi(q, 3, 0) = -0.25 * (1.0 + xi(q, 1));
301 m_dNxi(q, 0, 1) = -0.25 * (1.0 - xi(q, 0));
302 m_dNxi(q, 1, 1) = -0.25 * (1.0 + xi(q, 0));
303 m_dNxi(q, 2, 1) = +0.25 * (1.0 + xi(q, 0));
304 m_dNxi(q, 3, 1) = +0.25 * (1.0 - xi(q, 0));
314 m_dNx = xt::empty<double>({m_nelem, m_nip, s_nne, s_ndim});
317 this->compute_dN_impl();
324 template <
class T,
class R>
325 void interpQuad_vector_impl(
const T& elemvec, R& qvector)
const
332#pragma omp parallel for
333 for (
size_t e = 0; e < m_nelem; ++e) {
335 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
337 for (
size_t q = 0; q < m_nip; ++q) {
339 auto N = xt::adapt(&m_N(q, 0), xt::xshape<s_nne>());
340 auto ui = xt::adapt(&qvector(e, q, 0), xt::xshape<s_ndim>());
342 ui(0) = N(0) * u(0, 0) + N(1) * u(1, 0) + N(2) * u(2, 0) + N(3) * u(3, 0);
343 ui(1) = N(0) * u(0, 1) + N(1) * u(1, 1) + N(2) * u(2, 1) + N(3) * u(3, 1);
348 template <
class T,
class R>
349 void gradN_vector_impl(
const T& elemvec, R& qtensor)
const
354#pragma omp parallel for
355 for (
size_t e = 0; e < m_nelem; ++e) {
357 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
359 for (
size_t q = 0; q < m_nip; ++q) {
361 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
362 auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_ndim, s_ndim>());
365 gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) +
367 gradu(0, 1) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) +
369 gradu(1, 0) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) +
371 gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) +
377 template <
class T,
class R>
378 void gradN_vector_T_impl(
const T& elemvec, R& qtensor)
const
383#pragma omp parallel for
384 for (
size_t e = 0; e < m_nelem; ++e) {
386 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
388 for (
size_t q = 0; q < m_nip; ++q) {
390 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
391 auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_ndim, s_ndim>());
394 gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) +
396 gradu(1, 0) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) +
398 gradu(0, 1) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) +
400 gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) +
406 template <
class T,
class R>
407 void symGradN_vector_impl(
const T& elemvec, R& qtensor)
const
412#pragma omp parallel for
413 for (
size_t e = 0; e < m_nelem; ++e) {
415 auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
417 for (
size_t q = 0; q < m_nip; ++q) {
419 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
420 auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_ndim, s_ndim>());
424 eps(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) +
426 eps(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) +
428 eps(0, 1) = 0.5 * (dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) +
429 dNx(3, 0) * u(3, 1) + dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) +
430 dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0));
431 eps(1, 0) = eps(0, 1);
436 template <
class T,
class R>
437 void int_N_scalar_NT_dV_impl(
const T& qscalar, R& elemmat)
const
444#pragma omp parallel for
445 for (
size_t e = 0; e < m_nelem; ++e) {
447 auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape<s_nne * s_ndim, s_nne * s_ndim>());
449 for (
size_t q = 0; q < m_nip; ++q) {
451 auto N = xt::adapt(&m_N(q, 0), xt::xshape<s_nne>());
452 auto& vol = m_vol(e, q);
453 auto& rho = qscalar(e, q);
456 for (
size_t m = 0; m < s_nne; ++m) {
457 for (
size_t n = 0; n < s_nne; ++n) {
458 M(m * s_ndim + 0, n * s_ndim + 0) += N(m) * rho * N(n) * vol;
459 M(m * s_ndim + 1, n * s_ndim + 1) += N(m) * rho * N(n) * vol;
466 template <
class T,
class R>
467 void int_gradN_dot_tensor2_dV_impl(
const T& qtensor, R& elemvec)
const
474#pragma omp parallel for
475 for (
size_t e = 0; e < m_nelem; ++e) {
477 auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
479 for (
size_t q = 0; q < m_nip; ++q) {
481 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
482 auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_ndim, s_ndim>());
483 auto& vol = m_vol(e, q);
485 for (
size_t m = 0; m < s_nne; ++m) {
486 f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0)) * vol;
487 f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1)) * vol;
493 void compute_dN_impl()
497 array_type::tensor<double, 2> J = xt::empty<double>({2, 2});
498 array_type::tensor<double, 2> Jinv = xt::empty<double>({2, 2});
501 for (
size_t e = 0; e < m_nelem; ++e) {
503 auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape<s_nne, s_ndim>());
505 for (
size_t q = 0; q < m_nip; ++q) {
507 auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape<s_nne, s_ndim>());
508 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
511 J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) +
512 dNxi(3, 0) * x(3, 0);
513 J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) +
514 dNxi(3, 0) * x(3, 1);
515 J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) +
516 dNxi(3, 1) * x(3, 0);
517 J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) +
518 dNxi(3, 1) * x(3, 1);
520 double Jdet = detail::tensor<2>::inv(J, Jinv);
523 for (
size_t m = 0; m < s_nne; ++m) {
524 dNx(m, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1);
525 dNx(m, 1) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1);
528 m_vol(e, q) = m_w(q) * Jdet;
534 constexpr static size_t s_nne = 4;
535 constexpr static size_t s_ndim = 2;
536 constexpr static size_t s_tdim = 2;
540 array_type::tensor<double, 3> m_x;
541 array_type::tensor<double, 1> m_w;
542 array_type::tensor<double, 2> m_xi;
543 array_type::tensor<double, 2> m_N;
544 array_type::tensor<double, 3> m_dNxi;
545 array_type::tensor<double, 4> m_dNx;
546 array_type::tensor<double, 2> m_vol;
Convenience methods for integration point data.
Interpolation and quadrature.
Quadrature(const T &x, const X &xi, const W &w)
Constructor with custom integration.
Quadrature(const T &x)
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_qvector() const -> std::array< size_t, 3 >
Get the shape of a "qvector" (a "qtensor" of rank 1)
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::
array_type::tensor< double, 1 > w()
Integration point weights.
size_t nip()
Number of integration points:
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
size_t nip()
Number of integration points::
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
array_type::tensor< double, 1 > w()
Integration point weights.
size_t nip()
Number of integration points::
array_type::tensor< double, 1 > w()
Integration point weights.
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Toolbox to perform finite element computations.