Standalone tools for beam-based finite element method (FEM) No dependency with SubDyn types and representation. More...

Data Types
interface	findloci

Functions/Subroutines
subroutine	eigensolve (K, M, N, bCheckSingularity, EigVect, Omega, ErrStat, ErrMsg)
	Return eigenvalues, Omega, and eigenvectors. More...

pure subroutine	sort_in_place (a, key)

real(feki) function	determinant (A, ErrStat, ErrMsg)
	Compute the determinant of a real matrix using an LU factorization. More...

subroutine	chessboard (M, blackVal, whiteVal, nSpace, diagVal)
	Create a chessboard-like matrix with `valBlack` on the "black" cases, starting with black at (1,1) As a generalization, "black" values may be spaced every `nSpace` squares For instance, blackVal=9, whiteVal=0, nSpace=2 [9 0 0 9 0 0 9] [0 9 0 0 9 0 0] [0 0 9 0 0 9 0] Diagonal values may be overriden by `diagVal` Matrix M does not need to be square. More...

subroutine	breaksysmtrx (MM, KK, IDR, IDL, nR, nL, MRR, MLL, MRL, KRR, KLL, KRL, FG, FGR, FGL, CC, CRR, CLL, CRL)
	Partition matrices and vectors into Boundary (R) and internal (L) nodes M = [ MRR, MRL ] [ sym, MLL ] MRR = M(IDR, IDR), KRR = M(IDR, IDR), FR = F(IDR) MLL = M(IDL, IDL), KRR = K(IDL, IDL), FL = F(IDL) MRL = M(IDR, IDL), KRR = K(IDR, IDL) NOTE: generic code. More...

subroutine	craigbamptonreduction (MM, KK, IDR, nR, IDL, nL, nM, nM_Out, MBB, MBM, KBB, PhiL, PhiR, OmegaL, ErrStat, ErrMsg, FG, FGR, FGL, FGB, FGM, CC, CBB, CBM, CMM)
	Performs Craig-Bampton reduction of M and K matrices and optional Force vector TODO: (Damping optional) Convention is: "R": leader DOF -> "B": reduced leader DOF "L": interior DOF -> "M": reduced interior DOF (CB-modes) NOTE: More...

subroutine	craigbamptonreduction_frompartition (MRR, MLL, MRL, KRR, KLL, KRL, nR, nL, nM, nM_Out, MBB, MBM, KBB, PhiL, PhiR, OmegaL, ErrStat, ErrMsg, CRR, CLL, CRL, CBB, CBM, CMM)
	Performs Craig-Bampton reduction based on partitioned matrices M and K Convention is: "R": leader DOF -> "B": reduced leader DOF "L": interior DOF -> "M": reduced interior DOF (CB-modes) NOTE: More...

subroutine	eigensolvewrap (K, M, nDOF, NOmega, bCheckSingularity, EigVect, Omega, ErrStat, ErrMsg, bDOF)
	Wrapper function for eigen value analyses, for two cases: Case1: K and M are taken "as is", this is used for the "LL" part of the matrix Case2: K and M contain some constraints lines, and they need to be removed from the Mass/Stiffness matrix. More...

subroutine	removedof (A, bDOF, Ared, ErrStat, ErrMsg)
	Remove degrees of freedom from a matrix (lines and rows)

subroutine	insertdofrows (Ared, bDOF, DefaultVal, A, ErrStat, ErrMsg)
	Expand a matrix to includes rows where bDOF is False (inverse behavior as RemoveDOF)

integer(intki) function	findloci_reki (Array, Val)
	Returns index of val in Array (val is an integer!) More...

integer(intki) function	findloci_intki (Array, Val)
	Returns index of val in Array (val is an integer!) More...

subroutine	rigidtransformationline (dx, dy, dz, iLine, Line)

subroutine	getrigidtransformation (Pj, Pk, TRigid, ErrStat, ErrMsg)
	Rigid transformation matrix between DOFs of node j and k where node j is the leader node. More...

