#include <planej.h>

Inheritance diagram for AxialRotationJoint:

Collaboration diagram for AxialRotationJoint:

Public Member Functions
	AxialRotationJoint (unsigned int uL, const DofOwner pDO, const StructNode pN1, const StructNode pN2, const Vec3 &dTmp1, const Vec3 &dTmp2, const Mat3x3 &R1hTmp, const Mat3x3 &R2hTmp, const DriveCaller pDC, const OrientationDescription &od, flag fOut, const doublereal rr=0., const doublereal pref=0., BasicShapeCoefficient const sh=0, BasicFriction const f=0)

	~AxialRotationJoint (void)

virtual Joint::Type	GetJointType (void) const

virtual std::ostream &	Restart (std::ostream &out) const

virtual unsigned int	iGetNumDof (void) const

virtual std::ostream &	DescribeDof (std::ostream &out, const char *prefix="", bool bInitial=false) const

virtual void	DescribeDof (std::vector< std::string > &desc, bool bInitial=false, int i=-1) const

virtual std::ostream &	DescribeEq (std::ostream &out, const char *prefix="", bool bInitial=false) const

virtual void	DescribeEq (std::vector< std::string > &desc, bool bInitial=false, int i=-1) const

DofOrder::Order	GetDofType (unsigned int i) const

virtual void	SetValue (DataManager pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints ph=0)

virtual Hint *	ParseHint (DataManager pDM, const char s) const

virtual void	AfterConvergence (const VectorHandler &X, const VectorHandler &XP)

void	WorkSpaceDim (integer piNumRows, integer piNumCols) const

VariableSubMatrixHandler &	AssJac (VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)

SubVectorHandler &	AssRes (SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)

DofOrder::Order	GetEqType (unsigned int i) const

void	OutputPrepare (OutputHandler &OH)

void	Output (OutputHandler &OH) const

virtual unsigned int	iGetInitialNumDof (void) const

virtual void	InitialWorkSpaceDim (integer piNumRows, integer piNumCols) const

VariableSubMatrixHandler &	InitialAssJac (VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)

SubVectorHandler &	InitialAssRes (SubVectorHandler &WorkVec, const VectorHandler &XCurr)

virtual unsigned int	iGetNumPrivData (void) const

virtual unsigned int	iGetPrivDataIdx (const char *s) const

virtual doublereal	dGetPrivData (unsigned int i) const

virtual void	GetConnectedNodes (std::vector< const Node * > &connectedNodes) const

Public Member Functions inherited from Elem
	Elem (unsigned int uL, flag fOut)

virtual	~Elem (void)

virtual void	AssMats (VariableSubMatrixHandler &WorkMatA, VariableSubMatrixHandler &WorkMatB, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)

virtual bool	bInverseDynamics (void) const

void	SetInverseDynamicsFlags (unsigned uIDF)

unsigned	GetInverseDynamicsFlags (void) const

bool	bIsErgonomy (void) const

bool	bIsRightHandSide (void) const

virtual VariableSubMatrixHandler &	AssJac (VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)

virtual SubVectorHandler &	AssRes (SubVectorHandler &WorkVec, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr, const VectorHandler &XPrimePrimeCurr, InverseDynamics::Order iOrder=InverseDynamics::INVERSE_DYNAMICS)

virtual int	GetNumConnectedNodes (void) const

Public Member Functions inherited from WithLabel
	WithLabel (unsigned int uL=0, const std::string &sN="")

virtual	~WithLabel (void)

void	PutLabel (unsigned int uL)

void	PutName (const std::string &sN)

unsigned int	GetLabel (void) const

const std::string &	GetName (void) const

Public Member Functions inherited from SimulationEntity
	SimulationEntity (void)

virtual	~SimulationEntity (void)

virtual bool	bIsValidIndex (unsigned int i) const

virtual void	BeforePredict (VectorHandler &, VectorHandler &, VectorHandler &, VectorHandler &) const

virtual void	AfterPredict (VectorHandler &X, VectorHandler &XP)

virtual void	Update (const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)

virtual void	DerivativesUpdate (const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)

virtual void	AfterConvergence (const VectorHandler &X, const VectorHandler &XP, const VectorHandler &XPP)

virtual std::ostream &	OutputAppend (std::ostream &out) const

virtual void	ReadInitialState (MBDynParser &HP)

Public Member Functions inherited from ToBeOutput
	ToBeOutput (flag fOut=fDefaultOut)

virtual	~ToBeOutput (void)

virtual void	Output (OutputHandler &OH, const VectorHandler &X, const VectorHandler &XP) const

virtual flag	fToBeOutput (void) const

virtual bool	bToBeOutput (void) const

virtual void	SetOutputFlag (flag f=flag(1))

Public Member Functions inherited from Joint
	Joint (unsigned int uL, const DofOwner *pD, flag fOut)

virtual	~Joint (void)

virtual Elem::Type	GetElemType (void) const

std::ostream &	Output (std::ostream &out, const char *sJointName, unsigned int uLabel, const Vec3 &FLocal, const Vec3 &MLocal, const Vec3 &FGlobal, const Vec3 &MGlobal) const

virtual void	SetInitialValue (VectorHandler &)

virtual void	Update (const VectorHandler &XCurr, InverseDynamics::Order iOrder=InverseDynamics::INVERSE_DYNAMICS)

bool	bIsPrescribedMotion (void) const

bool	bIsTorque (void) const

Public Member Functions inherited from ElemGravityOwner
	ElemGravityOwner (unsigned int uL, flag fOut)

virtual	~ElemGravityOwner (void)

virtual doublereal	dGetM (void) const

Vec3	GetS (void) const

Mat3x3	GetJ (void) const

Vec3	GetB (void) const

Vec3	GetG (void) const

Public Member Functions inherited from GravityOwner
	GravityOwner (void)

virtual	~GravityOwner (void)

void	PutGravity (const Gravity *pG)

virtual bool	bGetGravity (const Vec3 &X, Vec3 &Acc) const

Public Member Functions inherited from ElemWithDofs
	ElemWithDofs (unsigned int uL, const DofOwner *pDO, flag fOut)

virtual	~ElemWithDofs (void)

Public Member Functions inherited from DofOwnerOwner
	DofOwnerOwner (const DofOwner *pDO)

virtual	~DofOwnerOwner ()

virtual const DofOwner *	pGetDofOwner (void) const

virtual integer	iGetFirstIndex (void) const

Public Member Functions inherited from InitialAssemblyElem
	InitialAssemblyElem (unsigned int uL, flag fOut)

virtual	~InitialAssemblyElem (void)

Public Member Functions inherited from SubjectToInitialAssembly
	SubjectToInitialAssembly (void)

virtual	~SubjectToInitialAssembly (void)

Public Member Functions inherited from DriveOwner
	DriveOwner (const DriveCaller *pDC=0)

	DriveOwner (const DriveOwner &drive)

virtual	~DriveOwner (void)

void	Set (const DriveCaller *pDC)

DriveCaller *	pGetDriveCaller (void) const

doublereal	dGet (const doublereal &dVar) const

doublereal	dGet (void) const

bool	bIsDifferentiable (void) const

doublereal	dGetP (const doublereal &dVar) const

doublereal	dGetP (void) const

Protected Attributes
OrientationDescription	od

Protected Attributes inherited from WithLabel
unsigned int	uLabel

std::string	sName

Protected Attributes inherited from ToBeOutput
flag	fOutput

Protected Attributes inherited from GravityOwner
Gravity *	pGravity

Protected Attributes inherited from DriveOwner
DriveCaller *	pDriveCaller

Private Attributes
const StructNode *	pNode1

const StructNode *	pNode2

Vec3	d1

Mat3x3	R1h

Vec3	d2

Mat3x3	R2h

Vec3	F

Vec3	M

int	NTheta

doublereal	dTheta

doublereal	dThetaWrapped

BasicShapeCoefficient *const	Sh_c

BasicFriction *const	fc

const doublereal	preF

const doublereal	r

doublereal	M3

Static Private Attributes
static const unsigned int	NumSelfDof

static const unsigned int	NumDof

Additional Inherited Members
Public Types inherited from Elem
enum	Type { UNKNOWN = -1, AIRPROPERTIES = 0, INDUCEDVELOCITY, AUTOMATICSTRUCTURAL, GRAVITY, BODY, JOINT, JOINT_REGULARIZATION, BEAM, PLATE, FORCE, INERTIA, ELECTRICBULK, ELECTRIC, THERMAL, HYDRAULIC, BULK, LOADABLE, DRIVEN, EXTERNAL, AEROMODAL, AERODYNAMIC, GENEL, SOCKETSTREAM_OUTPUT, RTAI_OUTPUT = SOCKETSTREAM_OUTPUT, LASTELEMTYPE }

Public Types inherited from SimulationEntity
typedef std::vector< Hint * >	Hints

Public Types inherited from ToBeOutput
enum	{ OUTPUT = 0x1U, OUTPUT_MASK = 0xFU, OUTPUT_PRIVATE = 0x10U, OUTPUT_PRIVATE_MASK = ~OUTPUT_MASK }

Public Types inherited from Joint
enum	Type { UNKNOWN = -1, DISTANCE = 0, DISTANCEWITHOFFSET, CLAMP, SPHERICALHINGE, PIN, UNIVERSALHINGE, UNIVERSALROTATION, UNIVERSALPIN, PLANEHINGE, PLANEROTATION, PLANEPIN, AXIALROTATION, PLANEDISP, PLANEDISPPIN, INPLANE, INPLANECONTACT, J_INLINE, ROD, RODBEZIER, DEFORMABLEHINGE, DEFORMABLEDISPJOINT, DEFORMABLEJOINT, DEFORMABLEAXIALJOINT, VISCOUSBODY, LINEARVELOCITY, ANGULARVELOCITY, LINEARACCELERATION, ANGULARACCELERATION, PRISMATIC, DRIVEHINGE, DRIVEDISP, DRIVEDISPPIN, IMPOSEDORIENTATION, IMPOSEDDISP, IMPOSEDDISPPIN, IMPOSEDKINEMATICS, BEAMSLIDER, BRAKE, GIMBAL, POINT_SURFACE_CONTACT, TOTALJOINT, TOTALPINJOINT, TOTALEQUATION, TOTALREACTION, MODAL, SCREWJOINT, LASTJOINTTYPE }

