Shell elements not positive definite exception

[rubsy92]

Hi, I'm trying to test the shell elements, but each time I run an analysis I get the "not positive definite" exception. I noticed that you always constrain all nodes in your examples, what's the reason for that? Correct me if I'm wrong but I would assume that you don't get any displacement when you constrain everything?

[epsi1on]

Constraining a node means set Node.Constraints = Constraints.Fixed. If set all nodes are constrainted as fix then there is no nodal displacement at all. Yes you are right.

The reason you do get not pos def error, is that there is at least a DoF on a node that is not fixed to ground nor connected to any element.

Which shell element you are using?

[rubsy92]

I used the full shell implementation (so all 3 components). I noticed that it does calculate when I change the solver type to ConjugateGradient (which you commented out in the code)... But there I noticed that the displacements are not correct (with a factor of my shell thickness). I'll try to look into it this weekend.

[epsi1on]

You mean use TriangleElement with behavior of FullThinShell? if you was not able to fix it, then pls serialize your model or a model with same error and send it here to check it out...,

[rubsy92]

Hi,

I used "TriangleElement" as the element type and " triElm.Behavior = TriangleElementBehaviours.Shell;" as the behaviour. The steps I took:

• Downloaded the latest version from github
• I copied your "Generate3DTriangleElementGrid" method and I replaced some things so that all nodes on the x=0 are fully constrained and that all nodes on the other side have a force in the positive x-direction. See below for the code - I'm sorry but I have no experience with serialisation.
• I made a model with the method mentioned above and only 1 row of elements, which comes down to a flat plane. Since the default solver was not working I changed the solver to "ConjugateGradient" and then it did ran, but with incorrect displacements... Since the difference with my benchmark software is exactly 1/(element thickness) I guess something goes wrong there. The displacements also don't change when I change the thickness so that's already a big pointer to a potential error in the element thickness...

########CODE####### MAIN METHOD: var nx = 6; var ny = 6; var nz = 1; var grd = StructureGenerator.Generate3DTriangleElementGridTest(nx, ny, nz); model = grd; grd.Solve_MPC();

METHOD FOR MODEL: public static Model Generate3DTriangleElementGridTest(int m, int n, int l) { var buf = new Model();

        var dx = 0.15;
        var dy = 0.15;
        var dz = 0.15;

        var nodes = new Node[m, n, l];

        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
            {
                for (int k = 0; k < l; k++)
                {
                    var pos = new Point(j * dx, i * dy, k * dz);
                    var nde = new Node() { Location = pos };
                    buf.Nodes.Add(nde);
                    nodes[i, j, k] = nde;
                }
            }
        }

        //elements parallel to XZ
        for (int j = 0; j < n; j++)
        {
            for (int i = 0; i < m - 1; i++)
            {
                for (int k = 0; k < l - 1; k++)
                {
                    {
                        var elm2 = new TriangleElement();
                        elm2.Nodes[0] = nodes[i, j, k];
                        elm2.Nodes[1] = nodes[i, j, k + 1];
                        elm2.Nodes[2] = nodes[i + 1, j, k];

                        buf.Elements.Add(elm2);
                    }

                    {
                        var elm2 = new TriangleElement();
                        elm2.Nodes[0] = nodes[i + 1, j, k + 1];
                        elm2.Nodes[1] = nodes[i, j, k + 1];
                        elm2.Nodes[2] = nodes[i + 1, j, k];

                        buf.Elements.Add(elm2);
                    }

                }
            }
        }

        //elements parallel to YZ
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n - 1; j++)
            {
                for (int k = 0; k < l - 1; k++)
                {
                    {
                        var elm2 = new TriangleElement();
                        elm2.Nodes[0] = nodes[i, j, k];
                        elm2.Nodes[1] = nodes[i, j, k + 1];
                        elm2.Nodes[2] = nodes[i, j + 1, k];

                        buf.Elements.Add(elm2);
                    }

                    {
                        var elm2 = new TriangleElement();
                        elm2.Nodes[0] = nodes[i, j + 1, k + 1];
                        elm2.Nodes[1] = nodes[i, j, k + 1];
                        elm2.Nodes[2] = nodes[i, j + 1, k];

                        buf.Elements.Add(elm2);
                    }

                }
            }
        }

        //elements parallel to XY
        for (int k = 0; k < l; k++)
        {
            for (int j = 0; j < n - 1; j++)
            {
                for (int i = 0; i < m - 1; i++)
                {
                    {
                        var elm2 = new TriangleElement();
                        elm2.Nodes[0] = nodes[i, j, k];
                        elm2.Nodes[1] = nodes[i + 1, j, k];
                        elm2.Nodes[2] = nodes[i, j + 1, k];

                        buf.Elements.Add(elm2);
                    }

                    {
                        var elm2 = new TriangleElement();
                        elm2.Nodes[0] = nodes[i + 1, j + 1, k];
                        elm2.Nodes[1] = nodes[i + 1, j, k];
                        elm2.Nodes[2] = nodes[i, j + 1, k];

                        buf.Elements.Add(elm2);
                    }

                }
            }
        }

        foreach (var elm in buf.Elements)
        {
            var triElm = elm as TriangleElement;

            if (triElm == null)
                continue;

            triElm.Behavior = TriangleElementBehaviours.Shell;

            var h = 0.03;
            var w = 0.003;

            var e = 210e9;
            var nu = 0.3;

            var sec = (Sections.UniformParametric2DSection)(triElm.Section = new Sections.UniformParametric2DSection());

            var mat = triElm.Material = Materials.UniformIsotropicMaterial.CreateFromYoungPoisson(e, nu);

            sec.T = 0.003;

        }

        for (int i = 0; i < m; i++)
        {

                nodes[i, 0, 0].Constraints = Constraint.Fixed;

        }
        for (int i = 0; i < m; i++)
        {
            nodes[i, n-1, 0].Loads.Add(new NodalLoad(new Force(1500,0,0,0,0,0), LoadCase.DefaultLoadCase));
        }
        return buf;
    }

[epsi1on]

Thanks for details, Let me check and investigate the code

[epsi1on]

there is a minor problem in DKT code, which is fixed, (forgot to put Math.Abs() around jacobian determinant!) Code is updated, and your code should not give that error anymore. Also thanks for reporting the bug! cheers