Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78664262
FE_Interpolate.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
Thu, Aug 22, 06:50
Size
34 KB
Mime Type
text/x-c
Expires
Sat, Aug 24, 06:50 (2 d)
Engine
blob
Format
Raw Data
Handle
20082966
Attached To
rLAMMPS lammps
FE_Interpolate.cpp
View Options
// ATC header files
#include "ATC_Error.h"
#include "FE_Element.h"
#include "FE_Interpolate.h"
#include "FE_Quadrature.h"
// Other headers
#include "math.h"
using
std
::
map
;
using
std
::
vector
;
namespace
ATC
{
FE_Interpolate
::
FE_Interpolate
(
FE_Element
*
feElement
)
:
feElement_
(
feElement
),
nSD_
(
feElement
->
num_dims
())
{
// Nothing to do here
}
FE_Interpolate
::~
FE_Interpolate
()
{
if
(
!
feQuadList_
.
empty
())
{
map
<
FeIntQuadrature
,
FE_Quadrature
*>::
iterator
qit
;
for
(
qit
=
feQuadList_
.
begin
();
qit
!=
feQuadList_
.
end
();
++
qit
)
{
delete
(
qit
->
second
);
}
}
}
void
FE_Interpolate
::
set_quadrature
(
FeEltGeometry
geo
,
FeIntQuadrature
quad
)
{
if
(
feQuadList_
.
count
(
quad
)
==
0
)
{
feQuad_
=
new
FE_Quadrature
(
geo
,
quad
);
feQuadList_
[
quad
]
=
feQuad_
;
}
else
{
feQuad_
=
feQuadList_
[
quad
];
}
precalculate_shape_functions
();
}
void
FE_Interpolate
::
precalculate_shape_functions
()
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
int
numFaces
=
feElement_
->
num_faces
();
int
numFaceNodes
=
feElement_
->
num_face_nodes
();
int
numIPs
=
feQuad_
->
numIPs
;
DENS_MAT
&
ipCoords
=
feQuad_
->
ipCoords
;
int
numFaceIPs
=
feQuad_
->
numFaceIPs
;
vector
<
DENS_MAT
>
&
ipFaceCoords
=
feQuad_
->
ipFaceCoords
;
DENS_MAT
&
ipFace2DCoords
=
feQuad_
->
ipFace2DCoords
;
// Compute elemental shape functions at ips
N_
.
reset
(
numIPs
,
numEltNodes
);
dNdr_
.
assign
(
numIPs
,
DENS_MAT
(
nSD_
,
numEltNodes
));
for
(
int
ip
=
0
;
ip
<
numIPs
;
ip
++
)
{
CLON_VEC
thisIP
=
column
(
ipCoords
,
ip
);
CLON_VEC
thisN
=
row
(
N_
,
ip
);
DENS_MAT
&
thisdNdr
=
dNdr_
[
ip
];
compute_N
(
thisIP
,
thisN
);
compute_N_dNdr
(
thisIP
,
thisN
,
thisdNdr
);
}
// Compute face shape functions at ip's
NFace_
.
assign
(
numFaces
,
DENS_MAT
(
numFaceIPs
,
numEltNodes
));
dNdrFace_
.
resize
(
numFaces
);
for
(
int
f
=
0
;
f
<
numFaces
;
f
++
)
{
dNdrFace_
[
f
].
assign
(
numIPs
,
DENS_MAT
(
nSD_
,
numEltNodes
));
}
for
(
int
f
=
0
;
f
<
numFaces
;
f
++
)
{
for
(
int
ip
=
0
;
ip
<
numFaceIPs
;
ip
++
)
{
CLON_VEC
thisIP
=
column
(
ipFaceCoords
[
f
],
ip
);
CLON_VEC
thisN
=
row
(
NFace_
[
f
],
ip
);
DENS_MAT
&
thisdNdr
=
dNdrFace_
[
f
][
ip
];
compute_N_dNdr
(
thisIP
,
thisN
,
thisdNdr
);
}
}
// Compute 2D face shape function derivatives
dNdrFace2D_
.
assign
(
numFaceIPs
,
DENS_MAT
(
nSD_
-
1
,
numFaceNodes
));
for
(
int
ip
=
0
;
ip
<
numFaceIPs
;
ip
++
)
{
CLON_VEC
thisIP
=
column
(
ipFace2DCoords
,
ip
);
DENS_MAT
&
thisdNdr
=
dNdrFace2D_
[
ip
];
compute_dNdr
(
thisIP
,
numFaceNodes
,
nSD_
-
1
,
feElement_
->
face_area
(),
thisdNdr
);
}
}
//-----------------------------------------------------------------
// shape function value at a particular point given local coordinates
//-----------------------------------------------------------------
void
FE_Interpolate
::
shape_function
(
const
VECTOR
&
xi
,
DENS_VEC
&
N
)
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
N
.
resize
(
numEltNodes
);
compute_N
(
xi
,
N
);
}
void
FE_Interpolate
::
shape_function
(
const
DENS_MAT
&
eltCoords
,
const
VECTOR
&
xi
,
DENS_VEC
&
N
,
DENS_MAT
&
dNdx
)
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
N
.
resize
(
numEltNodes
);
DENS_MAT
dNdr
(
nSD_
,
numEltNodes
,
false
);
compute_N_dNdr
(
xi
,
N
,
dNdr
);
DENS_MAT
eltCoordsT
=
transpose
(
eltCoords
);
DENS_MAT
dxdr
,
drdx
;
dxdr
=
dNdr
*
eltCoordsT
;
drdx
=
inv
(
dxdr
);
dNdx
=
drdx
*
dNdr
;
}
void
FE_Interpolate
::
shape_function_derivatives
(
const
DENS_MAT
&
eltCoords
,
const
VECTOR
&
xi
,
DENS_MAT
&
dNdx
)
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
DENS_MAT
dNdr
(
nSD_
,
numEltNodes
,
false
);
DENS_VEC
N
(
numEltNodes
);
compute_N_dNdr
(
xi
,
N
,
dNdr
);
DENS_MAT
eltCoordsT
=
transpose
(
eltCoords
);
DENS_MAT
dxdr
,
drdx
;
dxdr
=
dNdr
*
eltCoordsT
;
// tangents or Jacobian matrix
drdx
=
inv
(
dxdr
);
dNdx
=
drdx
*
dNdr
;
// dN/dx = dN/dxi (dx/dxi)^-1
}
void
FE_Interpolate
::
tangents
(
const
DENS_MAT
&
eltCoords
,
const
VECTOR
&
xi
,
DENS_MAT
&
dxdr
)
const
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
DENS_MAT
dNdr
(
nSD_
,
numEltNodes
,
false
);
DENS_VEC
N
(
numEltNodes
);
compute_N_dNdr
(
xi
,
N
,
dNdr
);
//dNdr.print("dNdr");
DENS_MAT
eltCoordsT
=
transpose
(
eltCoords
);
//eltCoordsT.print("elt coords");
dxdr
=
dNdr
*
eltCoordsT
;
//dxdr.print("dxdr");
}
void
FE_Interpolate
::
tangents
(
const
DENS_MAT
&
eltCoords
,
const
VECTOR
&
xi
,
vector
<
DENS_VEC
>
&
dxdxi
,
const
bool
normalize
)
const
{
DENS_MAT
dxdr
;
tangents
(
eltCoords
,
xi
,
dxdr
);
//dxdr.print("dxdr-post");
dxdxi
.
resize
(
nSD_
);
//for (int i = 0; i < nSD_; ++i) dxdxi[i] = CLON_VEC(dxdr,CLONE_COL,i);
for
(
int
i
=
0
;
i
<
nSD_
;
++
i
)
{
dxdxi
[
i
].
resize
(
nSD_
);
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
dxdxi
[
i
](
j
)
=
dxdr
(
i
,
j
);
}
}
//dxdxi[0].print("t1");
//dxdxi[1].print("t2");
//dxdxi[2].print("t3");
if
(
normalize
)
{
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
double
norm
=
0
;
VECTOR
&
t
=
dxdxi
[
j
];
for
(
int
i
=
0
;
i
<
nSD_
;
++
i
)
norm
+=
t
(
i
)
*
t
(
i
);
norm
=
1.
/
sqrt
(
norm
);
for
(
int
i
=
0
;
i
<
nSD_
;
++
i
)
t
(
i
)
*=
norm
;
}
}
}
// -------------------------------------------------------------
// shape_function values at nodes
// -------------------------------------------------------------
void
FE_Interpolate
::
shape_function
(
const
DENS_MAT
&
eltCoords
,
DENS_MAT
&
N
,
vector
<
DENS_MAT
>
&
dN
,
DIAG_MAT
&
weights
)
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
// Transpose eltCoords
DENS_MAT
eltCoordsT
(
transpose
(
eltCoords
));
// 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
(
feQuad_
->
numIPs
,
numEltNodes
);
weights
.
resize
(
feQuad_
->
numIPs
,
feQuad_
->
numIPs
);
// Create some temporary matrices:
// Jacobian matrix: [dx/dr dy/ds dz/dt | dx/ds ... ]
DENS_MAT
dxdr
,
drdx
,
dNdx
;
// Loop over integration points
for
(
int
ip
=
0
;
ip
<
feQuad_
->
numIPs
;
ip
++
)
{
// Compute dx/dxi matrix
dxdr
=
dNdr_
[
ip
]
*
eltCoordsT
;
drdx
=
inv
(
dxdr
);
// Compute dNdx and fill dN matrix
dNdx
=
drdx
*
dNdr_
[
ip
];
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
)
=
feQuad_
->
ipWeights
(
ip
)
*
J
;
}
}
//-----------------------------------------------------------------
// shape functions on a given face
//-----------------------------------------------------------------
void
FE_Interpolate
::
face_shape_function
(
const
DENS_MAT
&
eltCoords
,
const
DENS_MAT
&
faceCoords
,
const
int
faceID
,
DENS_MAT
&
N
,
DENS_MAT
&
n
,
DIAG_MAT
&
weights
)
{
int
numFaceIPs
=
feQuad_
->
numFaceIPs
;
// Transpose eltCoords
DENS_MAT
eltCoordsT
=
transpose
(
eltCoords
);
// Shape functions are simply the canonical element values
N
=
NFace_
[
faceID
];
// Create some temporary matrices:
// Jacobian matrix: [dx/dr dy/ds dz/dt | dx/ds ... ]
DENS_MAT
dxdr
,
drdx
,
dNdx
;
// Loop over integration points
DENS_VEC
normal
(
nSD_
);
n
.
resize
(
nSD_
,
numFaceIPs
);
weights
.
resize
(
numFaceIPs
,
numFaceIPs
);
for
(
int
ip
=
0
;
ip
<
numFaceIPs
;
ip
++
)
{
// Compute 2d jacobian determinant of dxdr at this ip
double
J
=
face_normal
(
faceCoords
,
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
)
=
feQuad_
->
ipFaceWeights
(
ip
)
*
J
;
}
}
void
FE_Interpolate
::
face_shape_function
(
const
DENS_MAT
&
eltCoords
,
const
DENS_MAT
&
faceCoords
,
const
int
faceID
,
DENS_MAT
&
N
,
vector
<
DENS_MAT
>
&
dN
,
vector
<
DENS_MAT
>
&
Nn
,
DIAG_MAT
&
weights
)
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
int
numFaceIPs
=
feQuad_
->
numFaceIPs
;
// Transpose eltCoords
DENS_MAT
eltCoordsT
=
transpose
(
eltCoords
);
// Shape functions are simply the canonical element values
N
=
NFace_
[
faceID
];
// Set sizes of matrices and vectors
if
((
int
)
dN
.
size
()
!=
nSD_
)
dN
.
resize
(
nSD_
);
if
((
int
)
Nn
.
size
()
!=
nSD_
)
Nn
.
resize
(
nSD_
);
for
(
int
isd
=
0
;
isd
<
nSD_
;
isd
++
)
{
dN
[
isd
].
resize
(
numFaceIPs
,
numEltNodes
);
Nn
[
isd
].
resize
(
numFaceIPs
,
numEltNodes
);
}
weights
.
resize
(
numFaceIPs
,
numFaceIPs
);
// Create some temporary matrices:
// Jacobian matrix: [dx/dr dy/ds dz/dt | dx/ds ... ]
DENS_MAT
dxdr
,
drdx
,
dNdx
;
DENS_VEC
normal
(
nSD_
);
// Loop over integration points
for
(
int
ip
=
0
;
ip
<
numFaceIPs
;
ip
++
)
{
// Compute dx/dxi matrix
dxdr
=
dNdrFace_
[
faceID
][
ip
]
*
eltCoordsT
;
drdx
=
inv
(
dxdr
);
// Compute 2d jacobian determinant of dxdr at this ip
double
J
=
face_normal
(
faceCoords
,
ip
,
normal
);
// Compute dNdx and fill dN matrix
dNdx
=
drdx
*
dNdrFace_
[
faceID
][
ip
];
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
)
=
feQuad_
->
ipFaceWeights
(
ip
)
*
J
;
}
}
// -------------------------------------------------------------
// face normal
// -------------------------------------------------------------
double
FE_Interpolate
::
face_normal
(
const
DENS_MAT
&
faceCoords
,
int
ip
,
DENS_VEC
&
normal
)
{
// Compute dx/dr for deformed element
DENS_MAT
faceCoordsT
=
transpose
(
faceCoords
);
DENS_MAT
dxdr
=
dNdrFace2D_
[
ip
]
*
faceCoordsT
;
// Normal vector from cross product, hardcoded for 3D, sad
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
);
// Jacobian (3D)
double
J
=
sqrt
(
normal
(
0
)
*
normal
(
0
)
+
normal
(
1
)
*
normal
(
1
)
+
normal
(
2
)
*
normal
(
2
));
double
invJ
=
1.0
/
J
;
normal
(
0
)
*=
invJ
;
normal
(
1
)
*=
invJ
;
normal
(
2
)
*=
invJ
;
return
J
;
}
int
FE_Interpolate
::
num_ips
()
const
{
return
feQuad_
->
numIPs
;
}
int
FE_Interpolate
::
num_face_ips
()
const
{
return
feQuad_
->
numFaceIPs
;
}
/*********************************************************
* Class FE_InterpolateCartLagrange
*
* For computing Lagrange shape functions using Cartesian
* coordinate systems (all quads/hexes fall under this
* category, and any elements derived by degenerating
* them). Not to be used for serendipity elements, which
* should be implemented for SPEED.
*
*********************************************************/
FE_InterpolateCartLagrange
::
FE_InterpolateCartLagrange
(
FE_Element
*
feElement
)
:
FE_Interpolate
(
feElement
)
{
set_quadrature
(
HEXA
,
GAUSS2
);
}
FE_InterpolateCartLagrange
::~
FE_InterpolateCartLagrange
()
{
// Handled by base class
}
void
FE_InterpolateCartLagrange
::
compute_N
(
const
VECTOR
&
point
,
VECTOR
&
N
)
{
// *** see comments for compute_N_dNdr ***
const
DENS_VEC
&
localCoords1d
=
feElement_
->
local_coords_1d
();
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
int
numEltNodes1d
=
feElement_
->
num_elt_nodes_1d
();
DENS_MAT
lagrangeTerms
(
nSD_
,
numEltNodes1d
);
DENS_MAT
lagrangeDenom
(
nSD_
,
numEltNodes1d
);
lagrangeTerms
=
1.0
;
lagrangeDenom
=
1.0
;
for
(
int
iSD
=
0
;
iSD
<
nSD_
;
++
iSD
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes1d
;
++
inode
)
{
for
(
int
icont
=
0
;
icont
<
numEltNodes1d
;
++
icont
)
{
if
(
inode
!=
icont
)
{
lagrangeDenom
(
iSD
,
inode
)
*=
(
localCoords1d
(
inode
)
-
localCoords1d
(
icont
));
lagrangeTerms
(
iSD
,
inode
)
*=
(
point
(
iSD
)
-
localCoords1d
(
icont
));
}
}
}
}
for
(
int
iSD
=
0
;
iSD
<
nSD_
;
++
iSD
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes1d
;
++
inode
)
{
lagrangeTerms
(
iSD
,
inode
)
/=
lagrangeDenom
(
iSD
,
inode
);
}
}
N
=
1.0
;
vector
<
int
>
mapping
(
nSD_
);
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
feElement_
->
mapping
(
inode
,
mapping
);
for
(
int
iSD
=
0
;
iSD
<
nSD_
;
++
iSD
)
{
N
(
inode
)
*=
lagrangeTerms
(
iSD
,
mapping
[
iSD
]);
}
}
}
// Sort of a test-ride for a generic version that can be used for
// faces too. The only thing that's not "generic" is localCoords,
// which very magically works in both cases.
void
FE_InterpolateCartLagrange
::
compute_dNdr
(
const
VECTOR
&
point
,
const
int
numNodes
,
const
int
nD
,
const
double
,
DENS_MAT
&
dNdr
)
{
// *** see comments for compute_N_dNdr ***
const
DENS_VEC
&
localCoords1d
=
feElement_
->
local_coords_1d
();
int
numEltNodes1d
=
feElement_
->
num_elt_nodes_1d
();
DENS_MAT
lagrangeTerms
(
nD
,
numEltNodes1d
);
DENS_MAT
lagrangeDenom
(
nD
,
numEltNodes1d
);
DENS_MAT
lagrangeDeriv
(
nD
,
numEltNodes1d
);
lagrangeDenom
=
1.0
;
lagrangeTerms
=
1.0
;
lagrangeDeriv
=
0.0
;
DENS_VEC
productRuleVec
(
numEltNodes1d
);
productRuleVec
=
1.0
;
for
(
int
iSD
=
0
;
iSD
<
nD
;
++
iSD
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes1d
;
++
inode
)
{
for
(
int
icont
=
0
;
icont
<
numEltNodes1d
;
++
icont
)
{
if
(
inode
!=
icont
)
{
lagrangeTerms
(
iSD
,
inode
)
*=
(
point
(
iSD
)
-
localCoords1d
(
icont
));
lagrangeDenom
(
iSD
,
inode
)
*=
(
localCoords1d
(
inode
)
-
localCoords1d
(
icont
));
for
(
int
dcont
=
0
;
dcont
<
numEltNodes1d
;
++
dcont
)
{
if
(
inode
==
dcont
)
{
productRuleVec
(
dcont
)
=
0.0
;
}
else
if
(
icont
==
dcont
)
{
}
else
{
productRuleVec
(
dcont
)
*=
(
point
(
iSD
)
-
localCoords1d
(
icont
));
}
}
}
}
for
(
int
dcont
=
0
;
dcont
<
numEltNodes1d
;
++
dcont
)
{
lagrangeDeriv
(
iSD
,
inode
)
+=
productRuleVec
(
dcont
);
}
productRuleVec
=
1.0
;
}
}
for
(
int
iSD
=
0
;
iSD
<
nD
;
++
iSD
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes1d
;
++
inode
)
{
lagrangeTerms
(
iSD
,
inode
)
/=
lagrangeDenom
(
iSD
,
inode
);
lagrangeDeriv
(
iSD
,
inode
)
/=
lagrangeDenom
(
iSD
,
inode
);
}
}
dNdr
=
1.0
;
vector
<
int
>
mapping
(
nD
);
for
(
int
inode
=
0
;
inode
<
numNodes
;
++
inode
)
{
feElement_
->
mapping
(
inode
,
mapping
);
for
(
int
iSD
=
0
;
iSD
<
nD
;
++
iSD
)
{
for
(
int
dSD
=
0
;
dSD
<
nD
;
++
dSD
)
{
if
(
iSD
==
dSD
)
{
dNdr
(
dSD
,
inode
)
*=
lagrangeDeriv
(
iSD
,
mapping
[
iSD
]);
}
else
{
dNdr
(
dSD
,
inode
)
*=
lagrangeTerms
(
iSD
,
mapping
[
iSD
]);
}
}
}
}
}
void
FE_InterpolateCartLagrange
::
compute_N_dNdr
(
const
VECTOR
&
point
,
VECTOR
&
N
,
DENS_MAT
&
dNdr
)
const
{
// Required data from element class
const
DENS_VEC
&
localCoords1d
=
feElement_
->
local_coords_1d
();
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
int
numEltNodes1d
=
feElement_
->
num_elt_nodes_1d
();
// lagrangeTerms stores the numerator for the various Lagrange polynomials
// in one dimension, that will be used to produce the three dimensional
// shape functions
DENS_MAT
lagrangeTerms
(
nSD_
,
numEltNodes1d
);
// lagrangeDenom stores the denominator. Stored separately to reduce
// redundancy, because it will be used for the shape functions and derivs
DENS_MAT
lagrangeDenom
(
nSD_
,
numEltNodes1d
);
// lagrangeDeriv stores the numerator for the derivative of the Lagrange
// polynomials
DENS_MAT
lagrangeDeriv
(
nSD_
,
numEltNodes1d
);
// Terms/Denom are products, Deriv will be a sum, so initialize as such:
lagrangeTerms
=
1.0
;
lagrangeDenom
=
1.0
;
lagrangeDeriv
=
0.0
;
// the derivative requires use of the product rule; to store the prodcuts
// which make up the terms produced by the product rule, we'll use this
// vector
DENS_VEC
productRuleVec
(
numEltNodes1d
);
productRuleVec
=
1.0
;
for
(
int
iSD
=
0
;
iSD
<
nSD_
;
++
iSD
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes1d
;
++
inode
)
{
for
(
int
icont
=
0
;
icont
<
numEltNodes1d
;
++
icont
)
{
if
(
inode
!=
icont
)
{
// each dimension and each 1d node per dimension has a
// contribution from all nodes besides the current node
lagrangeTerms
(
iSD
,
inode
)
*=
(
point
(
iSD
)
-
localCoords1d
(
icont
));
lagrangeDenom
(
iSD
,
inode
)
*=
(
localCoords1d
(
inode
)
-
localCoords1d
(
icont
));
// complciated; each sum produced by the product rule has one
// "derivative", and the rest are just identical to the terms
// above
for
(
int
dcont
=
0
;
dcont
<
numEltNodes1d
;
++
dcont
)
{
if
(
inode
==
dcont
)
{
// skip this term, derivative is 0
productRuleVec
(
dcont
)
=
0.0
;
}
else
if
(
icont
==
dcont
)
{
// no numerator contribution, derivative is 1
}
else
{
// part of the "constant"
productRuleVec
(
dcont
)
*=
(
point
(
iSD
)
-
localCoords1d
(
icont
));
}
}
}
}
// sum the terms produced by the product rule and store in Deriv
for
(
int
dcont
=
0
;
dcont
<
numEltNodes1d
;
++
dcont
)
{
lagrangeDeriv
(
iSD
,
inode
)
+=
productRuleVec
(
dcont
);
}
productRuleVec
=
1.0
;
}
}
// divide by denom
for
(
int
iSD
=
0
;
iSD
<
nSD_
;
++
iSD
)
{
for
(
int
inode
=
0
;
inode
<
numEltNodes1d
;
++
inode
)
{
lagrangeTerms
(
iSD
,
inode
)
/=
lagrangeDenom
(
iSD
,
inode
);
lagrangeDeriv
(
iSD
,
inode
)
/=
lagrangeDenom
(
iSD
,
inode
);
}
}
N
=
1.0
;
dNdr
=
1.0
;
// mapping returns the 1d nodes in each dimension that sould be multiplied
// to achieve the shape functions in 3d
vector
<
int
>
mapping
(
nSD_
);
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
feElement_
->
mapping
(
inode
,
mapping
);
for
(
int
iSD
=
0
;
iSD
<
nSD_
;
++
iSD
)
{
N
(
inode
)
*=
lagrangeTerms
(
iSD
,
mapping
[
iSD
]);
for
(
int
dSD
=
0
;
dSD
<
nSD_
;
++
dSD
)
{
// only use Deriv for the dimension in which we're taking the
// derivative, because the rest is essentially a "constant"
if
(
iSD
==
dSD
)
{
dNdr
(
dSD
,
inode
)
*=
lagrangeDeriv
(
iSD
,
mapping
[
iSD
]);
}
else
{
dNdr
(
dSD
,
inode
)
*=
lagrangeTerms
(
iSD
,
mapping
[
iSD
]);
}
}
}
}
}
/*********************************************************
* Class FE_InterpolateCartLin
*
* For computing linear shape functions using Cartesian
* coordinate systems (all quads/hexes fall under this
* category, and any elements derived by degenerating
* them).
*
*********************************************************/
FE_InterpolateCartLin
::
FE_InterpolateCartLin
(
FE_Element
*
feElement
)
:
FE_Interpolate
(
feElement
)
{
set_quadrature
(
HEXA
,
GAUSS2
);
}
FE_InterpolateCartLin
::~
FE_InterpolateCartLin
()
{
// Handled by base class
}
void
FE_InterpolateCartLin
::
compute_N
(
const
VECTOR
&
point
,
VECTOR
&
N
)
{
// *** see comments for compute_N_dNdr ***
const
DENS_MAT
&
localCoords
=
feElement_
->
local_coords
();
double
invVol
=
1.0
/
(
feElement_
->
vol
());
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
N
(
inode
)
=
invVol
;
for
(
int
isd
=
0
;
isd
<
nSD_
;
++
isd
)
{
N
(
inode
)
*=
(
1.0
+
point
(
isd
)
*
localCoords
(
isd
,
inode
));
}
}
}
// Sort of a test-ride for a generic version that can be used for
// faces too. The only thing that's not "generic" is localCoords,
// which very magically works in both cases.
void
FE_InterpolateCartLin
::
compute_dNdr
(
const
VECTOR
&
point
,
const
int
numNodes
,
const
int
nD
,
const
double
vol
,
DENS_MAT
&
dNdr
)
{
// *** see comments for compute_N_dNdr ***
const
DENS_MAT
&
localCoords
=
feElement_
->
local_coords
();
double
invVol
=
1.0
/
vol
;
for
(
int
inode
=
0
;
inode
<
numNodes
;
++
inode
)
{
for
(
int
idr
=
0
;
idr
<
nD
;
++
idr
)
{
dNdr
(
idr
,
inode
)
=
invVol
;
}
for
(
int
id
=
0
;
id
<
nD
;
++
id
)
{
for
(
int
idr
=
0
;
idr
<
nD
;
++
idr
)
{
if
(
id
==
idr
)
dNdr
(
idr
,
inode
)
*=
localCoords
(
id
,
inode
);
else
dNdr
(
idr
,
inode
)
*=
1.0
+
point
(
id
)
*
localCoords
(
id
,
inode
);
}
}
}
}
void
FE_InterpolateCartLin
::
compute_N_dNdr
(
const
VECTOR
&
point
,
VECTOR
&
N
,
DENS_MAT
&
dNdr
)
const
{
// Required data from element class
const
DENS_MAT
&
localCoords
=
feElement_
->
local_coords
();
double
invVol
=
1.0
/
(
feElement_
->
vol
());
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
// Fill in for each node
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
// Intiialize shape function and derivatives
N
(
inode
)
=
invVol
;
for
(
int
idr
=
0
;
idr
<
nSD_
;
++
idr
)
{
dNdr
(
idr
,
inode
)
=
invVol
;
}
for
(
int
isd
=
0
;
isd
<
nSD_
;
++
isd
)
{
// One term for each dimension
N
(
inode
)
*=
(
1.0
+
point
(
isd
)
*
localCoords
(
isd
,
inode
));
// One term for each dimension, only deriv in deriv's dimension
for
(
int
idr
=
0
;
idr
<
nSD_
;
++
idr
)
{
if
(
isd
==
idr
)
dNdr
(
idr
,
inode
)
*=
localCoords
(
isd
,
inode
);
else
dNdr
(
idr
,
inode
)
*=
1.0
+
point
(
isd
)
*
localCoords
(
isd
,
inode
);
}
}
}
}
/*********************************************************
* Class FE_InterpolateCartSerendipity
*
* For computing shape functions for quadratic serendipity
* elements, implemented for SPEED.
*
*********************************************************/
FE_InterpolateCartSerendipity
::
FE_InterpolateCartSerendipity
(
FE_Element
*
feElement
)
:
FE_Interpolate
(
feElement
)
{
set_quadrature
(
HEXA
,
GAUSS2
);
}
FE_InterpolateCartSerendipity
::~
FE_InterpolateCartSerendipity
()
{
// Handled by base class
}
void
FE_InterpolateCartSerendipity
::
compute_N
(
const
VECTOR
&
point
,
VECTOR
&
N
)
{
// *** see comments for compute_N_dNdr ***
const
DENS_MAT
&
localCoords
=
feElement_
->
local_coords
();
double
invVol
=
1.0
/
(
feElement_
->
vol
());
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
N
(
inode
)
=
invVol
;
for
(
int
isd
=
0
;
isd
<
nSD_
;
++
isd
)
{
if
(((
inode
==
8
||
inode
==
10
||
inode
==
16
||
inode
==
18
)
&&
(
isd
==
0
))
||
((
inode
==
9
||
inode
==
11
||
inode
==
17
||
inode
==
19
)
&&
(
isd
==
1
))
||
((
inode
==
12
||
inode
==
13
||
inode
==
14
||
inode
==
15
)
&&
(
isd
==
2
)))
{
N
(
inode
)
*=
(
1.0
-
pow
(
point
(
isd
),
2
))
*
2
;
}
else
{
N
(
inode
)
*=
(
1.0
+
point
(
isd
)
*
localCoords
(
isd
,
inode
));
}
}
if
(
inode
<
8
)
{
N
(
inode
)
*=
(
point
(
0
)
*
localCoords
(
0
,
inode
)
+
point
(
1
)
*
localCoords
(
1
,
inode
)
+
point
(
2
)
*
localCoords
(
2
,
inode
)
-
2
);
}
}
}
// Sort of a test-ride for a generic version that can be used for
// faces too. The only thing that's not "generic" is localCoords,
// which very magically works in both cases.
void
FE_InterpolateCartSerendipity
::
compute_dNdr
(
const
VECTOR
&
point
,
const
int
numNodes
,
const
int
nD
,
const
double
vol
,
DENS_MAT
&
dNdr
)
{
// *** see comments for compute_N_dNdr ***
const
DENS_MAT
&
localCoords
=
feElement_
->
local_coords
();
double
invVol
=
1.0
/
vol
;
bool
serendipityNode
=
false
;
double
productRule1
=
0.0
;
double
productRule2
=
0.0
;
if
(
nD
!=
3
&&
nD
!=
2
)
{
ATC_Error
(
"Serendipity dNdr calculations are too hard-wired to do "
"what you want them to. Only 2D and 3D currently work."
);
}
for
(
int
inode
=
0
;
inode
<
numNodes
;
++
inode
)
{
for
(
int
idr
=
0
;
idr
<
nD
;
++
idr
)
{
dNdr
(
idr
,
inode
)
=
invVol
;
}
for
(
int
id
=
0
;
id
<
nD
;
++
id
)
{
for
(
int
idr
=
0
;
idr
<
nD
;
++
idr
)
{
// identify nodes/dims differently if 3d or 2d case
if
(
nD
==
3
)
{
serendipityNode
=
(((
inode
==
8
||
inode
==
10
||
inode
==
16
||
inode
==
18
)
&&
(
id
==
0
))
||
((
inode
==
9
||
inode
==
11
||
inode
==
17
||
inode
==
19
)
&&
(
id
==
1
))
||
((
inode
==
12
||
inode
==
13
||
inode
==
14
||
inode
==
15
)
&&
(
id
==
2
)));
}
else
if
(
nD
==
2
)
{
serendipityNode
=
(((
inode
==
4
||
inode
==
6
)
&&
(
id
==
0
))
||
((
inode
==
5
||
inode
==
7
)
&&
(
id
==
1
)));
}
if
(
serendipityNode
)
{
if
(
id
==
idr
)
{
dNdr
(
idr
,
inode
)
*=
point
(
id
)
*
(
-
4
);
}
else
{
dNdr
(
idr
,
inode
)
*=
(
1.0
-
pow
(
point
(
id
),
2
))
*
2
;
}
}
else
{
if
(
id
==
idr
)
{
dNdr
(
idr
,
inode
)
*=
localCoords
(
id
,
inode
);
}
else
{
dNdr
(
idr
,
inode
)
*=
(
1.0
+
point
(
id
)
*
localCoords
(
id
,
inode
));
}
}
}
}
for
(
int
idr
=
0
;
idr
<
nD
;
++
idr
)
{
if
(
inode
<
8
)
{
// final corner contribution slightly different for 3d and 2d cases
if
(
nD
==
3
)
{
productRule2
=
(
point
(
0
)
*
localCoords
(
0
,
inode
)
+
point
(
1
)
*
localCoords
(
1
,
inode
)
+
point
(
2
)
*
localCoords
(
2
,
inode
)
-
2
);
}
else
if
(
nD
==
2
)
{
productRule2
=
(
point
(
0
)
*
localCoords
(
0
,
inode
)
+
point
(
1
)
*
localCoords
(
1
,
inode
)
-
1
);
}
productRule1
=
dNdr
(
idr
,
inode
)
*
(
1
+
point
(
idr
)
*
localCoords
(
idr
,
inode
));
productRule2
*=
dNdr
(
idr
,
inode
);
dNdr
(
idr
,
inode
)
=
productRule1
+
productRule2
;
}
}
}
}
void
FE_InterpolateCartSerendipity
::
compute_N_dNdr
(
const
VECTOR
&
point
,
VECTOR
&
N
,
DENS_MAT
&
dNdr
)
const
{
// Required data from element class
const
DENS_MAT
&
localCoords
=
feElement_
->
local_coords
();
double
invVol
=
1.0
/
(
feElement_
->
vol
());
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
// Will store terms for product rule derivative for dNdr
double
productRule1
;
double
productRule2
;
// Fill in for each node
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
// Initialize shape functions and derivatives
N
(
inode
)
=
invVol
;
for
(
int
idr
=
0
;
idr
<
nSD_
;
++
idr
)
{
dNdr
(
idr
,
inode
)
=
invVol
;
}
// Add components from each dimension
for
(
int
isd
=
0
;
isd
<
nSD_
;
++
isd
)
{
for
(
int
idr
=
0
;
idr
<
nSD_
;
++
idr
)
{
// Check to see if the node is NOT a corner node, and if its
// "0-coordinate" is in the same dimension as the one we're currently
// iterating over. If that's the case, we want to contribute to its
// shape functions and derivatives in a modified way:
if
(((
inode
==
8
||
inode
==
10
||
inode
==
16
||
inode
==
18
)
&&
(
isd
==
0
))
||
((
inode
==
9
||
inode
==
11
||
inode
==
17
||
inode
==
19
)
&&
(
isd
==
1
))
||
((
inode
==
12
||
inode
==
13
||
inode
==
14
||
inode
==
15
)
&&
(
isd
==
2
)))
{
// If the 1d shape function dimension matches the derivative
// dimension...
if
(
isd
==
idr
)
{
// contribute to N; sloppy, but this is the easiest way to get
// N to work right without adding extra, arguably unnecessary
// loops, while also computing the shape functions
N
(
inode
)
*=
(
1.0
-
pow
(
point
(
isd
),
2
))
*
2
;
// contribute to dNdr with the derivative of this shape function
// contribution
dNdr
(
idr
,
inode
)
*=
point
(
isd
)
*
(
-
4
);
}
else
{
// otherwise, just use the "constant" contribution to the deriv
dNdr
(
idr
,
inode
)
*=
(
1.0
-
pow
(
point
(
isd
),
2
))
*
2
;
}
}
else
{
// non-serendipity style contributions
if
(
isd
==
idr
)
{
N
(
inode
)
*=
(
1.0
+
point
(
isd
)
*
localCoords
(
isd
,
inode
));
dNdr
(
idr
,
inode
)
*=
localCoords
(
isd
,
inode
);
}
else
{
dNdr
(
idr
,
inode
)
*=
(
1.0
+
point
(
isd
)
*
localCoords
(
isd
,
inode
));
}
}
}
}
// serendipity corner nodes require more extra handling
if
(
inode
<
8
)
{
N
(
inode
)
*=
(
point
(
0
)
*
localCoords
(
0
,
inode
)
+
point
(
1
)
*
localCoords
(
1
,
inode
)
+
point
(
2
)
*
localCoords
(
2
,
inode
)
-
2
);
}
for
(
int
idr
=
0
;
idr
<
nSD_
;
++
idr
)
{
if
(
inode
<
8
)
{
productRule1
=
dNdr
(
idr
,
inode
)
*
(
1
+
point
(
idr
)
*
localCoords
(
idr
,
inode
));
productRule2
=
dNdr
(
idr
,
inode
)
*
(
point
(
0
)
*
localCoords
(
0
,
inode
)
+
point
(
1
)
*
localCoords
(
1
,
inode
)
+
point
(
2
)
*
localCoords
(
2
,
inode
)
-
2
);
dNdr
(
idr
,
inode
)
=
productRule1
+
productRule2
;
}
}
}
}
/*********************************************************
* Class FE_InterpolateSimpLin
*
* For computing linear shape functions of simplices,
* which are rather different from those computed
* in Cartesian coordinates.
*
* Note: degenerating quads/hexes can yield simplices
* as well, but this class is for computing these
* shape fucntions _natively_, in their own
* triangular/tetrahedral coordinate systems.
*
*********************************************************/
FE_InterpolateSimpLin
::
FE_InterpolateSimpLin
(
FE_Element
*
feElement
)
:
FE_Interpolate
(
feElement
)
{
set_quadrature
(
TETRA
,
GAUSS2
);
}
FE_InterpolateSimpLin
::~
FE_InterpolateSimpLin
()
{
// Handled by base class
}
void
FE_InterpolateSimpLin
::
compute_N
(
const
VECTOR
&
point
,
VECTOR
&
N
)
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
// Fill in for each node
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
if
(
inode
==
0
)
{
// Fill N...the ips are serving as proxies for "dimensions"
// since we're in tetrahedral coordinates, except that
// 0th node = 3rd "dimension" (u or O_o)
// 1st node = 0th "dimension" (x or r)
// 2nd node = 1st "dimension" (y or s)
// 3rd node = 3nd "dimension" (z or t)
// and remember that u = 1 - r - s - t for tet coords
N
(
inode
)
=
1
;
for
(
int
icont
=
0
;
icont
<
nSD_
;
++
icont
)
{
N
(
inode
)
-=
point
(
icont
);
}
}
else
{
N
(
inode
)
=
point
(
inode
-
1
);
}
}
}
void
FE_InterpolateSimpLin
::
compute_dNdr
(
const
VECTOR
&
,
const
int
numNodes
,
const
int
nD
,
const
double
,
DENS_MAT
&
dNdr
)
{
// Fill in for each node
for
(
int
inode
=
0
;
inode
<
numNodes
;
++
inode
)
{
// Fill dNdr_; we want 1 if the dimension of derivative
// and variable within N correspond. That is, if N == r,
// we want the 0th dimension to contain (d/dr)r = 1. Of
// course, (d/di)r = 0 forall i != r, so we need that as
// well. This is a bit elusively complicated. Also, the 0th
// integration point contains the term u = 1 - r - s - t.
// (which map to x, y, and z). Therefore, the derivative in
// each dimension are -1.
//
// The idea is similar for 2 dimensions, which this can
// handle as well.
for
(
int
idr
=
0
;
idr
<
nD
;
++
idr
)
{
if
(
inode
==
0
)
{
dNdr
(
idr
,
inode
)
=
-
1
;
}
else
{
dNdr
(
idr
,
inode
)
=
(
inode
==
(
idr
+
1
))
?
1
:
0
;
}
}
}
}
void
FE_InterpolateSimpLin
::
compute_N_dNdr
(
const
VECTOR
&
point
,
VECTOR
&
N
,
DENS_MAT
&
dNdr
)
const
{
int
numEltNodes
=
feElement_
->
num_elt_nodes
();
// Fill in for each node
for
(
int
inode
=
0
;
inode
<
numEltNodes
;
++
inode
)
{
// Fill N...
if
(
inode
==
0
)
{
N
(
inode
)
=
1
;
for
(
int
icont
=
0
;
icont
<
nSD_
;
++
icont
)
{
N
(
inode
)
-=
point
(
icont
);
}
}
else
{
N
(
inode
)
=
point
(
inode
-
1
);
}
// Fill dNdr...
for
(
int
idr
=
0
;
idr
<
nSD_
;
++
idr
)
{
if
(
inode
==
0
)
{
dNdr
(
idr
,
inode
)
=
-
1
;
}
else
{
dNdr
(
idr
,
inode
)
=
(
inode
==
(
idr
+
1
))
?
1
:
0
;
}
}
}
}
}
// namespace ATC
Event Timeline
Log In to Comment