using FeMM.Common.Models;
using FeMM.Grasshopper.ComponentAttributes;
using FeMM.Grasshopper.DataTypes.FeMM;
using FeMM.Grasshopper.Helpers;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Parameters;
using Grasshopper.Kernel.Types;
using Rhino;
using Rhino.Geometry;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.Versioning;
using TensileLib.FormFinding;
using TensileLib.Geometry;
namespace FeMM.Grasshopper.Components.Analysis
{
#if NETCOREAPP
[SupportedOSPlatform("windows")]
#endif
public class FormFindComponent : GH_Component
{
private bool _run;
private GH_Model _model;
private GH_Mesh _mesh;
private int _runCount;
protected RhinoDoc _doc;
private List<string> _log;
/// <summary>
/// Initializes a new instance of the FormFindComponent class.
/// </summary>
public FormFindComponent()
: base("Form Find", "FF", "Calculate the net cable form", CategoryNameConstants.CATEGORY_FEMM, CategoryNameConstants.SUBCATEGORY_ANALYSIS)
{
_run = false;
_model = new GH_Model();
_mesh = new GH_Mesh();
_runCount = 0;
_log = [];
}
public override void CreateAttributes()
{
var runBtn = new ComponentAttributes.ComponentNButtonsAttributes.ButtonBase("Run");
var btns = new List<ComponentNButtonsAttributes.ButtonBase>() { runBtn };
var attr = new ComponentNButtonsAttributes(this, btns);
attr.ButtonPressed += (i) =>
{
if (i == 0)
{
_run = true;
ExpireSolution(true);
}
};
m_attributes = attr;
}
/// <summary>
/// Registers all the input parameters for this component.
/// </summary>
protected override void RegisterInputParams(GH_InputParamManager pManager)
{
pManager.AddGenericParameter("Model", "M", "The FEM model", GH_ParamAccess.item);
pManager.AddNumberParameter("Gravity", "G", "The gravity acceleration for self-weight calculation purpose", GH_ParamAccess.item, 0);
pManager.AddNumberParameter("Pressure", "P", "The normal pressure applied to the model", GH_ParamAccess.item, 0);
int i = pManager.AddGenericParameter("Pressure Load Case", "LC", "The load case used for the resulting pressure forces", GH_ParamAccess.item);
pManager[i].Optional = true;
i = pManager.AddIntegerParameter("Solver Type", "ST", "The type of solver to be used", GH_ParamAccess.item, 0);
Param_Integer mtParam = pManager[i] as Param_Integer;
mtParam.AddNamedValue("Linear", 0);
mtParam.AddNamedValue("Linear step by step", 1);
mtParam.AddNamedValue("Non-Linear with unit weights", 2);
mtParam.AddNamedValue("Non-Linear with negligible weights", 3);
pManager.AddIntegerParameter("Max Iter", "MI", "Max number of iterations in solution (only for Linear step by step and Non-Linear solution)", GH_ParamAccess.item, 20);
pManager.AddNumberParameter("Tolerance", "T", "Tolerance on relative error durign solution", GH_ParamAccess.item, 1.0E-4);
pManager.AddIntegerParameter("Max Iter,s", "MI,s", "Max number of iterations for the algebric solver", GH_ParamAccess.item, 0);
pManager.AddNumberParameter("ToleranceA,s", "T,s", "Tolerance for the algebric solver", GH_ParamAccess.item, 0);
pManager.AddNumberParameter("ToleranceB,s", "T,s", "Tolerance for the algebric solver", GH_ParamAccess.item, 0);
pManager.AddBooleanParameter("Create mesh", "CM", "Automatically create the plates and the mesh", GH_ParamAccess.item, true);
i = pManager.AddGenericParameter("Plate Property", "PP", "The property to assign to the generated plates. Create mesh must be set to true", GH_ParamAccess.item);
pManager[i].Optional = true;
pManager.AddPlaneParameter("Plane", "P", "Plane of mesh generation", GH_ParamAccess.item, Plane.WorldXY);
pManager.AddBooleanParameter("Auto Run", "AR", "Activate or disable the live update", GH_ParamAccess.item, false);
}
/// <summary>
/// Registers all the output parameters for this component.
/// </summary>
protected override void RegisterOutputParams(GH_OutputParamManager pManager)
{
pManager.AddGenericParameter("Model", "M", "The FEM model", GH_ParamAccess.item);
pManager.AddMeshParameter("Mesh", "M", "The resulting mesh usefull to geometric operations", GH_ParamAccess.item);
pManager.AddTextParameter("Log", "L", "The analysis log", GH_ParamAccess.list);
}
/// <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)
{
Message = "";
GH_Model model = null;
double gravity = 0;
double pressure = 0;
GH_Case lc_pressure = null;
int solverType = 0;
int maxIterNumber = 0;
double tolerance = 0;
int maxIterNumberSolver = 0;
double toleranceSolverA = 0;
double toleranceSolverB = 0;
bool makeMesh = true;
GH_FunctionDefinition pp = null;
bool autoRun = false;
Plane plane = Plane.Unset;
if (!DA.GetData("Model", ref model))
return;
DA.GetData("Gravity", ref gravity);
DA.GetData("Pressure", ref pressure);
DA.GetData("Pressure Load Case", ref lc_pressure);
DA.GetData(4, ref solverType);
DA.GetData(5, ref maxIterNumber);
DA.GetData(6, ref tolerance);
DA.GetData(7, ref maxIterNumberSolver);
DA.GetData(8, ref toleranceSolverA);
DA.GetData(9, ref toleranceSolverB);
DA.GetData(10, ref makeMesh);
DA.GetData(11, ref pp);
DA.GetData(12, ref plane);
if (!DA.GetData("Auto Run", ref autoRun))
return;
// If the user gives the pressures it must give also the load case
if (pressure != 0 && lc_pressure == null)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, "Give the Load Case to assign at the pressure forces");
return;
}
if (_run || autoRun)
{
try
{
_log = [];
Model ff_model = Common.Helpers.FormFindingHelper.AssemblyModel(model.Value, plane, solverType, gravity, pressure, makeMesh, out List<string> warnings);
if (solverType == 1)
{
try
{
List<Beam> fixedAxialForceBeam = [.. ff_model.Beams.Where(j => j.HasNLFixedAxialForce)];
List<Beam> fixedLengthBeam = [.. ff_model.Beams.Where(j => j.HasNLFixedLength)];
bool exit = true;
bool converged = true;
int count = 1;
ForceDensitySolver solver = new LinearForceDensitySolver(ff_model, maxIterNumber, tolerance, toleranceSolverA, toleranceSolverB, maxIterNumberSolver);
do
{
_log.Add($"STEP N. {count}");
exit = true;
converged = true;
Solver.ResultCodes ret = solver.SolveForceDensity();
if (fixedAxialForceBeam.Count == 0 && fixedLengthBeam.Count == 0)
{
_log.Add("No target force beam or target length beam in the model. Linear solution solved");
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, "No target force beam or target length beam in the model. Linear solution solved.");
_log.AddRange(solver.Log);
exit = true;
}
else
{
var fixedLengthError = new List<double>();
var fixedForceError = new List<double>();
for (int i = 0; i < fixedLengthBeam.Count; i++)
{
Beam bb = fixedLengthBeam[i];
double error = Math.Abs((bb.Length - bb.NLFixedLength) / bb.NLFixedLength);
fixedLengthError.Add(error);
if (error > tolerance)
{
exit = false;
converged = false;
bb.TargetForceDensity = bb.TargetForceDensity * bb.Length / bb.NLFixedLength;
}
}
for (int i = 0; i < fixedAxialForceBeam.Count; i++)
{
Beam bb = fixedAxialForceBeam[i];
double error = Math.Abs((bb.AxialForce - bb.TargetAxialForce) / bb.TargetAxialForce);
fixedForceError.Add(error);
if (error > tolerance)
{
exit = false;
converged = false;
bb.TargetForceDensity = bb.TargetForceDensity * bb.TargetAxialForce / bb.AxialForce;
}
}
count++;
_log.AddRange(solver.Log);
double a = fixedLengthError.Count == 0 ? 0 : Math.Round(fixedLengthError.Average(), 5);
double b = fixedLengthError.Count == 0 ? 0 : Math.Round(fixedLengthError.Max(), 5);
double c = fixedForceError.Count == 0 ? 0 : Math.Round(fixedForceError.Average(), 5);
double d = fixedForceError.Count == 0 ? 0 : Math.Round(fixedForceError.Max(), 5);
_log.Add($"RESULTS: Average Target Length relative error: {a}");
_log.Add($"RESULTS: Max Target Length relative error: {b}");
_log.Add($"RESULTS: Average Target Force relative error: {c}");
_log.Add($"RESULTS: Max Target Force relative error: {d}");
if (count > maxIterNumber)
{
_log.Add("Not converged. The maximum number of iteration has been reached.");
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, "Not converged. The maximum number of iteration has been reached.");
exit = true;
converged = false;
}
}
} while (!exit);
if (converged)
_log.Add("TARGET FORCE AND TARGET LENGTH CONVERGED");
}
catch (Exception e)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, $"{e.Message}");
}
}
else
{
try
{
ForceDensitySolver solver;
if (solverType == 0)
solver = new LinearForceDensitySolver(ff_model, maxIterNumber, tolerance, toleranceSolverA, toleranceSolverB, maxIterNumberSolver);
else if (solverType == 2)
solver = new NonLinearForceDensitySolver(ff_model, maxIterNumber, tolerance, toleranceSolverA, toleranceSolverB, maxIterNumberSolver)
{ UseUnitCoefficients = true };
else if (solverType == 3)
solver = new NonLinearForceDensitySolver(ff_model, maxIterNumber, tolerance, toleranceSolverA, toleranceSolverB, maxIterNumberSolver)
{ UseUnitCoefficients = false };
else
return;
Solver.ResultCodes ret = solver.SolveForceDensity();
_log = solver.Log;
if (ret == Solver.ResultCodes.NotconvergedAtMaxIterationNumber)
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, "Not converged. The maximum number of iteration has been reached.");
}
catch (Exception e)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, $"{e.Message}");
}
}
// Create a deep copy of the model where to add the results
_model = (GH_Model)model.Duplicate();
// Update the model with the FF results
NodeModel[] nodeArray = [.. _model.Value.Nodes];
for (int i = 0; i < nodeArray.Length; i++)
{
NodeModel node = nodeArray[i];
Node ff_node = ff_model.Nodes.ElementAt(i);
node.Position = new Point3d(ff_node.X, ff_node.Y, ff_node.Z);
var p = new Vector3d(ff_node.PressureX, ff_node.PressureY, ff_node.PressureZ);
if (p.Length != 0)
node.Loads.Add(new NodeForceLoadModel((LoadCaseModel)lc_pressure.Value, p));
}
BeamModel[] beamsModel = [.. _model.Value.Beams];
for (int i = 0; i < beamsModel.Length; i++)
{
BeamModel beam = beamsModel[i];
Beam ff_beam = ff_model.Beams.ElementAt(i);
beam.PointFrom = new Point3d(ff_beam.Node1.X, ff_beam.Node1.Y, ff_beam.Node1.Z);
beam.PointTo = new Point3d(ff_beam.Node2.X, ff_beam.Node2.Y, ff_beam.Node2.Z);
ForceDensityModel fd = (ForceDensityModel)beam.Loads.Where(l => l.GetType() == typeof(ForceDensityModel)).First();
if (fd != null)
{
LoadCaseModel lc = fd.LoadCase; // Use the same Load Case of the force density
for (int k = 0; k < beam.Loads.Count; k++)
{
if (beam.Loads[k].GetType() == typeof(BeamPreLoadModel))
{
beam.Loads.Remove(beam.Loads[k]);
k--;
}
}
beam.Loads.Add(new BeamPreLoadModel(lc, ff_beam.AxialForce, BeamPreLoadModel.LoadType.Tension));
}
}
_mesh = new GH_Mesh();
if (makeMesh) // Create the mesh if required
{
// Generate the plates
try
{
if (ff_model.Plates == null || ff_model.Plates.Count == 0)
ff_model.GeneratePlatesFromBeams();
}
catch (Exception)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, "Fail to generate the mesh");
}
var mesh = new Mesh();
var idRhinoIndexMap = new Dictionary<int, int>();
for (int i = 0; i < ff_model.Nodes.Count; i++)
{
Node ff_node = ff_model.Nodes.ElementAt(i);
int tag = mesh.Vertices.Add(ff_node.X, ff_node.Y, ff_node.Z);
idRhinoIndexMap.Add(ff_node.Id, tag);
}
for (int i = 0; i < ff_model.Plates.Count; i++)
{
Plate ff_plate = ff_model.Plates.ElementAt(i);
int i1 = idRhinoIndexMap[ff_plate.Node1.Id];
int i2 = idRhinoIndexMap[ff_plate.Node2.Id];
int i3 = idRhinoIndexMap[ff_plate.Node3.Id];
int i4 = -1;
if (ff_plate.PlateType == Plate.PlateTypes.Quad4)
i4 = idRhinoIndexMap[ff_plate.Node4.Id];
if (ff_plate.PlateType == Plate.PlateTypes.Quad4)
mesh.Faces.AddFace(i1, i2, i3, i4);
else if (ff_plate.PlateType == Plate.PlateTypes.Tri3)
mesh.Faces.AddFace(i1, i2, i3);
}
mesh.Compact();
mesh.RebuildNormals();
_mesh = new GH_Mesh(mesh);
//_model.Value.Plates.Clear();
if (_model.Value.Plates.Count > 0)
{
for (int i = 0; i < ff_model.Plates.Count; i++)
{
Plate ff_plate = ff_model.Plates.ElementAt(i);
PlateModel femmPlate = _model.Value.Plates.ElementAt(i);
if (int.TryParse(femmPlate.PlateId, out int femmPlateId))
{
if (femmPlateId == ff_plate.Id)
{
Point3d[] points = null;
if (ff_plate.PlateType == Plate.PlateTypes.Quad4)
{
points =
[
new Point3d(ff_plate.Node1.X, ff_plate.Node1.Y, ff_plate.Node1.Z),
new Point3d(ff_plate.Node2.X, ff_plate.Node2.Y, ff_plate.Node2.Z),
new Point3d(ff_plate.Node3.X, ff_plate.Node3.Y, ff_plate.Node3.Z),
new Point3d(ff_plate.Node4.X, ff_plate.Node4.Y, ff_plate.Node4.Z)
];
}
else
{
points =
[
new Point3d(ff_plate.Node1.X, ff_plate.Node1.Y, ff_plate.Node1.Z),
new Point3d(ff_plate.Node2.X, ff_plate.Node2.Y, ff_plate.Node2.Z),
new Point3d(ff_plate.Node3.X, ff_plate.Node3.Y, ff_plate.Node3.Z)
];
}
femmPlate.Points = new System.Collections.ObjectModel.ObservableCollection<Point3d>(points);
femmPlate.CreateShapes();
var gH_Plate = new GH_Plate(femmPlate);
_model.Guids.Remove(femmPlate);
_model.Guids.Add(femmPlate, gH_Plate.Guid);
//var plate = new PlateModel(points, pp.Value)
//{
// PlateId = ff_plate.Id.ToString()
//};
//plates[i] = plate;
//var gH_Plate = new GH_Plate(plate);
//_model.AddGuid(plate, gH_Plate.Guid);
}
//_model.Value.AddElements(plates);
}
}
}
else
{
if (pp != null) // Add the plates to the model if it's required
{
PlateModel[] plates = new PlateModel[ff_model.Plates.Count];
for (int i = 0; i < ff_model.Plates.Count; i++)
{
Plate ff_plate = ff_model.Plates.ElementAt(i);
Point3d[] points = null;
if (ff_plate.PlateType == Plate.PlateTypes.Quad4)
{
points =
[
new Point3d(ff_plate.Node1.X, ff_plate.Node1.Y, ff_plate.Node1.Z),
new Point3d(ff_plate.Node2.X, ff_plate.Node2.Y, ff_plate.Node2.Z),
new Point3d(ff_plate.Node3.X, ff_plate.Node3.Y, ff_plate.Node3.Z),
new Point3d(ff_plate.Node4.X, ff_plate.Node4.Y, ff_plate.Node4.Z)
];
}
else
{
points =
[
new Point3d(ff_plate.Node1.X, ff_plate.Node1.Y, ff_plate.Node1.Z),
new Point3d(ff_plate.Node2.X, ff_plate.Node2.Y, ff_plate.Node2.Z),
new Point3d(ff_plate.Node3.X, ff_plate.Node3.Y, ff_plate.Node3.Z)
];
}
if (pp.Value is PlatePropertyModel platePropertyModel)
{
var plate = new PlateModel(points, platePropertyModel)
{
PlateId = ff_plate.Id.ToString()
};
plates[i] = plate;
var gH_Plate = new GH_Plate(plate);
_model.Guids.Add(plate, gH_Plate.Guid);
}
}
_model.Value.AddElements(plates);
}
}
}
Message = "Done";
}
catch (Exception e)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, e.Message);
Message = "Error";
}
_runCount--;
if (_runCount == 0)
_run = false;
}
DA.SetData(0, _model);
DA.SetData(1, _mesh);
DA.SetDataList(2, _log);
}
protected override void BeforeSolveInstance()
{
_runCount = Params.Input[0].VolatileDataCount;
base.BeforeSolveInstance();
}
/// <summary>
/// Provides an Icon for the component.
/// </summary>
protected override System.Drawing.Bitmap Icon => Properties.Resources.FormFindIcon;
/// <summary>
/// Gets the unique ID for this component. Do not change this ID after release.
/// </summary>
public override Guid ComponentGuid => new("83be34af-f862-457b-8c03-109fadcd679e");
public override GH_Exposure Exposure => GH_Exposure.primary;
}
}