quick questions

[David Mashburn]

Hello again,

I am performing elasticity simulations with boundary forces. I can tell
sfepy is capable of this, and starting with the linear_elasticity.py
example has already provided a close approximation of what I need, but I
am not clear on a couple of things.

My main two related questions right now are:

1. Can sfepy apply forces to boundaries instead of to volumes
(effectively a von Neumann boundary condition)? What is the prescribed
way to do this? I did notice that there is a surface term similar to the
volume term, but I was not sure how to use it...

2. Also, what exactly are the "regions" in sfepy? Are they collections
of individual mesh elements? If so, is there any way of specifying
borders instead? Also, is there a good way way to visualize regions
and/or get at the raw data that makes them up? (I am having trouble
getting just the)

Thanks,
-David

BTW, I was not as familiar with quadrature, so I made a small visualizer
for the different quadrature tables (I know the formatting will probably
get eaten, but you should still get the idea):

import numpy as np
import matplotlib.pyplot as plt
from mayavi import mlab
import sfepy.fem
quadType='3_8'
order=5

qp=sfepy.fem.quadratures.quadrature_tables[quadType][order]
c=qp.coors
w=qp.weights
print c
print w

plt.cla()
mlab.clf()

if quadType=='1_2':
plt.plot(qp.bounds,[0,0])
for i in range(len(c)):
plt.plot(c[i],0,'k.',ms=w[i]*100)

elif quadType=='2_3':
plt.plot([0,1,0,0],[0,0,1,0])
for i in range(len(c)):
plt.plot(c[i][0],c[i][1],'k.',ms=w[i]*100)

elif quadType=='2_4':
plt.plot([0,1,1,0,0],[0,0,1,1,0])
for i in range(len(c)):
plt.plot(c[i][0],c[i][1],'k.',ms=w[i]*100)

elif quadType=='3_4':
mlab.plot3d([0,0,0], [1,0,0], [0,0,1], color=(0,0,1),tube_radius=0.025)
mlab.plot3d([0,1,0,0,1], [0,0,1,0,0], [0,0,0,1,0], color=(0,0,1),
tube_radius=0.025)
for i in range(len(c)):
#plt.plot(c[i][0],c[i][1],'k.',ms=w[i]*100)
print [c[i][0]],[c[i][1]],[c[i][2]],w[i]

mlab.points3d([c[i][0]],[c[i][1]],[c[i][2]],[w[i]],color=(1,1,1),scale_factor=2)

elif quadType=='3_8':
pts = [[i,j,k] for i in [0,1] for j in [0,1] for k in [0,1]]
lines = [[i,j] for i in pts for j in pts if
(pts.index(i)<pts.index(j) and sum(abs(np.array(i)-j))==1)]
for l in lines:
l=np.transpose(l)
mlab.plot3d(l[0], l[1], l[2], color=(0,0,1),tube_radius=0.025)
for i in range(len(c)):
#plt.plot(c[i][0],c[i][1],'k.',ms=w[i]*100)
print [c[i][0]],[c[i][1]],[c[i][2]],w[i]

mlab.points3d([c[i][0]],[c[i][1]],[c[i][2]],[w[i]],color=(1,1,1),scale_factor=2)

[Robert Cimrman]

On Fri, 2 Mar 2012, David Mashburn wrote:

> Hello again,
>
> I am performing elasticity simulations with boundary forces. I can tell
> sfepy is capable of this, and starting with the linear_elasticity.py example
> has already provided a close approximation of what I need, but I am not clear
> on a couple of things.
>
> My main two related questions right now are:
>
> 1. Can sfepy apply forces to boundaries instead of to volumes (effectively a
> von Neumann boundary condition)? What is the prescribed way to do this? I did
> notice that there is a surface term similar to the volume term, but I was not
> sure how to use it...

Look at linear_elastic_tractions.py example in the same directory as the
linear elastic one.

> 2. Also, what exactly are the "regions" in sfepy? Are they collections of
> individual mesh elements? If so, is there any way of specifying borders
> instead? Also, is there a good way way to visualize regions and/or get at the
> raw data that makes them up? (I am having trouble getting just the)

The regions can define any mesh entity - nodes, elements, etc. if it
contains e.g. only surface nodes, it defines the surface faces and can be
used for traction forces. Check the user`s guide section about problem
description file.

> Thanks,
> -David
>
> BTW, I was not as familiar with quadrature, so I made a small visualizer for
> the different quadrature tables (I know the formatting will probably get
> eaten, but you should still get the idea):

Thanks, I will check it when I get to a computer.

Cheers,
r.

[Robert Cimrman]

Hi again!

On 03/03/2012 04:20 PM, Robert Cimrman wrote:
> On Fri, 2 Mar 2012, David Mashburn wrote:
>
>> Hello again,
>>
>> I am performing elasticity simulations with boundary forces. I can tell sfepy
>> is capable of this, and starting with the linear_elasticity.py example has
>> already provided a close approximation of what I need, but I am not clear on
>> a couple of things.
>>
>> My main two related questions right now are:
>>
>> 1. Can sfepy apply forces to boundaries instead of to volumes (effectively a
>> von Neumann boundary condition)? What is the prescribed way to do this? I did
>> notice that there is a surface term similar to the volume term, but I was not
>> sure how to use it...
>
> Look at linear_elastic_tractions.py example in the same directory as the linear
> elastic one.

More specifically, see [1] and [2] .

>> 2. Also, what exactly are the "regions" in sfepy? Are they collections of
>> individual mesh elements? If so, is there any way of specifying borders
>> instead? Also, is there a good way way to visualize regions and/or get at the
>> raw data that makes them up? (I am having trouble getting just the)
>
> The regions can define any mesh entity - nodes, elements, etc. if it contains
> e.g. only surface nodes, it defines the surface faces and can be used for
> traction forces. Check the user`s guide section about problem description file.

