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.