Operator Error in problem formulation

General Discussion
1

Hello,
I am solving a topology optimization problem using a density method.
I need to compute displacements of a 2D linear elastic structure, as well as the solution to the adjoint equation.
Below is the code I have been trying to run :

// Boundary conditions
int free = 0, loaded = 1, fixed = 2;
int[int] Labels = [0, 1, 0, 2];
mesh Th = square(200, 100, [2*x, 1*y], label=Labels);
plot(Th);

// Convergence Parameters
real p = 3;             // [/]      Stiffness Penalization power
real q = 0.5;           // [/]      Stress Penalization Power
// Material Parameters
real E1 = 200e9;        // [Pa]     Young's Modulus (isotropic scenario)
real E0 = 1e-9;         // [Pa]     Young's Modulus (void)
real nu = 0.3;          // [/]      Poisson's modulus
// Boudnary Conditions Parameters
real fx = 0;            // [N]      Horizontal force
real fy = -2e6;         // [N]      Vertical forc
// Secondary Parameters
real lambda = nu/((1+nu)*(1-2*nu));     // []       Lamé coefficient
real mu = 1/(2*(1+nu));                 // []       Lamé coefficient

// Computation of equivalent von Mises Stresses from given stress tensor components
macro EquivalentStresses(tx, ty, txy) (sqrt(((tx-ty)^2 + tx^2 + ty^2)/2.0 + 3*txy^2)) // EOM
// Computation of strain form displacement
macro e (x1,x2) [dx(x1), dy(x2), (dx(x2) + dy(x1))] // EOM
// Deviatoric stresses
macro s (x1,x2,x3) [(2*x1 - x2)/3, (2*x2 - x1)/3, x3] // EOM
// Element stiffness matrix
macro D(E) ([[E*(2*mu + lambda), E*lambda, 0], [E*lambda, E*(2*mu + lambda), 0], [0, 0, E*mu]]) // EOM

// Linear element-wise fields
fespace Vh(Th, [P1, P1]);
Vh [ux, uy], [vx, vy], [px, py];

// Constant element-wise fields
fespace Vhsig(Th, [P0, P0, P0]);
Vhsig [sigx, sigy, sigxy], [epsx, epsy, epsxy];

fespace Vh0 (Th, P0);
Vh0 E, theta, sigeq;

// Linear finite element analysis
problem Elasticity([ux,uy],[vx,vy])
    = int2d(Th) (E*((e(ux, uy)'*D(E1)*e(vx, vy))))
    - int1d(Th, loaded) (fx*vx + fy*vy)
    + on(fixed, ux=0, uy=0);

theta = 1;
E = theta^p*(1-E0) + E0;
Elasticity;
[epsx, epsy, epsxy] = e(ux, uy);
[sigx, sigy, sigxy] = sqrt(theta)*D(E1)*[epsx, epsy, epsxy];
sigeq = EquivalentStresses(sigx, sigy, sigxy);
plot(sigeq, fill=1, value=1, wait=1);

// Adjoint problem solving
problem Adjoint([px,py],[vx,vy])
    = int2d(Th) ((theta^q)*(e(vx, vy)'*D(E1)*(e(px, py) - sigeq*s(sigx, sigy, sigxy))))
    + on(fixed, px=0, py=0);

I get an operator error in the definition of the adjoint problem :

error operator -  <10LinearCombI7MGauche4C_F0E>, <d>

Hence, I have the following questions :

  • Is there a function to check the type of a variable ? When debugging, I couldn’t find which element was of type , is there a way to know ?
  • I was able to compute e(px, py) - sigeq*s(sigx, sigy, sigxy)) outside of the problem formulation, why is it causing an issue here ?
  • Any advice on how to fix it ?

Thank you in advance for your time and answers,
Best regards,
Hugo

2

The problem seems to be here the standard one:
e(px, py) depends on the unknown [px,py], whereas sigeq*s(sigx, sigy, sigxy) does not. Thus you have to put them in two different int2d.
This is a requirement of FF+ that we easily forget!

3

Thank you very much
Following my initial post, is there a way in freefem to check data type ?

4

Up to my knowledge there is no, as far as I understand your question.

1 Like