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.