subroutine	getdircos (P1, P2, DirCos, L_out, ErrStat, ErrMsg)
	Computes directional cosine matrix DirCos Transforms from element to global coordinates: xg = DC.xe, Kg = DC.Ke.DC^t Assumes that the element main direction is along ze. More...

subroutine	getorthvectors (e1, e2, e3, ErrStat, ErrMsg)
	Returns two vectors orthonormal to the input vector.

subroutine	elemk_beam (A, L, Ixx, Iyy, Jzz, Shear, kappa, E, G, DirCos, K)
	Element stiffness matrix for classical beam elements shear is true – non-tapered Timoshenko beam shear is false – non-tapered Euler-Bernoulli beam. More...

subroutine	elemk_cable (A, L, E, T0, DirCos, K)
	Element stiffness matrix for pretension cable Element coordinate system: z along the cable! More...

subroutine	elemm_beam (A, L, Ixx, Iyy, Jzz, rho, DirCos, M)
	Element mass matrix for classical beam elements. More...

subroutine	elemm_cable (A, L, rho, DirCos, M)
	Element stiffness matrix for pretension cable. More...

subroutine	elemg (A, L, rho, DirCos, F, g)
	calculates the lumped forces and moments due to gravity on a given element: the element has two nodes, with the loads for both elements stored in array F. More...

subroutine	elemf_cable (T0, DirCos, F)

subroutine	lumpforces (Area1, Area2, crat, L, rho, g, DirCos, F)
	Calculates the lumped gravity forces at the nodes given the element geometry It assumes a linear variation of the dimensions from node 1 to node 2, thus the area may be quadratically varying if crat<>1 bjj: note this routine is a work in progress, intended for future version of SubDyn. More...

subroutine	pseudoinverse (A, Ainv, ErrStat, ErrMsg)
	Method 1: pinv_A = A \ eye(m) (matlab) call _GELSS to solve A.X=B pinv(A) = B(1:n,1:m) Method 2: [U,S,V] = svd(A); pinv_A = ( V / S ) * U'; (matlab)

Variables
integer, parameter	feki = R8Ki

integer, parameter	laki = R8Ki

Detailed Description

Standalone tools for beam-based finite element method (FEM) No dependency with SubDyn types and representation.

Function/Subroutine Documentation

◆ breaksysmtrx()

subroutine fem::breaksysmtrx	(	real(feki), dimension(:,:), intent(in)	MM,
		real(feki), dimension(:,:), intent(in)	KK,
		integer(intki), dimension(nr), intent(in)	IDR,
		integer(intki), dimension(nl), intent(in)	IDL,
		integer(intki), intent(in)	nR,
		integer(intki), intent(in)	nL,
		real(feki), dimension(nr, nr), intent(out)	MRR,
		real(feki), dimension(nl, nl), intent(out)	MLL,
		real(feki), dimension(nr, nl), intent(out)	MRL,
		real(feki), dimension(nr, nr), intent(out)	KRR,
		real(feki), dimension(nl, nl), intent(out)	KLL,
		real(feki), dimension(nr, nl), intent(out)	KRL,
		real(feki), dimension(:), intent(in), optional	FG,
		real(feki), dimension(nr), intent(out), optional	FGR,
		real(feki), dimension(nl), intent(out), optional	FGL,
		real(feki), dimension(:,:), intent(in), optional	CC,
		real(feki), dimension(nr, nr), intent(out), optional	CRR,
		real(feki), dimension(nl, nl), intent(out), optional	CLL,
		real(feki), dimension(nr, nl), intent(out), optional	CRL
	)

Partition matrices and vectors into Boundary (R) and internal (L) nodes M = [ MRR, MRL ] [ sym, MLL ] MRR = M(IDR, IDR), KRR = M(IDR, IDR), FR = F(IDR) MLL = M(IDL, IDL), KRR = K(IDL, IDL), FL = F(IDL) MRL = M(IDR, IDL), KRR = K(IDR, IDL) NOTE: generic code.

