From b29fe66be244a7bba8eac968f7faa66f3ed40f19 Mon Sep 17 00:00:00 2001 From: pmetsis Date: Fri, 14 Oct 2016 08:11:45 +0300 Subject: [PATCH] Quad4 element and simple example Good morning everyone! I am uploading the Quad4 plane stress element. I also had to make some changes to quadrature and material. You can check the element by solving the simple example provided. Since I use Xamarin for Mac (and I haven't figured out the git extension for it yet), I have manually uploaded these files, if something else is missing please notify. Stresses/strains are not yet implemented but I couldn't comment them out for consistency reasons. I will also upload the solution of the simple example by hand or using some Matlab routines of my own in order to verify the code (I have already checked that it works). --- ElasticMaterial2D.cs | 116 +++++++ GaussQuadrature.cs | 167 ++++++++++ IFiniteElementMaterial2D.cs | 18 ++ Quad4.cs | 619 ++++++++++++++++++++++++++++++++++++ Quad4Example.cs | 125 ++++++++ 5 files changed, 1045 insertions(+) create mode 100644 ElasticMaterial2D.cs create mode 100644 GaussQuadrature.cs create mode 100644 IFiniteElementMaterial2D.cs create mode 100644 Quad4.cs create mode 100644 Quad4Example.cs diff --git a/ElasticMaterial2D.cs b/ElasticMaterial2D.cs new file mode 100644 index 000000000..8a2c719a1 --- /dev/null +++ b/ElasticMaterial2D.cs @@ -0,0 +1,116 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using ISAAR.MSolve.PreProcessor.Interfaces; +using ISAAR.MSolve.Matrices.Interfaces; +using ISAAR.MSolve.Matrices; + +namespace ISAAR.MSolve.PreProcessor.Materials +{ + public class ElasticMaterial2D : IFiniteElementMaterial2D + { + private readonly double[] strains = new double[3]; + private readonly double[] stresses = new double[3]; + private double[,] constitutiveMatrix = null; + public double YoungModulus { get; set; } + public double PoissonRatio { get; set; } + public double[] Coordinates { get; set; } + + private double[,] GetConstitutiveMatrix() + { + + //Panos Plane Stress + // + // [ 1 v 0 ] + // [D] = E/(1-v^2) [ v 1 0 ] + // [ 0 0 (1-v)/2 ] + // + + double fE1 = YoungModulus / (1 - PoissonRatio*PoissonRatio); + double fE2 = (1 - PoissonRatio)/2; + double[,] afE = new double[3, 3]; + afE[0, 0] = fE1; + afE[0, 1] = fE1 * PoissonRatio; + //afE[0, 2] = 0; + afE[1, 0] = fE1 * PoissonRatio; + afE[1, 1] = fE1; + //afE[1, 2] = 0; + //afE[2, 0] = 0; + //afE[2, 1] = 0; + afE[2, 2] = fE1*fE2; + + + Vector s = (new Matrix2D(afE)) * (new Vector(strains)); + s.Data.CopyTo(stresses, 0); + + return afE; + } + + #region IFiniteElementMaterial Members + + public int ID + { + get { return 1; } + } + + public bool Modified + { + get { return false; } + } + + public void ResetModified() + { + } + + #endregion + + #region IFiniteElementMaterial3D Members + + public double[] Stresses { get { return stresses; } } + + public IMatrix2D ConstitutiveMatrix + { + get + { + if (constitutiveMatrix == null) UpdateMaterial(new double[3]); + return new Matrix2D(constitutiveMatrix); + } + } + + public void UpdateMaterial(double[] strains) + { + //throw new NotImplementedException(); + + strains.CopyTo(this.strains, 0); + constitutiveMatrix = GetConstitutiveMatrix(); + } + + public void ClearState() + { + //throw new NotImplementedException(); + } + + public void SaveState() + { + //throw new NotImplementedException(); + } + + public void ClearStresses() + { + //throw new NotImplementedException(); + } + + #endregion + + #region ICloneable Members + + public object Clone() + { + return new ElasticMaterial2D() { YoungModulus = this.YoungModulus, PoissonRatio = this.PoissonRatio }; + } + + #endregion + + } +} diff --git a/GaussQuadrature.cs b/GaussQuadrature.cs new file mode 100644 index 000000000..9b6ffa9da --- /dev/null +++ b/GaussQuadrature.cs @@ -0,0 +1,167 @@ +namespace ISAAR.MSolve.PreProcessor.Elements.SupportiveClasses +{ + #region + + using System; + + #endregion + + public class GaussLegendrePoint1D + { + #region Properties + + public double Coordinate { get; set; } + + public double WeightFactor { get; set; } + + #endregion + } + + //Panos Input Start + public class GaussLegendrePoint2D + { + #region Constructors and Destructors + public GaussLegendrePoint2D( + double xi, double eta, double[,] deformationMatrix, double weightFactor) + { + this.Xi = xi; + this.Eta = eta; + this.DeformationMatrix = deformationMatrix; + this.WeightFactor = weightFactor; + } + #endregion + #region Properties + // Panos - I rearranged the appearce order - If you dont like it, change it... + public double Xi { get; private set; } + public double Eta { get; private set; } + public double WeightFactor { get; private set; } + public double[,] DeformationMatrix { get; private set; } + // + #endregion + } + //Panos Input End + public class GaussLegendrePoint3D + { + #region Constructors and Destructors + + public GaussLegendrePoint3D( + double xi, double eta, double zeta, double[,] deformationMatrix, double weightFactor) + { + this.Xi = xi; + this.Eta = eta; + this.Zeta = zeta; + this.DeformationMatrix = deformationMatrix; + this.WeightFactor = weightFactor; + } + + #endregion + + #region Properties + // Panos - I rearranged the appearce order - If you dont like it, change it... + public double Xi { get; private set; } + public double Eta { get; private set; } + public double Zeta { get; private set; } + public double WeightFactor { get; private set; } + public double[,] DeformationMatrix { get; private set; } + // + #endregion + } + + public class GaussQuadrature + { + /* For point coordinates, we encounter the following constants: + * 0.5773502691896 = 1 / Square Root 3 + * 0.7745966692415 = (Square Root 15)/ 5 + * 0.8611363115941 = Square Root( (3 + 2*sqrt(6/5))/7) + * 0.3399810435849 = Square Root( (3 - 2*sqrt(6/5))/7) + * + * For the weights, we encounter the followings constants: + * 0.5555555555556 = 5/9 + * 0.8888888888889 = 8/9 + * 0.3478548451375 = (18 - sqrt30)/36 + * 0.6521451548625 = (18 + sqrt30)/36 + */ + #region Constants and Fields + + private static readonly GaussLegendrePoint1D GaussLegendrePoint1 = new GaussLegendrePoint1D + { + Coordinate = 0.0, WeightFactor = 2.0 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint2A = new GaussLegendrePoint1D + { + Coordinate = -0.5773502691896, WeightFactor = 1.0 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint2B = new GaussLegendrePoint1D + { + Coordinate = 0.5773502691896, WeightFactor = 1.0 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint3A = new GaussLegendrePoint1D + { + Coordinate = -0.7745966692415, WeightFactor = 0.5555555555556 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint3B = new GaussLegendrePoint1D + { + Coordinate = 0.0, WeightFactor = 0.8888888888889 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint3C = new GaussLegendrePoint1D + { + Coordinate = 0.7745966692415, WeightFactor = 0.5555555555556 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint4A = new GaussLegendrePoint1D + { + Coordinate = -0.86113631159416, WeightFactor = 0.3478548451375 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint4B = new GaussLegendrePoint1D + { + Coordinate = -0.3399810435849, WeightFactor = 0.6521451548625 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint4C = new GaussLegendrePoint1D + { + Coordinate = 0.3399810435849, WeightFactor = 0.6521451548625 + }; + + private static readonly GaussLegendrePoint1D GaussLegendrePoint4D = new GaussLegendrePoint1D + { + Coordinate = 0.86113631159416, WeightFactor = 0.3478548451375 + }; + + #endregion + + #region Public Methods + + public static GaussLegendrePoint1D[] GetGaussLegendrePoints(int integrationDegree) + { + if (integrationDegree < 1) + { + throw new InvalidOperationException("Integration Degree must be greater or equal to 1. "); + } + + switch (integrationDegree) + { + case 1: + return new[] { GaussLegendrePoint1 }; + case 2: + return new[] { GaussLegendrePoint2A, GaussLegendrePoint2B }; + case 3: + return new[] { GaussLegendrePoint3A, GaussLegendrePoint3B, GaussLegendrePoint3C }; + case 4: + return new[] + { + GaussLegendrePoint4A, GaussLegendrePoint4B, GaussLegendrePoint4C, GaussLegendrePoint4D + }; + default: + throw new NotImplementedException("Integration Degree higher than 4 is not implemented yet. "); + } + } + + #endregion + } +} \ No newline at end of file diff --git a/IFiniteElementMaterial2D.cs b/IFiniteElementMaterial2D.cs new file mode 100644 index 000000000..20b5f98de --- /dev/null +++ b/IFiniteElementMaterial2D.cs @@ -0,0 +1,18 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using ISAAR.MSolve.Matrices.Interfaces; + +namespace ISAAR.MSolve.PreProcessor.Interfaces +{ + public interface IFiniteElementMaterial2D : IFiniteElementMaterial + { + double[] Stresses { get; } + IMatrix2D ConstitutiveMatrix { get; } + void UpdateMaterial(double[] strains); + void ClearState(); + void SaveState(); + void ClearStresses(); + } +} diff --git a/Quad4.cs b/Quad4.cs new file mode 100644 index 000000000..a81f8ecc1 --- /dev/null +++ b/Quad4.cs @@ -0,0 +1,619 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using ISAAR.MSolve.PreProcessor.Embedding; +using ISAAR.MSolve.PreProcessor.Interfaces; +using ISAAR.MSolve.Matrices.Interfaces; +using System.Runtime.InteropServices; +using ISAAR.MSolve.Matrices; +using ISAAR.MSolve.PreProcessor.Elements.SupportiveClasses; + +namespace ISAAR.MSolve.PreProcessor.Elements +{ + public class Quad4 : IStructuralFiniteElement, IEmbeddedHostElement + { + protected static double determinantTolerance = 0.00000001; + + protected static int iInt = 2; + protected static int iInt2 = iInt * iInt; + protected static int iInt3 = iInt2 * iInt; + protected readonly static DOFType[] nodalDOFTypes = new DOFType[] { DOFType.X, DOFType.Y}; //Panos - 2 DoF per node + protected readonly static DOFType[][] dofTypes = new DOFType[][] { nodalDOFTypes, nodalDOFTypes, + nodalDOFTypes,nodalDOFTypes}; //Panos - 4 nodes only + + //protected readonly IFiniteElementMaterial3D[] materialsAtGaussPoints; //Panos - What about 2D Material? Implemented + protected readonly IFiniteElementMaterial2D[] materialsAtGaussPoints; + + protected IFiniteElementDOFEnumerator dofEnumerator = new GenericDOFEnumerator(); + protected Quad4() + { + } + public Quad4(IFiniteElementMaterial2D material) + { + materialsAtGaussPoints = new IFiniteElementMaterial2D[iInt2]; + for (int i = 0; i < iInt2; i++) + materialsAtGaussPoints[i] = (IFiniteElementMaterial2D)material.Clone(); + } + public Quad4(IFiniteElementMaterial2D material, IFiniteElementDOFEnumerator dofEnumerator) + : this(material) + { + this.dofEnumerator = dofEnumerator; + } + public IFiniteElementDOFEnumerator DOFEnumerator + { + get { return dofEnumerator; } + set { dofEnumerator = value; } + } + //Panos - I dont think I need these + public double Density { get; set; } + public double RayleighAlpha { get; set; } + public double RayleighBeta { get; set; } + // + protected double[,] GetCoordinates(Element element) + { + //double[,] faXYZ = new double[dofTypes.Length, 3]; + double[,] faXY = new double[dofTypes.Length, 2]; + for (int i = 0; i < dofTypes.Length; i++) + { + faXY[i, 0] = element.Nodes[i].X; + faXY[i, 1] = element.Nodes[i].Y; + } + return faXY; + } + protected double[,] GetCoordinatesTranspose(Element element) + { + double[,] faXY = new double[2, dofTypes.Length]; + for (int i = 0; i < dofTypes.Length; i++) + { + faXY[0, i] = element.Nodes[i].X; + faXY[1, i] = element.Nodes[i].Y; + } + return faXY; + } + + #region IElementType Members + + public int ID + { + //Change Element ID to random 14? I am not sure why is this needed + get { return 14; } + } + + public ElementDimensions ElementDimensions + { + get { return ElementDimensions.TwoD; } + } + + public virtual IList> GetElementDOFTypes(Element element) + { + return dofTypes; + } + public IList GetNodesForMatrixAssembly(Element element) + { + return element.Nodes; + } + + private double[] CalcQ4Shape(double fXi, double fEta) + { + // QUAD4 ELEMENT + // + // Eta + // ^ + //#4(-1,1) | #3(1,1) + // ._____|_____. + // | | | + // | | | + // | |_____|____>Xi (or psi - the same) + // | | + // | | + // .___________. + // + //#1(-1,-1) #2(1,-1) + // + const double fSq050 = 0.50; //Panos - I need 0.25 so 0.50*0.50=0.25 + double fXiP = (1.0 + fXi) * fSq050; //(1+psi) - P means Plus + double fEtaP = (1.0 + fEta) * fSq050; //(1+eta) - P means Plus + double fXiM = (1.0 - fXi) * fSq050; //(1-psi) - M means Minus + double fEtaM = (1.0 - fEta) * fSq050; //(1-eta) - M means Minus + return new double[] + { + fXiM * fEtaM, //N1 = .25*(1-psi).*(1-eta); + fXiP * fEtaM, //N2 = .25*(1+psi).*(1-eta); + fXiP * fEtaP, //N3 = .25*(1+psi).*(1+eta); + fXiM * fEtaP, //N4 = .25*(1-psi).*(1+eta); + }; + } + + private double[] CalcQ4ShapeFunctionDerivatives(double fXi, double fEta) //Panos - Shape functions Derivatives + { + + const double fN025 = 0.25; //fSqC125 = 0.5; Panos - Should I Keep the same names? I dont multiply here + double fXiP = (1.0 + fXi) * fN025; //(1+psi) - P means Plus + double fEtaP = (1.0 + fEta) * fN025; //(1+eta) - P means Plus + double fXiM = (1.0 - fXi) * fN025; //(1-psi) - M means Minus + double fEtaM = (1.0 - fEta) * fN025; //(1-eta) - M means Minus + + double[] faDS = new double[8]; + //Panos - I put first with respect to psi(fXi here) and then eta + + faDS[0] = -fEtaM; //N1,psi = -.25*(1-eta); + faDS[1] = fEtaM; //N2,psi = .25*(1-eta); + faDS[2] = fEtaP; //N3,psi = .25*(1+eta); + faDS[3] = -fEtaP; //N4,psi = -.25*(1+eta); + + faDS[4] = -fXiM; //N1,eta = -.25*(1-psi); + faDS[5] = -fXiP; //N2,eta = -.25*(1+psi); + faDS[6] = fXiP; //N3,eta = .25*(1+psi); + faDS[7] = fXiM; //N4,eta = .25*(1-psi); + + return faDS; + } + + private double[,] CalcQ4J(double[,] faXY, double[] faDS, double fXi, double fEta) //Panos - Calculate the Jacobian matrix J + { + // Jacobian Matrix J + // [x1 y1] + // [ N1,psi N2,psi N3,psi N4,psi] [x2 y2] + //[J] = [ N1,eta N2,eta N3,eta N4,eta] [x3 y3] + // [x4 y4] + // + // [ ] multiplied by [ ] equals [ ] + // 2x4 4x2 2x2 + + // REMINDER faXY[i, 0] = element.Nodes[i].X; + // faXY[i, 1] = element.Nodes[i].Y; + + double[,] faJ = new double[2, 2]; + faJ[0, 0] = faDS[0] * faXY[0, 0] + faDS[1] * faXY[1, 0] + faDS[2] * faXY[2, 0] + faDS[3] * faXY[3, 0]; + faJ[0, 1] = faDS[0] * faXY[0, 1] + faDS[1] * faXY[1, 1] + faDS[2] * faXY[2, 1] + faDS[3] * faXY[3, 1]; + faJ[1, 0] = faDS[4] * faXY[0, 0] + faDS[5] * faXY[1, 0] + faDS[6] * faXY[2, 0] + faDS[7] * faXY[3, 0]; + faJ[1, 1] = faDS[4] * faXY[0, 1] + faDS[5] * faXY[1, 1] + faDS[6] * faXY[2, 1] + faDS[7] * faXY[3, 1]; + + return faJ; + } + private double CalcQ4JDetJ(double[,] faXY, double fXi, double fEta) //Calculate the determinant of the Jacobian Matrix + { + //|J| or det[j] = 1/8 [x1 x2 x3 x4] [ 0 1-eta eta-psi psi-1 ] [y1] + // [ eta-1 0 psi+1 -psi-eta] [y2] + // [psi-eta -psi-1 0 eta+1 ] [y3] + // [ 1-psi psi+eta -eta-1 0 ] [y4] + + const double fN0125 = 0.125; //=1/8 + + double fDetJ = (faXY[0, 0] * ( (1 - fEta) * faXY[1, 1] + (fEta - fXi) * faXY[2, 1] + (fXi - 1) * faXY[3, 1])) + + (faXY[1, 0] * ((fEta - 1) * faXY[0, 1] + (fXi + 1) * faXY[2, 1] + (-fXi - fEta) * faXY[3, 1])) + + (faXY[2, 0] * ((fXi - fEta) * faXY[0, 1] + (-fXi - 1) * faXY[1, 1] + (fEta + 1) * faXY[3, 1])) + + (faXY[3, 0] * ((1 - fXi) * faXY[0, 1] + (fXi + fEta) * faXY[1, 1] + (-fEta - 1) * faXY[2, 1])); + + fDetJ = fDetJ * fN0125; + + if (fDetJ < determinantTolerance) + throw new ArgumentException(String.Format("Jacobian determinant is negative or under tolerance ({0} < {1})." + + "Check the order of nodes or the element geometry.", fDetJ, determinantTolerance)); + return fDetJ; + } + private double[,] CalcQ4JInv(double fDetJ, double[,] faJ) //Calculate the full inverse of the Jacobian Matrix + { + // + //[J]^-1 = 1/det[J] [ J22 -J12] + // 1[-J21 J11] + + double fDetInv = 1.0 / fDetJ; + + double[,] faJInv = new double[2, 2]; + faJInv[0, 0] = (faJ[1, 1]) * fDetInv; //Panos - prosoxi stin arithmisi to J11 einai faJ[0,0]; + faJInv[0, 1] = (-faJ[0, 1]) * fDetInv; + faJInv[1, 0] = (-faJ[1, 0]) * fDetInv; + faJInv[1, 1] = (faJ[0, 0]) * fDetInv; + + return faJInv; + } + + private double[,] CalculateDeformationMatrix(double[,] faXY, double[] faDS, double fXi, double fEta, double fDetJ) //Panos - Calculate Deformation matrix B + { + + // [B] = 1/|J| [B1 B2 B3 B4] + // + // + // where [Bi] = [ a*Ni,psi - b*Ni,eta 0 ] + // + // [ 0 c*Ni,eta - d*Ni,psi ] + // + // [ c*Ni,eta - d*Ni,psi a*Ni,psi - b*Ni,eta] + // + // + // where a=1/4 {y1*(psi-1) + y2*(-1-psi) + y3*(1+psi) + y4*(1-psi)} + // b=1/4 {y1*(eta-1) + y2*(1-eta) + y3*(1+eta) + y4*(-1-eta)} + // c=1/4 {x1*(eta-1) + x2*(1-eta) + x3*(1+eta) + x4*(-1-eta)} + // d=1/4 {x1*(psi-1) + x2*(-1-psi) + x3*(1+psi) + x4*(1-psi)} + // + // + // REMINDER + // faDS[0] = N1,psi faDS[4] = N1,eta + // faDS[1] = N2,psi faDS[5] = N2,eta + // faDS[2] = N3,psi faDS[6] = N3,eta + // faDS[3] = N4,psi faDS[7] = N4,eta + + double fDetInv = 1.0 / fDetJ; //Panos - SOS, edo einai mono to klasma 1/det[j] kai oxi o full antistrofos pinakas; + double Aparameter; + double Bparameter; + double Cparameter; + double Dparameter; + + Aparameter = 0.25 * (faXY[0, 1] * (fXi - 1) + faXY[1, 1] * (-1 - fXi) + faXY[2, 1] * (1 + fXi) + faXY[3, 1] * (1 - fXi)); + Bparameter = 0.25 * (faXY[0, 1] * (fEta - 1) + faXY[1, 1] * (1 - fEta) + faXY[2, 1] * (1 + fEta) + faXY[3, 1] * (-1 - fEta)); + Cparameter = 0.25 * (faXY[0, 0] * (fEta - 1) + faXY[1, 0] * (1 - fEta) + faXY[2, 0] * (1 + fEta) + faXY[3, 0] * (-1 - fEta)); + Dparameter = 0.25 * (faXY[0, 0] * (fXi - 1) + faXY[1, 0] * (-1 - fXi) + faXY[2, 0] * (1 + fXi) + faXY[3, 0] * (1 - fXi)); + + double[,] Bmatrix = new double[3, 8]; + + Bmatrix[0, 0] = fDetInv * (Aparameter * faDS[0] - Bparameter * faDS[4]); + Bmatrix[1, 0] = 0; + Bmatrix[2, 0] = fDetInv * (Cparameter * faDS[4] - Dparameter * faDS[0]); + Bmatrix[0, 1] = 0; + Bmatrix[1, 1] = fDetInv * (Cparameter * faDS[4] - Dparameter * faDS[0]); + Bmatrix[2, 1] = fDetInv * (Aparameter * faDS[0] - Bparameter * faDS[4]); + + Bmatrix[0, 2] = fDetInv * (Aparameter * faDS[1] - Bparameter * faDS[5]); + Bmatrix[1, 2] = 0; + Bmatrix[2, 2] = fDetInv * (Cparameter * faDS[5] - Dparameter * faDS[1]); + Bmatrix[0, 3] = 0; + Bmatrix[1, 3] = fDetInv * (Cparameter * faDS[5] - Dparameter * faDS[1]); + Bmatrix[2, 3] = fDetInv * (Aparameter * faDS[1] - Bparameter * faDS[5]); + + Bmatrix[0, 4] = fDetInv * (Aparameter * faDS[2] - Bparameter * faDS[6]); + Bmatrix[1, 4] = 0; + Bmatrix[2, 4] = fDetInv * (Cparameter * faDS[6] - Dparameter*faDS[2]); + Bmatrix[0,5]= 0 ; + Bmatrix[1,5]= fDetInv * (Cparameter*faDS[6] - Dparameter*faDS[2]); + Bmatrix[2,5]= fDetInv * (Aparameter*faDS[2] - Bparameter*faDS[6]); + + Bmatrix[0,6]= fDetInv * (Aparameter*faDS[3] - Bparameter*faDS[7]); + Bmatrix[1,6]= 0 ; + Bmatrix[2,6]= fDetInv * (Cparameter*faDS[7] - Dparameter*faDS[3]); + Bmatrix[0,7]= 0 ; + Bmatrix[1,7]= fDetInv * (Cparameter*faDS[7] - Dparameter*faDS[3]); + Bmatrix[2,7]= fDetInv * (Aparameter*faDS[3] - Bparameter*faDS[7]); + + return Bmatrix; + //Bmatrix=fDetInv * Bmatrix; //Panos - den yparxei o operator fDetInv. * mhtroo??? - To kano analytika... + } + + private GaussLegendrePoint2D[] CalculateGaussMatrices(double[,] nodeCoordinates) + { + GaussLegendrePoint1D[] integrationPointsPerAxis = + GaussQuadrature.GetGaussLegendrePoints(iInt); + int totalSamplingPoints = (int)Math.Pow(integrationPointsPerAxis.Length, 2); + + GaussLegendrePoint2D[] integrationPoints = new GaussLegendrePoint2D[totalSamplingPoints]; + + int counter = -1; + foreach (GaussLegendrePoint1D pointXi in integrationPointsPerAxis) + { + foreach (GaussLegendrePoint1D pointEta in integrationPointsPerAxis) + { + counter += 1; + double xi = pointXi.Coordinate; + double eta = pointEta.Coordinate; + //double[] ShapeFunctions = this.CalcQ4Shape(xi,eta); //Panos - Uncommnent to check Shapefunctions + double[] faDS = this.CalcQ4ShapeFunctionDerivatives(xi, eta); + double [,] faJ = this.CalcQ4J( nodeCoordinates, faDS, xi, eta); + double fDetJ = this.CalcQ4JDetJ(nodeCoordinates, xi, eta); + //double[,] faJInv = this.CalcQ4JInv(fDetJ, faJ); //Panos - Uncommnent to check JInv + double [,] deformationMatrix = this.CalculateDeformationMatrix(nodeCoordinates,faDS, xi, eta, fDetJ); + double weightFactor = pointXi.WeightFactor * pointEta.WeightFactor * fDetJ; //Panos - we should also insert the thickness of the element t. Now it is assumed t=1; + integrationPoints[counter] = new GaussLegendrePoint2D(xi, eta, deformationMatrix, weightFactor); + } + } + return integrationPoints; + } + + public virtual IMatrix2D StiffnessMatrix(Element element) + { + double[,] coordinates = this.GetCoordinates(element); + GaussLegendrePoint2D[] integrationPoints = this.CalculateGaussMatrices(coordinates); + SymmetricMatrix2D stiffnessMatrix = new SymmetricMatrix2D(8); + int pointId = -1; + foreach (GaussLegendrePoint2D intPoint in integrationPoints) + { + pointId++; + IMatrix2D constitutiveMatrix = materialsAtGaussPoints[pointId].ConstitutiveMatrix; + double[,] b = intPoint.DeformationMatrix; + for (int i = 0; i < 8; i++) + { + double[] eb = new double[3]; + for (int iE = 0; iE < 3; iE++) + { + eb[iE] = (constitutiveMatrix[iE, 0] * b[0, i]) + (constitutiveMatrix[iE, 1] * b[1, i]) + + (constitutiveMatrix[iE, 2] * b[2, i]) ; + } + + for (int j = i; j < 8; j++) + { + double stiffness = (b[0, j] * eb[0]) + (b[1, j] * eb[1]) + (b[2, j] * eb[2]) ; + stiffnessMatrix[i, j] += stiffness * intPoint.WeightFactor; + } + } + } + + return stiffnessMatrix; + } + + public IMatrix2D CalculateConsistentMass(Element element) + { + double[,] coordinates = this.GetCoordinates(element); + GaussLegendrePoint2D[] integrationPoints = this.CalculateGaussMatrices(coordinates); + + SymmetricMatrix2D consistentMass = new SymmetricMatrix2D(24); + + foreach (GaussLegendrePoint2D intPoint in integrationPoints) + { + double[] shapeFunctionValues = this.CalcQ4Shape(intPoint.Xi, intPoint.Eta); + double weightDensity = intPoint.WeightFactor * this.Density; + for (int shapeFunctionI = 0; shapeFunctionI < shapeFunctionValues.Length; shapeFunctionI++) + { + for (int shapeFunctionJ = shapeFunctionI; shapeFunctionJ < shapeFunctionValues.Length; shapeFunctionJ++) + { + consistentMass[3 * shapeFunctionI, 3 * shapeFunctionJ] += shapeFunctionValues[shapeFunctionI] * + shapeFunctionValues[shapeFunctionJ] * + weightDensity; + } + + for (int shapeFunctionJ = shapeFunctionI; shapeFunctionJ < shapeFunctionValues.Length; shapeFunctionJ++) + { + consistentMass[(3 * shapeFunctionI) + 1, (3 * shapeFunctionJ) + 1] = + consistentMass[3 * shapeFunctionI, 3 * shapeFunctionJ]; + + consistentMass[(3 * shapeFunctionI) + 2, (3 * shapeFunctionJ) + 2] = + consistentMass[3 * shapeFunctionI, 3 * shapeFunctionJ]; + } + } + } + + return consistentMass; + } + + #endregion + + public virtual IMatrix2D MassMatrix(Element element) + { + return CalculateConsistentMass(element); + } + + public virtual IMatrix2D DampingMatrix(Element element) + { + var m = MassMatrix(element); + m.LinearCombination(new double[] { RayleighAlpha, RayleighBeta }, new IMatrix2D[] { MassMatrix(element), StiffnessMatrix(element) }); + return m; + //double[] faD = new double[300]; + //return new SymmetricMatrix2D(faD); + } + + public Tuple CalculateStresses(Element element, double[] localDisplacements, double[] localdDisplacements) + { + double[,] faXY = GetCoordinates(element); + double[,] faDS = new double[iInt3, 24]; + double[,] faS = new double[iInt3, 8]; + double[,,] faB = new double[iInt3, 24, 6]; + double[] faDetJ = new double[iInt3]; + double[,,] faJ = new double[iInt3, 3, 3]; + double[] faWeight = new double[iInt3]; + double[,] fadStrains = new double[iInt3, 6]; + double[,] faStrains = new double[iInt3, 6]; + // CalcQ4GaussMatrices(ref iInt, faXY, faWeight, faS, faDS, faJ, faDetJ, faB); + // CalcQ4Strains(ref iInt, faB, localDisplacements, faStrains); + // CalcQ4Strains(ref iInt, faB, localdDisplacements, fadStrains); + + double[] dStrains = new double[6]; + double[] strains = new double[6]; + for (int i = 0; i < materialsAtGaussPoints.Length; i++) + { + for (int j = 0; j < 6; j++) dStrains[j] = fadStrains[i, j]; + for (int j = 0; j < 6; j++) strains[j] = faStrains[i, j]; + materialsAtGaussPoints[i].UpdateMaterial(dStrains); + } + + return new Tuple(strains, materialsAtGaussPoints[materialsAtGaussPoints.Length - 1].Stresses); + } + + public double[] CalculateForcesForLogging(Element element, double[] localDisplacements) + { + return CalculateForces(element, localDisplacements, new double[localDisplacements.Length]); + } + + public double[] CalculateForces(Element element, double[] localTotalDisplacements, double[] localdDisplacements) + { + //Vector d = new Vector(localdDisplacements.Length); + //for (int i = 0; i < localdDisplacements.Length; i++) + // //d[i] = localdDisplacements[i] + localTotalDisplacements[i]; + // d[i] = localTotalDisplacements[i]; + //double[] faForces = new double[24]; + //StiffnessMatrix(element).Multiply(d, faForces); + + double[,] faStresses = new double[iInt3, 6]; + for (int i = 0; i < materialsAtGaussPoints.Length; i++) + for (int j = 0; j < 6; j++) faStresses[i, j] = materialsAtGaussPoints[i].Stresses[j]; + + double[,] faXYZ = GetCoordinates(element); + double[,] faDS = new double[iInt3, 24]; + double[,] faS = new double[iInt3, 8]; + double[,,] faB = new double[iInt3, 24, 6]; + double[] faDetJ = new double[iInt3]; + double[,,] faJ = new double[iInt3, 3, 3]; + double[] faWeight = new double[iInt3]; + double[] faForces = new double[24]; + // CalcQ4GaussMatrices(ref iInt, faXYZ, faWeight, faS, faDS, faJ, faDetJ, faB); + // CalcQ4Forces(ref iInt, faB, faWeight, faStresses, faForces); + + return faForces; + } + + public double[] CalculateAccelerationForces(Element element, IList loads) + { + Vector accelerations = new Vector(24); + IMatrix2D massMatrix = MassMatrix(element); + + foreach (MassAccelerationLoad load in loads) + { + int index = 0; + foreach (DOFType[] nodalDOFTypes in dofTypes) + foreach (DOFType dofType in nodalDOFTypes) + { + if (dofType == load.DOF) accelerations[index] += load.Amount; + index++; + } + } + double[] forces = new double[24]; + massMatrix.Multiply(accelerations, forces); + return forces; + } + + public void ClearMaterialState() + { + foreach (IFiniteElementMaterial2D m in materialsAtGaussPoints) m.ClearState(); + } + + public void SaveMaterialState() + { + foreach (IFiniteElementMaterial2D m in materialsAtGaussPoints) m.SaveState(); + } + + public void ClearMaterialStresses() + { + foreach (IFiniteElementMaterial2D m in materialsAtGaussPoints) m.ClearStresses(); + } + + #region IStructuralFiniteElement Members + + public bool MaterialModified + { + get + { + foreach (IFiniteElementMaterial2D material in materialsAtGaussPoints) + if (material.Modified) return true; + return false; + } + } + + public void ResetMaterialModified() + { + foreach (IFiniteElementMaterial2D material in materialsAtGaussPoints) material.ResetModified(); + } + + #endregion + + #region IEmbeddedHostElement Members + + public EmbeddedNode BuildHostElementEmbeddedNode(Element element, Node node, IEmbeddedDOFInHostTransformationVector transformationVector) + { + var points = GetNaturalCoordinates(element, node); + if (points.Length == 0) return null; + + element.EmbeddedNodes.Add(node); + var embeddedNode = new EmbeddedNode(node, element, transformationVector.GetDependentDOFTypes); + for (int i = 0; i < points.Length; i++) + embeddedNode.Coordinates.Add(points[i]); + return embeddedNode; + } + + public double[] GetShapeFunctionsForNode(Element element, EmbeddedNode node) + { + double[,] elementCoordinates = GetCoordinatesTranspose(element); + var shapeFunctions = CalcQ4Shape(node.Coordinates[0], node.Coordinates[1]); + var ShapeFunctionDerivatives = CalcQ4ShapeFunctionDerivatives(node.Coordinates[0], node.Coordinates[1]); + var jacobian = CalcQ4J(elementCoordinates, ShapeFunctionDerivatives, node.Coordinates[0], node.Coordinates[1]); + + return new double[] + { + shapeFunctions[0], shapeFunctions[1], shapeFunctions[2], shapeFunctions[3], + ShapeFunctionDerivatives[0], ShapeFunctionDerivatives[1], ShapeFunctionDerivatives[2], ShapeFunctionDerivatives[3], ShapeFunctionDerivatives[4], ShapeFunctionDerivatives[5], ShapeFunctionDerivatives[6], ShapeFunctionDerivatives[7], + jacobian[0, 0], jacobian[0, 1], jacobian[1, 0], jacobian[1, 1] + }; + + //double[] coords = GetNaturalCoordinates(element, node); + //return CalcH8Shape(coords[0], coords[1], coords[2]); + + //double fXiP = (1.0 + coords[0]) * 0.5; + //double fEtaP = (1.0 + coords[1]) * 0.5; + //double fZetaP = (1.0 + coords[2]) * 0.5; + //double fXiM = (1.0 - coords[0]) * 0.5; + //double fEtaM = (1.0 - coords[1]) * 0.5; + //double fZetaM = (1.0 - coords[2]) * 0.5; + + //return new double[] { fXiM * fEtaM * fZetaM, + // fXiP * fEtaM * fZetaM, + // fXiP * fEtaP * fZetaM, + // fXiM * fEtaP * fZetaM, + // fXiM * fEtaM * fZetaP, + // fXiP * fEtaM * fZetaP, + // fXiP * fEtaP * fZetaP, + // fXiM * fEtaP * fZetaP }; + } + + private double[] GetNaturalCoordinates(Element element, Node node) + { + double[] mins = new double[] { element.Nodes[0].X, element.Nodes[0].Y}; + double[] maxes = new double[] { element.Nodes[0].X, element.Nodes[0].Y}; + 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]; + 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]; + } + //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]; + 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 = CalcQ4Shape(naturalCoordinates[0], naturalCoordinates[1]); + double[] coordinateDifferences = new double[] { 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[0] = node.X - coordinateDifferences[0]; + coordinateDifferences[1] = node.Y - coordinateDifferences[1]; + + double[,] faXY = GetCoordinatesTranspose(element); + double[] ShapeFunctionDerivatives = CalcQ4ShapeFunctionDerivatives(naturalCoordinates[0], naturalCoordinates[1]); + //SOS PANOS - THIS IS WRONG, I SHOULD IMPLEMENT INVERSE JACOBIAN + var inverseJacobian = CalcQ4J(faXY, ShapeFunctionDerivatives,naturalCoordinates[0], naturalCoordinates[1]); + //SOS PANOS END + 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; + } + + #endregion + } +} + diff --git a/Quad4Example.cs b/Quad4Example.cs new file mode 100644 index 000000000..a344b077d --- /dev/null +++ b/Quad4Example.cs @@ -0,0 +1,125 @@ +using System; +using System.Collections.Generic; +using ISAAR.MSolve.PreProcessor; +using ISAAR.MSolve.Solvers.Skyline; +using ISAAR.MSolve.Problems; +using ISAAR.MSolve.Analyzers; +using ISAAR.MSolve.Logging; +using ISAAR.MSolve.PreProcessor.Materials; +using ISAAR.MSolve.PreProcessor.Elements; +using ISAAR.MSolve.Matrices; + +namespace ISAAR.MSolve.Problems.Quad4Example +{ + class Quad4Example + { + private static IList CreateNodes() + { + IList nodes = new List(); + Node node1 = new Node { ID = 1, X = 0.0, Y = 0.0, Z = 0.0 }; + Node node2 = new Node { ID = 2, X = 1.0, Y = 0.0, Z = 0.0 }; + Node node3 = new Node { ID = 3, X = 1.0, Y = 1.0, Z = 0.0 }; + Node node4 = new Node { ID = 4, X = 0.0, Y = 1.0, Z = 0.0 }; + + nodes.Add(node1); + nodes.Add(node2); + nodes.Add(node3); + nodes.Add(node4); + + return nodes; + } + + static void Main(string[] args) + { + VectorExtensions.AssignTotalAffinityCount(); + double youngModulus = 3.0e07; + double poissonRatio = 0.3; + double nodalLoad = 1000; + + // Create a new elastic 2D material + ElasticMaterial2D material = new ElasticMaterial2D() + { + YoungModulus = youngModulus, + PoissonRatio = poissonRatio + }; + + // Node creation + IList nodes = CreateNodes(); + + // Model creation + Model model = new Model(); + + // Add a single subdomain to the model + model.SubdomainsDictionary.Add(1, new Subdomain() { ID = 1 }); + + // Add nodes to the nodes dictonary of the model + for (int i = 0; i < nodes.Count; ++i) + { + model.NodesDictionary.Add(i + 1, nodes[i]); + } + + // Constrain left nodes of the model + model.NodesDictionary[1].Constraints.Add(DOFType.X); + model.NodesDictionary[1].Constraints.Add(DOFType.Y); + model.NodesDictionary[1].Constraints.Add(DOFType.Z); //Panos - I leave Z here + model.NodesDictionary[4].Constraints.Add(DOFType.X); + model.NodesDictionary[4].Constraints.Add(DOFType.Y); + model.NodesDictionary[4].Constraints.Add(DOFType.Z); //Panos - Also here + + + // Create a new Quad4 element + var element = new Element() + { + ID = 1, + ElementType = new Quad4(material) + }; + + // Add nodes to the created element + element.AddNode(model.NodesDictionary[1]); + element.AddNode(model.NodesDictionary[2]); + element.AddNode(model.NodesDictionary[3]); + element.AddNode(model.NodesDictionary[4]); + + + // Add Quad4 element to the element and subdomains dictionary of the model + model.ElementsDictionary.Add(element.ID, element); + model.SubdomainsDictionary[1].ElementsDictionary.Add(element.ID, element); + + // Add nodal load values to node 3 + model.Loads.Add(new Load() { Amount = nodalLoad, Node = model.NodesDictionary[3], DOF = DOFType.X }); + + // Needed in order to make all the required data structures + model.ConnectDataStructures(); + + // Choose linear equation system solver + SolverSkyline solver = new SolverSkyline(model); + + // Choose the provider of the problem -> here a structural problem + ProblemStructural provider = new ProblemStructural(model, solver.SubdomainsDictionary); + + // Choose parent and child analyzers -> Parent: Static, Child: Linear + Analyzers.LinearAnalyzer childAnalyzer = new LinearAnalyzer(solver, solver.SubdomainsDictionary); + StaticAnalyzer parentAnalyzer = new StaticAnalyzer(provider, childAnalyzer, solver.SubdomainsDictionary); + + // Choose dof types X, Y, Z to log for node 2 + childAnalyzer.LogFactories[1] = new LinearAnalyzerLogFactory(new int[] { + model.NodalDOFsDictionary[2][DOFType.X], + model.NodalDOFsDictionary[2][DOFType.Y], + // Choose dof types X, Y, Z to log for node 3 + model.NodalDOFsDictionary[3][DOFType.X], + model.NodalDOFsDictionary[3][DOFType.Y]}); + + // Analyze the problem + parentAnalyzer.BuildMatrices(); + parentAnalyzer.Initialize(); + parentAnalyzer.Solve(); + + // Write results to console + Console.WriteLine("Writing results for node 2 and node 3"); + Console.WriteLine("Dof and Values for Displacement X, Y"); + Console.WriteLine(childAnalyzer.Logs[1][0]); + + + } + } +} \ No newline at end of file