On 06/04/2012 01:46 PM, Steve Spicklemire wrote:

Hi Robert,

On Jun 4, 2012, at 1:44 AM, Robert Cimrman wrote:

...

Hi Steve,

On 06/04/2012 06:01 AM, steve wrote:

...

Hi Folks,

I'm extremely new to finite elements in general and sfepy in particular, but I have a couple of problems I'm interested in, and I'm hoping that sfepy may provide an easy way to get started. Both of these problems (1 acoustic normal modes in cavities with cylindrical symmetry and 2 electrostatic fields in cavities bounded by various equipotential surfaces) have in common a axisymmetric or cylindrical symmetry. I found this discussion

<http://groups.google.com/group/sfepy-devel/msg/17364ae1fcc57259>

but I didn't see any real follow up. Has there been more done on this?

I am slowly working on the plate/shell aspect, but not on the symmetry.

...

also.. I've spent the last couple of days skimming Hughes finite element analysis book to get my head screwed on straight WRT FEM and I found this interesting tidbit:

<http://books.google.com/books? id=cHH2n_qBK0IC&printsec=frontcover&dq=finite%20element%20method %20hughes&hl=en&sa=X&ei=ywHMT- bKBcfK2AWowNzZCw&output=reader&pg=GBS.PA101>

The link jumps to the front page - could you tell us the page?

Sorry.. it's page 101:

<http://books.google.com/books?id=cHH2n_qBK0IC&printsec=frontcover&dq=finite%20element%20method%20hughes&hl=en&sa=X&ei=ywHMT-bKBcfK2AWowNzZCw&output=reader&pg=GBS.PA101>

I used the google web page mailer to send... bad wrapping apparently!

Or my browser sucks. Anyway, I have found [1], where one can browse all the pages...

[1] http://seeinside.doverpublications.com/dover/0486411818

...

...

Since heat conduction is extremely similar to electrostatics (and only a term or so different from acoustics) I hoped that it might be possible to get away with a 2D mesh for my problem(s) rather than running a full 3D mesh of the highly symmetric cases. So.. I guess my question is, what's the easiest way to implement this feature in sfepy? Can anyone suggest a path to follow?

There are at least two paths:

a) Because the terms in symmetrical formulations are often just like the terms in the cartesian grid but multiplied by a coordinate-dependent factor (e.g. 1/r), the easiest way IMHO is to incorporate those factors to a material parameter of each term. The material parameters can be given as general functions of coordinates, and you get the coordinates, so all the information is there. Not all terms, however, would be possible to reuse in this way. The Laplacian operator, for example, would have to be reimplemented as a new term for the cylindrical symmetry.

This sounds like a reasonable approach.... so I guess that means I'd need to create a subclass of terms.termsLaplace.LaplaceTerm, maybe like "CylindricalLaplaceTerm" or something?

Or DiffusionTerm (LaplaceTerm is just a very simple subclass of that). The name could be DiffusionCylindricalTerm (dw_diffusion_cyl ?) so that it stays next to the original term in the term table. Anyway, it would be great if you could give it a shot. Here are some comments:

• Do you already have the equations written in the weak form? If not, start with that, I would like to see them. Maybe more than one term would be needed/better.

• I would pass in the values of r as a material parameter. In that way a user just writes a function for computing it given the physical quadrature point coordinates. It seems to me a little bit easier, than computing r in the term (using "qps = self.get_physical_qps()", self is a term instance).

• Writing a new term is not that simple, as it often requires going into C+Cython - the Python class then just dispatches the data. It is also probably under-documented (and changes), so ask in case of troubles.

I suppose there would have to be some kind of convention about the meaning of different dimensions, like x -> r, y -> z or vice-versa?

I would suggest the mapping used in Hughes: x -> r, y -> z, z -> phi.

r.

Would be nice to solve this in a general way.. but first things first. ;-)

...

b) Support somehow the symmetry internally, that is, do the above under the hood. This would be more difficult, and I am not sure if it is worth it. It depends on the number of new terms/reused terms that would have to implemented/or not in the case (a). Essentially, (b) == (a) with automatic use of the right terms given the symmetry.

Yes.. I agree.. that sounds like a lot more work. ;-)

thanks for the input! -steve

...

So solving your problems should be possible, with some work. I think I raised more questions than answered, so feel free to ask!

Also, if others have some remarks/opinions, do not hesitate to write them here.

Cheers, r.

On 06/04/2012 04:48 PM, steve wrote:

Thanks Robert,

Or DiffusionTerm (LaplaceTerm is just a very simple subclass of that). The

...