Protected Member Functions inherited from Joint
virtual void	OutputPrepare_int (const std::string &type, OutputHandler &OH, std::string &name)

Protected Member Functions inherited from ElemGravityOwner
virtual Vec3	GetS_int (void) const

virtual Mat3x3	GetJ_int (void) const

virtual Vec3	GetB_int (void) const

virtual Vec3	GetG_int (void) const

Detailed Description

Definition at line 344 of file planej.h.

Constructor & Destructor Documentation

AxialRotationJoint::AxialRotationJoint	(	unsigned int	uL,
		const DofOwner *	pDO,
		const StructNode *	pN1,
		const StructNode *	pN2,
		const Vec3 &	dTmp1,
		const Vec3 &	dTmp2,
		const Mat3x3 &	R1hTmp,
		const Mat3x3 &	R2hTmp,
		const DriveCaller *	pDC,
		const OrientationDescription &	od,
		flag	fOut,
		const doublereal	rr = `0.`,
		const doublereal	pref = `0.`,
		BasicShapeCoefficient *const	sh = `0`,
		BasicFriction *const	f = `0`
	)

Definition at line 2521 of file planej.cc.

References NO_OP.

 : Elem(uL, fOut), 
 Joint(uL, pDO, fOut), 
 DriveOwner(pDC), 
 pNode1(pN1), pNode2(pN2), 
 d1(dTmp1), R1h(R1hTmp), d2(dTmp2), R2h(R2hTmp), F(Zero3), M(Zero3),
 NTheta(0), dTheta(0.), dThetaWrapped(0.),
 #ifdef USE_NETCDF
 Var_Phi(0),
 Var_Omega(0),
 //Var_MFR(0),
 //Var_MU(0),
 #endif // USE_NETCDF
 Sh_c(sh), fc(f), preF(pref), r(rr),
 od(od)
 {
         NO_OP;
 }

AxialRotationJoint::~AxialRotationJoint ( void )

Definition at line 2554 of file planej.cc.

References fc, and Sh_c.

 {
         if (Sh_c) {
                 delete Sh_c;
         }
 
         if (fc) {
                 delete fc;
         }
 }

Member Function Documentation

void AxialRotationJoint::AfterConvergence	(	const VectorHandler &	X,
		const VectorHandler &	XP
	)

virtual

Reimplemented from SimulationEntity.

Definition at line 2952 of file planej.cc.

References BasicFriction::AfterConvergence(), grad::Dot(), dTheta, dThetaWrapped, F, fc, StructNode::GetRCurr(), StructNode::GetWCurr(), DofOwnerOwner::iGetFirstIndex(), M_PI, Vec3::Norm(), NTheta, NumSelfDof, pNode1, pNode2, preF, r, R1h, R2h, and RotManip::VecRot().

 {
         Vec3 v(RotManip::VecRot((pNode1->GetRCurr()*R1h).MulTM(pNode2->GetRCurr()*R2h)));
         doublereal dThetaTmp(v(3));
 
         // unwrap
         if (dThetaTmp - dThetaWrapped < -M_PI) {
                 NTheta++;
         }
 
         if (dThetaTmp - dThetaWrapped > M_PI) {
                 NTheta--;
         }
 
         // save new wrapped angle
         dThetaWrapped = dThetaTmp;
 
         // compute new unwrapped angle
         dTheta = 2*M_PI*NTheta + dThetaWrapped;
 
         if (fc) {
                 const Mat3x3& R1(pNode1->GetRCurr());
                 Mat3x3 R1hTmp(R1*R1h);
                 Vec3 e3a(R1hTmp.GetVec(3));
                 const Vec3& Omega1(pNode1->GetWCurr());
                 const Vec3& Omega2(pNode2->GetWCurr());
                 //relative velocity
                 doublereal v = (Omega1-Omega2).Dot(e3a)*r;
                 //reaction norm
                 doublereal modF = std::max(F.Norm(), preF);;
                 fc->AfterConvergence(modF,v,X,XP,iGetFirstIndex()+NumSelfDof);
         }
 }

Here is the call graph for this function:

VariableSubMatrixHandler & AxialRotationJoint::AssJac	(	VariableSubMatrixHandler &	WorkMat,
		doublereal	dCoef,
		const VectorHandler &	XCurr,
		const VectorHandler &	XPrimeCurr
	)

virtual

Implements Elem.

Definition at line 3007 of file planej.cc.

