Logan Sorenson wrote:
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
coupling. You can find the code at . 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 .
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
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?
- the integral_2 is a surface integral - its kind should be 's3' (3-vertex
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
logan@phoenix:~/projects/sfepy$ ./simple.py input/electrostatic.py
sfepy: left over: ['omega_squared', 'x_l', 'x_h',
'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',
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: ...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 )
sfepy: equation "2":
sfepy: dw_diffusion.i1.Y2( air.dielectric, psi, phi )
+ dw_electrostatic_structural_coupling.i2.Beam( u, psi )
Traceback (most recent call last):
File "./simple.py", line 100, in <module>
File "./simple.py", line 93, in main
app = SimpleApp( conf, options, output_prefix )
line 50, in __init__
File "/home/logan/projects/sfepy/sfepy/fem/problemDef.py", line 90,
obj.set_equations( conf.equations )
File "/home/logan/projects/sfepy/sfepy/fem/problemDef.py", line 170,
self.materials, user )
File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 202,
eq.setup_terms( regions, variables, materials, self.caches, user )
File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 372,
self.setup_term_args( variables, materials, user )
File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 322,
setup_term_args( self.terms, variables, materials, user )
File "/home/logan/projects/sfepy/sfepy/fem/equations.py", line 101,
ValueError: dw_electrostatic_structural_coupling: incompatible
regions: (term, field) ([0, 1](v) in (displacement)
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:
1. both u, phi are continuous over the whole domain, there is only the coupling
term on the interface, or
2. there is a discontinuity on the interface, so the coupling term is something
like a jump condition(?)
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
Also, would I have to implement some c code in
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.
Thanks, and best regards!
I hope that helps a bit, continue asking :)