Re: Cylindrical coordinates
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...@spvi.com> To: <sfepy...@googlegroups.com> Cc: "Steve Spicklemire" <st...@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...@spvi.com> To: <sfepy...@googlegroups.com> Cc: "Steve Spicklemire" <st...@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.
-- 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.
-- 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.
-- 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.
Hi Robert,
So.. I broke down and "template"ized a set of boundary condition functions in a separate file that can be imported into the problem description. It works pretty well I guess. My app can render the template and run the sfepy process to get the fields. First result attached. Noticed the "needle point" boundary which was the one I was hoping to handle with group ids.. but it's pretty easy to cook up a python function to find the nodes.
Now. I've got a vtk file full of lovely field data. In my relaxation version I interpolate over my rectangular grid to find the electric field and compute particle trajectories. So... I've found the vtk module and can read in data from the output file. I'm playing with various ways to interpolate and estimate the field at various positions in the domain, but I'm guessing this problem has been solved before. Is it best to re-map the data onto a rectangular grid, or is there some way to quickly find the nodes nearest some particular spatial coordinates and interpolate directly from the irregular grid data?
thanks! -steve
On 06/11/2012 05:17 PM, steve wrote:
Hi Robert,
So.. I broke down and "template"ized a set of boundary condition functions in a separate file that can be imported into the problem description. It works pretty well I guess. My app can render the template and run the sfepy process to get the fields. First result attached. Noticed the "needle point" boundary which was the one I was hoping to handle with group ids.. but it's pretty easy to cook up a python function to find the nodes.
Great!
Now. I've got a vtk file full of lovely field data. In my relaxation version I interpolate over my rectangular grid to find the electric field and compute particle trajectories. So... I've found the vtk module and can read in data from the output file. I'm playing with various ways to interpolate and estimate the field at various positions in the domain, but I'm guessing this problem has been solved before. Is it best to re-map the data onto a rectangular grid, or is there some way to quickly find the nodes nearest some particular spatial coordinates and interpolate directly from the irregular grid data?
Yes, check out [1]. I would use PointsProbe to get the data in a grid, instead of LineProbe mentioned in the Primer text. The text show how to make matplotlib plots of data over given lines. If you need more control, you will need to manually create the probe as in gen_lines(), and then use the code from get_it() in probe_hook() to sample the data. For this a ProblemDefinition instance is needed (called "problem" there), but I guess you have it in your code already. Anyway, ask in case of problems.
Cheers, r.
[1] http://sfepy.org/doc-devel/primer.html#probing
thanks! -steve
Awesome. Thanks!
-steve
On Jun 11, 2012, at 9:34 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
On 06/11/2012 05:17 PM, steve wrote:
Hi Robert,
So.. I broke down and "template"ized a set of boundary condition functions in a separate file that can be imported into the problem description. It works pretty well I guess. My app can render the template and run the sfepy process to get the fields. First result attached. Noticed the "needle point" boundary which was the one I was hoping to handle with group ids.. but it's pretty easy to cook up a python function to find the nodes.
Great!
Now. I've got a vtk file full of lovely field data. In my relaxation version I interpolate over my rectangular grid to find the electric field and compute particle trajectories. So... I've found the vtk module and can read in data from the output file. I'm playing with various ways to interpolate and estimate the field at various positions in the domain, but I'm guessing this problem has been solved before. Is it best to re-map the data onto a rectangular grid, or is there some way to quickly find the nodes nearest some particular spatial coordinates and interpolate directly from the irregular grid data?
Yes, check out [1]. I would use PointsProbe to get the data in a grid, instead of LineProbe mentioned in the Primer text. The text show how to make matplotlib plots of data over given lines. If you need more control, you will need to manually create the probe as in gen_lines(), and then use the code from get_it() in probe_hook() to sample the data. For this a ProblemDefinition instance is needed (called "problem" there), but I guess you have it in your code already. Anyway, ask in case of problems.
Cheers, r.
[1] http://sfepy.org/doc-devel/primer.html#probing
thanks! -steve
-- 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.
So.. I've discovered that I can do something like
z = array([x,y])-m.coors zz = (z*z).sum(axis=1) np.where(zz<0.1)
to get the index of the point closest to the field point (with sorting or something if I get more than one.)
Next I can figure out from m.conns which nodes are connected to that one, and go on from there?
is this a reasonable path?
thanks, -steve
Is it best to re-map the data onto a rectangular grid, or is there some way
to quickly find the nodes nearest some particular spatial coordinates and interpolate directly from the irregular grid data?
participants (3)
-
Robert Cimrman -
steve -
Steve Spicklemire