SfePy June 2012

sfepy@python.org

2 participants
9 discussions

Re: Torque
by Robert Cimrman 09 Mar '18

09 Mar '18

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.

4 7

1D line elements in SfePy
by Nimish 22 Jan '16

22 Jan '16

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

4 9

Modeling an adjustable X-ray optic
by Tom Aldcroft 21 May '13

21 May '13

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

3 6

Re: Cylindrical coordinates
by Robert Cimrman 11 Jun '12

11 Jun '12

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.

3 4

Re: Cylindrical coordinates
by Robert Cimrman 10 Jun '12

10 Jun '12

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.

2 2

Re: Cylindrical coordinates
by Robert Cimrman 09 Jun '12

09 Jun '12

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.

2 1

Re: Cylindrical coordinates
by Robert Cimrman 08 Jun '12

08 Jun '12

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.

2 5

Cylindrical coordinates
by steve 06 Jun '12

06 Jun '12

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

3 16

jdsym in motion...
by steve 05 Jun '12

05 Jun '12

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

3 4

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

SfePy June 2012