name could be DiffusionCylindricalTerm (dw_diffusion_cyl ?) so that it stays next to the original term in the term table. Anyway, it would be great if you could give it a shot. Here are some comments:

• Do you already have the equations written in the weak form? If not, start with that, I would like to see them. Maybe more than one term would be needed/better.

Well... the equations are embarrassingly simple, but as I'm completely new to the FEM I've zero experience converting to weak form, but I'll give it a go. For the electrostatics part it's just Laplace's equation in the interior with Dirichlet BCs on the domain boundary. For the electrostatics then it's just:

\nabla^2 u = 0

subject to u(r)=v(r) on the boundary. u is the scaler potential and v(r) is just it's value on the boundary of the cavity.

Which I guess in weak form would look like:

\int_\Omega \vec{\nabla} u \cdot \vec{\nabla} s \, dV + \int_\Gamma s v(\Gamma) d\Gamma = 0 ;\, \forall s

The acoustics problem is an eigenvalue problem (looking for normal modes) in a cylindrically symmetric enclosure, so it's more like:

\nabla^2 u + k^2 u = 0

subject to \hat{n}\cdot \vec{\nabla}u = 0, on the boundary. u in this case is the velocity potential and k is the wavenumber. So.. in my naive world the weak form of this one would be something like:

\int_\Omega \vec{\nabla} u \cdot \vec{\nabla} s \, dV + k^2 \int_\Omega s u dV + BT = 0 \forall s

Where "BT" is some kind of boundary term, that I haven't sorted out yet. Since it's a Neumann BC I'm not sure how it gets incorporated into the weak form statement.

Yes, those are simple and standard. But I meant the weak form of the equations in the cylindrical (and other) coordinates :). (With all those 1/r and other stuff.)

As far as C+Cython, I've used those before with numpy in other contexts, so that's not too scary, I'll see what I can do!

Cool, I am glad to hear that!

r.

On 06/04/2012 05:33 PM, steve wrote:

Hi Robert,

...

Where "BT" is some kind of boundary term, that I haven't sorted out yet.

...

Since it's a Neumann BC I'm not sure how it gets incorporated into the weak form statement.

Yes, those are simple and standard. But I meant the weak form of the equations in the cylindrical (and other) coordinates :). (With all those 1/r and other stuff.)

Ah... good point. So assuming u has no theta dependence (which I realize as I'm typing this is not necessarily true in the acoustics case.. but more on that later) the volume integrals reduce to:

2\pi \int_\Omega (u_{,r} s_{,r} + u_{,z} s_{,z}) r\,dr\,dz

and

2\pi k^2 \int_\Omega (s u ) r\,dr\,dz (for the acoustics case only)

OK. So it seems that for the above, you could use directly dw_laplace and dw_volume_dot terms, with material arguments equal to 2\pi r and 2\pi k^2 r, respectively. Or just r, and put the known constants directly into the equations string.

Where \Omega is now the 2-D interior of the 'r-z' plane. \Gamma is reduced to a 1D boundary, but d\Gamma is going to get a factor of 2\pi r as well to take into account the 'distance' around in the theta direction.

Since the weak form only depends on the gradient, and '1/r' only shows up in the theta part of the gradient, it doesn't appear at this point. The only 'r' factor is in the volume element itself.

I thought (googled) that the 'r' part of the cylindrical Laplacian was (d = \partial):

\frac{1}{r} \frac{d}{dr} (k r \frac{du}{dr}), which is

k \frac{d^2 u}{dr^2} + \frac{k}{r} \frac{du}{dr}.

so in weak form it comes to a laplace-like term with different coefficients by the 'u' and 'z' parts, and a gradient-like term w.r.t. 'r' (the second one in the line above). Just curious...

Now.. about the theta dependence of u. In the acoustics case there will be modes that have theta dependence, but they will have a continuity condition that will require the theta part go like exp(+mj theta) where m is an integer, so this will add a term to the laplacian. But I can worry about that after I get the basic concept working I think.

Good.

Cheers, r.

On 06/04/2012 07:14 PM, steve wrote:

...

OK. So it seems that for the above, you could use directly dw_laplace and dw_volume_dot terms, with material arguments equal to 2\pi r and 2\pi k^2 r, respectively. Or just r, and put the known constants directly into the equations string.

Hmm.. OK. I'll try to wrap my head around that. ;-)

I will try to support you, don't worry.

...

...

Where \Omega is now the 2-D interior of the 'r-z' plane. \Gamma is reduced to a 1D boundary, but d\Gamma is going to get a factor of 2\pi r as well to take into account the 'distance' around in the theta direction.

Since the weak form only depends on the gradient, and '1/r' only shows up in the theta part of the gradient, it doesn't appear at this point. The only 'r' factor is in the volume element itself.

I thought (googled) that the 'r' part of the cylindrical Laplacian was (d

\partial):

\frac{1}{r} \frac{d}{dr} (k r \frac{du}{dr}), which is

k \frac{d^2 u}{dr^2} + \frac{k}{r} \frac{du}{dr}.

so in weak form it comes to a laplace-like term with different coefficients by the 'u' and 'z' parts, and a gradient-like term w.r.t. 'r' (the second one in the line above). Just curious...

So let's see if I can work this out. To convert to the weak formulation you'd start with the laplacian and integrate against an arbitrary function. If we just focus on the 'r' part it would go like so:

\int s \frac{1}{r} \frac{d}{dr} (k r \frac{du}{dr}) dV

but the r part of dV is r dr (forgetting about the 2\pi for now) so that's

\int s \frac{1}{r} \frac{d}{dr} (k r \frac{du}{dr}) r dr

which looks like:

k \int s \frac{d}{dr} (r \frac{du}{dr}) dr

But if that gets integrated by parts it's works out to be

k (s r \frac{du}{dr})\big\vert_{\Gamma} - k \int \frac{ds}{dr} \frac{du}{dr} r dr

Which is I think more or less what I had before. Have I got that right?

Probably :) (would have to compile the LaTeX) The difference is, that I thought about first making the d/dr derivative by parts, and then doing the weak form from the resulting two terms, while you did the opposite. Not sure which way is better...

