Hi,
I was wondering how to implement the electrostatic-elastic coupling in SfePy. Basically, I think this coupling is described by the Maxwell stress tensor [1] which gives the stress due to the electrostatic field (setting the magnetic field to 0 in the quasistatic approximation) on the surface of a conductor.
I think the coupling can be handled in a similar manner with the piezoelectric coupling. For example, in input/piezo.py, the equations are:
equations = { '1' : """- %f * dw_mass_vector.i1.Y( inclusion.density, v, u ) + dw_lin_elastic_iso.i1.Y( inclusion.lam, inclusion.mu, v, u ) - dw_piezo_coupling.i1.Y2( inclusion.coupling, v, phi ) = 0""" % omega_squared, '2' : """dw_diffusion.i1.Y( inclusion.dielectric, psi, phi ) + dw_piezo_coupling.i1.Y2( inclusion.coupling, u, psi ) = 0""", }
I think there should be a similar set of equations for the electrostatic-elastic coupling case, but the integration should be done over the interface of the elastic material and the vacuum regions. I think that can be handled by adding a 2D integral in the problem definition. What I'm missing is how to define the coupling term in SfePy. I think the form of the term is basically the same as dw_surface_ltr in sfepy_manual.pdf. Also, it looks like the gradient of the electric potential can be obtained from cachesBasic.
Is it possible to use the gradient of the electric potential information to calculate the Maxwell stress tensor using appropriate vector operations (divergence and dot product), along with the dielectric tensor? Would this be doable through the dw_surface_ltr term? Or is it better to implement a new coupling term (based off dw_surface_ltr) to handle this (and how would one go about that :) )?
Hope everything is clear! Any advice is greatly appreciated!
Thanks, Logan
Hi Logan,
I am on a very thin ice with anything "electric" in its name, so do not take me too seriously... Let me ask some questions first:
- do we need the edge (Nedelec) elements for this (in 3D)?
- is the model similar to what is described in [2]? (I can access the fulltext.)
- do you have somewhere a complete weak formulation of the problem?
SfePy does not have the edge elements now, so if you need them, it will be a nontrivial implementation effort. Otherwise it should be easy.
Logan Sorenson wrote:
Hi,
I was wondering how to implement the electrostatic-elastic coupling in SfePy. Basically, I think this coupling is described by the Maxwell stress tensor [1] which gives the stress due to the electrostatic field (setting the magnetic field to 0 in the quasistatic approximation) on the surface of a conductor.
I think the coupling can be handled in a similar manner with the piezoelectric coupling. For example, in input/piezo.py, the equations are:
equations = { '1' : """- %f * dw_mass_vector.i1.Y( inclusion.density, v, u ) + dw_lin_elastic_iso.i1.Y( inclusion.lam, inclusion.mu, v, u ) - dw_piezo_coupling.i1.Y2( inclusion.coupling, v, phi ) = 0""" % omega_squared, '2' : """dw_diffusion.i1.Y( inclusion.dielectric, psi, phi ) + dw_piezo_coupling.i1.Y2( inclusion.coupling, u, psi ) = 0""", }
I think there should be a similar set of equations for the electrostatic-elastic coupling case, but the integration should be done over the interface of the elastic material and the vacuum regions. I think that can be handled by adding a 2D integral in the problem definition. What I'm missing is how to define the coupling term in SfePy. I think the form of the term is basically the same as dw_surface_ltr in sfepy_manual.pdf. Also, it looks like the gradient of the electric potential can be obtained from cachesBasic.
Is it possible to use the gradient of the electric potential information to calculate the Maxwell stress tensor using appropriate vector operations (divergence and dot product), along with the dielectric tensor? Would this be doable through the dw_surface_ltr term? Or is it better to implement a new coupling term (based off dw_surface_ltr) to handle this (and how would one go about that :) )?
I would probably implement a new term (called e.g. dw_electric_coupling?) by following the pattern of dw_piezo_coupling (sfepy/terms/termsPiezo.py), but working with surface integrals. Notice how dw_piezo_coupling treats its threefold nature - it can be called in the "grad" (as in eq. '1'), "div" (eq/ '2') and "eval" (both v, phi known) forms. The last one is not needed to solve problems. The grad, div terminology comes from the Stokes problem.
Hope everything is clear! Any advice is greatly appreciated!
I will try to assist you with the practical implementation, of course, but not much with the theory - you are the expert :).
cheers, r.
[1] http://en.wikipedia.org/wiki/Maxwell_stress_tensor [2] doi:10.1016/j.na.2005.01.055
Hi Robert,
On Wed, Oct 28, 2009 at 4:23 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi Logan,
I am on a very thin ice with anything "electric" in its name, so do not take me too seriously... Let me ask some questions first:
Don't worry, I'm on thin ice with lots of things! :)
- do we need the edge (Nedelec) elements for this (in 3D)?
I haven't encountered the Nedelec elements before, so I'll do some more research to see if they are necessary. I'm hoping not :).
- is the model similar to what is described in [2]? (I can access the fulltext.)
This is actually a really good paper. Thanks for pointing it out to me! [3] by the same author (found through the references of [2]) has some more detailed derivations. Also, her PhD thesis [4] is quite good. I think this is a better description of the problem than I could provide.
- do you have somewhere a complete weak formulation of the problem?
I haven't worked out the weak formulation yet. I think I need to consider some of the issues in [2-4], such as handling the changes in the electrical domain due to the displacements. I think the original construction below may work in the linear small displacement regime. By the way, do you know of a good reference which explains the weak form derivations? I got the basics from the sfepy_manual.pdf, but I think I need more practice for general things.
SfePy does not have the edge elements now, so if you need them, it will be a nontrivial implementation effort. Otherwise it should be easy.
Logan Sorenson wrote:
Hi,
I was wondering how to implement the electrostatic-elastic coupling in SfePy. Basically, I think this coupling is described by the Maxwell stress tensor [1] which gives the stress due to the electrostatic field (setting the magnetic field to 0 in the quasistatic approximation) on the surface of a conductor.
I think the coupling can be handled in a similar manner with the piezoelectric coupling. For example, in input/piezo.py, the equations are:
equations = { '1' : """- %f * dw_mass_vector.i1.Y( inclusion.density, v, u ) + dw_lin_elastic_iso.i1.Y( inclusion.lam, inclusion.mu, v, u ) - dw_piezo_coupling.i1.Y2( inclusion.coupling, v, phi ) = 0""" % omega_squared, '2' : """dw_diffusion.i1.Y( inclusion.dielectric, psi, phi ) + dw_piezo_coupling.i1.Y2( inclusion.coupling, u, psi ) = 0""", }
I think there should be a similar set of equations for the electrostatic-elastic coupling case, but the integration should be done over the interface of the elastic material and the vacuum regions. I think that can be handled by adding a 2D integral in the problem definition. What I'm missing is how to define the coupling term in SfePy. I think the form of the term is basically the same as dw_surface_ltr in sfepy_manual.pdf. Also, it looks like the gradient of the electric potential can be obtained from cachesBasic.
Is it possible to use the gradient of the electric potential information to calculate the Maxwell stress tensor using appropriate vector operations (divergence and dot product), along with the dielectric tensor? Would this be doable through the dw_surface_ltr term? Or is it better to implement a new coupling term (based off dw_surface_ltr) to handle this (and how would one go about that :) )?
I would probably implement a new term (called e.g. dw_electric_coupling?) by following the pattern of dw_piezo_coupling (sfepy/terms/termsPiezo.py), but working with surface integrals. Notice how dw_piezo_coupling treats its threefold nature - it can be called in the "grad" (as in eq. '1'), "div" (eq/ '2') and "eval" (both v, phi known) forms. The last one is not needed to solve problems. The grad, div terminology comes from the Stokes problem.
Ok, this makes sense and confirms my initial hunch. Is there any reference on the weak form derivation of the piezoelectric coupling term?
Hope everything is clear! Any advice is greatly appreciated!
I will try to assist you with the practical implementation, of course, but not much with the theory - you are the expert :).
Thank you very much as always. I know you have a lot of stuff going on, so it's pretty amazing and kind of you to answer my questions. :) I have some questions taking shape in my head, but for now I'll get some sleep and see if things make more sense tomorrow.
cheers, r.
Best, Logan
[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
Logan Sorenson wrote:
Hi Robert,
On Wed, Oct 28, 2009 at 4:23 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi Logan,
I am on a very thin ice with anything "electric" in its name, so do not take me too seriously... Let me ask some questions first:
Don't worry, I'm on thin ice with lots of things! :)
- do we need the edge (Nedelec) elements for this (in 3D)?
I haven't encountered the Nedelec elements before, so I'll do some more research to see if they are necessary. I'm hoping not :).
Fine! We are going to need them anyway, one day, but for the moment it's not on my TODO list.
- is the model similar to what is described in [2]? (I can access the fulltext.)
This is actually a really good paper. Thanks for pointing it out to me! [3] by the same author (found through the references of [2]) has some more detailed derivations. Also, her PhD thesis [4] is quite good. I think this is a better description of the problem than I could provide.
Heh, I just typed "electrostatic coupling FEM" to google, it was the first thing to jump up.
- do you have somewhere a complete weak formulation of the problem?
I haven't worked out the weak formulation yet. I think I need to consider some of the issues in [2-4], such as handling the changes in the electrical domain due to the displacements. I think the original construction below may work in the linear small displacement regime. By the way, do you know of a good reference which explains the weak form derivations? I got the basics from the sfepy_manual.pdf, but I think I need more practice for general things.
The basics are really that simple: you multiply your PDEs by some test functions, integrate over the whole domain, and use integration by parts on terms with the second order derivatives, if you do not like them ;) The black magic is in how to choose the right function spaces for the unknown and test functions, and how the weak solution relates to the strong one, but this was done already for the common problems. More can be found e.g. in [5].
Hi,
I was wondering how to implement the electrostatic-elastic coupling in SfePy. Basically, I think this coupling is described by the Maxwell stress tensor [1] which gives the stress due to the electrostatic field (setting the magnetic field to 0 in the quasistatic approximation) on the surface of a conductor.
I think the coupling can be handled in a similar manner with the piezoelectric coupling. For example, in input/piezo.py, the equations are:
equations = { '1' : """- %f * dw_mass_vector.i1.Y( inclusion.density, v, u ) + dw_lin_elastic_iso.i1.Y( inclusion.lam, inclusion.mu, v, u ) - dw_piezo_coupling.i1.Y2( inclusion.coupling, v, phi ) = 0""" % omega_squared, '2' : """dw_diffusion.i1.Y( inclusion.dielectric, psi, phi ) + dw_piezo_coupling.i1.Y2( inclusion.coupling, u, psi ) = 0""", }
I think there should be a similar set of equations for the electrostatic-elastic coupling case, but the integration should be done over the interface of the elastic material and the vacuum regions. I think that can be handled by adding a 2D integral in the problem definition. What I'm missing is how to define the coupling term in SfePy. I think the form of the term is basically the same as dw_surface_ltr in sfepy_manual.pdf. Also, it looks like the gradient of the electric potential can be obtained from cachesBasic.
Is it possible to use the gradient of the electric potential information to calculate the Maxwell stress tensor using appropriate vector operations (divergence and dot product), along with the dielectric tensor? Would this be doable through the dw_surface_ltr term? Or is it better to implement a new coupling term (based off dw_surface_ltr) to handle this (and how would one go about that :) )? I would probably implement a new term (called e.g. dw_electric_coupling?) by following the pattern of dw_piezo_coupling (sfepy/terms/termsPiezo.py), but working with surface integrals. Notice how dw_piezo_coupling treats its
Logan Sorenson wrote: threefold nature - it can be called in the "grad" (as in eq. '1'), "div" (eq/ '2') and "eval" (both v, phi known) forms. The last one is not needed to solve problems. The grad, div terminology comes from the Stokes problem.
Ok, this makes sense and confirms my initial hunch. Is there any reference on the weak form derivation of the piezoelectric coupling term?
You can find it e.g. in [6] (an article of my colleague E. Rohan), section 2, namely 2.4. It is within a broader setting (homogenization), so ignore everything with domain decomposition and epsilon in it :)
Hope everything is clear! Any advice is greatly appreciated! I will try to assist you with the practical implementation, of course, but not much with the theory - you are the expert :).
Thank you very much as always. I know you have a lot of stuff going on, so it's pretty amazing and kind of you to answer my questions. :) I have some questions taking shape in my head, but for now I'll get some sleep and see if things make more sense tomorrow.
You are welcome! It's good to have a break from the "must do" stuff fo a while.
cheers, r.
[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
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. :)
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].
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...
Also, would I have to implement some c code in sfepy/terms/extmods?
Thanks, and best regards! Logan
On Thu, Oct 29, 2009 at 2:24 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Logan Sorenson wrote:
Hi Robert,
On Wed, Oct 28, 2009 at 4:23 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi Logan,
I am on a very thin ice with anything "electric" in its name, so do not take me too seriously... Let me ask some questions first:
Don't worry, I'm on thin ice with lots of things! :)
- do we need the edge (Nedelec) elements for this (in 3D)? > I haven't encountered the Nedelec elements before, so I'll do some more research to see if they are necessary. I'm hoping not :).
Fine! We are going to need them anyway, one day, but for the moment it's not on my TODO list.
- is the model similar to what is described in [2]? (I can access the fulltext.)
This is actually a really good paper. Thanks for pointing it out to me! [3] by the same author (found through the references of [2]) has some more detailed derivations. Also, her PhD thesis [4] is quite good. I think this is a better description of the problem than I could provide.
Heh, I just typed "electrostatic coupling FEM" to google, it was the first thing to jump up.
- do you have somewhere a complete weak formulation of the problem?
I haven't worked out the weak formulation yet. I think I need to consider some of the issues in [2-4], such as handling the changes in the electrical domain due to the displacements. I think the original construction below may work in the linear small displacement regime. By the way, do you know of a good reference which explains the weak form derivations? I got the basics from the sfepy_manual.pdf, but I think I need more practice for general things.
The basics are really that simple: you multiply your PDEs by some test functions, integrate over the whole domain, and use integration by parts on terms with the second order derivatives, if you do not like them ;) The black magic is in how to choose the right function spaces for the unknown and test functions, and how the weak solution relates to the strong one, but this was done already for the common problems. More can be found e.g. in [5].
Hi,
I was wondering how to implement the electrostatic-elastic coupling in SfePy. Basically, I think this coupling is described by the Maxwell stress tensor [1] which gives the stress due to the electrostatic field (setting the magnetic field to 0 in the quasistatic approximation) on the surface of a conductor.
I think the coupling can be handled in a similar manner with the piezoelectric coupling. For example, in input/piezo.py, the equations are:
equations = { '1' : """- %f * dw_mass_vector.i1.Y( inclusion.density, v, u ) + dw_lin_elastic_iso.i1.Y( inclusion.lam, inclusion.mu, v, u ) - dw_piezo_coupling.i1.Y2( inclusion.coupling, v, phi ) = 0""" % omega_squared, '2' : """dw_diffusion.i1.Y( inclusion.dielectric, psi, phi ) + dw_piezo_coupling.i1.Y2( inclusion.coupling, u, psi ) = 0""", }
I think there should be a similar set of equations for the electrostatic-elastic coupling case, but the integration should be done over the interface of the elastic material and the vacuum regions. I think that can be handled by adding a 2D integral in the problem definition. What I'm missing is how to define the coupling term in SfePy. I think the form of the term is basically the same as dw_surface_ltr in sfepy_manual.pdf. Also, it looks like the gradient of the electric potential can be obtained from cachesBasic.
Is it possible to use the gradient of the electric potential information to calculate the Maxwell stress tensor using appropriate vector operations (divergence and dot product), along with the dielectric tensor? Would this be doable through the dw_surface_ltr term? Or is it better to implement a new coupling term (based off dw_surface_ltr) to handle this (and how would one go about that :) )? I would probably implement a new term (called e.g. dw_electric_coupling?) by following the pattern of dw_piezo_coupling (sfepy/terms/termsPiezo.py), but working with surface integrals. Notice how dw_piezo_coupling treats its
Logan Sorenson wrote: threefold nature - it can be called in the "grad" (as in eq. '1'), "div" (eq/ '2') and "eval" (both v, phi known) forms. The last one is not needed to solve problems. The grad, div terminology comes from the Stokes problem.
Ok, this makes sense and confirms my initial hunch. Is there any reference on the weak form derivation of the piezoelectric coupling term?
You can find it e.g. in [6] (an article of my colleague E. Rohan), section 2, namely 2.4. It is within a broader setting (homogenization), so ignore everything with domain decomposition and epsilon in it :)
Hope everything is clear! Any advice is greatly appreciated! I will try to assist you with the practical implementation, of course, but not much with the theory - you are the expert :).
Thank you very much as always. I know you have a lot of stuff going on, so it's pretty amazing and kind of you to answer my questions. :) I have some questions taking shape in my head, but for now I'll get some sleep and see if things make more sense tomorrow.
You are welcome! It's good to have a break from the "must do" stuff fo a while.
cheers, r.
[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
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
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 faces).
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".
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...
Thanks, and best regards! Logan
I hope that helps a bit, continue asking :)
r.
[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
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
Logan Sorenson wrote:
Hi Robert,
On Wed, Nov 4, 2009 at 10:49 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi Logan
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
Logan Sorenson wrote: 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. :)
Fine! Just a note: if you just test new things, try to use a really small mesh, as things are obviously much faster then. But I guess your mesh is ok in thus respect :) Also, you could maybe use the script/blockgen.py to generate the mesh:
./script/blockgen.py -d '[0.6, 0.1, 0.0015]' -s '[60, 10, 3]' -c '[0, 0, -0.00025]' -o out.{mesh or vtk}
The element groups of the layers would have to be hand-edited, though.
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.
I see. That should be doable, IMHO.
- 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?
It is mentioned on the other site [8]. That documentation got out-of-sync with the google site one. Now the only place should be the sphinx docs.
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.
I see. So you solve two problems on two subdomains with a common interface, where a coupling term plays.
I hope that helps a bit, continue asking :)
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!
So do I :)
r.
[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 [8] http://sfepy.kme.zcu.cz/index.cgi/Sfepy input file#Integrals
participants (2)
-
Logan Sorenson -
Robert Cimrman