Hi Robert,
On Wed, Nov 4, 2009 at 10:49 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote:
Hi Logan
Logan Sorenson wrote:
Hi Robert,
Hope you're making good progress toward your deadlines! As such, I understand if you can't reply immediately to the following, so please take your time. :)
Two are over (well I guess), the third approaching :). So I will respond only partially for the moment, but it may point you to a right direction.
I made an initial try at implementing the electrostatic-structural coupling. You can find the code at [7]. I added a term in sfepy/terms/termsSurface.py by copying the piezoelectric terms over and modifying everything to surfaces instead of volumes. Is this a good place for the term, or should I create a source file? Also, I created input/electrostatic.py and database/electrostatic.mesh. My inspiration for this test case was the micro-bridge found in [2].
Quick first-glance notes:
- The mesh is essentially 2D - is there any reason you use a 3D mesh? Oh, I see, you have different layers and an interface in between. OK! Discovered using:
$ ./simple.py --solve-not --save-regions input/electrostatic.py
Yes, only because eventually I want to be able to handle 3D structures. :)
unfortunately this saves the regions to the medit format only...
So you have the elastic displacements in one layer, and the potential in the other?
Yes, that's right, the layer with elastic displacements should be some kind of electrical conductor (e.g., metal/semiconductor) and the region with potential should be air or vacuum where the electric field exists.
- the integral_2 is a surface integral - its kind should be 's3' (3-vertex faces).
Oh, ok, I didn't really know what to put for the kind parameter of the integral. Since I saw the other integrals marked with 'v', I thought it was related to the virtual variable also called 'v'. Now I'm guessing that 'v' actually stands for 'volume'. Are the options documented somewhere?
I don't really expect the coupling term to work yet since I haven't implemented the Maxwell stress tensor This exercise was more for me to learn how to hook terms into SfePy. Right now, I'm getting the following exception:
logan@phoenix:~/projects/sfepy$ ./simple.py input/electrostatic.py sfepy: left over: ['omega_squared', 'x_l', 'x_h', 'frequency', 'io', 'nm', 'conf_dir', 'omega', 'get_silicon_pars', 'y_l', 'geom', 'y_h', 'z_h', 'z_l', '__builtins__', 'epsilon', '__file__', 'MeshIO', 'bbox', '__name__', 'dim', '__doc__', 'get_air_pars', '_filename', 'output', 'os'] sfepy: reading mesh (database/electrostatic.mesh)... sfepy: ...done in 0.64 s sfepy: setting up domain edges... sfepy: ...done in 0.81 s sfepy: setting up domain faces... sfepy: ...done in 0.63 s sfepy: creating regions... sfepy: Y1 sfepy: Beam sfepy: Y2 sfepy: Anchor sfepy: Ground sfepy: ...done in 0.93 s sfepy: equation "1": sfepy: - 3947.841760 * dw_mass_vector.i1.Y1( silicon.density, v, u ) + dw_lin_elastic_iso.i1.Y1( silicon.lam, silicon.mu, v, u ) - dw_electrostatic_structural_coupling.i2.Beam( v, phi ) = 0 sfepy: equation "2": sfepy: dw_diffusion.i1.Y2( air.dielectric, psi, phi ) + dw_electrostatic_structural_coupling.i2.Beam( u, psi ) = 0 Traceback (most recent call last): File "./simple.py", line 100, in <module> main() File "./simple.py", line 93, in main app = SimpleApp( conf, options, output_prefix ) File "/home/logan/projects/sfepy/sfepy/applications/simple_app.py", line 50, in __init__ **kwargs ) File "/home/logan/projects/sfepy/sfepy/fem/problemDef.py", line 90, in from_conf obj.set_equations( conf.equations ) File "/home/logan/projects/sfepy/sfepy/fem/problemDef.py", line 170, in set_equations self.materials, user ) File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 202, in setup_terms eq.setup_terms( regions, variables, materials, self.caches, user ) File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 372, in setup_terms self.setup_term_args( variables, materials, user ) File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 322, in setup_term_args setup_term_args( self.terms, variables, materials, user ) File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 101, in setup_term_args raise ValueError(msg) ValueError: dw_electrostatic_structural_coupling: incompatible regions: (term, field) (0, 1 in 0
I think this makes sense since when I step through with winpdb I can see that the igs for the 'v' virtual variable get picked up from elements on both sides of the interface, but the field only has one type of element and therefore one ig. This is making me wonder about my assumption for the weakform. If you look in sfepy/terms/termsSurface.py (or generate the sfepy_manual.pdf), I'm would like to treat the coupling as in the linear traction (dw_surface_ltr) case where the stress tensor is giving rise to the surface traction is calculated from the solution for the electric field in region 2. I'm not too clear on the meaning of 'field' in SfePy, but is it possible for the field to have more than one ig? Does it make sense the way I have defined the weakform for the coupling? I'm still trying to understand the papers referenced below, but I wanted to try this way out since it makes sense in my head...
Now it depends what you need: like a jump condition(?)
- both u, phi are continuous over the whole domain, there is only the coupling term on the interface, or
- there is a discontinuity on the interface, so the coupling term is something
Both are possible :)
I guess you aim for (2.), so have a look at input/subdomains.py and the dw_jump term, that can take a boundary trace argument of a field defined "from the other side".
Ahh, yes, (2.) was more along the lines of what I was thinking. Of course (1.) may be more accurate in reality, but for most cases I think it is best to treat the conductor as ideal so that the electric field inside it goes to 0.
Also, would I have to implement some c code in sfepy/terms/extmods?
Not necessarily - you can prototype easily just using Python, only if it's too slow you go C. There are terms that have no special C functions, e.g. dw_surface_integrate_variable...
This is a good example, thanks!
Thanks, and best regards! Logan
I hope that helps a bit, continue asking :)
r.
It's a big help! I need to do more poking and prodding in the code in the places you suggested, so I'll get back with more questions later on. ;) Hope your third project can be completed soon!
Best, Logan
[snip]
[1] http://en.wikipedia.org/wiki/Maxwell_stress_tensor [2] doi:10.1016/j.na.2005.01.055 [3] doi:10.1117/12.482820 [4] http://www.ltas-vis.ulg.ac.be/cmsms/uploads/File/PHD_Rochus2006.pdf [5] http://en.wikipedia.org/wiki/Weak_formulation [6] doi:10.1016/j.jmps.2005.05.006
[7] git://github.com/logansorenson/electrostatic_structural_test.git