References ExpandableRowVector::Add(), FullSubMatrixHandler::Add(), BasicFriction::AssJac(), Vec3::Cross(), d1, d2, DEBUGCOUT, Vec3::dGet(), grad::Dot(), BasicShapeCoefficient::dSh_c(), F, BasicFriction::fc(), fc, StructNode::GetRRef(), Mat3x3::GetVec(), StructNode::GetWCurr(), StructNode::GetWRef(), DofOwnerOwner::iGetFirstIndex(), StructDispNode::iGetFirstMomentumIndex(), StructDispNode::iGetFirstPositionIndex(), iGetNumDof(), ExpandableRowVector::Link(), M, MatCross, MatCrossCross, Vec3::Norm(), NumSelfDof, pNode1, pNode2, preF, FullSubMatrixHandler::PutCoef(), FullSubMatrixHandler::PutColIndex(), FullSubMatrixHandler::PutRowIndex(), r, R1h, R2h, ExpandableRowVector::ReDim(), FullSubMatrixHandler::ResizeReset(), ExpandableRowVector::Set(), VariableSubMatrixHandler::SetFull(), BasicShapeCoefficient::Sh_c(), Sh_c, ExpandableRowVector::Sub(), FullSubMatrixHandler::Sub(), and WorkSpaceDim().

 {
    DEBUGCOUT("Entering AxialRotationJoint::AssJac()" << std::endl);
    
    /* Setta la sottomatrice come piena (e' un po' dispersivo, ma lo jacobiano 
     * e' complicato */                                  
    FullSubMatrixHandler& WM = WorkMat.SetFull();
    
    /* Ridimensiona la sottomatrice in base alle esigenze */
    integer iNumRows = 0;
    integer iNumCols = 0;
    WorkSpaceDim(&iNumRows, &iNumCols);
    WM.ResizeReset(iNumRows, iNumCols);
    
    /* Recupera gli indici delle varie incognite */
    integer iNode1FirstPosIndex = pNode1->iGetFirstPositionIndex();
    integer iNode1FirstMomIndex = pNode1->iGetFirstMomentumIndex();
    integer iNode2FirstPosIndex = pNode2->iGetFirstPositionIndex();
    integer iNode2FirstMomIndex = pNode2->iGetFirstMomentumIndex();
    integer iFirstReactionIndex = iGetFirstIndex();
 
    /* Recupera i dati che servono */
    const Mat3x3& R1(pNode1->GetRRef());
    const Mat3x3& R2(pNode2->GetRRef());   
    const Vec3& Omega2(pNode2->GetWRef()); /* Serve dopo */
    Vec3 d1Tmp(R1*d1);
    Vec3 d2Tmp(R2*d2);
    Mat3x3 R1hTmp(R1*R1h);
    Mat3x3 R2hTmp(R2*R2h);
    /* Suppongo che le reazioni F, M siano gia' state aggiornate da AssRes;
     * ricordo che la forza F e' nel sistema globale, mentre la coppia M
     * e' nel sistema locale ed il terzo termine, M(3), e' nullo in quanto
     * diretto come l'asse attorno al quale la rotazione e' consentita */
    
       
    /* 
     * La cerniera piana ha le prime 3 equazioni uguali alla cerniera sferica;
     * inoltre ha due equazioni che affermano la coincidenza dell'asse 3 del
     * riferimento solidale con la cerniera visto dai due nodi.
     * 
     *      (R1 * R1h * e1)^T * (R2 * R2h * e3) = 0
     *      (R1 * R1h * e2)^T * (R2 * R2h * e3) = 0
     * 
     * A queste equazioni corrisponde una reazione di coppia agente attorno 
     * agli assi 1 e 2 del riferimento della cerniera. La coppia attorno 
     * all'asse 3 e' nulla per definizione. Quindi la coppia nel sistema 
     * globale e':
     *      -R1 * R1h * M       per il nodo 1,
     *       R2 * R2h * M       per il nodo 2
     * 
     * 
     *       xa   ga                   xb   gb                     F     M 
     * Qa |  0    0                     0    0                     I     0  | | xa |   | -F           |
     * Ga |  0    c*(F/\da/\-(Sa*M)/\)  0    0                     da/\  Sa | | ga |   | -da/\F-Sa*M |
     * Qb |  0    0                     0    0                    -I     0  | | xb | = |  F           |
     * Gb |  0    0                     0   -c*(F/\db/\-(Sb*M)/\) -db/\ -Sb | | gb |   |  db/\F+Sb*M |
     * F  | -c*I  c*da/\                c*I -c*db/\                0     0  | | F  |   |  xa+da-xb-db |
     * M1 |  0    c*Tb1^T*Ta3/\         0    c*Ta3^T*Tb1/\         0     0  | | M  |   |  Sb^T*Ta3    |
     * M2 |  0    c*Tb2^T*Ta3/\         0    c*Ta3^T*Tb2/\         0     0  | 
     * 
     * con da = R1*d01, db = R2*d02, c = dCoef,
     * Sa = R1*R1h*[e1,e2], Sb = R2*R2h*[e1, e2],
     * Ta3 = R1*R1h*e3, Tb1 = R2*R2h*e1, Tb2 = R2*R2h*e2
     */
 
    /* Setta gli indici delle equazioni */
    for (int iCnt = 1; iCnt <= 6; iCnt++) {
       WM.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
       WM.PutColIndex(iCnt, iNode1FirstPosIndex+iCnt);
       WM.PutRowIndex(6+iCnt, iNode2FirstMomIndex+iCnt);
       WM.PutColIndex(6+iCnt, iNode2FirstPosIndex+iCnt);
    }
    
    for (unsigned int iCnt = 1; iCnt <= iGetNumDof(); iCnt++) {
       WM.PutRowIndex(12+iCnt, iFirstReactionIndex+iCnt);
       WM.PutColIndex(12+iCnt, iFirstReactionIndex+iCnt);
    }
    
    /* Contributo della forza alle equazioni di equilibrio dei due nodi */
    for (int iCnt = 1; iCnt <= 3; iCnt++) {      
       WM.PutCoef(iCnt, 12+iCnt, 1.);
       WM.PutCoef(6+iCnt, 12+iCnt, -1.);
    }
    
    WM.Add(4, 13, Mat3x3(MatCross, d1Tmp));
    WM.Sub(10, 13, Mat3x3(MatCross, d2Tmp));   
    
    /* Moltiplica la forza ed il momento per il coefficiente
     * del metodo */
    Vec3 FTmp = F*dCoef;
    Vec3 MTmp = M*dCoef;
    
    Vec3 e3a(R1hTmp.GetVec(3));
    Vec3 e1b(R2hTmp.GetVec(1));
    Vec3 e2b(R2hTmp.GetVec(2));
    MTmp = e2b*MTmp.dGet(1) - e1b*MTmp.dGet(2);
    
    Mat3x3 MWedgee3aWedge(MatCrossCross, MTmp, e3a);
    Mat3x3 e3aWedgeMWedge(MatCrossCross, e3a, MTmp);
    
    WM.Add(4, 4, Mat3x3(MatCrossCross, FTmp, d1Tmp) - MWedgee3aWedge - Mat3x3(MatCross, e3a*M(3)));
    WM.Add(4, 10, e3aWedgeMWedge);
    
    WM.Add(10, 4, MWedgee3aWedge);
    WM.Sub(10, 10, Mat3x3(MatCrossCross, FTmp, d2Tmp) + e3aWedgeMWedge
                    - Mat3x3(MatCross, e3a*M(3)));
  
    /* Contributo del momento alle equazioni di equilibrio dei nodi */
    Vec3 Tmp1(e2b.Cross(e3a));
    Vec3 Tmp2(e3a.Cross(e1b));
    
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       doublereal d = Tmp1(iCnt);
       WM.PutCoef(3+iCnt, 16, d);
       WM.PutCoef(9+iCnt, 16, -d);
       
       d = Tmp2(iCnt);
       WM.PutCoef(3+iCnt, 17, d);
       WM.PutCoef(9+iCnt, 17, -d);
       
       d = e3a(iCnt);
       WM.PutCoef(3+iCnt, 18, d);        
       WM.PutCoef(9+iCnt, 18, -d);       
    }         
    
    /* Modifica: divido le equazioni di vincolo per dCoef */
    
    /* Equazioni di vincolo degli spostamenti */
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       WM.PutCoef(12+iCnt, iCnt, -1.);
       WM.PutCoef(12+iCnt, 6+iCnt, 1.);
    }
    
    WM.Add(13, 4, Mat3x3(MatCross, d1Tmp));
    WM.Sub(13, 10, Mat3x3(MatCross, d2Tmp));
    
    /* Equazione di vincolo del momento
     * 
     * Attenzione: bisogna scrivere il vettore trasposto
     *   (Sb[1]^T*(Sa[3]/\))*dCoef
     * Questo pero' e' uguale a:
     *   (-Sa[3]/\*Sb[1])^T*dCoef,
     * che puo' essere ulteriormente semplificato:
     *   (Sa[3].Cross(Sb[1])*(-dCoef))^T;
     */
    
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       doublereal d = Tmp1(iCnt);
       WM.PutCoef(16, 3+iCnt, d);
       WM.PutCoef(16, 9+iCnt, -d);
       
       d = Tmp2.dGet(iCnt);
       WM.PutCoef(17, 3+iCnt, -d);
       WM.PutCoef(17, 9+iCnt, d);
    }
    
    /* Questa equazione non viene divisa per dCoef */
    
    /* Equazione di imposizione della velocita' di rotazione: 
     * (e3a+c(wb/\e3a))^T*(Delta_gPb-Delta_gPa) */
    Vec3 Tmp = e3a + Omega2.Cross(e3a*dCoef);
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       doublereal d = Tmp(iCnt);
       WM.PutCoef(18, 3+iCnt, -d);
       WM.PutCoef(18, 9+iCnt, d);
    }   
 
    if (fc) {
       //retrive
           //friction coef
       doublereal f = fc->fc();
           //shape function
       doublereal shc = Sh_c->Sh_c();
           //omega and omega rif
       const Vec3& Omega1(pNode1->GetWCurr());
       const Vec3& Omega2(pNode2->GetWCurr());
       // const Vec3& Omega1r(pNode1->GetWRef());
       // const Vec3& Omega2r(pNode2->GetWRef());   
       //compute 
           //relative velocity
       doublereal v = (Omega1-Omega2).Dot(e3a)*r;
           //reaction norm
       doublereal modF = std::max(F.Norm(), preF);
           //reaction moment
       //doublereal M3 = shc*modF*f*r;
       
       ExpandableRowVector dfc;
       ExpandableRowVector dF;
       ExpandableRowVector dv;
           //variation of reaction force
       dF.ReDim(3);
       if ((modF == 0.) or (F.Norm() > preF)) {
           dF.Set(0.,1,12+1);
           dF.Set(0.,2,12+2);
           dF.Set(0.,3,12+3);
       } else {
           dF.Set(F.dGet(1)/modF,1,12+1);
           dF.Set(F.dGet(2)/modF,2,12+2);
           dF.Set(F.dGet(3)/modF,3,12+3);
       }
           //variation of relative velocity
       dv.ReDim(6);
       
 /* (approximate: assume constant triads orientations) 
  * relative velocity linearization 
 */
       dv.Set(e3a*r,1, 0+4);
       dv.Set(-e3a*r,4, 6+4);
 
       //assemble friction states
       fc->AssJac(WM,dfc,12+NumSelfDof,iFirstReactionIndex+NumSelfDof,dCoef,modF,v,
                 XCurr,XPrimeCurr,dF,dv);
       ExpandableRowVector dM3;
       ExpandableRowVector dShc;
       //compute 
           //variation of shape function
       Sh_c->dSh_c(dShc,f,modF,v,dfc,dF,dv);
           //variation of moment component
       dM3.ReDim(2);
       dM3.Set(shc * r,1); dM3.Link(1,&dF);
       dM3.Set(modF * r,2); dM3.Link(2,&dShc);
       //assemble first node
           //variation of moment component
       dM3.Add(WM,0+4,e3a.dGet(1));
       dM3.Add(WM,0+5,e3a.dGet(2));
       dM3.Add(WM,0+6,e3a.dGet(3));
       //assemble second node
           //variation of moment component
       dM3.Sub(WM,6+4,e3a.dGet(1));
       dM3.Sub(WM,6+5,e3a.dGet(2));
       dM3.Sub(WM,6+6,e3a.dGet(3));
    }
    
    return WorkMat;
 }

Here is the call graph for this function:

SubVectorHandler & AxialRotationJoint::AssRes	(	SubVectorHandler &	WorkVec,
		doublereal	dCoef,
		const VectorHandler &	XCurr,
		const VectorHandler &	XPrimeCurr
	)

virtual

!!!!!!!!!!!!!

Implements Elem.

Definition at line 3248 of file planej.cc.

References VectorHandler::Add(), ASSERT, BasicFriction::AssRes(), Vec3::Cross(), d1, d2, DEBUGCOUT, Vec3::dGet(), DriveCaller::dGet(), Vec3::Dot(), grad::Dot(), F, BasicFriction::fc(), fc, StructNode::GetRCurr(), Mat3x3::GetVec(), StructNode::GetWCurr(), StructDispNode::GetXCurr(), DofOwnerOwner::iGetFirstIndex(), StructDispNode::iGetFirstMomentumIndex(), iGetNumDof(), M, M3, MBDYN_EXCEPT_ARGS, Vec3::Norm(), NumSelfDof, DriveOwner::pGetDriveCaller(), pNode1, pNode2, preF, VectorHandler::PutCoef(), SubVectorHandler::PutRowIndex(), r, R1h, R2h, VectorHandler::ResizeReset(), BasicShapeCoefficient::Sh_c(), Sh_c, VectorHandler::Sub(), and WorkSpaceDim().

 {
    DEBUGCOUT("Entering AxialRotationJoint::AssRes()" << std::endl);
 
    /* Dimensiona e resetta la matrice di lavoro */
    integer iNumRows = 0;
    integer iNumCols = 0;
    WorkSpaceDim(&iNumRows, &iNumCols);
    WorkVec.ResizeReset(iNumRows);
  
    /* Indici */
    integer iNode1FirstMomIndex = pNode1->iGetFirstMomentumIndex();
    integer iNode2FirstMomIndex = pNode2->iGetFirstMomentumIndex();
    integer iFirstReactionIndex = iGetFirstIndex();
    
    /* Indici dei nodi */
    for (int iCnt = 1; iCnt <= 6; iCnt++) {      
       WorkVec.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
       WorkVec.PutRowIndex(6+iCnt, iNode2FirstMomIndex+iCnt);
    }
    /* Indici del vincolo */
    for (unsigned int iCnt = 1; iCnt <= iGetNumDof(); iCnt++) {
       WorkVec.PutRowIndex(12+iCnt, iFirstReactionIndex+iCnt);
    }
    
    /* Aggiorna i dati propri */
    F = Vec3(XCurr, iFirstReactionIndex+1);
    M = Vec3(XCurr, iFirstReactionIndex+4);
    
    /* Recupera i dati */
    const Vec3& x1(pNode1->GetXCurr());
    const Vec3& x2(pNode2->GetXCurr());
    const Mat3x3& R1(pNode1->GetRCurr());
    const Mat3x3& R2(pNode2->GetRCurr());
    
    /* Costruisce i dati propri nella configurazione corrente */
    Vec3 dTmp1(R1*d1);
    Vec3 dTmp2(R2*d2);
    Mat3x3 R1hTmp(R1*R1h);
    Mat3x3 R2hTmp(R2*R2h);
    
    Vec3 e3a(R1hTmp.GetVec(3));
    Vec3 e1b(R2hTmp.GetVec(1));
    Vec3 e2b(R2hTmp.GetVec(2));
    
    Vec3 MTmp(e2b.Cross(e3a)*M.dGet(1)+e3a.Cross(e1b)*M.dGet(2));
    
    /* Equazioni di equilibrio, nodo 1 */
    WorkVec.Sub(1, F);
    WorkVec.Add(4, F.Cross(dTmp1)-MTmp-e3a*M.dGet(3)); /* Sfrutto  F/\d = -d/\F */
    
    /* Equazioni di equilibrio, nodo 2 */
    WorkVec.Add(7, F);
    WorkVec.Add(10, dTmp2.Cross(F)+MTmp+e3a*M.dGet(3));
 
    /* Modifica: divido le equazioni di vincolo per dCoef */
    ASSERT(dCoef != 0.);
       
    /* Equazione di vincolo di posizione */
    WorkVec.Add(13, (x1+dTmp1-x2-dTmp2)/dCoef);
       
    /* Equazioni di vincolo di rotazione */
    Vec3 Tmp = Vec3(R1hTmp.GetVec(3));
    WorkVec.PutCoef(16, Tmp.Dot(R2hTmp.GetVec(2))/dCoef);
    WorkVec.PutCoef(17, Tmp.Dot(R2hTmp.GetVec(1))/dCoef);
    
    /* Questa equazione non viene divisa per dCoef */
    
    /* Equazione di vincolo di velocita' di rotazione */
    const Vec3& Omega1(pNode1->GetWCurr());
    const Vec3& Omega2(pNode2->GetWCurr());
    doublereal dOmega0 = pGetDriveCaller()->dGet();
    WorkVec.PutCoef(18, dOmega0-e3a.Dot(Omega2-Omega1));
 
    if (fc) {
       bool ChangeJac(false);
       doublereal v = (Omega1-Omega2).Dot(e3a)*r;
       doublereal modF = std::max(F.Norm(), preF);
       try {
           fc->AssRes(WorkVec,12+NumSelfDof,iFirstReactionIndex+NumSelfDof,modF,v,XCurr,XPrimeCurr);
       }
       catch (Elem::ChangedEquationStructure) {
           ChangeJac = true;
       }
       doublereal f = fc->fc();
       doublereal shc = Sh_c->Sh_c(f,modF,v);
       M3 = shc*modF*r;
       WorkVec.Sub(4,e3a*M3);
       WorkVec.Add(10,e3a*M3);
 //      M += e3a*M3;
       if (ChangeJac) {
           throw Elem::ChangedEquationStructure(MBDYN_EXCEPT_ARGS);
       }
    }
    
    return WorkVec;
 }

