ISAAR.MSolve.IGA/Elements/TSplineKirchhoffLoveShellElementMaterial.cs
/// <summary>
/// Calculates the stiffness matrix of the element.
/// </summary>
/// <param name="element">An element of type <see cref="TSplineKirchhoffLoveShellElementMaterial"/>.</param>
/// <returns>An <see cref="IMatrix"/> containing the stiffness matrix of an <see cref="TSplineKirchhoffLoveShellElementMaterial"/>.</returns>
Calculates the stiffness matrix of the element.
Calculates the stiffness matrix of the element .
public IMatrix StiffnessMatrix(IElement element)
		{
			var shellElement = (TSplineKirchhoffLoveShellElementMaterial)element;
			var gaussPoints = _materialsAtThicknessGp.Keys.ToArray();
			var elementControlPoints = shellElement.ControlPoints.ToArray();
			var bRows = 3;
			var bCols = elementControlPoints.Length * 3;
			var stiffnessMatrix = new double[bCols, bCols];
			var BmTranspose = new double[bCols, bRows];
			var BbTranspose = new double[bCols, bRows];
			var BmTransposeMultStiffness = new double[bCols, bRows];
			var BbTransposeMultStiffness = new double[bCols, bRows];
			var BmbTransposeMultStiffness = new double[bCols, bRows];
			var BbmTransposeMultStiffness = new double[bCols, bRows];
			for (int j = 0; j < gaussPoints.Length; j++)
			{
				var jacobianMatrix = CalculateJacobian(elementControlPoints, _tsplines, j);
				var hessianMatrix = CalculateHessian(elementControlPoints, _tsplines, j);
				var surfaceBasisVector1 = CalculateSurfaceBasisVector1(jacobianMatrix, 0);
				var surfaceBasisVector2 = CalculateSurfaceBasisVector1(jacobianMatrix, 1);
				var surfaceBasisVector3 = new[]
				{
					surfaceBasisVector1[1] * surfaceBasisVector2[2] - surfaceBasisVector1[2] * surfaceBasisVector2[1],
					surfaceBasisVector1[2] * surfaceBasisVector2[0] - surfaceBasisVector1[0] * surfaceBasisVector2[2],
					surfaceBasisVector1[0] * surfaceBasisVector2[1] - surfaceBasisVector1[1] * surfaceBasisVector2[0]
				};
				foreach (var integrationPointMaterial in _materialsAtThicknessGp[gaussPoints[j]].Values)
				{
					integrationPointMaterial.TangentVectorV1 = surfaceBasisVector1;
					integrationPointMaterial.TangentVectorV2 = surfaceBasisVector2;
					integrationPointMaterial.NormalVectorV3 = surfaceBasisVector3;
				}
				double norm = 0;
				for (int i = 0; i < surfaceBasisVector3.Length; i++)
					norm += surfaceBasisVector3[i] * surfaceBasisVector3[i];
				var J1 = Math.Sqrt(norm);
				for (int i = 0; i < surfaceBasisVector3.Length; i++)
					surfaceBasisVector3[i] = surfaceBasisVector3[i] / J1;
				var surfaceBasisVectorDerivative1 = CalculateSurfaceBasisVector1(hessianMatrix, 0);
				var surfaceBasisVectorDerivative2 = CalculateSurfaceBasisVector1(hessianMatrix, 1);
				var surfaceBasisVectorDerivative12 = CalculateSurfaceBasisVector1(hessianMatrix, 2);
				var Bmembrane = CalculateMembraneDeformationMatrix(_tsplines, j, surfaceBasisVector1,
					surfaceBasisVector2, elementControlPoints);
				var Bbending = CalculateBendingDeformationMatrix(surfaceBasisVector3, _tsplines, j, surfaceBasisVector2,
					surfaceBasisVectorDerivative1, surfaceBasisVector1, J1, surfaceBasisVectorDerivative2,
					surfaceBasisVectorDerivative12, elementControlPoints);
				var (MembraneConstitutiveMatrix, BendingConstitutiveMatrix, CouplingConstitutiveMatrix) =
					IntegratedConstitutiveOverThickness(gaussPoints, j);
				double wFactor = J1 * gaussPoints[j].WeightFactor;
				double tempb = 0;
				double tempm = 0;
				Array.Clear(BmTranspose, 0, bRows * bCols);
				Array.Clear(BbTranspose, 0, bRows * bCols);
				for (int i = 0; i < bRows; i++)
				{
					for (int k = 0; k < bCols; k++)
					{
						BmTranspose[k, i] = Bmembrane[i, k] * wFactor;
						BbTranspose[k, i] = Bbending[i, k] * wFactor;
					}
				}
				double tempcm = 0;
				double tempcb = 0;
				double tempcc = 0;
				Array.Clear(BmTransposeMultStiffness, 0, bRows * bCols);
				Array.Clear(BbTransposeMultStiffness, 0, bRows * bCols);
				Array.Clear(BmbTransposeMultStiffness, 0, bRows * bCols);
				Array.Clear(BbmTransposeMultStiffness, 0, bRows * bCols);
				for (int i = 0; i < bCols; i++)
				{
					for (int k = 0; k < bRows; k++)
					{
						tempm = BmTranspose[i, k];
						tempb = BbTranspose[i, k];
						for (int m = 0; m < bRows; m++)
						{
							tempcm = MembraneConstitutiveMatrix[k, m];
							tempcb = BendingConstitutiveMatrix[k, m];
							tempcc = CouplingConstitutiveMatrix[k, m];
							BmTransposeMultStiffness[i, m] += tempm * tempcm;
							BbTransposeMultStiffness[i, m] += tempb * tempcb;
							BmbTransposeMultStiffness[i, m] += tempm * tempcc;
							BbmTransposeMultStiffness[i, m] += tempb * tempcc;
						}
					}
				}
				double tempmb = 0;
				double tempbm = 0;
				double mem = 0;
				double ben = 0;
				for (int i = 0; i < bCols; i++)
				{
					for (int k = 0; k < bRows; k++)
					{
						tempm = BmTransposeMultStiffness[i, k];
						tempb = BbTransposeMultStiffness[i, k];
						tempmb = BmbTransposeMultStiffness[i, k];
						tempbm = BbmTransposeMultStiffness[i, k];
						for (int m = 0; m < bCols; m++)
						{
							mem = Bmembrane[k, m];
							ben = Bbending[k, m];
							stiffnessMatrix[i, m] += tempm * mem + tempb * ben + tempmb * ben + tempbm * mem;
						}
					}
				}
			}
			return Matrix.CreateFromArray(stiffnessMatrix);
		}
