OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added: Framework for co-rotational beam elements

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2714 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2014-02-13 14:58:05 +00:00
committed by Knut Morten Okstad
parent 078db406a6
commit e202ea3478
12 changed files with 386 additions and 177 deletions
Download Patch File Download Diff File Expand all files Collapse all files

184
src/ASM/ASMbase.C
View File

@@ -270,6 +270,9 @@ public:
case 1: return bc.CX == 0;
case 2: return bc.CY == 0;
case 3: return bc.CZ == 0;
case 4: return bc.RX == 0;
case 5: return bc.RY == 0;
case 6: return bc.RZ == 0;
}
return false;
}
@@ -353,15 +356,6 @@ bool ASMbase::add3PC (int slave, int dir, int master1, int master2, int code)
}
bool ASMbase::addSPC (int node, int dir, int code)
{
if (dir < 1 || dir > nf) return true;
MPC* mpc = new MPC(node,dir);
return this->addMPC(mpc,code);
}
MPC* ASMbase::findMPC (int node, int dof) const
{
MPC slave(node,dof);
@@ -394,41 +388,12 @@ bool ASMbase::addPeriodicity (size_t master, size_t slave, int dir)
void ASMbase::makePeriodic (size_t master, size_t slave, int dirs)
{
switch (dirs)
{
case 1:
case 2:
case 3:
this->addPeriodicity(master,slave,dirs);
break;
case 12:
case 21:
for (int dir = 1; dir <= 2; dir++)
this->addPeriodicity(master,slave,dir);
break;
case 13:
case 31:
for (int dir = 1; dir <= 3; dir += 2)
this->addPeriodicity(master,slave,dir);
break;
case 23:
case 32:
for (int dir = 2; dir <= 3; dir++)
this->addPeriodicity(master,slave,dir);
break;
default:
std::cerr <<" ** ASMbase::makePeriodic: Invalid DOF code "<< dirs
<<" replaced by 123"<< std::endl;
case 0:
case 123:
case 132:
case 213:
case 231:
case 312:
case 321:
// If all DOFs are going to be coupled, assign a common global node number
ASMbase::collapseNodes(*this,master,*this,slave);
}
std::set<int> dofs(utl::getDigits(dirs));
if (dofs.size() == nf && *dofs.rbegin() == nf)
// If all DOFs are going to be coupled, assign a common global node number
ASMbase::collapseNodes(*this,master,*this,slave);
else for (std::set<int>::iterator it = dofs.begin(); it != dofs.end(); ++it)
this->addPeriodicity(master,slave,*it);
}
@@ -439,48 +404,43 @@ void ASMbase::constrainPatch (int dof, int code)
}
void ASMbase::prescribe (size_t inod, int dirs, int code)
int ASMbase::prescribe (size_t inod, int dirs, int code)
{
if (code == 0 && fixHomogeneousDirichlet)
return this->fix(inod,dirs);
int node = this->getNodeID(inod);
if (node < 1) return;
if (node < 1 || dirs < 1) return dirs;
switch (dirs)
int ignoredDirs = 0;
std::set<int> dofs(utl::getDigits(dirs));
for (std::set<int>::const_iterator it = dofs.begin(); it != dofs.end(); ++it)
if (*it <= nf)
{
case 1:
case 2:
case 3:
this->addSPC(node,dirs,code);
break;
case 12:
case 21:
for (int dir = 1; dir <= 2; dir++)
this->addSPC(node,dir,code);
break;
case 13:
case 31:
for (int dir = 1; dir <= 3; dir += 2)
this->addSPC(node,dir,code);
break;
case 23:
case 32:
for (int dir = 2; dir <= 3; dir++)
this->addSPC(node,dir,code);
break;
default:
std::cerr <<" ** ASMbase::prescribe: Invalid DOF code "<< dirs
<<" replaced by 123"<< std::endl;
case 123:
case 132:
case 213:
case 231:
case 312:
case 321:
for (int dir = 1; dir <= 3; dir++)
this->addSPC(node,dir,code);
MPC* mpc = new MPC(node,*it);
if (!this->addMPC(mpc,code))
ignoredDirs = 10*ignoredDirs + *it;
}
else
{
ignoredDirs = 10*ignoredDirs + *it;
std::cerr <<" ** ASMbase::prescribe: Ignoring invalid DOF code "<< *it
<< std::endl;
}
return ignoredDirs;
}
/*!
\brief Equality operator for BC objects comparing node numbers and DOF codes.
*/
bool operator== (const ASMbase::BC& rhs, const ASMbase::BC& lhs)
{
return (rhs.node == lhs.node &&
rhs.CX == lhs.CX && rhs.CY == lhs.CY && rhs.CZ == lhs.CZ &&
rhs.RX == lhs.RX && rhs.RY == lhs.RY && rhs.RZ == lhs.RZ);
}
@@ -494,10 +454,10 @@ bool operator== (const ASMbase::BC& rhs, const int& lhs)
}
void ASMbase::fix (size_t inod, int dirs)
int ASMbase::fix (size_t inod, int dirs)
{
int node = this->getNodeID(inod);
if (node < 1) return;
if (node < 1 || dirs < 1) return dirs;
BCVec::iterator bit = std::find(BCode.begin(),BCode.end(),node);
if (bit == BCode.end())
@@ -506,52 +466,34 @@ void ASMbase::fix (size_t inod, int dirs)
bit = BCode.end()-1;
}
switch (dirs)
#if SP_DEBUG > 1
BC old = *bit;
#endif
int invalidDOFs = 0;
std::set<int> dofs(utl::getDigits(dirs));
for (std::set<int>::const_iterator it = dofs.begin(); it != dofs.end(); ++it)
if (*it <= nf)
switch (*it) {
case 1: bit->CX = 0; break;
case 2: bit->CY = 0; break;
case 3: bit->CZ = 0; break;
case 4: bit->RX = 0; break;
case 5: bit->RY = 0; break;
case 6: bit->RZ = 0; break;
}
else
{
case 1:
if (bit->CX == 0) return;
bit->CX = 0;
break;
case 2:
if (bit->CY == 0) return;
bit->CY = 0;
break;
case 3:
if (bit->CZ == 0) return;
bit->CZ = 0;
break;
case 12:
case 21:
if (bit->CX + bit->CY == 0) return;
bit->CX = bit->CY = 0;
break;
case 13:
case 31:
if (bit->CX + bit->CZ == 0) return;
bit->CX = bit->CZ = 0;
break;
case 23:
case 32:
if (bit->CY + bit->CZ == 0) return;
bit->CY = bit->CZ = 0;
break;
default:
std::cerr <<" ** ASMbase::fix: Invalid DOF code "<< dirs
<<" replaced by 123"<< std::endl;
dirs = 123;
case 123:
case 132:
case 213:
case 231:
case 312:
case 321:
if (bit->CX + bit->CY + bit->CZ == 0) return;
bit->CX = bit->CY = bit->CZ = 0;
invalidDOFs = 10*invalidDOFs + *it;
std::cerr <<" ** ASMbase::fix: Ignoring invalid DOF code "<< *it
<< std::endl;
}
#if SP_DEBUG > 1
std::cout <<"\tFixed node: "<< node <<" "<< dirs << std::endl;
if (!(old == *bit))
std::cout <<"\tFixed node: "<< node <<" "<< dirs << std::endl;
#endif
return invalidDOFs;
}