Parameters

[in]	mm	Mass Matrix
[in]	kk	Stiffness matrix
[in]	idr	Indices of leader DOFs
[in]	idl	Indices of interior DOFs
[in]	fg	Force vector
[in]	cc	Stiffness matrix

◆ chessboard()

subroutine fem::chessboard	(	real(reki), dimension(:,:), intent(out)	M,
		real(reki), intent(in)	blackVal,
		real(reki), intent(in)	whiteVal,
		integer(intki), intent(in), optional	nSpace,
		real(reki), intent(in), optional	diagVal
	)

Create a chessboard-like matrix with valBlack on the "black" cases, starting with black at (1,1) As a generalization, "black" values may be spaced every nSpace squares For instance, blackVal=9, whiteVal=0, nSpace=2 [9 0 0 9 0 0 9] [0 9 0 0 9 0 0] [0 0 9 0 0 9 0] Diagonal values may be overriden by diagVal Matrix M does not need to be square.

Parameters

[out]	m	Output matrix
[in]	blackval	value for black squares
[in]	whiteval	value for white squre
[in]	nspace	spacing between black values, default 1
[in]	diagval	Value to override diagonal

◆ craigbamptonreduction()

subroutine fem::craigbamptonreduction	(	real(feki), dimension(:, :), intent(in)	MM,
		real(feki), dimension(:, :), intent(in)	KK,
		integer(intki), dimension(nr), intent(in)	IDR,
		integer(intki), intent(in)	nR,
		integer(intki), dimension(nl), intent(in)	IDL,
		integer(intki), intent(in)	nL,
		integer(intki), intent(in)	nM,
		integer(intki), intent(in)	nM_Out,
		real(feki), dimension( nr, nr), intent(out)	MBB,
		real(feki), dimension( nr, nm), intent(out)	MBM,
		real(feki), dimension( nr, nr), intent(out)	KBB,
		real(feki), dimension(nl, nm_out), intent(out)	PhiL,
		real(feki), dimension(nl, nr), intent(out)	PhiR,
		real(feki), dimension(nm_out), intent(out)	OmegaL,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg,
		real(feki), dimension(:), intent(in), optional	FG,
		real(feki), dimension(nr), intent(out), optional	FGR,
		real(feki), dimension(nl), intent(out), optional	FGL,
		real(feki), dimension(nr), intent(out), optional	FGB,
		real(feki), dimension(nm), intent(out), optional	FGM,
		real(feki), dimension(:, :), intent(in), optional	CC,
		real(feki), dimension(nr, nr), intent(out), optional	CBB,
		real(feki), dimension(nr, nm), intent(out), optional	CBM,
		real(feki), dimension(nm, nm), intent(out), optional	CMM
	)

Performs Craig-Bampton reduction of M and K matrices and optional Force vector TODO: (Damping optional) Convention is: "R": leader DOF -> "B": reduced leader DOF "L": interior DOF -> "M": reduced interior DOF (CB-modes) NOTE:

M_MM = Identity and K_MM = Omega*2 hence these matrices are not returned
Possibility to get more CB modes using the input nM_Out>nM

NOTE: generic code

Parameters

[in]	mm	Mass matrix
[in]	kk	Stiffness matrix
[in]	idr	Indices of leader DOFs
[in]	idl	Indices of interior DOFs
[in]	nm	Number of CB modes
[in]	nm_out	Number of modes returned for PhiL & OmegaL
[out]	mbb	Reduced Guyan Mass Matrix
[out]	kbb	Reduced Guyan Stiffness matrix
[out]	mbm	Cross term
[out]	phir	Guyan Modes
[out]	phil	Craig-Bampton modes
[out]	omegal	Eigenvalues
[in]	fg	Force vector (typically a constant force, like gravity)
[out]	fgr	Force vector partitioned for R DOFs (TODO remove me)
[out]	fgl	Force vector partitioned for L DOFs (TODO somehow for Static improvment..)
[out]	fgb	Force vector in Guyan modes = FR+PhiR^t FL
[out]	fgm	Force vector in CB modes = PhiM^t FL
[in]	cc	Damping matrix
[out]	cbb	Guyan Damping matrix
[out]	cbm	Coupling Damping matrix
[out]	cmm	Craig-Bampton Damping matrix