Here is the call graph for this function:

std::ostream & AxialRotationJoint::DescribeDof	(	std::ostream &	out,
		const char *	prefix = `""`,
		bool	bInitial = `false`
	)		const

virtual

Reimplemented from Elem.

Definition at line 2567 of file planej.cc.

References SimulationEntity::DescribeDof(), fc, DofOwnerOwner::iGetFirstIndex(), SimulationEntity::iGetNumDof(), and NumSelfDof.

 {
         integer iIndex = iGetFirstIndex();
 
         out
                 << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
                         "reaction forces [Fx,Fy,Fz]" << std::endl
                 << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
                         "reaction couples [mx,my,mz]" << std::endl;
 
         if (bInitial) {
                 iIndex += NumSelfDof;
                 out
                         << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
                                 "reaction force derivatives [FPx,FPy,FPz]" << std::endl
                         << prefix << iIndex + 4 << "->" << iIndex + 5 << ": "
                                 "reaction couple derivatives [mPx,mPy]" << std::endl;
         }
         
         iIndex += NumSelfDof;
         if (fc) {
                 integer iFCDofs = fc->iGetNumDof();
                 if (iFCDofs > 0) {
                         out << prefix << iIndex + 1;
                         if (iFCDofs > 1) {
                                 out << "->" << iIndex + iFCDofs;
                         }
                         out << ": friction dof(s)" << std::endl
                                 << "        ", fc->DescribeDof(out, prefix, bInitial);
                 }
         }
 
         return out;
 }

Here is the call graph for this function:

void AxialRotationJoint::DescribeDof	(	std::vector< std::string > &	desc,
		bool	bInitial = `false`,
		int	i = `-1`
	)		const

virtual

Reimplemented from Elem.

Definition at line 2603 of file planej.cc.

