Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86661859
DynamicsDiagonalQuad4.h
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, Oct 7, 20:57
Size
11 KB
Mime Type
text/x-c++
Expires
Wed, Oct 9, 20:57 (2 d)
Engine
blob
Format
Raw Data
Handle
21465411
Attached To
rGOOSEFEM GooseFEM
DynamicsDiagonalQuad4.h
View Options
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_QUAD4_H
#define GOOSEFEM_DYNAMICS_DIAGONAL_QUAD4_H
#include "Macros.h"
#include <cppmat/cppmat.h>
// -------------------------------------------------------------------------------------------------
namespace
GooseFEM
{
namespace
Dynamics
{
namespace
Diagonal
{
// -------------------------------------------------------------------------------------------------
using
vec
=
cppmat
::
cartesian2d
::
vector
<
double
>
;
using
T2
=
cppmat
::
cartesian2d
::
tensor2
<
double
>
;
using
T2s
=
cppmat
::
cartesian2d
::
tensor2s
<
double
>
;
using
T2d
=
cppmat
::
cartesian2d
::
tensor2d
<
double
>
;
// ========================== N.B. most loops are unrolled for efficiency ==========================
template
<
class
QuadraturePoint
>
class
Quad4
{
public:
// variables
// ---------
// arrays / matrices
cppmat
::
tiny
::
matrix2
<
double
,
4
,
2
>
xe
,
ue
,
ve
,
xi
,
xi_n
,
dNdxi
,
dNdx
;
cppmat
::
tiny
::
vector
<
double
,
4
>
w
,
w_n
;
cppmat
::
tiny
::
matrix2
<
double
,
8
,
8
>
M
;
cppmat
::
tiny
::
vector
<
double
,
8
>
fu
,
fv
;
// tensors
cppmat
::
cartesian2d
::
tensor2
<
double
>
J
,
Jinv
,
gradu
,
gradv
;
// scalars
double
Jdet
,
V
;
// sizes (nodes per element, dimensions, quadrature points)
size_t
nne
=
4
,
ndim
=
2
,
nk
=
4
;
// quadrature point
std
::
shared_ptr
<
QuadraturePoint
>
quad
;
// constructor
// -----------
Quad4
(
std
::
shared_ptr
<
QuadraturePoint
>
quad
);
// functions
// ---------
void
computeM
(
size_t
elem
);
// mass matrix <- quad->density
void
computeFu
(
size_t
elem
);
// displacement dependent forces <- quad->stressStrain
void
computeFv
(
size_t
elem
);
// displacement dependent forces <- quad->stressStrainRate
void
post
(
size_t
elem
);
// post-process <- quad->stressStrain(Rate)
};
// =================================================================================================
template
<
class
QuadraturePoint
>
Quad4
<
QuadraturePoint
>::
Quad4
(
std
::
shared_ptr
<
QuadraturePoint
>
_quad
)
{
// quadrature point routines
quad
=
_quad
;
// integration point coordinates/weights: normal Gauss integration
xi
(
0
,
0
)
=
-
1.
/
std
::
sqrt
(
3.
);
xi
(
0
,
1
)
=
-
1.
/
std
::
sqrt
(
3.
);
w
(
0
)
=
1.
;
xi
(
1
,
0
)
=
+
1.
/
std
::
sqrt
(
3.
);
xi
(
1
,
1
)
=
-
1.
/
std
::
sqrt
(
3.
);
w
(
1
)
=
1.
;
xi
(
2
,
0
)
=
+
1.
/
std
::
sqrt
(
3.
);
xi
(
2
,
1
)
=
+
1.
/
std
::
sqrt
(
3.
);
w
(
2
)
=
1.
;
xi
(
3
,
0
)
=
-
1.
/
std
::
sqrt
(
3.
);
xi
(
3
,
1
)
=
+
1.
/
std
::
sqrt
(
3.
);
w
(
3
)
=
1.
;
// integration point coordinates/weights: integration at the nodes
xi_n
(
0
,
0
)
=
-
1.
;
xi_n
(
0
,
1
)
=
-
1.
;
w_n
(
0
)
=
1.
;
xi_n
(
1
,
0
)
=
+
1.
;
xi_n
(
1
,
1
)
=
-
1.
;
w_n
(
1
)
=
1.
;
xi_n
(
2
,
0
)
=
+
1.
;
xi_n
(
2
,
1
)
=
+
1.
;
w_n
(
2
)
=
1.
;
xi_n
(
3
,
0
)
=
-
1.
;
xi_n
(
3
,
1
)
=
+
1.
;
w_n
(
3
)
=
1.
;
}
// =================================================================================================
template
<
class
QuadraturePoint
>
void
Quad4
<
QuadraturePoint
>::
computeM
(
size_t
elem
)
{
// zero-initialize element mass matrix
M
.
zeros
();
// loop over integration points (coincide with the nodes to get a diagonal mass matrix)
for
(
size_t
k
=
0
;
k
<
nne
;
++
k
)
{
// - shape function gradients (local coordinates)
dNdxi
(
0
,
0
)
=
-
.25
*
(
1.
-
xi_n
(
k
,
1
));
dNdxi
(
0
,
1
)
=
-
.25
*
(
1.
-
xi_n
(
k
,
0
));
dNdxi
(
1
,
0
)
=
+
.25
*
(
1.
-
xi_n
(
k
,
1
));
dNdxi
(
1
,
1
)
=
-
.25
*
(
1.
+
xi_n
(
k
,
0
));
dNdxi
(
2
,
0
)
=
+
.25
*
(
1.
+
xi_n
(
k
,
1
));
dNdxi
(
2
,
1
)
=
+
.25
*
(
1.
+
xi_n
(
k
,
0
));
dNdxi
(
3
,
0
)
=
-
.25
*
(
1.
+
xi_n
(
k
,
1
));
dNdxi
(
3
,
1
)
=
+
.25
*
(
1.
-
xi_n
(
k
,
0
));
// - Jacobian
// J(i,j) += dNdxi(m,i) * xe(m,j)
J
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
J
(
0
,
0
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
0
);
J
(
0
,
1
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
1
);
J
(
1
,
0
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
0
);
J
(
1
,
1
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
1
);
}
// - determinant and inverse of the Jacobian
Jdet
=
J
.
det
();
Jinv
=
J
.
inv
();
// - integration point volume (== volume associated with the node, in this element)
V
=
w_n
(
k
)
*
Jdet
;
// - assemble to element mass matrix (use the delta properties of the shape functions)
// M(m+i,n+i) = N(m) * rho * V * N(n);
M
(
k
*
2
,
k
*
2
)
=
quad
->
density
(
elem
,
k
,
V
)
*
V
;
M
(
k
*
2
+
1
,
k
*
2
+
1
)
=
quad
->
density
(
elem
,
k
,
V
)
*
V
;
}
}
// =================================================================================================
template
<
class
QuadraturePoint
>
void
Quad4
<
QuadraturePoint
>::
computeFu
(
size_t
elem
)
{
// zero-initialize element force
fu
.
zeros
();
// loop over integration points
for
(
size_t
k
=
0
;
k
<
nk
;
++
k
)
{
// - shape function gradients (local coordinates)
dNdxi
(
0
,
0
)
=
-
.25
*
(
1.
-
xi
(
k
,
1
));
dNdxi
(
0
,
1
)
=
-
.25
*
(
1.
-
xi
(
k
,
0
));
dNdxi
(
1
,
0
)
=
+
.25
*
(
1.
-
xi
(
k
,
1
));
dNdxi
(
1
,
1
)
=
-
.25
*
(
1.
+
xi
(
k
,
0
));
dNdxi
(
2
,
0
)
=
+
.25
*
(
1.
+
xi
(
k
,
1
));
dNdxi
(
2
,
1
)
=
+
.25
*
(
1.
+
xi
(
k
,
0
));
dNdxi
(
3
,
0
)
=
-
.25
*
(
1.
+
xi
(
k
,
1
));
dNdxi
(
3
,
1
)
=
+
.25
*
(
1.
-
xi
(
k
,
0
));
// - Jacobian
// J(i,j) += dNdxi(m,i) * xe(m,j)
J
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
J
(
0
,
0
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
0
);
J
(
0
,
1
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
1
);
J
(
1
,
0
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
0
);
J
(
1
,
1
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
1
);
}
// - determinant and inverse of the Jacobian
Jdet
=
J
.
det
();
Jinv
=
J
.
inv
();
// - integration point volume
V
=
w
(
k
)
*
Jdet
;
// - shape function gradients (global coordinates)
// dNdx(m,i) += Jinv(i,j) * dNdxi(m,j)
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
dNdx
(
m
,
0
)
=
Jinv
(
0
,
0
)
*
dNdxi
(
m
,
0
)
+
Jinv
(
0
,
1
)
*
dNdxi
(
m
,
1
);
dNdx
(
m
,
1
)
=
Jinv
(
1
,
0
)
*
dNdxi
(
m
,
0
)
+
Jinv
(
1
,
1
)
*
dNdxi
(
m
,
1
);
}
// - displacement gradient
// gradu(i,j) += dNdx(m,i) * ue(m,j)
gradu
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
gradu
(
0
,
0
)
+=
dNdx
(
m
,
0
)
*
ue
(
m
,
0
);
gradu
(
0
,
1
)
+=
dNdx
(
m
,
0
)
*
ue
(
m
,
1
);
gradu
(
1
,
0
)
+=
dNdx
(
m
,
1
)
*
ue
(
m
,
0
);
gradu
(
1
,
1
)
+=
dNdx
(
m
,
1
)
*
ue
(
m
,
1
);
}
// - strain tensor (symmetric part of "gradu")
// quad->eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) )
quad
->
eps
(
0
,
0
)
=
gradu
(
0
,
0
);
quad
->
eps
(
0
,
1
)
=
.5
*
(
gradu
(
0
,
1
)
+
gradu
(
1
,
0
)
);
quad
->
eps
(
1
,
0
)
=
quad
->
eps
(
0
,
1
);
quad
->
eps
(
1
,
1
)
=
gradu
(
1
,
1
);
// - constitutive response
quad
->
stressStrain
(
elem
,
k
,
V
);
// - assemble to element force
// fu(m*ndim+j) += dNdx(m,i) * quad->sig(i,j) * V;
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
fu
(
m
*
ndim
+
0
)
+=
dNdx
(
m
,
0
)
*
quad
->
sig
(
0
,
0
)
*
V
+
dNdx
(
m
,
1
)
*
quad
->
sig
(
1
,
0
)
*
V
;
fu
(
m
*
ndim
+
1
)
+=
dNdx
(
m
,
0
)
*
quad
->
sig
(
0
,
1
)
*
V
+
dNdx
(
m
,
1
)
*
quad
->
sig
(
1
,
1
)
*
V
;
}
}
}
// =================================================================================================
template
<
class
QuadraturePoint
>
void
Quad4
<
QuadraturePoint
>::
computeFv
(
size_t
elem
)
{
// zero-initialize element force
fv
.
zeros
();
// loop over integration points
for
(
size_t
k
=
0
;
k
<
nk
;
++
k
)
{
// - shape function gradients (local coordinates)
dNdxi
(
0
,
0
)
=
-
.25
*
(
1.
-
xi
(
k
,
1
));
dNdxi
(
0
,
1
)
=
-
.25
*
(
1.
-
xi
(
k
,
0
));
dNdxi
(
1
,
0
)
=
+
.25
*
(
1.
-
xi
(
k
,
1
));
dNdxi
(
1
,
1
)
=
-
.25
*
(
1.
+
xi
(
k
,
0
));
dNdxi
(
2
,
0
)
=
+
.25
*
(
1.
+
xi
(
k
,
1
));
dNdxi
(
2
,
1
)
=
+
.25
*
(
1.
+
xi
(
k
,
0
));
dNdxi
(
3
,
0
)
=
-
.25
*
(
1.
+
xi
(
k
,
1
));
dNdxi
(
3
,
1
)
=
+
.25
*
(
1.
-
xi
(
k
,
0
));
// - Jacobian
// J(i,j) += dNdxi(m,i) * xe(m,j)
J
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
J
(
0
,
0
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
0
);
J
(
0
,
1
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
1
);
J
(
1
,
0
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
0
);
J
(
1
,
1
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
1
);
}
// - determinant and inverse of the Jacobian
Jdet
=
J
.
det
();
Jinv
=
J
.
inv
();
// - integration point volume
V
=
w
(
k
)
*
Jdet
;
// - shape function gradients (global coordinates)
// dNdx(m,i) += Jinv(i,j) * dNdxi(m,j)
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
dNdx
(
m
,
0
)
=
Jinv
(
0
,
0
)
*
dNdxi
(
m
,
0
)
+
Jinv
(
0
,
1
)
*
dNdxi
(
m
,
1
);
dNdx
(
m
,
1
)
=
Jinv
(
1
,
0
)
*
dNdxi
(
m
,
0
)
+
Jinv
(
1
,
1
)
*
dNdxi
(
m
,
1
);
}
// - velocity gradient
// gradv(i,j) += dNdx(m,i) * ve(m,j)
gradv
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
gradv
(
0
,
0
)
+=
dNdx
(
m
,
0
)
*
ve
(
m
,
0
);
gradv
(
0
,
1
)
+=
dNdx
(
m
,
0
)
*
ve
(
m
,
1
);
gradv
(
1
,
0
)
+=
dNdx
(
m
,
1
)
*
ve
(
m
,
0
);
gradv
(
1
,
1
)
+=
dNdx
(
m
,
1
)
*
ve
(
m
,
1
);
}
// - strain-rate tensor (symmetric part of "gradv")
// quad->epsdot(i,j) = .5 * ( gradv(i,j) + gradv(j,i) )
quad
->
epsdot
(
0
,
0
)
=
gradv
(
0
,
0
);
quad
->
epsdot
(
0
,
1
)
=
.5
*
(
gradv
(
0
,
1
)
+
gradv
(
1
,
0
)
);
quad
->
epsdot
(
1
,
0
)
=
quad
->
epsdot
(
0
,
1
);
quad
->
epsdot
(
1
,
1
)
=
gradv
(
1
,
1
);
// - constitutive response
quad
->
stressStrainRate
(
elem
,
k
,
V
);
// - assemble to element force
// fv(m*ndim+j) += dNdx(m,i) * quad->sig(i,j) * V;
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
fv
(
m
*
2
+
0
)
+=
dNdx
(
m
,
0
)
*
quad
->
sig
(
0
,
0
)
*
V
+
dNdx
(
m
,
1
)
*
quad
->
sig
(
1
,
0
)
*
V
;
fv
(
m
*
2
+
1
)
+=
dNdx
(
m
,
0
)
*
quad
->
sig
(
0
,
1
)
*
V
+
dNdx
(
m
,
1
)
*
quad
->
sig
(
1
,
1
)
*
V
;
}
}
}
// =================================================================================================
template
<
class
QuadraturePoint
>
void
Quad4
<
QuadraturePoint
>::
post
(
size_t
elem
)
{
// loop over integration points
for
(
size_t
k
=
0
;
k
<
nk
;
++
k
)
{
// - shape function gradients (local coordinates)
dNdxi
(
0
,
0
)
=
-
.25
*
(
1.
-
xi
(
k
,
1
));
dNdxi
(
0
,
1
)
=
-
.25
*
(
1.
-
xi
(
k
,
0
));
dNdxi
(
1
,
0
)
=
+
.25
*
(
1.
-
xi
(
k
,
1
));
dNdxi
(
1
,
1
)
=
-
.25
*
(
1.
+
xi
(
k
,
0
));
dNdxi
(
2
,
0
)
=
+
.25
*
(
1.
+
xi
(
k
,
1
));
dNdxi
(
2
,
1
)
=
+
.25
*
(
1.
+
xi
(
k
,
0
));
dNdxi
(
3
,
0
)
=
-
.25
*
(
1.
+
xi
(
k
,
1
));
dNdxi
(
3
,
1
)
=
+
.25
*
(
1.
-
xi
(
k
,
0
));
// - Jacobian
// J(i,j) += dNdxi(m,i) * xe(m,j)
J
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
J
(
0
,
0
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
0
);
J
(
0
,
1
)
+=
dNdxi
(
m
,
0
)
*
xe
(
m
,
1
);
J
(
1
,
0
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
0
);
J
(
1
,
1
)
+=
dNdxi
(
m
,
1
)
*
xe
(
m
,
1
);
}
// - determinant and inverse of the Jacobian
Jdet
=
J
.
det
();
Jinv
=
J
.
inv
();
// - integration point volume
V
=
w
(
k
)
*
Jdet
;
// - shape function gradients (global coordinates)
// dNdx(m,i) += Jinv(i,j) * dNdxi(m,j)
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
dNdx
(
m
,
0
)
=
Jinv
(
0
,
0
)
*
dNdxi
(
m
,
0
)
+
Jinv
(
0
,
1
)
*
dNdxi
(
m
,
1
);
dNdx
(
m
,
1
)
=
Jinv
(
1
,
0
)
*
dNdxi
(
m
,
0
)
+
Jinv
(
1
,
1
)
*
dNdxi
(
m
,
1
);
}
// - displacement gradient
// gradu(i,j) += dNdx(m,i) * ue(m,j)
gradu
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
gradu
(
0
,
0
)
+=
dNdx
(
m
,
0
)
*
ue
(
m
,
0
);
gradu
(
0
,
1
)
+=
dNdx
(
m
,
0
)
*
ue
(
m
,
1
);
gradu
(
1
,
0
)
+=
dNdx
(
m
,
1
)
*
ue
(
m
,
0
);
gradu
(
1
,
1
)
+=
dNdx
(
m
,
1
)
*
ue
(
m
,
1
);
}
// - strain tensor (symmetric part of "gradu")
// quad->eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) )
quad
->
eps
(
0
,
0
)
=
gradu
(
0
,
0
);
quad
->
eps
(
0
,
1
)
=
.5
*
(
gradu
(
0
,
1
)
+
gradu
(
1
,
0
)
);
quad
->
eps
(
1
,
0
)
=
quad
->
eps
(
0
,
1
);
quad
->
eps
(
1
,
1
)
=
gradu
(
1
,
1
);
// - velocity gradient
// gradv(i,j) += dNdx(m,i) * ve(m,j)
gradv
.
zeros
();
for
(
size_t
m
=
0
;
m
<
nne
;
++
m
)
{
gradv
(
0
,
0
)
+=
dNdx
(
m
,
0
)
*
ve
(
m
,
0
);
gradv
(
0
,
1
)
+=
dNdx
(
m
,
0
)
*
ve
(
m
,
1
);
gradv
(
1
,
0
)
+=
dNdx
(
m
,
1
)
*
ve
(
m
,
0
);
gradv
(
1
,
1
)
+=
dNdx
(
m
,
1
)
*
ve
(
m
,
1
);
}
// - strain-rate tensor (symmetric part of "gradv")
// quad->epsdot(i,j) = .5 * ( gradv(i,j) + gradv(j,i) )
quad
->
epsdot
(
0
,
0
)
=
gradv
(
0
,
0
);
quad
->
epsdot
(
0
,
1
)
=
.5
*
(
gradv
(
0
,
1
)
+
gradv
(
1
,
0
)
);
quad
->
epsdot
(
1
,
0
)
=
quad
->
epsdot
(
0
,
1
);
quad
->
epsdot
(
1
,
1
)
=
gradv
(
1
,
1
);
// - constitutive response
quad
->
stressStrainPost
(
elem
,
k
,
V
);
quad
->
stressStrainRatePost
(
elem
,
k
,
V
);
}
}
// =================================================================================================
}}}
// namespace GooseFEM::Dynamics::Diagonal
#endif
Event Timeline
Log In to Comment