Im new to FreeFem, can someone clear up what my program is actually doing?

General Discussion
1

Basically as the title says Im learning FreeFem for a university project and I want to better understand what the program is actually doing. I have the following program modified from https://modules.freefem.org/modules/elasticity/:

load "msh3"
load "gmsh"

//* mesh
int Clamped = 27;		//Beam fixed 
mesh3 Th = readmesh3("rails");
plot(Th, wait=true);

//* Some parameters defined here...
real E = 210e9;         //* Young's Modulus
real Nu = 0.27;         //* poisson's ratio
real Rho = 7800;        //* density
real Gravity = -9.81;   //* gravitational acceleration

//* fespace
func Pk = P1;
fespace Uh(Th, [Pk, Pk, Pk]);
Uh [ux, uy, uz];
Uh [vx, vy, vz];
Uh [uxp, uyp, uzp];
Uh [uxpp, uypp, uzpp];

//* Macro
real sqrt2 = sqrt(2.);
macro Epsilon(ux, uy, uz) [dx(ux), dy(uy), dz(uz),
	(dz(uy)+dy(uz))/sqrt2,
	(dz(ux)+dx(uz))/sqrt2,
	(dy(ux)+dx(uy))/sqrt2] //
macro Divergence(ux, uy, uz) (dx(ux) + dy(uy) + dz(uz)) //


//* Problem
real Mu = E/(2.*(1.+Nu));
real Lambda = E*Nu/((1.+Nu)*(1.-2.*Nu));


varf vElasticity ([ux, uy, uz], [vx, vy, vz])
	= int3d(Th)(
		  Lambda * Divergence(vx, vy, vz) * Divergence(ux, uy, uz)
		+ 2. * Mu * (
			  Epsilon(vx, vy, vz)' * Epsilon(ux, uy, uz)
		)
	)
	+ int3d(Th)(
		  Rho * Gravity * vz
	)
    - int2d(Th, 30)(
        30000*vy
    )
	+ on(Clamped, ux=0, uy=0, uz=0)
        + on(26, ux=0, uy=0, uz=0)
	;

matrix Elasticity = vElasticity(Uh, Uh, solver=sparsesolver);
real[int] ElasticityBoundary = vElasticity(0, Uh);
ux[] = Elasticity^-1 * ElasticityBoundary;

//*Movemesh
real dmax = ux[].max;
cout << "max displacement = " << dmax << endl;
real coef = 50;
//Th = movemesh(Th, [x+ux, y+uy, z+uz]);
//[ux, uy, uz] = [ux, uy, uz];

mesh3 Thm = movemesh3(Th, transfo=[x+ux*coef, y+uy*coef, z+uz*coef]);
Thm = change(Thm);

//*Plot
//plot([ux, uy, uz], value=true);
plot(Th, Thm, value=true, wait=true, cmm="coef amplification = "+coef);

A few of the points Im confused about are:

  1. under fespace, what does Uh [uxp, uyp, uzp] and Uh [uxpp, uypp, uzpp] do/represent? These aren’t used elsewhere so does not sure what it does.

  2. when solving the matrix vector equation Ax=b, why do the results get written to ux[]? I have applied the forces in x,y and z direction but it doesn’t seem to matter? so why the x in ux[]?

  3. When applying external forces I am using int2d(Th, No.). What does int2d represent? I also found int3d, int1d and int0d. Would int0d for example be a point force? And the extra number I have to pass, does that correspond to a line, point or surface on my mesh? (Im preparing in gmsh fyi).

  4. The program I have can display the displacements and not the stresses. I am actually mostly interested in these stresses. What modifications would I need to calculate the stresses in the structure?

Sorry if my code is a bit bad, Im very new to this software and am trying to understand it better. Any help is greatly appreciated!

p.s. This is my result so far:

original vs deflected structure
original vs deflected structure1093Ă—623 64.5 KB

2
  1. under fespace, what does Uh [uxp, uyp, uzp] and Uh [uxpp, uypp, uzpp] do/represent? These aren’t used elsewhere so does not sure what it does.

I don’t think these are used… Does the code still run with those lines removed?

  1. when solving the matrix vector equation Ax=b, why do the results get written to ux[]? I have applied the forces in x,y and z direction but it doesn’t seem to matter? so why the x in ux[]?

Yes, this is confusing when I started using FreeFEM as well. Even though the FEspace consists of ux,uy,uz, the name of the vector that stores all of this information is just the first element of the array, ie ux[].

  1. When applying external forces I am using int2d(Th, No.). What does int2d represent? I also found int3d, int1d and int0d. Would int0d for example be a point force? And the extra number I have to pass, does that correspond to a line, point or surface on my mesh? (Im preparing in gmsh fyi).

int2D(Th, No.) means 2d integration over the part of mesh Th numbered No. In 2d problems where you integrate over the whole mesh you dont need to provide the No argument.

  1. The program I have can display the displacements and not the stresses. I am actually mostly interested in these stresses. What modifications would I need to calculate the stresses in the structure?

I’m guessing that this is what the [uxp,uyp,uzp] and [uxpp,uypp,uzpp] arrays were originally used for. (I’m guessing the p stands for prime?) Just define these arrays based on the stress relationship and your known displacements.

Hope that helps!

2 Likes
3

Hi Chris,

thanks a lot for your answer. As you said commenting out the uxp and uxpp doesn’t affect the program so I guess its not important in that sense.

I did manage to get the von mises stresses from my code by just using the equations from Lamé.

thank you for your help!

1 Like