understanding my first problem
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
I assume 'Omega' is a domain that includes the whole problem, but that I can also name domains and some how use different domains for different purposes. Is this correct?
What does the '3_4_P1' mean?
Does this section of the setup file determine what the output is that is dumped to the vtk output files? Is the variable 'u' what is ultimately being sought, reported, and visualized?
What are the units on 'u'?
======================
The next section I have questions about is this:
equations = { 'balance_of_forces' : """dw_lin_elastic_iso.i1.Omega( solid.lame, v, u ) = dw_surface_ltr.isurf.Top( traction.val, v )""", }
material_1 = { 'name' : 'solid', 'mode' : 'here', 'region' : 'Omega', 'lame' : {'lambda' : 1e1, 'mu' : 1e0}, # Lame coefficients. }
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it?
I think that is enough for now. Thanks for your ongoing help and patience.
Ryan
On Sat, Jul 12, 2008 at 11:03 PM, Ryan Krauss ryan...@gmail.com wrote:
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
I assume 'Omega' is a domain that includes the whole problem, but that I can also name domains and some how use different domains for different purposes. Is this correct?
Robert will answer the above when he comes to work tomorrow.
What does the '3_4_P1' mean?
3D element with 4 nodes and P1 is the polynomial approximation, linear in this case.
Does this section of the setup file determine what the output is that is dumped to the vtk output files? Is the variable 'u' what is
In my case (schroedinger.py), this is done by the actual script schroedinger.py. In your case, simple.py, I cannot find it though. I don't have time now to dig into that, so I'll leave this for Robert to answer.
ultimately being sought, reported, and visualized?
What are the units on 'u'?
======================
The next section I have questions about is this:
equations = { 'balance_of_forces' : """dw_lin_elastic_iso.i1.Omega( solid.lame, v, u ) = dw_surface_ltr.isurf.Top( traction.val, v )""", }
material_1 = { 'name' : 'solid', 'mode' : 'here', 'region' : 'Omega', 'lame' : {'lambda' : 1e1, 'mu' : 1e0}, # Lame coefficients. }
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it?
I think that is enough for now. Thanks for your ongoing help and patience.
Again couple more questions for Robert.
Robert, when you answer all those questions, could you please prepare a nice rst file (that we can put into the sphinx docs) that explains all of the above (assuming it's going to stay the current way in the future)? Because I was asking exactly the same questions too.
In my case, rather than answering the questions, I'd prefer pointers how I can easily find the answers myself, e.g. with the lambda and mu question, for me a pointer in the source code + couple comments would be enough. It think this will greatly improve our docs.
Thanks a lot, Ondrej
Ondrej Certik wrote:
Robert, when you answer all those questions, could you please prepare a nice rst file (that we can put into the sphinx docs) that explains all of the above (assuming it's going to stay the current way in the future)? Because I was asking exactly the same questions too.
There are many rst-like mark-ups. Have you some summary of the one we should use (i.e. others use with sphinx)?
In my case, rather than answering the questions, I'd prefer pointers how I can easily find the answers myself, e.g. with the lambda and mu question, for me a pointer in the source code + couple comments would be enough. It think this will greatly improve our docs.
doc/sfepy_manual.pdf + search the relevant term. All the basic ones should be already documented, if not, file a bugreport.
Concerning sources, all the terms are under sfepy/terms/, the files are named as in the manual.
r.
On Sat, 2008-07-12 at 16:03 -0500, Ryan Krauss wrote:
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Just to cut down on Robert's work, I'll try to reply to questions I can answer
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
u is the displacement sought after. v is the "dual" of u. Nitchze(?). Usually used to do a posteriori estimate of the error. This is interesting in that you can use this error info to do automatic refining of the mesh until the error is below some preset value. I would like to use this for navier-stokes solver to solve without any turbulence modeling.
I assume 'Omega' is a domain that includes the whole problem, but that I can also name domains and some how use different domains for different purposes. Is this correct?
yes
What does the '3_4_P1' mean?
Does this section of the setup file determine what the output is that is dumped to the vtk output files? Is the variable 'u' what is ultimately being sought, reported, and visualized?
What are the units on 'u'?
usually units in FE are anything that are "consistent".
======================
The next section I have questions about is this:
equations = { 'balance_of_forces' : """dw_lin_elastic_iso.i1.Omega( solid.lame, v, u ) = dw_surface_ltr.isurf.Top( traction.val, v )""", }
material_1 = { 'name' : 'solid', 'mode' : 'here', 'region' : 'Omega', 'lame' : {'lambda' : 1e1, 'mu' : 1e0}, # Lame coefficients. }
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
those sound like lame's constants, are related to E and nu. You can look them in any elasticity theory book.
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it?
traction loads are being applied to a portion of the boundary of Omega. They would be tangent to the surface/line. Units again whatever you choose but they hav eto be consistent with other units you have used.
Hope this helps.
-osman
osman wrote:
On Sat, 2008-07-12 at 16:03 -0500, Ryan Krauss wrote:
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Just to cut down on Robert's work, I'll try to reply to questions I can answer
Hey, great!
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
u is the displacement sought after. v is the "dual" of u. Nitchze(?). Usually used to do a posteriori estimate of the error. This is interesting in that you can use this error info to do automatic refining of the mesh until the error is below some preset value. I would like to use this for navier-stokes solver to solve without any turbulence modeling.
This is interesting - how you would use a test variable for a posteriori estimates? I know almost nothing about this, but am interested to have some error estimates in SfePy, preferably such that do not depend on the equations solved (i.e. discretization-error based). Have you some link?
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
those sound like lame's constants, are related to E and nu. You can look them in any elasticity theory book.
exactly.
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it?
traction loads are being applied to a portion of the boundary of Omega. They would be tangent to the surface/line. Units again whatever you choose but they hav eto be consistent with other units you have used.
...the portion of the boundary defined by the region 'Top'.
Hope this helps.
thanks for joining the discussion!
r.
That is funny. I have never taken a theory of elasticity class (actually I was in one for a day and dropped it) and had never heard of Lame (I don't know how to put the accent mark on the e). So, I thought from the comments in the original setup file that the coeffiecents were lame in this sense:
http://dictionary.reference.com/search?q=lame&x=0&y=0
namely:
- weak; inadequate; unsatisfactory; clumsy: a lame excuse. or
- Slang. out of touch with modern fads or trends; unsophisticated. –verb (used with object)
Ryan
On Mon, Jul 14, 2008 at 5:18 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote:
osman wrote:
On Sat, 2008-07-12 at 16:03 -0500, Ryan Krauss wrote:
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Just to cut down on Robert's work, I'll try to reply to questions I can answer
Hey, great!
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
u is the displacement sought after. v is the "dual" of u. Nitchze(?). Usually used to do a posteriori estimate of the error. This is interesting in that you can use this error info to do automatic refining of the mesh until the error is below some preset value. I would like to use this for navier-stokes solver to solve without any turbulence modeling.
This is interesting - how you would use a test variable for a posteriori estimates? I know almost nothing about this, but am interested to have some error estimates in SfePy, preferably such that do not depend on the equations solved (i.e. discretization-error based). Have you some link?
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
those sound like lame's constants, are related to E and nu. You can look them in any elasticity theory book.
exactly.
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it?
traction loads are being applied to a portion of the boundary of Omega. They would be tangent to the surface/line. Units again whatever you choose but they hav eto be consistent with other units you have used.
...the portion of the boundary defined by the region 'Top'.
Hope this helps.
thanks for joining the discussion!
r.
Ryan Krauss wrote:
That is funny. I have never taken a theory of elasticity class (actually I was in one for a day and dropped it) and had never heard of Lame (I don't know how to put the accent mark on the e). So, I thought from the comments in the original setup file that the coeffiecents were lame in this sense:
http://dictionary.reference.com/search?q=lame&x=0&y=0
namely:
- weak; inadequate; unsatisfactory; clumsy: a lame excuse. or
- Slang. out of touch with modern fads or trends; unsophisticated. –verb (used with object)
:-) I thought you might have thought something like that, but was lazy to write with diacritics. But now you know one more French guy!
http://en.wikipedia.org/wiki/Lam%C3%A9_parameters
r.
On Mon, Jul 14, 2008 at 3:37 PM, Robert Cimrman cimr...@ntc.zcu.cz wrote:
Ryan Krauss wrote:
That is funny. I have never taken a theory of elasticity class (actually I was in one for a day and dropped it) and had never heard of Lame (I don't know how to put the accent mark on the e). So, I thought from the comments in the original setup file that the coeffiecents were lame in this sense:
http://dictionary.reference.com/search?q=lame&x=0&y=0
namely:
- weak; inadequate; unsatisfactory; clumsy: a lame excuse. or
- Slang. out of touch with modern fads or trends; unsophisticated. –verb (used with object)
:-) I thought you might have thought something like that, but was lazy to write with diacritics. But now you know one more French guy!
BTW, you can use copy & paste to get the "é" letter from the wikipedia.
Ondřej
Ondrej Certik wrote:
On Mon, Jul 14, 2008 at 3:37 PM, Robert Cimrman cimr...@ntc.zcu.cz wrote:
Ryan Krauss wrote:
That is funny. I have never taken a theory of elasticity class (actually I was in one for a day and dropped it) and had never heard of Lame (I don't know how to put the accent mark on the e). So, I thought from the comments in the original setup file that the coeffiecents were lame in this sense:
http://dictionary.reference.com/search?q=lame&x=0&y=0
namely:
- weak; inadequate; unsatisfactory; clumsy: a lame excuse. or
- Slang. out of touch with modern fads or trends; unsophisticated. –verb (used with object) :-) I thought you might have thought something like that, but was lazy to write with diacritics. But now you know one more French guy!
BTW, you can use copy & paste to get the "é" letter from the wikipedia.
I can even type it easily, but I am so used to omit the accents (even in Czech, where there are tons of them) that I do not bother... My name has none, unlike yours :)
On Mon, 2008-07-14 at 12:18 +0200, Robert Cimrman wrote:
osman wrote:
On Sat, 2008-07-12 at 16:03 -0500, Ryan Krauss wrote:
This is interesting - how you would use a test variable for a posteriori estimates? I know almost nothing about this, but am interested to have some error estimates in SfePy, preferably such that do not depend on the equations solved (i.e. discretization-error based). Have you some link?
hmm.. I guess I was wrong, wrote that at 3am local time! Just replied without reading it closely. When I saw "dual" I thought you were talking about "adjoint" that is in one of your presentations (closed channel shape optimization). I was also looking at one of the Fenics documents, and they use dual/adjoint for estimating the error in some integrated quantity (in this case they were looking at lift/drag coeffs). And use it for auto refinements. I would like to do the same using sfepy.
br, -osman
osman wrote:
On Mon, 2008-07-14 at 12:18 +0200, Robert Cimrman wrote:
osman wrote:
On Sat, 2008-07-12 at 16:03 -0500, Ryan Krauss wrote:
This is interesting - how you would use a test variable for a posteriori estimates? I know almost nothing about this, but am interested to have some error estimates in SfePy, preferably such that do not depend on the equations solved (i.e. discretization-error based). Have you some link?
hmm.. I guess I was wrong, wrote that at 3am local time! Just replied without reading it closely. When I saw "dual" I thought you were talking about "adjoint" that is in one of your presentations (closed channel shape optimization). I was also looking at one of the Fenics documents, and they use dual/adjoint for estimating the error in some integrated quantity (in this case they were looking at lift/drag coeffs). And use it for auto refinements. I would like to do the same using sfepy.
I see. Yes, we use the adjoint problem to a weak formulation of Navier-Stokes equations to perform a shape sensitivity analysis.
It would be great to have error estimators in sfepy. If you have some specific wishes, submit an issue, please, at 'http://code.google.com/p/sfepy/issues/list'.
r.
Hi Ryan,
Ryan Krauss wrote:
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Yes a tutorial would certainly be very useful. I have to update the docs at http://ui505p06-mbs.ntc.zcu.cz/sfe/Documentation, and I have already updated the google page to point there, too.
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
In solving the elasticity problem, your only unknown are the displacements 'u'. However, in finite element context problems are usually written using the weak formulation, where you multiply the original equations by so-called test or virtual functions - this is 'v' in your file. A very basic intro can be found at http://en.wikipedia.org/wiki/Weak_formulation or http://magnet.atp.tuwien.ac.at/scholz/projects/diss/html/node6.html (a simple google search returned that).
Now having a weak but continuous formulation of your problem, the problem has to be discretized (i.e. the continuous functions, or fields approximated by a finite-dimensional approximation) so that it can solved on computer. This is what the finite element method does - it replaces the original continuous functions 'u' and 'v' by their finite-dimensional counterparts 'u_h' and 'v_h' (usual notation, h indicates dependence on a mesh element size). There are many ways how to do this discretization in the finite element (FE) context. In SfePy we use a simple piecewise-polynomial approximation, e.g. P1 (linear polynomials), P2 (second order (quadratic) polynomials) on a simplex and Q1 (bi-, or tri-linear polynomials) on rectangular elements. Now in SfePy you can have many variables, unknown, test, or other, that can use different (or share the same) approximations. In your example both 'u' and 'v' are approximated in the same way, by P1 polynomials on tetrahedral (4-node) elements in 3 space dimensions, hence '3_4_P1' in the 'field' entry (thinking about this now, I think a better name for 'field' would be 'function_space'...).
To conclude, I will comment the code below:
field_1 = { 'name' : 'displacement', # Name of the field approximation, referred to below in the variables (can be arbitrary). 'dim' : (3,1), # dimension: (1,1) for scalar fields, (3,1) for vector fields in 3D, (2,1) in 2D. 'domain' : 'Omega', # Domain where the field (and subsequently of all variables of this field) is defined. 'bases' : {'Omega' : '3_4_P1'} # FE function spaces to use on subdomains of Omega, usually only one is used on the whole domain, just like here. }
variable_1 = { 'name' : 'u', # Name to use in 'equations' (can be arbitrary). 'kind' : 'unknown field', # Type: this is the unknown we wish to compute. 'field' : 'displacement', # It uses the FE approximation defined in the field called 'displacement'. 'order' : 0, # Order in the vector where all the unknown are concatenated together. Ignored here, as there is only one unknown. }
variable_2 = { 'name' : 'v', # As above. 'kind' : 'test field', # Test field (not an unkown). 'field' : 'displacement', # Uses the same FE approximation as 'u', i.e. 'displacement'. 'dual' : 'u', # It is a test field related to 'u'. }
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
In a weak formulation of a PDE, the equation must hold for all possible test functions from a given function space, see the links above and your books on finite elements.
I assume 'Omega' is a domain that includes the whole problem, but that
Yes.
I can also name domains and some how use different domains for different purposes. Is this correct?
Yes. It is the name in 'region*' entries.
What does the '3_4_P1' mean?
See above.
Does this section of the setup file determine what the output is that is dumped to the vtk output files? Is the variable 'u' what is ultimately being sought, reported, and visualized?
The variable name ('u') is the one you see in Paraview, yes. By default, 'simple.py' names the output file by the mesh file you use. It can be changed by '-o' option of simple.py, see '$ ./simple.py --help'
What are the units on 'u'?
As osman answered, it is up to you how you interpret you data: e.g. if you inputs are in: length: mm time: s mass: kg then forces are in mN and stresses (and pressure) in kPa, so that the unit set is consistent.
======================
The next section I have questions about is this:
equations = { 'balance_of_forces' : """dw_lin_elastic_iso.i1.Omega( solid.lame, v, u ) = dw_surface_ltr.isurf.Top( traction.val, v )""", }
material_1 = { 'name' : 'solid', 'mode' : 'here', 'region' : 'Omega', 'lame' : {'lambda' : 1e1, 'mu' : 1e0}, # Lame coefficients. }
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
Those are the Lame coefficients, see http://en.wikipedia.org/wiki/Lam%C3%A9_coefficients Density is not specified, as the input file defines a quasistatic problem only (no inertial forces).
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
Yes.
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it
It can be an arbitrary function of space and time, and is defined in the region specified in 'dw_surface_ltr.isurf.Top', i.e. 'Top'. So it is not defined at all the nodes, but only at the nodes in 'Top'.
Units: consistent with other units, see the remark above. Google says this: http://hyperphysics.phy-astr.gsu.edu/Hbase/units.html
compressive or tensile? What direction is associated with it?
see 'dw surface ltr' in doc/sfepy_manual.pdf. As tractionLoad() returns a scalar at each node of Top, is corresponds to specifying a pressure traction on the top surface. About the sign, I am never sure - look at the resulting deformation :)
I think that is enough for now. Thanks for your ongoing help and patience.
It is good to know what should go into the documentaion :-)
r.
So, I think I am finally understanding this problem. I worked out a very long winded solution, using FEA and Rayleigh-Ritz as well as an analytical solution that I put here if anyone wants to check it out and tell me if I am really getting this: http://www.siue.edu/~rkrauss/sfepy/axial_rod_2_out.pdf
Attached is a png showing my 1D analytical solution vs. the z axis displacement from SfePy. I think they look pretty close (SfePy is a little stiffer, but I think that is reasonable). The comparison is also between a 3D and a 1D model, so I think perfect agreement is unreasonable to expect.
If you have any input, please let me know. Otherwise I will now focus on understanding how the linear material model is implemented in SfePy in preparation to create a nonlinear material.
Ryan
On Sat, Jul 12, 2008 at 4:03 PM, Ryan Krauss ryan...@gmail.com wrote:
So, everything is correctly installed and thanks to Robert I can now work in the directory of my choice. Basically, things are working. Now I need to understand what the code is actually doing.
My intention is to ask some questions about the attached setup file, increase my understanding, add some comments to, and possibly have a good tutorial problem when it's done.
Referring to the attached file, which gets passed to simple.py as an input, I have the following questions:
The file starts with
field_1 = { 'name' : 'displacement', 'dim' : (3,1), 'domain' : 'Omega', 'bases' : {'Omega' : '3_4_P1'} }
variable_1 = { 'name' : 'u', 'kind' : 'unknown field', 'field' : 'displacement', 'order' : 0, } variable_2 = { 'name' : 'v', 'kind' : 'test field', 'field' : 'displacement', 'dual' : 'u', }
All 3 of these are displacements, I am trying to understand how they are related or different. I assume a field is a vector. Maybe the variables are scalars. If 'u' is one component of the vector field displacement, where is the degree of freedom specified (i.e. which component of 'displacement' does 'u' refer to)?
From the description, it sounds like 'v' is a test field of 'u'. Does that mean that when 'v' converges, 'u' is the result? Or something else?
I assume 'Omega' is a domain that includes the whole problem, but that I can also name domains and some how use different domains for different purposes. Is this correct?
What does the '3_4_P1' mean?
Does this section of the setup file determine what the output is that is dumped to the vtk output files? Is the variable 'u' what is ultimately being sought, reported, and visualized?
What are the units on 'u'?
======================
The next section I have questions about is this:
equations = { 'balance_of_forces' : """dw_lin_elastic_iso.i1.Omega( solid.lame, v, u ) = dw_surface_ltr.isurf.Top( traction.val, v )""", }
material_1 = { 'name' : 'solid', 'mode' : 'here', 'region' : 'Omega', 'lame' : {'lambda' : 1e1, 'mu' : 1e0}, # Lame coefficients. }
What are lambda and mu and what are their units? Youngs modulus and mass density per unit length? Something else? I assume it is something mass related and something stiffness related.
It seems like equations is setting up a force balance on the top of the material specimen where the load is being applied. Is this correct?
It seems like this function:
def tractionLoad( ts, coor, region, ig ): """ts = TimeStepper, coor = coordinates of field nodes in region.""" import numpy as nm from sfepy.base.base import pause, debug nt = ts.nt
val = nm.zeros( (coor.shape[0],), dtype = nm.float64 ) val.fill( 1e-1 * nt )
print val return {'val' : val}
defines a traction load at all nodes (or maybe on some infinitesimal volume or area) that is the same at all points and ramping up in time with a slope of 0.1 per time step. Can this force be an arbitrary function of time? What are the units of the traction load? Is it compressive or tensile? What direction is associated with it?
I think that is enough for now. Thanks for your ongoing help and patience.
Ryan
On Thu, Jul 31, 2008 at 8:42 PM, Ryan Krauss ryan...@gmail.com wrote:
So, I think I am finally understanding this problem. I worked out a very long winded solution, using FEA and Rayleigh-Ritz as well as an analytical solution that I put here if anyone wants to check it out and tell me if I am really getting this: http://www.siue.edu/~rkrauss/sfepy/axial_rod_2_out.pdf
Wow, what a long file! Will you put it contribute it to some sfepy tutorial?
For me it's unintuitive those funny {E} symbols, but I guess that's the standard notation in FEM. :)
Ondrej
I am glad to make it into some sort of tutorial, if someone can verify that it is right and agree that it is useful. I will also gladly add pictures and additional explanations or remove stuff.
On Thu, Jul 31, 2008 at 2:12 PM, Ondrej Certik ond...@certik.cz wrote:
On Thu, Jul 31, 2008 at 8:42 PM, Ryan Krauss ryan...@gmail.com wrote:
So, I think I am finally understanding this problem. I worked out a very long winded solution, using FEA and Rayleigh-Ritz as well as an analytical solution that I put here if anyone wants to check it out and tell me if I am really getting this: http://www.siue.edu/~rkrauss/sfepy/axial_rod_2_out.pdf
Wow, what a long file! Will you put it contribute it to some sfepy tutorial?
For me it's unintuitive those funny {E} symbols, but I guess that's the standard notation in FEM. :)
Ondrej
On Thu, Jul 31, 2008 at 9:29 PM, Ryan Krauss ryan...@gmail.com wrote:
I am glad to make it into some sort of tutorial, if someone can verify that it is right and agree that it is useful. I will also gladly add pictures and additional explanations or remove stuff.
I can verify that it is useful. Robert can check that it is also right. :)
Ondrej
And we can modify the symbols and terms to be whatever we want. I was mainly trying to follow the conventions of the book by Cook et.al. that I am using.
On Thu, Jul 31, 2008 at 3:23 PM, Ondrej Certik ond...@certik.cz wrote:
On Thu, Jul 31, 2008 at 9:29 PM, Ryan Krauss ryan...@gmail.com wrote:
I am glad to make it into some sort of tutorial, if someone can verify that it is right and agree that it is useful. I will also gladly add pictures and additional explanations or remove stuff.
I can verify that it is useful. Robert can check that it is also right. :)
Ondrej
On Thu, Jul 31, 2008 at 10:29 PM, Ryan Krauss ryan...@gmail.com wrote:
And we can modify the symbols and terms to be whatever we want. I was mainly trying to follow the conventions of the book by Cook et.al. that I am using.
robert should say what he likes. My background is in theoretical physics, so I am used to using different kind of symbols.
Ondrej
Hi Ryan,
Ondrej Certik wrote:
On Thu, Jul 31, 2008 at 10:29 PM, Ryan Krauss ryan...@gmail.com wrote:
And we can modify the symbols and terms to be whatever we want. I was mainly trying to follow the conventions of the book by Cook et.al. that I am using.
robert should say what he likes. My background is in theoretical physics, so I am used to using different kind of symbols.
I will print it and go through it during this weekend, ok?
On the first glance it looks ok. It uses one of the engineering notations (matrices denoted by [], vectors by {}, right?).
I would prefer the notation used in doc/sfepy_manual.pdf, both indicial and tensor are ok.
For example, the first term in your (1) would look like:
continous:
1/2 \int E_{ijkl} \varepsilon_{ij} \varepsilon_{kl} or 1/2 \int E_(4) : \ull{\varepsilon} \ull{\varepsilon} here \ull is a shortcut for double underline (denoting a second-order tensor) and (4) denotes a fourth-order tensor (4 underlines are ugly)
discretized: 1/2 \bm{u}^T \bm{K} \bm{u} (u = displacements)
\bm{K} is global stiffness matrix, \bm{K}|_e contribution of element e
(4) would be: \ul{u} \approx \bm{\varphi} \bm{u} \nabla \ul{u} \approx \bm{G} \bm{u}, i.e. G denotes \nabla \varphi (gradients of the FE base funtions)
Anyway, it looks fine as it is, I have proposed the above just in case you have enough time...
Some figure would be nice, though.
cheers, r.
ps: an interesting link I have just found: http://www.sv.vt.edu/classes/ESM4714/methods/EEG.html
We can work out the notations details and make sure it agrees with other SfePy docs for consistency.
matrices denoted by [], vectors by {}
Yes.
But I think a matrix is a second order tensor and a vector is a first order tensor so that the first equation should be:
1/2 \int E_{ij} \varepsilon_{i} \varepsilon_{j}
or something like that - at least for the 1D problem.
Am I thinking correctly?
Ryan
On Fri, Aug 1, 2008 at 4:33 AM, Robert Cimrman cimr...@ntc.zcu.cz wrote:
Hi Ryan,
Ondrej Certik wrote:
On Thu, Jul 31, 2008 at 10:29 PM, Ryan Krauss ryan...@gmail.com wrote:
And we can modify the symbols and terms to be whatever we want. I was mainly trying to follow the conventions of the book by Cook et.al. that I am using.
robert should say what he likes. My background is in theoretical physics, so I am used to using different kind of symbols.
I will print it and go through it during this weekend, ok?
On the first glance it looks ok. It uses one of the engineering notations (matrices denoted by [], vectors by {}, right?).
I would prefer the notation used in doc/sfepy_manual.pdf, both indicial and tensor are ok.
For example, the first term in your (1) would look like:
continous:
1/2 \int E_{ijkl} \varepsilon_{ij} \varepsilon_{kl} or 1/2 \int E_(4) : \ull{\varepsilon} \ull{\varepsilon} here \ull is a shortcut for double underline (denoting a second-order tensor) and (4) denotes a fourth-order tensor (4 underlines are ugly)
discretized: 1/2 \bm{u}^T \bm{K} \bm{u} (u = displacements)
\bm{K} is global stiffness matrix, \bm{K}|_e contribution of element e
(4) would be: \ul{u} \approx \bm{\varphi} \bm{u} \nabla \ul{u} \approx \bm{G} \bm{u}, i.e. G denotes \nabla \varphi (gradients of the FE base funtions)
Anyway, it looks fine as it is, I have proposed the above just in case you have enough time...
Some figure would be nice, though.
cheers, r.
ps: an interesting link I have just found: http://www.sv.vt.edu/classes/ESM4714/methods/EEG.html
Ryan Krauss wrote:
We can work out the notations details and make sure it agrees with other SfePy docs for consistency.
matrices denoted by [], vectors by {}
Yes.
But I think a matrix is a second order tensor and a vector is a first order tensor so that the first equation should be:
1/2 \int E_{ij} \varepsilon_{i} \varepsilon_{j}
or something like that - at least for the 1D problem.
Am I thinking correctly?
There are two possible approaches: engineering (using just vectors and matrices) and "mathematical". The engineering approach employs usually symmetry, consider (in a pseudo-latex):
e_ij = 1/2 (du_i/dx_j + du_j/dx_i) - a symmetric second order tensor can be stored as {e} = [e_11, e_22, e_33, e_12, e_13, e_23] - a vector.
I would prefer working with proper tensors in the docs. The storage is just an implementation detail (yes, in sfepy we use this storage too).
r.
participants (4)
-
Ondrej Certik -
osman -
Robert Cimrman -
Ryan Krauss