8
$\begingroup$

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.

<< NDSolve`FEM`
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"] 
 mesh1 = ToQuadMesh[mesh];
mesh["Wireframe"]
mesh1["Wireframe"]

enter image description here

enter image description here

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?

enter image description here

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

2 Answers 2

10
$\begingroup$

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

Needs["NDSolve`FEM`"];
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"]
 ]

enter image description here

Edit

This should work now for any function f.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – LejlaS
    Aug 12, 2019 at 14:14
8
$\begingroup$

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];

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.