Deformed mesh plot displaced

General Discussion
1

I’ve computed a solution for a solid wave in a rectangular bar. This particular case has an analytical solution, and line plots confirm that it’s been found. However, plots of the deformed mesh have unexpected creases, as shown (for real and imaginary parts of the solution) here:

defmshr
defmshr1404×804 13.8 KB

defmshi
defmshi1424×602 12.5 KB

The code that produces these (and the line plots) is listed below. Quiver and contour plots are also generated, in case they’re helpful.

// Waves on a bar of rectangular cross-section.
// Developed from gempro-code-bareigFF.cpp

load “msh3”

// 1. Geometric and Material Parameters
real Lx = 2.0 * sqrt(2); real nu = 0.3; // Used for Lame modes
real Ly = 2.0;
real mu = 1.0; // Shear modulus
real lam = 2nu/(1 - 2nu); // Lame’s first parameter (lambda)
real rho = 1.0; // Density
real k = 3.1416; // Given real wavenumber (M&F zeta; T&I gamma)

int nx = 40;
int ny = 30;

// 2. Mesh Generation
mesh Th = square(nx, ny, [xLx, yLy]);
plot(Th, wait=true);

// 3. Finite Element Space
// P2 elements for all three displacement components, complex-valued
fespace Vh(Th, [P2, P2, P2]);
Vh [u1, u2, u3];
Vh [v1, v2, v3];

// 5. Variational Formulations
// Stiffness Matrix (Strain Energy)
varf vK([u1, u2, u3], [v1, v2, v3]) =
int2d(Th)(
lam * (dx(u1)+dy(u2)-1iku3) * (dx(v1)+dy(v2)+1ikv3)

  • 2.0 * mu * (
    dx(u1) * dx(v1)
  • dy(u2) * dy(v2)
  • k^2 * u3 * v3)
  • mu * (
    (dy(u1)+ dx(u2)) * (dy(v1)+ dx(v2))
  • (-1iku2 + dy(u3)) * (1ikv2 + dy(v3))
  • (dx(u3) - 1iku1) * (dx(v3) + 1ikv1)
    )
    );

// Mass Matrix (Kinetic Energy)
varf vM([u1, u2, u3], [v1, v2, v3]) =
int2d(Th)(
rho * (u1v1 + u2v2 + u3*v3)
);

// Assemble Matrices
matrix K = vK(Vh, Vh);
matrix M = vM(Vh, Vh);

// 6. Eigenvalue Solution via ARPACK
int nev = 10; // Number of requested eigenvalues
complex[int] ev(nev); // Array to store eigenvalues
Vh[int] eVec1, eVec2, eVec3; // Eigenvectors

// Solve generalized eigenvalue problem A * x = lambda * B * x
// where lambda = omega^2
real sig = 20.0;
matrix Kshft = K - sig*M;
set(Kshft, solver=sparsesolver);

int kfound = EigenValue(Kshft, M, sigma=sig, value=ev, vector=eVec1, nev=nev);

// 8. Plot eigenvectors
fespace Vhplt(Th,P2);
int kplt = 2; // mode to plot, counting from one
Vhplt uxr = real(eVec1[kplt-1]);
Vhplt uxi = imag(eVec1[kplt-1]);
Vhplt uyr = real(eVec2[kplt-1]);
Vhplt uyi = imag(eVec2[kplt-1]);
Vhplt uzr = real(eVec3[kplt-1]);
Vhplt uzi = imag(eVec3[kplt-1]);
meshS ThNew = movemesh23(Th, transfo=[x + uxr, y + uyr, uzr]);
plot(ThNew, wait=true, value=true, ps=“defmshr.eps”, cmm=“Deformed mesh, real displacements”);
meshS ThNewi = movemesh23(Th, transfo=[x + uxi, y + uyi, uzi]);
plot(ThNewi, wait=true, value=true, ps=“defmshi.eps”, cmm=“Deformed mesh, imag displacements”);

// — Define Slice Parameters —
real xSlc = 0.5*Lx;
real yMin = 0.0, yMax = Ly;
int Ny = 50; // Grid resolution for the slice
real[int] ySlc(Ny+1), uyrsl(Ny+1), uyisl(Ny+1), uzrsl(Ny+1), uzisl(Ny+1);

real deey = (yMax - yMin) / Ny;

for (int j = 0; j <= Ny; j++) {

ySlc[j] = yMin + j * deey;

// FreeFEM automatically interpolates the FE functions
uyrsl[j] = uyr(xSlc, ySlc[j]);
uyisl[j] = uyi(xSlc, ySlc[j]);
uzrsl[j] = uzr(xSlc, ySlc[j]);
uzisl[j] = uzi(xSlc, ySlc[j]);

}

plot([uxr,uyr],value=true,cmm=“Re(ux,uy)”,wait=true);
plot([ySlc,uyrsl],value=true,cmm=“Re(uy) vs y”,wait=true);

plot(uzr,cmm=“Re(uz)”,value=true, wait=true);
plot([ySlc,uzrsl],value=true,cmm=“Re(uz) vs y”,wait=true);

plot([uxi,uyi],value=true,cmm=“Im(ux,uy)”,wait=true);
plot([ySlc,uyisl],value=true,cmm=“Im(uy) vs y”,wait=true);

plot(uzi,value=true,cmm=“Im(uz)”,wait=true);
plot([ySlc,uzisl],value=true,cmm=“Im(uz) vs y”,wait=true);