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...@spvi.com> To: <sfepy...@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...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
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.
participants (2)
-
Robert Cimrman -
Steve Spicklemire