In examples/large_deformation/hyperelastic.py a rotation by displacements is applied. By using a similar function the vectors defining the force couples could be defined for dw_surface_ltr (IMHO). Does it make sense?
r.
----- Reply message -----
From: "Andre Smit" <freev...(a)gmail.com>
To: <sfepy...(a)googlegroups.com>
Subject: Torque
Date: Sat, Dec 18, 2010 05:10
What is the best way to apply a torque load to a model?
--
Andre
--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
I am currrently looking for FEM packages to help me solve a system of
beams and columns, basically a collection of 1D bernoulli/timoshenko
line elements.
I started reading SfePy docs and i am getting the idea that doing the
above is not really possible here, am i right?
Are only 2D area elements permitted in SfePy?
Or is there any direct support for solving 1D line elements too..
Cheers
Nimish
I'm working on modeling a next-generation X-ray mirror for which the
shape can be actively controlled by use of many thin piezo-electric
actuators mounted on the mirror surface. The mirror is basically a
glass conical paraboloid with a 1 meter radius and 200 micron
thickness (e.g. http://en.wikipedia.org/wiki/X-ray_optics). Our
project is currently using a proprietary FEA package, but the model
setup and turnaround time is slow, in part because there is only one
part-time engineer who can run it.
SfePy looks like a great package and we're hoping that it could be
used to automate running a large number of different cases. I've
spent some time reading the documentation but I have a few questions
that I hope can be answered before going too much further. I want to
apologize in advance if some of my wording is imprecise, I have a
physics background but this topic is a bit outside my realm...
- Is SfePy appropriate for this problem?
- If a specify a grid with about 800 x 400 points (azimuthal, axial)
and about 10 boundary conditions (corresponding to mount points), what
is the rough order of magnitude of time to compute the solution? Is
it seconds, minutes, hours, or days?
- The linear elastic examples show a problem with a specified
displacement. How do I specify an input force? The piezo essentially
provides a tensile force along the surface.
- Is there a way to specify the problem and solve in cylindrical
coordinates? This is the natural coordinate system.
- How do I specify 6-DOF constraints which correspond to the mirror
mounts?
Thanks in advance for any help!
Tom Aldcroft
This selects a single node in an element group. It's useful if, for example, pressure is determined up to a constant, and you want to fix it.
r.
----- Reply message -----
From: "Steve Spicklemire" <st...(a)spvi.com>
To: <sfepy...(a)googlegroups.com>
Cc: "Steve Spicklemire" <st...(a)spvi.com>
Subject: Cylindrical coordinates
Date: Sun, Jun 10, 2012 15:14
Do you know.. how do most folks use 'node of group x' syntax?
thanks,
-steve
On Jun 10, 2012, at 7:06 AM, Steve Spicklemire wrote:
> OK.. no prob. I can write equations for each boundary... just thought there might be an easier way. ;-)
>
> -steve
>
> On Jun 10, 2012, at 6:39 AM, Robert Cimrman wrote:
>
>> Hi Steve,
>>
>> I cannotbhelp you much with gmsh - try asking the gmsh authors directly, if it is possible at all to export nodal group ids.
>>
>> r.
>>
>> ----- Reply message -----
>> From: "Steve Spicklemire" <st...(a)spvi.com>
>> To: <sfepy...(a)googlegroups.com>
>> Cc: "Steve Spicklemire" <st...(a)spvi.com>
>> Subject: Cylindrical coordinates
>> Date: Sat, Jun 9, 2012 15:01
>>
>>
>> Hi Richard,
>>
>> I don't *think* that's the problem.
>>
>> What follows is a ridiculously long winded explanation of the fact that I don't know how to specify groups of nodes in a .geo file. I've played around with various possibilities (Physical Line, Physical Point, etc.), but none seems to do the trick. My current plan is to start with a .geo template file that my app will render with particular values depending on user-prescribed variables. Once a .geo file is saved on disk, I can run gmsh/mesh_to_vtk.py and then sfepy. Anyway.. sorry for the verbosity.. but here it goes!
>>
>> -------------
>>
>> When I specify the BCs this way (case 1) everything works:
>>
>> regions = {
>> 'Omega' : ('all', {}),
>> 'Omega_Source' : ('nodes in (y<0.15)', {}),
>> 'Omega_Target' : ('nodes in (y>0.25)', {}),
>> }
>>
>> But this is awkward since I'll have some potentially complicated geometries. I'd rather "tag" some surfaces/lines as belonging to a particular group and using the syntax (case 2):
>>
>> regions = {
>> 'Omega' : ('all', {}),
>> 'Omega_Source' : ('nodes of group 1', {}),
>> 'Omega_Target' : ('nodes of group 2', {}),
>> }
>>
>> Where group 1 and group 2 are set up in my original .geo file.
>>
>> (BTW I think the socket api for gmsh is only for connecting a solver or a postprocessor using gmsh as the UI, not for creating meshes from outside applications)
>>
>> Let's say I start with a drastically simple .geo file like so:
>>
>> lc = 0.1;
>>
>> Point(1) = {0, 0, 0, lc};
>> Point(2) = {.1, 0, 0, lc} ;
>> Point(3) = {.1, .3, 0, lc} ;
>> Point(4) = {0, .3, 0, lc} ;
>>
>> Line(1) = {1,2} ;
>> Line(2) = {3,2} ;
>> Line(3) = {3,4} ;
>> Line(4) = {4,1} ;
>>
>> Line Loop(5) = {4,1,-2,3} ;
>>
>> Plane Surface(6) = {5} ;
>>
>> Physical Point(1) = {1,2} ;
>> Physical Point(2) = {3,4} ;
>>
>> Physical Surface(7) = {6};
>>
>> (I'm hoping that "Physical Point(1)" might correspond to "group 1" in the region definition (case 2) and "Physical Point(2)" may be "group 2". No luck on that I'm afraid.)
>>
>> I use gmsh to convert it to a mesh:
>>
>> aluminum:sfepy steve$ gmsh -2 t1.geo -format mesh
>> Info : Running 'gmsh -2 t1.geo -format mesh' [1 node(s), max. 1 thread(s)]
>> Info : Started on Sat Jun 9 06:33:30 2012
>> Info : Reading 't1.geo'...
>> Info : Done reading 't1.geo'
>> Info : Meshing 1D...
>> Info : Meshing curve 1 (Line)
>> Info : Meshing curve 2 (Line)
>> Info : Meshing curve 3 (Line)
>> Info : Meshing curve 4 (Line)
>> Info : Done meshing 1D (0.000122 s)
>> Info : Meshing 2D...
>> Info : Meshing surface 6 (Plane, MeshAdapt)
>> Info : Done meshing 2D (0.001485 s)
>> Info : 6 vertices 14 elements
>> Info : Writing 't1.mesh'...
>> Info : Done writing 't1.mesh'
>> Info : Stopped on Sat Jun 9 06:33:30 2012
>>
>> Now I've got the following, but no hints of group 1 or group 2 ;-(.
>>
>> MeshVersionFormatted 1
>> Dimension
>> 3
>> Vertices
>> 6
>> 0 0 0 1
>> 0.1 0 0 2
>> 0.1 0.3 0 3
>> 0 0.3 0 4
>> 0.1 0.15000000000017 0 2
>> 0 0.15000000000017 0 4
>> Triangles
>> 4
>> 1 2 6 6
>> 2 5 6 6
>> 3 4 6 6
>> 5 3 6 6
>> End
>>
>> Sadly, sfepy see's this as 3D.
>>
>> aluminum:sfepy steve$ ./dump-info.py t1.mesh
>> sfepy: reading mesh (t1.mesh)...
>> sfepy: ...done in 0.00 s
>> Mesh:t1
>> conns:
>> list: [array([[0, 1, 5],
>> [1, 4, 5],
>> [2, 3, 5],
>> [4, 2, 5]])]
>> coors:
>> (6, 3) array of float64
>> descs:
>> list: ['2_3']
>> dim:
>> 3
>> dims:
>> list: [2]
>> el_offsets:
>> (2,) array of int32
>> io:
>> None
>> mat_ids:
>> list: [array([6, 6, 6, 6])]
>> n_e_ps:
>> (1,) array of int32
>> n_el:
>> 4
>> n_els:
>> (1,) array of int32
>> n_nod:
>> 6
>> name:
>> t1
>> ngroups:
>> (6,) array of float64
>> nodal_bcs:
>> dict with keys: []
>> setup_done:
>> 0
>>
>> (BTW dump_info is just the following script:
>>
>> aluminum:sfepy steve$ cat dump-info.py
>> #!/usr/bin/env python
>>
>> from sfepy.fem import Mesh
>> import sys
>>
>> m = Mesh.from_file(sys.argv[1])
>>
>> print m
>> )
>>
>> I can use the same trick that schroedinger.py does and convert it a .vtk which gets interpreted as 2D:
>>
>> aluminum:sfepy steve$ ./script/mesh_to_vtk.py t1.mesh tmp/t1.vtk
>>
>> Now I've got:
>>
>> aluminum:sfepy steve$ cat tmp/t1.vtk
>> # vtk DataFile Version 2.0
>> generated by mesh_to_vtk.py
>> ASCII
>> DATASET UNSTRUCTURED_GRID
>> POINTS 6 float
>> 0 0 0
>> 0.1 0 0
>> 0.1 0.3 0
>> 0 0.3 0
>> 0.1 0.15000000000017 0
>> 0 0.15000000000017 0
>> CELLS 4 16
>> 3 0 1 5
>> 3 1 4 5
>> 3 2 3 5
>> 3 4 2 5
>> CELL_TYPES 4
>> 5
>> 5
>> 5
>> 5
>> CELL_DATA 4
>>
>> SCALARS mat_id int 1
>> LOOKUP_TABLE default
>> 6
>> 6
>> 6
>> 6
>>
>> aluminum:sfepy steve$ ./dump-info.py tmp/t1.vtk
>> sfepy: reading mesh (tmp/t1.vtk)...
>> sfepy: ...done in 0.00 s
>> Mesh:tmp/t1
>> conns:
>> list: [array([[0, 1, 5],
>> [1, 4, 5],
>> [2, 3, 5],
>> [4, 2, 5]])]
>> coors:
>> (6, 2) array of float64
>> descs:
>> list: ['2_3']
>> dim:
>> 2
>> dims:
>> list: [2]
>> el_offsets:
>> (2,) array of int32
>> io:
>> None
>> mat_ids:
>> list: [array([6, 6, 6, 6])]
>> n_e_ps:
>> (1,) array of int32
>> n_el:
>> 4
>> n_els:
>> (1,) array of int32
>> n_nod:
>> 6
>> name:
>> tmp/t1
>> ngroups:
>> (6,) array of int32
>> nodal_bcs:
>> dict with keys: []
>> setup_done:
>> 0
>>
>> Which get's correctly loaded up as 2D, but also no hint of group 1 and group 2.
>>
>> Now... if I run simple.py with regions defined as
>>
>> regions = {
>> 'Omega' : ('all', {}),
>> 'Omega_Source' : ('nodes in (y<0.15)', {}),
>> 'Omega_Target' : ('nodes in (y>0.25)', {}),
>> }
>>
>> everything is fine, and I get the correct result:
>>
>> aluminum:sfepy steve$ python simple.py es-foo.py
>> sfepy: left over: ['newregions', 'verbose', '__builtins__', 'get_r_dependence', '__doc__', '__name__', 'data_dir', 'nm', '__package__', '_filename', '__file__']
>> sfepy: reading mesh (/Users/steve/Development/sfepy/tmp/t1.vtk)...
>> sfepy: ...done in 0.00 s
>> sfepy: creating regions...
>> sfepy: Omega_Target
>> sfepy: Omega
>> sfepy: Omega_Source
>> sfepy: ...done in 0.01 s
>> sfepy: equation "Laplace equation":
>> sfepy: dw_laplace.i1.Omega( m.val, v, u ) = 0
>> sfepy: setting up dof connectivities...
>> sfepy: ...done in 0.00 s
>> sfepy: using solvers:
>> nls: newton
>> ls: ls
>> sfepy: updating variables...
>> sfepy: ...done
>> sfepy: matrix shape: (2, 2)
>> sfepy: assembling matrix graph...
>> sfepy: ...done in 0.00 s
>> sfepy: matrix structural nonzeros: 4 (1.00e+00% fill)
>> sfepy: updating materials...
>> sfepy: m
>> sfepy: ...done in 0.00 s
>> sfepy: nls: iter: 0, residual: 1.600781e+01 (rel: 1.000000e+00)
>> sfepy: rezidual: 0.00 [s]
>> sfepy: solve: 0.00 [s]
>> sfepy: matrix: 0.00 [s]
>> sfepy: nls: iter: 1, residual: 4.076197e-15 (rel: 2.546380e-16)
>>
>> results look good:
>>
>> aluminum:sfepy steve$ cat t1.vtk
>> # vtk DataFile Version 2.0
>> step 0 time 0.000000e+00 normalized time 0.000000e+00, generated by simple.py
>> ASCII
>> DATASET UNSTRUCTURED_GRID
>>
>> POINTS 6 float
>> 0.000000e+00 0.000000e+00 0.000000e+00
>> 1.000000e-01 0.000000e+00 0.000000e+00
>> 1.000000e-01 3.000000e-01 0.000000e+00
>> 0.000000e+00 3.000000e-01 0.000000e+00
>> 1.000000e-01 1.500000e-01 0.000000e+00
>> 0.000000e+00 1.500000e-01 0.000000e+00
>>
>> CELLS 4 16
>> 3 0 1 5
>> 3 1 4 5
>> 3 2 3 5
>> 3 4 2 5
>>
>> CELL_TYPES 4
>> 5
>> 5
>> 5
>> 5
>>
>> POINT_DATA 6
>>
>> SCALARS node_groups int 1
>> LOOKUP_TABLE default
>> 0
>> 0
>> 0
>> 0
>> 0
>> 0
>>
>> SCALARS u float 1
>> LOOKUP_TABLE default
>> 0.000000e+00
>> 0.000000e+00
>> 1.500000e+02
>> 1.500000e+02
>> 1.102273e+02
>> 1.147727e+02
>>
>> CELL_DATA 4
>> SCALARS mat_id int 1
>> LOOKUP_TABLE default
>> 6
>> 6
>> 6
>> 6
>>
>> But if I try the regions defined by group I get:
>>
>> aluminum:sfepy steve$ python simple.py es-foo.py
>> sfepy: left over: ['verbose', '__builtins__', 'get_r_dependence', '__doc__', '__name__', 'data_dir', 'nm', '__package__', 'oldregions', '_filename', '__file__']
>> sfepy: reading mesh (/Users/steve/Development/sfepy/tmp/t1.vtk)...
>> sfepy: ...done in 0.00 s
>> sfepy: creating regions...
>> sfepy: Omega_Target
>> sfepy: Omega
>> sfepy: Omega_Source
>> sfepy: ...done in 0.01 s
>> sfepy: equation "Laplace equation":
>> sfepy: dw_laplace.i1.Omega( m.val, v, u ) = 0
>> sfepy: setting up dof connectivities...
>> sfepy: ...done in 0.00 s
>> sfepy: using solvers:
>> nls: newton
>> ls: ls
>> sfepy: updating variables...
>> sfepy: ...done
>> sfepy: matrix shape: (6, 6)
>> sfepy: assembling matrix graph...
>> sfepy: ...done in 0.00 s
>> sfepy: matrix structural nonzeros: 24 (6.67e-01% fill)
>> sfepy: updating materials...
>> sfepy: m
>> sfepy: ...done in 0.00 s
>> sfepy: nls: iter: 0, residual: 0.000000e+00 (rel: 0.000000e+00)
>>
>> aluminum:sfepy steve$ cat t1.vtk
>> # vtk DataFile Version 2.0
>> step 0 time 0.000000e+00 normalized time 0.000000e+00, generated by simple.py
>> ASCII
>> DATASET UNSTRUCTURED_GRID
>>
>> POINTS 6 float
>> 0.000000e+00 0.000000e+00 0.000000e+00
>> 1.000000e-01 0.000000e+00 0.000000e+00
>> 1.000000e-01 3.000000e-01 0.000000e+00
>> 0.000000e+00 3.000000e-01 0.000000e+00
>> 1.000000e-01 1.500000e-01 0.000000e+00
>> 0.000000e+00 1.500000e-01 0.000000e+00
>>
>> CELLS 4 16
>> 3 0 1 5
>> 3 1 4 5
>> 3 2 3 5
>> 3 4 2 5
>>
>> CELL_TYPES 4
>> 5
>> 5
>> 5
>> 5
>>
>> POINT_DATA 6
>>
>> SCALARS node_groups int 1
>> LOOKUP_TABLE default
>> 0
>> 0
>> 0
>> 0
>> 0
>> 0
>>
>> SCALARS u float 1
>> LOOKUP_TABLE default
>> 0.000000e+00
>> 0.000000e+00
>> 0.000000e+00
>> 0.000000e+00
>> 0.000000e+00
>> 0.000000e+00
>>
>> CELL_DATA 4
>> SCALARS mat_id int 1
>> LOOKUP_TABLE default
>> 6
>> 6
>> 6
>> 6
>>
>> Everything is just zero. ;-(
>>
>> Just for completeness here is my test 'es-foo.py' file:
>>
>> r"""
>> Laplace equation with variable material in 2D
>>
>> """
>> import numpy as nm
>> from sfepy import data_dir
>>
>> #filename_mesh = data_dir + '/es-lens.mesh'
>> #filename_mesh = data_dir + '/tmp/es-lens.vtk'
>> filename_mesh = data_dir + '/tmp/t1.vtk'
>>
>> options = {
>> 'nls' : 'newton',
>> 'ls' : 'ls',
>> }
>>
>> materials = {
>> 'm' : 'get_r_dependence',
>> }
>>
>> oldregions = {
>> 'Omega' : ('all', {}),
>> 'Omega_Source' : ('nodes in (y<0.15)', {}),
>> 'Omega_Target' : ('nodes in (y>0.25)', {}),
>> }
>>
>> regions = {
>> 'Omega' : ('all', {}),
>> 'Omega_Source' : ('nodes of group 1', {}),
>> 'Omega_Target' : ('nodes of group 2', {}),
>> }
>>
>> fields = {
>> 'potential' : ('real', 1, 'Omega', 1),
>> }
>>
>> variables = {
>> 'u' : ('unknown field', 'potential', 0),
>> 'v' : ('test field', 'potential', 'u'),
>> }
>>
>> ebcs = {
>> 'u1' : ('Omega_Source', {'u.0' : 0.0}),
>> 'u2' : ('Omega_Target', {'u.0' : 150.0}),
>> }
>>
>> integrals = {
>> 'i1' : ('v', 'gauss_o1_d2'),
>> }
>>
>> equations = {
>> 'Laplace equation' :
>> """dw_laplace.i1.Omega( m.val, v, u ) = 0 """
>> }
>>
>> solvers = {
>> 'ls' : ('ls.scipy_direct', {}),
>> 'newton' : ('nls.newton',
>> {'i_max' : 1,
>> 'eps_a' : 1e-10,
>> }),
>> }
>>
>> def get_r_dependence(ts, coors, mode=None, **kwargs):
>> """
>> We want to add an r factor to the laplacian integral.
>>
>> For scalar parameters, the shape has to be set to `(coors.shape[0], 1, 1)`.
>> """
>> if mode == 'qp':
>> x = coors[:,1]
>> val = x.copy()
>> val.shape = (coors.shape[0], 1, 1)
>> return {'val' : val}
>>
>> functions = {
>> 'get_r_dependence' : (get_r_dependence,),
>> }
>>
>> Anyway... sorry for the extremely long winded explanation!
>>
>> Thanks for any insight.
>> -steve
>>
>>
>> On Jun 9, 2012, at 2:08 AM, Robert Cimrman wrote:
>>
>>> The ebcs appear to be given correctly. They should be applied automatically, when using simple.py. In your eigenvalue problem solver you should take care of that by calling problem.time_update(), as in schroedinger_app.py, line 104.
>>>
>>> Does it help?
>>>
>>> r.
>>>
>>
>> --
>> You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
>> To post to this group, send email to sfepy...(a)googlegroups.com.
>> To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
>> For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
>>
>>
>> --
>> You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
>> To post to this group, send email to sfepy...(a)googlegroups.com.
>> To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
>> For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
>
> --
> You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
> To post to this group, send email to sfepy...(a)googlegroups.com.
> To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
>
--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
Hi Steve,
I cannotbhelp you much with gmsh - try asking the gmsh authors directly, if it is possible at all to export nodal group ids.
r.
----- Reply message -----
From: "Steve Spicklemire" <st...(a)spvi.com>
To: <sfepy...(a)googlegroups.com>
Cc: "Steve Spicklemire" <st...(a)spvi.com>
Subject: Cylindrical coordinates
Date: Sat, Jun 9, 2012 15:01
Hi Richard,
I don't *think* that's the problem.
What follows is a ridiculously long winded explanation of the fact that I don't know how to specify groups of nodes in a .geo file. I've played around with various possibilities (Physical Line, Physical Point, etc.), but none seems to do the trick. My current plan is to start with a .geo template file that my app will render with particular values depending on user-prescribed variables. Once a .geo file is saved on disk, I can run gmsh/mesh_to_vtk.py and then sfepy. Anyway.. sorry for the verbosity.. but here it goes!
-------------
When I specify the BCs this way (case 1) everything works:
regions = {
'Omega' : ('all', {}),
'Omega_Source' : ('nodes in (y<0.15)', {}),
'Omega_Target' : ('nodes in (y>0.25)', {}),
}
But this is awkward since I'll have some potentially complicated geometries. I'd rather "tag" some surfaces/lines as belonging to a particular group and using the syntax (case 2):
regions = {
'Omega' : ('all', {}),
'Omega_Source' : ('nodes of group 1', {}),
'Omega_Target' : ('nodes of group 2', {}),
}
Where group 1 and group 2 are set up in my original .geo file.
(BTW I think the socket api for gmsh is only for connecting a solver or a postprocessor using gmsh as the UI, not for creating meshes from outside applications)
Let's say I start with a drastically simple .geo file like so:
lc = 0.1;
Point(1) = {0, 0, 0, lc};
Point(2) = {.1, 0, 0, lc} ;
Point(3) = {.1, .3, 0, lc} ;
Point(4) = {0, .3, 0, lc} ;
Line(1) = {1,2} ;
Line(2) = {3,2} ;
Line(3) = {3,4} ;
Line(4) = {4,1} ;
Line Loop(5) = {4,1,-2,3} ;
Plane Surface(6) = {5} ;
Physical Point(1) = {1,2} ;
Physical Point(2) = {3,4} ;
Physical Surface(7) = {6};
(I'm hoping that "Physical Point(1)" might correspond to "group 1" in the region definition (case 2) and "Physical Point(2)" may be "group 2". No luck on that I'm afraid.)
I use gmsh to convert it to a mesh:
aluminum:sfepy steve$ gmsh -2 t1.geo -format mesh
Info : Running 'gmsh -2 t1.geo -format mesh' [1 node(s), max. 1 thread(s)]
Info : Started on Sat Jun 9 06:33:30 2012
Info : Reading 't1.geo'...
Info : Done reading 't1.geo'
Info : Meshing 1D...
Info : Meshing curve 1 (Line)
Info : Meshing curve 2 (Line)
Info : Meshing curve 3 (Line)
Info : Meshing curve 4 (Line)
Info : Done meshing 1D (0.000122 s)
Info : Meshing 2D...
Info : Meshing surface 6 (Plane, MeshAdapt)
Info : Done meshing 2D (0.001485 s)
Info : 6 vertices 14 elements
Info : Writing 't1.mesh'...
Info : Done writing 't1.mesh'
Info : Stopped on Sat Jun 9 06:33:30 2012
Now I've got the following, but no hints of group 1 or group 2 ;-(.
MeshVersionFormatted 1
Dimension
3
Vertices
6
0 0 0 1
0.1 0 0 2
0.1 0.3 0 3
0 0.3 0 4
0.1 0.15000000000017 0 2
0 0.15000000000017 0 4
Triangles
4
1 2 6 6
2 5 6 6
3 4 6 6
5 3 6 6
End
Sadly, sfepy see's this as 3D.
aluminum:sfepy steve$ ./dump-info.py t1.mesh
sfepy: reading mesh (t1.mesh)...
sfepy: ...done in 0.00 s
Mesh:t1
conns:
list: [array([[0, 1, 5],
[1, 4, 5],
[2, 3, 5],
[4, 2, 5]])]
coors:
(6, 3) array of float64
descs:
list: ['2_3']
dim:
3
dims:
list: [2]
el_offsets:
(2,) array of int32
io:
None
mat_ids:
list: [array([6, 6, 6, 6])]
n_e_ps:
(1,) array of int32
n_el:
4
n_els:
(1,) array of int32
n_nod:
6
name:
t1
ngroups:
(6,) array of float64
nodal_bcs:
dict with keys: []
setup_done:
0
(BTW dump_info is just the following script:
aluminum:sfepy steve$ cat dump-info.py
#!/usr/bin/env python
from sfepy.fem import Mesh
import sys
m = Mesh.from_file(sys.argv[1])
print m
)
I can use the same trick that schroedinger.py does and convert it a .vtk which gets interpreted as 2D:
aluminum:sfepy steve$ ./script/mesh_to_vtk.py t1.mesh tmp/t1.vtk
Now I've got:
aluminum:sfepy steve$ cat tmp/t1.vtk
# vtk DataFile Version 2.0
generated by mesh_to_vtk.py
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 6 float
0 0 0
0.1 0 0
0.1 0.3 0
0 0.3 0
0.1 0.15000000000017 0
0 0.15000000000017 0
CELLS 4 16
3 0 1 5
3 1 4 5
3 2 3 5
3 4 2 5
CELL_TYPES 4
5
5
5
5
CELL_DATA 4
SCALARS mat_id int 1
LOOKUP_TABLE default
6
6
6
6
aluminum:sfepy steve$ ./dump-info.py tmp/t1.vtk
sfepy: reading mesh (tmp/t1.vtk)...
sfepy: ...done in 0.00 s
Mesh:tmp/t1
conns:
list: [array([[0, 1, 5],
[1, 4, 5],
[2, 3, 5],
[4, 2, 5]])]
coors:
(6, 2) array of float64
descs:
list: ['2_3']
dim:
2
dims:
list: [2]
el_offsets:
(2,) array of int32
io:
None
mat_ids:
list: [array([6, 6, 6, 6])]
n_e_ps:
(1,) array of int32
n_el:
4
n_els:
(1,) array of int32
n_nod:
6
name:
tmp/t1
ngroups:
(6,) array of int32
nodal_bcs:
dict with keys: []
setup_done:
0
Which get's correctly loaded up as 2D, but also no hint of group 1 and group 2.
Now... if I run simple.py with regions defined as
regions = {
'Omega' : ('all', {}),
'Omega_Source' : ('nodes in (y<0.15)', {}),
'Omega_Target' : ('nodes in (y>0.25)', {}),
}
everything is fine, and I get the correct result:
aluminum:sfepy steve$ python simple.py es-foo.py
sfepy: left over: ['newregions', 'verbose', '__builtins__', 'get_r_dependence', '__doc__', '__name__', 'data_dir', 'nm', '__package__', '_filename', '__file__']
sfepy: reading mesh (/Users/steve/Development/sfepy/tmp/t1.vtk)...
sfepy: ...done in 0.00 s
sfepy: creating regions...
sfepy: Omega_Target
sfepy: Omega
sfepy: Omega_Source
sfepy: ...done in 0.01 s
sfepy: equation "Laplace equation":
sfepy: dw_laplace.i1.Omega( m.val, v, u ) = 0
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: using solvers:
nls: newton
ls: ls
sfepy: updating variables...
sfepy: ...done
sfepy: matrix shape: (2, 2)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 4 (1.00e+00% fill)
sfepy: updating materials...
sfepy: m
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 1.600781e+01 (rel: 1.000000e+00)
sfepy: rezidual: 0.00 [s]
sfepy: solve: 0.00 [s]
sfepy: matrix: 0.00 [s]
sfepy: nls: iter: 1, residual: 4.076197e-15 (rel: 2.546380e-16)
results look good:
aluminum:sfepy steve$ cat t1.vtk
# vtk DataFile Version 2.0
step 0 time 0.000000e+00 normalized time 0.000000e+00, generated by simple.py
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 6 float
0.000000e+00 0.000000e+00 0.000000e+00
1.000000e-01 0.000000e+00 0.000000e+00
1.000000e-01 3.000000e-01 0.000000e+00
0.000000e+00 3.000000e-01 0.000000e+00
1.000000e-01 1.500000e-01 0.000000e+00
0.000000e+00 1.500000e-01 0.000000e+00
CELLS 4 16
3 0 1 5
3 1 4 5
3 2 3 5
3 4 2 5
CELL_TYPES 4
5
5
5
5
POINT_DATA 6
SCALARS node_groups int 1
LOOKUP_TABLE default
0
0
0
0
0
0
SCALARS u float 1
LOOKUP_TABLE default
0.000000e+00
0.000000e+00
1.500000e+02
1.500000e+02
1.102273e+02
1.147727e+02
CELL_DATA 4
SCALARS mat_id int 1
LOOKUP_TABLE default
6
6
6
6
But if I try the regions defined by group I get:
aluminum:sfepy steve$ python simple.py es-foo.py
sfepy: left over: ['verbose', '__builtins__', 'get_r_dependence', '__doc__', '__name__', 'data_dir', 'nm', '__package__', 'oldregions', '_filename', '__file__']
sfepy: reading mesh (/Users/steve/Development/sfepy/tmp/t1.vtk)...
sfepy: ...done in 0.00 s
sfepy: creating regions...
sfepy: Omega_Target
sfepy: Omega
sfepy: Omega_Source
sfepy: ...done in 0.01 s
sfepy: equation "Laplace equation":
sfepy: dw_laplace.i1.Omega( m.val, v, u ) = 0
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: using solvers:
nls: newton
ls: ls
sfepy: updating variables...
sfepy: ...done
sfepy: matrix shape: (6, 6)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 24 (6.67e-01% fill)
sfepy: updating materials...
sfepy: m
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 0.000000e+00 (rel: 0.000000e+00)
aluminum:sfepy steve$ cat t1.vtk
# vtk DataFile Version 2.0
step 0 time 0.000000e+00 normalized time 0.000000e+00, generated by simple.py
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 6 float
0.000000e+00 0.000000e+00 0.000000e+00
1.000000e-01 0.000000e+00 0.000000e+00
1.000000e-01 3.000000e-01 0.000000e+00
0.000000e+00 3.000000e-01 0.000000e+00
1.000000e-01 1.500000e-01 0.000000e+00
0.000000e+00 1.500000e-01 0.000000e+00
CELLS 4 16
3 0 1 5
3 1 4 5
3 2 3 5
3 4 2 5
CELL_TYPES 4
5
5
5
5
POINT_DATA 6
SCALARS node_groups int 1
LOOKUP_TABLE default
0
0
0
0
0
0
SCALARS u float 1
LOOKUP_TABLE default
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
CELL_DATA 4
SCALARS mat_id int 1
LOOKUP_TABLE default
6
6
6
6
Everything is just zero. ;-(
Just for completeness here is my test 'es-foo.py' file:
r"""
Laplace equation with variable material in 2D
"""
import numpy as nm
from sfepy import data_dir
#filename_mesh = data_dir + '/es-lens.mesh'
#filename_mesh = data_dir + '/tmp/es-lens.vtk'
filename_mesh = data_dir + '/tmp/t1.vtk'
options = {
'nls' : 'newton',
'ls' : 'ls',
}
materials = {
'm' : 'get_r_dependence',
}
oldregions = {
'Omega' : ('all', {}),
'Omega_Source' : ('nodes in (y<0.15)', {}),
'Omega_Target' : ('nodes in (y>0.25)', {}),
}
regions = {
'Omega' : ('all', {}),
'Omega_Source' : ('nodes of group 1', {}),
'Omega_Target' : ('nodes of group 2', {}),
}
fields = {
'potential' : ('real', 1, 'Omega', 1),
}
variables = {
'u' : ('unknown field', 'potential', 0),
'v' : ('test field', 'potential', 'u'),
}
ebcs = {
'u1' : ('Omega_Source', {'u.0' : 0.0}),
'u2' : ('Omega_Target', {'u.0' : 150.0}),
}
integrals = {
'i1' : ('v', 'gauss_o1_d2'),
}
equations = {
'Laplace equation' :
"""dw_laplace.i1.Omega( m.val, v, u ) = 0 """
}
solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton',
{'i_max' : 1,
'eps_a' : 1e-10,
}),
}
def get_r_dependence(ts, coors, mode=None, **kwargs):
"""
We want to add an r factor to the laplacian integral.
For scalar parameters, the shape has to be set to `(coors.shape[0], 1, 1)`.
"""
if mode == 'qp':
x = coors[:,1]
val = x.copy()
val.shape = (coors.shape[0], 1, 1)
return {'val' : val}
functions = {
'get_r_dependence' : (get_r_dependence,),
}
Anyway... sorry for the extremely long winded explanation!
Thanks for any insight.
-steve
On Jun 9, 2012, at 2:08 AM, Robert Cimrman wrote:
> The ebcs appear to be given correctly. They should be applied automatically, when using simple.py. In your eigenvalue problem solver you should take care of that by calling problem.time_update(), as in schroedinger_app.py, line 104.
>
> Does it help?
>
> r.
>
--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
The ebcs appear to be given correctly. They should be applied automatically, when using simple.py. In your eigenvalue problem solver you should take care of that by calling problem.time_update(), as in schroedinger_app.py, line 104.
Does it help?
r.
----- Reply message -----
From: "steve" <st...(a)spvi.com>
To: <sfepy...(a)googlegroups.com>
Subject: Cylindrical coordinates
Date: Fri, Jun 8, 2012 14:54
I've simplified this to a case the works a little better (no crash) but still not there.
I've started with the gmsh tutorial file 1: t1.geo and modified it a bit
------------------------------------------lc = 0.009;
Point(1) = {0, 0, 0, lc};
Point(2) = {.1, 0, 0, lc} ;
Point(3) = {.1, .3, 0, lc} ;
Point(4) = {0, .3, 0, lc} ;
Line(1) = {1,2} ;
Line(2) = {3,2} ;
Line(3) = {3,4} ;
Line(4) = {4,1} ;
Line Loop(5) = {4,1,-2,3} ;
Plane Surface(6) = {5} ;
Physical Point(1) = {1,2} ;
Physical Point(2) = {3,4} ;
Physical Surface(7) = {6};-------------------------------
I'd like to set the voltage along line 1 to 0V, and along line 3 to 150V say. I tried 'Physical Line', but that added '2_2' elements to the vtk file that sfepy didn't like. So.. I tried 'Physical Point' and that seemed to allow my code to run, but its still not setting up the region correctly. I used:
regions = {
'Omega' : ('all', {}),
'Omega_Source' : ('nodes of group 1', {}),
'Omega_Target' : ('nodes of group 2', {}),
}
which ran "OK", but didn't set the voltages with:
ebcs = {
'u1' : ('Omega_Source', {'u.0' : 0.0}),
'u2' : ('Omega_Target', {'u.0' : 150.0}),
}
I'll keep looking. Thanks for any insight.
-steve
--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To view this discussion on the web visit https://groups.google.com/d/msg/sfepy-devel/-/XDGjDglwIKwJ.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
The code checks that certain keywords are not missing. Try using ebcs = {} in the input file.
r.
----- Reply message -----
From: "steve" <st...(a)spvi.com>
To: <sfepy...(a)googlegroups.com>
Subject: Cylindrical coordinates
Date: Wed, Jun 6, 2012 23:33
Velocity potential is pretty well defined here:
<http://en.wikipedia.org/wiki/Potential_flow#Description_and_characteristics>
When I try to run with no boundary conditions I'm getting the error below. Does it lead to a singular matrix or something with no essential boundary conditions?
thanks!-steve
Traceback (most recent call last): File "acoustics.py", line 113, in <module> main() File "acoustics.py", line 104, in main conf = ProblemConf.from_file_and_options(filename_in, options, required, other) File "/Users/steve/Development/sfepy/sfepy/base/conf.py", line 311, in from_file_and_options override=override) File "/Users/steve/Development/sfepy/sfepy/base/conf.py", line 290, in from_file required, other, verbose, override) File "/Users/steve/Development/sfepy/sfepy/base/conf.py", line 359, in __init__ required=required, other=other) File "/Users/steve/Development/sfepy/sfepy/base/conf.py", line 369, in setup other_missing = self.validate( required = required, other = other ) File "/Users/steve/Development/sfepy/sfepy/base/conf.py", line 421, in validate raise ValueError('required missing: %s' % required_missing)ValueError: required missing: ['ebc_[0-9]+|ebcs']
On Wednesday, June 6, 2012 3:10:58 PM UTC-6, Robert Cimrman wrote:On 06/06/2012 09:40 PM, steve wrote:
> OK.. on to the eigenvalue problem.
>
> I started with the quantum program... and stripped out all the quantum
> stuff. ;-)
>
> Basically the acoustics problem is exactly the square well potential (V=0
> everywhere) BUT with different boundary conditions.
>
> In quantum_common.py we've got:
>
> ebc_1 = {
> 'name' : 'ZeroSurface',
> 'region' : 'Surface',
> 'dofs' : {'Psi.0' : 0.0},
> }
>
> Which forces psi to be zero on the surface of Omega. Basically Dirichlet
> conditions. For acoustics the velocity potential needs to have no normal
> gradient at the surfaces.. more or less Neumann conditions. Is there an
> easy way to implement that?
What is the exact definition of the velocity potential? If it's really a
Neumann-like boundary integral, having it zero = omitting the term in
equations + no Dirichlet boundary conditions (ebcs) at all.
r.
--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To view this discussion on the web visit https://groups.google.com/d/msg/sfepy-devel/-/e9Qjtk3ai88J.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
Hi Folks,
I'm extremely new to finite elements in general and sfepy in
particular, but I have a couple of problems I'm interested in, and I'm
hoping that sfepy may provide an easy way to get started. Both of
these problems (1 acoustic normal modes in cavities with cylindrical
symmetry and 2 electrostatic fields in cavities bounded by various
equipotential surfaces) have in common a axisymmetric or cylindrical
symmetry. I found this discussion
<http://groups.google.com/group/sfepy-devel/msg/17364ae1fcc57259>
but I didn't see any real follow up. Has there been more done on this?
also.. I've spent the last couple of days skimming Hughes finite
element analysis book to get my head screwed on straight WRT FEM and I
found this interesting tidbit:
<http://books.google.com/books?
id=cHH2n_qBK0IC&printsec=frontcover&dq=finite%20element%20method
%20hughes&hl=en&sa=X&ei=ywHMT-
bKBcfK2AWowNzZCw&output=reader&pg=GBS.PA101>
Since heat conduction is extremely similar to electrostatics (and only
a term or so different from acoustics) I hoped that it might be
possible to get away with a 2D mesh for my problem(s) rather than
running a full 3D mesh of the highly symmetric cases. So.. I guess my
question is, what's the easiest way to implement this feature in
sfepy? Can anyone suggest a path to follow?
thanks!
-steve
Hi Folks,
Trying to get familiar with finding eigenvalues etc and noticed the
schroedinger.py app does almost exactly what I need, except in the quantum
domain. Anyway when I tried it run it I get:
Traceback (most recent call last):
File "schroedinger.py", line 215, in <module>
main()
File "schroedinger.py", line 212, in main
app()
File
"/Users/steve/Downloads/sfepy-2012.2-epd/sfepy/applications/application.py",
line 29, in call_basic
return self.call(**kwargs)
File
"/Users/steve/Downloads/sfepy-2012.2-epd/sfepy/physics/schroedinger_app.py",
line 86, in call
evp = self.solve_eigen_problem()
File
"/Users/steve/Downloads/sfepy-2012.2-epd/sfepy/physics/schroedinger_app.py",
line 122, in solve_eigen_problem
eigs, mtx_s_phi = eig(mtx_a, mtx_b, n_eigs)
File "/Users/steve/Downloads/sfepy-2012.2-epd/sfepy/solvers/eigen.py",
line 43, in _standard_call
**kwargs)
File "/Users/steve/Downloads/sfepy-2012.2-epd/sfepy/solvers/eigen.py",
line 283, in __call__
from pysparse import jdsym, itsolvers, precon
ImportError: cannot import name jdsym
I see that stuff has moved around in pysparse. Snooping a bit I found that
this seems to work:
aluminum:sfepy-2012.2-epd steve$ diff
/Users/steve/Downloads/sfepy-2012.2-epd/sfepy/solvers/eigen.py
../sfepy-2012.2-py2.7/sfepy/solvers/eigen.py
283c283,284
< from pysparse import jdsym, itsolvers, precon
---
> from pysparse.eigen import jdsym
> from pysparse import itsolvers, precon
Looks like I installed pysparse 1.2-dev224. What version do you guys
recommend?
thanks!
-steve
