From a195eb894bf19c25ae461922e3c7c2622ee3f504 Mon Sep 17 00:00:00 2001
From: Cecilie Brandt-Olsen <c.brandtolsen@gmail.com>
Date: Sun, 11 Feb 2018 22:39:41 +0100
Subject: [PATCH 1/4] Added pressure load goal with output

---
 .../K2Engineering/DataTypes/PressureData.cs   |  42 +++
 .../K2Engineering/K2Engineering.csproj        |   9 +
 K2Engineering/K2Engineering/Pressure.cs       | 259 ++++++++++++++++++
 K2Engineering/K2Engineering/PressureOutput.cs | 101 +++++++
 .../Properties/Resources.Designer.cs          |  20 ++
 .../K2Engineering/Properties/Resources.resx   |   6 +
 .../K2Engineering/Resources/Pressure.png      | Bin 0 -> 1260 bytes
 .../Resources/PressureOutput.png              | Bin 0 -> 1298 bytes
 8 files changed, 437 insertions(+)
 create mode 100644 K2Engineering/K2Engineering/DataTypes/PressureData.cs
 create mode 100644 K2Engineering/K2Engineering/Pressure.cs
 create mode 100644 K2Engineering/K2Engineering/PressureOutput.cs
 create mode 100644 K2Engineering/K2Engineering/Resources/Pressure.png
 create mode 100644 K2Engineering/K2Engineering/Resources/PressureOutput.png

diff --git a/K2Engineering/K2Engineering/DataTypes/PressureData.cs b/K2Engineering/K2Engineering/DataTypes/PressureData.cs
new file mode 100644
index 0000000..56047d7
--- /dev/null
+++ b/K2Engineering/K2Engineering/DataTypes/PressureData.cs
@@ -0,0 +1,42 @@
+using Rhino.Geometry;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace K2Engineering.DataTypes
+{
+    public class PressureData
+    {
+        public Point3d[] Locations;
+        public Vector3d[] Loads;
+        public double PresStart;
+        public double PresEnd;
+        public double VolStart;
+        public double VolEnd;
+        public int MolStart;
+        public int MolEnd;
+
+        public PressureData() { }
+
+        public PressureData(Point3d[] locations, Vector3d[] loads, double presStart, double presEnd, double volStart, double volEnd, int molStart, int molEnd)
+        {
+            Locations = new Point3d[locations.Length];
+            Loads = new Vector3d[loads.Length];
+
+            for(int i=0; i<locations.Length; i++)
+            {
+                Locations[i] = locations[i];
+                Loads[i] = loads[i];
+            }
+
+            PresStart = presStart;
+            PresEnd = presEnd;
+            VolStart = volStart;
+            VolEnd = volEnd;
+            MolStart = molStart;
+            MolEnd = molEnd;
+        }
+    }
+}
diff --git a/K2Engineering/K2Engineering/K2Engineering.csproj b/K2Engineering/K2Engineering/K2Engineering.csproj
index 6aa5043..0d82586 100644
--- a/K2Engineering/K2Engineering/K2Engineering.csproj
+++ b/K2Engineering/K2Engineering/K2Engineering.csproj
@@ -78,7 +78,10 @@
     <Compile Include="BucklingAnalysis.cs" />
     <Compile Include="Cable.cs" />
     <Compile Include="BarOutput.cs" />
+    <Compile Include="DataTypes\PressureData.cs" />
     <Compile Include="LoadOutput.cs" />
+    <Compile Include="Pressure.cs" />
+    <Compile Include="PressureOutput.cs" />
     <Compile Include="RodOutput.cs" />
     <Compile Include="SupportOutput.cs" />
     <Compile Include="CircularCS.cs" />
@@ -198,6 +201,12 @@
   <ItemGroup>
     <None Include="Resources\Pretension.png" />
   </ItemGroup>
+  <ItemGroup>
+    <None Include="Resources\Pressure.png" />
+  </ItemGroup>
+  <ItemGroup>
+    <None Include="Resources\PressureOutput.png" />
+  </ItemGroup>
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
        Other similar extension points exist, see Microsoft.Common.targets.