r.

On 06/06/2012 04:33 PM, steve wrote:

Hi Robert,

...

...

OK. So it seems that for the above, you could use directly dw_laplace and

...

dw_volume_dot terms, with material arguments equal to 2\pi r and 2\pi k^2 r, respectively. Or just r, and put the known constants directly into the equations string.

Hmm.. OK. I'll try to wrap my head around that. ;-)

I will try to support you, don't worry.

Well... I've had some success here, I think!

Glad to hear that!

here's the setup file 'es-test.py' which is trying to model a cylindrically symmetric electrostatic lens. The results compare favorably with a relaxation code with the same setup. Now of course I have a couple of issues that I hope you can help me with.

Nice! Could you send me also the mesh? (off-list, if it's big...) So this just solves the direct problem, not the eigenvalue problem, right?

1. In the relaxation version I have to impose the du/dr=0 boundary condition at r=0 at each time step. I made no mention of that BC in the setup file, but the solution seems to respect it more or less automatically. I remember in the examples seeing mention of "zero flux" over the boundary. Is this just a consequence of leaving out the boundary integral in the weak formulation?

Yes, by leaving out the boundary integral you impose the zero flux.

2. I've developed a simple wxpython app that allows a user to specify voltages on, and position/geometry of the lenses and then compute particle trajectories through the system. My relaxation algorithm is currently limited to a uniform rectangular mesh. Part of my motivation for pursuing sfepy is the possibility of using variable/arbitrary mesh density. However, to run this test I had to generate a mesh in gmsh and then hand-edit the output to get a .mesh file that was acceptable to sfepy. I'd like to increase the density of nodes around the corners of the conductors in the problem (especially around the needle where the particles are introduced). Any thoughts on dynamically building meshes in code to accomplish this?

What had to be changed in the gmsh-generated mesh? INRIA medit format (.mesh) should work - if it does not, it should be fixed.

As sfepy cannot do the meshing/adaptive refinement, so you will need to call an external mesher (e.g. gmsh) to do that. Skimming the gmsh docs seems to suggest that it's possible to drive it from Python. If you choose to go that path, I would really like to see the result - having a MeshIO class taking a gmsh mesh description file and returning the mesh would be great.

Anyway.. thanks for the help! I think it's working pretty well after I finally groked how variable materials work in the setup file.

Hth, r.

That's ok - this can happen, when the initial residual is high. Try increasing 'lin_red' option of the Newton solver to 1e1.

More details: Usually, one wants the linear solver precision match the nonlinear solver precision, so sfepy warns you when the linear system solution rezidual norm is > than eps_a * lin_red of the nonlinear solver, and default values are eps_a = 1e-10, lin_red = 1. In this case the relative rezidual is 2.823174e-15 which is about as good as it can get in double precision.

Conclusion: ignore the message :)

r.

On 06/06/2012 05:35 PM, steve wrote:

Hmm.. also I just noticed that simple.py complains that all is not well. The output really *looks* reasonable though!