References SimulationEntity::DescribeDof(), fc, WithLabel::GetLabel(), MBDYN_EXCEPT_ARGS, NumSelfDof, and xyz.

 {
         std::ostringstream os;
         os << "AxialRotationJoint(" << GetLabel() << ")";
 
         unsigned short nself = NumSelfDof;
         if (bInitial) {
                 nself += NumSelfDof - 1;
         }
         if (fc && (i == -1 || i >= nself)) {
                 fc->DescribeDof(desc, bInitial, i - nself);
                 if (i != -1) {
                         desc[0] = os.str() + ": " + desc[0];
                         return;
                 }
         }
 
         if (i == -1) {
                 // move fc desc to the end
                 unsigned short nfc = 0;
                 if (fc) {
                         nfc = desc.size();
                 }
                 desc.resize(nfc + nself);
                 for (unsigned i = nfc; i-- > 0; ) {
                         desc[nself + i] = os.str() + ": " + desc[nfc];
                 }
 
                 std::string name = os.str();
 
                 for (unsigned i = 0; i < 3; i++) {
                         os.str(name);
                         os.seekp(0, std::ios_base::end);
                         os << ": reaction force f" << xyz[i];
                         desc[i] = os.str();
                 }
 
                 for (unsigned i = 0; i < 3; i++) {
                         os.str(name);
                         os.seekp(0, std::ios_base::end);
                         os << ": reaction couple m" << xyz[i];
                         desc[3 + i] = os.str();
                 }
 
                 if (bInitial) {
                         for (unsigned i = 0; i < 3; i++) {
                                 os.str(name);
                                 os.seekp(0, std::ios_base::end);
                                 os << ": reaction force derivative fP" << xyz[i];
                                 desc[6 + i] = os.str();
                         }
         
                         for (unsigned i = 0; i < 2; i++) {
                                 os.str(name);
                                 os.seekp(0, std::ios_base::end);
                                 os << ": reaction couple derivative mP" << xyz[i];
                                 desc[9 + i] = os.str();
                         }
                 }
 
         } else {
                 if (i < -1) {
                         // error
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
 
                 if (i >= nself) {
                         // error
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
 
                 desc.resize(1);
 
                 switch (i) {
                 case 0:
                 case 1:
                 case 2:
                         os << ": reaction force f" << xyz[i];
                         break;
 
                 case 3:
                 case 4:
                 case 5:
                         os << ": reaction couple m" << xyz[i - 3];
                         break;
 
                 case 6:
                 case 7:
                 case 8:
                         os << ": reaction force derivative fP" << xyz[i - 6];
                         break;
 
                 case 9:
                 case 10:
                         os << ": reaction couple derivative mP" << xyz[i - 9];
                         break;
                 }
                 desc[0] = os.str();
         }
 }

Here is the call graph for this function:

std::ostream & AxialRotationJoint::DescribeEq	(	std::ostream &	out,
		const char *	prefix = `""`,
		bool	bInitial = `false`
	)		const

virtual

Reimplemented from Elem.

Definition at line 2705 of file planej.cc.

References SimulationEntity::DescribeEq(), fc, DofOwnerOwner::iGetFirstIndex(), SimulationEntity::iGetNumDof(), and NumSelfDof.

 {
         integer iIndex = iGetFirstIndex();
 
         out
                 << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
                         "position constraints [Px1=Px2,Py1=Py2,Pz1=Pz2]" << std::endl
                 << prefix << iIndex + 4 << "->" << iIndex + 5 << ": "
                         "orientation constraints [gx1=gx2,gy1=gy2]" << std::endl
                 << prefix << iIndex + 6 << ": "
                         "angular velocity constraint wz2-wz1=Omega]" << std::endl;
 
         if (bInitial) {
                 iIndex += NumSelfDof;
                 out
                         << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
                                 "velocity constraints [vx1=vx2,vy1=vy2,vz1=vz2]" << std::endl
                         << prefix << iIndex + 4 << "->" << iIndex + 5 << ": "
                                 "reaction couple derivatives [wx1=wx2,wy1=wy2]" << std::endl;
         }
         
         iIndex += NumSelfDof;
         if (fc) {
                 integer iFCDofs = fc->iGetNumDof();
                 if (iFCDofs > 0) {
                         out << prefix << iIndex + 1;
                         if (iFCDofs > 1) {
                                 out << "->" << iIndex + iFCDofs;
                         }
                         out << ": friction equation(s)" << std::endl
                                 << "        ", fc->DescribeEq(out, prefix, bInitial);
                 }
         }
 
         return out;
 }

Here is the call graph for this function:

void AxialRotationJoint::DescribeEq	(	std::vector< std::string > &	desc,
		bool	bInitial = `false`,
		int	i = `-1`
	)		const

virtual

Reimplemented from Elem.

Definition at line 2743 of file planej.cc.

References SimulationEntity::DescribeEq(), fc, WithLabel::GetLabel(), MBDYN_EXCEPT_ARGS, NumSelfDof, and xyz.

 {
         std::ostringstream os;
         os << "AxialRotationJoint(" << GetLabel() << ")";
 
         unsigned short nself = NumSelfDof;
         if (bInitial) {
                 nself += NumSelfDof - 1;
         }
         if (fc && (i == -1 || i >= nself)) {
                 fc->DescribeEq(desc, bInitial, i - nself);
                 if (i != -1) {
                         desc[0] = os.str() + ": " + desc[0];
                         return;
                 }
         }
 
         if (i == -1) {
                 // move fc desc to the end
                 unsigned short nfc = 0;
                 if (fc) {
                         nfc = desc.size();
                 }
                 desc.resize(nfc + nself);
                 for (unsigned i = nfc; i-- > 0; ) {
                         desc[nself + i] = os.str() + ": " + desc[nfc];
                 }
 
                 std::string name = os.str();
 
                 for (unsigned i = 0; i < 3; i++) {
                         os.str(name);
                         os.seekp(0, std::ios_base::end);
                         os << ": position constraint P" << xyz[i];
                         desc[i] = os.str();
                 }
 
                 for (unsigned i = 0; i < 2; i++) {
                         os.str(name);
                         os.seekp(0, std::ios_base::end);
                         os << ": orientation constraint g" << xyz[i];
                         desc[3 + i] = os.str();
                 }
 
                 os.str(name);
                 os.seekp(0, std::ios_base::end);
                 os << ": angular velocity constraint wz";
                 desc[5] = os.str();
 
                 if (bInitial) {
                         for (unsigned i = 0; i < 3; i++) {
                                 os.str(name);
                                 os.seekp(0, std::ios_base::end);
                                 os << ": position constraint derivative v" << xyz[i];
                                 desc[6 + i] = os.str();
                         }
         
                         for (unsigned i = 0; i < 2; i++) {
                                 os.str(name);
                                 os.seekp(0, std::ios_base::end);
                                 os << ": orientation constraint derivative w" << xyz[i];
                                 desc[9 + i] = os.str();
                         }
                 }
 
         } else {
                 if (i < -1) {
                         // error
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
 
                 if (i >= nself) {
                         // error
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
 
                 desc.resize(1);
 
                 switch (i) {
                 case 0:
                 case 1:
                 case 2:
                         os << ": position constraint P" << xyz[i];
                         break;
 
                 case 3:
                 case 4:
                         os << ": orientation constraint g" << xyz[i - 3];
                         break;
 
                 case 5:
                         os << ": angular velocity constraint wz";
                         break;
 
                 case 6:
                 case 7:
                 case 8:
                         os << ": position constraint derivative v" << xyz[i - 6];
                         break;
 
                 case 9:
                 case 10:
                         os << ": orientation constraint derivative w" << xyz[i - 9];
                         break;
                 }
                 desc[0] = os.str();
         }
 }

Here is the call graph for this function:

doublereal AxialRotationJoint::dGetPrivData ( unsigned int i ) const

virtual

Reimplemented from SimulationEntity.

Definition at line 3930 of file planej.cc.

References ASSERT, DriveOwner::dGet(), dThetaWrapped, F, WithLabel::GetLabel(), StructNode::GetRCurr(), iGetNumPrivData(), M, M_PI, MBDYN_EXCEPT_ARGS, NTheta, pNode1, pNode2, R1h, R2h, and RotManip::VecRot().

 {
         ASSERT(i >= 1 && i <= iGetNumPrivData());
    
         switch (i) {
         case 1: {
                 Vec3 v(RotManip::VecRot((pNode1->GetRCurr()*R1h).MulTM(pNode2->GetRCurr()*R2h)));
                 doublereal dThetaTmp(v(3));
 
                 int n = 0;
 
                 if (dThetaTmp - dThetaWrapped < -M_PI) {
                         n++;
                 }
 
                 if (dThetaTmp - dThetaWrapped > M_PI) {
                         n--;
                 }
 
                 return 2*M_PI*(NTheta + n) + dThetaTmp;
         }
       
         case 2: 
                 return dGet();
       
         case 3:
         case 4:
         case 5:
                 return F(i - 2);
       
         case 6:
         case 7:
         case 8:
                 return M(i - 5);
         }
 
         silent_cerr("AxialRotationJoint(" << GetLabel() << "): "
                 "illegal private data " << i << std::endl);
         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 }

Here is the call graph for this function:

virtual void AxialRotationJoint::GetConnectedNodes ( std::vector< const Node * > & connectedNodes ) const

inlinevirtual

Reimplemented from Elem.

Definition at line 506 of file planej.h.

References pNode1, and pNode2.

                                                                                  {
      connectedNodes.resize(2);
      connectedNodes[0] = pNode1;
      connectedNodes[1] = pNode2;
    };

DofOrder::Order AxialRotationJoint::GetDofType ( unsigned int i ) const

inlinevirtual

Reimplemented from Elem.

Definition at line 435 of file planej.h.

References DofOrder::ALGEBRAIC, ASSERT, fc, SimulationEntity::GetDofType(), iGetNumDof(), and NumSelfDof.

                                                   {
       ASSERT(i >= 0 && i < iGetNumDof());
       if (i<NumSelfDof) {
           return DofOrder::ALGEBRAIC; 
       } else {
           return fc->GetDofType(i-NumSelfDof);
       }
    };

Here is the call graph for this function:

DofOrder::Order AxialRotationJoint::GetEqType ( unsigned int i ) const

virtual

Reimplemented from SimulationEntity.

Definition at line 3351 of file planej.cc.

References DofOrder::ALGEBRAIC, ASSERTMSGBREAK, DofOrder::DIFFERENTIAL, fc, SimulationEntity::GetEqType(), iGetNumDof(), and NumSelfDof.

 {
         ASSERTMSGBREAK(i < iGetNumDof(),
                 "INDEX ERROR in AxialRotationJoint::GetEqType");
         if (i == 5) {
                 return DofOrder::DIFFERENTIAL;
         } else if (i < NumSelfDof) {
                 return DofOrder::ALGEBRAIC;
         } else {
                 return fc->GetEqType(i-NumSelfDof);
         }
 }

Here is the call graph for this function:

virtual Joint::Type AxialRotationJoint::GetJointType ( void ) const

inlinevirtual

Implements Joint.

Definition at line 406 of file planej.h.

References Joint::AXIALROTATION.

                                               { 
       return Joint::AXIALROTATION; 
    };

virtual unsigned int AxialRotationJoint::iGetInitialNumDof ( void ) const

inlinevirtual

Implements SubjectToInitialAssembly.

Definition at line 481 of file planej.h.

                                                       { 
       return 11;
    };

virtual unsigned int AxialRotationJoint::iGetNumDof ( void ) const

inlinevirtual

Reimplemented from Elem.

Definition at line 413 of file planej.h.

References fc, SimulationEntity::iGetNumDof(), and NumSelfDof.

Referenced by AssJac(), AssRes(), GetDofType(), and GetEqType().

                                                { 
        unsigned int i = NumSelfDof;
        if (fc) {
            i+=fc->iGetNumDof();
        } 
        return i;
    };

Here is the call graph for this function:

unsigned int AxialRotationJoint::iGetNumPrivData ( void ) const

virtual

Reimplemented from SimulationEntity.

Definition at line 3887 of file planej.cc.

Referenced by dGetPrivData().

 {
         return 8;
 }

unsigned int AxialRotationJoint::iGetPrivDataIdx ( const char * s ) const

virtual

Reimplemented from SimulationEntity.

Definition at line 3893 of file planej.cc.

References ASSERT.

 {
         ASSERT(s != NULL);
 
         unsigned int idx = 0;
 
         switch (s[0]) {
         case 'w':
                 idx++;
         case 'r':
                 idx++;
                 if (s[1] == 'z') {
                         return idx;
                 }
                 break;
 
         case 'M':
                 idx += 3;
         case 'F':
                 idx += 2;
 
                 switch (s[1]) {
                 case 'x':
                         return idx + 1;
 
                 case 'y':
                         return idx + 2;
 
                 case 'z':
                         return idx + 3;
                 }
         }
 
         return 0;
 }

VariableSubMatrixHandler & AxialRotationJoint::InitialAssJac	(	VariableSubMatrixHandler &	WorkMat,
		const VectorHandler &	XCurr
	)

virtual

Implements SubjectToInitialAssembly.

Definition at line 3494 of file planej.cc.

References FullSubMatrixHandler::Add(), Vec3::Cross(), d1, d2, Vec3::dGet(), F, WithLabel::GetLabel(), StructNode::GetRRef(), Mat3x3::GetVec(), StructNode::GetWRef(), DofOwnerOwner::iGetFirstIndex(), StructDispNode::iGetFirstPositionIndex(), InitialWorkSpaceDim(), M, MatCross, MatCrossCross, MBDYN_EXCEPT_ARGS, pNode1, pNode2, FullSubMatrixHandler::PutCoef(), FullSubMatrixHandler::PutColIndex(), FullSubMatrixHandler::PutRowIndex(), R1h, R2h, FullSubMatrixHandler::ResizeReset(), VariableSubMatrixHandler::SetFull(), and FullSubMatrixHandler::Sub().

 {
    /* Per ora usa la matrice piena; eventualmente si puo' 
     * passare a quella sparsa quando si ottimizza */
    FullSubMatrixHandler& WM = WorkMat.SetFull();
    
    /* Dimensiona e resetta la matrice di lavoro */
    integer iNumRows = 0;
    integer iNumCols = 0;
    InitialWorkSpaceDim(&iNumRows, &iNumCols);
    WM.ResizeReset(iNumRows, iNumCols);
 
    /* Equazioni: vedi joints.dvi */
    
    /*       equazioni ed incognite
     * F1                                     Delta_x1         0+1 =  1
     * M1                                     Delta_g1         3+1 =  4
     * FP1                                    Delta_xP1        6+1 =  7
     * MP1                                    Delta_w1         9+1 = 10
     * F2                                     Delta_x2        12+1 = 13
     * M2                                     Delta_g2        15+1 = 16
     * FP2                                    Delta_xP2       18+1 = 19
     * MP2                                    Delta_w2        21+1 = 22
     * vincolo spostamento                    Delta_F         24+1 = 25
     * vincolo rotazione                      Delta_M         27+1 = 28
     * derivata vincolo spostamento           Delta_FP        29+1 = 30
     * derivata vincolo rotazione             Delta_MP        32+1 = 33
     * vincolo di velocita' di rotazione
     */
         
     
    /* Indici */
    integer iNode1FirstPosIndex = pNode1->iGetFirstPositionIndex();
    integer iNode1FirstVelIndex = iNode1FirstPosIndex+6;
    integer iNode2FirstPosIndex = pNode2->iGetFirstPositionIndex();
    integer iNode2FirstVelIndex = iNode2FirstPosIndex+6;
    integer iFirstReactionIndex = iGetFirstIndex();
    integer iReactionPrimeIndex = iFirstReactionIndex+5;
    
    /* Nota: le reazioni vincolari sono: 
     * Forza,       3 incognite, riferimento globale, 
     * Momento,     2 incognite, riferimento locale
     */
 
    /* Setta gli indici dei nodi */
    for (int iCnt = 1; iCnt <= 6; iCnt++) {
       WM.PutRowIndex(iCnt, iNode1FirstPosIndex+iCnt);
       WM.PutColIndex(iCnt, iNode1FirstPosIndex+iCnt);
       WM.PutRowIndex(6+iCnt, iNode1FirstVelIndex+iCnt);
       WM.PutColIndex(6+iCnt, iNode1FirstVelIndex+iCnt);
       WM.PutRowIndex(12+iCnt, iNode2FirstPosIndex+iCnt);
       WM.PutColIndex(12+iCnt, iNode2FirstPosIndex+iCnt);
       WM.PutRowIndex(18+iCnt, iNode2FirstVelIndex+iCnt);
       WM.PutColIndex(18+iCnt, iNode2FirstVelIndex+iCnt);
    }
    
    /* Setta gli indici delle reazioni */
    for (int iCnt = 1; iCnt <= 5; iCnt++) {
       WM.PutRowIndex(24+iCnt, iFirstReactionIndex+iCnt);
       WM.PutColIndex(24+iCnt, iFirstReactionIndex+iCnt);        
    }   
    
    for (int iCnt = 1; iCnt <= 6; iCnt++) {
       WM.PutRowIndex(29+iCnt, iReactionPrimeIndex+iCnt);
       WM.PutColIndex(29+iCnt, iReactionPrimeIndex+iCnt);        
    }   
    
    /* Recupera i dati */
    const Mat3x3& R1(pNode1->GetRRef());
    const Mat3x3& R2(pNode2->GetRRef());
    const Vec3& Omega1(pNode1->GetWRef());
    const Vec3& Omega2(pNode2->GetWRef());   
    /* F ed M sono gia' state aggiornate da InitialAssRes */
    Vec3 FPrime(XCurr, iReactionPrimeIndex+1);
    Vec3 MPrime(XCurr, iReactionPrimeIndex+4);
 
    /* Matrici identita' */   
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       /* Contributo di forza all'equazione della forza, nodo 1 */
       WM.PutCoef(iCnt, 24+iCnt, 1.);
       
       /* Contrib. di der. di forza all'eq. della der. della forza, nodo 1 */
       WM.PutCoef(6+iCnt, 29+iCnt, 1.);
       
       /* Contributo di forza all'equazione della forza, nodo 2 */
       WM.PutCoef(12+iCnt, 24+iCnt, -1.);
       
       /* Contrib. di der. di forza all'eq. della der. della forza, nodo 2 */
       WM.PutCoef(18+iCnt, 29+iCnt, -1.);
       
       /* Equazione di vincolo, nodo 1 */
       WM.PutCoef(24+iCnt, iCnt, -1.);
       
       /* Derivata dell'equazione di vincolo, nodo 1 */
       WM.PutCoef(29+iCnt, 6+iCnt, -1.);
       
       /* Equazione di vincolo, nodo 2 */
       WM.PutCoef(24+iCnt, 12+iCnt, 1.);
       
       /* Derivata dell'equazione di vincolo, nodo 2 */
       WM.PutCoef(29+iCnt, 18+iCnt, 1.);
    }
       
    /* Distanze e matrici di rotazione dai nodi alla cerniera 
     * nel sistema globale */
    Vec3 d1Tmp(R1*d1);
    Vec3 d2Tmp(R2*d2);
    Mat3x3 R1hTmp(R1*R1h);
    Mat3x3 R2hTmp(R2*R2h);
    
    Vec3 e3a(R1hTmp.GetVec(3));
    Vec3 e1b(R2hTmp.GetVec(1));
    Vec3 e2b(R2hTmp.GetVec(2));
    
    /* Ruota il momento e la sua derivata con le matrici della cerniera 
     * rispetto ai nodi */
    Vec3 MTmp(e2b*M.dGet(1)-e1b*M.dGet(2));
    Vec3 MPrimeTmp(e2b*MPrime.dGet(1)-e1b*MPrime.dGet(2));
 
    /* Matrici F/\d1/\, -F/\d2/\ */
    Mat3x3 FWedged1Wedge(MatCrossCross, F, d1Tmp);
    Mat3x3 FWedged2Wedge(MatCrossCross, F, -d2Tmp);
    
    /* Matrici (omega1/\d1)/\, -(omega2/\d2)/\ */
    Mat3x3 O1Wedged1Wedge(MatCross, Omega1.Cross(d1Tmp));
    Mat3x3 O2Wedged2Wedge(MatCross, d2Tmp.Cross(Omega2));
    
    Mat3x3 MDeltag1((Mat3x3(MatCross, Omega2.Cross(MTmp) + MPrimeTmp)
                            + Mat3x3(MatCrossCross, MTmp, Omega1))*Mat3x3(MatCross, e3a));
    Mat3x3 MDeltag2(Mat3x3(MatCrossCross, Omega1.Cross(e3a), MTmp)
                    + Mat3x3(MatCrossCross, e3a, MPrimeTmp)
                    + e3a.Cross(Mat3x3(MatCrossCross, Omega2, MTmp)));
 
    /* Vettori temporanei */
    Vec3 Tmp1(e2b.Cross(e3a));   
    Vec3 Tmp2(e3a.Cross(e1b));
    
    /* Prodotto vettore tra il versore 3 della cerniera secondo il nodo 1
     * ed il versore 1 della cerniera secondo il nodo 2. A convergenza
     * devono essere ortogonali, quindi il loro prodotto vettore deve essere 
     * unitario */
 
    /* Error handling: il programma si ferma, pero' segnala dov'e' l'errore */
    if (Tmp1.Dot() <= std::numeric_limits<doublereal>::epsilon() || Tmp2.Dot() <= std::numeric_limits<doublereal>::epsilon()) {
       silent_cerr("AxialRotationJoint(" << GetLabel() << "): "
               "first and second node hinge axes are (nearly) orthogonal" 
               << std::endl);
       throw ErrNullNorm(MBDYN_EXCEPT_ARGS);
    }   
    
    Vec3 TmpPrime1(e2b.Cross(Omega1.Cross(e3a))-e3a.Cross(Omega2.Cross(e2b)));
    Vec3 TmpPrime2(e3a.Cross(Omega2.Cross(e1b))-e1b.Cross(Omega1.Cross(e3a)));
    
    /* Equazione di momento, nodo 1 */
    WM.Add(4, 4, FWedged1Wedge - Mat3x3(MatCrossCross, MTmp, e3a));
    WM.Add(4, 16, Mat3x3(MatCrossCross, e3a, MTmp));
    WM.Add(4, 25, Mat3x3(MatCross, d1Tmp));
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       WM.PutCoef(3+iCnt, 28, Tmp1(iCnt));
       WM.PutCoef(3+iCnt, 29, Tmp2(iCnt));
    }
    
    /* Equazione di momento, nodo 2 */
    WM.Add(4, 16, Mat3x3(MatCrossCross, MTmp, e3a));
    WM.Add(16, 16, FWedged2Wedge - Mat3x3(MatCrossCross, e3a, MTmp));
    WM.Sub(16, 25, Mat3x3(MatCross, d2Tmp));
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       WM.PutCoef(15+iCnt, 28, -Tmp1(iCnt));
       WM.PutCoef(15+iCnt, 29, -Tmp2(iCnt));     
    }
  
    /* Derivata dell'equazione di momento, nodo 1 */
    WM.Add(10, 4, (Mat3x3(MatCross, FPrime) + Mat3x3(MatCrossCross, F, Omega1))*Mat3x3(MatCross, d1Tmp)
                    - MDeltag1
                    - Mat3x3(MatCross, e3a*MPrime(3)));
    WM.Add(10, 10, FWedged1Wedge
                    - Mat3x3(MatCrossCross, MTmp, e3a)
                    - Mat3x3(MatCross, e3a*MPrime(3)));
    WM.Add(10, 16, MDeltag2);
    WM.Add(10, 22, Mat3x3(MatCrossCross, e3a, MTmp));
    WM.Add(10, 25, O1Wedged1Wedge);
    WM.Add(10, 30, Mat3x3(MatCross, d1Tmp));
  
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       WM.PutCoef(9+iCnt, 28, TmpPrime1(iCnt));
       WM.PutCoef(9+iCnt, 29, TmpPrime2(iCnt));
       WM.PutCoef(9+iCnt, 33, Tmp1(iCnt));
       WM.PutCoef(9+iCnt, 34, Tmp2(iCnt));       
       
       /* Contributo del momento attorno all'asse 3, dovuto alla velocita' 
        * imposta */
       WM.PutCoef(9+iCnt, 35, e3a(iCnt));        
    }
    
    /* Derivata dell'equazione di momento, nodo 2 */
    WM.Add(22, 4, MDeltag1 + Mat3x3(MatCross, e3a*MPrime(3)));
    WM.Add(22, 10, Mat3x3(MatCrossCross, MTmp, e3a) + Mat3x3(MatCross, e3a*MPrime(3)));
    WM.Sub(22, 16, (Mat3x3(MatCross, FPrime) + Mat3x3(MatCrossCross, F, Omega2))*Mat3x3(MatCross, d2Tmp) + MDeltag2);
    WM.Add(22, 22, FWedged2Wedge - Mat3x3(MatCrossCross, e3a, MTmp));
    WM.Add(22, 25, O2Wedged2Wedge);
    WM.Sub(22, 30, Mat3x3(MatCross, d2Tmp));
 
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       WM.PutCoef(21+iCnt, 28, -TmpPrime1(iCnt));
       WM.PutCoef(21+iCnt, 29, -TmpPrime2(iCnt));
       WM.PutCoef(21+iCnt, 33, -Tmp1(iCnt));
       WM.PutCoef(21+iCnt, 34, -Tmp2(iCnt));
       
       /* Contributo del momento attorno all'asse 3, dovuto alla velocita' 
        * imposta */
       WM.PutCoef(21+iCnt, 35, -e3a(iCnt));
    }
    
    /* Equazione di vincolo di posizione */
    WM.Add(25, 4, Mat3x3(MatCross, d1Tmp));
    WM.Sub(25, 16, Mat3x3(MatCross, d2Tmp));
    
    /* Derivata dell'equazione di vincolo di posizione */
    WM.Add(30, 4, O1Wedged1Wedge);
    WM.Add(30, 10, Mat3x3(MatCross, d1Tmp));
    WM.Add(30, 16, O2Wedged2Wedge);
    WM.Sub(30, 22, Mat3x3(MatCross, d2Tmp));
    
    /* Equazioni di vincolo di rotazione: e1b~e3a, e2b~e3a */
             
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       doublereal d = Tmp1(iCnt);
       WM.PutCoef(28, 3+iCnt, d);
       WM.PutCoef(28, 15+iCnt, -d);
       
       /* Queste sono per la derivata dell'equazione, sono qui solo per 
        * ottimizzazione */
       WM.PutCoef(33, 9+iCnt, d);
       WM.PutCoef(33, 21+iCnt, -d);
       
       d = Tmp2.dGet(iCnt);
       WM.PutCoef(29, 3+iCnt, -d);
       WM.PutCoef(29, 15+iCnt, d);
       
       /* Queste sono per la derivata dell'equazione, sono qui solo per 
        * ottimizzazione */
       WM.PutCoef(34, 9+iCnt, -d);
       WM.PutCoef(34, 21+iCnt, d);
    }   
    
    /* Derivate delle equazioni di vincolo di rotazione: e1b~e3a, e2b~e3a */
    Vec3 O1mO2(Omega1 - Omega2);
    TmpPrime1 = e3a.Cross(O1mO2.Cross(e2b));   
    TmpPrime2 = e2b.Cross(e3a.Cross(O1mO2));
    for (int iCnt = 1; iCnt <= 3; iCnt++) {      
       WM.PutCoef(33, 3+iCnt, TmpPrime1(iCnt));
       WM.PutCoef(33, 15+iCnt, TmpPrime2(iCnt));
    }
    
    TmpPrime1 = e3a.Cross(O1mO2.Cross(e1b));
    TmpPrime2 = e1b.Cross(e3a.Cross(O1mO2));
    for (int iCnt = 1; iCnt <= 3; iCnt++) {      
       WM.PutCoef(34, 3+iCnt, TmpPrime1(iCnt));
       WM.PutCoef(34, 15+iCnt, TmpPrime2(iCnt));
    }
    
    /* Equazione di vincolo di velocita' di rotazione; viene posta qui perche'
     * a questo numero di equazione corrisponde il numero della 
     * relativa incognita */
    Vec3 Tmp = O1mO2.Cross(e3a);
    for (int iCnt = 1; iCnt <= 3; iCnt++) {
       WM.PutCoef(35, 3+iCnt, Tmp(iCnt));
       doublereal d = e3a(iCnt);
       WM.PutCoef(35, 9+iCnt, -d);
       WM.PutCoef(35, 21+iCnt, d);
    }
    
    return WorkMat;
 }

