12:41 a.m.
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:
--
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:
--
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.
--
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.
8:47 p.m.
New subject: Cylindrical coordinates
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
9:04 p.m.
New subject: Cylindrical coordinates
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
9:06 p.m.
New subject: Cylindrical coordinates
Awesome. Thanks!
-steve
On Jun 11, 2012, at 9:34 AM, Robert Cimrman <cimr...(a)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...(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.
9:04 p.m.
New subject: Cylindrical coordinates
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?
2533
days inactive
2534
days old
4 comments
3 participants
participants (3)
-
Robert Cimrman -
steve -
Steve Spicklemire