sfepy: reading mesh (/Users/steve/Development/sfepy/es-lens.mesh)... sfepy: ...done in 0.10 s sfepy: warning: bad element orientation, trying to correct... sfepy: ...corrected sfepy: creating regions... sfepy: Omega_Source sfepy: Omega_Lens2 sfepy: Omega_Lens1 sfepy: Omega_Target sfepy: Omega sfepy: Omega_Needle sfepy: ...done in 0.04 s sfepy: equation "Laplace equation": sfepy: dw_laplace.i1.Omega( m.val, v, u ) = 0 sfepy: setting up dof connectivities... sfepy: ...done in 0.00 s sfepy: using solvers: nls: newton ls: ls sfepy: updating variables... sfepy: ...done sfepy: matrix shape: (4465, 4465) sfepy: assembling matrix graph... sfepy: ...done in 0.00 s sfepy: matrix structural nonzeros: 30483 (1.53e-03% fill) sfepy: updating materials... sfepy: m sfepy: ...done in 0.00 s sfepy: nls: iter: 0, residual: 4.936809e+04 (rel: 1.000000e+00) sfepy: rezidual: 0.00 [s] sfepy: solve: 0.03 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.335327e-10< 1.000000e-10) sfepy: nls: iter: 1, residual: 1.393747e-10 (rel: 2.823174e-15)

-steve

On Wednesday, June 6, 2012 9:04:39 AM UTC-6, Robert Cimrman wrote:

...

...

Hi Robert,

...

...

OK. So it seems that for the above, you could use directly dw_laplace and

...

dw_volume_dot terms, with material arguments equal to 2\pi r and 2\pi k^2 r, respectively. Or just r, and put the known constants directly into

On 06/06/2012 04:33 PM, steve wrote: the

...

...

...

...

equations string.

Hmm.. OK. I'll try to wrap my head around that. ;-)

I will try to support you, don't worry.

Well... I've had some success here, I think!

Glad to hear that!

...

here's the setup file 'es-test.py' which is trying to model a cylindrically symmetric electrostatic lens. The results compare favorably with a relaxation code with the same setup. Now of course I have a couple of issues that I hope you can help me with.

Nice! Could you send me also the mesh? (off-list, if it's big...) So this just solves the direct problem, not the eigenvalue problem, right?

...

1. In the relaxation version I have to impose the du/dr=0 boundary condition at r=0 at each time step. I made no mention of that BC in the setup file, but the solution seems to respect it more or less automatically. I remember in the examples seeing mention of "zero flux" over the boundary. Is this just a consequence of leaving out the boundary integral in the weak formulation?

Yes, by leaving out the boundary integral you impose the zero flux.

...

2. I've developed a simple wxpython app that allows a user to specify voltages on, and position/geometry of the lenses and then compute particle trajectories through the system. My relaxation algorithm is currently limited to a uniform rectangular mesh. Part of my motivation for pursuing sfepy is the possibility of using variable/arbitrary mesh density. However, to run this test I had to generate a mesh in gmsh and then hand-edit the output to get a .mesh file that was acceptable to sfepy. I'd like to increase the density of nodes around the corners of the conductors in the problem (especially around the needle where the particles are introduced). Any thoughts on dynamically building meshes in code to accomplish this?

What had to be changed in the gmsh-generated mesh? INRIA medit format (.mesh) should work - if it does not, it should be fixed.

As sfepy cannot do the meshing/adaptive refinement, so you will need to call an external mesher (e.g. gmsh) to do that. Skimming the gmsh docs seems to suggest that it's possible to drive it from Python. If you choose to go that path, I would really like to see the result - having a MeshIO class taking a gmsh mesh description file and returning the mesh would be great.

...

Anyway.. thanks for the help! I think it's working pretty well after I finally groked how variable materials work in the setup file.

Hth, r.

On 06/06/2012 09:40 PM, steve wrote:

OK.. on to the eigenvalue problem.

I started with the quantum program... and stripped out all the quantum stuff. ;-)

Basically the acoustics problem is exactly the square well potential (V=0 everywhere) BUT with different boundary conditions.

In quantum_common.py we've got:

 ebc_1 = {
     'name' : 'ZeroSurface',
     'region' : 'Surface',
     'dofs' : {'Psi.0' : 0.0},
 }

Which forces psi to be zero on the surface of Omega. Basically Dirichlet conditions. For acoustics the velocity potential needs to have no normal gradient at the surfaces.. more or less Neumann conditions. Is there an easy way to implement that?

What is the exact definition of the velocity potential? If it's really a Neumann-like boundary integral, having it zero = omitting the term in equations + no Dirichlet boundary conditions (ebcs) at all.

r.