From eeb11011f75814eaac95e0c0768cc4520ca99b57 Mon Sep 17 00:00:00 2001 From: alexios <106801979+alexios725@users.noreply.github.com> Date: Fri, 7 Oct 2022 20:06:50 +0300 Subject: [PATCH] Implemented IEmbeddedHostElement in ContinuumElement3D. Added the following classes in FEM.Structural.Embedding: BeamElementEmbedder, BeamElementEmbedderBase, CohesiveBeam3DToBeam3D, EmbeddedBeam3DGrouping Added the following classes in FEM.Structural.Line: SupportiveCohesiveBeam3DCorotationalAbstract, SupportiveCohesiveBeam3DCorotationalQuaternion Added a unit test in FEM.Structural.Tests: CohesiveBeam3DToBeam3DTest --- .../Continuum/ContinuumElement3D.cs | 142 +++- .../Embedding/BeamElementEmbedder.cs | 26 + .../Embedding/BeamElementEmbedderBase.cs | 344 +++++++++ .../Embedding/CohesiveBeam3DToBeam3D.cs | 686 ++++++++++++++++++ .../Embedding/EmbeddedBeam3DGrouping.cs | 91 +++ ...rtiveCohesiveBeam3DCorotationalAbstract.cs | 636 ++++++++++++++++ ...iveCohesiveBeam3DCorotationalQuaternion.cs | 211 ++++++ .../Elements/CohesiveBeam3DToBeam3DTest.cs | 91 +++ 8 files changed, 2226 insertions(+), 1 deletion(-) create mode 100644 src/MGroup.FEM.Structural/Embedding/BeamElementEmbedder.cs create mode 100644 src/MGroup.FEM.Structural/Embedding/BeamElementEmbedderBase.cs create mode 100644 src/MGroup.FEM.Structural/Embedding/CohesiveBeam3DToBeam3D.cs create mode 100644 src/MGroup.FEM.Structural/Embedding/EmbeddedBeam3DGrouping.cs create mode 100644 src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalAbstract.cs create mode 100644 src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalQuaternion.cs create mode 100644 tests/MGroup.FEM.Structural.Tests/Elements/CohesiveBeam3DToBeam3DTest.cs diff --git a/src/MGroup.FEM.Structural/Continuum/ContinuumElement3D.cs b/src/MGroup.FEM.Structural/Continuum/ContinuumElement3D.cs index f5d2bb1..1dd1c66 100644 --- a/src/MGroup.FEM.Structural/Continuum/ContinuumElement3D.cs +++ b/src/MGroup.FEM.Structural/Continuum/ContinuumElement3D.cs @@ -15,6 +15,8 @@ using MGroup.MSolve.Discretization.Meshes; using MGroup.MSolve.Geometry.Coordinates; using MGroup.LinearAlgebra.Vectors; +using MGroup.MSolve.Discretization.Embedding; +using System.Linq; namespace MGroup.FEM.Structural.Continuum { @@ -23,7 +25,7 @@ namespace MGroup.FEM.Structural.Continuum /// the appropriate , etc. strategies. /// Authors: Dimitris Tsapetis /// - public class ContinuumElement3D : IStructuralElementType, ICell + public class ContinuumElement3D : IStructuralElementType, ICell, IEmbeddedHostElement { private readonly static IDofType[] nodalDOFTypes = new IDofType[] { @@ -396,5 +398,143 @@ private Matrix BuildShapeFunctionMatrix(double[] shapeFunctions) return shapeFunctionMatrix; } + public EmbeddedNode BuildHostElementEmbeddedNode(IElementType element, INode node, + IEmbeddedDOFInHostTransformationVector transformation) + { + var points = GetNaturalCoordinates(element, (Node)node); + if (points.Length == 0) return null; + + //element.EmbeddedNodes.Add(node); + var embeddedNode = new EmbeddedNode(node, element, transformation.GetDependentDOFTypes); + for (int i = 0; i < points.Length; i++) embeddedNode.Coordinates.Add(points[i]); + return embeddedNode; + } + + private double[] GetNaturalCoordinates(IElementType element, Node node) + { + double[] mins = new double[] { element.Nodes[0].X, element.Nodes[0].Y, element.Nodes[0].Z }; + double[] maxes = new double[] { element.Nodes[0].X, element.Nodes[0].Y, element.Nodes[0].Z }; + for (int i = 0; i < element.Nodes.Count; i++) + { + mins[0] = mins[0] > element.Nodes[i].X ? element.Nodes[i].X : mins[0]; + mins[1] = mins[1] > element.Nodes[i].Y ? element.Nodes[i].Y : mins[1]; + mins[2] = mins[2] > element.Nodes[i].Z ? element.Nodes[i].Z : mins[2]; + maxes[0] = maxes[0] < element.Nodes[i].X ? element.Nodes[i].X : maxes[0]; + maxes[1] = maxes[1] < element.Nodes[i].Y ? element.Nodes[i].Y : maxes[1]; + maxes[2] = maxes[2] < element.Nodes[i].Z ? element.Nodes[i].Z : maxes[2]; + } + //return new double[] { (node.X - mins[0]) / ((maxes[0] - mins[0]) / 2) - 1, + // (node.Y - mins[1]) / ((maxes[1] - mins[1]) / 2) - 1, + // (node.Z - mins[2]) / ((maxes[2] - mins[2]) / 2) - 1 }; + + bool maybeInsideElement = node.X <= maxes[0] && node.X >= mins[0] && + node.Y <= maxes[1] && node.Y >= mins[1] && + node.Z <= maxes[2] && node.Z >= mins[2]; + if (maybeInsideElement == false) return new double[0]; + + const int jacobianSize = 3; + const int maxIterations = 1000; + const double tolerance = 1e-10; + int iterations = 0; + double deltaNaturalCoordinatesNormSquare = 100; + double[] naturalCoordinates = new double[] { 0, 0, 0 }; + const double toleranceSquare = tolerance * tolerance; + + while (deltaNaturalCoordinatesNormSquare > toleranceSquare && iterations < maxIterations) + { + iterations++; + var shapeFunctions = Interpolation.EvaluateFunctionsAt(new NaturalPoint(naturalCoordinates)); + double[] coordinateDifferences = new double[] { 0, 0, 0 }; + for (int i = 0; i < shapeFunctions.Length; i++) + { + coordinateDifferences[0] += shapeFunctions[i] * element.Nodes[i].X; + coordinateDifferences[1] += shapeFunctions[i] * element.Nodes[i].Y; + coordinateDifferences[2] += shapeFunctions[i] * element.Nodes[i].Z; + } + coordinateDifferences[0] = node.X - coordinateDifferences[0]; + coordinateDifferences[1] = node.Y - coordinateDifferences[1]; + coordinateDifferences[2] = node.Z - coordinateDifferences[2]; + + double[,] faXYZ = GetCoordinatesTranspose(element); + var nablaShapeFunctions = Interpolation.EvaluateNaturalGradientsAt(new NaturalPoint(naturalCoordinates)); + var inverseJacobian = new IsoparametricJacobian3D(element.Nodes.ToArray(), nablaShapeFunctions).InverseMatrix; + + double[] deltaNaturalCoordinates = new double[] { 0, 0, 0 }; + for (int i = 0; i < jacobianSize; i++) + for (int j = 0; j < jacobianSize; j++) + deltaNaturalCoordinates[i] += inverseJacobian[j, i] * coordinateDifferences[j]; + for (int i = 0; i < 3; i++) + naturalCoordinates[i] += deltaNaturalCoordinates[i]; + + deltaNaturalCoordinatesNormSquare = 0; + for (int i = 0; i < 3; i++) + deltaNaturalCoordinatesNormSquare += deltaNaturalCoordinates[i] * deltaNaturalCoordinates[i]; + //deltaNaturalCoordinatesNormSquare = Math.Sqrt(deltaNaturalCoordinatesNormSquare); + } + + return naturalCoordinates.Count(x => Math.Abs(x) - 1.0 > tolerance) > 0 ? new double[0] : naturalCoordinates; + } + + protected double[,] GetCoordinatesTranspose(IElementType element) + { + double[,] nodeCoordinatesXYZ = new double[3, dofTypes.Length]; + for (int i = 0; i < dofTypes.Length; i++) + { + nodeCoordinatesXYZ[0, i] = element.Nodes[i].X; + nodeCoordinatesXYZ[1, i] = element.Nodes[i].Y; + nodeCoordinatesXYZ[2, i] = element.Nodes[i].Z; + } + return nodeCoordinatesXYZ; + } + + double[] IEmbeddedHostElement.GetShapeFunctionsForNode(IElementType element, EmbeddedNode node) + { + var naturalPoint = new NaturalPoint(node.Coordinates.ToArray()); + var shapeFunctions = Interpolation.EvaluateFunctionsAt(naturalPoint); + var shapeFunctionGradients = Interpolation.EvaluateNaturalGradientsAt(naturalPoint); + var jacobian = new IsoparametricJacobian3D(element.Nodes.ToArray(), shapeFunctionGradients); + + var shapeFunctionGradients2DArray = shapeFunctionGradients.CopyToArray2D(); + var shapeFunctionGradients1DArray = new double[shapeFunctionGradients2DArray.GetLength(0) * shapeFunctionGradients2DArray.GetLength(1)]; + for (int j = 0; j < shapeFunctionGradients2DArray.GetLength(1); j++) + { + for (int i = 0; i < shapeFunctionGradients2DArray.GetLength(0); i++) + { + shapeFunctionGradients1DArray[(shapeFunctionGradients2DArray.GetLength(0) * j) + i] = shapeFunctionGradients2DArray[i, j]; + } + } + + var jacobianDirect2DArray = jacobian.DirectMatrix.CopyToArray2D(); + var jacobianInverse2DArray = jacobian.InverseMatrix.CopyToArray2D(); + var jacobianDirect1DArray = new double[jacobianDirect2DArray.GetLength(0) * jacobianDirect2DArray.GetLength(1)]; + var jacobianInverse1DArray = new double[jacobianInverse2DArray.GetLength(0) * jacobianInverse2DArray.GetLength(1)]; + for (int i = 0; i < jacobianDirect2DArray.GetLength(0); i++) + { + for (int j = 0; j < jacobianDirect2DArray.GetLength(1); j++) + { + jacobianDirect1DArray[(jacobianDirect2DArray.GetLength(1) * i) + j] = jacobianDirect2DArray[i, j]; + jacobianInverse1DArray[(jacobianInverse2DArray.GetLength(1) * i) + j] = jacobianInverse2DArray[i, j]; + } + } + var returnValueList = new List(); + foreach (double shapeFunction in shapeFunctions) + { + returnValueList.Add(shapeFunction); + } + foreach (double shapeFunctionGradient in shapeFunctionGradients1DArray) + { + returnValueList.Add(shapeFunctionGradient); + } + foreach (double value in jacobianDirect1DArray) + { + returnValueList.Add(value); + } + foreach (double value in jacobianInverse1DArray) + { + returnValueList.Add(value); + } + + return returnValueList.ToArray(); + } } } diff --git a/src/MGroup.FEM.Structural/Embedding/BeamElementEmbedder.cs b/src/MGroup.FEM.Structural/Embedding/BeamElementEmbedder.cs new file mode 100644 index 0000000..52702fb --- /dev/null +++ b/src/MGroup.FEM.Structural/Embedding/BeamElementEmbedder.cs @@ -0,0 +1,26 @@ +using System.Collections.Generic; +using System.Linq; + +//using ISAAR.MSolve.FEM.Interfaces; + +using MGroup.MSolve.Discretization; +using MGroup.MSolve.Discretization.Embedding; +using MGroup.MSolve.Discretization.Entities; + +namespace MGroup.FEM.Structural.Embedding +{ + public class BeamElementEmbedder : BeamElementEmbedderBase + { + public BeamElementEmbedder(Model model, IElementType embeddedElement, IEmbeddedDOFInHostTransformationVector transformation) + : base(embeddedElement, transformation) + { + } + + protected void CalculateTransformationMatrix() + { + //var e = (IEmbeddedBeamElement)(embeddedElement.ElementType); + base.CalculateTransformationMatrix(); + //transformationMatrix = e.CalculateRotationMatrix().MultiplyRight(transformationMatrix,true); + } + } +} diff --git a/src/MGroup.FEM.Structural/Embedding/BeamElementEmbedderBase.cs b/src/MGroup.FEM.Structural/Embedding/BeamElementEmbedderBase.cs new file mode 100644 index 0000000..2f51422 --- /dev/null +++ b/src/MGroup.FEM.Structural/Embedding/BeamElementEmbedderBase.cs @@ -0,0 +1,344 @@ +using System.Collections.Generic; +using System.Linq; + +using MGroup.LinearAlgebra.Matrices; +using MGroup.MSolve.DataStructures; +using MGroup.MSolve.Discretization; +using MGroup.MSolve.Discretization.Dofs; +using MGroup.MSolve.Discretization.Embedding; +using MGroup.MSolve.Discretization.Entities; + +namespace MGroup.FEM.Structural.Embedding +{ + public class BeamElementEmbedderBase : IElementDofEnumerator + { + private readonly IElementType embeddedElement; + private readonly IEmbeddedDOFInHostTransformationVector transformation; + private readonly Dictionary superElementMap = new Dictionary(); + private readonly Dictionary> dofToHostMapping = new Dictionary>(); + private Matrix transformationMatrix; + //private bool isElementEmbedded = false; + + public BeamElementEmbedderBase(IElementType embeddedElement, IEmbeddedDOFInHostTransformationVector transformation) + { + this.embeddedElement = embeddedElement; + this.transformation = transformation; + Initialize(); + } + + private void InitializeMappings() + { + var e = (IEmbeddedElement)(embeddedElement); + var nodesList = embeddedElement.Nodes.ToList(); + superElementMap.Clear(); + int index = 0; + foreach (var embeddedNode in e.EmbeddedNodes) + { + int nodeOrderInEmbeddedElement = embeddedElement.Nodes.FindFirstIndex(embeddedNode.Node); //TODO: Fix that. It was embeddedElement.Nodes.IndexOf(embeddedNode.Node) + var currentEmbeddedNodeDOFs = embeddedElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedElement)[nodeOrderInEmbeddedElement]; + //var currentNodeDOFs = currentEmbeddedNodeDOFs.Intersect(embeddedNode.DependentDOFs); + var independentEmbeddedDOFs = currentEmbeddedNodeDOFs.Except(embeddedNode.DependentDOFs); + + // TODO: Optimization to exclude host DOFs that embedded node does not depend on. + for (int i = 0; i < embeddedNode.EmbeddedInElement.Nodes.Count; i++) + { + var currentNodeDOFs = embeddedNode.EmbeddedInElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedNode.EmbeddedInElement)[i]; + foreach (var dof in currentNodeDOFs) + { + var superElementDOF = new SuperElementDof() { DOF = dof, EmbeddedNode = embeddedNode.Node, HostNode = embeddedNode.EmbeddedInElement.Nodes[i], Element = embeddedNode.EmbeddedInElement }; + if (!superElementMap.ContainsKey(superElementDOF)) + { + superElementMap.Add(superElementDOF, index); + index++; + } + } + } + + //var independentEmbeddedDOFs = model.NodalDOFsDictionary[embeddedNode.Node.ID].Select(x => x.Key).Except(embeddedNode.DependentDOFs); + //int nodeOrderInEmbeddedElement = embeddedElement.Nodes.IndexOf(embeddedNode.Node); + + //var independentEmbeddedDOFs = embeddedElement.ElementType.DOFEnumerator.GetDOFTypes(embeddedElement)[nodeOrderInEmbeddedElement].Except(embeddedNode.DependentDOFs); + + foreach (var dof in independentEmbeddedDOFs) + { + var superElementDOF = new SuperElementDof() { DOF = dof, EmbeddedNode = embeddedNode.Node, HostNode = null, Element = null }; + if (!superElementMap.ContainsKey(superElementDOF)) + { + superElementMap.Add(superElementDOF, index); + index++; + } + } + } + + foreach (var node in embeddedElement.Nodes.Except(e.EmbeddedNodes.Select(x => x.Node))) + { + int nodeOrderInEmbeddedElement = embeddedElement.Nodes.FindFirstIndex(node); //TODO: Fix that. It was embeddedElement.Nodes.IndexOf(embeddedNode.Node) + var currentNodeDOFs = embeddedElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedElement)[nodeOrderInEmbeddedElement]; + foreach (var dof in currentNodeDOFs) + { + var superElementDOF = new SuperElementDof() { DOF = dof, EmbeddedNode = node, HostNode = null, Element = null }; + if (!superElementMap.ContainsKey(superElementDOF)) + { + superElementMap.Add(superElementDOF, index); + index++; + } + } + } + } + + public void CalculateTransformationMatrix() + { + var e = (IEmbeddedElement)(embeddedElement); + var nodesList = embeddedElement.Nodes.ToList(); + int row = 0; + int col = 0; + int totalRows = embeddedElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedElement).SelectMany(x => x).Count(); + var matrix = new double[totalRows, superElementMap.Count]; + + foreach (var embeddedNode in e.EmbeddedNodes) + { + var localTransformationMatrix = transformation.GetTransformationVector(embeddedNode); + var localHostDOFs = transformation.GetDOFTypesOfHost(embeddedNode); + int nodeOrderInEmbeddedElement = nodesList.FindIndex(x => x.ID == embeddedNode.Node.ID); //TODO: Fix that. It was embeddedElement.Nodes.IndexOf(embeddedNode.Node) + var embeddedNodeDOFQuantity = embeddedElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedElement)[nodeOrderInEmbeddedElement].Count; + int dependentDOFs = transformation.GetDependentDOFTypes.Count; + + for (int i = 0; i < dependentDOFs; i++) + { + col = 0; + for (int j = 0; j < localHostDOFs.Count; j++) + { + for (int k = 0; k < localHostDOFs[j].Count; k++) + { + var superelement = new SuperElementDof() { DOF = localHostDOFs[j][k], Element = embeddedNode.EmbeddedInElement, EmbeddedNode = embeddedNode.Node, HostNode = embeddedNode.EmbeddedInElement.Nodes[j] }; + matrix[row + i, superElementMap[superelement]] = localTransformationMatrix[i][col]; + col++; + } + } + } + row += dependentDOFs; + + var independentEmbeddedDOFs = embeddedElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedElement)[nodeOrderInEmbeddedElement].Except(embeddedNode.DependentDOFs).ToArray(); + for (int j = 0; j < independentEmbeddedDOFs.Length; j++) + { + var superelement = new SuperElementDof() { DOF = independentEmbeddedDOFs[j], Element = null, HostNode = null, EmbeddedNode = embeddedNode.Node }; + matrix[row, superElementMap[superelement]] = 1; + row++; + } + } + + foreach (var node in embeddedElement.Nodes.Except(e.EmbeddedNodes.Select(x => x.Node))) + { + int nodeOrderInEmbeddedElement = embeddedElement.Nodes.FindFirstIndex(node); + //int nodeOrderInEmbeddedElement = node.ID; //TODO: Fix that. It was embeddedElement.Nodes.IndexOf(embeddedNode.Node) + var currentNodeDOFs = embeddedElement.DofEnumerator.GetDofTypesForMatrixAssembly(embeddedElement)[nodeOrderInEmbeddedElement]; + for (int j = 0; j < currentNodeDOFs.Count; j++) + { + var superelement = new SuperElementDof() { DOF = currentNodeDOFs[j], Element = null, HostNode = null, EmbeddedNode = node }; + matrix[row, superElementMap[superelement]] = 1; + row++; + } + } + + //StreamWriter sw = File.CreateText(String.Format("TransformationMatrix{0}.txt", embeddedElement.ID)); + //for (int i = 0; i < totalRows; i++) + //{ + // var line = String.Empty; + // for (int j = 0; j < superElementMap.Count; j++) + // line += matrix[i,j].ToString() + ";"; + // sw.WriteLine(line); + //} + //sw.Close(); + transformationMatrix = Matrix.CreateFromArray(matrix); + } + + private void Initialize() + { + var e = embeddedElement as IEmbeddedElement; + if (e == null) return; + if (e.EmbeddedNodes.Count == 0) return; + + InitializeMappings(); + CalculateTransformationMatrix(); + } + + public IMatrix GetTransformedMatrix(IMatrix matrix) + { + var e = embeddedElement as IEmbeddedElement; + //if (e == null || !isElementEmbedded) return matrix; + if (e == null) return matrix; + if (e.EmbeddedNodes.Count == 0) return matrix; + + return transformationMatrix.ThisTransposeTimesOtherTimesThis(matrix); + } + + public double[] GetTransformedDisplacementsVector(double[] vector) + { + var e = embeddedElement as IEmbeddedElement; + //if (e == null || !isElementEmbedded) return matrix; + if (e == null) return vector; + if (e.EmbeddedNodes.Count == 0) return vector; + + return transformationMatrix.Multiply(vector, false); + } + + public double[] GetTransformedForcesVector(double[] vector) //compa prosthiki msolve + { + var e = embeddedElement as IEmbeddedElement; + //if (e == null || !isElementEmbedded) return matrix; + if (e == null) return vector; + if (e.EmbeddedNodes.Count == 0) return vector; + + return transformationMatrix.Multiply(vector, true); + } + + + //public IReadOnlyList> GetDofTypesForMatrixAssembly(IElementType element) + //{ + // //return element.ElementType.GetElementDOFTypes(element); + + // var dofs = new List>(); + // INode currentNode = null; + // List nodeDOFs = null; + + // foreach (var superElement in superElementMap) + // { + // INode node = superElement.Key.HostNode == null ? superElement.Key.EmbeddedNode : superElement.Key.HostNode; + + // if (currentNode != node) + // { + // if (nodeDOFs != null) + // dofs.Add(nodeDOFs); + // currentNode = node; + // nodeDOFs = new List(); + // } + // nodeDOFs.Add(superElement.Key.DOF); + // } + // if (nodeDOFs != null) + // dofs.Add(nodeDOFs); + + // return dofs; + //} + + public IReadOnlyList> GetDofTypesForMatrixAssembly(IElementType element) + { + //return element.ElementType.GetElementDOFTypes(element); + + var dofs = new List>(); + INode currentNode = null; + INode currentEmbeddedNode = null; + List nodeDOFs = null; + + foreach (var superElement in superElementMap) + { + INode node = superElement.Key.HostNode == null ? superElement.Key.EmbeddedNode : superElement.Key.HostNode; + INode embeddedNode = superElement.Key.EmbeddedNode; + if (currentNode != node || currentEmbeddedNode != embeddedNode) + { + if (nodeDOFs != null) + dofs.Add(nodeDOFs); + currentNode = node; + currentEmbeddedNode = embeddedNode; + nodeDOFs = new List(); + } + nodeDOFs.Add(superElement.Key.DOF); + } + if (nodeDOFs != null) + dofs.Add(nodeDOFs); + + return dofs; + } + + public IReadOnlyList> GetDofTypesForDofEnumeration(IElementType element) + { + //if (embeddedElement != element) throw new ArgumentException(); + + var nodesDictionary = new Dictionary(); + int index = 0; + foreach (var node in element.Nodes) + { + nodesDictionary.Add(node, index); + index++; + } + + var dofs = new List>(); + for (int i = 0; i < element.Nodes.Count; i++) + dofs.Add(new List()); + + INode currentNode = null; + List nodeDOFs = null; + + foreach (var superElement in superElementMap) + { + if (superElement.Key.HostNode != null) continue; + INode node = superElement.Key.EmbeddedNode; + //Node node = superElement.Key.HostNode == null ? superElement.Key.EmbeddedNode : superElement.Key.HostNode; + + if (currentNode != node) + { + if (nodeDOFs != null) + dofs[nodesDictionary[currentNode]] = nodeDOFs; + currentNode = node; + nodeDOFs = new List(); + } + nodeDOFs.Add(superElement.Key.DOF); + } + if (nodeDOFs != null) + dofs[nodesDictionary[currentNode]] = nodeDOFs; + //dofs.Add(nodeDOFs); + + return dofs; + } + + //public IReadOnlyList GetNodesForMatrixAssembly(IElementType element) + //{ + // var nodes = new List(); + // INode currentNode = null; + // foreach (var superElement in superElementMap) + // { + // INode node = superElement.Key.HostNode == null ? superElement.Key.EmbeddedNode : superElement.Key.HostNode; + // if (currentNode != node) + // { + // if (currentNode != null) + // nodes.Add(currentNode); + // currentNode = node; + // } + // //if (nodes.IndexOf(node) < 0) + // // nodes.Add(node); + // } + // if (currentNode != null) + // nodes.Add(currentNode); + + // return nodes; + //} + + public IReadOnlyList GetNodesForMatrixAssembly(IElementType element) + { + var nodes = new List(); + INode currentNode = null; + INode currentEmbeddedNode = null; + foreach (var superElement in superElementMap) + { + INode node = superElement.Key.HostNode == null ? superElement.Key.EmbeddedNode : superElement.Key.HostNode; + INode embeddedNode = superElement.Key.EmbeddedNode; + if (currentNode != node || currentEmbeddedNode != embeddedNode) + { + if (currentNode != null) + { + nodes.Add(currentNode); + } + currentNode = node; + currentEmbeddedNode = embeddedNode; + } + //if (nodes.IndexOf(node) < 0) + // nodes.Add(node); + } + if (currentNode != null) + nodes.Add(currentNode); + + return nodes; + } + + + } +} diff --git a/src/MGroup.FEM.Structural/Embedding/CohesiveBeam3DToBeam3D.cs b/src/MGroup.FEM.Structural/Embedding/CohesiveBeam3DToBeam3D.cs new file mode 100644 index 0000000..39720cc --- /dev/null +++ b/src/MGroup.FEM.Structural/Embedding/CohesiveBeam3DToBeam3D.cs @@ -0,0 +1,686 @@ +using System; +using System.Collections.Generic; +using System.Xml.Linq; + +using MGroup.Constitutive.Structural; +using MGroup.Constitutive.Structural.Cohesive; +using MGroup.Constitutive.Structural.Continuum; +using MGroup.Constitutive.Structural.Line; +using MGroup.FEM.Structural.Line; +using MGroup.LinearAlgebra.Matrices; +using MGroup.LinearAlgebra.Vectors; +using MGroup.MSolve.Discretization; +using MGroup.MSolve.Discretization.Dofs; +using MGroup.MSolve.Discretization.Embedding; +using MGroup.MSolve.Discretization.Entities; +using MGroup.MSolve.Numerics.Integration.Quadratures; +using MGroup.MSolve.Numerics.Interpolation; + +namespace MGroup.FEM.Structural.Embedding +{ + public class CohesiveBeam3DToBeam3D : IStructuralElementType, IEmbeddedElement + { + private List nodes = new List(); + protected readonly static IDofType[] nodalDOFTypes = new IDofType[] { StructuralDof.TranslationX, StructuralDof.TranslationY, StructuralDof.TranslationZ, StructuralDof.RotationX, StructuralDof.RotationY, StructuralDof.RotationZ }; + protected readonly static IDofType[][] dofTypes = new IDofType[][] { nodalDOFTypes, nodalDOFTypes, nodalDOFTypes, + nodalDOFTypes }; + protected readonly ICohesiveZoneMaterial[] materialsAtGaussPoints; + private readonly InterpolationTruss1D interpolation = InterpolationTruss1D.UniqueInstance; + private int nGaussPoints; + protected IElementDofEnumerator dofEnumerator = new GenericDofEnumerator(); + private SupportiveCohesiveBeam3DCorotationalQuaternion supportive_beam; + private SupportiveCohesiveBeam3DCorotationalQuaternion supportive_clone; + private readonly double perimeter; + + /// + /// Initial nodel coordinates of 4 node inner cohesive element + /// + private double[][] initialNodalCoordinates; // ox_i + + /// + /// Unrolled current nodal coordinates of 4 node inner cohesive element + /// + private double[] currentNodalCoordinates; // x_local + + public bool MatrixIsNotInitialized = true; + public double[,] k_cohesive_element_total {get; private set;} + private double[][] tractions; + private double[] localDisplacements; + private double[] localDisplacementsLastConverged; + //protected CohesiveBeam3DToBeam3D() + //{ + //} + + public CohesiveBeam3DToBeam3D(List nodes, ICohesiveZoneMaterial material, IQuadrature1D quadratureForStiffness, IList nodes_beam, + IList nodes_clone, IIsotropicContinuumMaterial3D material_beam, double density, BeamSection3D beamSection, double perimeter) + { + this.QuadratureForStiffness = quadratureForStiffness; + this.nGaussPoints = quadratureForStiffness.IntegrationPoints.Count; + materialsAtGaussPoints = new ICohesiveZoneMaterial[nGaussPoints]; + for (int i = 0; i < nGaussPoints; i++) + { + var materialCopy = (ICloneable)material; + materialsAtGaussPoints[i] = (ICohesiveZoneMaterial)materialCopy.Clone(); + } + this.supportive_beam = new SupportiveCohesiveBeam3DCorotationalQuaternion(nodes_beam, material_beam, density, beamSection); + this.supportive_clone = new SupportiveCohesiveBeam3DCorotationalQuaternion(nodes_clone, material_beam, density, beamSection); + this.perimeter = perimeter; + this.nodes = nodes; + } + + public IQuadrature1D QuadratureForStiffness { get; } + + private void GetInitialGeometricDataAndInitializeMatrices(IElementType element) + { + initialNodalCoordinates = new double[4][]; + + //TODOgsoim: why 0 to 4?? + for (int j = 0; j < 4; j++) + { initialNodalCoordinates[j] = new double[] { element.Nodes[j].X, element.Nodes[j].Y, element.Nodes[j].Z, }; } + + currentNodalCoordinates = new double[12]; //12:without rotations, 24:with rotations + } + + private double[][] UpdateCoordinateDataAndCalculateDisplacementVector(double[] localdisplacements) + { + //Matrix shapeFunctionDerivatives = interpolation.EvaluateGradientsAt(); + IReadOnlyList N1 = interpolation.EvaluateFunctionsAtGaussPoints(QuadratureForStiffness); + //IReadOnlyList N3 = interpolation.EvaluateN3ShapeFunctionsReorganized(QuadratureForStiffness); //Shape functions matrix [N_beam] + + double[,] u_prok = new double[3, 2];// [3d-axes, nodes] - of the mid-surface(#i_m, #j_m) + double[,] x_bar = new double[3, 2]; + double[] u = new double[3]; + + double[][] Delta = new double[nGaussPoints][]; + for (int j = 0; j < nGaussPoints; j++) + { + Delta[j] = new double[3]; + } + + //double[][,] R = new double[nGaussPoints][,]; //TODO: maybe cache R + //for (int j = 0; j < nGaussPoints; j++) + //{ + // R[j] = new double[3, 3]; + // R[j][0, 0] = 0.5 * (supportive_beam.currentRotationMatrix[0, 0] + supportive_clone.currentRotationMatrix[0, 0]); + // R[j][0, 1] = 0.5 * (supportive_beam.currentRotationMatrix[0, 1] + supportive_clone.currentRotationMatrix[0, 1]); + // R[j][0, 2] = 0.5 * (supportive_beam.currentRotationMatrix[0, 2] + supportive_clone.currentRotationMatrix[0, 2]); + // R[j][1, 0] = 0.5 * (supportive_beam.currentRotationMatrix[1, 0] + supportive_clone.currentRotationMatrix[1, 0]); + // R[j][1, 1] = 0.5 * (supportive_beam.currentRotationMatrix[1, 1] + supportive_clone.currentRotationMatrix[1, 1]); + // R[j][1, 2] = 0.5 * (supportive_beam.currentRotationMatrix[1, 2] + supportive_clone.currentRotationMatrix[1, 2]); + // R[j][2, 0] = 0.5 * (supportive_beam.currentRotationMatrix[2, 0] + supportive_clone.currentRotationMatrix[2, 0]); + // R[j][2, 1] = 0.5 * (supportive_beam.currentRotationMatrix[2, 1] + supportive_clone.currentRotationMatrix[2, 1]); + // R[j][2, 2] = 0.5 * (supportive_beam.currentRotationMatrix[2, 2] + supportive_clone.currentRotationMatrix[2, 2]); + //} + + Matrix R = CalculateRotationMatrix(); + + // Update x_local + for (int j = 0; j < 4; j++) + { + for (int k = 0; k < 3; k++) + { + currentNodalCoordinates[3 * j + k] = initialNodalCoordinates[j][k] + localdisplacements[6 * j + k]; + } + } + for (int j = 0; j < 2; j++) + { + for (int k = 0; k < 3; k++) + { + u_prok[k, j] = currentNodalCoordinates[k + 3 * j] - currentNodalCoordinates[6 + k + 3 * j]; + x_bar[k, j] = currentNodalCoordinates[k + 3 * j] + currentNodalCoordinates[6 + k + 3 * j]; + } + } + + //Calculate Delta for all GPs + for (int npoint1 = 0; npoint1 < nGaussPoints; npoint1++) + { + for (int l = 0; l < 3; l++) + { u[l] = 0; } + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 2; m++) // must be 2 for beam cohesive 8-nodes + { + u[l] += u_prok[l, m] * N1[npoint1][m]; + } + } + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 3; m++) + { + Delta[npoint1][l] += R[m, l] * u[m]; + } + } + } + return Delta; + } + + private Tuple CalculateNecessaryMatricesForStiffnessMatrixAndForcesVectorCalculations() + { + //IReadOnlyList N1 = interpolation.EvaluateFunctionsAtGaussPoints(QuadratureForStiffness); + IReadOnlyList N3 = EvaluateN3ShapeFunctionsReorganized(QuadratureForStiffness); //Shape functions matrix [N_beam] + //Matrix shapeFunctionDerivatives = interpolation.EvaluateGradientsAt(); + + double[] integrationsCoeffs = new double[nGaussPoints]; + Matrix[] RtN3 = new Matrix[nGaussPoints]; + //Matrix[] R = new Matrix[nGaussPoints]; //TODO: perhaps cache matrices in InitializeMatrices() where RtN3 is calculated + //for (int j = 0; j < nGaussPoints; j++) + //{ R[j] = Matrix.CreateZero(3, 3); } + + // Calculate Delta for all GPs + Matrix R = CalculateRotationMatrix(); + for (int npoint1 = 0; npoint1 < nGaussPoints; npoint1++) + { + //double[,] u_prok = new double[3, 2];// [3d-axes, nodes] - of the mid-surface(#i_m, #j_m) + //double[] u = new double[3]; + + double[][] Delta = new double[nGaussPoints][]; + for (int j = 0; j < nGaussPoints; j++) + { + Delta[j] = new double[3]; + } + + //R[npoint1][0, 0] = 0.5 * (supportive_beam.currentRotationMatrix[0, 0] + supportive_clone.currentRotationMatrix[0, 0]); + //R[npoint1][0, 1] = 0.5 * (supportive_beam.currentRotationMatrix[0, 1] + supportive_clone.currentRotationMatrix[0, 1]); + //R[npoint1][0, 2] = 0.5 * (supportive_beam.currentRotationMatrix[0, 2] + supportive_clone.currentRotationMatrix[0, 2]); + //R[npoint1][1, 0] = 0.5 * (supportive_beam.currentRotationMatrix[1, 0] + supportive_clone.currentRotationMatrix[1, 0]); + //R[npoint1][1, 1] = 0.5 * (supportive_beam.currentRotationMatrix[1, 1] + supportive_clone.currentRotationMatrix[1, 1]); + //R[npoint1][1, 2] = 0.5 * (supportive_beam.currentRotationMatrix[1, 2] + supportive_clone.currentRotationMatrix[1, 2]); + //R[npoint1][2, 0] = 0.5 * (supportive_beam.currentRotationMatrix[2, 0] + supportive_clone.currentRotationMatrix[2, 0]); + //R[npoint1][2, 1] = 0.5 * (supportive_beam.currentRotationMatrix[2, 1] + supportive_clone.currentRotationMatrix[2, 1]); + //R[npoint1][2, 2] = 0.5 * (supportive_beam.currentRotationMatrix[2, 2] + supportive_clone.currentRotationMatrix[2, 2]); + + integrationsCoeffs[npoint1] = ((supportive_beam.currentLength + supportive_clone.currentLength) / 4.0) * QuadratureForStiffness.IntegrationPoints[npoint1].Weight; + + + // Calculate RtN3 here instead of in InitializeRN3() and then in UpdateForces() + RtN3[npoint1] = R.Transpose() * N3[npoint1]; + } + return new Tuple(RtN3, integrationsCoeffs); + } + + private double[] UpdateForces(/*IElementType element, */Matrix[] RtN3, double[] integrationCoeffs, double[] localTotalDisplacements) + { + double[] fxk1_coh = new double[24]; // Beam3D: 4 nodes, 6 dofs/node (translations + rotations) + + for (int npoint1 = 0; npoint1 < nGaussPoints; npoint1++) + { + double[] T_int_integration_coeffs = new double[3]; + for (int l = 0; l < 3; l++) + { + T_int_integration_coeffs[l] = /*materialsAtGaussPoints[npoint1].Tractions[l]*/Tractions[npoint1][l] * integrationCoeffs[npoint1] * perimeter; + } + double[] r_int_1 = new double[6]; + for (int l = 0; l < 6; l++) + { + for (int m = 0; m < 3; m++) + { r_int_1[l] += RtN3[npoint1][m, l] * T_int_integration_coeffs[m]; } + } + for (int l = 0; l < 3; l++) + { + fxk1_coh[l] += r_int_1[l]; + fxk1_coh[12 + l] += (-r_int_1[l]); + } + for (int l = 6; l < 9; l++) + { + fxk1_coh[l] += r_int_1[l - 3]; + fxk1_coh[12 + l] += (-r_int_1[l - 3]); + } + + Matrix R = CalculateRotationMatrix(); + Matrix ConstaintsRotations = Matrix.CreateZero(3, 3); + ConstaintsRotations[0, 0] = 1000.0; + ConstaintsRotations[1, 1] = 1000.0; + ConstaintsRotations[2, 2] = 1000.0; + Matrix Constr_R = ConstaintsRotations * R; + Matrix M2 = R.Transpose() * Constr_R; + var r_int_2a = M2 * Vector.CreateFromArray(new double[3] {localTotalDisplacements[3]-localTotalDisplacements[15], + localTotalDisplacements[4] - localTotalDisplacements[16], localTotalDisplacements[5]-localTotalDisplacements[17] }) * integrationCoeffs[npoint1]; + var r_int_2b = M2 * Vector.CreateFromArray(new double[3] {localTotalDisplacements[9]-localTotalDisplacements[21], + localTotalDisplacements[10] - localTotalDisplacements[22], localTotalDisplacements[11]-localTotalDisplacements[23] }) * integrationCoeffs[npoint1]; + + for (int i = 0; i < 3; i++) + { + fxk1_coh[i + 3] += r_int_2a[i]; + fxk1_coh[i + 9] += r_int_2b[i]; + + fxk1_coh[i + 15] += -r_int_2a[i]; + fxk1_coh[i + 21] += -r_int_2b[i]; + } + } + + //Matrix R = CalculateRotationMatrix(); + //Matrix ConstaintsRotations = Matrix.CreateZero(3, 3); + //ConstaintsRotations[0, 0] = 1.0; // 10.0; // + //Matrix Constr_R = ConstaintsRotations * R; + //Matrix M2 = R.Transpose() * Constr_R; + //var r_int_2a = M2 * Vector.CreateFromArray(new double[3] {localTotalDisplacements[3]-localTotalDisplacements[15], + // localTotalDisplacements[4] - localTotalDisplacements[16], localTotalDisplacements[5]-localTotalDisplacements[17] }); + //var r_int_2b = M2 * Vector.CreateFromArray(new double[3] {localTotalDisplacements[9]-localTotalDisplacements[21], + // localTotalDisplacements[10] - localTotalDisplacements[22], localTotalDisplacements[11]-localTotalDisplacements[23] }); + + //for (int i = 0; i < 3; i++) + //{ + // fxk1_coh[i + 3] = r_int_2a[i]; + // fxk1_coh[i + 9] = r_int_2b[i]; + + // fxk1_coh[i + 15] = -r_int_2a[i]; + // fxk1_coh[i + 21] = -r_int_2b[i]; + //} + + + return fxk1_coh; + } + + private double[,] UpdateKmatrices(IElementType element, Matrix[] RtN3, double[] integrationCoeffs) + { + k_cohesive_element_total = new double[24, 24]; + //double[,] k_cohesive_element = new double[12, 12]; + + for (int npoint1 = 0; npoint1 < nGaussPoints; npoint1++) + { + Matrix D_tan_sunt_ol = Matrix.CreateZero(3, 3); + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 3; m++) + { + D_tan_sunt_ol[l, m] = materialsAtGaussPoints[npoint1].ConstitutiveMatrix[l, m] * integrationCoeffs[npoint1]; + } + } + Matrix D_RtN3_sunt_ol = D_tan_sunt_ol * RtN3[npoint1]; + Matrix M = RtN3[npoint1].Transpose() * D_RtN3_sunt_ol * perimeter; //TODOgsoim: is perimeter wright here?? + + //k_cohesive_element_total - Contrains rotations, that create zero-elements in the diagonal of the total stiffness matrix. + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 3; m++) + { + k_cohesive_element_total[l, m] += M[l, m]; + k_cohesive_element_total[l, 12 + m] += -M[l, m]; + k_cohesive_element_total[12 + l, m] += -M[l, m]; + k_cohesive_element_total[12 + l, 12 + m] += M[l, m]; + } + } + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 3; m++) + { + k_cohesive_element_total[l + 6, m] += M[l + 3, m]; + k_cohesive_element_total[l + 6, 12 + m] += -M[l + 3, m]; + k_cohesive_element_total[l + 18, m] += -M[l + 3, m]; + k_cohesive_element_total[l + 18, 12 + m] += M[l + 3, m]; + } + } + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 3; m++) + { + k_cohesive_element_total[l, m + 6] += M[l, m + 3]; + k_cohesive_element_total[l, m + 18] += -M[l, m + 3]; + k_cohesive_element_total[l + 12, m + 6] += -M[l, m + 3]; + k_cohesive_element_total[l + 12, m + 18] += M[l, m + 3]; + } + } + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 3; m++) + { + k_cohesive_element_total[l + 6, m + 6] += M[l + 3, m + 3]; + k_cohesive_element_total[l + 6, m + 18] += -M[l + 3, m + 3]; + k_cohesive_element_total[l + 18, m + 6] += -M[l + 3, m + 3]; + k_cohesive_element_total[l + 18, m + 18] += M[l + 3, m + 3]; + } + } + + // ***NEW*** + Matrix R = CalculateRotationMatrix(); + Matrix ConstaintsRotations = Matrix.CreateZero(3, 3); + ConstaintsRotations[0, 0] = 1000.0; + ConstaintsRotations[1, 1] = 1000.0; + ConstaintsRotations[2, 2] = 1000.0; + Matrix Constr_R = ConstaintsRotations * R; + Matrix M2 = R.Transpose() * Constr_R * integrationCoeffs[npoint1] * perimeter; //TODOgsoim: should perimeter be here?? + + for (int ii = 0; ii < 3; ii++) + { + for (int jj = 0; jj < 3; jj++) + { + k_cohesive_element_total[ii + 3, jj + 3] += M2[ii, jj]; // + k_cohesive_element_total[ii + 3, jj + 15] += -M2[ii, jj]; + k_cohesive_element_total[ii + 9, jj + 9] += M2[ii, jj]; // + k_cohesive_element_total[ii + 9, jj + 21] += -M2[ii, jj]; + k_cohesive_element_total[ii + 15, jj + 3] += -M2[ii, jj]; + k_cohesive_element_total[ii + 15, jj + 15] += M2[ii, jj]; // + k_cohesive_element_total[ii + 21, jj + 9] += -M2[ii, jj]; + k_cohesive_element_total[ii + 21, jj + 21] += M2[ii, jj]; // + } + } + } + return k_cohesive_element_total; //k_cohesive_element; // + } + + public Tuple CalculateResponse(double[] localTotalDisplacementsSuperElement) + { + if (this.localDisplacements == null) + { + this.localDisplacements = new double[localTotalDisplacementsSuperElement.Length]; + } + if (localDisplacementsLastConverged == null) + { + localDisplacementsLastConverged = new double[localTotalDisplacementsSuperElement.Length]; + } + double[] localdDisplacementsSuperElement = new double[this.localDisplacements.Length]; + + for (int i = 0; i < this.localDisplacements.Length; i++) + { + localdDisplacementsSuperElement[i] = localTotalDisplacementsSuperElement[i] - localDisplacementsLastConverged[i]; + } + this.localDisplacements = localTotalDisplacementsSuperElement; + double[][] Delta = new double[nGaussPoints][]; + double[] localTotalDisplacements = dofEnumerator.GetTransformedDisplacementsVector(localTotalDisplacementsSuperElement); + double[] localTotaldDisplacements = dofEnumerator.GetTransformedDisplacementsVector(localdDisplacementsSuperElement); + double[] localTotaldDisplacements_beam = new double[12]; + double[] localTotaldDisplacements_clone = new double[12]; + + for (int i1 = 0; i1 < 12; i1++) + { + localTotaldDisplacements_beam[i1] = localTotaldDisplacements[12 + i1]; + localTotaldDisplacements_clone[i1] = localTotaldDisplacements[i1]; + } + + supportive_beam.CalculateStresses(localTotaldDisplacements_beam); + supportive_clone.CalculateStresses(localTotaldDisplacements_clone); + Delta = this.UpdateCoordinateDataAndCalculateDisplacementVector(localTotalDisplacements); + + double[][] tractions = new double[materialsAtGaussPoints.Length][]; + for (int i = 0; i < materialsAtGaussPoints.Length; i++) + { + tractions[i] = materialsAtGaussPoints[i].UpdateConstitutiveMatrixAndEvaluateResponse(Delta[i]); + } + this.tractions = tractions; + return new Tuple(Delta[materialsAtGaussPoints.Length - 1], tractions[materialsAtGaussPoints.Length - 1]); + } + + public double[] CalculateResponseIntegral() + { + double[] localTotalDisplacementsSuperElement = this.localDisplacements; + double[] localTotalDisplacements = dofEnumerator.GetTransformedDisplacementsVector(localTotalDisplacementsSuperElement); + double[] fxk2_coh; + Tuple RtN3AndIntegrationCoeffs; + RtN3AndIntegrationCoeffs = CalculateNecessaryMatricesForStiffnessMatrixAndForcesVectorCalculations(); // Rt * Nbeam + Matrix[] RtN3; + RtN3 = RtN3AndIntegrationCoeffs.Item1; + double[] integrationCoeffs; + integrationCoeffs = RtN3AndIntegrationCoeffs.Item2; + fxk2_coh = this.UpdateForces(/*element, */RtN3, integrationCoeffs, localTotalDisplacements); // sxesh 18 k 19 + return dofEnumerator.GetTransformedForcesVector(fxk2_coh);// embedding + } + + public double[] CalculateResponseIntegralForLogging(double[] localDisplacements) + { + return CalculateResponseIntegral(); + } + + public IMatrix StiffnessMatrix(IElementType element) + { + double[,] k_stoixeiou_coh2; + if (MatrixIsNotInitialized) + { + this.GetInitialGeometricDataAndInitializeMatrices(element); + this.UpdateCoordinateDataAndCalculateDisplacementVector(new double[24]); //returns Delta that can't be used for the initial material state + MatrixIsNotInitialized = false; + } + + Tuple RtN3AndIntegrationCoeffs; + RtN3AndIntegrationCoeffs = CalculateNecessaryMatricesForStiffnessMatrixAndForcesVectorCalculations(); //Rt*Nbeam + Matrix[] RtN3; + RtN3 = RtN3AndIntegrationCoeffs.Item1; + double[] integrationCoeffs; + integrationCoeffs = RtN3AndIntegrationCoeffs.Item2; + k_stoixeiou_coh2 = this.UpdateKmatrices(element, RtN3, integrationCoeffs); + IMatrix element_stiffnessMatrix = Matrix.CreateFromArray(k_stoixeiou_coh2); + return dofEnumerator.GetTransformedMatrix(element_stiffnessMatrix); + } + + //public bool MaterialModified + //{ + // get + // { + // foreach (ICohesiveZoneMaterial3D material in materialsAtGaussPoints) + // if (material.Modified) return true; + // return false; + // } + //} + + //public void ResetMaterialModified() + //{ + // foreach (ICohesiveZoneMaterial3D material in materialsAtGaussPoints) material.ResetModified(); + //} + + //public void ClearMaterialState() + //{ + // foreach (ICohesiveZoneMaterial3D m in materialsAtGaussPoints) m.ClearState(); + //} + + //public void SaveMaterialState() + //{ + // foreach (ICohesiveZoneMaterial3D m in materialsAtGaussPoints) m.SaveState(); + // supportive_beam.SaveMaterialState(); + // supportive_clone.SaveMaterialState(); + //} + + //public void ClearMaterialStresses() + //{ + // foreach (ICohesiveZoneMaterial3D m in materialsAtGaussPoints) m.ClearTractions(); + //} + + //public int ID + //{ + // get { return 13; } + //} + public ElementDimensions ElementDimensions + { + get { return ElementDimensions.ThreeD; } + } + + public IElementDofEnumerator DOFEnumerator + { + get { return dofEnumerator; } + set { dofEnumerator = value; } + } + + //public virtual IList> GetElementDOFTypes(IElementType element) + //{ + // return dofTypes; + //} + + //public double[] CalculateAccelerationForces(IElementType element, IList loads) + //{ + // return new double[64]; + //} + + public virtual IMatrix MassMatrix(IElementType element) + { + return Matrix.CreateZero(64, 64); + } + + public virtual IMatrix DampingMatrix(IElementType element) + { + + return Matrix.CreateZero(64, 64); + } + + #region EMBEDDED + private readonly IReadOnlyList embeddedNodes = new List(); + public IList EmbeddedNodes { get { return (IList)embeddedNodes; } } + + public IElementDofEnumerator DofEnumerator + { + get { return dofEnumerator; } + set { dofEnumerator = value; } + } + + public Dictionary GetInternalNodalDOFs(IElementType element, INode node)// + { + int index = 0; + foreach (var elementNode in element.Nodes) + { + if (node.ID == elementNode.ID) + break; + index++; + } + if (index >= 16) + throw new ArgumentException(String.Format("GetInternalNodalDOFs: Node {0} not found in element {1}.", node.ID, element.ID)); + + if (index >= 8) + { + int index2 = index - 8; + return new Dictionary() { { StructuralDof.TranslationX, 39 + 3 * index2 + 1 }, { StructuralDof.TranslationY, 39 + 3 * index2 + 2 }, { StructuralDof.TranslationZ, 39 + 3 * index2 + 3 } }; + } + else + { + return new Dictionary() { { StructuralDof.TranslationX, + 5 * index + 0 }, { StructuralDof.TranslationY, + 5 * index + 1 }, { StructuralDof.TranslationZ, + 5 * index + 2 }, + { StructuralDof.RotationX, + 5 * index + 3 }, { StructuralDof.RotationY, + 5 * index + 4 }}; + } + } + + public double[] GetLocalDOFValues(IElementType hostElement, double[] hostDOFValues) // omoiws Beam3D + { + //if (transformation == null) + // throw new InvalidOperationException("Requested embedded node values for element that has no embedded nodes."); + //if (hostElementList == null) + // throw new InvalidOperationException("Requested host element list for element that has no embedded nodes."); + //int index = hostElementList.IndexOf(hostElement); + //if (index < 0) + // throw new ArgumentException("Requested host element is not inside host element list."); + + //double[] values = new double[transformation.Columns]; + //int multiplier = hostElement.ElementType.DOFEnumerator.GetDOFTypes(hostElement).SelectMany(d => d).Count(); + //int vectorIndex = 0; + //for (int i = 0; i < index; i++) + // vectorIndex += isNodeEmbedded[i] ? 3 : multiplier; + //Array.Copy(hostDOFValues, 0, values, vectorIndex, multiplier); + + //return (transformation * new Vector(values)).Data; + return dofEnumerator.GetTransformedDisplacementsVector(hostDOFValues); + } + + //public double[] CalculateAccelerationForces(Element element, IList loads) + //{ + // throw new NotImplementedException(); + //} + + //IMatrix IElementType.StiffnessMatrix(IElement element) + //{ + // throw new NotImplementedException(); + //} + + //public IMatrix MassMatrix(IElement element) + //{ + // throw new NotImplementedException(); + //} + + //public IMatrix DampingMatrix(IElement element) + //{ + // throw new NotImplementedException(); + //} + + //public virtual IReadOnlyList> GetElementDofTypes(IElementType element) => dofTypes; + + //Dictionary IEmbeddedElement.GetInternalNodalDOFs(Element element, Node node) + //{ + // throw new NotImplementedException(); + //} + + //double[] IEmbeddedElement.GetLocalDOFValues(Element hostElement, double[] hostDOFValues) + //{ + // throw new NotImplementedException(); + //} + #endregion + + public Matrix CalculateRotationMatrix() + { + Matrix R = Matrix.CreateZero(3, 3); + // R = 0.5*(supportive_beam.currentRotationMatrix + supportive_clone.currentRotationMatrix) + R[0, 1] = 0.5 * (supportive_beam.currentRotationMatrix[0, 1] + supportive_clone.currentRotationMatrix[0, 1]); + R[0, 2] = 0.5 * (supportive_beam.currentRotationMatrix[0, 2] + supportive_clone.currentRotationMatrix[0, 2]); + R[1, 0] = 0.5 * (supportive_beam.currentRotationMatrix[1, 0] + supportive_clone.currentRotationMatrix[1, 0]); + R[0, 0] = 0.5 * (supportive_beam.currentRotationMatrix[0, 0] + supportive_clone.currentRotationMatrix[0, 0]); + R[1, 1] = 0.5 * (supportive_beam.currentRotationMatrix[1, 1] + supportive_clone.currentRotationMatrix[1, 1]); + R[1, 2] = 0.5 * (supportive_beam.currentRotationMatrix[1, 2] + supportive_clone.currentRotationMatrix[1, 2]); + R[2, 0] = 0.5 * (supportive_beam.currentRotationMatrix[2, 0] + supportive_clone.currentRotationMatrix[2, 0]); + R[2, 1] = 0.5 * (supportive_beam.currentRotationMatrix[2, 1] + supportive_clone.currentRotationMatrix[2, 1]); + R[2, 2] = 0.5 * (supportive_beam.currentRotationMatrix[2, 2] + supportive_clone.currentRotationMatrix[2, 2]); + return R; + } + + public double[][] Tractions => tractions; + public IMatrix StiffnessMatrix() => StiffnessMatrix(this); + public IMatrix MassMatrix() => throw new NotImplementedException(); + public IMatrix DampingMatrix() => throw new NotImplementedException(); + public void SaveConstitutiveLawState() + { + localDisplacementsLastConverged = new double[localDisplacements.Length]; + for (int i = 0; i < localDisplacements.Length; i++) + { + localDisplacementsLastConverged[i] = localDisplacements[i]; + } + foreach (ICohesiveZoneMaterial m in materialsAtGaussPoints) m.CreateState(); + supportive_beam.SaveMaterialState(); + supportive_clone.SaveMaterialState(); + } + //public Tuple CalculateResponse(double[] localDisplacements) => throw new NotImplementedException(); + //public double[] CalculateResponseIntegral() => throw new NotImplementedException(); + //public double[] CalculateResponseIntegralForLogging(double[] localDisplacements) => throw new NotImplementedException(); + public IMatrix PhysicsMatrix() + { + return StiffnessMatrix(); + } + public IReadOnlyList> GetElementDofTypes() => dofTypes; + + //public IReadOnlyList Materials { get; } + + public CellType CellType { get; } + int IElementType.ID { get; set; } + + public IReadOnlyList Nodes{ get => nodes; } + + public int SubdomainID { get; set; } + + public IReadOnlyList EvaluateN3ShapeFunctionsReorganized(IQuadrature1D quadrature) + { + var cachedN3AtGPs = new Dictionary>(); + bool isCached = cachedN3AtGPs.TryGetValue(quadrature, + out IReadOnlyList N3AtGPs); + if (isCached) + { + return N3AtGPs; + } + else + { + IReadOnlyList N1 = this.interpolation.EvaluateFunctionsAtGaussPoints(quadrature); + N3AtGPs = GetN3ShapeFunctionsReorganized(quadrature, N1); + cachedN3AtGPs.Add(quadrature, N3AtGPs); + return N3AtGPs; + } + } + + private IReadOnlyList GetN3ShapeFunctionsReorganized(IQuadrature1D quadrature, IReadOnlyList N1) + { + //TODO reorganize cohesive shell to use only N1 (not reorganised) + + int nGaussPoints = quadrature.IntegrationPoints.Count; + var N3 = new Matrix[nGaussPoints]; // shapeFunctionsgpData + for (int npoint = 0; npoint < nGaussPoints; npoint++) + { + double ksi = quadrature.IntegrationPoints[npoint].Xi; + var N3gp = Matrix.CreateZero(3, 6); //8=nShapeFunctions; + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < 2; m++) N3gp[l, l + 3 * m] = N1[npoint][m]; + } + N3[npoint] = N3gp; + } + return N3; + } + } +} diff --git a/src/MGroup.FEM.Structural/Embedding/EmbeddedBeam3DGrouping.cs b/src/MGroup.FEM.Structural/Embedding/EmbeddedBeam3DGrouping.cs new file mode 100644 index 0000000..b3ae29b --- /dev/null +++ b/src/MGroup.FEM.Structural/Embedding/EmbeddedBeam3DGrouping.cs @@ -0,0 +1,91 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MGroup.FEM.Structural.Embedding; +using MGroup.MSolve.Discretization; +using MGroup.MSolve.Discretization.Embedding; +using MGroup.MSolve.Discretization.Entities; + +namespace MGroup.FEM.Structural.Embedding +{ + public class EmbeddedBeam3DGrouping + { + private readonly Model model; + private readonly IEnumerable hostGroup; + private readonly IEnumerable embeddedGroup; + private readonly bool hasEmbeddedRotations = false; + private readonly int skip; + public IEnumerable HostGroup { get { return hostGroup; } } + public IEnumerable EmbeddedGroup { get { return embeddedGroup; } } + + public static EmbeddedBeam3DGrouping CreateFullyBonded(Model model, IEnumerable hostGroup, IEnumerable embeddedGroup, bool hasEmbeddedRotations) + { + return new EmbeddedBeam3DGrouping(model, hostGroup, embeddedGroup, hasEmbeddedRotations, 0); + } + + public static EmbeddedBeam3DGrouping CreateCohesive(Model model, IEnumerable hostGroup, IEnumerable embeddedGroup, bool hasEmbeddedRotations) + { + return new EmbeddedBeam3DGrouping(model, hostGroup, embeddedGroup, hasEmbeddedRotations, 2); + } + + private EmbeddedBeam3DGrouping(Model model, IEnumerable hostGroup, IEnumerable embeddedGroup, bool hasEmbeddedRotations, int skip) + { + this.model = model; + this.hostGroup = hostGroup; + this.embeddedGroup = embeddedGroup; + this.hasEmbeddedRotations = hasEmbeddedRotations; + this.skip = skip; + hostGroup.Select(e => e).Distinct().ToList().ForEach(et => + { + if (!(et is IEmbeddedHostElement)) + throw new ArgumentException("EmbeddedGrouping: One or more elements of host group does NOT implement IEmbeddedHostElement."); + }); + embeddedGroup.Select(e => e).Distinct().ToList().ForEach(et => + { + if (!(et is IEmbeddedElement)) + throw new ArgumentException("EmbeddedGrouping: One or more elements of embedded group does NOT implement IEmbeddedElement."); + }); + UpdateNodesBelongingToEmbeddedElements(); + } + + private void UpdateNodesBelongingToEmbeddedElements() + { + IEmbeddedDOFInHostTransformationVector transformer; + if (hasEmbeddedRotations) + transformer = new Hexa8TranslationAndRotationTransformationVector(); + else + transformer = new Hexa8LAndNLTranslationTransformationVector(); + + foreach (var embeddedElement in embeddedGroup) + { + var elType = (IEmbeddedElement)embeddedElement; + foreach (var node in embeddedElement.Nodes.Skip(skip)) + { + var embeddedNodes = hostGroup + .Select(e => ((IEmbeddedHostElement)e).BuildHostElementEmbeddedNode(e, node, transformer)) + .Where(e => e != null); + foreach (var embeddedNode in embeddedNodes) + { + if (elType.EmbeddedNodes.Count(x => x.Node == embeddedNode.Node) == 0) + elType.EmbeddedNodes.Add(embeddedNode); + + // Update embedded node information for elements that are not inside the embedded group but contain an embedded node. + foreach (var element in model.ElementsDictionary.Values.Except(embeddedGroup)) + if (element is IEmbeddedElement && element.Nodes.Contains(embeddedNode.Node)) + { + var currentElementType = (IEmbeddedElement)element; + if (!currentElementType.EmbeddedNodes.Contains(embeddedNode)) + { + currentElementType.EmbeddedNodes.Add(embeddedNode); + element.DofEnumerator = new BeamElementEmbedder(model, element, transformer); + } + } + } + } + embeddedElement.DofEnumerator = new BeamElementEmbedder(model, embeddedElement, transformer); + } + } + } +} diff --git a/src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalAbstract.cs b/src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalAbstract.cs new file mode 100644 index 0000000..511689a --- /dev/null +++ b/src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalAbstract.cs @@ -0,0 +1,636 @@ +using System; +using System.Collections.Generic; + +using MGroup.Constitutive.Structural; +using MGroup.Constitutive.Structural.Continuum; +using MGroup.Constitutive.Structural.Line; +using MGroup.FEM.Structural.Line; +using MGroup.LinearAlgebra.Matrices; +using MGroup.MSolve.Discretization; +using MGroup.MSolve.Discretization.Dofs; +using MGroup.MSolve.Discretization.Embedding; +using MGroup.MSolve.Discretization.Entities; + +namespace MGroup.FEM.Structural.Line +{ + public abstract class SupportiveCohesiveBeam3DCorotationalAbstract //: IFiniteElement, IEmbeddedElement + { + protected static readonly int NATURAL_DEFORMATION_COUNT = 6; + protected static readonly int FREEDOM_DEGREE_COUNT = 12; + protected static readonly int AXIS_COUNT = 3; + protected static readonly int NODE_COUNT = 2; + + protected readonly static IDofType[] nodalDOFTypes = new IDofType[] { StructuralDof.TranslationX, StructuralDof.TranslationY, StructuralDof.TranslationZ, StructuralDof.RotationX, StructuralDof.RotationY, StructuralDof.RotationZ }; + protected readonly static IDofType[][] dofTypes = new IDofType[][] { nodalDOFTypes, nodalDOFTypes }; + private static readonly IDofType[][] dofs = new IDofType[][] { nodalDOFTypes, nodalDOFTypes }; + //protected static final List> FREEDOM_DEGREE_TYPES = + // Collections.nCopies(NODE_COUNT, FreedomDegreeTypeSets.X_Y_Z_ROTX_ROTY_ROTZ); + protected IElementDofEnumerator dofEnumerator = new GenericDofEnumerator(); + protected readonly IIsotropicContinuumMaterial3D material; + protected readonly IList nodes; + protected readonly double density; + protected BeamSection3D beamSection; + protected readonly double initialLength; + public double currentLength; + public Matrix currentRotationMatrix; + protected double[] naturalDeformations; + protected double[] beamAxisX; + protected double[] beamAxisY; + protected double[] beamAxisZ; + private readonly List embeddedNodes = new List(); + + public double RayleighAlpha { get; set; } + public double RayleighBeta { get; set; } + + protected SupportiveCohesiveBeam3DCorotationalAbstract(IList nodes, IIsotropicContinuumMaterial3D material, double density, + BeamSection3D beamSection) + { + this.nodes = nodes; + this.material = material; + this.density = density; + this.beamSection = beamSection; + this.initialLength = Math.Sqrt(Math.Pow(nodes[0].X - nodes[1].X, 2) + Math.Pow(nodes[0].Y - nodes[1].Y, 2) + + Math.Pow(nodes[0].Z - nodes[1].Z, 2)); + this.currentLength = this.initialLength; + this.currentRotationMatrix = Matrix.CreateZero(AXIS_COUNT, AXIS_COUNT); + this.naturalDeformations = new double[NATURAL_DEFORMATION_COUNT]; + this.beamAxisX = new double[AXIS_COUNT]; + this.beamAxisY = new double[AXIS_COUNT]; + this.beamAxisZ = new double[AXIS_COUNT]; + } + + public int ID => 100; + public ElementDimensions ElementDimensions => ElementDimensions.ThreeD; + public IList EmbeddedNodes { get { return embeddedNodes; } } + //public bool MaterialModified => material.Modified; + public IElementDofEnumerator DofEnumerator + { + get { return dofEnumerator; } + set { dofEnumerator = value; } + } + + public abstract void SaveGeometryState(); + public abstract void UpdateState(double[] incrementalNodeDisplacements); + + private Matrix CalculateBlockRotationMatrix() + { + Matrix blockRotationMatrix = Matrix.CreateZero(FREEDOM_DEGREE_COUNT, FREEDOM_DEGREE_COUNT); + int totalBlocks = 4; + int blockSize = 3; + double R11 = this.currentRotationMatrix[0, 0]; + double R12 = this.currentRotationMatrix[0, 1]; + double R13 = this.currentRotationMatrix[0, 2]; + double R21 = this.currentRotationMatrix[1, 0]; + double R22 = this.currentRotationMatrix[1, 1]; + double R23 = this.currentRotationMatrix[1, 2]; + double R31 = this.currentRotationMatrix[2, 0]; + double R32 = this.currentRotationMatrix[2, 1]; + double R33 = this.currentRotationMatrix[2, 2]; + + for (int block = 0; block < totalBlocks; block++) + { + int cellDistanceCount = block * blockSize; + + blockRotationMatrix[cellDistanceCount, cellDistanceCount] = R11; + blockRotationMatrix[cellDistanceCount, cellDistanceCount + 1] = R12; + blockRotationMatrix[cellDistanceCount, cellDistanceCount + 2] = R13; + + blockRotationMatrix[cellDistanceCount + 1, cellDistanceCount] = R21; + blockRotationMatrix[cellDistanceCount + 1, cellDistanceCount + 1] = R22; + blockRotationMatrix[cellDistanceCount + 1, cellDistanceCount + 2] = R23; + + blockRotationMatrix[cellDistanceCount + 2, cellDistanceCount] = R31; + blockRotationMatrix[cellDistanceCount + 2, cellDistanceCount + 1] = R32; + blockRotationMatrix[cellDistanceCount + 2, cellDistanceCount + 2] = R33; + } + return blockRotationMatrix; + } + + /** + * Calculates the constitutive stiffness of the element. + * + * @return The constitutive stiffness + */ + private Matrix CalculateConstitutiveStiffness() + { + var constitutiveStiffness = SymmetricMatrix.CreateZero(FREEDOM_DEGREE_COUNT); + var materialCopy = (ElasticMaterial3D)this.material; + double E = materialCopy.YoungModulus; + double G = E / (2d * (1d + materialCopy.PoissonRatio)); + double Iy = this.beamSection.InertiaY; + double Iz = this.beamSection.InertiaZ; + double J = this.beamSection.TorsionalInertia; + double A = this.beamSection.Area; + double Ay = this.beamSection.EffectiveAreaY; + double Az = this.beamSection.EffectiveAreaZ; + double l = this.currentLength; + double lSqruared = l * l; + double lCubed = l * l * l; + double phiY = (12.0 * E * Iy) / (l * l * G * Az); + double phiZ = (12.0 * E * Iz) / (l * l * G * Ay); + double psiY = 1.0 / (1.0 + phiY); + double psiZ = 1.0 / (1.0 + phiZ); + double EAOverL = (E * A) / l; + double GJOverL = (G * J) / l; + double psiZ_E_Iz_12_OverlCubed = (12.0 * psiZ * E * Iz) / lCubed; + double psiZ_E_Iz_6_OverlSquared = (6.0 * psiZ * E * Iz) / lSqruared; + double psiY_E_12_Iy_OverlCubed = (12.0 * psiY * E * Iy) / lCubed; + double psiY_E_6_Iy_OverlSquared = (6.0 * psiY * E * Iy) / lSqruared; + double psiZ_3_Plus_1_E_Iz_Overl = (((3.0 * psiZ) + 1.0) * E * Iz) / l; + double psiZ_3_Minus_1_E_Iz_Overl = (((3.0 * psiZ) - 1.0) * E * Iz) / l; + double psiY_3_Plus_1_E_Iy_Overl = (((3.0 * psiY) + 1.0) * E * Iy) / l; + double psiY_3_Minus_1_E_Iy_Overl = (((3.0 * psiY) - 1.0) * E * Iy) / l; + + constitutiveStiffness[0, 0] = EAOverL; + constitutiveStiffness[0, 6] = -EAOverL; + + constitutiveStiffness[1, 1] = psiZ_E_Iz_12_OverlCubed; + constitutiveStiffness[1, 5] = psiZ_E_Iz_6_OverlSquared; + constitutiveStiffness[1, 7] = -psiZ_E_Iz_12_OverlCubed; + constitutiveStiffness[1, 11] = psiZ_E_Iz_6_OverlSquared; + + constitutiveStiffness[2, 2] = psiY_E_12_Iy_OverlCubed; + constitutiveStiffness[2, 4] = -psiY_E_6_Iy_OverlSquared; + constitutiveStiffness[2, 8] = -psiY_E_12_Iy_OverlCubed; + constitutiveStiffness[2, 10] = -psiY_E_6_Iy_OverlSquared; + + constitutiveStiffness[3, 3] = GJOverL; + constitutiveStiffness[3, 9] = -GJOverL; + + constitutiveStiffness[4, 4] = psiY_3_Plus_1_E_Iy_Overl; + constitutiveStiffness[4, 8] = psiZ_E_Iz_6_OverlSquared; + constitutiveStiffness[4, 10] = psiY_3_Minus_1_E_Iy_Overl; + + constitutiveStiffness[5, 5] = psiZ_3_Plus_1_E_Iz_Overl; + constitutiveStiffness[5, 7] = -psiY_E_6_Iy_OverlSquared; + constitutiveStiffness[5, 11] = psiZ_3_Minus_1_E_Iz_Overl; + + constitutiveStiffness[6, 6] = EAOverL; + + constitutiveStiffness[7, 7] = psiZ_E_Iz_12_OverlCubed; + constitutiveStiffness[7, 11] = -psiZ_E_Iz_6_OverlSquared; + + constitutiveStiffness[8, 8] = psiY_E_12_Iy_OverlCubed; + constitutiveStiffness[8, 10] = psiY_E_6_Iy_OverlSquared; + + constitutiveStiffness[9, 9] = GJOverL; + + constitutiveStiffness[10, 10] = psiY_3_Plus_1_E_Iy_Overl; + + constitutiveStiffness[11, 11] = psiZ_3_Plus_1_E_Iz_Overl; + + return constitutiveStiffness.CopyToFullMatrix(); + } + + /** + * Calculates the forces in the global coordinate system. + * + * @return The forces in the global coordinate system + */ + private double[] CalculateForcesInGlobalSystem() + { + Matrix transformationMatrix = this.CalculateNaturalToGlobalTransormMatrix(); + double[] forcesNatural = this.CalculateForcesInNaturalSystem(); + double[] forcesGlobal = transformationMatrix.Multiply(forcesNatural); + return forcesGlobal; + } + + /** + * Calculates the forces in the local coordinate system. + * + * @return The forces in the local coordinate system + */ + private double[] CalculateForcesInLocalSystem() + { + Matrix naturalToLocal = this.CalculateNaturalToLocalTranformMatrix(); + double[] naturalForces = this.CalculateForcesInNaturalSystem(); + double[] forcesLocal = naturalToLocal.Multiply(naturalForces); + return forcesLocal; + } + + /** + * Calculates forces in the natural coordinate system. + * + * @return The forces in the natural coordinate system + */ + private double[] CalculateForcesInNaturalSystem() + { + var forcesNatural = new double[NATURAL_DEFORMATION_COUNT]; + var materialCopy = (ElasticMaterial3D)this.material; + double E = materialCopy.YoungModulus; + double G = E / (2d * (1d + materialCopy.PoissonRatio)); + double Iy = this.beamSection.InertiaY; + double Iz = this.beamSection.InertiaZ; + double J = this.beamSection.TorsionalInertia; + double A = this.beamSection.Area; + double Ay = this.beamSection.EffectiveAreaY; + double Az = this.beamSection.EffectiveAreaZ; + double l = this.currentLength; + double phiY = (12.0 * E * Iy) / (l * l * G * Az); + double phiZ = (12.0 * E * Iz) / (l * l * G * Ay); + double psiY = 1.0 / (1.0 + phiY); + double psiZ = 1.0 / (1.0 + phiZ); + forcesNatural[NaturalDeformationMode3D.TORSION] = + (G * J * this.naturalDeformations[NaturalDeformationMode3D.TORSION]) / l; + forcesNatural[NaturalDeformationMode3D.SYMMETRIC_BENDING_Y] = + (E * Iy * this.naturalDeformations[NaturalDeformationMode3D.SYMMETRIC_BENDING_Y]) / l; + forcesNatural[NaturalDeformationMode3D.SYMMETRIC_BENDING_Z] = + (E * Iz * this.naturalDeformations[NaturalDeformationMode3D.SYMMETRIC_BENDING_Z]) / l; + forcesNatural[NaturalDeformationMode3D.EXTENSION] = + (E * A * this.naturalDeformations[NaturalDeformationMode3D.EXTENSION]) / l; + forcesNatural[NaturalDeformationMode3D.ANTISYMMETRIC_BENDING_Y] = + (3.0 * psiY * E * Iy * this.naturalDeformations[NaturalDeformationMode3D.ANTISYMMETRIC_BENDING_Y]) / l; + forcesNatural[NaturalDeformationMode3D.ANTISYMMETRIC_BENDING_Z] = + (3.0 * psiZ * E * Iz * this.naturalDeformations[NaturalDeformationMode3D.ANTISYMMETRIC_BENDING_Z]) / l; + return forcesNatural; + } + + /** + * Calculates the geometric stiffness of the element. + * + * @return The geometric stiffness + */ + private Matrix CalculateGeometricStiffness() + { + var geometricStiffness = SymmetricMatrix.CreateZero(FREEDOM_DEGREE_COUNT); + var forcesInNaturalSystem = this.CalculateForcesInNaturalSystem(); + var forcesInLocalSystem = this.CalculateForcesInLocalSystem(); + double torsionalMoment = forcesInNaturalSystem[NaturalDeformationMode3D.TORSION]; + double axialForce = forcesInNaturalSystem[NaturalDeformationMode3D.EXTENSION]; + double momentY_A = forcesInLocalSystem[4]; + double momentZ_A = forcesInLocalSystem[5]; + double momentY_B = forcesInLocalSystem[10]; + double momentZ_B = forcesInLocalSystem[11]; + double length = this.currentLength; + double Qy = -(momentZ_A + momentZ_B) / length; + double Qz = (momentY_A + momentY_B) / length; + double Qy_Over_L = Qy / length; + double Qz_Over_L = Qz / length; + double Qy_L_Over_6 = (Qy * length) / 6.0; + double Qz_L_Over_6 = (Qz * length) / 6.0; + double N_6_Over_5_L = (6.0 * axialForce) / (5.0 * length); + double N_Over_10 = axialForce / 10.0; + double N_L_4_Over_30 = (4.0 * axialForce * length) / 30.0; + double N_L_Over_30 = (axialForce * length) / 30.0; + double MyA_Over_L = momentY_A / length; + double MyB_Over_L = momentY_B / length; + double MzA_Over_L = momentZ_A / length; + double MzB_Over_L = momentZ_B / length; + double M_Over_L = torsionalMoment / length; + double M_3_Over_6 = (3.0 * torsionalMoment) / 6.0; + + geometricStiffness[0, 1] = -Qy_Over_L; + geometricStiffness[0, 2] = -Qz_Over_L; + geometricStiffness[0, 7] = Qy_Over_L; + geometricStiffness[0, 8] = Qz_Over_L; + + geometricStiffness[1, 1] = N_6_Over_5_L; + geometricStiffness[1, 3] = MyA_Over_L; + geometricStiffness[1, 4] = M_Over_L; + geometricStiffness[1, 5] = N_Over_10; + geometricStiffness[1, 6] = Qy_Over_L; + geometricStiffness[1, 7] = -N_6_Over_5_L; + geometricStiffness[1, 9] = MyB_Over_L; + geometricStiffness[1, 10] = -M_Over_L; + geometricStiffness[1, 11] = N_Over_10; + + geometricStiffness[2, 2] = N_6_Over_5_L; + geometricStiffness[2, 3] = MzA_Over_L; + geometricStiffness[2, 4] = -N_Over_10; + geometricStiffness[2, 5] = M_Over_L; + geometricStiffness[2, 6] = Qz_Over_L; + geometricStiffness[2, 8] = -N_6_Over_5_L; + geometricStiffness[2, 9] = MzB_Over_L; + geometricStiffness[2, 10] = -N_Over_10; + geometricStiffness[2, 11] = -M_Over_L; + + geometricStiffness[3, 4] = ((-2.0 * momentZ_A) + momentZ_B) / 6.0; + geometricStiffness[3, 5] = ((2.0 * momentY_A) - momentY_B) / 6.0; + geometricStiffness[3, 7] = -MyA_Over_L; + geometricStiffness[3, 8] = -MzA_Over_L; + geometricStiffness[3, 10] = Qy_L_Over_6; + geometricStiffness[3, 11] = Qz_L_Over_6; + + geometricStiffness[4, 4] = N_L_4_Over_30; + geometricStiffness[4, 7] = -M_Over_L; + geometricStiffness[4, 8] = +N_Over_10; + geometricStiffness[4, 9] = Qy_L_Over_6; + geometricStiffness[4, 10] = -N_L_Over_30; + geometricStiffness[4, 11] = M_3_Over_6; + + geometricStiffness[5, 5] = N_L_4_Over_30; + geometricStiffness[5, 7] = -N_Over_10; + geometricStiffness[5, 8] = -M_Over_L; + geometricStiffness[5, 9] = Qz_L_Over_6; + geometricStiffness[5, 10] = -M_3_Over_6; + geometricStiffness[5, 11] = -N_L_Over_30; + + geometricStiffness[6, 7] = -Qy_Over_L; + geometricStiffness[6, 8] = -Qz_Over_L; + + geometricStiffness[7, 7] = N_6_Over_5_L; + geometricStiffness[7, 9] = -MyB_Over_L; + geometricStiffness[7, 10] = M_Over_L; + geometricStiffness[7, 11] = -N_Over_10; + + geometricStiffness[8, 8] = N_6_Over_5_L; + geometricStiffness[8, 9] = -MzB_Over_L; + geometricStiffness[8, 10] = N_Over_10; + geometricStiffness[8, 11] = M_Over_L; + + geometricStiffness[9, 10] = ((-2.0 * momentZ_B) + momentZ_A) / 6.0; + geometricStiffness[9, 11] = ((+2.0 * momentY_B) - momentY_A) / 6.0; + + geometricStiffness[10, 10] = N_L_4_Over_30; + + geometricStiffness[11, 11] = N_L_4_Over_30; + + return geometricStiffness.CopyToFullMatrix(); + } + + /** + * Calculates the stiffness matrix in the local coordinate system. + * + * @return The stiffness matrix in the local coordinate system. + */ + private Matrix CalculateLocalStiffnessMatrix() + { + Matrix constitutivePart = this.CalculateConstitutiveStiffness(); + Matrix geometricPart = this.CalculateGeometricStiffness(); + constitutivePart.AddIntoThis(geometricPart); + return constitutivePart; + } + + /** + * Calculates the transformation matrix from natural to local coordinate system. + * + * @return The natural to local transformation matrix + */ + private Matrix CalculateNaturalToGlobalTransormMatrix() + { + var transformMatrix = Matrix.CreateZero(FREEDOM_DEGREE_COUNT, NATURAL_DEFORMATION_COUNT); + double L = this.currentLength; + double nx1 = this.beamAxisX[0]; + double nx2 = this.beamAxisX[1]; + double nx3 = this.beamAxisX[2]; + double ny1 = this.beamAxisY[0]; + double ny2 = this.beamAxisY[1]; + double ny3 = this.beamAxisY[2]; + double nz1 = this.beamAxisZ[0]; + double nz2 = this.beamAxisZ[1]; + double nz3 = this.beamAxisZ[2]; + + transformMatrix[0, 3] = -nx1; + transformMatrix[0, 4] = (-2.0 * nz1) / L; + transformMatrix[0, 5] = (2.0 * ny1) / L; + + transformMatrix[1, 3] = -nx2; + transformMatrix[1, 4] = (-2.0 * nz2) / L; + transformMatrix[1, 5] = (2.0 * ny2) / L; + + transformMatrix[2, 3] = -nx3; + transformMatrix[2, 4] = (-2.0 * nz3) / L; + transformMatrix[2, 5] = (2.0 * ny3) / L; + + transformMatrix[3, 0] = -nx1; + transformMatrix[3, 1] = -ny1; + transformMatrix[3, 2] = -nz1; + transformMatrix[3, 4] = ny1; + transformMatrix[3, 5] = nz1; + + transformMatrix[4, 0] = -nx2; + transformMatrix[4, 1] = -ny2; + transformMatrix[4, 2] = -nz2; + transformMatrix[4, 4] = ny2; + transformMatrix[4, 5] = nz2; + + transformMatrix[5, 0] = -nx3; + transformMatrix[5, 1] = -ny3; + transformMatrix[5, 2] = -nz3; + transformMatrix[5, 4] = ny3; + transformMatrix[5, 5] = nz3; + + transformMatrix[6, 3] = nx1; + transformMatrix[6, 4] = (2.0 * nz1) / L; + transformMatrix[6, 5] = (-2.0 * ny1) / L; + + transformMatrix[7, 3] = nx2; + transformMatrix[7, 4] = (2.0 * nz2) / L; + transformMatrix[7, 5] = (-2.0 * ny2) / L; + + transformMatrix[8, 3] = nx3; + transformMatrix[8, 4] = (2.0 * nz3) / L; + transformMatrix[8, 5] = (-2.0 * ny3) / L; + + transformMatrix[9, 0] = nx1; + transformMatrix[9, 1] = ny1; + transformMatrix[9, 2] = nz1; + transformMatrix[9, 4] = ny1; + transformMatrix[9, 5] = nz1; + + transformMatrix[10, 0] = nx2; + transformMatrix[10, 1] = ny2; + transformMatrix[10, 2] = nz2; + transformMatrix[10, 4] = ny2; + transformMatrix[10, 5] = nz2; + + transformMatrix[11, 0] = nx3; + transformMatrix[11, 1] = ny3; + transformMatrix[11, 2] = nz3; + transformMatrix[11, 4] = ny3; + transformMatrix[11, 5] = nz3; + + return transformMatrix; + } + + /** + * Calculates the transformation matrix from natural to local coordinate system. + * + * @return The natural to local transformation matrix + */ + private Matrix CalculateNaturalToLocalTranformMatrix() + { + var transformMatrix = Matrix.CreateZero(FREEDOM_DEGREE_COUNT, NATURAL_DEFORMATION_COUNT); + double L = this.currentLength; + double nx1 = 1.0; + double ny2 = 1.0; + double nz3 = 1.0; + + transformMatrix[0, 3] = -nx1; + + transformMatrix[1, 5] = (2.0 * ny2) / L; + + transformMatrix[2, 4] = (-2.0 * nz3) / L; + + transformMatrix[3, 0] = -nx1; + + transformMatrix[4, 1] = -ny2; + transformMatrix[4, 4] = ny2; + + transformMatrix[5, 2] = -nz3; + transformMatrix[5, 5] = nz3; + + transformMatrix[6, 3] = nx1; + + transformMatrix[7, 5] = (-2.0 * ny2) / L; + + transformMatrix[8, 4] = (2.0 * nz3) / L; + + transformMatrix[9, 0] = nx1; + + transformMatrix[10, 1] = ny2; + transformMatrix[10, 4] = ny2; + + transformMatrix[11, 2] = nz3; + transformMatrix[11, 5] = nz3; + + return transformMatrix; + } + + public IList> GetElementDOFTypes(IElementType element) => dofTypes; + + public IMatrix StiffnessMatrix(IElementType element) + { + Matrix rotationMatrixBlock = this.CalculateBlockRotationMatrix(); + Matrix localStiffnessMatrix = this.CalculateLocalStiffnessMatrix(); + Matrix s = rotationMatrixBlock.MultiplyRight(localStiffnessMatrix).MultiplyRight(rotationMatrixBlock, false, true); + return dofEnumerator.GetTransformedMatrix(s); + } + + public IMatrix MassMatrix(IElementType element) + { + //throw new NotImplementedException(); + double area = beamSection.Area; + double inertiaY = beamSection.InertiaY; + double inertiaZ = beamSection.InertiaZ; + double x2 = Math.Pow(element.Nodes[1].X - element.Nodes[0].X, 2); + double y2 = Math.Pow(element.Nodes[1].Y - element.Nodes[0].Y, 2); + double z2 = Math.Pow(element.Nodes[1].Z - element.Nodes[0].Z, 2); + double L = Math.Sqrt(x2 + y2 + z2); + double fullMass = density * area * L; + + var massMatrix = Matrix.CreateZero(FREEDOM_DEGREE_COUNT, FREEDOM_DEGREE_COUNT); + massMatrix[0, 0] = (1.0 / 3.0) * fullMass; + massMatrix[0, 6] = (1.0 / 6.0) * fullMass; + + massMatrix[1, 1] = (13.0 / 35.0) * fullMass; + massMatrix[1, 5] = (11.0 * L / 210.0) * fullMass; + massMatrix[1, 7] = (9.0 / 70.0) * fullMass; + massMatrix[1, 11] = -(13.0 * L / 420.0) * fullMass; + + massMatrix[2, 2] = (13.0 / 35.0) * fullMass; + massMatrix[2, 4] = -(11.0 * L / 210.0) * fullMass; + massMatrix[2, 8] = (9.0 / 70.0) * fullMass; + massMatrix[2, 10] = -(13.0 * L / 420.0) * fullMass; + + massMatrix[3, 3] = ((inertiaY + inertiaZ) / (3.0 * area)) * fullMass; + massMatrix[3, 9] = ((inertiaY + inertiaZ) / (6.0 * area)) * fullMass; + + massMatrix[4, 4] = ((L * L) / 105.0) * fullMass; + massMatrix[4, 8] = -(13.0 * L / 420.0) * fullMass; + massMatrix[4, 10] = -((L * L) / 105.0) * fullMass; + + massMatrix[5, 5] = ((L * L) / 105.0) * fullMass; + massMatrix[5, 7] = (13.0 * L / 420.0) * fullMass; + massMatrix[5, 11] = -((L * L) / 105.0) * fullMass; + + massMatrix[6, 6] = (1.0 / 3.0) * fullMass; + + massMatrix[7, 7] = (13.0 / 35.0) * fullMass; + massMatrix[7, 11] = -(11.0 * L / 210.0) * fullMass; + + massMatrix[8, 8] = (13.0 / 35.0) * fullMass; + massMatrix[8, 10] = (11.0 * L / 210.0) * fullMass; + + massMatrix[9, 9] = ((inertiaY + inertiaZ) / (3.0 * area)) * fullMass; + + massMatrix[10, 10] = ((L * L) / 105.0) * fullMass; + + massMatrix[11, 11] = ((L * L) / 105.0) * fullMass; + + massMatrix[6, 0] = (1.0 / 6.0) * fullMass; + massMatrix[5, 1] = (11.0 * L / 210.0) * fullMass; + massMatrix[7, 1] = (9.0 / 70.0) * fullMass; + massMatrix[11, 1] = -(13.0 * L / 420.0) * fullMass; + massMatrix[4, 2] = -(11.0 * L / 210.0) * fullMass; + massMatrix[8, 2] = (9.0 / 70.0) * fullMass; + massMatrix[10, 2] = -(13.0 * L / 420.0) * fullMass; + massMatrix[9, 3] = ((inertiaY + inertiaZ) / (6.0 * area)) * fullMass; + massMatrix[8, 4] = -(13.0 * L / 420.0) * fullMass; + massMatrix[10, 4] = -((L * L) / 105.0) * fullMass; + massMatrix[7, 5] = (13.0 * L / 420.0) * fullMass; + massMatrix[11, 5] = -((L * L) / 105.0) * fullMass; + massMatrix[11, 7] = -(11.0 * L / 210.0) * fullMass; + massMatrix[10, 8] = (11.0 * L / 210.0) * fullMass; + + return dofEnumerator.GetTransformedMatrix(massMatrix); + } + + public IMatrix DampingMatrix(IElementType element) + { + IMatrix k = StiffnessMatrix(element); + IMatrix m = MassMatrix(element); + k.LinearCombinationIntoThis(RayleighBeta, m, RayleighAlpha); + return dofEnumerator.GetTransformedMatrix(k); + } + + //public void ResetMaterialModified() => material.ResetModified(); + + public Tuple CalculateStresses(double[] localdDisplacements) + { + UpdateState(dofEnumerator.GetTransformedDisplacementsVector(localdDisplacements)); + //TODO: Should calculate strains and update material as well + //material.UpdateMaterial(strains); + //TODO: Should calculate stresses as well + return new Tuple(new double[FREEDOM_DEGREE_COUNT], new double[FREEDOM_DEGREE_COUNT]); + } + + public double[] CalculateForces(IElementType element, double[] localDisplacements, double[] localdDisplacements) + { + double[] internalForces = this.CalculateForcesInGlobalSystem(); + return dofEnumerator.GetTransformedForcesVector(internalForces); + } + + public double[] CalculateForcesForLogging(IElementType element, double[] localDisplacements) + { + throw new NotImplementedException(); + } + + //public double[] CalculateAccelerationForces(IElementType element, IList loads) + //{ + // var accelerations = new double[6]; + // IMatrix massMatrix = MassMatrix(element); + + // int index = 0; + // foreach (MassAccelerationLoad load in loads) + // foreach (IDofType[] nodalDOFTypes in dofs) + // foreach (IDofType dofType in nodalDOFTypes) + // { + // if (dofType == load.DOF) accelerations[index] += load.Amount; + // index++; + // } + + // return massMatrix.Multiply(accelerations); + //} + + public void SaveMaterialState() + { + SaveGeometryState(); + material.CreateState(); + } + + //public void ClearMaterialState() => material.ClearState(); + + //public void ClearMaterialStresses() => material.ClearStresses(); + + public Dictionary GetInternalNodalDOFs(IElementType element, Node node) + { + throw new NotImplementedException(); + } + + public double[] GetLocalDOFValues(IElementType hostElement, double[] hostDOFValues) + { + throw new NotImplementedException(); + } + } +} diff --git a/src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalQuaternion.cs b/src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalQuaternion.cs new file mode 100644 index 0000000..c32e296 --- /dev/null +++ b/src/MGroup.FEM.Structural/Line/SupportiveCohesiveBeam3DCorotationalQuaternion.cs @@ -0,0 +1,211 @@ +using System; +using System.Collections.Generic; + +using MGroup.Constitutive.Structural.Continuum; +using MGroup.Constitutive.Structural.Line; +using MGroup.LinearAlgebra.Matrices; +using MGroup.MSolve.Discretization.Entities; +using MGroup.MSolve.Geometry; + +namespace MGroup.FEM.Structural.Line +{ + public class SupportiveCohesiveBeam3DCorotationalQuaternion : SupportiveCohesiveBeam3DCorotationalAbstract + { + private readonly double[] lastDisplacements; + private readonly double[] currentDisplacements; + private readonly double[] displacementsOfCurrentIncrement; + private readonly double[] initialAxisY; + private readonly double[] initialAxisZ; + private readonly double[] initialAxisX; + private Quaternion quaternionCurrentNodeA; + private Quaternion quaternionCurrentNodeB; + private readonly double[] currentBeamAxis; + private Quaternion quaternionLastNodeA; + private Quaternion quaternionLastNodeB; + + /** + * Creates a new instance of {@link Beam3DCorotationalIncremental} class. + * + * @param nodes + * The element nodes + * @param material + * The element material + * @param density + * The element density + * @param beamSection + * The beam section. + */ + public SupportiveCohesiveBeam3DCorotationalQuaternion(IList nodes, IIsotropicContinuumMaterial3D material, double density, + BeamSection3D beamSection) + : base(nodes, material, density, beamSection) + { + this.displacementsOfCurrentIncrement = new double[FREEDOM_DEGREE_COUNT]; + this.lastDisplacements = new double[FREEDOM_DEGREE_COUNT]; + this.currentDisplacements = new double[FREEDOM_DEGREE_COUNT]; + this.quaternionLastNodeA = Quaternion.OfZeroAngle(); + this.quaternionLastNodeB = Quaternion.OfZeroAngle(); + this.quaternionCurrentNodeA = Quaternion.OfZeroAngle(); + this.quaternionCurrentNodeB = Quaternion.OfZeroAngle(); + this.initialAxisX = new double[AXIS_COUNT]; + this.initialAxisY = new double[AXIS_COUNT]; + this.initialAxisZ = new double[AXIS_COUNT]; + this.currentBeamAxis = new double[AXIS_COUNT]; + this.InitializeElementAxes(); + } + + public override void SaveGeometryState() + { + displacementsOfCurrentIncrement.Clear(); + //displacementsOfCurrentIncrement.ScaleIntoThis(0d); + lastDisplacements.CopyFrom(currentDisplacements); + + double[] qA = quaternionCurrentNodeA.VectorPart.Copy(); + double[] qB = quaternionCurrentNodeB.VectorPart.Copy(); + this.quaternionLastNodeA = new Quaternion(quaternionCurrentNodeA.ScalarPart, qA); + this.quaternionLastNodeB = new Quaternion(quaternionCurrentNodeB.ScalarPart, qB); + } + + public override void UpdateState(double[] incrementalNodeDisplacements) + { + //displacementsOfCurrentIncrement.Add(new Vector(incrementalNodeDisplacements)); + //lastDisplacements.CopyTo(currentDisplacements.Data, 0); + //currentDisplacements.Add(displacementsOfCurrentIncrement); + + //this.currentDisplacements.addIntoThis(this.lastDisplacements, incrementalNodeDisplacements); + currentDisplacements.CopyFrom(lastDisplacements); + currentDisplacements.AddIntoThis(incrementalNodeDisplacements); + displacementsOfCurrentIncrement.CopyFrom(incrementalNodeDisplacements); + + var incrementalRotationsA = new double[AXIS_COUNT]; + var incrementalRotationsB = new double[AXIS_COUNT]; + + incrementalRotationsA[0] = displacementsOfCurrentIncrement[3]; + incrementalRotationsA[1] = displacementsOfCurrentIncrement[4]; + incrementalRotationsA[2] = displacementsOfCurrentIncrement[5]; + incrementalRotationsB[0] = displacementsOfCurrentIncrement[9]; + incrementalRotationsB[1] = displacementsOfCurrentIncrement[10]; + incrementalRotationsB[2] = displacementsOfCurrentIncrement[11]; + + double[] qA = quaternionLastNodeA.VectorPart.Copy(); + double[] qB = quaternionLastNodeB.VectorPart.Copy(); + this.quaternionCurrentNodeA = new Quaternion(quaternionLastNodeA.ScalarPart, qA); + this.quaternionCurrentNodeB = new Quaternion(quaternionLastNodeB.ScalarPart, qB); + this.quaternionCurrentNodeA.ApplyIncrementalRotation(incrementalRotationsA); + this.quaternionCurrentNodeB.ApplyIncrementalRotation(incrementalRotationsB); + + double scalarA = this.quaternionCurrentNodeA.ScalarPart; + double scalarB = this.quaternionCurrentNodeB.ScalarPart; + var vectorPartA = this.quaternionCurrentNodeA.VectorPart; + var vectorPartB = this.quaternionCurrentNodeB.VectorPart; + double[] sumOfVectorParts = vectorPartA.Copy(); + sumOfVectorParts.AddIntoThis(vectorPartB); + double sumOfVectorPartsNorm = sumOfVectorParts.Norm2(); + double scalarPartDifference = 0.5 * Math.Sqrt(((scalarA + scalarB) * (scalarA + scalarB)) + (sumOfVectorPartsNorm * sumOfVectorPartsNorm)); + double meanRotationScalarPart = (0.5 * (scalarA + scalarB)) / scalarPartDifference; + double[] meanRotationVectorPart = vectorPartA.Copy(); + meanRotationVectorPart.AddIntoThis(vectorPartB); + meanRotationVectorPart.ScaleIntoThis(1d / (2d * scalarPartDifference)); + double[] vectorPartDifference = vectorPartB.Copy(); + vectorPartDifference.ScaleIntoThis(scalarA); + vectorPartDifference.AddIntoThis(vectorPartA.Scale(-scalarB)); + //vectorPartDifference.doPointwise(vectorPartA, DoubleBinaryOps.alphaPlusScaledBeta(-scalarB)); + vectorPartDifference.AddIntoThis(vectorPartA.CrossProduct(vectorPartB)); + vectorPartDifference.ScaleIntoThis(1d / (2d * scalarPartDifference)); + var meanRotationQuaternion = Quaternion.CreateFromIndividualParts(meanRotationScalarPart, meanRotationVectorPart); + this.currentRotationMatrix = meanRotationQuaternion.GetRotationMatrix(); + this.CalculateUpdatedBeamAxis(); + this.UpdateRotationMatrix(); + this.UpdateNaturalDeformations(vectorPartDifference); + + //SaveGeometryState(); + } + + private void CalculateUpdatedBeamAxis() + { + currentRotationMatrix.MultiplyIntoResult(initialAxisX, beamAxisX); + currentRotationMatrix.MultiplyIntoResult(initialAxisY, beamAxisY); + currentRotationMatrix.MultiplyIntoResult(initialAxisZ, beamAxisZ); + + double dX = ((nodes[1].X - nodes[0].X) + this.currentDisplacements[6]) - this.currentDisplacements[0]; + double dY = ((nodes[1].Y - nodes[0].Y) + this.currentDisplacements[7]) - this.currentDisplacements[1]; + double dZ = ((nodes[1].Z - nodes[0].Z) + this.currentDisplacements[8]) - this.currentDisplacements[2]; + this.currentLength = Math.Sqrt((dX * dX) + (dY * dY) + (dZ * dZ)); + //var delta = new[] { dX, dY, dZ }; + this.currentBeamAxis[0] = dX / currentLength + beamAxisX[0]; + this.currentBeamAxis[1] = dY / currentLength + beamAxisX[1]; + this.currentBeamAxis[2] = dZ / currentLength + beamAxisX[2]; + this.currentBeamAxis.ScaleIntoThis(1d / this.currentBeamAxis.Norm2()); + + //this.currentBeamAxis.divideAllIntoThis(delta, this.currentLength); + //this.currentBeamAxis.add(this.beamAxisX); + //this.currentBeamAxis.divideAll(this.currentBeamAxis.normEntrywise(EntrywiseNorms.L2)); + } + + private void InitializeElementAxes() + { + var globalVectorX = new double[AXIS_COUNT]; + var globalVectorY = new double[AXIS_COUNT]; + var globalVectorZ = new double[AXIS_COUNT]; + globalVectorX[0] = 1d; + globalVectorY[1] = 1d; + globalVectorZ[2] = 1d; + double deltaX = nodes[1].X - nodes[0].X; + double deltaY = nodes[1].Y - nodes[0].Y; + double deltaZ = nodes[1].Z - nodes[0].Z; + var elementAxis = new double[3]; + elementAxis[0] = deltaX; + elementAxis[1] = deltaY; + elementAxis[2] = deltaZ; + elementAxis.ScaleIntoThis(1d / elementAxis.Norm2()); + + currentRotationMatrix = RotationMatrix.CalculateRotationMatrix(globalVectorX, elementAxis); + currentRotationMatrix.MultiplyIntoResult(globalVectorX, initialAxisX); + currentRotationMatrix.MultiplyIntoResult(globalVectorY, initialAxisY); + currentRotationMatrix.MultiplyIntoResult(globalVectorZ, initialAxisZ); + beamAxisX.CopyFrom(initialAxisX); + beamAxisY.CopyFrom(initialAxisY); + beamAxisZ.CopyFrom(initialAxisZ); + } + + private void UpdateNaturalDeformations(double[] vectorPartDifference) + { + double extension = this.currentLength - this.initialLength; + double[] symmetricRotation = currentRotationMatrix.Multiply(vectorPartDifference, true); + double[] antisymmetricRotation = + currentRotationMatrix.Multiply(this.beamAxisX.CrossProduct(this.currentBeamAxis), true); + + //currentRotationMatrix.Transpose().Multiply(vectorPartDifference, symmetricRotation.Data); + //currentRotationMatrix.Multiply(vectorPartDifference, symmetricRotation.Data); + //currentRotationMatrix.Multiply(this.beamAxisX ^ this.currentBeamAxis, antisymmetricRotation.Data); + //currentRotationMatrix.Transpose().Multiply(this.beamAxisX ^ this.currentBeamAxis, antisymmetricRotation.Data); + //currentRotationMatrix.MultiplyTranspose2(this.beamAxisX ^ this.currentBeamAxis, antisymmetricRotation.Data); + //final VectorView symmetricRotation = vectorPartDifference.leftMultiplyWithMatrix(this.currentRotationMatrix); + //final VectorView antisymmetricRotation = + // GeometricUtils.crossProduct(this.beamAxisX, this.currentBeamAxis) + // .leftMultiplyWithMatrix(this.currentRotationMatrix); + + this.naturalDeformations[0] = symmetricRotation[0] * 4.0; + this.naturalDeformations[1] = symmetricRotation[1] * 4.0; + this.naturalDeformations[2] = symmetricRotation[2] * 4.0; + this.naturalDeformations[3] = extension; + this.naturalDeformations[4] = antisymmetricRotation[1] * 4.0; + this.naturalDeformations[5] = antisymmetricRotation[2] * 4.0; + } + + private void UpdateRotationMatrix() + { + var rotationMatrixLocal = Matrix.CreateIdentity(AXIS_COUNT); + rotationMatrixLocal.AxpyIntoThis(currentBeamAxis.TensorProduct(currentBeamAxis), -2d); + //rotationMatrixLocal.LinearCombinationGOAT(new[] { -2d }, new[] { Matrix2D.FromVector(currentBeamAxis.Data) * Matrix2D.FromVectorTranspose(currentBeamAxis.Data) }); + + double[] minusBeamAxisX = beamAxisX.Scale(-1.0); + var tempBeamAxisY = beamAxisY.Copy(); + var tempBeamAxisZ = beamAxisZ.Copy(); + rotationMatrixLocal.MultiplyIntoResult(minusBeamAxisX, beamAxisX); + rotationMatrixLocal.MultiplyIntoResult(tempBeamAxisY, beamAxisY); + rotationMatrixLocal.MultiplyIntoResult(tempBeamAxisZ, beamAxisZ); + this.currentRotationMatrix = RotationMatrix.CalculateFromOrthonormalVectors(this.beamAxisX, this.beamAxisY, this.beamAxisZ); + } + + } +} diff --git a/tests/MGroup.FEM.Structural.Tests/Elements/CohesiveBeam3DToBeam3DTest.cs b/tests/MGroup.FEM.Structural.Tests/Elements/CohesiveBeam3DToBeam3DTest.cs new file mode 100644 index 0000000..cf6a981 --- /dev/null +++ b/tests/MGroup.FEM.Structural.Tests/Elements/CohesiveBeam3DToBeam3DTest.cs @@ -0,0 +1,91 @@ +namespace MGroup.FEM.Structural.Tests.Elements +{ + using System; + using System.Collections.Generic; + using System.Linq; + using System.Text; + using System.Threading.Tasks; + using DotNumerics.Optimization.TN; + + using MGroup.Constitutive.Structural.Cohesive; + using MGroup.Constitutive.Structural.Continuum; + using MGroup.Constitutive.Structural.Line; + using MGroup.FEM.Structural.Embedding; + using MGroup.MSolve.Discretization.Entities; + using MGroup.MSolve.Numerics.Integration.Quadratures; + + using Xunit; + + public class CohesiveBeam3DToBeam3DTest + { + [Fact] + private static void RunTest() + { + var nodeSet = new List() + { + new Node(0, -11.261562, 21.900388, 1.716324), + new Node(1, -42.483085, -17.153164, 1.525223), + new Node(2, -11.261562, 21.900388, 1.716324), + new Node(3, -42.483085, -17.153164, 1.525223), + }; + var elementNodesBeam = new List() + { + new Node(0, -11.261562, 21.900388, 1.716324), + new Node(1, -42.483085, -17.153164, 1.525223), + }; + + IIsotropicContinuumMaterial3D matrixMaterial = new ElasticMaterial3D(youngModulus: 3.5, poissonRatio: 0.4); + var K_el = 10; + var K_pl = 1.0; + var T_max = 0.10; + var cohesiveMaterial = new BondSlipMaterial(K_el, K_pl, 100.0, T_max, new double[2], new double[2], 1e-3); + + double area = 694.77; + double inertiaY = 100.18; + double inertiaZ = 100.18; + double torsionalInertia = 68.77; + var beamSection = new BeamSection3D(area, inertiaY, inertiaZ, torsionalInertia, 0, 0); + + double mi = 8.0; + double ni = 8.0; + double thickness_CNT = 0.34; + double a = 0.241; + double diameter_CNT = (a / Math.PI) * Math.Sqrt(Math.Pow(ni, 2) + ni * mi + Math.Pow(mi, 2)); + double radius_CNT = diameter_CNT / 2.0; + double radius_CNT_outer = radius_CNT + (thickness_CNT / 2); + double CntPerimeter = 2.0 * Math.PI * radius_CNT_outer; + var element = new CohesiveBeam3DToBeam3D(nodeSet, cohesiveMaterial, GaussLegendre1D.GetQuadratureWithOrder(2), elementNodesBeam, elementNodesBeam, matrixMaterial, 1, beamSection, CntPerimeter); + + var elementStiffnessMatrix = element.StiffnessMatrix(); + + double[] localDisplacements = + { + -0.22817359334790202, -0.093157552807095936, -0.00032526368182175715, 1.3688506028227088E-05, 0.000259653566909844, + 4.9007290948375906E-05, -0.8453951256523183, 0.067674668134558316, -0.0095639398099483467, -3.1091507972934296E-05, + 0.00035233238898999089, -0.00012733831428290391, -0.23387978688156119, -0.10029516508465312, -0.00036019020640142024, + 1.3688021032242603E-05, 0.00025965431209030315, 4.9005687219101045E-05, -0.83968893211865792, 0.074812280412115614, + -0.0095290132853686309, -3.1091022976958875E-05, 0.00035233164380952571, -0.00012733671055360516, + }; + + Tuple results = element.CalculateResponse(localDisplacements); + + var forces = element.CalculateResponseIntegral(); + + + var expectedStrains = new double[] { 0.0052756089927531976, 6.11479059092248E-05, 1.3659736733557528E-06 }; + var expectedStresses = new double[] { 0.052756089927531974, 0.000611479059092248, 0.00013659736733557528 }; + var expectedForces = new double[] + { + 2.1064736285205248, 2.6349929239668777, 0.0050331374077453254, 2.4373019262327917E-05, -3.7448346514700273E-05, + 8.0593913678409413E-05, -2.1064736285244603, -2.6349929239668777, -0.0050331374078813225, -2.4373018806862645E-05, + 3.7448346817095162E-05, -8.0593914878794527E-05, -2.1064736285205248, -2.6349929239668777, -0.0050331374077453254, + -2.4373019262327917E-05, 3.7448346514700273E-05, -8.0593913678409413E-05, 2.1064736285244603, 2.6349929239668777, + 0.0050331374078813225, 2.4373018806862645E-05, -3.7448346817095162E-05, 8.0593914878794527E-05, + }; + + Assert.Equal(results.Item1, expectedStrains); + Assert.Equal(results.Item2, expectedStresses); + Assert.Equal(forces, expectedForces); + } + } +}