See [3]. To visualize the regions, the best option is to use:

./simple.py <your_input> --save-regions-as-groups --solve-not

and view the created file <mesh name>_regions.vtk with postproc.py.
The '--solve-not' option is useful if the input is not finished (e.g. missing
equations, conditions, ...) and you just need to verify that the regions are
correctly defined.

To see the actual data, the easiest thing is to set a breakpoint in simple.py,
line 123 and explore 'print app.problem.domain.regions'. See also [4].

>> Thanks,
>> -David
>>
>> BTW, I was not as familiar with quadrature, so I made a small visualizer for
>> the different quadrature tables (I know the formatting will probably get
>> eaten, but you should still get the idea):
>
> Thanks, I will check it when I get to a computer.

Nice, do you think it would be useful to have the images in the docs?

r.

[1]
http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tractions.html
[2]
http://sfepy.org/doc-devel/src/sfepy/terms/termsSurface.html#module-sfepy.terms.termsSurface
[3] http://sfepy.org/doc-devel/users_guide.html#regions
[4] http://sfepy.org/doc-devel/tutorial.html#interactive-example-linear-elasticity

[David Mashburn]

Thanks Robert.

On 03/05/2012 05:33 AM, Robert Cimrman wrote:
> Hi again!
>
> On 03/03/2012 04:20 PM, Robert Cimrman wrote:
>> On Fri, 2 Mar 2012, David Mashburn wrote:
>>
>>> Hello again,
>>>
>>> I am performing elasticity simulations with boundary forces. I can
>>> tell sfepy
>>> is capable of this, and starting with the linear_elasticity.py
>>> example has
>>> already provided a close approximation of what I need, but I am not
>>> clear on
>>> a couple of things.
>>>
>>> My main two related questions right now are:
>>>
>>> 1. Can sfepy apply forces to boundaries instead of to volumes
>>> (effectively a
>>> von Neumann boundary condition)? What is the prescribed way to do
>>> this? I did
>>> notice that there is a surface term similar to the volume term, but
>>> I was not
>>> sure how to use it...
>>
>> Look at linear_elastic_tractions.py example in the same directory as
>> the linear
>> elastic one.
>
> More specifically, see [1] and [2] .

Ok, great!

>
>
>>> 2. Also, what exactly are the "regions" in sfepy? Are they
>>> collections of
>>> individual mesh elements? If so, is there any way of specifying borders
>>> instead? Also, is there a good way way to visualize regions and/or
>>> get at the
>>> raw data that makes them up? (I am having trouble getting just the)
>>
>> The regions can define any mesh entity - nodes, elements, etc. if it
>> contains
>> e.g. only surface nodes, it defines the surface faces and can be used
>> for
>> traction forces. Check the user`s guide section about problem
>> description file.
>
> See [3]. To visualize the regions, the best option is to use:
>
> ./simple.py <your_input> --save-regions-as-groups --solve-not
>
> and view the created file <mesh name>_regions.vtk with postproc.py.
> The '--solve-not' option is useful if the input is not finished (e.g.
> missing equations, conditions, ...) and you just need to verify that
> the regions are correctly defined.
>
> To see the actual data, the easiest thing is to set a breakpoint in
> simple.py, line 123 and explore 'print app.problem.domain.regions'.
> See also [4].
>

Excellent, these are exactly what I was looking for, I'll try them out.

>>> Thanks,
>>> -David
>>>
>>> BTW, I was not as familiar with quadrature, so I made a small
>>> visualizer for
>>> the different quadrature tables (I know the formatting will probably
>>> get
>>> eaten, but you should still get the idea):
>>
>> Thanks, I will check it when I get to a computer.
>
> Nice, do you think it would be useful to have the images in the docs?
>

???
You know your audience better than I do; the only advantage I see (other
than for people rusty on quadrature) is for people to quickly see which
element types are supported, but you already show that.
Feel free to use code and/or images if you think it would help, though.

[Robert Cimrman]

On 03/06/2012 01:28 AM, David Mashburn wrote:
>>>> BTW, I was not as familiar with quadrature, so I made a small visualizer for
>>>> the different quadrature tables (I know the formatting will probably get
>>>> eaten, but you should still get the idea):
>>>
>>> Thanks, I will check it when I get to a computer.
>>
>> Nice, do you think it would be useful to have the images in the docs?
>>
> ???
> You know your audience better than I do; the only advantage I see (other than
> for people rusty on quadrature) is for people to quickly see which element
> types are supported, but you already show that.

Ok. I thought that if it was interesting for you, it might be interesting for
others.

> Feel free to use code and/or images if you think it would help, though.

I would like to add it to a new script/plot_quadrature.py - but for that the
quad_type and order should be given on command line (using optparse) and I
would prefer using matplotlib even for 3D plots [1], as Mayavi does not
interact well with matplotlib on some systems (mine included), and for drawing
a bunch of points it's an overkill. I do not have time to do it, and it seems
like an easy task for anyone, so patches welcome!

r.

[1] http://matplotlib.sourceforge.net/mpl_toolkits/mplot3d/tutorial.html

[David Mashburn]

I'm attaching a start to this script... I'm not really familiar with
OptionParser, but wanted to do something with this while it was fresh on
my mind. Either I or someone else can clean it up later...

Thanks,
-David

[Robert Cimrman]

Thanks, David.

I will create an enhancement issue and attach your script to it. It's a good start.

r.