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:

-- 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.