◆ craigbamptonreduction_frompartition()

subroutine fem::craigbamptonreduction_frompartition	(	real(feki), dimension( nr, nr), intent(in)	MRR,
		real(feki), dimension( nl, nl), intent(in)	MLL,
		real(feki), dimension( nr, nl), intent(in)	MRL,
		real(feki), dimension( nr, nr), intent(in)	KRR,
		real(feki), dimension( nl, nl), intent(inout)	KLL,
		real(feki), dimension( nr, nl), intent(in)	KRL,
		integer(intki), intent(in)	nR,
		integer(intki), intent(in)	nL,
		integer(intki), intent(in)	nM,
		integer(intki), intent(in)	nM_Out,
		real(feki), dimension( nr, nr), intent(out)	MBB,
		real(feki), dimension( nr, nm), intent(out)	MBM,
		real(feki), dimension( nr, nr), intent(out)	KBB,
		real(feki), dimension(nl, nm_out), intent(out)	PhiL,
		real(feki), dimension(nl, nr), intent(out)	PhiR,
		real(feki), dimension(nm_out), intent(out)	OmegaL,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg,
		real(feki), dimension( nr, nr), intent(in), optional	CRR,
		real(feki), dimension( nl, nl), intent(in), optional	CLL,
		real(feki), dimension( nr, nl), intent(in), optional	CRL,
		real(feki), dimension( nr, nr), intent(out), optional	CBB,
		real(feki), dimension( nr, nm), intent(out), optional	CBM,
		real(feki), dimension( nm, nm), intent(out), optional	CMM
	)

Performs Craig-Bampton reduction based on partitioned matrices M and K Convention is: "R": leader DOF -> "B": reduced leader DOF "L": interior DOF -> "M": reduced interior DOF (CB-modes) NOTE:

M_MM = Identity and K_MM = Omega*2 hence these matrices are not returned
Possibility to get more CB modes using the input nM_Out>nM (e.g. for static improvement)

NOTE: generic code

Parameters

[in]	mrr	Partitioned mass and stiffness matrices
[out]	phir	Guyan Modes
[out]	phil	Craig-Bampton modes
[out]	omegal	Eigenvalues
[out]	errstat	Error status of the operation
[out]	errmsg	Error message if ErrStat /= ErrID_None
[in]	crr	Partitioned damping matrices
[out]	cbb	Guyan damping matrix
[out]	cbm	Coupling damping matrix
[out]	cmm	CB damping matrix

◆ determinant()

real(feki) function fem::determinant	(	real(feki), dimension(:, :), intent(in)	A,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg
	)

Compute the determinant of a real matrix using an LU factorization.

Parameters

[in]	a	Input matrix, no side effect
[out]	errstat	Error status of the operation
[out]	errmsg	Error message if ErrStat /= ErrID_None

Returns: May easily overflow

◆ eigensolve()

subroutine fem::eigensolve	(	real(laki), dimension(n, n), intent(inout)	K,
		real(laki), dimension(n, n), intent(inout)	M,
		integer, intent(in)	N,
		logical, intent(in)	bCheckSingularity,
		real(laki), dimension(n, n), intent(inout)	EigVect,
		real(laki), dimension(n), intent(inout)	Omega,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg
	)

Return eigenvalues, Omega, and eigenvectors.

Parameters

[in]	n	Number of degrees of freedom, size of M and K
[in,out]	k	Stiffness matrix
[in,out]	m	Mass matrix
[in,out]	eigvect	Returned Eigenvectors
[in,out]	omega	Returned Eigenvalues
[out]	errstat	Error status of the operation
[out]	errmsg	Error message if ErrStat /= ErrID_None

