Implementation of new element and material for OpenSees in .NET environment (C# and MATLAB)


var staticLoading = false;
var theDomain = new OpenSees.Components.DomainWrapper();
var node1 = new OpenSees.Components.NodeWrapper(1, 2, 0, 0);
var node2 = new OpenSees.Components.NodeWrapper(2, 2, 144.0, 0);
var node3 = new OpenSees.Components.NodeWrapper(3, 2, 168.0, 0);
var node4 = new OpenSees.Components.NodeWrapper(4, 2, 72.0, 96.0);

theDomain.AddNode(new OpenSees.Components.NodeWrapper[] { node1, node2, node3, node4 });

var inclDamping = true;
var e = 3000;
var rho = 0.005;

//var theMaterial = new OpenSees.Materials.Uniaxials.ElasticMaterialWrapper(1, e, 0);
var theMaterial = new ElasticExternalUniaxialMaterial(1, e, 0, rho);

var truss1 = new OpenSees.Elements.TrussWrapper(1, 2, 1, 4, theMaterial, 10.0, 10 * rho, inclDamping ? 1 : 0, 0);
var truss2 = new OpenSees.Elements.TrussWrapper(2, 2, 2, 4, theMaterial, 5.0, 5 * rho, inclDamping ? 1 : 0, 0);
var truss3 = new OpenSees.Elements.TrussWrapper(3, 2, 3, 4, theMaterial, 5.0, 5 * rho, inclDamping ? 1 : 0, 0);


//var truss1 = new TestExternals.Models.Truss2d_2ExternalElement(1, 1, 4, 10.0, theMaterial, inclDamping);
//var truss2 = new TestExternals.Models.Truss2d_2ExternalElement(2, 2, 4, 5.0, theMaterial, inclDamping);
//var truss3 = new TestExternals.Models.Truss2d_2ExternalElement(3, 3, 4, 5.0, theMaterial, inclDamping);


theDomain.AddElement(truss1);
theDomain.AddElement(truss2);
theDomain.AddElement(truss3);

theDomain.SetMass(4, new OpenSees.MatrixWrapper(new double[,] { { 1, 0 }, { 0, 2 } }));


var sp1 = new OpenSees.Components.Constraints.SP_ConstraintWrapper(1, 0, 0.0, true);
var sp2 = new OpenSees.Components.Constraints.SP_ConstraintWrapper(1, 1, 0.0, true);
var sp3 = new OpenSees.Components.Constraints.SP_ConstraintWrapper(2, 0, 0.0, true);
var sp4 = new OpenSees.Components.Constraints.SP_ConstraintWrapper(2, 1, 0.0, true);
var sp5 = new OpenSees.Components.Constraints.SP_ConstraintWrapper(3, 0, 0.0, true);
var sp6 = new OpenSees.Components.Constraints.SP_ConstraintWrapper(3, 1, 0.0, true);
theDomain.AddSP_Constraint(new OpenSees.Components.Constraints.SP_ConstraintWrapper[] { sp1, sp2, sp4, sp5, sp6, sp3 });

if (staticLoading){
    var theSeries = new OpenSees.Components.Timeseries.LinearSeriesWrapper();
    var theLoadPattern = new OpenSees.Components.LoadPatterns.LoadPatternWrapper(1);
    theLoadPattern.SetTimeSeries(theSeries);
    theDomain.AddLoadPattern(theLoadPattern);

    var theLoadValues = new OpenSees.VectorWrapper(2);
    theLoadValues[0] = 100;
    theLoadValues[1] = -50;

    var nodalLoad = new OpenSees.Components.Loads.NodalLoadWrapper(1, 2, theLoadValues, false);
    theDomain.AddNodalLoad(nodalLoad, 1);
}
else
{

    var accelSeries = new OpenSees.Components.Timeseries.TriangleSeriesWrapper(0,4,2,0,1000,0);
    var groundMotion = new OpenSees.Components.GroundMotions.GroundMotionWrapper(null, null, accelSeries, null, 0.01, 1);
    var uniformExcitationPattern = new OpenSees.Components.LoadPatterns.UniformExcitationWrapper(2, groundMotion, 0, 0, 1);
    theDomain.AddLoadPattern(uniformExcitationPattern);
    theDomain.SetRayleighDampingFactors(0.1, 0.2, 0.3, 0.15);
}

{
    var opsstream = new OpenSees.Handlers.DataFileStreamWrapper(System.Environment.CurrentDirectory + $"\\FCol.out");
    var recorder = new OpenSees.Recorders.ElementRecorderWrapper(
        new OpenSees.IDWrapper(new int[] { 1 }),
        new string[] { "globalForce" },
        true,
        theDomain,
        opsstream
    );
    theDomain.AddRecorder(recorder);
}

var theModel = new OpenSees.AnalysisModelWrapper();
var theSolnAlgo = new OpenSees.Algorithms.NewtonRaphsonWrapper();
//var theIntegrator = new OpenSees.Integrators.Static.LoadControlWrapper(0.5, 1, 0.5, 0.5);
var theIntegrator = new OpenSees.Integrators.Transient.NewmarkWrapper(0.5, 0.25);
var theHandler = new OpenSees.Handlers.PlainHandlerWrapper();
var theRCM = new OpenSees.GraphNumberers.RCMWrapper(false);
var theNumberer = new OpenSees.Numberers.DOF_NumbererWrapper(theRCM);
var theSolver = new OpenSees.Systems.Linears.BandSPDLinLapackSolverWrapper();
var theSOE = new OpenSees.Systems.Linears.BandSPDLinSOEWrapper(theSolver);
var theTest = new OpenSees.ConvergenceTests.CTestNormDispIncrWrapper(1e-8, 6, 0, 2, 1.0e10);
var theAnalysis = new OpenSees.Analysis.DirectIntegrationAnalysisWrapper(
    theDomain,
    theHandler,
    theNumberer,
    theModel,
    theSolnAlgo,
    theSOE,
    theIntegrator,
    theTest);

double dt = 0.02;
double t = 10.0;
int nstep = (int)(t/dt);

for (var i = 0; i < nstep; i++)
{
    var ret = theAnalysis.Analyze(1, dt);
    Console.WriteLine(dt * i + dt);
    if (ret != 0)
        return;
}


var commitDisp = node4.GetCommitDisp();
Console.WriteLine(string.Join(",", commitDisp));
Console.ReadKey();

                        

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;


namespace TestExternals.Models
{
    interface IExternalUniaxialMaterial
    {
        OpenSees.Materials.Uniaxials.MDSetTrialStrain UMSetTrialStrain { get; }
        OpenSees.Materials.Uniaxials.MDGetStress UMGetStress { get; }
        OpenSees.Materials.Uniaxials.MDGetTangent UMGetTangent { get; }
        OpenSees.Materials.Uniaxials.MDGetInitialTangent UMGetInitialTangent { get; }
        OpenSees.Materials.Uniaxials.MDGetStrain UMGetStrain { get; }
        OpenSees.Materials.Uniaxials.MDCommitState UMCommitState { get; }
        OpenSees.Materials.Uniaxials.MDRevertToLastCommit UMRevertToLastCommit { get; }
        OpenSees.Materials.Uniaxials.MDRevertToStart UMRevertToStart { get; }
        OpenSees.Materials.Uniaxials.MDGetCopy UMGetCopy { get; }
        OpenSees.Materials.Uniaxials.MDPrint UMPrint { get; }
        OpenSees.Materials.Uniaxials.MDGetDampTangent UMDGetDampTangent { get; }
        OpenSees.Materials.Uniaxials.MDGetRho UMDGetRho { get; }
        OpenSees.Materials.Uniaxials.MDGetStrainRate UMDGetStrainRate { get; }

        void SetLinks();
    }

    //[System.Runtime.InteropServices.ComVisible(true)]
    public class ElasticExternalUniaxialMaterial : 
        OpenSees.Materials.Uniaxials.ExternalUniaxialMaterialWrapper, IExternalUniaxialMaterial
    {
        private int tag { get; set; }
        private double e { get; set; }
        private double eta { get; set; }
        private double rho { get; set; }
        private double e_commit { get; set; }
        private double e_trial { get; set; }
        private double s_commit { get; set; }
        private double s_trial { get; set; }
        public ElasticExternalUniaxialMaterial(int tag, double e, double eta, double rho)
            : base(tag)
        {
            this.tag = tag;
            this.e = e;
            this.eta = eta;
            this.rho = rho;
            e_commit = e_trial = s_commit = s_trial = 0;
            this.SetLinks();
        }

        public OpenSees.Materials.Uniaxials.MDSetTrialStrain UMSetTrialStrain => setTrialStrain;

        private int setTrialStrain(double strain, double strainRate)
        {
            e_trial = strain;
            s_trial = e * strain;
            return 0;
        }

        public OpenSees.Materials.Uniaxials.MDGetStress UMGetStress => getStress;
        private double getStress()
        {
            return s_trial;
        }

        public OpenSees.Materials.Uniaxials.MDGetTangent UMGetTangent => getTangent;
        private double getTangent()
        {
            return e;
        }

        public OpenSees.Materials.Uniaxials.MDGetInitialTangent UMGetInitialTangent => getInitialTangent;
        private double getInitialTangent()
        {
            return e;
        }

        public OpenSees.Materials.Uniaxials.MDGetStrain UMGetStrain => getStrain;
        private double getStrain()
        {
            return e_trial;
        }

        public OpenSees.Materials.Uniaxials.MDCommitState UMCommitState => commitState;
        private int commitState()
        {
            e_commit = e_trial;
            s_commit = e_trial;
            return 0;
        }

        public OpenSees.Materials.Uniaxials.MDRevertToLastCommit UMRevertToLastCommit => revertToLastCommit;
        private int revertToLastCommit()
        {
            e_trial = e_commit;
            e_trial = s_commit;
            return 0;
        }

        public OpenSees.Materials.Uniaxials.MDRevertToStart UMRevertToStart => revertToStart;
        private int revertToStart()
        {
            e_commit = e_trial = s_commit = s_trial = 0;
            return 0;
        }

        public OpenSees.Materials.Uniaxials.MDGetCopy UMGetCopy => getCopy;
        private OpenSees.Materials.Uniaxials.UniaxialMaterialWrapper getCopy()
        {
            var mat = new ElasticExternalUniaxialMaterial(tag, e, eta, rho);
            mat.SetLinks();
            return mat;
        }


        public void SetLinks()
        {
            this.SetLinks(this.UMSetTrialStrain, this.UMGetStress, this.UMGetTangent, this.UMGetInitialTangent,
                this.UMDGetDampTangent, this.UMGetStrain, this.UMDGetStrainRate,
                this.UMDGetRho, this.UMCommitState, this.UMRevertToLastCommit,
                 this.UMRevertToStart, this.UMGetCopy, this.UMPrint);
        }

        public OpenSees.Materials.Uniaxials.MDPrint UMPrint => print;

        private void print()
        {
            Console.WriteLine($"ElasticExternalUniaxialMaterial :");
            Console.WriteLine($"Elastic Modulus :{e}");
            Console.WriteLine($"Current Stress t/c :{s_trial} {s_commit}");
            Console.WriteLine($"Current Strain t/c :{e_trial} {e_commit}");
        }

        public OpenSees.Materials.Uniaxials.MDGetDampTangent UMDGetDampTangent => () => { return eta; };

        public OpenSees.Materials.Uniaxials.MDGetRho UMDGetRho => () => { return rho; };

        public OpenSees.Materials.Uniaxials.MDGetStrainRate UMDGetStrainRate => () => { return 0; };

    }
}

                        

using System;
using System.Linq;
using OpenSees.Elements;
using OpenSees.Materials.Uniaxials;

namespace TestExternals.Models
{

    public class Truss2d_2ExternalElement : ExternalElementWrapper, IExternalElement
    {
        int Node1 { get; set; }
        int Node2 { get; set; }
        UniaxialMaterialWrapper theMaterial;
        double a { get; set; }
        double c, s, length;
        bool doRayleighDamping { get; set; }
        OpenSees.Components.DomainWrapper theDomain { get; set; }
        OpenSees.Components.NodeWrapper ptrNode1;
        OpenSees.Components.NodeWrapper ptrNode2;
        static OpenSees.MatrixWrapper theMatrix;
        static OpenSees.VectorWrapper theVector;
        OpenSees.VectorWrapper theLoad;
        static int dimension = 2;
        public Truss2d_2ExternalElement(int tag, int node1, int node2, double a, UniaxialMaterialWrapper theMaterial, bool doRayleighDamping = false) : base(tag)
        {
            Node1 = node1;
            Node2 = node2;
            this.theMaterial = theMaterial.GetCopy();
            this.a = a;
            this.doRayleighDamping = doRayleighDamping;

            if (theMatrix == null)
                theMatrix = new OpenSees.MatrixWrapper(4, 4);
            if (theVector == null)
                theVector = new OpenSees.VectorWrapper(4);
            if (theLoad == null)
                theLoad = new OpenSees.VectorWrapper(4);

            SetLinks(DEle_GetClassType, DEle_GetNumExternalNodes, DEle_GetNodePtrs, DEle_GetExternalNodes,
                DEle_GetNumDOF, DEle_SetDomain, DEle_CommitState, DEle_RevertToLastCommit, DEle_RevertToStart, DEle_Update,
                DEle_GetTangentStiff, DEle_GetInitialStiff, DEle_GetDamp, DEle_GetMass, DEle_ZeroLoad, DEle_AddLoad, DEle_AddInertiaLoadToUnbalance,
                DEle_GetResistingForce, DEle_GetResistingForceIncInertia, DEle_Print, DEle_SetResponse, DEle_GetResponse);
        }

        public DEle_GetClassType DEle_GetClassType => () =>
        {
            return "Truss2dExternalElement";
        };

        public DEle_GetNumExternalNodes DEle_GetNumExternalNodes => () =>
        {
            return 2;
        };

        public DEle_GetExternalNodes DEle_GetExternalNodes => () =>
        {
            return new OpenSees.IDWrapper(new int[] { Node1, Node2 });
        };

        public DEle_GetNodePtrs DEle_GetNodePtrs => () =>
        {
            return new OpenSees.Components.NodeWrapper[] { ptrNode1, ptrNode2 };
        };

        public DEle_GetNumDOF DEle_GetNumDOF => () => { return 4; };

        public DEle_SetDomain DEle_SetDomain => (OpenSees.Components.DomainWrapper theDomain) =>
        {
            this.theDomain = theDomain;
            ptrNode1 = theDomain.GetNode(Node1);
            ptrNode2 = theDomain.GetNode(Node2);

            if (ptrNode1.GetNumberDOF() != ptrNode2.GetNumberDOF() || ptrNode1.GetNumberDOF() != 2)
            {
                Console.WriteLine($"WARNING Truss2dExternalElement::setDomain(): nodes {Node1} and {Node2}" +
                $" have differing dof at ends for truss2d {this.GetTag()} ");
                System.Environment.Exit(-1);
            }

            var crd1 = ptrNode1.GetCrds();
            var crd2 = ptrNode2.GetCrds();


            length = Math.Sqrt((crd2[0] - crd1[0]) * (crd2[0] - crd1[0]) + (crd2[1] - crd1[1]) * (crd2[1] - crd1[1]));
            s = (crd2[1] - crd1[1]) / length;
            c = (crd2[0] - crd1[0]) / length;
        };

        public DEle_CommitState DEle_CommitState => () =>
        {
            return theMaterial.CommitState();
        };

        public DEle_RevertToLastCommit DEle_RevertToLastCommit => () =>
        {
            return theMaterial.RevertToLastCommit();
        };

        public DEle_RevertToStart DEle_RevertToStart => () =>
        {
            return theMaterial.RevertToStart();
        };

        public DEle_Update DEle_Update => () =>
        {
            return 0;
        };

        public DEle_GetTangentStiff DEle_GetTangentStiff => () =>
        {
            return getTangent(theMaterial.GetTangent());
        };

        public DEle_GetInitialStiff DEle_GetInitialStiff => () =>
        {
            return getTangent(theMaterial.GetInitialTangent());
        };

        private OpenSees.MatrixWrapper getTangent(double e)
        {
            var st = e * a / length;
            theMatrix[0, 0] = c * c * st;
            theMatrix[0, 1] = s * c * st;
            theMatrix[0, 2] = -c * c * st;
            theMatrix[0, 3] = -s * c * st;

            theMatrix[1, 1] = s * s * st;
            theMatrix[1, 2] = -s * c * st;
            theMatrix[1, 3] = -s * s * st;

            theMatrix[2, 2] = c * c * st;
            theMatrix[2, 3] = c * s * st;

            theMatrix[3, 3] = s * s * st;

            for (var i = 0; i < 4; i++)
                for (var j = 0; j < i; j++)
                    theMatrix[i, j] = theMatrix[j, i];

            return theMatrix;
        }

        public DEle_AddLoad DEle_AddLoad => (e, f) =>
        {
            return 0;
        };

        public DEle_AddInertiaLoadToUnbalance DEle_AddInertiaLoadToUnbalance => (OpenSees.VectorWrapper accel) =>
        {
            double rho = theMaterial.GetRho() * a;
            if (length == 0 || rho == 0) return 0;
            var raccel1 = ptrNode1.GetRVVector(accel);
            var raccel2 = ptrNode2.GetRVVector(accel);
            int nodalDOF = GetNumDOF() / 2;

            // lumped mass
            double m = 0.5 * rho * length;
            for (int i = 0; i < dimension; i++)
            {
                theLoad[i] -= m * raccel1[i];
                theLoad[i + nodalDOF] -= m * raccel2[i];
            }

            return 0;
        };

        public DEle_GetResistingForce DEle_GetResistingForce => () =>
        {
            return getResistingForce();
        };

        private OpenSees.VectorWrapper getResistingForce()
        {
            var st = theMaterial.GetTangent() * a / length;
            var disp1 = ptrNode1.GetTrialDisp();
            var disp2 = ptrNode2.GetTrialDisp();
            double t = st * ((disp2[0] - disp1[0]) * c + (disp2[1] - disp1[1]) * s);
            theVector[0] = -t * c;
            theVector[1] = -t * s;
            theVector[2] = t * c;
            theVector[3] = t * s;

            return theVector;
        }

        public DEle_GetResistingForceIncInertia DEle_GetResistingForceIncInertia => () =>
        {
            double rho = theMaterial.GetRho() * a;
            getResistingForce(); // set resisting force to theVector
            
            theVector.Add(1, theLoad, -1);
            
            if (length != 0 && rho != 0)
            {
                var accel1 = ptrNode1.GetTrialAccel();
                var accel2 = ptrNode2.GetTrialAccel();
                int numDOF2 = GetNumDOF() / 2;
                double m = 0.5 * rho * length;
                for (int i = 0; i < dimension; i++)
                {
                    theVector[i] += m * accel1[i];
                    theVector[i + numDOF2] += m * accel2[i];
                }
            }
           
            if (doRayleighDamping)
                theVector.Add(1, this.GetRayleighDampingForcesVectorWrapper(), 1);
            
            return theVector;
        };

        public DEle_Print DEle_Print => (int flag) =>
        {
            Console.WriteLine($"{GetClassType()} ({GetTag()}) Nodes {Node1},{Node2}");
        };

        public DEle_GetDamp DEle_GetDamp => () =>
        {
            theMatrix.Zero();
            
            if (doRayleighDamping)
                theMatrix.Add(0, this.GetBaseElementDamp(), 1);

            double eta = theMaterial.GetDampTangent();
            if (eta == 0)
                return theMatrix;
            int numDOF2 = GetNumDOF() / 2;
            double temp;
            double etaAoverL = eta * a / length;
            var cosX = new double[] { c, s };
            for (int i = 0; i < dimension; i++)
            {
                for (int j = 0; j < dimension; j++)
                {
                    temp = cosX[i] * cosX[j] * etaAoverL;
                    theMatrix[i, j] += temp;
                    theMatrix[i + numDOF2, j] += -temp;
                    theMatrix[i, j + numDOF2] += -temp;
                    theMatrix[i + numDOF2, j + numDOF2] += temp;
                }
            }

            return theMatrix;
        };

        public DEle_GetMass DEle_GetMass => () =>
        {
            theMatrix.Zero();
            double rho = theMaterial.GetRho() * a;
            if (rho == 0 || length == 0)
                return theMatrix;
            double m = 0.5 * rho * length;
            int numDOF2 = this.GetNumDOF() / 2;

            for (int i = 0; i < dimension; i++)
            {
                theMatrix[i, i] = m;
                theMatrix[i + numDOF2, i + numDOF2] = m;
            }

            return theMatrix;
        };

        public DEle_ZeroLoad DEle_ZeroLoad => () =>
        {
            theLoad.Zero();
        };

        enum Truss2dResponseType { globalForce = 1, localForce = 2, deformation = 3, material = 4 }
        public DEle_SetResponse DEle_SetResponse => (string[] argv, OpenSees.Handlers.OPS_StreamWrapper output) =>
        {

            OpenSees.Recorders.ResponseWrapper response = null;
            output.Tag("ElementOutput");
            output.Attr("eleType", "Truss2d");
            output.Attr("eleTag", this.GetTag());
            output.Attr("node1", this.Node1);
            output.Attr("node2", this.Node2);

            if (argv[0] == "globalForce")
            {
                output.Tag("ResponseType", "P1_1");
                output.Tag("ResponseType", "P1_2");
                output.Tag("ResponseType", "P2_1");
                output.Tag("ResponseType", "P2_2");

                response = new OpenSees.Recorders.ElementResponseWrapper(this, (int)Truss2dResponseType.globalForce, new OpenSees.VectorWrapper(4));
            }
            else if (argv[0] == "localForce")
            {
                response = new OpenSees.Recorders.ElementResponseWrapper(this, (int)Truss2dResponseType.localForce, 0.0);
            }
            else if (argv[0] == "deformation")
            {
                response = new OpenSees.Recorders.ElementResponseWrapper(this, (int)Truss2dResponseType.deformation, 0.0);
            }
            else if (argv[0] == "material" && argv.Length > 1)
            {
                var _argv = argv.ToList();
                _argv.RemoveAt(0);
                response = theMaterial.SetResponse(_argv.ToArray(), output);
            }

            output.EndTag();
            return response;
        };

        public DEle_GetResponse DEle_GetResponse => (int id, OpenSees.Recorders.InformationWrapper info)=>
        {
            
            switch ((Truss2dResponseType)id)
            {
                case Truss2dResponseType.globalForce:
                    return info.SetVector(getResistingForce());
                case Truss2dResponseType.localForce:
                    var stress = theMaterial.GetStress();
                    return info.SetDouble(stress * this.a);
                case Truss2dResponseType.deformation:
                    var strian = theMaterial.GetStrain();
                    return info.SetDouble(strian * length);
                case Truss2dResponseType.material:
                    return theMaterial.GetResponse(id, info);
                default:
                    break;
            }
            return 0;
        };
    }
}
                            

B: MATLAB


function ExternalTest()
NET.addAssembly('System');
NET.addAssembly('System.IO');
NET.addAssembly(strcat(pwd,'\OpenSees.NET.X64.dll'));
staticLoading = false;
inclDamping = true;
inclDampingInt = 1;
rho = 0.005;
e = 3000;
theDomain = OpenSees.Components.DomainWrapper();
node1 = OpenSees.Components.NodeWrapper(1, 2, 0, 0);
node2 = OpenSees.Components.NodeWrapper(2, 2, 144.0, 0);
node3 = OpenSees.Components.NodeWrapper(3, 2, 168.0, 0);
node4 = OpenSees.Components.NodeWrapper(4, 2, 72.0, 96.0);

theDomain.AddNode(node1);
theDomain.AddNode(node2);
theDomain.AddNode(node3);
theDomain.AddNode(node4);

theMaterial = OpenSees.Materials.Uniaxials.ElasticMaterialWrapper(1, e, 0);
%theMaterial = ElasticExternalUniaxialMaterial(1,e,0,rho).material;

% truss1 = Truss2dExternalElement(1, 1, 4, 10.0, theMaterial, inclDamping).theElement;
% truss2 = Truss2dExternalElement(2, 2, 4, 5.0, theMaterial, inclDamping).theElement;
% truss3 = Truss2dExternalElement(3, 3, 4, 5.0, theMaterial, inclDamping).theElement;

truss1 = OpenSees.Elements.TrussWrapper(1, 2, 1, 4, theMaterial, 10.0, 10.0*rho, inclDampingInt, 0);
truss2 = OpenSees.Elements.TrussWrapper(2, 2, 2, 4, theMaterial, 5.0, 5.0*rho, inclDampingInt, 0);
truss3 = OpenSees.Elements.TrussWrapper(3, 2, 3, 4, theMaterial, 5.0, 5.0*rho, inclDampingInt, 0);

theDomain.AddElement(truss1);
theDomain.AddElement(truss2);
theDomain.AddElement(truss3);

theDomain.SetMass(4, OpenSees.MatrixWrapper([1,0;0,2]));

sp1 = OpenSees.Components.Constraints.SP_ConstraintWrapper(1, 0, 0.0, true);
sp2 = OpenSees.Components.Constraints.SP_ConstraintWrapper(1, 1, 0.0, true);
sp3 = OpenSees.Components.Constraints.SP_ConstraintWrapper(2, 0, 0.0, true);
sp4 = OpenSees.Components.Constraints.SP_ConstraintWrapper(2, 1, 0.0, true);
sp5 = OpenSees.Components.Constraints.SP_ConstraintWrapper(3, 0, 0.0, true);
sp6 = OpenSees.Components.Constraints.SP_ConstraintWrapper(3, 1, 0.0, true);

theDomain.AddSP_Constraint(sp1);
theDomain.AddSP_Constraint(sp2);
theDomain.AddSP_Constraint(sp3);
theDomain.AddSP_Constraint(sp4);
theDomain.AddSP_Constraint(sp5);
theDomain.AddSP_Constraint(sp6);

theModel = OpenSees.AnalysisModelWrapper();
theSolnAlgo = OpenSees.Algorithms.NewtonRaphsonWrapper();
theHandler = OpenSees.Handlers.PlainHandlerWrapper();
theRCM = OpenSees.GraphNumberers.RCMWrapper(false);
theNumberer = OpenSees.Numberers.DOF_NumbererWrapper(theRCM);
theSolver = OpenSees.Systems.Linears.BandSPDLinLapackSolverWrapper();
theSOE = OpenSees.Systems.Linears.BandSPDLinSOEWrapper(theSolver);
test = OpenSees.ConvergenceTests.CTestNormDispIncrWrapper(1e-8, 6, 0, 2, 1.0e10);

if staticLoading 
    % static
    theSeries = OpenSees.Components.Timeseries.LinearSeriesWrapper();
    theLoadPattern = OpenSees.Components.LoadPatterns.LoadPatternWrapper(1);
    theLoadPattern.SetTimeSeries(theSeries);
    theDomain.AddLoadPattern(theLoadPattern);
    theLoadValues = OpenSees.VectorWrapper([100,-50]);
    nodalLoad = OpenSees.Components.Loads.NodalLoadWrapper(1, 4, theLoadValues, false);
    theDomain.AddNodalLoad(nodalLoad, 1);
    theIntegrator = OpenSees.Integrators.Static.LoadControlWrapper(1.0, 1, 1.0, 1.0);
    
    theAnalysis = OpenSees.Analysis.StaticAnalysisWrapper(...
    theDomain,...
    theHandler,...
    theNumberer,...
    theModel,...
    theSolnAlgo,...
    theSOE,...
    theIntegrator,...
    test);

    numSteps = 1;
    theAnalysis.Analyze(numSteps);
else
    % dynamic
    accelSeries = OpenSees.Components.Timeseries.TriangleSeriesWrapper(0,4,2,0,1000,0);
    groundMotion = OpenSees.Components.GroundMotions.GroundMotionWrapper(...
        [], ...
        [], ...
        accelSeries, ...
        [], ...
        0.01, 1);
    uniformExcitationPattern = OpenSees.Components.LoadPatterns.UniformExcitationWrapper(2, ...
        groundMotion, 0, 0, 1);
    theDomain.AddLoadPattern(uniformExcitationPattern);
    theDomain.SetRayleighDampingFactors(0.1, 0.2, 0.3, 0.15);
    theIntegrator = OpenSees.Integrators.Transient.NewmarkWrapper(0.5, 0.25);
    theAnalysis = OpenSees.Analysis.DirectIntegrationAnalysisWrapper(...
    theDomain,...
    theHandler,...
    theNumberer,...
    theModel,...
    theSolnAlgo,...
    theSOE,...
    theIntegrator,...
    test);
    
    dt = 0.02;
    t = 10.0;
    numSteps = int32(t/dt);
    for i=1:numSteps
        theAnalysis.Analyze(1,dt);
    end
end
disp(double(node4.GetCommitDisp()));
end

classdef ElasticExternalUniaxialMaterial < handle
    
    properties
        tag,e,eps_trial,eps_commit,str_trial,str_commit;
        eta,rho;
        material;
    end
    
    methods
        function this = ElasticExternalUniaxialMaterial(tag,e,eta,rho)
            this.tag = tag;
            this.e = e;
            this.eta = eta;
            this.rho = rho;
            
            this.eps_trial=0;
            this.eps_commit=0;
            this.str_trial=0;
            this.str_commit=0;
            
            this.material = OpenSees.Materials.Uniaxials.ExternalUniaxialMaterialWrapper(this.tag);
            this.material.SetLinks(@this.SetTrialStrain,@this.GetStress, ...
                @this.GetTangent,@this.GetInitialTangent,@this.GetDampTangent,...
                @this.GetStrain,@this.GetStrainRate,@this.GetRho, ...
                @this.CommitState,@this.RevertToLastCommit,@this.RevertToStart, ...
                @this.GetCopy,@this.Print);
        end
        
        function r = SetTrialStrain(this,strain, ~)
            this.eps_trial = strain;
            this.str_trial = this.eps_trial*this.e;
            r = 0;
        end
        
        function r = GetStress(this)
            r = this.str_trial;
        end
        
        function r = GetTangent(this)
            r = this.e;
        end
        
        function r = GetInitialTangent(this)
            r = this.e;
        end
        
        function r = GetStrain(this)
            r = this.eps_trial;
        end
        
        function r = CommitState(this)
            this.eps_commit = this.eps_trial;
            this.str_commit = this.str_trial;
            r = 0;
        end
        
        function r = RevertToLastCommit(this)
            this.eps_trial = this.eps_commit;
            this.str_trial = this.str_commit;
            r = 0;
        end
        
        function r = RevertToStart(this)
            this.eps_trial = 0;
            this.eps_commit = 0;
            this.str_trial = 0;
            this.str_commit = 0;
            r = 0;
        end
        
        function r = GetCopy(this)
            copy = ElasticExternalUniaxialMaterial(this.tag,this.e,this.eta,this.rho);
            r = copy.material;
        end
        
        function Print(this)
            disp(strcat('ElasticExternalUniaxialMaterial Tag :',this.tag));
            disp(strcat('Elastic Modulus :',this.e));
            disp(strcat('Current Stress t/c :',this.str_trial,' ',this.str_commit));
            disp(strcat('Current Strain t/c :',this.eps_trial,' ',this.eps_commit));
        end
        
        function r = GetDampTangent(this)
            r= this.eta;
        end
        
        function r = GetRho(this)
            r= this.rho;
        end
        
        function r = GetStrainRate(~)
            r= 0;
        end
        
            
        
    end
end

classdef Truss2dExternalElement < handle
    properties (Access = private)
        tag;
        node1,node2;
        a;
        doRayleighDamping;
        c,s, length;
        
        theDomain;
        ptrNode1,ptrNode2;
        theMaterial;  
        theLoad;
    end
    
    properties (Access = public)     
        theElement;
    end
    
    methods (Access = private, Static = true)
        function val = theMatrix(newval)
            persistent currentval;
            if nargin >= 1
                currentval = newval;
            end
            val = currentval;
        end
        function val = theVector(newval)
            persistent currentval;
            if nargin >= 1
                currentval = newval;
            end
            val = currentval;
        end
     end
    methods (Access = private)
        function computeTangent(this, e)
            st = e*this.a/this.length;
            
            this.theMatrix.Set(0,0,this.c*this.c*st);
            this.theMatrix.Set(0,1,this.s*this.c*st);
            this.theMatrix.Set(0,2,-this.c*this.c*st);
            this.theMatrix.Set(0,3,-this.s*this.c*st);
            
            this.theMatrix.Set(1,1,this.s*this.s*st);
            this.theMatrix.Set(1,2,-this.s*this.c*st);
            this.theMatrix.Set(1,3,-this.s*this.c*st);
            
            this.theMatrix.Set(2,2,this.c*this.c*st);
            this.theMatrix.Set(2,3,this.s*this.c*st);
            
            this.theMatrix.Set(3,3,this.s*this.s*st);
            
            for i=0:3
                for j=0:i-1
                   this.theMatrix.Set(i,j,this.theMatrix.Get(j,i)); 
                end
            end
        end
        
        function computeResistingForce(this, e)
            st = e*this.a/this.length;
            d1 = this.ptrNode1.GetTrialDisp();
            d2 = this.ptrNode2.GetTrialDisp();
            t = st*((d2(1)-d1(1))*this.c+((d2(2)-d1(2))*this.s));
            this.theVector.Set(0,-t*this.c);
            this.theVector.Set(1,-t*this.s);
            this.theVector.Set(2,t*this.c);
            this.theVector.Set(3,t*this.s);
            
        end
    end
    methods
        function this = Truss2dExternalElement(tag,node1,node2,a,theMaterial,doRayleighDamping)
            this.tag = tag;
            this.node1 = node1;
            this.node2 = node2;
            this.theMaterial = theMaterial.GetCopy();
            this.a = a;
            this.doRayleighDamping = doRayleighDamping;
            
            
            % init this.theMatrix,this.this.theVector, theLoad
            this.theLoad = OpenSees.VectorWrapper(zeros(1,4));
            if isempty(this.theMatrix)
                this.theMatrix(OpenSees.MatrixWrapper(zeros(4,4)));
            end
            if isempty(this.theVector)
                this.theVector(OpenSees.VectorWrapper(zeros(1,4)));
            end
            
            this.theElement = OpenSees.Elements.ExternalElementWrapper(this.tag);
            this.theElement.SetLinks(@this.GetClassType,@this.GetNumExternalNodes,...
                @this.GetNodePtrs,@this.GetExternalNodes,...
                @this.GetNumDOF,@this.SetDomain,@this.CommitState, ...
                @this.RevertToLastCommit,@this.RevertToStart,@this.Update,...
                @this.GetTangentStiff,@this.GetInitialStiff,@this.GetDamp,...
                @this.GetMass,@this.ZeroLoad,@this.AddLoad,...
                @this.AddInertiaLoadToUnbalance,@this.GetResistingForce,@this.GetResistingForceIncInertia,...
                @this.Print,@this.SetResponse, @this.GetResponse);
        end
        
        function r = GetClassType(~)
            r = "Truss2dExternalElement";
        end
        
        function r = GetNumExternalNodes(~)
            r = 2;
        end
        
        function r = GetExternalNodes(this)
            r = OpenSees.IDWrapper([this.node1,this.node2]);
        end
        
        function r = GetNodePtrs(this)
            list = NET.createArray('OpenSees.Components.NodeWrapper',2);
            list(1)= this.ptrNode1;
            list(2)= this.ptrNode2;
            r = list;
        end
        
        function r = GetNumDOF(~)
            r = 4;
        end
        
        function SetDomain(this,theDomain)
            this.theDomain = theDomain;
            this.ptrNode1 = this.theDomain.GetNode(this.node1);
            this.ptrNode2 = this.theDomain.GetNode(this.node2);
            if (this.ptrNode1.GetNumberDOF() ~= 2) && (this.ptrNode2.GetNumberDOF() ~= 2)
                disp('WARNING Truss2dExternalElement::setDomain(): nodes dof not compatible');
            end 
            crd1 = this.ptrNode1.GetCrds();
            crd2 = this.ptrNode2.GetCrds();
            this.length = sqrt((crd2(1)-crd1(1))^2+(crd2(2)-crd1(2))^2);
            this.s = (crd2(2)-crd1(2))/this.length;
            this.c = (crd2(1)-crd1(1))/this.length;
        end
        
        function r = CommitState(this)
            r = this.theMaterial.CommitState();
        end
        
        function r = RevertToLastCommit(this)
            r = this.theMaterial.RevertToLastCommit();
        end
        
        function r = RevertToStart(this)
            r = this.theMaterial.RevertToStart();
        end
        
        function r = Update(~)
            r = 0;
        end
        
        function r = GetTangentStiff(this)
            e = this.theMaterial.GetTangent();
            this.computeTangent(e);
            r = this.theMatrix;
        end
        
        function r = GetInitialStiff(this)
            e0 = this.theMaterial.GetInitialTangent();
            this.computeTangent(e0);
            r = this.theMatrix;
        end
        
        function r = AddLoad(~,~,~) %ElementalLoadWrapper theLoad, double loadFactor
            r = 0;
        end
        
        function r = AddInertiaLoadToUnbalance(this, accel)
            rho = this.theMaterial.GetRho()*this.a;
            if this.length == 0 || rho == 0
                r = 0;
                return;
            end
            raccel1 = this.ptrNode1.GetRV(accel);
            raccel2 = this.ptrNode2.GetRV(accel);
            nodalDOF = this.theElement.GetNumDOF() / 2;
            m = 0.5 *rho * this.length;
            data = [0,0,0,0];
            for i=1:nodalDOF
                data(i) = -m*raccel1(i);
                data(i+nodalDOF) = -m*raccel2(i);
            end
            this.theLoad.SetData(data,1,1);
            r = 0;
        end
        
        function r = GetResistingForce(this)
            e = this.theMaterial.GetTangent();
            this.computeResistingForce(e);
            r = this.theVector;
        end
        
        function r = GetResistingForceIncInertia(this)
            e = this.theMaterial.GetTangent();
            this.computeResistingForce(e);
            this.theVector.Add(1,this.theLoad,-1);
            rho = this.theMaterial.GetRho()*this.a;
            if rho ~= 0 && this.length ~= 0
                accel1 = this.ptrNode1.GetTrialAccel();
                accel2 = this.ptrNode2.GetTrialAccel();
                numDOF2 = this.theElement.GetNumDOF() / 2;
                m = 0.5 * rho * this.length;
                data = zeros(4,1);
                for i = 1:2
                    data(i) = m * accel1(i);
                    data(i+numDOF2) = m * accel2(i);
                end
                this.theVector.SetData(data,1,1);
            end          
            if this.doRayleighDamping
                this.theVector.Add(1,this.theElement.GetRayleighDampingForcesVectorWrapper(),1);
            end
            r = this.theVector;
        end
        
        function Print(~) % flag
            
        end
        
        function r = GetDamp(this)
            this.theMatrix.Zero();
            if this.doRayleighDamping
                this.theMatrix.Add(0.0,this.theElement.GetBaseElementDamp(),1);
            end
            eta = this.theMaterial.GetDampTangent();
            if eta == 0
                r = this.theMatrix;
                return;
            end
            numDOF2 = this.theElement.GetNumDOF()/2;
            etaAoverL = eta * this.a / this.length;
            cosX = [this.c,this.s];
            data = zeros(4,4);
            for i=1:2
                for j=1:2
                    temp = cosX(i)*cosX(j)*etaAoverL;
                    data(i,j) = data(i,j) + temp;                                         
                    data(i+numDOF2,j) = data(i+numDOF2,j) - temp;
                    data(i,j+numDOF2) = data(i,j+numDOF2) - temp;
                    data(i+numDOF2,j+numDOF2) = data(i+numDOF2,j+numDOF2) + temp;
                end
            end
            this.theMatrix.SetData(data,1,1);            
            r= this.theMatrix;
        end
        
        function r = GetMass(this)
            this.theMatrix.Zero();
            rho = this.theMaterial.GetRho()*this.a;
            if rho == 0 && this.length == 0
                r = this.theMatrix;
                return;
            end
            m = 0.5 *rho * this.length;
            numDOF2 = this.theElement.GetNumDOF()/2;
            data = zeros(4,4);
            for i=1:2
               data(i,i) = m;
               data(i+ numDOF2,i+ numDOF2) = m;
            end
            this.theMatrix.SetData(data,1,1);
            r= this.theMatrix;
        end
        
        function ZeroLoad(this)
            this.theLoad.Zero();
        end
        
        function r = SetResponse(~,~,~) % this,(string[] argv, OpenSees.Handlers.OPS_StreamWrapper output)
            r = []; % ResponseWrapper
        end
        
        function r = GetResponse(~,~) % this,int id, OpenSees.Recorders.InformationWrapper info
            r = 0;
        end
    end
end