27
src/ASM/ASMbase.h
View File

@@ -60,8 +60,11 @@ public:
char CX; //!< Boundary condition code for X-translation
char CY; //!< Boundary condition code for Y-translation
char CZ; //!< Boundary condition code for Z-translation
char RX; //!< Boundary condition code for X-rotation
char RY; //!< Boundary condition code for Y-rotation
char RZ; //!< Boundary condition code for Z-rotation
//! \brief Constructor initializing a BC instance.
BC(int n) : node(n), CX(1), CY(1), CZ(1) {}
BC(int n) : node(n), CX(1), CY(1), CZ(1), RX(1), RY(1), RZ(1) {}
};
typedef std::vector<BC> BCVec; //!< Nodal boundary condition container
@@ -127,6 +130,11 @@ public:
//! \param[in] ng Number of Gauss points in each parameter direction
void setGauss(int ng) { nGauss = ng; }
//! \brief Defines the number of solution fields \a nf in the patch.
//! \brief This method is to be used by simulators where \a nf is not known
//! when the patch is constructed, e.g., it depends on the input file content.
//! It must be invoked only before the SIMbase::preprocess is invoked.
void setNoFields(unsigned char n) { nf = n; }
// Service methods for query of various model data
// ===============================================
@@ -223,9 +231,9 @@ public:
// Various preprocessing methods
// =============================
//! \brief Copy the parameter domain from another ASM
//! \param other The ASM to copy parameter domain from
//! \brief Copies the parameter domain from another patch.
//! \param[in] other The patch to copy parameter domain from
virtual void copyParameterDomain(const ASMbase* other) {}
//! \brief Merges a given node in this patch with a given global node.
@@ -528,11 +536,6 @@ protected:
//! \param[in] master2 Global node number of 2nd master node of the constraint
//! \param[in] code Identifier for inhomogeneous Dirichlet condition field
bool add3PC(int slave, int dir, int master1, int master2, int code = 0);
//! \brief Creates and adds a single-point constraint to this patch.
//! \param[in] node Global node number of the node to constrain
//! \param[in] dir Which local DOF to constrain (1, 2, 3)
//! \param[in] code Identifier for inhomogeneous Dirichlet condition field
bool addSPC(int node, int dir, int code);
//! \brief Creates and adds a periodicity constraint to this patch.
//! \param[in] master 1-based local index of the master node
//! \param[in] slave 1-based local index of the slave node to constrain
@@ -556,11 +559,13 @@ public:
//! \param[in] inod 1-based node index local to current patch
//! \param[in] dirs Which local DOFs to constrain (1, 2, 3, 12, 23, 123)
//! \param[in] code Identifier for inhomogeneous Dirichlet condition field
void prescribe(size_t inod, int dirs, int code);
//! \return Invalid local DOFs, or DOFs that already are constrained
int prescribe(size_t inod, int dirs, int code);
//! \brief Constrains DOFs in the given node to zero.
//! \param[in] inod 1-based node index local to current patch
//! \param[in] dirs Which local DOFs to constrain (1, 2, 3, 12, 23, 123)
void fix(size_t inod, int dirs = 123);
//! \return Invalid local DOFs
int fix(size_t inod, int dirs = 123);
//! \brief Checks whether given DOFs are fixed or not.
//! \param[in] node Global node number of the DOF to check
//! \param[in] dof Local indices of the DOFs to check