diff --git a/K2Engineering/K2Engineering/Pressure.cs b/K2Engineering/K2Engineering/Pressure.cs
new file mode 100644
index 0000000..2f06f99
--- /dev/null
+++ b/K2Engineering/K2Engineering/Pressure.cs
@@ -0,0 +1,259 @@
+using System;
+using System.Collections.Generic;
+using Grasshopper.Kernel;
+using Rhino.Geometry;
+using Plankton;
+using KangarooSolver;
+
+namespace K2Engineering
+{
+    public class Pressure : GH_Component
+    {
+        /// <summary>
+        /// Initializes a new instance of the Pressure class.
+        /// </summary>
+        public Pressure()
+          : base("Pressure", "Pressure",
+              "A pressure load goal following the principle of the Ideal Gas Law and with forces perpendicular to the mesh faces",
+              "K2Eng", "2 Load")
+        {
+        }
+
+        /// <summary>
+        /// Registers all the input parameters for this component.
+        /// </summary>
+        protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager)
+        {
+            pManager.AddGenericParameter("PMesh", "pMesh", "A plankton mesh in [m]. The mesh has to be triangulated", GH_ParamAccess.item);
+            pManager.AddNumberParameter("Pressure", "p", "The pressure in [kN/m2]", GH_ParamAccess.item);
+            pManager.AddBooleanParameter("PressureOption", "opt", "If true, the pressure is constant. If false, the amount of molecules is constant but the pressure varies according to the change in volume", GH_ParamAccess.item, true);
+        }
+
+        /// <summary>
+        /// Registers all the output parameters for this component.
+        /// </summary>
+        protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager)
+        {
+            pManager.AddGenericParameter("PressureGoal", "P", "The pressure goal", GH_ParamAccess.item);
+        }
+
+        /// <summary>
+        /// This is the method that actually does the work.
+        /// </summary>
+        /// <param name="DA">The DA object is used to retrieve from inputs and store in outputs.</param>
+        protected override void SolveInstance(IGH_DataAccess DA)
+        {
+            //Input
+            PlanktonMesh pMesh = new PlanktonMesh();
+            DA.GetData(0, ref pMesh);
+
+            double pressure = 0.0;
+            DA.GetData(1, ref pressure);
+            pressure *= 1e3;   //From [kN/m2] to [N/m2]
+
+            bool constantP = true;
+            DA.GetData(2, ref constantP);
+
+
+            //Calculate
+            GoalObject pressureGoal = new PressureGoal(pMesh, pressure, constantP);
+
+
+            //Output
+            DA.SetData(0, pressureGoal);
+
+        }
+
+
+        public class PressureGoal : GoalObject
+        {
+            //Global properties
+            PlanktonMesh pMesh;
+
+            double pres0;
+            double pres1;
+
+            double vol0;
+            double vol1;
+
+            double mol0;
+            double mol1;
+
+            double gasConst;
+            double temp;
+
+            bool opt;
+
+            public PressureGoal(PlanktonMesh mesh, double p, bool option)
+            {
+                pMesh = (PlanktonMesh) mesh;
+
+                //Pressure
+                pres0 = p;      // Unit: [N/m2] / [Pa]
+                pres1 = pres0;
+
+                //Calculate starting volume
+                vol0 = 0.0;         // Unit: [m3]
+
+                for (int j = 0; j < pMesh.Faces.Count; j++)
+                {
+                    int[] faceVertices = pMesh.Faces.GetFaceVertices(j);
+
+                    Point3d pA = new Point3d(pMesh.Vertices[faceVertices[0]].X, pMesh.Vertices[faceVertices[0]].Y, pMesh.Vertices[faceVertices[0]].Z);
+                    Point3d pB = new Point3d(pMesh.Vertices[faceVertices[1]].X, pMesh.Vertices[faceVertices[1]].Y, pMesh.Vertices[faceVertices[1]].Z);
+                    Point3d pC = new Point3d(pMesh.Vertices[faceVertices[2]].X, pMesh.Vertices[faceVertices[2]].Y, pMesh.Vertices[faceVertices[2]].Z);
+
+                    Point3d pD = new Point3d(0, 0, 0);
+
+                    Vector3d AB = pB - pA;
+                    Vector3d AC = pC - pA;
+                    Vector3d AD = pD - pA;
+
+                    Vector3d cross = Vector3d.CrossProduct(AB, -AC);
+                    double dot = Vector3d.Multiply(cross, AD);
+
+                    double v = (1.0 / 6.0) * dot;
+
+                    vol0 += v;
+                }
+
+                vol1 = vol0;
+
+                //Gas constant and temperature
+                gasConst = 8.314472;            // Unit: [J/(K*mol)]
+                temp = 273 + 20;                // Unit: [Kelvins] Set to 20 degrees by default
+
+                //Molecules
+                mol0 = (pres0 * vol0) / (gasConst * temp);      // Unit: [moles]
+                mol1 = mol0;
+
+                opt = option;
+
+
+                PPos = new Point3d[pMesh.Vertices.Count];
+                Move = new Vector3d[pMesh.Vertices.Count];
+                Weighting = new double[pMesh.Vertices.Count];
+
+                for (int i = 0; i < pMesh.Vertices.Count; i++)
+                {
+                    PPos[i] = new Point3d(pMesh.Vertices[i].X, pMesh.Vertices[i].Y, pMesh.Vertices[i].Z);
+                    Move[i] = new Vector3d(0, 0, 0);
+                    Weighting[i] = 1.0;
+                }
+
+            }
+
+
+            public override void Calculate(List<KangarooSolver.Particle> p)
+            {
+                //Lists of vertex normals and areas
+                Vector3d[] normals = new Vector3d[pMesh.Vertices.Count];
+                double[] vertexAreas = new double[pMesh.Vertices.Count];
+
+                for (int j = 0; j < pMesh.Vertices.Count; j++)
+                {
+                    normals[j] = new Vector3d(0, 0, 0);
+                    vertexAreas[j] = 0.0;
+                }
+
+                vol1 = 0.0;
+
+
+                //Run through all mesh faces to calculate vertex normals, areas and current volume
+                for (int i = 0; i < pMesh.Faces.Count; i++)
+                {
+
+                    //Get face vertices
+                    int[] faceVertices = pMesh.Faces.GetFaceVertices(i);
+
+                    Point3d PA = p[PIndex[faceVertices[0]]].Position;
+                    Point3d PB = p[PIndex[faceVertices[1]]].Position;
+                    Point3d PC = p[PIndex[faceVertices[2]]].Position;
+
+                    Vector3d AB = PB - PA;
+                    Vector3d BC = PC - PB;
+                    Vector3d CA = PA - PC;
+
+                    Vector3d normal = Vector3d.CrossProduct(AB, BC);
+                    double nodalArea = normal.Length / 6.0;
+                    normal.Unitize();
+
+                    normals[faceVertices[0]] += normal;
+                    normals[faceVertices[1]] += normal;
+                    normals[faceVertices[2]] += normal;
+
+                    vertexAreas[faceVertices[0]] += nodalArea;
+                    vertexAreas[faceVertices[1]] += nodalArea;
+                    vertexAreas[faceVertices[2]] += nodalArea;
+
+
+                    //Current volume
+                    Point3d PD = new Point3d(0, 0, 0);
+                    Vector3d AD = PD - PA;
+                    Vector3d cross = Vector3d.CrossProduct(AB, CA);
+                    double dot = Vector3d.Multiply(cross, AD);
+                    vol1 += (1.0 / 6.0) * dot;
+                }
+
+                //constant molecules -> pressure adjusts
+                if (opt == false)
+                {
+                    pres1 = (mol0 * gasConst * temp) / vol1;
+                }
+                //constant pressure -> molecules adjust
+                else
+                {
+                    mol1 = (pres0 * vol1) / (gasConst * temp);
+                }
+
+
+                //Calculate pressure vector acting in each vertex in [N]
+                for (int k = 0; k < pMesh.Vertices.Count; k++)
+                {
+                    Vector3d n = new Vector3d(normals[k]);
+                    n.Unitize();
+                    Move[k] = pres1 * vertexAreas[k] * n;
+                }
+
+            }
+
+
+            public override object Output(List<KangarooSolver.Particle> p)
+            {
+                Point3d[] vertices = new Point3d[pMesh.Vertices.Count];
+                Vector3d[] pForces = new Vector3d[pMesh.Vertices.Count];
+
+                for (int i = 0; i < pMesh.Vertices.Count; i++)
+                {
+                    vertices[i] = p[PIndex[i]].Position;
+                    pForces[i] = Move[i] * Weighting[i] * 1e-3;       // Unit: [kN]
+                }
+
+                //Create pressure data object to store output information
+                DataTypes.PressureData presData = new DataTypes.PressureData(vertices, pForces, Math.Round(pres0*1e-3,3), Math.Round(pres1*1e-3,3), Math.Round(vol0,3), Math.Round(vol1,3), Convert.ToInt32(mol0), Convert.ToInt32(mol1));
+                return presData;
+            }
+
+        }
+
+
+        /// <summary>
+        /// Provides an Icon for the component.
+        /// </summary>
+        protected override System.Drawing.Bitmap Icon
+        {
+            get
+            {
+                return Properties.Resources.Pressure;
+            }
+        }
+
+        /// <summary>
+        /// Gets the unique ID for this component. Do not change this ID after release.
+        /// </summary>
+        public override Guid ComponentGuid
+        {
+            get { return new Guid("dfe7573e-5d67-481b-8226-15fe9b07d3fd"); }
+        }
+    }
+}
\ No newline at end of file
diff --git a/K2Engineering/K2Engineering/PressureOutput.cs b/K2Engineering/K2Engineering/PressureOutput.cs
new file mode 100644
index 0000000..eeb084f
--- /dev/null
+++ b/K2Engineering/K2Engineering/PressureOutput.cs
@@ -0,0 +1,101 @@
+using System;
+using System.Collections.Generic;
+using Grasshopper.Kernel;
+using Rhino.Geometry;
+
+namespace K2Engineering
+{
+    public class PressureOutput : GH_Component
+    {
+        /// <summary>
+        /// Initializes a new instance of the PressureOutput class.
+        /// </summary>
+        public PressureOutput()
+          : base("PressureOutput", "PressureOutput",
+              "Extract the output of the pressure goal",
+              "K2Eng", "6 Utility")
+        {
+        }
+
+        /// <summary>
+        /// Registers all the input parameters for this component.
+        /// </summary>
+        protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager)
+        {
+            pManager.AddGenericParameter("PressureData", "PD", "The PressureData from the output of the Pressure goal", GH_ParamAccess.item);
+        }
+
+        /// <summary>
+        /// Registers all the output parameters for this component.
+        /// </summary>
+        protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager)
+        {
+            pManager.AddPointParameter("Positions", "pos", "The positions of the load", GH_ParamAccess.list);
+            pManager.AddVectorParameter("Loads", "loads", "The nodal loads in [kN]", GH_ParamAccess.list);
+            pManager.AddNumberParameter("PressureStart", "presStart", "The pressure at the start in [kN/m2]", GH_ParamAccess.item);
+            pManager.AddNumberParameter("PressureEnd", "presEnd", "The pressure at the end in [kN/m2]", GH_ParamAccess.item);
+            pManager.AddNumberParameter("VolumeStart", "volStart", "The volume at the start in [m3]", GH_ParamAccess.item);
+            pManager.AddNumberParameter("VolumeEnd", "volEnd", "The volume at the end in [m3]", GH_ParamAccess.item);
+            pManager.AddIntegerParameter("MoleculesStart", "molStart", "The amount of molecules at the start in [mole]", GH_ParamAccess.item);
+            pManager.AddIntegerParameter("MoleculesEnd", "molEnd", "The amount of molecules at the end in [mole]", GH_ParamAccess.item);
+        }
+
+        /// <summary>
+        /// This is the method that actually does the work.
+        /// </summary>
+        /// <param name="DA">The DA object is used to retrieve from inputs and store in outputs.</param>
+        protected override void SolveInstance(IGH_DataAccess DA)
+        {
+            //Input
+            DataTypes.PressureData presData = new DataTypes.PressureData();
+            DA.GetData(0, ref presData);
+
+            //Extract properties
+            Point3d[] pts = new Point3d[presData.Locations.Length];
+            Vector3d[] forces = new Vector3d[presData.Loads.Length];
+
+            for(int i=0; i<presData.Locations.Length; i++)
+            {
+                pts[i] = presData.Locations[i];
+                forces[i] = presData.Loads[i];
+            }
+
+            double p0 = presData.PresStart;
+            double p1 = presData.PresEnd;
+            double v0 = presData.VolStart;
+            double v1 = presData.VolEnd;
+            int mol0 = presData.MolStart;
+            int mol1 = presData.MolEnd;
+
+
+            //Output
+            DA.SetDataList(0, pts);
+            DA.SetDataList(1, forces);
+            DA.SetData(2, p0);
+            DA.SetData(3, p1);
+            DA.SetData(4, v0);
+            DA.SetData(5, v1);
+            DA.SetData(6, mol0);
+            DA.SetData(7, mol1);
+        }
+
+        /// <summary>
+        /// Provides an Icon for the component.
+        /// </summary>
+        protected override System.Drawing.Bitmap Icon
+        {
+            get
+            {
+                return Properties.Resources.PressureOutput;
+            }
+        }
+
+        /// <summary>
+        /// Gets the unique ID for this component. Do not change this ID after release.
+        /// </summary>
+        public override Guid ComponentGuid
+        {
+            get { return new Guid("a2cfee5b-fb59-4092-b25b-dea8203ea8a5"); }
+        }
+    }
+}
\ No newline at end of file
diff --git a/K2Engineering/K2Engineering/Properties/Resources.Designer.cs b/K2Engineering/K2Engineering/Properties/Resources.Designer.cs
index 25a8748..5842e78 100644
--- a/K2Engineering/K2Engineering/Properties/Resources.Designer.cs
+++ b/K2Engineering/K2Engineering/Properties/Resources.Designer.cs
@@ -220,6 +220,26 @@ internal static System.Drawing.Bitmap MeshVertexArea {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized resource of type System.Drawing.Bitmap.
+        /// </summary>
+        internal static System.Drawing.Bitmap Pressure {
+            get {
+                object obj = ResourceManager.GetObject("Pressure", resourceCulture);
+                return ((System.Drawing.Bitmap)(obj));
+            }
+        }
+        
+        /// <summary>
+        ///   Looks up a localized resource of type System.Drawing.Bitmap.
+        /// </summary>
+        internal static System.Drawing.Bitmap PressureOutput {
+            get {
+                object obj = ResourceManager.GetObject("PressureOutput", resourceCulture);
+                return ((System.Drawing.Bitmap)(obj));
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized resource of type System.Drawing.Bitmap.
         /// </summary>
diff --git a/K2Engineering/K2Engineering/Properties/Resources.resx b/K2Engineering/K2Engineering/Properties/Resources.resx
index 8957ac7..1552697 100644
--- a/K2Engineering/K2Engineering/Properties/Resources.resx
+++ b/K2Engineering/K2Engineering/Properties/Resources.resx
@@ -199,4 +199,10 @@
   <data name="Pretension" type="System.Resources.ResXFileRef, System.Windows.Forms">
     <value>..\Resources\Pretension.png;System.Drawing.Bitmap, System.Drawing, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a</value>
   </data>
+  <data name="Pressure" type="System.Resources.ResXFileRef, System.Windows.Forms">
+    <value>..\Resources\Pressure.png;System.Drawing.Bitmap, System.Drawing, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a</value>
+  </data>
+  <data name="PressureOutput" type="System.Resources.ResXFileRef, System.Windows.Forms">
+    <value>..\Resources\PressureOutput.png;System.Drawing.Bitmap, System.Drawing, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a</value>
+  </data>
 </root>
\ No newline at end of file
diff --git a/K2Engineering/K2Engineering/Resources/Pressure.png b/K2Engineering/K2Engineering/Resources/Pressure.png
new file mode 100644
index 0000000000000000000000000000000000000000..9a81a325977bc9e78f69d4b04ad0d25420890097
GIT binary patch
literal 1260
zcmV<I1QYv-P)<h;3K|Lk000e1NJLTq000;O000;W1^@s6;CDUv00004b3#c}2nYxW
zd<bNS00009a7bBm000ie000ie0hKEb8vp<R8FWQhbW?9;ba!ELWdL_~cP?peYja~^
zaAhuUa%Y?FJQ@H11Zhb`K~zYIwUu3LlXVovf9Gx25r#k;5!Q{ltZ%!_;gEO%a<!on
zvyh?4=7^CiH5$D_Fb0VZM2TO41Q#`uEJi-G8}h9zNJ0W=)Ql(zgq3#=0&N(K3Zs;*
zd;1<QSW8!6#0$>d|D5MM&vVZI{eb%%GFGv=y4o`|G*o98MuVz8s;X0gLJ^6HNRQ(<
zdw~6xWexma0V5a;E(P8LPKwBG;P>lZ|K-BDFMEwWEskLrGgP$!s1=dcL?W>v7K<hC
zDNtTsKH2SdZyh<Vuq^9DrsnU^FI1@$;rZBNIR#KvRaI=;b_+09RiCvit0$Mwm3t^F
zEAx6hp2H&2X<625@p$}VmVBE|$Hkya8FU0!BE0%IK!1OKA|8+L_4$1LBC^%z^L57K
z@k_T&P+ME;jz**VMdV=O!|sy^Z=?JcYn0^~%#^TB@CYlTjj}Q;+^+Kx8gM8ndN;HQ
zST7=TBauiyfRO>7KY#vJ;A$iiStpPNensfL(*=!+aU21P#Nf*5)2H{T>So7rKFX90
zfSQ_`A`y8<L{<Q(&EI144NB7Ap>N!c+rLN0GeoP)W-Si;00svKKNXRug27-_Rv?*7
z)&oZ)k;p#)u_baY+ANV+qqN;kB$xbLv-Cupq+?_gH-N8H^`)$Uh%8jq9d~jWYm_$N
zMoGIaxCA#uXhT&PWMkT4&ADdj$?r+0(>p|DVKxyA29Je8p^E9>YoMfEm)>z7s-qt1
z@-`ip4=LHMb^fqEJ;HezaJfe&C{8whIh1zM3GDWF=<3`?vDp$8iR9)IJ_`Kp+A_)m
zAW&EGgL+vN>Z`1*oa(-5P95)G9m)8td56;u#|ly$%x#>KAQ8=T+uuOKzn3+m1RNew
zhth&!4ghOYsZ<&O(=`9AsHg}6OxvdA`J<P&={WNzEzHgN!&;e-oBE?ZR2^0iKqwT7
z14ibOj~j+D2jIdIIg?MGnn3Hd1b%?BE$RgjgPH!YuE@pqEs~2Tnj{qn1ZJq}8Ne|{
zR<sLPl<%lIte(^aHag1ZeNEE)_}5x6T*xYb=wb<r!P$I~5yOjAbtl|<lRHIZrfHhh
z0HqyT|L}G#{#R&ql=T;y<YYcczDP;CzJ;oora4hWUU9ix-(>|X%NkPE_eEp_U^vPg
zJI%*qk!Pad8oNapP^kB9;Q_-iT2wXM-QC@rI{}tuZ3eEIruhkQ9Rm!Zek^hTorqUq
zvIA!x{A~7fs`{kM<!Z^UD;Mnd`^Og*75$>BU6UqFTGrLoHJqtU3+t2uhjQCPv)Js9
zHOYb8-cTs?rm8NpZM&|wxA$x&Z{&hvu~>3=cz7N#abRHJsA-zDnc4-2r(L`*Hnqqt
zW5HmsDijLsRMi)3+kPUK{I&^lBM=BYV;IIdRkcNAm#Y4eE}L~FJ%9OZ&!!iTi^yZZ
z0>D(&4=l_2mYbjXQ8DIEAHX!tSt7Cks8Q8o5pe^3z$sPzxwN$O$e}}r?0fOtxA+ew
Wi}~-0c;>eN0000<MNUMnLSTZ*DnvH`

literal 0
HcmV?d00001

diff --git a/K2Engineering/K2Engineering/Resources/PressureOutput.png b/K2Engineering/K2Engineering/Resources/PressureOutput.png
new file mode 100644
index 0000000000000000000000000000000000000000..aa1f12dc246dfc8b6a381bd26419deaf9734364a
GIT binary patch
literal 1298
zcmV+t1?~EYP)<h;3K|Lk000e1NJLTq000;O000;W1^@s6;CDUv00004b3#c}2nYxW
zd<bNS00009a7bBm000ie000ie0hKEb8vp<R8FWQhbW?9;ba!ELWdL_~cP?peYja~^
zaAhuUa%Y?FJQ@H11dmBXK~zYIy_8FARCN@_f9KA0T4{r{Nr6tGG@S?5Vu5HPHZEFF
z)Jlp$TOcS4QG<~s2qIBHw1^sr*r?H@8y?aQ@>Uu`VuC=_gdhn7hg%?ljy$F#)R}3g
zz4sm$bV@t67#DhW=k<N(eE$dh$0qkmR#a4YQmIspVHkC)`jD!Q1F}UVA|maM;~W4E
zo2J?IZwZV*Ag~yC2RJ7p`+=Vl`QFa#s#o%j18W__Fs7(#9Z)SIP4ReqOC%CW{v$zY
zY3UfZ+r6XTU16H$*>uhu($7?>5#jmB3vwQyyu7@~vaGd0m8w2xnr8c8Jl9~Sq@*O@
z<MA95kv7vbm&RhT8w21wbs{c$TuP54xZ*AA&HzNC(ReHtJ6KRq5EYRf1qB6du~_Wp
z-!`bOu6B2HbQ~6uqw)8DJBRQl%E?H*EKk=IwrHM5Slv-Cs|Sc%bv8mB4#l}|2fqe3
zi%3;C9F77Q>B6g5uf7J{4u`{y0yc00;nHviG{kcF5y%yTTNf@|IHaoE9LM<}oi+d}
zD=TwFWR-}l08oox#o{ZJq&K8*jD+nyprc69<g(a+!y$m4o}Q0H<cUBaP(DDAOeSlA
z<Kb}lB0yxJT<KUK;YhtSk3=#U{K^7p?`V)vzY^WR7pnTw0D*|iRn<MiO&O_|W}v&c
zRp;D<ZV{SM6?)ieJFLI5K-x2VY}?)=B69}<fk5C?Fc>VG_>KCCTXpfUdZ-S$(xuIs
zlUY){TWh>6`eZ+5C*X4TZ%~x1|2$>8Xan|pL%MEoqgZV35D5>C#S6&tw(8m;GXMfL
z#k(}0HNoqXCr=*dzB8OU`KLq4*z}plY==`>y&N6fSQsbXkuhzwj=1+A>xU3HJg5#6
zvidjzZ0POnwE^(^{lAozl?4FC@6^)F)k~W-$N!yXXT~#a+_yTmqD4Ib!C){380kYk
zV;IH^fNKloawd3h9!-D7@dA`=*ZfT2NxQYED5Osr$^)r9#OL!(QPs<UWAsb3515zP
zSka=M-aNKC$|u(w<kG`mYF1x1Ywl{I;tlE6K8KT{%*SSBGC|m`tX-z6d*SX&ZWEEI
ze!ss0U_wZ1C+*gvKZGVn*?g@*&J8Dtx0Qv!tQ=)(q+VYBy+NKO(YdMbnN>@4bm5Kx
z0@E~8s`{>oYyk{MnPJ&{crS1Ogt-ErUTl!t=}_Lr$9<ZU+nri{V$|RbFimqCaNF<q
zZvqk+U<kD%g5R++70k5b_xl&AZO>9q?rGb>H-rGtL?W?RM5;!;6__x-@sTBwdigk=
z9N(gaH|}Q|2BaSo5EGj*0Ks5znTWh+Syr7G?$HR*Oada2NV2c5ZzeFhtE=m{-|w$Z
z=dMA_cJaDcR3rW0yl_eroi7H1!M&<F-?FU7tj3FO1fNSed#7Xs^Z9&F8-~%Ss+NfC
zQ`MjClIgeX*~=gGe7)$5h)mNl-g(_iTdE{_=}2Em<vLf^ePdBJqU?$^$ohK{4AOPF
zh|B>hRkcV&+`x6<ysCcRz48ZB^NU|anTik<hr{Wg^?$Ly0Ip95L-Pl{r2qf`07*qo
IM6N<$f@FhTtN;K2

literal 0
HcmV?d00001


From 5cbe5265c13c1ab563da2484c843ff6961bf4693 Mon Sep 17 00:00:00 2001
From: Cecilie Brandt-Olsen <c.brandtolsen@gmail.com>
Date: Sun, 11 Feb 2018 23:08:23 +0100
Subject: [PATCH 2/4] Fixed negative pressure bug

---
 K2Engineering/K2Engineering/Pressure.cs | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/K2Engineering/K2Engineering/Pressure.cs b/K2Engineering/K2Engineering/Pressure.cs
index 2f06f99..fd0e3ea 100644
--- a/K2Engineering/K2Engineering/Pressure.cs
+++ b/K2Engineering/K2Engineering/Pressure.cs
@@ -83,14 +83,23 @@ public class PressureGoal : GoalObject
             double temp;
 
             bool opt;
+            bool negativePres;
 
             public PressureGoal(PlanktonMesh mesh, double p, bool option)
             {
                 pMesh = (PlanktonMesh) mesh;
 
                 //Pressure
-                pres0 = p;      // Unit: [N/m2] / [Pa]
+                pres0 = Math.Abs(p);      // Unit: [N/m2] / [Pa]
                 pres1 = pres0;
+                if (p < 0.0)
+                {
+                    negativePres = true;
+                }
+                else
+                {
+                    negativePres = false;
+                }
 
                 //Calculate starting volume
                 vol0 = 0.0;         // Unit: [m3]
@@ -212,6 +221,10 @@ public override void Calculate(List<KangarooSolver.Particle> p)
                 {
                     Vector3d n = new Vector3d(normals[k]);
                     n.Unitize();
+                    if (negativePres)
+                    {
+                        n.Reverse();
+                    }
                     Move[k] = pres1 * vertexAreas[k] * n;
                 }
 
@@ -230,7 +243,15 @@ public override object Output(List<KangarooSolver.Particle> p)
                 }
 
                 //Create pressure data object to store output information
-                DataTypes.PressureData presData = new DataTypes.PressureData(vertices, pForces, Math.Round(pres0*1e-3,3), Math.Round(pres1*1e-3,3), Math.Round(vol0,3), Math.Round(vol1,3), Convert.ToInt32(mol0), Convert.ToInt32(mol1));
+                double p0 = pres0 * 1e-3;       //Unit: [kN/m2]
+                double p1 = pres1 * 1e-3;       //Unit: [kN/m2]
+                if (negativePres)
+                {
+                    p0 *= -1.0;
+                    p1 *= -1.0;
+                }
+
+                DataTypes.PressureData presData = new DataTypes.PressureData(vertices, pForces, Math.Round(p0,3), Math.Round(p1,3), Math.Round(vol0,3), Math.Round(vol1,3), Convert.ToInt32(mol0), Convert.ToInt32(mol1));
                 return presData;
             }
 

From 558e48fc358f3805c49c48d5e359ab423544e700 Mon Sep 17 00:00:00 2001
From: habibamoubarak <48219758+habibamoubarak@users.noreply.github.com>
Date: Wed, 20 Mar 2019 16:30:34 +0100
Subject: [PATCH 3/4] Updated beam / kinda working

---
 BeamHabiba.cs | 349 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 349 insertions(+)
 create mode 100644 BeamHabiba.cs

diff --git a/BeamHabiba.cs b/BeamHabiba.cs
new file mode 100644
index 0000000..e7cf793
--- /dev/null
+++ b/BeamHabiba.cs
@@ -0,0 +1,349 @@
+using System;
+using System.Collections.Generic;
+using Grasshopper.Kernel;
+using Rhino.Geometry;
+using KangarooSolver;
+
+namespace K2Engineering
+{
+     public class Beam : GH_Component
+     {
+         /// <summary>
+         /// Initializes a new instance of the Beam class.
+         /// </summary>
+         public Beam()
+           : base("Beam", "Beam",
+               "A goal that represents a beam element with biaxial bending and torsion behaviour",
+               "K2Eng", "0 Elements")
+         {
+         }
+
+         /// <summary>
+         /// Registers all the input parameters for this component.
+         /// </summary>
+         protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager)
+         {
+             pManager.AddPlaneParameter("StartPlane", "startPln", "The start plane of the beam (model in [m])", GH_ParamAccess.item);
+             pManager.AddPlaneParameter("EndPlane", "endPln", "The end plane of the beam (model in [m])", GH_ParamAccess.item);
+             pManager.AddNumberParameter("E-modulus", "E", "Young's modulus in [MPa]", GH_ParamAccess.item);
+             pManager.AddNumberParameter("G-modulus", "G", "The shear modulus in [MPa]", GH_ParamAccess.item);
+             pManager.AddNumberParameter("A", "A", "The cross section area in [mm2]", GH_ParamAccess.item);
+             pManager.AddNumberParameter("Iy", "Iy", "The moment of inertia about the cross section y-axis in [mm4]", GH_ParamAccess.item);
+             pManager.AddNumberParameter("Iz", "Iz", "The moment of inertia about the cross section z-axis in [mm4]", GH_ParamAccess.item);
+             pManager.AddNumberParameter("It", "It", "The torsional moment of inertia in [mm4]", GH_ParamAccess.item);
+         }
+
+         /// <summary>
+         /// Registers all the output parameters for this component.
+         /// </summary>
+         protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager)
+         {
+             pManager.AddGenericParameter("Beam", "Beam", "The beam goal", GH_ParamAccess.item);
+         }
+
+         /// <summary>
+         /// This is the method that actually does the work.
+         /// </summary>
+         /// <param name="DA">The DA object is used to retrieve from inputs and store in outputs.</param>
+         protected override void SolveInstance(IGH_DataAccess DA)
+         {
+             //Input
+             Plane startPln = new Plane();
+             DA.GetData(0, ref startPln);
+
+             Plane endPln = new Plane();
+             DA.GetData(1, ref endPln);
+
+             double eModulus = 0.0;
+             DA.GetData(2, ref eModulus);
+
+             double gModulus = 0.0;
+             DA.GetData(3, ref gModulus);
+
+             double area = 0.0;
+             DA.GetData(4, ref area);
+
+             double inertiaY = 0.0;
+             DA.GetData(5, ref inertiaY);
+
+             double inertiaZ = 0.0;
+             DA.GetData(6, ref inertiaZ);
+
+             double inertiaT = 0.0;
+             DA.GetData(7, ref inertiaT);
+
+
+             //Calculate
+             GoalObject beamElement = new BeamGoal(startPln, endPln, eModulus, gModulus, area, inertiaY, inertiaZ, inertiaT);
+
+
+             //Output
+             DA.SetData(0, beamElement);
+         }
+
+         //Define beam goal
+         
+         public class BeamGoal : GoalObject
+         {
+             Plane P0;               //Start plane (local)
+             Plane P1;               //End plane (local)
+             Plane P0R;              //Updated start plane (local)
+             Plane P1R;              //Updated end plane (local)
+
+             double restLength;
+             double E, G, A, Iy, Iz, It;
+             double thetaY0, thetaX0, thetaY1, thetaX1, thetaX;
+
+             double N, MY0, MX0, MY1, MX1, MZ; //, SZ, SY;
+            // Plane nodeO;
+
+             public BeamGoal(Plane startPlane, Plane endPlane, double eModulus, double gModulus, double area, double inertiaY, double inertiaZ, double inertiaT)
+             {
+                 restLength = startPlane.Origin.DistanceTo(endPlane.Origin);
+                 E = eModulus;
+                 G = gModulus;
+                 A = area;
+                 Iy = inertiaY;
+                 Iz = inertiaZ;
+                 It = inertiaT;
+
+
+                 //stiffness properties used to set the weightings
+                 double axialStiffness = (E * A) / restLength;                               //Unit: N/m
+                 int axialDigits = axialStiffness.ToString().Split('.')[0].Length;      // this is just making the above 10 raised to number of digits
+
+                 //max bending stiffness
+                 /*
+                 double bsY = (E * Iy * 1e-6) / restLength;                              //Unit: N*m
+                 double bsZ = (E * Iz * 1e-6) / restLength;                              //Unit: N*m 
+
+                 double bendingStiffness = bsY;
+                 if(bsZ > bsY)
+                 {
+                     bendingStiffness = bsZ;
+                 }
+                 
+                 */
+                double bendingStiffness = E * Math.Max(Iy, Iz);
+                bendingStiffness *= 0.01;                                                   //Still needs some adjustment to achieve acceptable convergence speed
+                int bendingDigits = bendingStiffness.ToString().Split('.')[0].Length;
+
+                double bsY = E * Iy * 0.01;
+                int bDy = bsY.ToString().Split('.')[0].Length;
+                double bsZ = E * Iz * 0.01;
+                int bDz = bsZ.ToString().Split('.')[0].Length;
+                double bst = G * It * 0.01;
+                int bDt = bst.ToString().Split('.')[0].Length;
+
+
+                //K2 properties
+                PPos = new Point3d[2] {startPlane.Origin, endPlane.Origin};
+                Move = new Vector3d[2];
+                Weighting = new double[2] { Math.Pow(10, axialDigits), Math.Pow(10, axialDigits) };           
+
+                Torque = new Vector3d[2];
+                //TorqueWeighting = new double[2] { Math.Pow(10, bendingDigits), Math.Pow(10, bendingDigits) };
+                TorqueWeighting = new double[3] { Math.Pow(10, bDy), Math.Pow(10, bDz), Math.Pow(10, bDt) };
+               
+                //Trial 14/03/2019
+                Plane startGlobal = new Plane(startPlane.Origin, Vector3d.XAxis, Vector3d.YAxis);
+                Plane endGlobal = new Plane(endPlane.Origin, Vector3d.XAxis, Vector3d.YAxis);
+                InitialOrientation = new Plane[2] { startGlobal, endGlobal };       //Needs to be global because the calculated forces and rotations in nodes are calculated globally
+                //P0.Transform(Transform.ChangeBasis(Plane.WorldXY, startGlobal));
+               // P1.Transform(Transform.ChangeBasis(Plane.WorldXY, endGlobal));
+
+
+                //Local end planes
+                P0 = startPlane;
+                P1 = endPlane;
+                P0R = P0;
+                P1R = P1;
+            }
+
+        
+            public override void Calculate(List<KangarooSolver.Particle> p)
+            {
+                //Get the current positions/orientations of the nodes (global: related to InitialOrientation)... I think the output of the beam frame is local
+               Plane P0Current = p[PIndex[0]].Orientation;
+                Plane P1Current = p[PIndex[1]].Orientation;
+
+                Vector3d elementVec = new Vector3d(P1Current.Origin - P0Current.Origin);
+                double currentLength = elementVec.Length;
+                Vector3d elementDir = new Vector3d(elementVec);
+                elementDir.Unitize();
+
+
+                //this works
+
+
+                P0R.Transform(Transform.PlaneToPlane(P0R, P0Current));
+                P1R.Transform(Transform.PlaneToPlane(P1R, P1Current));
+
+                LineCurve cur = new LineCurve(P0Current.Origin, P1Current.Origin);
+   
+                Plane perFrame0;
+                Plane perFrame1;
+                    cur.PerpendicularFrameAt((double)0, out perFrame0);
+                    cur.PerpendicularFrameAt((double)1, out perFrame1);
+                    P0R.Transform(Transform.PlaneToPlane(P0R, perFrame0));
+                    P1R.Transform(Transform.PlaneToPlane(P1R, perFrame1)); 
+            
+                  // nodeO = P0Current;
+              
+                P0R = new Plane(P0Current.Origin, P0Current.YAxis, -P0Current.ZAxis);
+                P1R = new Plane(P1Current.Origin, P1Current.YAxis, -P1Current.ZAxis); 
+
+                
+
+                //Mine
+               
+                Vector3d Y0 = P0R.YAxis;
+                Vector3d X0 = P0R.XAxis;
+                Vector3d Y1 = P1R.YAxis;
+                Vector3d X1 = P1R.XAxis;
+                //Bending angle changes around local axes
+                thetaX0 = 1.0 * Vector3d.Multiply(Y0, elementVec) / elementVec.Length;
+                thetaY1 = Vector3d.Multiply(X1, elementVec) / elementVec.Length;
+
+                thetaY0 = 1.0 * Vector3d.Multiply(X0, elementVec) / elementVec.Length;     //NB! Negative sign /// right hand rule should in reality be theta x
+                
+                thetaX1 = 1.0 * Vector3d.Multiply(Y1, elementVec) / elementVec.Length;     //NB! Negative sign
+
+                //Twist angle change around element axis //which here we're taking to be the x
+                thetaX = (Vector3d.Multiply(X0, Y1) - Vector3d.Multiply(X1, Y0)) / 2.0;
+
+
+                //Axial
+
+                //double extA = (Math.Pow(currentLength, 2) - Math.Pow(restLength, 2)) / (2.0 * restLength);
+                //double extB = (restLength / 60.0) * (4.0 * (Math.Pow(thetaY0, 2) + Math.Pow(thetaZ0, 2)) - 2.0 * ((thetaY0 * thetaY1) + (thetaZ0 * thetaZ1)) + 4.0 * (Math.Pow(thetaY1, 2) + Math.Pow(thetaZ1, 2)));
+                //double extension = extA + extB;        //Unit: [m]
+
+                //For now let's use the most simple expression similar to the Bar Goal (the effect from bowing should be relatively small)
+                double extension = currentLength - restLength;        //Unit: [m]
+
+                //Element internal forces
+                N = ((E * A) / restLength) * extension;      // Unit: [N]
+                
+                MX0 = ( ((N * restLength) / 30.0) * ((4.0 * thetaX0) - thetaX1) ) + ( ((E * Iy * 1e-6) / restLength) * ((4.0 * thetaX0) + (2.0 * thetaX1)) );           //Unit: [Nm]
+                MY0 = ( ((N * restLength) / 30.0) * ((4.0 * thetaY0) - thetaY1) ) + ( ((E * Iz * 1e-6) / restLength) * ((4.0 * thetaY0) + (2.0 * thetaY1)) );           //Unit: [Nm]
+
+                MX1 = ( ((N * restLength) / 30.0) * ((4.0 * thetaX1) - thetaX0) ) + ( ((E * Iy * 1e-6) / restLength) * ((4.0 * thetaX1) + (2.0 * thetaX0)) );           //Unit: [Nm]
+                MY1 = ( ((N * restLength) / 30.0) * ((4.0 * thetaY1) - thetaY0) ) + ( ((E * Iz * 1e-6) / restLength) * ((4.0 * thetaY1) + (2.0 * thetaY0)) );           //Unit: [Nm]
+
+                MZ = ((G * It * 1e-6) / restLength) * thetaX;            //Unit: [Nm]
+
+                //Global forces
+                //Force start
+                double F0X = (1.0 / restLength) * ((N * elementVec.X) + (MX0 * Y0.X) - (MY0 * X0.X) + (MX1 * Y1.X) - (MY1 * X1.X));
+                double F0Y = (1.0 / restLength) * ((N * elementVec.Y) + (MX0 * Y0.Y) - (MY0 * X0.Y) + (MX1 * Y1.Y) - (MY1 * X1.Y));
+                double F0Z = (1.0 / restLength) * ((N * elementVec.Z) + (MX0 * Y0.Z) - (MY0 * X0.Z) + (MX1 * Y1.Z) - (MY1 * X1.Z)); 
+                Vector3d F0 = new Vector3d(F0X, F0Y, F0Z);          //Unit: [N]
+                //Vector3d F0 = new Vector3d(N, 0, 0);          //Unit: [N]
+
+                //Force end
+                Vector3d F1 = -1.0 * F0;
+
+
+                //Permutation symbol: Includes 6 non-zero components. Is a triple product of vectors (elementVec,y,z) of an orthogonal frame in a right-handed coordinate system
+
+                //Moment start
+                //i=1, j=2, k=3
+                double M0X_pos = (-1.0) * (((MX0 * elementVec.Z * Y0.Y) / restLength) - ((MY0 * elementVec.Z * X0.Y) / restLength) + ((MZ * ((X0.Y * Y1.Z) - (Y0.Y * X1.Z))) / 2.0));
+
+                //i=1, j=3, k=2
+                double M0X_neg = (1.0) * (((MX0 * elementVec.Y * Y0.Z) / restLength) - ((MY0 * elementVec.Y * X0.Z) / restLength) + ((MZ * ((X0.Z * Y1.Y) - (Y0.Z * X1.Y))) / 2.0));
+
+                //i=2, j=3, k=1
+                double M0Y_pos = (-1.0) * (((MX0 * elementVec.X * Y0.Z) / restLength) - ((MY0 * elementVec.X * X0.Z) / restLength) + ((MZ * ((X0.Z * Y1.X) - (Y0.Z * X1.X))) / 2.0));
+
+                //i=2, j=1, k=3
+                double M0Y_neg = (1.0) * (((MX0 * elementVec.Z * Y0.X) / restLength) - ((MY0 * elementVec.Z * X0.X) / restLength) + ((MZ * ((X0.X * Y1.Z) - (Y0.X * X1.Z))) / 2.0));
+
+                //i=3, j=1, k=2
+                double M0Z_pos = (-1.0) * (((MX0 * elementVec.Y * Y0.X) / restLength) - ((MY0 * elementVec.Y * X0.X) / restLength) + ((MZ * ((X0.X * Y1.Y) - (Y0.X * X1.Y))) / 2.0));
+
+                //i=3, j=2, k=1
+                double M0Z_neg = (1.0) * (((MX0 * elementVec.X * Y0.Y) / restLength) - ((MY0 * elementVec.X * X0.Y) / restLength) + ((MZ * ((X0.Y * Y1.X) - (Y0.Y * X1.X))) / 2.0));
+
+                //Sum of components
+              Vector3d M0 = new Vector3d(M0X_pos + M0X_neg, M0Y_pos + M0Y_neg, M0Z_pos + M0Z_neg);          //Unit: [Nm]
+                
+               // Vector3d M0 = new Vector3d(M0X_pos, M0Y_pos , M0Z_pos );
+               
+
+                //Moment end
+                //i=1, j=2, k=3
+                double M1X_pos = (-1.0) * (((MX1 * elementVec.Z * Y1.Y) / restLength) - ((MY1 * elementVec.Z * X1.Y) / restLength) + ((MZ * ((X0.Y * Y1.Z) - (Y0.Y * X1.Z))) / 2.0));
+
+                //i=1, j=3, k=2
+                double M1X_neg = (1.0) * (((MX1 * elementVec.Y * Y1.Z) / restLength) - ((MY1 * elementVec.Y * X1.Z) / restLength) + ((MZ * ((X0.Z * Y1.Y) - (Y0.Z * X1.Y))) / 2.0));
+
+                //i=2, j=3, k=1
+                double M1Y_pos = (-1.0) * (((MX1 * elementVec.X * Y1.Z) / restLength) - ((MY1 * elementVec.X * X1.Z) / restLength) + ((MZ * ((X0.Z * Y1.X) - (Y0.Z * X1.X))) / 2.0));
+
+                //i=2, j=1, k=3
+               double M1Y_neg = (1.0) * (((MX1 * elementVec.Z * Y1.X) / restLength) - ((MY1 * elementVec.Z * X1.X) / restLength) + ((MZ * ((X0.X * Y1.Z) - (Y0.X * X1.Z))) / 2.0));
+
+                //i=3, j=1, k=2
+                double M1Z_pos = (-1.0) * (((MX1 * elementVec.Y * Y1.X) / restLength) - ((MY1 * elementVec.Y * X1.X) / restLength) + ((MZ * ((X0.X * Y1.Y) - (Y0.X * X1.Y))) / 2.0));
+
+                //i=3, j=2, k=1
+                double M1Z_neg = (1.0) * (((MX1 * elementVec.X * Y1.Y) / restLength) - ((MY1 * elementVec.X * X1.Y) / restLength) + ((MZ * ((X0.Y * Y1.X) - (Y0.Y * X1.X))) / 2.0));
+
+                //Sum of components
+                Vector3d M1 = new Vector3d(M1X_pos + M1X_neg, M1Y_pos + M1Y_neg, M1Z_pos + M1Z_neg);            //Unit: [Nm]
+                //Vector3d M1 = new Vector3d(M1X_pos, M1Y_pos, M1Z_pos);
+               
+
+                //Move and torque vectors
+                Move[0] = F0 / Weighting[0];
+                Move[1] = F1 / Weighting[1];
+                ///Move[0] = elementDir*(extension/2); //moving the start point by this much, but in every direction? NO the element dir VECTOR is what tell us where the node extension applies and how much it is... 
+                ///Move[1] = -elementDir * (extension / 2);
+
+                Torque[0] = M0 / TorqueWeighting[0];
+                Torque[1] = M1 / TorqueWeighting[1];
+                //Torque[0] = new Vector3d((M0X_pos) / TorqueWeighting[0], (M0Y_pos) / TorqueWeighting[1], (M0Z_pos) / TorqueWeighting[2]);
+                //Torque[1] = new Vector3d((M1X_pos) / TorqueWeighting[0], (M1Y_pos) / TorqueWeighting[1], (M1Z_pos) / TorqueWeighting[2]);
+                
+                
+                //above doesnt work.. will have to try with matrices
+                //Matrix matrix1 = new Matrix(3, 1);
+                // Torque[0] = matrix1[1,1], ;
+            }
+            
+
+            //Output moment in [kNm] and normal force/shear in [kN]
+            public override object Output(List<KangarooSolver.Particle> p)
+            {
+                DataTypes.BeamData beamData = new DataTypes.BeamData(P0R, P1R, Math.Round(N * 1e-3, 3), Math.Round(MZ * 1e-3, 3), Math.Round(MX0 * 1e-3, 3), Math.Round(MY0 * 1e-3, 3), Math.Round(MX1 * 1e-3, 3), Math.Round(MY1 * 1e-3, 3)/* nodeO*/);
+                return beamData;
+            }
+
+        } 
+
+
+
+        /// <summary>
+        /// Provides an Icon for the component.
+        /// </summary>
+        protected override System.Drawing.Bitmap Icon
+        {
+            get
+            {
+                //You can add image files to your project resources and access them like this:
+                // return Resources.IconForThisComponent;
+                return null;
+            }
+        }
+
+        /// <summary>
+        /// Gets the unique ID for this component. Do not change this ID after release.
+        /// </summary>
+        public override Guid ComponentGuid
+        {
+            get { return new Guid("f76cfabb-06ab-4869-a60d-c6c4865320ea"); }
+        }
+    }
+}
\ No newline at end of file

