Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80749052
ElementQuad4Axisymmetric.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 2, 09:56
Size
17 KB
Mime Type
text/x-c
Expires
Wed, Sep 4, 09:56 (2 d)
Engine
blob
Format
Raw Data
Handle
20427574
Attached To
rGOOSEFEM GooseFEM
ElementQuad4Axisymmetric.hpp
View Options
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_ELEMENTQUADAXISYMMETRIC_HPP
#define GOOSEFEM_ELEMENTQUADAXISYMMETRIC_HPP
#include "ElementQuad4Axisymmetric.h"
namespace
GooseFEM
{
namespace
Element
{
namespace
Quad4
{
inline
QuadratureAxisymmetric
::
QuadratureAxisymmetric
(
const
xt
::
xtensor
<
double
,
3
>&
x
)
:
QuadratureAxisymmetric
(
x
,
Gauss
::
xi
(),
Gauss
::
w
())
{
}
inline
QuadratureAxisymmetric
::
QuadratureAxisymmetric
(
const
xt
::
xtensor
<
double
,
3
>&
x
,
const
xt
::
xtensor
<
double
,
2
>&
xi
,
const
xt
::
xtensor
<
double
,
1
>&
w
)
:
m_x
(
x
),
m_w
(
w
),
m_xi
(
xi
)
{
this
->
initQuadratureBase
(
m_x
.
shape
(
0
),
m_w
.
size
());
GOOSEFEM_ASSERT
(
m_x
.
shape
(
1
)
==
m_nne
);
GOOSEFEM_ASSERT
(
m_x
.
shape
(
2
)
==
m_ndim
);
GOOSEFEM_ASSERT
(
m_xi
.
shape
(
0
)
==
m_nip
);
GOOSEFEM_ASSERT
(
m_xi
.
shape
(
1
)
==
m_ndim
);
GOOSEFEM_ASSERT
(
m_w
.
size
()
==
m_nip
);
m_N
=
xt
::
empty
<
double
>
({
m_nip
,
m_nne
});
m_dNxi
=
xt
::
empty
<
double
>
({
m_nip
,
m_nne
,
m_ndim
});
m_B
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nip
,
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
});
m_vol
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nip
});
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
));
}
compute_dN
();
}
inline
xt
::
xtensor
<
double
,
2
>
QuadratureAxisymmetric
::
dV
()
const
{
return
m_vol
;
}
inline
void
QuadratureAxisymmetric
::
update_x
(
const
xt
::
xtensor
<
double
,
3
>&
x
)
{
GOOSEFEM_ASSERT
(
x
.
shape
()
==
m_x
.
shape
());
xt
::
noalias
(
m_x
)
=
x
;
compute_dN
();
}
inline
void
QuadratureAxisymmetric
::
compute_dN
()
{
// most components remain zero, and are not written
m_B
.
fill
(
0.0
);
#pragma omp parallel
{
xt
::
xtensor
<
double
,
2
>
J
=
xt
::
empty
<
double
>
({
2
,
2
});
xt
::
xtensor
<
double
,
2
>
Jinv
=
xt
::
empty
<
double
>
({
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
<
m_nne
,
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
dNxi
=
xt
::
adapt
(
&
m_dNxi
(
q
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_ndim
>
());
auto
B
=
xt
::
adapt
(
&
m_B
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
N
=
xt
::
adapt
(
&
m_N
(
q
,
0
),
xt
::
xshape
<
m_nne
>
());
// 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
=
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
<
m_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
;
}
}
}
}
inline
void
QuadratureAxisymmetric
::
gradN_vector
(
const
xt
::
xtensor
<
double
,
3
>&
elemvec
,
xt
::
xtensor
<
double
,
4
>&
qtensor
)
const
{
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
elemvec
,
{
m_nelem
,
m_nne
,
m_ndim
}));
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
qtensor
,
{
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
}));
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
<
m_nne
,
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
B
=
xt
::
adapt
(
&
m_B
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
gradu
=
xt
::
adapt
(
&
qtensor
(
e
,
q
,
0
,
0
),
xt
::
xshape
<
m_tdim
,
m_tdim
>
());
// 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
);
}
}
}
inline
void
QuadratureAxisymmetric
::
gradN_vector_T
(
const
xt
::
xtensor
<
double
,
3
>&
elemvec
,
xt
::
xtensor
<
double
,
4
>&
qtensor
)
const
{
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
elemvec
,
{
m_nelem
,
m_nne
,
m_ndim
}));
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
qtensor
,
{
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
}));
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
<
m_nne
,
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
B
=
xt
::
adapt
(
&
m_B
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
gradu
=
xt
::
adapt
(
&
qtensor
(
e
,
q
,
0
,
0
),
xt
::
xshape
<
m_tdim
,
m_tdim
>
());
// 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
);
}
}
}
inline
void
QuadratureAxisymmetric
::
symGradN_vector
(
const
xt
::
xtensor
<
double
,
3
>&
elemvec
,
xt
::
xtensor
<
double
,
4
>&
qtensor
)
const
{
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
elemvec
,
{
m_nelem
,
m_nne
,
m_ndim
}));
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
qtensor
,
{
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
}));
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
<
m_nne
,
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
B
=
xt
::
adapt
(
&
m_B
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
eps
=
xt
::
adapt
(
&
qtensor
(
e
,
q
,
0
,
0
),
xt
::
xshape
<
m_tdim
,
m_tdim
>
());
// 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
);
}
}
}
inline
void
QuadratureAxisymmetric
::
int_N_scalar_NT_dV
(
const
xt
::
xtensor
<
double
,
2
>&
qscalar
,
xt
::
xtensor
<
double
,
3
>&
elemmat
)
const
{
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
qscalar
,
{
m_nelem
,
m_nip
}));
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
elemmat
,
{
m_nelem
,
m_nne
*
m_ndim
,
m_nne
*
m_ndim
}));
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
<
m_nne
*
m_ndim
,
m_nne
*
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
N
=
xt
::
adapt
(
&
m_N
(
q
,
0
),
xt
::
xshape
<
m_nne
>
());
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
<
m_nne
;
++
m
)
{
for
(
size_t
n
=
0
;
n
<
m_nne
;
++
n
)
{
M
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
0
)
+=
N
(
m
)
*
rho
*
N
(
n
)
*
vol
;
M
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
N
(
m
)
*
rho
*
N
(
n
)
*
vol
;
}
}
}
}
}
inline
void
QuadratureAxisymmetric
::
int_gradN_dot_tensor2_dV
(
const
xt
::
xtensor
<
double
,
4
>&
qtensor
,
xt
::
xtensor
<
double
,
3
>&
elemvec
)
const
{
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
qtensor
,
{
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
}));
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
elemvec
,
{
m_nelem
,
m_nne
,
m_ndim
}));
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
<
m_nne
,
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
B
=
xt
::
adapt
(
&
m_B
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
sig
=
xt
::
adapt
(
&
qtensor
(
e
,
q
,
0
,
0
),
xt
::
xshape
<
m_tdim
,
m_tdim
>
());
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
<
m_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
));
}
}
}
}
inline
void
QuadratureAxisymmetric
::
int_gradN_dot_tensor4_dot_gradNT_dV
(
const
xt
::
xtensor
<
double
,
6
>&
qtensor
,
xt
::
xtensor
<
double
,
3
>&
elemmat
)
const
{
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
qtensor
,
{
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
,
m_tdim
,
m_tdim
}));
GOOSEFEM_ASSERT
(
xt
::
has_shape
(
elemmat
,
{
m_nelem
,
m_nne
*
m_ndim
,
m_nne
*
m_ndim
}));
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
<
m_nne
*
m_ndim
,
m_nne
*
m_ndim
>
());
for
(
size_t
q
=
0
;
q
<
m_nip
;
++
q
)
{
auto
B
=
xt
::
adapt
(
&
m_B
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_nne
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
C
=
xt
::
adapt
(
&
qtensor
(
e
,
q
,
0
,
0
,
0
,
0
),
xt
::
xshape
<
m_tdim
,
m_tdim
,
m_tdim
,
m_tdim
>
());
auto
&
vol
=
m_vol
(
e
,
q
);
// K(m*m_ndim+perm(c), n*m_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
<
m_nne
;
++
m
)
{
for
(
size_t
n
=
0
;
n
<
m_nne
;
++
n
)
{
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
0
,
0
,
0
)
*
C
(
0
,
0
,
0
,
0
)
*
B
(
n
,
0
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
0
,
0
,
0
)
*
C
(
0
,
0
,
1
,
1
)
*
B
(
n
,
1
,
1
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
0
,
0
,
0
)
*
C
(
0
,
0
,
2
,
2
)
*
B
(
n
,
2
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
0
,
0
,
0
)
*
C
(
0
,
0
,
2
,
0
)
*
B
(
n
,
0
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
0
,
0
,
0
)
*
C
(
0
,
0
,
0
,
2
)
*
B
(
n
,
2
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
1
,
1
,
0
)
*
C
(
1
,
1
,
0
,
0
)
*
B
(
n
,
0
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
1
,
1
,
0
)
*
C
(
1
,
1
,
1
,
1
)
*
B
(
n
,
1
,
1
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
1
,
1
,
0
)
*
C
(
1
,
1
,
2
,
2
)
*
B
(
n
,
2
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
1
,
1
,
0
)
*
C
(
1
,
1
,
2
,
0
)
*
B
(
n
,
0
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
1
,
1
,
0
)
*
C
(
1
,
1
,
0
,
2
)
*
B
(
n
,
2
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
2
,
2
,
2
)
*
C
(
2
,
2
,
0
,
0
)
*
B
(
n
,
0
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
2
,
2
,
2
)
*
C
(
2
,
2
,
1
,
1
)
*
B
(
n
,
1
,
1
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
2
,
2
,
2
)
*
C
(
2
,
2
,
2
,
2
)
*
B
(
n
,
2
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
2
,
2
,
2
)
*
C
(
2
,
2
,
2
,
0
)
*
B
(
n
,
0
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
2
,
2
,
2
)
*
C
(
2
,
2
,
0
,
2
)
*
B
(
n
,
2
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
0
,
2
,
2
)
*
C
(
0
,
2
,
0
,
0
)
*
B
(
n
,
0
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
0
,
2
,
2
)
*
C
(
0
,
2
,
1
,
1
)
*
B
(
n
,
1
,
1
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
0
,
2
,
2
)
*
C
(
0
,
2
,
2
,
2
)
*
B
(
n
,
2
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
0
,
2
,
2
)
*
C
(
0
,
2
,
2
,
0
)
*
B
(
n
,
0
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
0
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
0
,
2
,
2
)
*
C
(
0
,
2
,
0
,
2
)
*
B
(
n
,
2
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
2
,
0
,
0
)
*
C
(
2
,
0
,
0
,
0
)
*
B
(
n
,
0
,
0
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
2
,
0
,
0
)
*
C
(
2
,
0
,
1
,
1
)
*
B
(
n
,
1
,
1
,
0
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
2
,
0
,
0
)
*
C
(
2
,
0
,
2
,
2
)
*
B
(
n
,
2
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
0
)
+=
B
(
m
,
2
,
0
,
0
)
*
C
(
2
,
0
,
2
,
0
)
*
B
(
n
,
0
,
2
,
2
)
*
vol
;
K
(
m
*
m_ndim
+
1
,
n
*
m_ndim
+
1
)
+=
B
(
m
,
2
,
0
,
0
)
*
C
(
2
,
0
,
0
,
2
)
*
B
(
n
,
2
,
0
,
0
)
*
vol
;
}
}
}
}
}
inline
xt
::
xtensor
<
double
,
4
>
QuadratureAxisymmetric
::
GradN_vector
(
const
xt
::
xtensor
<
double
,
3
>&
elemvec
)
const
{
xt
::
xtensor
<
double
,
4
>
qtensor
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
});
this
->
gradN_vector
(
elemvec
,
qtensor
);
return
qtensor
;
}
inline
xt
::
xtensor
<
double
,
4
>
QuadratureAxisymmetric
::
GradN_vector_T
(
const
xt
::
xtensor
<
double
,
3
>&
elemvec
)
const
{
xt
::
xtensor
<
double
,
4
>
qtensor
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
});
this
->
gradN_vector_T
(
elemvec
,
qtensor
);
return
qtensor
;
}
inline
xt
::
xtensor
<
double
,
4
>
QuadratureAxisymmetric
::
SymGradN_vector
(
const
xt
::
xtensor
<
double
,
3
>&
elemvec
)
const
{
xt
::
xtensor
<
double
,
4
>
qtensor
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nip
,
m_tdim
,
m_tdim
});
this
->
symGradN_vector
(
elemvec
,
qtensor
);
return
qtensor
;
}
inline
xt
::
xtensor
<
double
,
3
>
QuadratureAxisymmetric
::
Int_N_scalar_NT_dV
(
const
xt
::
xtensor
<
double
,
2
>&
qscalar
)
const
{
xt
::
xtensor
<
double
,
3
>
elemmat
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nne
*
m_ndim
,
m_nne
*
m_ndim
});
this
->
int_N_scalar_NT_dV
(
qscalar
,
elemmat
);
return
elemmat
;
}
inline
xt
::
xtensor
<
double
,
3
>
QuadratureAxisymmetric
::
Int_gradN_dot_tensor2_dV
(
const
xt
::
xtensor
<
double
,
4
>&
qtensor
)
const
{
xt
::
xtensor
<
double
,
3
>
elemvec
=
xt
::
empty
<
double
>
({
m_nelem
,
m_nne
,
m_ndim
});
this
->
int_gradN_dot_tensor2_dV
(
qtensor
,
elemvec
);
return
elemvec
;
}
inline
xt
::
xtensor
<
double
,
3
>
QuadratureAxisymmetric
::
Int_gradN_dot_tensor4_dot_gradNT_dV
(
const
xt
::
xtensor
<
double
,
6
>&
qtensor
)
const
{
xt
::
xtensor
<
double
,
3
>
elemmat
=
xt
::
empty
<
double
>
({
m_nelem
,
m_ndim
*
m_nne
,
m_ndim
*
m_nne
});
this
->
int_gradN_dot_tensor4_dot_gradNT_dV
(
qtensor
,
elemmat
);
return
elemmat
;
}
}
// namespace Quad4
}
// namespace Element
}
// namespace GooseFEM
#endif
Event Timeline
Log In to Comment