Here is the call graph for this function:

SubVectorHandler & AxialRotationJoint::InitialAssRes	(	SubVectorHandler &	WorkVec,
		const VectorHandler &	XCurr
	)

virtual

Implements SubjectToInitialAssembly.

Definition at line 3773 of file planej.cc.

References VectorHandler::Add(), Vec3::Cross(), grad::Cross(), d1, d2, DEBUGCOUT, Vec3::dGet(), DriveCaller::dGet(), F, StructNode::GetRCurr(), StructDispNode::GetVCurr(), Mat3x3::GetVec(), StructNode::GetWCurr(), StructDispNode::GetXCurr(), DofOwnerOwner::iGetFirstIndex(), StructDispNode::iGetFirstPositionIndex(), InitialWorkSpaceDim(), M, DriveOwner::pGetDriveCaller(), pNode1, pNode2, VectorHandler::PutCoef(), SubVectorHandler::PutRowIndex(), R1h, R2h, VectorHandler::ResizeReset(), and VectorHandler::Sub().

 {   
    DEBUGCOUT("Entering AxialRotationJoint::InitialAssRes()" << std::endl);
    
    /* Dimensiona e resetta la matrice di lavoro */
    integer iNumRows = 0;
    integer iNumCols = 0;
    InitialWorkSpaceDim(&iNumRows, &iNumCols);
    WorkVec.ResizeReset(iNumRows);
 
    /* Indici */
    integer iNode1FirstPosIndex = pNode1->iGetFirstPositionIndex();
    integer iNode1FirstVelIndex = iNode1FirstPosIndex+6;
    integer iNode2FirstPosIndex = pNode2->iGetFirstPositionIndex();
    integer iNode2FirstVelIndex = iNode2FirstPosIndex+6;
    integer iFirstReactionIndex = iGetFirstIndex();
    integer iReactionPrimeIndex = iFirstReactionIndex+5;
    
    /* Setta gli indici */
    for (int iCnt = 1; iCnt <= 6; iCnt++) {      
       WorkVec.PutRowIndex(iCnt, iNode1FirstPosIndex+iCnt);
       WorkVec.PutRowIndex(6+iCnt, iNode1FirstVelIndex+iCnt);
       WorkVec.PutRowIndex(12+iCnt, iNode2FirstPosIndex+iCnt);
       WorkVec.PutRowIndex(18+iCnt, iNode2FirstVelIndex+iCnt);
    }
    
    for (int iCnt = 1; iCnt <= 5; iCnt++) {
       WorkVec.PutRowIndex(24+iCnt, iFirstReactionIndex+iCnt);
    }
 
    for (int iCnt = 1; iCnt <= 6; iCnt++) {
       WorkVec.PutRowIndex(29+iCnt, iReactionPrimeIndex+iCnt);
    }
 
    /* Recupera i dati */
    const Vec3& x1(pNode1->GetXCurr());
    const Vec3& x2(pNode2->GetXCurr());
    const Vec3& v1(pNode1->GetVCurr());
    const Vec3& v2(pNode2->GetVCurr());
    const Mat3x3& R1(pNode1->GetRCurr());
    const Mat3x3& R2(pNode2->GetRCurr());
    const Vec3& Omega1(pNode1->GetWCurr());
    const Vec3& Omega2(pNode2->GetWCurr());
    
    /* Aggiorna F ed M, che restano anche per InitialAssJac */
    F = Vec3(XCurr, iFirstReactionIndex+1);
    M = Vec3(XCurr(iFirstReactionIndex+4),
             XCurr(iFirstReactionIndex+5),
             0.);
    Vec3 FPrime(XCurr, iReactionPrimeIndex+1);
    Vec3 MPrime(XCurr, iReactionPrimeIndex+4);
    
    /* Distanza nel sistema globale */
    Vec3 d1Tmp(R1*d1);
    Vec3 d2Tmp(R2*d2);
    Mat3x3 R1hTmp(R1*R1h);
    Mat3x3 R2hTmp(R2*R2h);
 
    /* Vettori omega1/\d1, -omega2/\d2 */
    Vec3 O1Wedged1(Omega1.Cross(d1Tmp));
    Vec3 O2Wedged2(Omega2.Cross(d2Tmp));
    
    Vec3 e3a(R1hTmp.GetVec(3));
    Vec3 e1b(R2hTmp.GetVec(1));
    Vec3 e2b(R2hTmp.GetVec(2));  
    
    /* Ruota il momento e la sua derivata con le matrici della cerniera 
     * rispetto ai nodi */
    Vec3 MTmp(e2b*M.dGet(1)-e1b*M.dGet(2));       
    Vec3 MPrimeTmp(e3a.Cross(MTmp.Cross(Omega2))+MTmp.Cross(Omega1.Cross(e3a))+
                   e2b.Cross(e3a)*MPrime.dGet(1)+e3a.Cross(e1b)*MPrime.dGet(2)); 
    
    /* Equazioni di equilibrio, nodo 1 */
    WorkVec.Sub(1, F);
    WorkVec.Add(4, F.Cross(d1Tmp)-MTmp.Cross(e3a)); /* Sfrutto il fatto che F/\d = -d/\F */
    
    /* Derivate delle equazioni di equilibrio, nodo 1 */
    WorkVec.Sub(7, FPrime);
    WorkVec.Add(10, FPrime.Cross(d1Tmp)-O1Wedged1.Cross(F)-MPrimeTmp-
                e3a*MPrime.dGet(3));
    
    /* Equazioni di equilibrio, nodo 2 */
    WorkVec.Add(13, F);
    WorkVec.Add(16, d2Tmp.Cross(F)+MTmp.Cross(e3a));
    
    /* Derivate delle equazioni di equilibrio, nodo 2 */
    WorkVec.Add(19, FPrime);
    WorkVec.Add(22, d2Tmp.Cross(FPrime)+O2Wedged2.Cross(F)+MPrimeTmp+
                e3a*MPrime.dGet(3));
    
    /* Equazione di vincolo di posizione */
    WorkVec.Add(25, x1+d1Tmp-x2-d2Tmp);
    
    /* Derivata dell'equazione di vincolo di posizione */
    WorkVec.Add(30, v1+O1Wedged1-v2-O2Wedged2);
 
    /* Equazioni di vincolo di rotazione */
    WorkVec.PutCoef(28, e2b.Dot(e3a));
    WorkVec.PutCoef(29, e1b.Dot(e3a));
    
    /* Derivate delle equazioni di vincolo di rotazione: e1b~e3a, e2b~e3a */
    Vec3 Tmp((Omega1-Omega2).Cross(e3a));
    WorkVec.PutCoef(33, e2b.Dot(Tmp));
    WorkVec.PutCoef(34, e1b.Dot(Tmp));
 
    /* Equazione di vincolo di velocita' di rotazione */
    doublereal Omega0 = pGetDriveCaller()->dGet();
    WorkVec.PutCoef(35, Omega0-e3a.Dot(Omega2-Omega1));
 
    return WorkVec;
 }

