Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81519134
FE_Element.cpp
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
Fri, Sep 6, 12:34
Size
17 KB
Mime Type
text/x-c
Expires
Sun, Sep 8, 12:34 (2 d)
Engine
blob
Format
Raw Data
Handle
20513524
Attached To
rLAMMPS lammps
FE_Element.cpp
View Options
// ATC header files
#include "ATC_Error.h"
#include "FE_Element.h"
#include "FE_Mesh.h"
// Other headers
#include "math.h"
namespace
ATC
{
// -------------------------------------------------------------
// -------------------------------------------------------------
// class FE_Element
// -------------------------------------------------------------
// -------------------------------------------------------------
FE_Element
::
FE_Element
(
FE_Mesh
*
feMesh
,
int
nSD
,
int
numEltNodes
,
int
numIPs
,
int
numFaces
,
int
numFaceNodes
,
int
numFaceIPs
)
:
feMesh_
(
feMesh
),
nSD_
(
nSD
),
numEltNodes_
(
numEltNodes
),
numIPs_
(
numIPs
),
numFaces_
(
numFaces
),
numFaceNodes_
(
numFaceNodes
),
numFaceIPs_
(
numFaceIPs
),
localCoords_
(
nSD
,
numEltNodes
),
ipCoords_
(
nSD
,
numIPs
),
N_
(
numIPs
,
numEltNodes
),
ipWeights_
(
numIPs
),
localFaceConn_
(
numFaces
,
numFaceNodes
),
ipFace2DCoords_
(
nSD
-
1
,
numFaceIPs
),
ipFaceCoords_
(
numFaces
),
// localFaceCoords_(nSD-1,numFaceNodes),
NFace_
(
numFaces
),
ipFaceWeights_
(
numFaceIPs
)
{
for
(
int
f
=
0
;
f
<
numFaces
;
f
++
)
{
ipFaceCoords_
[
f
].
reset
(
nSD
,
numFaceIPs
);
NFace_
[
f
].
reset
(
numFaceIPs
,
numEltNodes
);
}
}
FE_Element
::~
FE_Element
()
{
// Nothing else to do
}
// -------------------------------------------------------------
// -------------------------------------------------------------
// class FE_ElementHex
// -------------------------------------------------------------
// -------------------------------------------------------------
FE_ElementHex
::
FE_ElementHex
(
FE_Mesh
*
feMesh
)
:
FE_Element
(
feMesh
,
3
,
8
,
8
,
6
,
4
,
4
)
{
// Matrix of local nodal coordinates
localCoords_
(
0
,
0
)
=
-
1
;
localCoords_
(
1
,
0
)
=
-
1
;
localCoords_
(
2
,
0
)
=
-
1
;
localCoords_
(
0
,
1
)
=
+
1
;
localCoords_
(
1
,
1
)
=
-
1
;
localCoords_
(
2
,
1
)
=
-
1
;
localCoords_
(
0
,
2
)
=
+
1
;
localCoords_
(
1
,
2
)
=
+
1
;
localCoords_
(
2
,
2
)
=
-
1
;
localCoords_
(
0
,
3
)
=
-
1
;
localCoords_
(
1
,
3
)
=
+
1
;
localCoords_
(
2
,
3
)
=
-
1
;
localCoords_
(
0
,
4
)
=
-
1
;
localCoords_
(
1
,
4
)
=
-
1
;
localCoords_
(
2
,
4
)
=
+
1
;
localCoords_
(
0
,
5
)
=
+
1
;
localCoords_
(
1
,
5
)
=
-
1
;
localCoords_
(
2
,
5
)
=
+
1
;
localCoords_
(
0
,
6
)
=
+
1
;
localCoords_
(
1
,
6
)
=
+
1
;
localCoords_
(
2
,
6
)
=
+
1
;
localCoords_
(
0
,
7
)
=
-
1
;
localCoords_
(
1
,
7
)
=
+
1
;
localCoords_
(
2
,
7
)
=
+
1
;
// 3 --- 2
// /| /|
// / 0 --/ 1 y
// 7 --- 6 / |
// | |/ |
// 4 --- 5 ---> x
// /
// /
// z
// Matrix of local face connectivity
// -x
localFaceConn_
(
0
,
0
)
=
0
;
localFaceConn_
(
0
,
1
)
=
4
;
localFaceConn_
(
0
,
2
)
=
7
;
localFaceConn_
(
0
,
3
)
=
3
;
// +x
localFaceConn_
(
1
,
0
)
=
1
;
localFaceConn_
(
1
,
1
)
=
2
;
localFaceConn_
(
1
,
2
)
=
6
;
localFaceConn_
(
1
,
3
)
=
5
;
// -y
localFaceConn_
(
2
,
0
)
=
0
;
localFaceConn_
(
2
,
1
)
=
1
;
localFaceConn_
(
2
,
2
)
=
5
;
localFaceConn_
(
2
,
3
)
=
4
;
// +y
localFaceConn_
(
3
,
0
)
=
2
;
localFaceConn_
(
3
,
1
)
=
3
;
localFaceConn_
(
3
,
2
)
=
7
;
localFaceConn_
(
3
,
3
)
=
6
;
// -z
localFaceConn_
(
4
,
0
)
=
0
;
localFaceConn_
(
4
,
1
)
=
3
;
localFaceConn_
(
4
,
2
)
=
2
;
localFaceConn_
(
4
,
3
)
=
1
;
// +z
localFaceConn_
(
5
,
0
)
=
4
;
localFaceConn_
(
5
,
1
)
=
5
;
localFaceConn_
(
5
,
2
)
=
6
;
localFaceConn_
(
5
,
3
)
=
7
;
/*
localFaceCoords_(0,0) = -1; localFaceCoords_(1,0) = -1;
localFaceCoords_(0,1) = +1; localFaceCoords_(1,1) = -1;
localFaceCoords_(0,2) = +1; localFaceCoords_(1,2) = +1;
localFaceCoords_(0,3) = -1; localFaceCoords_(1,3) = +1;
*/
set_quadrature
(
GAUSSIAN_QUADRATURE
);
}
//-----------------------------------------------------------------
void
FE_ElementHex
::
set_quadrature
(
int
quadrature
)
{
double
a
=
1.
/
sqrt
(
3.
);
if
(
quadrature
==
NODAL_QUADRATURE
)
{
a
=
1.0
;
}
// Matrix of integration point locations (Gaussian) & follows local conn
ipCoords_
(
0
,
0
)
=
-
a
;
ipCoords_
(
1
,
0
)
=
-
a
;
ipCoords_
(
2
,
0
)
=
-
a
;
ipCoords_
(
0
,
1
)
=
+
a
;
ipCoords_
(
1
,
1
)
=
-
a
;
ipCoords_
(
2
,
1
)
=
-
a
;
ipCoords_
(
0
,
2
)
=
+
a
;
ipCoords_
(
1
,
2
)
=
+
a
;
ipCoords_
(
2
,
2
)
=
-
a
;
ipCoords_
(
0
,
3
)
=
-
a
;
ipCoords_
(
1
,
3
)
=
+
a
;
ipCoords_
(
2
,
3
)
=
-
a
;
ipCoords_
(
0
,
4
)
=
-
a
;
ipCoords_
(
1
,
4
)
=
-
a
;
ipCoords_
(
2
,
4
)
=
+
a
;
ipCoords_
(
0
,
5
)
=
+
a
;
ipCoords_
(
1
,
5
)
=
-
a
;
ipCoords_
(
2
,
5
)
=
+
a
;
ipCoords_
(
0
,
6
)
=
+
a
;
ipCoords_
(
1
,
6
)
=
+
a
;
ipCoords_
(
2
,
6
)
=
+
a
;
ipCoords_
(
0
,
7
)
=
-
a
;
ipCoords_
(
1
,
7
)
=
+
a
;
ipCoords_
(
2
,
7
)
=
+
a
;
// Compute shape functions at ip's
for
(
int
ip
=
0
;
ip
<
numIPs_
;
ip
++
)
{
double
r
=
ipCoords_
(
0
,
ip
);
double
s
=
ipCoords_
(
1
,
ip
);
double
t
=
ipCoords_
(
2
,
ip
);
for
(
int
iNode
=
0
;
iNode
<
numEltNodes_
;
iNode
++
)
{
double
rI
=
localCoords_
(
0
,
iNode
);
double
sI
=
localCoords_
(
1
,
iNode
);
double
tI
=
localCoords_
(
2
,
iNode
);
N_
(
ip
,
iNode
)
=
0.125
*
(
1
+
r
*
rI
)
*
(
1
+
s
*
sI
)
*
(
1
+
t
*
tI
);
}
}
// Integration point weights
ipWeights_
=
1.0
;
// integration points by face
ipFace2DCoords_
(
0
,
0
)
=-
a
;
ipFace2DCoords_
(
1
,
0
)
=-
a
;
ipFace2DCoords_
(
0
,
1
)
=
a
;
ipFace2DCoords_
(
1
,
1
)
=-
a
;
ipFace2DCoords_
(
0
,
2
)
=
a
;
ipFace2DCoords_
(
1
,
2
)
=
a
;
ipFace2DCoords_
(
0
,
3
)
=-
a
;
ipFace2DCoords_
(
1
,
3
)
=
a
;
ipFaceCoords_
[
0
](
0
,
0
)
=-
1
;
ipFaceCoords_
[
0
](
1
,
0
)
=-
a
;
ipFaceCoords_
[
0
](
2
,
0
)
=-
a
;
ipFaceCoords_
[
0
](
0
,
1
)
=-
1
;
ipFaceCoords_
[
0
](
1
,
1
)
=
a
;
ipFaceCoords_
[
0
](
2
,
1
)
=-
a
;
ipFaceCoords_
[
0
](
0
,
2
)
=-
1
;
ipFaceCoords_
[
0
](
1
,
2
)
=
a
;
ipFaceCoords_
[
0
](
2
,
2
)
=
a
;
ipFaceCoords_
[
0
](
0
,
3
)
=-
1
;
ipFaceCoords_
[
0
](
1
,
3
)
=-
a
;
ipFaceCoords_
[
0
](
2
,
3
)
=
a
;
ipFaceCoords_
[
1
](
0
,
0
)
=
1
;
ipFaceCoords_
[
1
](
1
,
0
)
=-
a
;
ipFaceCoords_
[
1
](
2
,
0
)
=-
a
;
ipFaceCoords_
[
1
](
0
,
1
)
=
1
;
ipFaceCoords_
[
1
](
1
,
1
)
=
a
;
ipFaceCoords_
[
1
](
2
,
1
)
=-
a
;
ipFaceCoords_
[
1
](
0
,
2
)
=
1
;
ipFaceCoords_
[
1
](
1
,
2
)
=
a
;
ipFaceCoords_
[
1
](
2
,
2
)
=
a
;
ipFaceCoords_
[
1
](
0
,
3
)
=
1
;
ipFaceCoords_
[
1
](
1
,
3
)
=-
a
;
ipFaceCoords_
[
1
](
2
,
3
)
=
a
;
ipFaceCoords_
[
2
](
0
,
0
)
=-
a
;
ipFaceCoords_
[
2
](
1
,
0
)
=-
1
;
ipFaceCoords_
[
2
](
2
,
0
)
=-
a
;
ipFaceCoords_
[
2
](
0
,
1
)
=-
a
;
ipFaceCoords_
[
2
](
1
,
1
)
=-
1
;
ipFaceCoords_
[
2
](
2
,
1
)
=
a
;
ipFaceCoords_
[
2
](
0
,
2
)
=
a
;
ipFaceCoords_
[
2
](
1
,
2
)
=-
1
;
ipFaceCoords_
[
2
](
2
,
2
)
=
a
;
ipFaceCoords_
[
2
](
0
,
3
)
=
a
;
ipFaceCoords_
[
2
](
1
,
3
)
=-
1
;
ipFaceCoords_
[
2
](
2
,
3
)
=-
a
;
ipFaceCoords_
[
3
](
0
,
0
)
=-
a
;
ipFaceCoords_
[
3
](
1
,
0
)
=
1
;
ipFaceCoords_
[
3
](
2
,
0
)
=-
a
;
ipFaceCoords_
[
3
](
0
,
1
)
=-
a
;
ipFaceCoords_
[
3
](
1
,
1
)
=
1
;
ipFaceCoords_
[
3
](
2
,
1
)
=
a
;
ipFaceCoords_
[
3
](
0
,
2
)
=
a
;
ipFaceCoords_
[
3
](
1
,
2
)
=
1
;
ipFaceCoords_
[
3
](
2
,
2
)
=
a
;
ipFaceCoords_
[
3
](
0
,
3
)
=
a
;
ipFaceCoords_
[
3
](
1
,
3
)
=
1
;
ipFaceCoords_
[
2
](
2
,
3
)
=-
a
;
ipFaceCoords_
[
4
](
0
,
0
)
=-
a
;
ipFaceCoords_
[
4
](
1
,
0
)
=-
a
;
ipFaceCoords_
[
4
](
2
,
0
)
=-
1
;
ipFaceCoords_
[
4
](
0
,
1
)
=
a
;
ipFaceCoords_
[
4
](
1
,
1
)
=-
a
;
ipFaceCoords_
[
4
](
2
,
1
)
=-
1
;
ipFaceCoords_
[
4
](
0
,
2
)
=
a
;
ipFaceCoords_
[
4
](
1
,
2
)
=
a
;
ipFaceCoords_
[
4
](
2
,
2
)
=-
1
;
ipFaceCoords_
[
4
](
0
,
3
)
=-
a
;
ipFaceCoords_
[
4
](
1
,
3
)
=
a
;
ipFaceCoords_
[
4
](
2
,
3
)
=-
1
;
ipFaceCoords_
[
5
](
0
,
0
)
=-
a
;
ipFaceCoords_
[
5
](
1
,
0
)
=-
a
;
ipFaceCoords_
[
5
](
2
,
0
)
=
1
;
ipFaceCoords_
[
5
](
0
,
1
)
=
a
;
ipFaceCoords_
[
5
](
1
,
1
)
=-
a
;
ipFaceCoords_
[
5
](
2
,
1
)
=
1
;
ipFaceCoords_
[
5
](
0
,
2
)
=
a
;
ipFaceCoords_
[
5
](
1
,
2
)
=
a
;
ipFaceCoords_
[
5
](
2
,
2
)
=
1
;
ipFaceCoords_
[
5
](
0
,
3
)
=-
a
;
ipFaceCoords_
[
5
](
1
,
3
)
=
a
;
ipFaceCoords_
[
5
](
2
,
3
)
=
1
;
// Compute shape functions at ip's
for
(
int
f
=
0
;
f
<
numFaces_
;
f
++
)
{
for
(
int
ip
=
0
;
ip
<
numFaceIPs_
;
ip
++
)
{
double
r
=
ipFaceCoords_
[
f
](
0
,
ip
);
double
s
=
ipFaceCoords_
[
f
](
1
,
ip
);
double
t
=
ipFaceCoords_
[
f
](
2
,
ip
);
for
(
int
iNode
=
0
;
iNode
<
numEltNodes_
;
iNode
++
)
{
double
rI
=
localCoords_
(
0
,
iNode
);
double
sI
=
localCoords_
(
1
,
iNode
);
double
tI
=
localCoords_
(
2
,
iNode
);
NFace_
[
f
](
ip
,
iNode
)
=
0.125
*
(
1
+
r
*
rI
)
*
(
1
+
s
*
sI
)
*
(
1
+
t
*
tI
);
}
}
}
// integration points
ipFaceWeights_
=
1.0
;
}
FE_ElementHex
::~
FE_ElementHex
()
{
// Nothing to do here
}
// -------------------------------------------------------------
// shape_function
// -------------------------------------------------------------
void
FE_ElementHex
::
shape_function
(
const
int
eltID
,
DENS_MAT
&
N
,
vector
<
DENS_MAT
>
&
dN
,
DIAG_MAT
&
weights
)
{
// Get element node coordinates from mesh
DENS_MAT
xCoords
;
feMesh_
->
element_coordinates
(
eltID
,
xCoords
);
// Transpose xCoords
DENS_MAT
xCoordsT
(
transpose
(
xCoords
));
// Shape functions are simply the canonical element values
N
=
N_
;
// Set sizes of matrices and vectors
if
((
int
)
dN
.
size
()
!=
nSD_
)
dN
.
resize
(
nSD_
);
for
(
int
isd
=
0
;
isd
<
nSD_
;
isd
++
)
dN
[
isd
].
resize
(
numIPs_
,
numEltNodes_
);
weights
.
resize
(
numIPs_
,
numIPs_
);
// Create some temporary matrices:
// Shape function deriv.'s w.r.t. element coords
DENS_MAT
dNdr
(
nSD_
,
numEltNodes_
,
false
);
// Jacobian matrix: [dx/dr dy/ds dz/dt | dx/ds ... ]
DENS_MAT
dxdr
,
drdx
,
dNdx
;
// Loop over integration points
for
(
int
ip
=
0
;
ip
<
numIPs_
;
ip
++
)
{
const
double
&
r
=
ipCoords_
(
0
,
ip
);
const
double
&
s
=
ipCoords_
(
1
,
ip
);
const
double
&
t
=
ipCoords_
(
2
,
ip
);
// Fill dNdr
for
(
int
inode
=
0
;
inode
<
numEltNodes_
;
inode
++
)
{
const
double
&
rI
=
localCoords_
(
0
,
inode
);
const
double
&
sI
=
localCoords_
(
1
,
inode
);
const
double
&
tI
=
localCoords_
(
2
,
inode
);
dNdr
(
0
,
inode
)
=
0.125
*
rI
*
(
1
+
s
*
sI
)
*
(
1
+
t
*
tI
);
dNdr
(
1
,
inode
)
=
0.125
*
sI
*
(
1
+
r
*
rI
)
*
(
1
+
t
*
tI
);
dNdr
(
2
,
inode
)
=
0.125
*
tI
*
(
1
+
r
*
rI
)
*
(
1
+
s
*
sI
);
}
// Compute dx/dxi matrix
dxdr
=
dNdr
*
xCoordsT
;
drdx
=
inv
(
dxdr
);
// Compute dNdx and fill dN matrix
dNdx
=
drdx
*
dNdr
;
for
(
int
isd
=
0
;
isd
<
nSD_
;
isd
++
)
for
(
int
inode
=
0
;
inode
<
numEltNodes_
;
inode
++
)
dN
[
isd
](
ip
,
inode
)
=
dNdx
(
isd
,
inode
);
// Compute jacobian determinant of dxdr at this ip
double
J
=
dxdr
(
0
,
0
)
*
(
dxdr
(
1
,
1
)
*
dxdr
(
2
,
2
)
-
dxdr
(
2
,
1
)
*
dxdr
(
1
,
2
)
)
-
dxdr
(
0
,
1
)
*
(
dxdr
(
1
,
0
)
*
dxdr
(
2
,
2
)
-
dxdr
(
1
,
2
)
*
dxdr
(
2
,
0
)
)
+
dxdr
(
0
,
2
)
*
(
dxdr
(
1
,
0
)
*
dxdr
(
2
,
1
)
-
dxdr
(
1
,
1
)
*
dxdr
(
2
,
0
)
);
// Compute ip weight
weights
(
ip
,
ip
)
=
ipWeights_
(
ip
)
*
J
;
}
}
//-----------------------------------------------------------------
void
FE_ElementHex
::
shape_function
(
const
int
eltID
,
const
VECTOR
&
xi
,
DENS_VEC
&
N
)
{
N
.
resize
(
numEltNodes_
);
for
(
int
inode
=
0
;
inode
<
numEltNodes_
;
++
inode
)
{
double
shp
=
0.125
;
for
(
int
i
=
0
;
i
<
nSD_
;
++
i
)
shp
*=
(
1
+
xi
(
i
)
*
localCoords_
(
i
,
inode
));
N
(
inode
)
=
shp
;
}
}
//-----------------------------------------------------------------
void
FE_ElementHex
::
shape_function
(
const
int
eltID
,
const
VECTOR
&
xi
,
DENS_VEC
&
N
,
DENS_MAT
&
dNdx
)
{
N
.
resize
(
numEltNodes_
);
// Get element node coordinates from mesh
DENS_MAT
xCoords
;
feMesh_
->
element_coordinates
(
eltID
,
xCoords
);
DENS_MAT
xCoordsT
=
transpose
(
xCoords
);
DENS_MAT
dNdr
(
3
,
8
,
false
),
dxdr
(
3
,
3
,
false
),
drdx
;
double
r
=
xi
(
0
);
double
s
=
xi
(
1
);
double
t
=
xi
(
2
);
for
(
int
inode
=
0
;
inode
<
numEltNodes_
;
++
inode
)
{
double
rI
=
localCoords_
(
0
,
inode
);
double
sI
=
localCoords_
(
1
,
inode
);
double
tI
=
localCoords_
(
2
,
inode
);
// Shape function
N
(
inode
)
=
0.125
*
(
1.0
+
r
*
rI
)
*
(
1.0
+
s
*
sI
)
*
(
1.0
+
t
*
tI
);
// Shape function derivative wrt (r,s,t)
dNdr
(
0
,
inode
)
=
0.125
*
rI
*
(
1.0
+
s
*
sI
)
*
(
1.0
+
t
*
tI
);
dNdr
(
1
,
inode
)
=
0.125
*
sI
*
(
1.0
+
r
*
rI
)
*
(
1.0
+
t
*
tI
);
dNdr
(
2
,
inode
)
=
0.125
*
tI
*
(
1.0
+
r
*
rI
)
*
(
1.0
+
s
*
sI
);
}
// Derivatives wrt (x,y,z)
dxdr
=
dNdr
*
xCoordsT
;
drdx
=
inv
(
dxdr
);
dNdx
=
drdx
*
dNdr
;
}
//-----------------------------------------------------------------
void
FE_ElementHex
::
face_shape_function
(
const
PAIR
&
face
,
DENS_MAT
&
N
,
vector
<
DENS_MAT
>
&
dN
,
vector
<
DENS_MAT
>
&
Nn
,
DIAG_MAT
&
weights
)
{
int
eltID
=
face
.
first
;
int
local_faceID
=
face
.
second
;
// const double * face_normal = face.normal();
// double * normal = face_normal;
DENS_VEC
normal
(
nSD_
);
// Get element node coordinates from mesh
DENS_MAT
xCoords
;
feMesh_
->
element_coordinates
(
eltID
,
xCoords
);
// Transpose xCoords
DENS_MAT
xCoordsT
=
transpose
(
xCoords
);
// Shape functions are simply the canonical element values
N
.
reset
(
numFaceIPs_
,
numEltNodes_
);
N
=
NFace_
[
local_faceID
];
// Set sizes of matrices and vectors
if
((
int
)
dN
.
size
()
!=
nSD_
)
dN
.
resize
(
nSD_
);
for
(
int
isd
=
0
;
isd
<
nSD_
;
isd
++
)
dN
[
isd
].
reset
(
numFaceIPs_
,
numEltNodes_
);
weights
.
reset
(
numFaceIPs_
,
numFaceIPs_
);
// Create some temporary matrices:
// Shape function deriv.'s w.r.t. element coords
DENS_MAT
dNdr
(
nSD_
,
numEltNodes_
);
// Jacobian matrix: [dx/dr dy/ds dz/dt | dx/ds ... ]
DENS_MAT
dxdr
,
drdx
,
dNdx
;
// Loop over integration points
for
(
int
ip
=
0
;
ip
<
numFaceIPs_
;
ip
++
)
{
double
r
=
ipFaceCoords_
[
local_faceID
](
0
,
ip
);
double
s
=
ipFaceCoords_
[
local_faceID
](
1
,
ip
);
double
t
=
ipFaceCoords_
[
local_faceID
](
2
,
ip
);
// Fill dNdr
for
(
int
inode
=
0
;
inode
<
numEltNodes_
;
inode
++
)
{
double
rI
=
localCoords_
(
0
,
inode
);
double
sI
=
localCoords_
(
1
,
inode
);
double
tI
=
localCoords_
(
2
,
inode
);
dNdr
(
0
,
inode
)
=
0.125
*
rI
*
(
1
+
s
*
sI
)
*
(
1
+
t
*
tI
);
dNdr
(
1
,
inode
)
=
0.125
*
sI
*
(
1
+
r
*
rI
)
*
(
1
+
t
*
tI
);
dNdr
(
2
,
inode
)
=
0.125
*
tI
*
(
1
+
r
*
rI
)
*
(
1
+
s
*
sI
);
}
// Compute dx/dxi matrix
dxdr
=
dNdr
*
xCoordsT
;
drdx
=
inv
(
dxdr
);
// Compute 2d jacobian determinant of dxdr at this ip
double
J
=
face_normal
(
face
,
ip
,
normal
);
// Compute dNdx and fill dN matrix
dNdx
=
drdx
*
dNdr
;
for
(
int
isd
=
0
;
isd
<
nSD_
;
isd
++
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes_
;
inode
++
)
{
dN
[
isd
](
ip
,
inode
)
=
dNdx
(
isd
,
inode
);
Nn
[
isd
](
ip
,
inode
)
=
N
(
ip
,
inode
)
*
normal
(
isd
);
}
}
// Compute ip weight
weights
(
ip
,
ip
)
=
ipFaceWeights_
(
ip
)
*
J
;
}
}
//-----------------------------------------------------------------
void
FE_ElementHex
::
face_shape_function
(
const
PAIR
&
face
,
DENS_MAT
&
N
,
DENS_MAT
&
n
,
DIAG_MAT
&
weights
)
{
int
eltID
=
face
.
first
;
int
local_faceID
=
face
.
second
;
// const double * face_normal = face.normal();
// double * normal = face_normal;
DENS_VEC
normal
(
nSD_
);
// Get element node coordinates from mesh
DENS_MAT
xCoords
;
feMesh_
->
element_coordinates
(
eltID
,
xCoords
);
// Transpose xCoords
DENS_MAT
xCoordsT
=
transpose
(
xCoords
);
// Shape functions are simply the canonical element values
N
.
reset
(
numFaceIPs_
,
numEltNodes_
);
N
=
NFace_
[
local_faceID
];
weights
.
reset
(
numFaceIPs_
,
numFaceIPs_
);
// Create some temporary matrices:
// Shape function deriv.'s w.r.t. element coords
DENS_MAT
dNdr
(
nSD_
,
numEltNodes_
);
// Jacobian matrix: [dx/dr dy/ds dz/dt | dx/ds ... ]
DENS_MAT
dxdr
,
drdx
,
dNdx
;
// Loop over integration points
for
(
int
ip
=
0
;
ip
<
numFaceIPs_
;
ip
++
)
{
double
r
=
ipFaceCoords_
[
local_faceID
](
0
,
ip
);
double
s
=
ipFaceCoords_
[
local_faceID
](
1
,
ip
);
double
t
=
ipFaceCoords_
[
local_faceID
](
2
,
ip
);
// Compute 2d jacobian determinant of dxdr at this ip
double
J
=
face_normal
(
face
,
ip
,
normal
);
// Copy normal at integration point
for
(
int
isd
=
0
;
isd
<
nSD_
;
isd
++
)
{
n
(
isd
,
ip
)
=
normal
(
isd
);
}
// Compute ip weight
weights
(
ip
,
ip
)
=
ipFaceWeights_
(
ip
)
*
J
;
}
}
// -------------------------------------------------------------
// face_normal
// -------------------------------------------------------------
double
FE_ElementHex
::
face_normal
(
const
PAIR
&
face
,
int
ip
,
DENS_VEC
&
normal
)
{
// Get element node coordinates from mesh
DENS_MAT
xCoords
;
feMesh_
->
face_coordinates
(
face
,
xCoords
);
// Transpose xCoords
DENS_MAT
xCoordsT
=
transpose
(
xCoords
);
double
r
=
ipFace2DCoords_
(
0
,
ip
);
double
s
=
ipFace2DCoords_
(
1
,
ip
);
DENS_MAT
dNdr
(
nSD_
-
1
,
numFaceNodes_
);
for
(
int
inode
=
0
;
inode
<
numFaceNodes_
;
inode
++
)
{
double
rI
=
localCoords_
(
0
,
inode
);
// NOTE a little dangerous
double
sI
=
localCoords_
(
1
,
inode
);
dNdr
(
0
,
inode
)
=
0.25
*
rI
*
(
1
+
s
*
sI
);
dNdr
(
1
,
inode
)
=
0.25
*
sI
*
(
1
+
r
*
rI
);
}
DENS_MAT
dxdr
(
2
,
3
);
// Compute dx/dxi matrix
dxdr
=
dNdr
*
xCoordsT
;
normal
(
0
)
=
dxdr
(
0
,
1
)
*
dxdr
(
1
,
2
)
-
dxdr
(
0
,
2
)
*
dxdr
(
1
,
1
)
;
normal
(
1
)
=
dxdr
(
0
,
2
)
*
dxdr
(
1
,
0
)
-
dxdr
(
0
,
0
)
*
dxdr
(
1
,
2
)
;
normal
(
2
)
=
dxdr
(
0
,
0
)
*
dxdr
(
1
,
1
)
-
dxdr
(
0
,
1
)
*
dxdr
(
1
,
0
)
;
double
J
=
sqrt
(
normal
(
0
)
*
normal
(
0
)
+
normal
(
1
)
*
normal
(
1
)
+
normal
(
2
)
*
normal
(
2
));
// double inv_J = normal_sense(face)/J;
double
inv_J
=
1.0
/
J
;
normal
(
0
)
*=
inv_J
;
normal
(
1
)
*=
inv_J
;
normal
(
2
)
*=
inv_J
;
return
J
;
}
};
// namespace ATC_Transfer
Event Timeline
Log In to Comment