◆ eigensolvewrap()

subroutine fem::eigensolvewrap	(	real(feki), dimension(ndof, ndof), intent(in)	K,
		real(feki), dimension(ndof, ndof), intent(in)	M,
		integer, intent(in)	nDOF,
		integer, intent(in)	NOmega,
		logical, intent(in)	bCheckSingularity,
		real(feki), dimension(ndof, nomega), intent(out)	EigVect,
		real(feki), dimension(nomega), intent(out)	Omega,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg,
		logical, dimension(ndof), intent(in), optional	bDOF
	)

Wrapper function for eigen value analyses, for two cases: Case1: K and M are taken "as is", this is used for the "LL" part of the matrix Case2: K and M contain some constraints lines, and they need to be removed from the Mass/Stiffness matrix.

Used for full system

◆ elemf_cable()

subroutine fem::elemf_cable	(	real(reki), intent(in)	T0,
		real(feki), dimension(3,3), intent(in)	DirCos,
		real(feki), dimension(12), intent(out)	F
	)

Parameters

[in]	t0	Pretension load [N]
[in]	dircos	From element to global: xg = DC.xe, Kg = DC.Ke.DC^t
[out]	f	returned loads. 1-6 for node 1; 7-12 for node 2.

◆ elemg()

subroutine fem::elemg	(	real(reki), intent(in)	A,
		real(reki), intent(in)	L,
		real(reki), intent(in)	rho,
		real(feki), dimension(3,3), intent(in)	DirCos,
		real(feki), dimension(12), intent(out)	F,
		real(reki), intent(in)	g
	)

calculates the lumped forces and moments due to gravity on a given element: the element has two nodes, with the loads for both elements stored in array F.

Indexing of F is: Fx_n1=1,Fy_n1=2,Fz_n1=3,Mx_n1= 4,My_n1= 5,Mz_n1= 6, Fx_n2=7,Fy_n2=8,Fz_n2=9,Mx_n2=10,My_n2=11,Mz_n2=12

Parameters

[in]	a	area
[in]	l	element length
[in]	rho	density
[in]	dircos	From element to global: xg = DC.xe, Kg = DC.Ke.DC^t
[in]	g	gravity
[out]	f	returned loads. positions 1-6 are the loads for node 1 ; 7-12 are loads for node 2.

◆ elemk_beam()

subroutine fem::elemk_beam	(	real(reki), intent(in)	A,
		real(reki), intent(in)	L,
		real(reki), intent(in)	Ixx,
		real(reki), intent(in)	Iyy,
		real(reki), intent(in)	Jzz,
		logical, intent(in)	Shear,
		real(reki), intent(in)	kappa,
		real(reki), intent(in)	E,
		real(reki), intent(in)	G,
		real(feki), dimension(3,3), intent(in)	DirCos,
		real(feki), dimension(12, 12), intent(out)	K
	)

Element stiffness matrix for classical beam elements shear is true – non-tapered Timoshenko beam shear is false – non-tapered Euler-Bernoulli beam.

Parameters

[in] dircos From element to global: xg = DC.xe, Kg = DC.Ke.DC^t

◆ elemk_cable()

subroutine fem::elemk_cable	(	real(reki), intent(in)	A,
		real(reki), intent(in)	L,
		real(reki), intent(in)	E,
		real(reki), intent(in)	T0,
		real(feki), dimension(3,3), intent(in)	DirCos,
		real(feki), dimension(12, 12), intent(out)	K
	)

Element stiffness matrix for pretension cable Element coordinate system: z along the cable!

Parameters

[in] dircos From element to global: xg = DC.xe, Kg = DC.Ke.DC^t

◆ elemm_beam()