Here is the call graph for this function:

virtual void AxialRotationJoint::InitialWorkSpaceDim	(	integer *	piNumRows,
		integer *	piNumCols
	)		const

inlinevirtual

Implements SubjectToInitialAssembly.

Definition at line 484 of file planej.h.

Referenced by InitialAssJac(), and InitialAssRes().

                                                               {
       *piNumRows = 35; 
       *piNumCols = 35; 
    };

void AxialRotationJoint::Output ( OutputHandler & OH ) const

virtual

Reimplemented from ToBeOutput.

Definition at line 3391 of file planej.cc.

References ToBeOutput::bToBeOutput(), DriveOwner::dGet(), dRaDegr, EULER_123, EULER_313, EULER_321, F, BasicFriction::fc(), fc, OutputHandler::GetCurrentStep(), WithLabel::GetLabel(), StructNode::GetRCurr(), StructNode::GetWCurr(), OutputHandler::JOINTS, OutputHandler::Joints(), M, M3, MatR2EulerAngles123(), MatR2EulerAngles313(), MatR2EulerAngles321(), od, ORIENTATION_MATRIX, ORIENTATION_VECTOR, Joint::Output(), Vec3::pGetVec(), pNode1, pNode2, R1h, R2h, OutputHandler::UseNetCDF(), OutputHandler::UseText(), and RotManip::VecRot().

 {
    if (bToBeOutput()) {
       Mat3x3 R2Tmp(pNode2->GetRCurr()*R2h);
       Mat3x3 RTmp((pNode1->GetRCurr()*R1h).MulTM(R2Tmp));
                 Vec3 E;
                 switch (od) {
                 case EULER_123:
                         E = MatR2EulerAngles123(RTmp)*dRaDegr;
                         break;
 
                 case EULER_313:
                         E = MatR2EulerAngles313(RTmp)*dRaDegr;
                         break;
 
                 case EULER_321:
                         E = MatR2EulerAngles321(RTmp)*dRaDegr;
                         break;
 
                 case ORIENTATION_VECTOR:
                         E = RotManip::VecRot(RTmp);
                         break;
 
                 case ORIENTATION_MATRIX:
                         break;
 
                 default:
                         /* impossible */
                         break;
                 }
       Vec3 OmegaTmp(R2Tmp.MulTV(pNode2->GetWCurr()-pNode1->GetWCurr()));
       
 #ifdef USE_NETCDF
                 if (OH.UseNetCDF(OutputHandler::JOINTS)) {
                         Var_F_local->put_rec((R2Tmp.MulTV(F)).pGetVec(), OH.GetCurrentStep());
                         Var_M_local->put_rec(M.pGetVec(), OH.GetCurrentStep());
                         Var_F_global->put_rec(F.pGetVec(), OH.GetCurrentStep());
                         Var_M_global->put_rec((R2Tmp*M).pGetVec(), OH.GetCurrentStep());
                         switch (od) {
                         case EULER_123:
                         case EULER_313:
                         case EULER_321:
                         case ORIENTATION_VECTOR:
                                 Var_Phi->put_rec(E.pGetVec(), OH.GetCurrentStep());
                                 break;
 
                         case ORIENTATION_MATRIX:
                                 Var_Phi->put_rec(RTmp.pGetMat(), OH.GetCurrentStep());
                                 break;
 
                         default:
                                 /* impossible */
                                 break;
                         }
                         Var_Omega->put_rec(OmegaTmp.pGetVec(), OH.GetCurrentStep());
 /*
                         if (fc) {
                                         Var_MFR->put_rec(&M3, OH.GetCurrentStep());
                                         Var_MU->put_rec(fc->fc(), OH.GetCurrentStep());
                         }
                         else
                         {
                                 Var_MFR->put_rec(0, OH.GetCurrentStep());
                                 Var_MU->put_rec(0, OH.GetCurrentStep());
                         }
 */
                 }
 #endif // USE_NETCDF
                 if (OH.UseText(OutputHandler::JOINTS)) {
                   std::ostream &of = Joint::Output(OH.Joints(), "AxialRotation", GetLabel(),
                                 R2Tmp.MulTV(F), M, F, R2Tmp*M)
                   << " ";
 
                         switch (od) {
                         case EULER_123:
                         case EULER_313:
                         case EULER_321:
                         case ORIENTATION_VECTOR:
                                 OH.Joints() << E;
                                 break;
 
                         case ORIENTATION_MATRIX:
                                 OH.Joints() << RTmp;
                                 break;
 
                         default:
                                 /* impossible */
                                 break;
                         }
 
                         OH.Joints() << " " << dGet()
                   << " " << OmegaTmp;
                   if (fc) {
                           of << " " << M3 << " " << fc->fc();
                   }
                 of << std::endl;
                 }
    }
 }

