Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77522686
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
Wed, Aug 14, 23:37
Size
34 KB
Mime Type
text/x-c
Expires
Fri, Aug 16, 23:37 (2 d)
Engine
blob
Format
Raw Data
Handle
19764438
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