Is there a limit in the number of degrees of freedom with the lm_feasible algorithm? If so, what is the limit?

Question

I am developing a finite element software that minimizes the energy of a mechanical structure. Using octave and its optim package, I run into a strange issue: The lm_feasible algorithm doesn't calculate at all when I use more than 300 degrees of freedom (DoF). Another algorithm (sqp) performs the calculation but doesn't work well when I complexify the structure and are out of my test case.

Is there a limit in the number of DoF with lm_feasible algorithm?

If so, how many DoF are maximally possible?

To give an overview and general idea of how the code works:

[x,y] = geometryGenerator()

U = zeros(lenght(x)*2,1);
U(1:2:end-1) = x;
U(2:2:end) = y;

%Non geometric argument are not optimised, and fixed during calculation
fct =@(U)complexFunctionOfEnergyIWrap(U(1:2:end-1),U(2:2:end), variousMaterialPropertiesAndOtherArgs)

para = optimset("f_equc_idx",contEq,"lb",lb,"ub",ub,"objf_grad",dEne,"objf_hessian",d2Ene,"MaxIter",1000);
[U,eneFinale,cvg,outp] = nonlin_min(fct,U,para)

Full example:

clear

pkg load optim

function [x,y] = geometryGenerator(r,elts = 100)
  teta  = linspace(0,pi,elts = 100);
  x = r * cos(teta);
  y = r * sin(teta);
endfunction

function ene  = complexFunctionOfEnergyIWrap (x,y,E,P, X,Y)
  ene = 0;
  for i = 1:length(x)-1
    ene += E*(x(i)/X(i))^4+ E*(y(i)/Y(i))^4- P *(x(i)^2+(x(i+1)^2)-x(i)*x(i+1))*abs(y(i)-y(i+1));
  endfor
endfunction

[x,y] = geometryGenerator(5,100)

%Little distance from axis to avoid division by zero
x +=1e-6;
y +=1e-6;
%Saving initial geometry
X = x;
Y = y;

%Vectorisation of the function
%% Initial vector
U = zeros(length(x)*2,1);
U(1:2:end-1) = linspace(min(x),max(x),length(x));
U(2:2:end) = linspace(min(y),max(y),length(y));

%%Constraints
Aeq = zeros(3,length(U));
%%% Blocked bottom
    Aeq(1,1) = 1;
    Aeq(2,2) = 1;
%%% Sliding top    
    Aeq(3,end-1) = 1;
%%%Initial condition
    beq = zeros(3,1);
    beq(1) = U(1);
    beq(2) = U(2);
    beq(3) = U(end-1);

    contEq = @(U) Aeq * U - beq;

%Parameter
Mat = 0.2e9;
pressure = 50;

%% Vectorized function. Non geometric argument are not optimised, and fixed during calculation
fct =@(U)complexFunctionOfEnergyIWrap(U(1:2:end-1),U(2:2:end), Mat, pressure, X, Y)

para = optimset("Algorithm","lm_feasible","f_equc_idx",contEq,"MaxIter",1000);
[U,eneFinale,cvg,outp] = nonlin_min(fct,U,para)

xFinal = U(1:2:end-1);
yFinal = U(2:2:end);

plot(x,y,';Initial geo;',xFinal,yFinal,'--x;Final geo;')

Answer 1

Finite Element Method is typically formulated as the optimal criteria for the minimization problem, which is equivalent to the Virtual Work Principle (see books like Hughes of Bathe). The Virtual Work, represents a set of linear (or nonlinear) equations which can be solved more efficiently (with fsolve).

If for some motive you must solve the problem as an optimization problem, then, if you are considering linear elasticity, your strain energy is quadratic, thus you could use the qp Octave function.

To use sparse matrices could also be helpful.