public IMatrix StiffnessMatrix ( IElement element ) { var shellElement = ( TSplineKirchhoffLoveShellElementMaterial ) element ; var gaussPoints = _materialsAtThicknessGp . Keys . ToArray ( ) ; var elementControlPoints = shellElement . ControlPoints . ToArray ( ) ; var bRows = 3 ; var bCols = elementControlPoints . Length * 3 ; var stiffnessMatrix = new double [ bCols , bCols ] ; var BmTranspose = new double [ bCols , bRows ] ; var BbTranspose = new double [ bCols , bRows ] ; var BmTransposeMultStiffness = new double [ bCols , bRows ] ; var BbTransposeMultStiffness = new double [ bCols , bRows ] ; var BmbTransposeMultStiffness = new double [ bCols , bRows ] ; var BbmTransposeMultStiffness = new double [ bCols , bRows ] ; for ( int j = 0 ; j < gaussPoints . Length ; j ++ ) { var jacobianMatrix = CalculateJacobian ( elementControlPoints , _tsplines , j ) ; var hessianMatrix = CalculateHessian ( elementControlPoints , _tsplines , j ) ; var surfaceBasisVector1 = CalculateSurfaceBasisVector1 ( jacobianMatrix , 0 ) ; var surfaceBasisVector2 = CalculateSurfaceBasisVector1 ( jacobianMatrix , 1 ) ; var surfaceBasisVector3 = new [ ] { surfaceBasisVector1 [ 1 ] * surfaceBasisVector2 [ 2 ] - surfaceBasisVector1 [ 2 ] * surfaceBasisVector2 [ 1 ] , surfaceBasisVector1 [ 2 ] * surfaceBasisVector2 [ 0 ] - surfaceBasisVector1 [ 0 ] * surfaceBasisVector2 [ 2 ] , surfaceBasisVector1 [ 0 ] * surfaceBasisVector2 [ 1 ] - surfaceBasisVector1 [ 1 ] * surfaceBasisVector2 [ 0 ] } ; foreach ( var integrationPointMaterial in _materialsAtThicknessGp [ gaussPoints [ j ] ] . Values ) { integrationPointMaterial . TangentVectorV1 = surfaceBasisVector1 ; integrationPointMaterial . TangentVectorV2 = surfaceBasisVector2 ; integrationPointMaterial . NormalVectorV3 = surfaceBasisVector3 ; } double norm = 0 ; for ( int i = 0 ; i < surfaceBasisVector3 . Length ; i ++ ) norm += surfaceBasisVector3 [ i ] * surfaceBasisVector3 [ i ] ; var J1 = Math . Sqrt ( norm ) ; for ( int i = 0 ; i < surfaceBasisVector3 . Length ; i ++ ) surfaceBasisVector3 [ i ] = surfaceBasisVector3 [ i ] / J1 ; var surfaceBasisVectorDerivative1 = CalculateSurfaceBasisVector1 ( hessianMatrix , 0 ) ; var surfaceBasisVectorDerivative2 = CalculateSurfaceBasisVector1 ( hessianMatrix , 1 ) ; var surfaceBasisVectorDerivative12 = CalculateSurfaceBasisVector1 ( hessianMatrix , 2 ) ; var Bmembrane = CalculateMembraneDeformationMatrix ( _tsplines , j , surfaceBasisVector1 , surfaceBasisVector2 , elementControlPoints ) ; var Bbending = CalculateBendingDeformationMatrix ( surfaceBasisVector3 , _tsplines , j , surfaceBasisVector2 , surfaceBasisVectorDerivative1 , surfaceBasisVector1 , J1 , surfaceBasisVectorDerivative2 , surfaceBasisVectorDerivative12 , elementControlPoints ) ; var ( MembraneConstitutiveMatrix , BendingConstitutiveMatrix , CouplingConstitutiveMatrix ) = IntegratedConstitutiveOverThickness ( gaussPoints , j ) ; double wFactor = J1 * gaussPoints [ j ] . WeightFactor ; double tempb = 0 ; double tempm = 0 ; Array . Clear ( BmTranspose , 0 , bRows * bCols ) ; Array . Clear ( BbTranspose , 0 , bRows * bCols ) ; for ( int i = 0 ; i < bRows ; i ++ ) { for ( int k = 0 ; k < bCols ; k ++ ) { BmTranspose [ k , i ] = Bmembrane [ i , k ] * wFactor ; BbTranspose [ k , i ] = Bbending [ i , k ] * wFactor ; } } double tempcm = 0 ; double tempcb = 0 ; double tempcc = 0 ; Array . Clear ( BmTransposeMultStiffness , 0 , bRows * bCols ) ; Array . Clear ( BbTransposeMultStiffness , 0 , bRows * bCols ) ; Array . Clear ( BmbTransposeMultStiffness , 0 , bRows * bCols ) ; Array . Clear ( BbmTransposeMultStiffness , 0 , bRows * bCols ) ; for ( int i = 0 ; i < bCols ; i ++ ) { for ( int k = 0 ; k < bRows ; k ++ ) { tempm = BmTranspose [ i , k ] ; tempb = BbTranspose [ i , k ] ; for ( int m = 0 ; m < bRows ; m ++ ) { tempcm = MembraneConstitutiveMatrix [ k , m ] ; tempcb = BendingConstitutiveMatrix [ k , m ] ; tempcc = CouplingConstitutiveMatrix [ k , m ] ; BmTransposeMultStiffness [ i , m ] += tempm * tempcm ; BbTransposeMultStiffness [ i , m ] += tempb * tempcb ; BmbTransposeMultStiffness [ i , m ] += tempm * tempcc ; BbmTransposeMultStiffness [ i , m ] += tempb * tempcc ; } } } double tempmb = 0 ; double tempbm = 0 ; double mem = 0 ; double ben = 0 ; for ( int i = 0 ; i < bCols ; i ++ ) { for ( int k = 0 ; k < bRows ; k ++ ) { tempm = BmTransposeMultStiffness [ i , k ] ; tempb = BbTransposeMultStiffness [ i , k ] ; tempmb = BmbTransposeMultStiffness [ i , k ] ; tempbm = BbmTransposeMultStiffness [ i , k ] ; for ( int m = 0 ; m < bCols ; m ++ ) { mem = Bmembrane [ k , m ] ; ben = Bbending [ k , m ] ; stiffnessMatrix [ i , m ] += tempm * mem + tempb * ben + tempmb * ben + tempbm * mem ; } } } } return Matrix . CreateFromArray ( stiffnessMatrix ) ; }
Calculates the stiffness matrix of the element.
Calculates the stiffness matrix of the element .
An  containing the stiffness matrix of an .
An containing the stiffness matrix of an .
element
IElement
An element of type .
An element of type .