The latest commits in the master branch on github incorporate meshio
package [1] into SfePy. We use meshio to read the following file
formats: abaqus, exodus, gmsh, medit, nastran, vtk, vtu, med, xdmf and
tetgen. Besides these formats, we keep on maintaining own code for
reading: ansys_cdb, xyz, comsol, hmascii, gambit and mesh3d.
Bug fixes and improvements are welcome!
Regards,
Vladimír Lukeš
[1] https://github.com/nschloe/meshio
I want to define my own field. For example, pressure.
```python
fields = {
'temperature' : ('real', 1, 'Omega', 1),
'pressure' : ('real', 1, 'Omega', 1),
}
variables = {
't' : ('unknown field', 'temperature', 0),
's' : ('test field', 'temperature', 't'),
'p' : ( 'parameter field', 'pressure', {'setter' : 'get_pressure'} ),
}
def get_pressure(ts, coors, region=None):
z = coors[:,2]
return z
functions = {
'get_pressure': (get_pressure,),
}
```
The idea is to create very primitive field that depends on z coordinates.
Now I want to export the grad(p) to the output mesh. And so I add this in
postprocess function:
```python
p_grad = pb.evaluate('ev_grad.i.Omega(p)', mode='el_avg', verbose=False)
out['p'] = Struct(name='output_data', mode='cell', data=p_grad, dofs=None)
```
But I get the following error:
ValueError: argument p not found!
What I'm doing wrong?
And how can I export the pressure field itself (not its gradient)?
The answer from Robert Cimrman:
Hi,
In a post-processing function, the only variables that are defined are
those
that are used in the equations. If p is not used, you need to create it, as
shown below (also how to save p and its cell averages.
def post_process(out, pb, state, extend=False):
from sfepy.base.base import Struct
pvar = pb.create_variables(['p'])['p']
pvar.time_update(None, pb.functions)
out.update(pvar.create_output())
ap = pb.evaluate('ev_volume_integrate.i.Omega(p)', mode='el_avg',
p=pvar,
verbose=False)
out['ap'] = Struct(name='output_data', mode='cell', data=ap, dofs=None)
p_grad = pb.evaluate('ev_grad.i.Omega(p)', mode='el_avg', p=pvar,
verbose=False)
out['gp'] = Struct(name='output_data', mode='cell', data=p_grad,
dofs=None)
return out
r.
PS: This might be interesting to others, so, if you do not mind, you can
forward it to the list.
Yes, mmesh.point_data['mat_ids'] would work, but I did not see it in your mesh :)
r.
On 2/3/20 4:43 AM, FaTune wrote:
> Your little snipped has worked like magic! Thank you so much!
>
> Am I right that to also save 'mat_ids' I have to change the argument
> [np.zeros_like(conn[:, 0])] to mmesh.point_data['mat_ids']?
>
> ```python
> smesh = Mesh.from_data(
> 'meshio-mesh',
> mmesh.points,
> mmesh.point_data['node_groups'],
> [conn],
> mmesh.point_data['mat_ids'],
> ['3_4'],
> )
> ```
>
> On Mon, Feb 3, 2020 at 3:45 AM Robert Cimrman <cimrman3(a)ntc.zcu.cz> wrote:
>
>> I was thinking about meshio too. But because it does not write the old VTK
>> format, you could try the following:
>>
>> ```
>> import numpy as np
>> import meshio
>>
>> from sfepy.discrete.fem import Mesh
>>
>> mmesh = meshio.read('trian.vtk')
>>
>> conn = mmesh.cells['tetra']
>> smesh = Mesh.from_data(
>> 'meshio-mesh',
>> mmesh.points,
>> mmesh.point_data['node_groups'],
>> [conn],
>> [np.zeros_like(conn[:, 0])],
>> ['3_4'],
>> )
>>
>> smesh.write('trian2.vtk')
>> ```
>>
>> Then
>>
>> postproc.py --only-names=node_groups -b trian2.vtk --wire
>>
>> seems to show the regions correctly.
>>
>> r.
>>
>> On 2/2/20 1:14 PM, FaTune wrote:
>>> Thank you for help. I'm not sure I know how to convert the mesh into
>> older
>>> format. Is there any converter? I've tried to convert trian.vtk to
>>> trian.mesh using meshio-convert in ubuntu. But the resulting trian.mesh
>>> doesn't have any information about node_groups.
>>>
>>> On Sun, Feb 2, 2020 at 9:05 PM Robert Cimrman <cimrman3(a)ntc.zcu.cz>
>> wrote:
>>>
>>>> Can you try converting the mesh to "# vtk DataFile Version 2.0" format?
>>>> Your
>>>> mesh is in "# vtk DataFile Version 4.2", and the reader in sfepy does
>> not
>>>> read
>>>> it correctly.
>>>>
>>>> r.
>>>>
>>>> On 2/2/20 1:49 AM, FaTune wrote:
>>>>> I see the node_groups in paraview, but don't see it if I use your
>>>> snipped:
>>>>>
>>>>> python postproc.py trian.vtk --all -b
>>>>>
>>>>> The file trian.vtk is in the attachment. Clearly something is wrong in
>>>> this
>>>>> mesh.
>>>>>
>>>>> On Sun, Feb 2, 2020 at 10:49 AM Robert Cimrman <cimrman3(a)ntc.zcu.cz>
>>>> wrote:
>>>>>
>>>>>> On 2/2/20 12:28 AM, FaTune wrote:
>>>>>>> In addition to the previous message. I created vtk mesh with values
>> in
>>>>>>> field "node_groups" to be 1 or 2. When I use 'vertices of group 0'
>>>> sfepy
>>>>>>> works without errors, but set this region to every single vertex in
>> the
>>>>>>> mesh, and I don't have values 0 in 'node_groups' field. When I
>> change 0
>>>>>> to
>>>>>>> 1 or 2 it raises error. Apparently sfepy just ignores 'node_groups'
>>>>>> field.
>>>>>>
>>>>>> It is possible to display the vertex and cell groups in a mesh with
>> the
>>>>>> following (or use any VTK viewer, such as paraview):
>>>>>>
>>>>>> python postproc.py meshes/3d/acoustic_wg.vtk --all -b
>>>>>>
>>>>>> Does it show your regions correctly? If the mesh and the node_groups
>> are
>>>>>> well
>>>>>> defined, what you do should work.
>>>>>>
>>>>>> If the above does not help you finding what's wrong, send here a
>> minimal
>>>>>> example that demonstrates the problem.
>>>>>>
>>>>>> r.
>>>>>>
>>>>>>> On Sun, Feb 2, 2020 at 10:16 AM <fatune(a)gmail.com> wrote:
>>>>>>>
>>>>>>>> regions = {
>>>>>>>> 'Omega' : 'all',
>>>>>>>> 'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'),
>>>>>>>> 'Gamma_Right' : ('vertices in (x > 99.099999)', 'facet'),
>>>>>>>> 'GammaIn': ('vertices of group 0', 'edge'),
>>>>>>>> }
>>>>>>>>
>>>>>>>>
>>>>>>>> The group I am struggling with is GammaIn. I was using
>>>>>> vibro_acoustic3d.py
>>>>>>>> as an example when I was creating my vtk mesh. Am I right in my
>>>>>> assumption
>>>>>>>> that it is enough to add field "node_groups" to my vtk and sfepy
>> will
>>>>>>>> gather information of the group number from this field?
I was thinking about meshio too. But because it does not write the old VTK
format, you could try the following:
```
import numpy as np
import meshio
from sfepy.discrete.fem import Mesh
mmesh = meshio.read('trian.vtk')
conn = mmesh.cells['tetra']
smesh = Mesh.from_data(
'meshio-mesh',
mmesh.points,
mmesh.point_data['node_groups'],
[conn],
[np.zeros_like(conn[:, 0])],
['3_4'],
)
smesh.write('trian2.vtk')
```
Then
postproc.py --only-names=node_groups -b trian2.vtk --wire
seems to show the regions correctly.
r.
On 2/2/20 1:14 PM, FaTune wrote:
> Thank you for help. I'm not sure I know how to convert the mesh into older
> format. Is there any converter? I've tried to convert trian.vtk to
> trian.mesh using meshio-convert in ubuntu. But the resulting trian.mesh
> doesn't have any information about node_groups.
>
> On Sun, Feb 2, 2020 at 9:05 PM Robert Cimrman <cimrman3(a)ntc.zcu.cz> wrote:
>
>> Can you try converting the mesh to "# vtk DataFile Version 2.0" format?
>> Your
>> mesh is in "# vtk DataFile Version 4.2", and the reader in sfepy does not
>> read
>> it correctly.
>>
>> r.
>>
>> On 2/2/20 1:49 AM, FaTune wrote:
>>> I see the node_groups in paraview, but don't see it if I use your
>> snipped:
>>>
>>> python postproc.py trian.vtk --all -b
>>>
>>> The file trian.vtk is in the attachment. Clearly something is wrong in
>> this
>>> mesh.
>>>
>>> On Sun, Feb 2, 2020 at 10:49 AM Robert Cimrman <cimrman3(a)ntc.zcu.cz>
>> wrote:
>>>
>>>> On 2/2/20 12:28 AM, FaTune wrote:
>>>>> In addition to the previous message. I created vtk mesh with values in
>>>>> field "node_groups" to be 1 or 2. When I use 'vertices of group 0'
>> sfepy
>>>>> works without errors, but set this region to every single vertex in the
>>>>> mesh, and I don't have values 0 in 'node_groups' field. When I change 0
>>>> to
>>>>> 1 or 2 it raises error. Apparently sfepy just ignores 'node_groups'
>>>> field.
>>>>
>>>> It is possible to display the vertex and cell groups in a mesh with the
>>>> following (or use any VTK viewer, such as paraview):
>>>>
>>>> python postproc.py meshes/3d/acoustic_wg.vtk --all -b
>>>>
>>>> Does it show your regions correctly? If the mesh and the node_groups are
>>>> well
>>>> defined, what you do should work.
>>>>
>>>> If the above does not help you finding what's wrong, send here a minimal
>>>> example that demonstrates the problem.
>>>>
>>>> r.
>>>>
>>>>> On Sun, Feb 2, 2020 at 10:16 AM <fatune(a)gmail.com> wrote:
>>>>>
>>>>>> regions = {
>>>>>> 'Omega' : 'all',
>>>>>> 'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'),
>>>>>> 'Gamma_Right' : ('vertices in (x > 99.099999)', 'facet'),
>>>>>> 'GammaIn': ('vertices of group 0', 'edge'),
>>>>>> }
>>>>>>
>>>>>>
>>>>>> The group I am struggling with is GammaIn. I was using
>>>> vibro_acoustic3d.py
>>>>>> as an example when I was creating my vtk mesh. Am I right in my
>>>> assumption
>>>>>> that it is enough to add field "node_groups" to my vtk and sfepy will
>>>>>> gather information of the group number from this field?
>>>>>> _______________________________________________
>>>> _______________________________________________
>>>> SfePy mailing list -- sfepy(a)python.org
>>>> To unsubscribe send an email to sfepy-leave(a)python.org
>>>> https://mail.python.org/mailman3/lists/sfepy.python.org/
>>>>
>>>
>> _______________________________________________
>> SfePy mailing list -- sfepy(a)python.org
>> To unsubscribe send an email to sfepy-leave(a)python.org
>> https://mail.python.org/mailman3/lists/sfepy.python.org/
>>
>
Can you try converting the mesh to "# vtk DataFile Version 2.0" format? Your
mesh is in "# vtk DataFile Version 4.2", and the reader in sfepy does not read
it correctly.
r.
On 2/2/20 1:49 AM, FaTune wrote:
> I see the node_groups in paraview, but don't see it if I use your snipped:
>
> python postproc.py trian.vtk --all -b
>
> The file trian.vtk is in the attachment. Clearly something is wrong in this
> mesh.
>
> On Sun, Feb 2, 2020 at 10:49 AM Robert Cimrman <cimrman3(a)ntc.zcu.cz> wrote:
>
>> On 2/2/20 12:28 AM, FaTune wrote:
>>> In addition to the previous message. I created vtk mesh with values in
>>> field "node_groups" to be 1 or 2. When I use 'vertices of group 0' sfepy
>>> works without errors, but set this region to every single vertex in the
>>> mesh, and I don't have values 0 in 'node_groups' field. When I change 0
>> to
>>> 1 or 2 it raises error. Apparently sfepy just ignores 'node_groups'
>> field.
>>
>> It is possible to display the vertex and cell groups in a mesh with the
>> following (or use any VTK viewer, such as paraview):
>>
>> python postproc.py meshes/3d/acoustic_wg.vtk --all -b
>>
>> Does it show your regions correctly? If the mesh and the node_groups are
>> well
>> defined, what you do should work.
>>
>> If the above does not help you finding what's wrong, send here a minimal
>> example that demonstrates the problem.
>>
>> r.
>>
>>> On Sun, Feb 2, 2020 at 10:16 AM <fatune(a)gmail.com> wrote:
>>>
>>>> regions = {
>>>> 'Omega' : 'all',
>>>> 'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'),
>>>> 'Gamma_Right' : ('vertices in (x > 99.099999)', 'facet'),
>>>> 'GammaIn': ('vertices of group 0', 'edge'),
>>>> }
>>>>
>>>>
>>>> The group I am struggling with is GammaIn. I was using
>> vibro_acoustic3d.py
>>>> as an example when I was creating my vtk mesh. Am I right in my
>> assumption
>>>> that it is enough to add field "node_groups" to my vtk and sfepy will
>>>> gather information of the group number from this field?
>>>> _______________________________________________
>> _______________________________________________
>> SfePy mailing list -- sfepy(a)python.org
>> To unsubscribe send an email to sfepy-leave(a)python.org
>> https://mail.python.org/mailman3/lists/sfepy.python.org/
>>
>
Hello. Is there a way to somehow mark some vertices in vtk mesh so I could later use them in region definition? I've tried to create my vtk mesh with vertices field "node_groups". And some vertices have this field to set to 1, and other to 2. But it doesnt work. Sfepy tell me
>raise ValueError('region "%s" has no entities!' % self.name)
Thank you
