iainsproat / SharpFE

Lightweight, modular finite element solver for .net
BSD 3-Clause "New" or "Revised" License
37 stars 9 forks source link

Full3D model fail to solve. #1

Closed AndreasBak closed 11 years ago

AndreasBak commented 11 years ago

The following test throws a zero determinant exception.

    [Test]
    public void PortalFrameNegativeCoordinates()
    {
        FiniteElementModel model = new FiniteElementModel(ModelType.Full3D);
        FiniteElementNode node1 = model.NodeFactory.Create(-10,0,0); 
        model.ConstrainNode(node1, DegreeOfFreedom.X);
        model.ConstrainNode(node1, DegreeOfFreedom.Y);
        model.ConstrainNode(node1, DegreeOfFreedom.Z);
        model.ConstrainNode(node1, DegreeOfFreedom.XX);
        model.ConstrainNode(node1, DegreeOfFreedom.YY);
        model.ConstrainNode(node1, DegreeOfFreedom.ZZ);

        FiniteElementNode node2 = model.NodeFactory.Create(-10,0,10);

        FiniteElementNode node3 = model.NodeFactory.Create(0,0,14);
        FiniteElementNode node4 = model.NodeFactory.Create(10,0,10);
        FiniteElementNode node5 = model.NodeFactory.Create(10,0,0);
        model.ConstrainNode(node5, DegreeOfFreedom.X);
        model.ConstrainNode(node5, DegreeOfFreedom.Y);
        model.ConstrainNode(node5, DegreeOfFreedom.Z);

        model.ConstrainNode(node5, DegreeOfFreedom.XX);
        model.ConstrainNode(node5, DegreeOfFreedom.YY);
        model.ConstrainNode(node5, DegreeOfFreedom.ZZ);

        IMaterial material = new GenericElasticMaterial(100, 2000, 0.2, 1000);
        ICrossSection section = new SolidRectangle(0.5, 0.1);

        model.ElementFactory.CreateLinear3DBeam(node1, node2, material, section);
        model.ElementFactory.CreateLinear3DBeam(node2,node3,material,section);
        model.ElementFactory.CreateLinear3DBeam(node3,node4,material,section);
        model.ElementFactory.CreateLinear3DBeam(node4,node5,material,section);

        ForceVector force = model.ForceFactory.Create(10, 0, 0,0,0,0);
        model.ApplyForceToNode(force, node3);

        IFiniteElementSolver solver = new LinearSolver(model);
        FiniteElementResults results = solver.Solve();

        DisplacementVector displacement = results.GetDisplacement(node2);
        Assert.AreNotEqual(0.0, displacement.X);
        Assert.AreEqual(0.0, displacement.Y, 0.001);
    }
iainsproat commented 11 years ago

I've amended the example and created a similar example for the 1D beam.

It seems that the negative coordinates are not the problem. The 1D beam examples works OK, so the problem seems to be either: (a) some problem with not enough boundary conditions i.e. constraining nodes (b) implementation of the 3D beam stiffness matrix class

iainsproat commented 11 years ago

It works if you constrain either node2, node3 or node4 in the X-axis, however this modifies the intent of the model so is not a valid fix.

AndreasBak commented 11 years ago

I should think that a moment frame fully fixed at both supports (node1 and node5) should be constrained sufficiently.

AndreasBak commented 11 years ago

I'm having the same issue trying to implement the ConstantStrainTriangle in the Full3d model. I'll write a test later to show this issue. This however this suggest that the issue might not be in the beam stiffness matrix class.

AndreasBak commented 11 years ago

Great job.