subroutine fem::elemm_beam	(	real(reki), intent(in)	A,
		real(reki), intent(in)	L,
		real(reki), intent(in)	Ixx,
		real(reki), intent(in)	Iyy,
		real(reki), intent(in)	Jzz,
		real(reki), intent(in)	rho,
		real(feki), dimension(3,3), intent(in)	DirCos,
		real(feki), dimension(12, 12), intent(out)	M
	)

Element mass matrix for classical beam elements.

Parameters

[in] dircos From element to global: xg = DC.xe, Kg = DC.Ke.DC^t

◆ elemm_cable()

subroutine fem::elemm_cable	(	real(reki), intent(in)	A,
		real(feki), intent(in)	L,
		real(reki), intent(in)	rho,
		real(feki), dimension(3,3), intent(in)	DirCos,
		real(feki), dimension(12, 12), intent(out)	M
	)

Element stiffness matrix for pretension cable.

Parameters

[in] dircos From element to global: xg = DC.xe, Kg = DC.Ke.DC^t

◆ findloci_intki()

integer(intki) function fem::findloci_intki	(	integer(intki), dimension(:), intent(in)	Array,
		integer(intki), intent(in)	Val
	)

Returns index of val in Array (val is an integer!)

Parameters

[in] array Array to search in

Returns: Index of joint in joint table

◆ findloci_reki()

integer(intki) function fem::findloci_reki	(	real(reki), dimension(:), intent(in)	Array,
		integer(intki), intent(in)	Val
	)

Returns index of val in Array (val is an integer!)

Parameters

[in] array Array to search in

Returns: Index of joint in joint table

◆ getdircos()

subroutine fem::getdircos	(	real(reki), dimension(3), intent(in)	P1,
		real(reki), dimension(3), intent(in)	P2,
		real(feki), dimension(3, 3), intent(out)	DirCos,
		real(reki), intent(out)	L_out,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg
	)

Computes directional cosine matrix DirCos Transforms from element to global coordinates: xg = DC.xe, Kg = DC.Ke.DC^t Assumes that the element main direction is along ze.

bjj: note that this is the transpose of what is normally considered the Direction Cosine Matrix in the FAST framework.

◆ getrigidtransformation()

subroutine fem::getrigidtransformation	(	real(reki), dimension(3), intent(in)	Pj,
		real(reki), dimension(3), intent(in)	Pk,
		real(reki), dimension(6,6), intent(out)	TRigid,
		integer(intki), intent(out)	ErrStat,
		character(*), intent(out)	ErrMsg
	)

Rigid transformation matrix between DOFs of node j and k where node j is the leader node.

◆ lumpforces()

subroutine fem::lumpforces	(	real(reki), intent(in)	Area1,
		real(reki), intent(in)	Area2,
		real(reki), intent(in)	crat,
		real(reki), intent(in)	L,
		real(reki), intent(in)	rho,
		real(reki), intent(in)	g,
		real(reki), dimension(3,3), intent(in)	DirCos,
		real(reki), dimension(12), intent(out)	F
	)

Calculates the lumped gravity forces at the nodes given the element geometry It assumes a linear variation of the dimensions from node 1 to node 2, thus the area may be quadratically varying if crat<>1 bjj: note this routine is a work in progress, intended for future version of SubDyn.

Compare with ElemG.

Parameters

[in]	crat	X-sectional areas at node 1 and node 2, t2/t1 thickness ratio
[in]	g	gravity
[in]	l	Length of element
[in]	rho	density
[in]	dircos	From element to global: xg = DC.xe, Kg = DC.Ke.DC^t
[out]	f	Lumped forces

Data Types

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ breaksysmtrx()

◆ chessboard()

◆ craigbamptonreduction()

◆ craigbamptonreduction_frompartition()

◆ determinant()

◆ eigensolve()

◆ eigensolvewrap()

◆ elemf_cable()

◆ elemg()

◆ elemk_beam()

◆ elemk_cable()

◆ elemm_beam()

◆ elemm_cable()

◆ findloci_intki()

◆ findloci_reki()

◆ getdircos()

◆ getrigidtransformation()

◆ lumpforces()