# Make uniform mesh with quad elements

I have the following code that generates a finite element mesh with triangular elements which by using command ToQuadMesh I convert to the quadrilateral elements, but I am not able to make the uniform distribution.

<< NDSolveFEM
Needs["FEMAddOns"]
r0 = .5; Emod = 2*10^6; \[Nu] = 0.3; \[Rho] = 7860; g = 10; F = 100000; k = \[Rho]*g*.5^2*Pi/2/F;
reg = ImplicitRegion[0 <= r <= r0*Exp[k*z] && 0 <= z <= 4, {r, z}];
mesh = ToElementMesh[reg, MaxCellMeasure -> {"Length" -> 1.5}, "MaxBoundaryCellMeasure" -> .59, Method -> "Continuation"]
mesh["Wireframe"]
mesh1["Wireframe"]


How to generate the uniform quad element mesh like the one that I obtain in Ansys, so that i can compare the solution from Ndsolve to Ansys?

• Comments are not for extended discussion; this conversation has been moved to chat.
– Kuba
Sep 13, 2019 at 12:32
• I reverted your recent edit, if you have a follow up/a separate question please post a new one and link this one there.
– Kuba
Sep 13, 2019 at 12:34

In the following, the parameters m and n enable you to control the mesh resolution in each direction.

Needs["NDSolveFEM"];
r0 = .5; Emod = 2*10^6; ν = 0.3; ρ = 7860; g = 10; F = 100000; k = ρ* g*.5^2*Pi/2/F;
m = 5;
n = 18;
R = ToElementMesh[Rectangle[{0, 0}, {m, n}], MaxCellMeasure -> {"Length" -> 1}, "MeshOrder" -> 1];

f = z \[Function] r0 Exp[k z];
F = {r, z} \[Function] {(r/m) f[4 z/n], (4 z/n)};
R2 = ToElementMesh["Coordinates" -> F @@@ R["Coordinates"], "MeshElements" -> R["MeshElements"]];
Show[
RegionPlot[ImplicitRegion[0 <= r <= f[z] && 0 <= z <= 4, {r, z}]],
R2["Wireframe"]
]


# Edit

This should work now for any function f.

• Thank you again but it seems to me that it doesn't approximate on the etire domain, only from 0 to 0.2, this is probably because of the Mathematica version, i have 11.1. Aug 12, 2019 at 14:14

An alternative method would be to deform a rectangle:

mesh2 = ToElementMesh[ (* determines horizontal scale: r0*Exp[k*4]/2 *)
Rectangle @@ Transpose[{{0, r0*Exp[k*4]/2}, {0, 4}}]];

drfn = Function[{r, z}, r ( 2 E^(k (-4 + z)) - 1)]; (* r * (scale - 1) *)
dzfn = Function[{r, z}, 0];     (* 0 deformation of z coordinates *)

dr = ElementMeshInterpolation[{mesh2}, rfn @@@ mesh2["Coordinates"]];
dz = ElementMeshInterpolation[{mesh2}, zfn @@@ mesh2["Coordinates"]];

mesh3 = ElementMeshDeformation[mesh2, {dr, dz}];

Show[mesh3["Wireframe"], Frame -> True]


Show[
%,
mesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]],
Frame -> True]


Update: Since the OP is using the FEMAddOns package, I would also suggest the following:

nx = 5; ny = 18;
mesh3 = StructuredMesh[Reverse@Array[
Compile @@ {{r, z}, {r*r0*Exp[k*z], z}},
{2 nx + 1, 2 ny + 1}, {{0., 1.}, {0., 4.}}], {ny, nx},
InterpolationOrder -> 2];
`