Artifacts in solution after switch from tet to hex mesh?
Hello,
Thank you to everyone for all the effort that has been put into SfePy, it is a pleasure to use.
I am trying to solve a transient diffusion problem similar to the poisson_periodic_boundary_condition.py example but on a mesh of hexahedral elements. The example works well with the provided tet mesh. When I switch to a hex mesh the solution has some problems. The solution on gridpoints just outside the source region goes below zero at early times and the solution has some waves along line x=y=z.
It seem like the elements may be read in wrong but I cannot figure out the problem. I create the hex mesh, test_hex.mesh, with this Python code:
import numpy as np from sfepy.mesh.mesh_generators import gen_block_mesh from sfepy.discrete.fem.meshio import MeditMeshIO
mesh = gen_block_mesh(np.array((10,10,1)), np.array((11,11,2)), np.array((5.0,5.0,0.5)), mat_id=9, verbose=True) m = MeditMeshIO("dummy_name") m.write("test_hex.mesh", mesh)
The attached code is a simplification of the poisson_periodic_boundary_condition.py example, I removed the periodic boundary condition and the time-dependent material properties. with_tet_mesh.py gives good results and with_hex_mesh.py shows the problem.
I am running Python 2.7.13 and SfePy version 2017.2; I get the same result from the current SfePy github master. I run the code with the command python simple.py with_hex_mesh.py.
Can someone help me understand where I am going wrong?
Thanks, Jason
Hi Jason,
On 07/11/2017 10:07 PM, Jason Furtney wrote:
Hello,
Thank you to everyone for all the effort that has been put into SfePy, it is a pleasure to use.
Thanks!
I am trying to solve a transient diffusion problem similar to the poisson_periodic_boundary_condition.py example but on a mesh of hexahedral elements. The example works well with the provided tet mesh. When I switch to a hex mesh the solution has some problems. The solution on gridpoints just outside the source region goes below zero at early times and the solution has some waves along line x=y=z.
It seem like the elements may be read in wrong but I cannot figure out the problem. I create the hex mesh, test_hex.mesh, with this Python code:
import numpy as np from sfepy.mesh.mesh_generators import gen_block_mesh from sfepy.discrete.fem.meshio import MeditMeshIO
mesh = gen_block_mesh(np.array((10,10,1)), np.array((11,11,2)), np.array((5.0,5.0,0.5)), mat_id=9, verbose=True) m = MeditMeshIO("dummy_name") m.write("test_hex.mesh", mesh)
The attached code is a simplification of the poisson_periodic_boundary_condition.py example, I removed the periodic boundary condition and the time-dependent material properties. with_tet_mesh.py gives good results and with_hex_mesh.py shows the problem.
I am running Python 2.7.13 and SfePy version 2017.2; I get the same result from the current SfePy github master. I run the code with the command python simple.py with_hex_mesh.py.
Can someone help me understand where I am going wrong?
With the tet mesh, the approximation is linear. With the hex mesh, the functions are tri-linear - even cubic terms appear in the basis. So the problem with the hex mesh is under-integrated, and that causes the artifacts. The solution is easy: increase the order of 'i' to at least 2.
(The reasoning is a bit tricky: the mass matrix term contains a product of two tri-linear functions = order 2 * 3, the laplace term contains a product of two gradients of trilinear functions = order 2 * 2, but the integral order in the hex case means the order in one axis... so 'i' : 1 should be ok for a single tri-linear function, but not for a product of two such functions - there, 2 is needed. Note also that having the mass matrix term means that 'i' : 1 under-integrates even with the tet mesh...).
-> tldr; If you see something weird, try increasing the numerical integration order :)
r.
Dear Robert,
Thanks for the quick reply and for the explanation. Yes, increasing the integration order does the trick.
Thanks, Jason
On Tue, Jul 11, 2017 at 3:38 PM, Robert Cimrman <cimrman3@ntc.zcu.cz> wrote:
Hi Jason,
On 07/11/2017 10:07 PM, Jason Furtney wrote:
Hello,
Thank you to everyone for all the effort that has been put into SfePy, it is a pleasure to use.
Thanks!
I am trying to solve a transient diffusion problem similar to the poisson_periodic_boundary_condition.py example but on a mesh of hexahedral elements. The example works well with the provided tet mesh. When I switch to a hex mesh the solution has some problems. The solution on gridpoints just outside the source region goes below zero at early times and the solution has some waves along line x=y=z.
It seem like the elements may be read in wrong but I cannot figure out the problem. I create the hex mesh, test_hex.mesh, with this Python code:
import numpy as np from sfepy.mesh.mesh_generators import gen_block_mesh from sfepy.discrete.fem.meshio import MeditMeshIO
mesh = gen_block_mesh(np.array((10,10,1)), np.array((11,11,2)), np.array((5.0,5.0,0.5)), mat_id=9, verbose=True) m = MeditMeshIO("dummy_name") m.write("test_hex.mesh", mesh)
The attached code is a simplification of the poisson_periodic_boundary_condition.py example, I removed the periodic boundary condition and the time-dependent material properties. with_tet_mesh.py gives good results and with_hex_mesh.py shows the problem.
I am running Python 2.7.13 and SfePy version 2017.2; I get the same result from the current SfePy github master. I run the code with the command python simple.py with_hex_mesh.py.
Can someone help me understand where I am going wrong?
With the tet mesh, the approximation is linear. With the hex mesh, the functions are tri-linear - even cubic terms appear in the basis. So the problem with the hex mesh is under-integrated, and that causes the artifacts. The solution is easy: increase the order of 'i' to at least 2.
(The reasoning is a bit tricky: the mass matrix term contains a product of two tri-linear functions = order 2 * 3, the laplace term contains a product of two gradients of trilinear functions = order 2 * 2, but the integral order in the hex case means the order in one axis... so 'i' : 1 should be ok for a single tri-linear function, but not for a product of two such functions - there, 2 is needed. Note also that having the mass matrix term means that 'i' : 1 under-integrates even with the tet mesh...).
-> tldr; If you see something weird, try increasing the numerical integration order :)
r.
SfePy mailing list sfepy@python.org https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/
participants (3)
-
Jason Furtney -
Jason Furtney -
Robert Cimrman