From 63518147e3cdbe9931f6ec340f837c5f6811b7ad Mon Sep 17 00:00:00 2001
From: habibamoubarak <48219758+habibamoubarak@users.noreply.github.com>
Date: Thu, 21 Mar 2019 14:28:48 +0100
Subject: [PATCH 4/4] Delete BeamHabiba.cs

---
 BeamHabiba.cs | 349 --------------------------------------------------
 1 file changed, 349 deletions(-)
 delete mode 100644 BeamHabiba.cs

diff --git a/BeamHabiba.cs b/BeamHabiba.cs
deleted file mode 100644
index e7cf793..0000000
--- a/BeamHabiba.cs
+++ /dev/null
@@ -1,349 +0,0 @@
-using System;
-using System.Collections.Generic;
-using Grasshopper.Kernel;
-using Rhino.Geometry;
-using KangarooSolver;
-
-namespace K2Engineering
-{
-     public class Beam : GH_Component
-     {
-         /// <summary>
-         /// Initializes a new instance of the Beam class.
-         /// </summary>
-         public Beam()
-           : base("Beam", "Beam",
-               "A goal that represents a beam element with biaxial bending and torsion behaviour",
-               "K2Eng", "0 Elements")
-         {
-         }
-
-         /// <summary>
-         /// Registers all the input parameters for this component.
-         /// </summary>
-         protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager)
-         {
-             pManager.AddPlaneParameter("StartPlane", "startPln", "The start plane of the beam (model in [m])", GH_ParamAccess.item);
-             pManager.AddPlaneParameter("EndPlane", "endPln", "The end plane of the beam (model in [m])", GH_ParamAccess.item);
-             pManager.AddNumberParameter("E-modulus", "E", "Young's modulus in [MPa]", GH_ParamAccess.item);
-             pManager.AddNumberParameter("G-modulus", "G", "The shear modulus in [MPa]", GH_ParamAccess.item);
-             pManager.AddNumberParameter("A", "A", "The cross section area in [mm2]", GH_ParamAccess.item);
-             pManager.AddNumberParameter("Iy", "Iy", "The moment of inertia about the cross section y-axis in [mm4]", GH_ParamAccess.item);
-             pManager.AddNumberParameter("Iz", "Iz", "The moment of inertia about the cross section z-axis in [mm4]", GH_ParamAccess.item);
-             pManager.AddNumberParameter("It", "It", "The torsional moment of inertia in [mm4]", GH_ParamAccess.item);
-         }
-
-         /// <summary>
-         /// Registers all the output parameters for this component.
-         /// </summary>
-         protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager)
-         {
-             pManager.AddGenericParameter("Beam", "Beam", "The beam goal", GH_ParamAccess.item);
-         }
-
-         /// <summary>
-         /// This is the method that actually does the work.
-         /// </summary>
-         /// <param name="DA">The DA object is used to retrieve from inputs and store in outputs.</param>
-         protected override void SolveInstance(IGH_DataAccess DA)
-         {
-             //Input
-             Plane startPln = new Plane();
-             DA.GetData(0, ref startPln);
-
-             Plane endPln = new Plane();
-             DA.GetData(1, ref endPln);
-
-             double eModulus = 0.0;
-             DA.GetData(2, ref eModulus);
-
-             double gModulus = 0.0;
-             DA.GetData(3, ref gModulus);
-
-             double area = 0.0;
-             DA.GetData(4, ref area);
-
-             double inertiaY = 0.0;
-             DA.GetData(5, ref inertiaY);
-
-             double inertiaZ = 0.0;
-             DA.GetData(6, ref inertiaZ);
-
-             double inertiaT = 0.0;
-             DA.GetData(7, ref inertiaT);
-
-
-             //Calculate
-             GoalObject beamElement = new BeamGoal(startPln, endPln, eModulus, gModulus, area, inertiaY, inertiaZ, inertiaT);
-
-
-             //Output
-             DA.SetData(0, beamElement);
-         }
-
-         //Define beam goal
-         
-         public class BeamGoal : GoalObject
-         {
-             Plane P0;               //Start plane (local)
-             Plane P1;               //End plane (local)
-             Plane P0R;              //Updated start plane (local)
-             Plane P1R;              //Updated end plane (local)
-
-             double restLength;
-             double E, G, A, Iy, Iz, It;
-             double thetaY0, thetaX0, thetaY1, thetaX1, thetaX;
-
-             double N, MY0, MX0, MY1, MX1, MZ; //, SZ, SY;
-            // Plane nodeO;
-
-             public BeamGoal(Plane startPlane, Plane endPlane, double eModulus, double gModulus, double area, double inertiaY, double inertiaZ, double inertiaT)
-             {
-                 restLength = startPlane.Origin.DistanceTo(endPlane.Origin);
-                 E = eModulus;
-                 G = gModulus;
-                 A = area;
-                 Iy = inertiaY;
-                 Iz = inertiaZ;
-                 It = inertiaT;
-
-
-                 //stiffness properties used to set the weightings
-                 double axialStiffness = (E * A) / restLength;                               //Unit: N/m
-                 int axialDigits = axialStiffness.ToString().Split('.')[0].Length;      // this is just making the above 10 raised to number of digits
-
-                 //max bending stiffness
-                 /*
-                 double bsY = (E * Iy * 1e-6) / restLength;                              //Unit: N*m
-                 double bsZ = (E * Iz * 1e-6) / restLength;                              //Unit: N*m 
-
-                 double bendingStiffness = bsY;
-                 if(bsZ > bsY)
-                 {
-                     bendingStiffness = bsZ;
-                 }
-                 
-                 */
-                double bendingStiffness = E * Math.Max(Iy, Iz);
-                bendingStiffness *= 0.01;                                                   //Still needs some adjustment to achieve acceptable convergence speed
-                int bendingDigits = bendingStiffness.ToString().Split('.')[0].Length;
-
-                double bsY = E * Iy * 0.01;
-                int bDy = bsY.ToString().Split('.')[0].Length;
-                double bsZ = E * Iz * 0.01;
-                int bDz = bsZ.ToString().Split('.')[0].Length;
-                double bst = G * It * 0.01;
-                int bDt = bst.ToString().Split('.')[0].Length;
-
-
-                //K2 properties
-                PPos = new Point3d[2] {startPlane.Origin, endPlane.Origin};
-                Move = new Vector3d[2];
-                Weighting = new double[2] { Math.Pow(10, axialDigits), Math.Pow(10, axialDigits) };           
-
-                Torque = new Vector3d[2];
-                //TorqueWeighting = new double[2] { Math.Pow(10, bendingDigits), Math.Pow(10, bendingDigits) };
-                TorqueWeighting = new double[3] { Math.Pow(10, bDy), Math.Pow(10, bDz), Math.Pow(10, bDt) };
-               
-                //Trial 14/03/2019
-                Plane startGlobal = new Plane(startPlane.Origin, Vector3d.XAxis, Vector3d.YAxis);
-                Plane endGlobal = new Plane(endPlane.Origin, Vector3d.XAxis, Vector3d.YAxis);
-                InitialOrientation = new Plane[2] { startGlobal, endGlobal };       //Needs to be global because the calculated forces and rotations in nodes are calculated globally
-                //P0.Transform(Transform.ChangeBasis(Plane.WorldXY, startGlobal));
-               // P1.Transform(Transform.ChangeBasis(Plane.WorldXY, endGlobal));
-
-
-                //Local end planes
-                P0 = startPlane;
-                P1 = endPlane;
-                P0R = P0;
-                P1R = P1;
-            }
-
-        
-            public override void Calculate(List<KangarooSolver.Particle> p)
-            {
-                //Get the current positions/orientations of the nodes (global: related to InitialOrientation)... I think the output of the beam frame is local
-               Plane P0Current = p[PIndex[0]].Orientation;
-                Plane P1Current = p[PIndex[1]].Orientation;
-
-                Vector3d elementVec = new Vector3d(P1Current.Origin - P0Current.Origin);
-                double currentLength = elementVec.Length;
-                Vector3d elementDir = new Vector3d(elementVec);
-                elementDir.Unitize();
-
-
-                //this works
-
-
-                P0R.Transform(Transform.PlaneToPlane(P0R, P0Current));
-                P1R.Transform(Transform.PlaneToPlane(P1R, P1Current));
-
-                LineCurve cur = new LineCurve(P0Current.Origin, P1Current.Origin);
-   
-                Plane perFrame0;
-                Plane perFrame1;
-                    cur.PerpendicularFrameAt((double)0, out perFrame0);
-                    cur.PerpendicularFrameAt((double)1, out perFrame1);
-                    P0R.Transform(Transform.PlaneToPlane(P0R, perFrame0));
-                    P1R.Transform(Transform.PlaneToPlane(P1R, perFrame1)); 
-            
-                  // nodeO = P0Current;
-              
-                P0R = new Plane(P0Current.Origin, P0Current.YAxis, -P0Current.ZAxis);
-                P1R = new Plane(P1Current.Origin, P1Current.YAxis, -P1Current.ZAxis); 
-
-                
-
-                //Mine
-               
-                Vector3d Y0 = P0R.YAxis;
-                Vector3d X0 = P0R.XAxis;
-                Vector3d Y1 = P1R.YAxis;
-                Vector3d X1 = P1R.XAxis;
-                //Bending angle changes around local axes
-                thetaX0 = 1.0 * Vector3d.Multiply(Y0, elementVec) / elementVec.Length;
-                thetaY1 = Vector3d.Multiply(X1, elementVec) / elementVec.Length;
-
-                thetaY0 = 1.0 * Vector3d.Multiply(X0, elementVec) / elementVec.Length;     //NB! Negative sign /// right hand rule should in reality be theta x
-                
-                thetaX1 = 1.0 * Vector3d.Multiply(Y1, elementVec) / elementVec.Length;     //NB! Negative sign
-
-                //Twist angle change around element axis //which here we're taking to be the x
-                thetaX = (Vector3d.Multiply(X0, Y1) - Vector3d.Multiply(X1, Y0)) / 2.0;
-
-
-                //Axial
-
-                //double extA = (Math.Pow(currentLength, 2) - Math.Pow(restLength, 2)) / (2.0 * restLength);
-                //double extB = (restLength / 60.0) * (4.0 * (Math.Pow(thetaY0, 2) + Math.Pow(thetaZ0, 2)) - 2.0 * ((thetaY0 * thetaY1) + (thetaZ0 * thetaZ1)) + 4.0 * (Math.Pow(thetaY1, 2) + Math.Pow(thetaZ1, 2)));
-                //double extension = extA + extB;        //Unit: [m]
-
-                //For now let's use the most simple expression similar to the Bar Goal (the effect from bowing should be relatively small)
-                double extension = currentLength - restLength;        //Unit: [m]
-
-                //Element internal forces
-                N = ((E * A) / restLength) * extension;      // Unit: [N]
-                
-                MX0 = ( ((N * restLength) / 30.0) * ((4.0 * thetaX0) - thetaX1) ) + ( ((E * Iy * 1e-6) / restLength) * ((4.0 * thetaX0) + (2.0 * thetaX1)) );           //Unit: [Nm]
-                MY0 = ( ((N * restLength) / 30.0) * ((4.0 * thetaY0) - thetaY1) ) + ( ((E * Iz * 1e-6) / restLength) * ((4.0 * thetaY0) + (2.0 * thetaY1)) );           //Unit: [Nm]
-
-                MX1 = ( ((N * restLength) / 30.0) * ((4.0 * thetaX1) - thetaX0) ) + ( ((E * Iy * 1e-6) / restLength) * ((4.0 * thetaX1) + (2.0 * thetaX0)) );           //Unit: [Nm]
-                MY1 = ( ((N * restLength) / 30.0) * ((4.0 * thetaY1) - thetaY0) ) + ( ((E * Iz * 1e-6) / restLength) * ((4.0 * thetaY1) + (2.0 * thetaY0)) );           //Unit: [Nm]
-
-                MZ = ((G * It * 1e-6) / restLength) * thetaX;            //Unit: [Nm]
-
-                //Global forces
-                //Force start
-                double F0X = (1.0 / restLength) * ((N * elementVec.X) + (MX0 * Y0.X) - (MY0 * X0.X) + (MX1 * Y1.X) - (MY1 * X1.X));
-                double F0Y = (1.0 / restLength) * ((N * elementVec.Y) + (MX0 * Y0.Y) - (MY0 * X0.Y) + (MX1 * Y1.Y) - (MY1 * X1.Y));
-                double F0Z = (1.0 / restLength) * ((N * elementVec.Z) + (MX0 * Y0.Z) - (MY0 * X0.Z) + (MX1 * Y1.Z) - (MY1 * X1.Z)); 
-                Vector3d F0 = new Vector3d(F0X, F0Y, F0Z);          //Unit: [N]
-                //Vector3d F0 = new Vector3d(N, 0, 0);          //Unit: [N]
-
-                //Force end
-                Vector3d F1 = -1.0 * F0;
-
-
-                //Permutation symbol: Includes 6 non-zero components. Is a triple product of vectors (elementVec,y,z) of an orthogonal frame in a right-handed coordinate system
-
-                //Moment start
-                //i=1, j=2, k=3
-                double M0X_pos = (-1.0) * (((MX0 * elementVec.Z * Y0.Y) / restLength) - ((MY0 * elementVec.Z * X0.Y) / restLength) + ((MZ * ((X0.Y * Y1.Z) - (Y0.Y * X1.Z))) / 2.0));
-
-                //i=1, j=3, k=2
-                double M0X_neg = (1.0) * (((MX0 * elementVec.Y * Y0.Z) / restLength) - ((MY0 * elementVec.Y * X0.Z) / restLength) + ((MZ * ((X0.Z * Y1.Y) - (Y0.Z * X1.Y))) / 2.0));
-
-                //i=2, j=3, k=1
-                double M0Y_pos = (-1.0) * (((MX0 * elementVec.X * Y0.Z) / restLength) - ((MY0 * elementVec.X * X0.Z) / restLength) + ((MZ * ((X0.Z * Y1.X) - (Y0.Z * X1.X))) / 2.0));
-
-                //i=2, j=1, k=3
-                double M0Y_neg = (1.0) * (((MX0 * elementVec.Z * Y0.X) / restLength) - ((MY0 * elementVec.Z * X0.X) / restLength) + ((MZ * ((X0.X * Y1.Z) - (Y0.X * X1.Z))) / 2.0));
-
-                //i=3, j=1, k=2
-                double M0Z_pos = (-1.0) * (((MX0 * elementVec.Y * Y0.X) / restLength) - ((MY0 * elementVec.Y * X0.X) / restLength) + ((MZ * ((X0.X * Y1.Y) - (Y0.X * X1.Y))) / 2.0));
-
-                //i=3, j=2, k=1
-                double M0Z_neg = (1.0) * (((MX0 * elementVec.X * Y0.Y) / restLength) - ((MY0 * elementVec.X * X0.Y) / restLength) + ((MZ * ((X0.Y * Y1.X) - (Y0.Y * X1.X))) / 2.0));
-
-                //Sum of components
-              Vector3d M0 = new Vector3d(M0X_pos + M0X_neg, M0Y_pos + M0Y_neg, M0Z_pos + M0Z_neg);          //Unit: [Nm]
-                
-               // Vector3d M0 = new Vector3d(M0X_pos, M0Y_pos , M0Z_pos );
-               
-
-                //Moment end
-                //i=1, j=2, k=3
-                double M1X_pos = (-1.0) * (((MX1 * elementVec.Z * Y1.Y) / restLength) - ((MY1 * elementVec.Z * X1.Y) / restLength) + ((MZ * ((X0.Y * Y1.Z) - (Y0.Y * X1.Z))) / 2.0));
-
-                //i=1, j=3, k=2
-                double M1X_neg = (1.0) * (((MX1 * elementVec.Y * Y1.Z) / restLength) - ((MY1 * elementVec.Y * X1.Z) / restLength) + ((MZ * ((X0.Z * Y1.Y) - (Y0.Z * X1.Y))) / 2.0));
-
-                //i=2, j=3, k=1
-                double M1Y_pos = (-1.0) * (((MX1 * elementVec.X * Y1.Z) / restLength) - ((MY1 * elementVec.X * X1.Z) / restLength) + ((MZ * ((X0.Z * Y1.X) - (Y0.Z * X1.X))) / 2.0));
-
-                //i=2, j=1, k=3
-               double M1Y_neg = (1.0) * (((MX1 * elementVec.Z * Y1.X) / restLength) - ((MY1 * elementVec.Z * X1.X) / restLength) + ((MZ * ((X0.X * Y1.Z) - (Y0.X * X1.Z))) / 2.0));
-
-                //i=3, j=1, k=2
-                double M1Z_pos = (-1.0) * (((MX1 * elementVec.Y * Y1.X) / restLength) - ((MY1 * elementVec.Y * X1.X) / restLength) + ((MZ * ((X0.X * Y1.Y) - (Y0.X * X1.Y))) / 2.0));
-
-                //i=3, j=2, k=1
-                double M1Z_neg = (1.0) * (((MX1 * elementVec.X * Y1.Y) / restLength) - ((MY1 * elementVec.X * X1.Y) / restLength) + ((MZ * ((X0.Y * Y1.X) - (Y0.Y * X1.X))) / 2.0));
-
-                //Sum of components
-                Vector3d M1 = new Vector3d(M1X_pos + M1X_neg, M1Y_pos + M1Y_neg, M1Z_pos + M1Z_neg);            //Unit: [Nm]
-                //Vector3d M1 = new Vector3d(M1X_pos, M1Y_pos, M1Z_pos);
-               
-
-                //Move and torque vectors
-                Move[0] = F0 / Weighting[0];
-                Move[1] = F1 / Weighting[1];
-                ///Move[0] = elementDir*(extension/2); //moving the start point by this much, but in every direction? NO the element dir VECTOR is what tell us where the node extension applies and how much it is... 
-                ///Move[1] = -elementDir * (extension / 2);
-
-                Torque[0] = M0 / TorqueWeighting[0];
-                Torque[1] = M1 / TorqueWeighting[1];
-                //Torque[0] = new Vector3d((M0X_pos) / TorqueWeighting[0], (M0Y_pos) / TorqueWeighting[1], (M0Z_pos) / TorqueWeighting[2]);
-                //Torque[1] = new Vector3d((M1X_pos) / TorqueWeighting[0], (M1Y_pos) / TorqueWeighting[1], (M1Z_pos) / TorqueWeighting[2]);
-                
-                
-                //above doesnt work.. will have to try with matrices
-                //Matrix matrix1 = new Matrix(3, 1);
-                // Torque[0] = matrix1[1,1], ;
-            }
-            
-
-            //Output moment in [kNm] and normal force/shear in [kN]
-            public override object Output(List<KangarooSolver.Particle> p)
-            {
-                DataTypes.BeamData beamData = new DataTypes.BeamData(P0R, P1R, Math.Round(N * 1e-3, 3), Math.Round(MZ * 1e-3, 3), Math.Round(MX0 * 1e-3, 3), Math.Round(MY0 * 1e-3, 3), Math.Round(MX1 * 1e-3, 3), Math.Round(MY1 * 1e-3, 3)/* nodeO*/);
-                return beamData;
-            }
-
-        } 
-
-
-
-        /// <summary>
-        /// Provides an Icon for the component.
-        /// </summary>
-        protected override System.Drawing.Bitmap Icon
-        {
-            get
-            {
-                //You can add image files to your project resources and access them like this:
-                // return Resources.IconForThisComponent;
-                return null;
-            }
-        }
-
-        /// <summary>
-        /// Gets the unique ID for this component. Do not change this ID after release.
-        /// </summary>
-        public override Guid ComponentGuid
-        {
-            get { return new Guid("f76cfabb-06ab-4869-a60d-c6c4865320ea"); }
-        }
-    }
-}
\ No newline at end of file