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" <s...@spvi.com>
To: <sfep...@googlegroups.com>
Cc: "Steve Spicklemire" <s...@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...@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.