101
src/ASM/ASMs1D.C
View File

@@ -31,13 +31,13 @@
ASMs1D::ASMs1D (unsigned char n_s, unsigned char n_f)
: ASMstruct(1,n_s,n_f), curv(0)
: ASMstruct(1,n_s,n_f), curv(0), nodalT(myT), elmCS(myCS)
{
}
ASMs1D::ASMs1D (const ASMs1D& patch, unsigned char n_f)
: ASMstruct(patch,n_f), curv(patch.curv)
: ASMstruct(patch,n_f), curv(patch.curv), nodalT(patch.myT), elmCS(patch.myCS)
{
}
@@ -211,6 +211,12 @@ bool ASMs1D::generateFEMTopology ()
myMLGE.resize(n1-p1+1,0);
myMLGN.resize(n1);
myMNPC.resize(myMLGE.size());
if (nsd == 3 && nf == 6)
{
// This is a 3D beam problem, allocation rotation tensors
myT.resize(n1,Tensor(3,true));
myCS.resize(n1-p1+1,Tensor(3));
}
int iel = 0;
int inod = 0;
@@ -236,6 +242,15 @@ bool ASMs1D::generateFEMTopology ()
#ifdef SP_DEBUG
std::cout <<"NEL = "<< iel <<" NNOD = "<< inod << std::endl;
#endif
// Calculate local element axes for 3D beam elements
for (size_t i = 0; i < myCS.size(); i++)
{
Vec3 X1 = this->getCoord(1+MNPC[i].back());
Vec3 X2 = this->getCoord(1+MNPC[i].front());
myCS[i] = Tensor(X2-X1).shift();
}
return true;
}
@@ -414,6 +429,32 @@ bool ASMs1D::getElementCoordinates (Matrix& X, int iel) const
}
bool ASMs1D::getElementNodalRotations (TensorVec& T, int iel) const
{
#ifdef INDEX_CHECK
if (iel < 1 || (size_t)iel > MNPC.size())
{
std::cerr <<" *** ASMs1D::getElementNodalRotations: Element index "<< iel
<<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
return false;
}
#endif
T.clear();
if (nodalT.empty())
return true;
const IntVec& mnpc = MNPC[iel-1];
Tensor Tgl(elmCS[iel-1],true);
T.reserve(mnpc.size());
for (size_t i = 0; i < mnpc.size(); i++)
T.push_back(Tgl*nodalT[mnpc[i]]);
return true;
}
void ASMs1D::getNodalCoordinates (Matrix& X) const
{
const int n1 = curv->numCoefs();
@@ -448,6 +489,25 @@ bool ASMs1D::updateCoords (const Vector& displ)
}
bool ASMs1D::updateRotations (const Vector& displ)
{
if (shareFE || nf != 6) return true;
if (displ.size() != 6*myT.size())
{
std::cerr <<" *** ASMs1D::updateRotations: Invalid dimension "
<< displ.size() <<" on displ, should be "
<< 6*myT.size() << std::endl;
return false;
}
for (size_t i = 0; i < myT.size(); i++)
myT[i] *= Tensor(displ[6*i+3],displ[6*i+4],displ[6*i+5]);
return true;
}
void ASMs1D::getBoundaryNodes (int lIndex, IntVec& glbNodes) const
{
if (!curv) return; // silently ignore empty patches
@@ -563,9 +623,9 @@ void ASMs1D::extractBasis (double u, Vector& N, Matrix& dNdu,
d2Ndu2.resize(p1,1,1);
for (int i = 1; i <= p1; i++)
{
N(i) = bas[3*i-3];
dNdu(i,1) = bas[3*i-2];
d2Ndu2(i,1,1) = bas[3*i-1];
N(i) = bas[3*i-3];
dNdu(i,1) = bas[3*i-2];
d2Ndu2(i,1,1) = bas[3*i-1];
}
}
@@ -622,6 +682,9 @@ bool ASMs1D::integrate (Integrand& integrand,
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementEnds(i1-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::NODAL_ROTATIONS)
this->getElementNodalRotations(fe.T,iel);
// Initialize element matrices
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel);
bool ok = integrand.initElement(MNPC[iel-1],fe,X,nRed,*A);
@@ -673,11 +736,11 @@ bool ASMs1D::integrate (Integrand& integrand,
// Compute basis functions and derivatives
if (integrand.getIntegrandType() & Integrand::NO_DERIVATIVES)
curv->basis().computeBasisValues(fe.u,fe.N.ptr());
curv->basis().computeBasisValues(fe.u,fe.N.ptr());
else if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
this->extractBasis(fe.u,fe.N,dNdu,d2Ndu2);
this->extractBasis(fe.u,fe.N,dNdu,d2Ndu2);
else
this->extractBasis(fe.u,fe.N,dNdu);
this->extractBasis(fe.u,fe.N,dNdu);
if (!dNdu.empty())
{
@@ -698,7 +761,7 @@ bool ASMs1D::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dL*wg[i];
if (ok && !integrand.evalInt(*A,fe,time,X))
ok = false;
ok = false;
}
// Finalize the element quantities
@@ -950,11 +1013,11 @@ bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
for (size_t i = 0; i < nPoints; i++)
{
IntVec ip;
scatterInd(p1,curv->basis().lastKnotInterval(),ip);
switch (deriv) {
case 0: // Evaluate the solution
curv->basis().computeBasisValues(upar[i],&basis.front());
scatterInd(p1,curv->basis().lastKnotInterval(),ip);
utl::gather(ip,nComp,locSol,Xtmp);
Xtmp.multiply(basis,ptSol);
sField.fillColumn(1+i,ptSol);
@@ -962,6 +1025,7 @@ bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
case 1: // Evaluate first derivatives of the solution
this->extractBasis(upar[i],basis,dNdu);
scatterInd(p1,curv->basis().lastKnotInterval(),ip);
utl::gather(ip,nsd,Xnod,Xtmp);
utl::Jacobian(Jac,dNdX,Xtmp,dNdu);
utl::gather(ip,nComp,locSol,Xtmp);
@@ -971,6 +1035,7 @@ bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
case 2: // Evaluate second derivatives of the solution
this->extractBasis(upar[i],basis,dNdu,d2Ndu2);
scatterInd(p1,curv->basis().lastKnotInterval(),ip);
utl::gather(ip,nsd,Xnod,Xtmp);
utl::Jacobian(Jac,dNdX,Xtmp,dNdu);
utl::Hessian(Hess,d2NdX2,Jac,Xtmp,d2Ndu2,dNdX);
@@ -1132,15 +1197,23 @@ bool ASMs1D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
bool ASMs1D::globalL2projection (Matrix& sField,
const IntegrandBase& integrand) const
const IntegrandBase& integrand,
bool continuous) const
{
if (!curv) return true; // silently ignore empty patches
if (continuous)
{
std::cerr <<" *** ASMs1D::globalL2projection: Only discrete L2-projection"
<<" is available in this method.\n "
<<" Use ASMbase::L2projection instead."<< std::endl;
}
const int p1 = curv->order();
const int n1 = curv->numCoefs();
const int nel = n1 - p1 + 1;
// Get Gaussian quadrature points
// Get Gaussian quadrature point coordinates
const int ng1 = p1 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
if (!xg) return false;
@@ -1205,9 +1278,9 @@ bool ASMs1D::globalL2projection (Matrix& sField,
}
bool ASMs1D::getNoStructElms(int& n1, int& n2, int& n3) const
bool ASMs1D::getNoStructElms (int& n1, int& n2, int& n3) const
{
n1 = getNoElms();
n1 = curv->numCoefs() - curv->order() + 1;
n2 = n3 = 0;
return true;

33
src/ASM/ASMs1D.h
View File

@@ -15,6 +15,9 @@
#define _ASM_S1D_H
#include "ASMstruct.h"
#include "Tensor.h"
typedef std::vector<Tensor> TensorVec; //!< An array of non-symmetric tensors
namespace Go {
class SplineCurve;
@@ -77,6 +80,10 @@ public:
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Updates the nodal rotations for this patch.
//! \param[in] displ Incremental displacements to update the rotations with
bool updateRotations(const Vector& displ);
//! \brief Finds the global number of the node on a patch end.
//! \param[in] lIndex Local index of the end point
//! \param glbNodes Array of global boundary node numbers
@@ -211,12 +218,13 @@ public:
//! \param[in] nSegSpan Number of visualization segments over each knot-span
virtual bool getGridParameters(RealArray& prm, int nSegSpan) const;
using ASMstruct::globalL2projection;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand) const;
const IntegrandBase& integrand,
bool continuous = false) const;
protected:
@@ -273,6 +281,11 @@ protected:
//! \param[out] XC Coordinates of the element corners
void getElementEnds(int i, std::vector<Vec3>& XC) const;
//! \brief Returns nodal rotation matrices for an element, if any.
//! \param[out] T Array of nodal rotation matrices
//! \param[in] iel Element index
bool getElementNodalRotations(TensorVec& T, int iel) const;
public:
//! \brief Auxilliary function for computation of basis function indices.
static void scatterInd(int p1, int start, IntVec& index);
@@ -285,19 +298,25 @@ public:
//! \brief Returns the number of nodal points in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[out] n2 Number of nodes in second (v) direction (always zero)
//! \param[out] n3 Number of nodes in third (w) direction (always zero)
//! \param[in] basis Which basis to return size parameters for (mixed methods)
virtual bool getSize(int& n1, int& n2, int& n3, int basis) const;
//! \brief Returns the number of elements in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[out] n1 Number of elements in first (u) direction
//! \param[out] n2 Number of elements in second (v) direction (always zero)
//! \param[out] n3 Number of elements in third (w) direction (always zero)
virtual bool getNoStructElms(int& n1, int& n2, int& n3) const;
protected:
Go::SplineCurve* curv; //!< Pointer to the actual spline curve object
const TensorVec& nodalT; //!< Nodal rotation tensors (for 3D beams)
const TensorVec& elmCS; //!< Element coordinate systems (for 3D beams)
TensorVec myT; //!< The actual nodal rotation tensors
TensorVec myCS; //!< The actual element coordinate systems
};
#endif

6
src/ASM/ASMstruct.h
View File

@@ -56,9 +56,9 @@ public:
virtual bool getSize(int& n1, int& n2, int& n3, int basis = 0) const = 0;
//! \brief Returns the number of elements in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[out] n1 Number of elements in first (u) direction
//! \param[out] n2 Number of elements in second (v) direction
//! \param[out] n3 Number of elements in third (w) direction
virtual bool getNoStructElms(int& n1, int& n2, int& n3) const = 0;
using ASMbase::evalSolution;

7
src/ASM/FiniteElement.h
View File

@@ -16,6 +16,7 @@
#include "MatVec.h"
#include "Vec3.h"
#include "Tensor.h"
/*!
@@ -26,7 +27,7 @@ class FiniteElement
{
public:
//! \brief Default constructor.
FiniteElement(size_t n = 0, size_t i = 0) : iel(0), iGP(i), N(n), detJxW(1.0)
FiniteElement(size_t n = 0, size_t i = 0) : iel(0), iGP(i), detJxW(1.0), N(n)
{ u = v = w = xi = eta = zeta = 0.0; }
//! \brief Empty destructor.
@@ -40,13 +41,15 @@ public:
double xi; //!< First local coordinate of current integration point
double eta; //!< Second local coordinate of current integration point
double zeta; //!< Third local coordinate of current integration point
double detJxW; //!< Weighted determinant of the coordinate mapping
Vector N; //!< Basis function values
Vector Navg; //!< Volume-averaged basis function values
Matrix dNdX; //!< First derivatives (gradient) of the basis functions
Matrix3D d2NdX2; //!< Second derivatives of the basis functions
Matrix G; //!< Matrix used for stabilized methods
Vec3Vec XC; //!< Matrix with element corner coordinates
double detJxW; //!< Weighted determinant of the coordinate mapping
std::vector<Tensor> T; //!< Array of nodal rotation matrices
};

3
src/ASM/Integrand.h
View File

@@ -141,7 +141,8 @@ public:
ELEMENT_CORNERS = 8, //!< Integrand wants element corner coordinates
ELEMENT_CENTER = 16,//!< Integrand wants element center coordinates
G_MATRIX = 32,//!< Integrand wants the G matrix
XO_ELEMENTS = 64 //!< Integrand is defined on extraordinary elements
NODAL_ROTATIONS = 64,//!< Integrand wants nodal rotation tesnros
XO_ELEMENTS = 128//!< Integrand is defined on extraordinary elements
};
//! \brief Defines which FE quantities are needed by the integrand.

3
src/ASM/SAMpatch.C
View File

@@ -119,6 +119,9 @@ bool SAMpatch::initNodeDofs (const ASMVec& model)
if (nndof > 0) msc[idof1-1] *= bit->CX;
if (nndof > 1) msc[idof1 ] *= bit->CY;
if (nndof > 2) msc[idof1+1] *= bit->CZ;
if (nndof > 3) msc[idof1+2] *= bit->RX;
if (nndof > 4) msc[idof1+3] *= bit->RY;
if (nndof > 5) msc[idof1+4] *= bit->RZ;
}
if (ierr == 0) return true;

162
src/Utility/Tensor.C
View File

@@ -13,6 +13,12 @@
#include "Tensor.h"
#include "Vec3.h"
#include <algorithm>
#ifndef epsZ
//! \brief Zero tolerance for the incremental rotations.
#define epsZ 1.0e-16
#endif
std::ostream& Tensor::print (std::ostream& os) const
@@ -33,6 +39,16 @@ std::ostream& Tensor::print (std::ostream& os) const
}
Tensor::Tensor (const t_ind nsd, bool identity) : n(nsd)
{
v.resize(n*n,Real(0));
if (identity)
for (t_ind i = 0; i < n*n; i += n+1)
v[i] = Real(1);
}
/*!
The provided vector \a vn is taken as the local Z-axis.
The local X- and Y-axes are then defined by projecting either the global X-
@@ -101,10 +117,43 @@ Tensor::Tensor (const Vec3& v1, const Vec3& v2, const Vec3& v3) : n(3)
}
Tensor::Tensor (const Tensor& T) : n(T.n)
/*!
This constructor computes an incremental rotation tensor from the given
rotation angles via a quaternion representation of the rotation.
The angles provided are those related to a Rodrigues parameterization.
*/
Tensor::Tensor (Real a1, Real a2, Real a3) : n(3)
{
Vec3 q(a1,a2,a3);
double theta = q.length();
q *= (theta < epsZ ? 0.5 : sin(0.5*theta)/theta);
Real q0 = cos(0.5*theta);
Real ql = sqrt(q0*q0 + q.x*q.x + q.y*q.y + q.z*q.z);
q0 /= ql;
q /= ql;
v[0] = 2.0*(q.x*q.x + q0*q0) - 1.0;
v[4] = 2.0*(q.y*q.y + q0*q0) - 1.0;
v[8] = 2.0*(q.z*q.z + q0*q0) - 1.0;
v[3] = 2.0*(q.x*q.y - q.z*q0);
v[6] = 2.0*(q.x*q.z + q.y*q0);
v[7] = 2.0*(q.y*q.z - q.x*q0);
v[1] = 2.0*(q.y*q.x + q.z*q0);
v[2] = 2.0*(q.z*q.x - q.y*q0);
v[4] = 2.0*(q.z*q.y + q.x*q0);
}
Tensor::Tensor (const Tensor& T, bool transpose) : n(T.n)
{
v.resize(n*n);
this->operator=(T);
std::copy(T.v.begin(),T.v.end(),v.begin());
if (transpose)
this->transpose();
}
@@ -120,6 +169,15 @@ void Tensor::define3Dtransform (const Vec3& v1, const Vec3& v2, const Vec3& v3)
}
Vec3 Tensor::operator[] (t_ind i) const
{
if (i >= n)
return Vec3();
return Vec3(&v[0]+n*i,n);
}
Tensor& Tensor::operator= (const Tensor& T)
{
if (v.size() == T.v.size())
@@ -238,6 +296,34 @@ Tensor& Tensor::operator*= (Real val)
}
Tensor& Tensor::operator*= (const Tensor& B)
{
switch (std::min(n,B.n)) {
case 1:
v[0] *= B.v[0];
break;
case 2:
{
Tensor A(*this);
for (int i = 1; i <= 2; i++)
for (int j = 1; j <= 2; j++)
v[this->index(i,j)] = A(i,1)*B(1,j) + A(i,2)*B(2,j);
break;
}
case 3:
{
Tensor A(*this);
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
v[this->index(i,j)] = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j);
break;
}
}
return *this;
}
Real Tensor::innerProd (const Tensor& T) const
{
Real value = Real(0);
@@ -302,6 +388,30 @@ Tensor& Tensor::symmetrize ()
}
Tensor& Tensor::shift (short int idx)
{
if (this->symmetric() || idx <= -n || (n+idx)%n == 0)
return *this;
if (n == 2)
{
std::swap(v[0],v[2]);
std::swap(v[1],v[3]);
}
else if (n == 3)
{
t_ind j = (3+idx)%3;
for (t_ind r = 0; r < 3; r++)
{
std::swap(v[r],v[r+3*j]);
std::swap(v[r],v[r-3*j+9]);
}
}
return *this;
}
Real Tensor::trace () const
{
if (n == 3)
@@ -378,7 +488,7 @@ Vec3 operator* (const Tensor& T, const Vec3& v)
{
switch (T.n) {
case 1:
return Vec3(T(1,1)*v.x, v.y, v.z);
return Vec3(T.v[0]*v.x, v.y, v.z);
case 2:
return Vec3(T(1,1)*v.x + T(1,2)*v.y,
T(2,1)*v.x + T(2,2)*v.y, v.z);
@@ -399,7 +509,7 @@ Vec3 operator* (const Vec3& v, const Tensor& T)
{
switch (T.n) {
case 1:
return Vec3(T(1,1)*v.x, v.y, v.z);
return Vec3(T.v[0]*v.x, v.y, v.z);
case 2:
return Vec3(T(1,1)*v.x + T(2,1)*v.y,
T(1,2)*v.x + T(2,2)*v.y, v.z);
@@ -412,6 +522,34 @@ Vec3 operator* (const Vec3& v, const Tensor& T)
}
/*!
\brief Multiplication between two Tensors.
*/
Tensor operator* (const Tensor& A, const Tensor& B)
{
Tensor C(std::min(A.n,B.n));
switch (C.n) {
case 1:
C.v[0] = A.v[0]*B.v[0];
break;
case 2:
for (int i = 1; i <= 2; i++)
for (int j = 1; j <= 2; j++)
C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j);
break;
case 3:
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j);
break;
}
return C;
}
std::ostream& SymmTensor::print (std::ostream& os) const
{
switch (n) {
@@ -703,10 +841,10 @@ void SymmTensor::principal (Vec3& p) const
Real b1 = this->trace() / Real(3);
Real s1 = v[0] - b1;
Real s2 = v[1] - b1;
Real s3 = v.size() > 3 ? v[2] - b1 : 0.0;
Real s3 = v.size() > 3 ? v[2] - b1 : Real(0);
Real s4 = n > 2 ? v[3] : v.back();
Real s5 = n > 2 ? v[4] : 0.0;
Real s6 = n > 2 ? v[5] : 0.0;
Real s5 = n > 2 ? v[4] : Real(0);
Real s6 = n > 2 ? v[5] : Real(0);
// Compute 2nd and 3rd invariants of deviator J_2 and J_3
@@ -752,7 +890,7 @@ SymmTensor operator+ (const SymmTensor& T, Real a)
{
SymmTensor S(T);
for (Tensor::t_ind i = 0; i < S.n; i++)
for (unsigned short int i = 0; i < S.n; i++)
S.v[i] += a;
if (S.v.size() == 4)
@@ -770,7 +908,7 @@ SymmTensor operator- (const SymmTensor& T, Real a)
{
SymmTensor S(T);
for (Tensor::t_ind i = 0; i < S.n; i++)
for (unsigned short int i = 0; i < S.n; i++)
S.v[i] -= a;
if (S.v.size() == 4)
@@ -788,7 +926,7 @@ SymmTensor operator* (Real a, const SymmTensor& T)
{
SymmTensor S(T.dim(), T.v.size() == 4);
for (Tensor::t_ind i = 0; i < T.v.size(); i++)
for (size_t i = 0; i < T.v.size(); i++)
S.v[i] = a*T.v[i];
return S;
@@ -815,13 +953,13 @@ SymmTensor4::SymmTensor4 (const std::vector<Real>& x, t_ind nsd)
const Real& SymmTensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l) const
{
return v[index(i,j)*m+index(k,l)];
return v[this->index(i,j)*m+this->index(k,l)];
}
Real& SymmTensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l)
{
return ptr[index(i,j)*m+index(k,l)];
return ptr[this->index(i,j)*m+this->index(k,l)];
}

22
src/Utility/Tensor.h
View File

@@ -26,9 +26,9 @@ class Vec3;
class Tensor
{
public:
typedef unsigned short int t_ind; //!< Tensor index type
protected:
typedef unsigned short int t_ind; //!< Tensor index type (for convenience)
const t_ind n; //!< Number of spatial dimensions for the tensor
std::vector<Real> v; //!< The actual tensor component values
@@ -44,16 +44,18 @@ private:
void define3Dtransform(const Vec3& v1, const Vec3& v2, const Vec3& v3);
public:
//! \brief Constructor creating a zero tensor.
Tensor(const t_ind nsd) : n(nsd) { v.resize(n*n,Real(0)); }
//! \brief Constructor creating a zero or identity tesnor.
Tensor(const t_ind nsd, bool identity = false);
//! \brief Constructor creating a transformation from a face normal vector.
Tensor(const Vec3& vn);
//! \brief Constructor creating a transformation from two tangent vectors.
Tensor(const Vec3& t1, const Vec3& t2, bool t1isZ = false);
//! \brief Constructor creating a transformation from three unit vectors.
Tensor(const Vec3& v1, const Vec3& v2, const Vec3& v3);
//! \brief Copy constructor.
Tensor(const Tensor& T);
//! \brief Constructor creating a transformation from three rotation angles.
Tensor(Real a1, Real a2, Real a3);
//! \brief Copy constructor, optionally creating the transpose of \b T.
Tensor(const Tensor& T, bool transpose = false);
//! \brief Sets \a this to the 0-tensor.
void zero() { std::fill(v.begin(),v.end(),Real(0)); }
@@ -70,6 +72,8 @@ public:
const Real& operator()(t_ind i, t_ind j) const { return v[this->index(i,j)]; }
//! \brief Index-1 based component access.
Real& operator()(t_ind i, t_ind j) { return v[this->index(i,j)]; }
//! \brief Index-0 based column reference.
Vec3 operator[](t_ind i) const;
//! \brief Assignment operator.
Tensor& operator=(const Tensor& T);
@@ -90,6 +94,8 @@ public:
//! \brief Scaling operator.
Tensor& operator*=(Real val);
//! \brief Post-multiplication with another Tensor.
Tensor& operator*=(const Tensor& B);
//! \brief Returns the inner-product of \a *this and the given tensor.
Real innerProd(const Tensor& T) const;
@@ -110,6 +116,8 @@ public:
virtual Tensor& transpose();
//! \brief Makes the tensor symmetric.
virtual Tensor& symmetrize();
//! \brief Performs a cyclic permutation of the tensor columns.
Tensor& shift(short int idx = 1);
//! \brief Returns the trace of the tensor.
virtual Real trace() const;
@@ -128,6 +136,8 @@ public:
friend Vec3 operator*(const Tensor& T, const Vec3& v);
//! \brief Multiplication between a point vector and transpose of a tensor.
friend Vec3 operator*(const Vec3& v, const Tensor& T);
//! \brief Multiplication between two tensors.
friend Tensor operator*(const Tensor& A, const Tensor& B);
//! \brief Output stream operator.
friend std::ostream& operator<<(std::ostream& os, const Tensor& T)

11
src/Utility/Utilities.C
View File

@@ -388,3 +388,14 @@ std::string utl::adjustRight (size_t width, const std::string& s,
static std::string blank(32,' ');
return blank.substr(0,width-s.size()) + s + suffix;
}
std::set<int> utl::getDigits (int num)
{
std::set<int> digits;
for (; num != 0; num /= 10)
digits.insert(num%10);
return digits;
}

4
src/Utility/Utilities.h
View File

@@ -20,6 +20,7 @@
#include <iostream>
#include <vector>
#include <map>
#include <set>
class TiXmlElement;
class TiXmlNode;
@@ -179,6 +180,9 @@ namespace utl
//! \brief Right-justifies the input string to the given total \a width.
std::string adjustRight(size_t width, const std::string& s,
const std::string& suffix = " : ");
//! \brief Splits an integer into its (unique) digits in ascending order.
std::set<int> getDigits(int value);
}
#endif