Here is the call graph for this function:

void AxialRotationJoint::OutputPrepare ( OutputHandler & OH )

virtual

Reimplemented from ToBeOutput.

Definition at line 3365 of file planej.cc.

References ToBeOutput::bToBeOutput(), OutputHandler::JOINTS, od, Joint::OutputPrepare_int(), and OutputHandler::UseNetCDF().

 {
         if (bToBeOutput()) {
 #ifdef USE_NETCDF
                 if (OH.UseNetCDF(OutputHandler::JOINTS)) {
                         std::string name;
                         OutputPrepare_int("axial rotation", OH, name);
 
                         Var_Phi = OH.CreateRotationVar(name, "", od, "global");
 
                         Var_Omega = OH.CreateVar<Vec3>(name + "Omega", "radian/s",
                                 "local relative angular velocity (x, y, z)");
 
 /* TODO
                         Var_MFR = OH.CreateVar<doublereal>(name + "MFR", "Nm",
                                 "friciton moment ");
 
                         Var_MU = OH.CreateVar<doublereal>(name + "MU", "--",
                                         "friction model specific data: friction coefficient?");
 */
                 }
 #endif // USE_NETCDF
         }
 }

Here is the call graph for this function:

Hint * AxialRotationJoint::ParseHint	(	DataManager *	pDM,
		const char *	s
	)		const

virtual

Reimplemented from SimulationEntity.

Definition at line 2912 of file planej.cc.

References fc, SimulationEntity::ParseHint(), and STRLENOF.

 {
         if (strncasecmp(s, "offset{" /*}*/, STRLENOF("offset{" /*}*/)) == 0) {
                 s += STRLENOF("offset{" /*}*/);
 
                 if (strcmp(&s[1], /*{*/ "}") != 0) {
                         return 0;
                 }
 
                 switch (s[0]) {
                 case '1':
                         return new Joint::OffsetHint<1>;
 
                 case '2':
                         return new Joint::OffsetHint<2>;
                 }
 
         } else if (strncasecmp(s, "hinge{" /*}*/, STRLENOF("hinge{" /*}*/)) == 0) {
                 s += STRLENOF("hinge{" /*}*/);
 
                 if (strcmp(&s[1], /*{*/ "}") != 0) {
                         return 0;
                 }
 
                 switch (s[0]) {
                 case '1':
                         return new Joint::HingeHint<1>;
 
                 case '2':
                         return new Joint::HingeHint<2>;
                 }
 
         } else if (fc) {
                 return fc->ParseHint(pDM, s);
         }
 
         return 0;
 }

Here is the call graph for this function:

std::ostream & AxialRotationJoint::Restart ( std::ostream & out ) const

virtual

Implements Elem.

Definition at line 2989 of file planej.cc.

References d1, d2, WithLabel::GetLabel(), Mat3x3::GetVec(), DriveOwner::pGetDriveCaller(), pNode1, pNode2, R1h, R2h, Joint::Restart(), DriveCaller::Restart(), Write(), and Vec3::Write().

 {
    Joint::Restart(out) << ", axial rotation, "
      << pNode1->GetLabel() 
      << ", reference, node, ", d1.Write(out, ", ")
      << ", hinge, reference, node, 1, ", (R1h.GetVec(1)).Write(out, ", ")
      << ", 2, ", (R1h.GetVec(2)).Write(out, ", ") << ", "
      << pNode2->GetLabel() 
      << ", reference, node, ", d2.Write(out, ", ") 
      << ", hinge, reference, node, 1, ", (R2h.GetVec(1)).Write(out, ", ")
      << ", 2, ", (R2h.GetVec(2)).Write(out, ", ") << ", ";
    
    return pGetDriveCaller()->Restart(out) << ';' << std::endl;
 }

Here is the call graph for this function:

void AxialRotationJoint::SetValue	(	DataManager *	pDM,
		VectorHandler &	X,
		VectorHandler &	XP,
		SimulationEntity::Hints *	ph = `0`
	)

virtual

Reimplemented from Joint.

Definition at line 2853 of file planej.cc.

References d1, d2, Vec3::dGet(), dTheta, dThetaWrapped, fc, StructNode::GetRCurr(), StructDispNode::GetXCurr(), DofOwnerOwner::iGetFirstIndex(), MatR2EulerAngles(), MBDYN_EXCEPT_ARGS, Mat3x3::MulTM(), NumSelfDof, DriveHint::pCreateDrive(), pNode1, pNode2, R1h, R2h, DriveOwner::Set(), and BasicFriction::SetValue().

 {
         if (ph) {
                 for (unsigned i = 0; i < ph->size(); i++) {
                         Joint::JointHint *pjh = dynamic_cast<Joint::JointHint *>((*ph)[i]);
 
                         if (pjh == 0) {
                                 DriveHint *pdh = dynamic_cast<DriveHint *>((*ph)[i]);
                                 if (pdh) {
                                         DriveCaller *pDC = pdh->pCreateDrive(pDM);
                                         if (pDC == 0) {
                                                 silent_cerr("AxialRotationJoint::SetValue: "
                                                         "unable to create drive after hint "
                                                         "#" << i << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                         DriveOwner::Set(pDC);
                                 }
                                 continue;
                         }
 
                         if (dynamic_cast<Joint::OffsetHint<1> *>(pjh)) {
                                 const Mat3x3& R1(pNode1->GetRCurr());
                                 Vec3 dTmp2(pNode2->GetRCurr()*d2);
    
                                 d1 = R1.MulTV(pNode2->GetXCurr() + dTmp2 - pNode1->GetXCurr());
 
                         } else if (dynamic_cast<Joint::OffsetHint<2> *>(pjh)) {
                                 const Mat3x3& R2(pNode2->GetRCurr());
                                 Vec3 dTmp1(pNode1->GetRCurr()*d1);
    
                                 d2 = R2.MulTV(pNode1->GetXCurr() + dTmp1 - pNode2->GetXCurr());
 
                         } else if (dynamic_cast<Joint::HingeHint<1> *>(pjh)) {
                                 R1h = pNode1->GetRCurr().MulTM(pNode2->GetRCurr()*R2h);
 
                         } else if (dynamic_cast<Joint::HingeHint<2> *>(pjh)) {
                                 R2h = pNode2->GetRCurr().MulTM(pNode1->GetRCurr()*R1h);
 
                         } else if (dynamic_cast<Joint::ReactionsHint *>(pjh)) {
                                 /* TODO */
                         }
                 }
         }
 
         Mat3x3 RTmp((pNode1->GetRCurr()*R1h).MulTM(pNode2->GetRCurr()*R2h));
         Vec3 v(MatR2EulerAngles(RTmp));
 
         dThetaWrapped = dTheta = v.dGet(3);
         
         if (fc) {
                 fc->SetValue(pDM, X, XP, ph, iGetFirstIndex() + NumSelfDof);
         }
 }

Here is the call graph for this function:

void AxialRotationJoint::WorkSpaceDim	(	integer *	piNumRows,
		integer *	piNumCols
	)		const

inlinevirtual

Implements Elem.

Definition at line 454 of file planej.h.

References fc, SimulationEntity::iGetNumDof(), and NumDof.

Referenced by AssJac(), and AssRes().

                                                                    { 
       *piNumRows = NumDof;
       *piNumCols = NumDof;
       if (fc) {
           *piNumRows += fc->iGetNumDof();
           *piNumCols += fc->iGetNumDof();
       } 
    };