Point data in high order meshes #706
-
How would one go about saving from skfem.io import from_meshio
from pygmsh import opencascade, generate_mesh as gm
import numpy as np
from skfem import *
def generateMesh():
R_out = 1.0e-2
geo = opencascade.Geometry
geom = geo(R_out / 30.0, R_out / 10.)
geom.add_raw_code("Mesh.ElementOrder=2;")
outer_ball = geom.add_ball([0.0, 0.0, 0.0], R_out)
msh = gm(geom)
return msh
mesh = from_meshio(generateMesh())
uelem = ElementVectorH1(ElementTetCCR())
elems = {
"u": uelem
}
basis = {
field: InteriorBasis(mesh, e, intorder=4)
for field, e in elems.items()
}
du = basis["u"].zeros()
mesh.save("filename.xdmf", point_data = {"u": du[basis["u"].nodal_dofs].T}) # would work if linear mesh raises Traceback (most recent call last):
File "D:\Research\runHyperelasticSphere.py", line 433, in <module>
mesh.save("filename.xdmf", point_data = {"u": du[basis["u"].nodal_dofs].T}) # would work if linear mesh
File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\mesh\mesh.py", line 424, in save
return to_file(self, filename, point_data, **kwargs)
File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\io\meshio.py", line 154, in to_file
meshio.write(filename, to_meshio(mesh, point_data), **kwargs)
File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\io\meshio.py", line 150, in to_meshio
return meshio.Mesh(mesh.p.T, cells, point_data)
File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\meshio\_mesh.py", line 61, in __init__
raise ValueError(
ValueError: len(points) = 30317, but len(point_data["u"]) = 4123
How would one go about incorporating |
Beta Was this translation helpful? Give feedback.
Replies: 12 comments
-
Actually I don't get the ValueError here... I get an XDMF of which
and which ParaView is happy to open and render. |
Beta Was this translation helpful? Give feedback.
-
Ah, O. K., after updating meshio I do get the ValueError, which makes more sense, since, as you say, the lengths of the points and point-data don't match. |
Beta Was this translation helpful? Give feedback.
-
General issue on serialization of meshes and solutions: #158 |
Beta Was this translation helpful? Give feedback.
-
I'd guess you need to provide one value per |
Beta Was this translation helpful? Give feedback.
-
If you use matching from skfem import *
from skfem.helpers import dot, grad
m = MeshTri2().init_circle()
e = ElementTriP2()
basis = CellBasis(m, e)
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))
@LinearForm
def rhs(v, _):
return 1.0 * v
A = asm(laplace, basis)
b = asm(rhs, basis)
A, b = enforce(A, b, D=basis.get_dofs())
x = solve(A, b)
m.save("hosave.vtk", point_data={'x': x}) If this is not the case, then it becomes more complicated. |
Beta Was this translation helpful? Give feedback.
-
If your finite element space has higher order than P2, then L2-projection is well defined and you can use |
Beta Was this translation helpful? Give feedback.
-
Otherwise, you can hack together something like this: from skfem import *
from skfem.helpers import dot, grad
m = MeshTri2().init_circle()
e = ElementTriP1()
basis = CellBasis(m, e, intorder=4)
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))
@LinearForm
def rhs(v, _):
return 1.0 * v
A = asm(laplace, basis)
b = asm(rhs, basis)
A, b = enforce(A, b, D=basis.get_dofs())
x = solve(A, b)
X = basis.interpolate(x).value
y = projection(lambda _: X, basis.with_element(ElementTriP2()), basis)
m.save("hosave.vtk", point_data={'y': y}) Maybe |
Beta Was this translation helpful? Give feedback.
-
Actually, this seems to work also: from skfem import *
from skfem.helpers import dot, grad
m = MeshTri2().init_circle()
e = ElementTriP1()
basis = CellBasis(m, e, intorder=4)
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))
@LinearForm
def rhs(v, _):
return 1.0 * v
A = asm(laplace, basis)
b = asm(rhs, basis)
A, b = enforce(A, b, D=basis.get_dofs())
x = solve(A, b)
y = projection(x, basis.with_element(ElementTriP2()), basis)
m.save("hosave.vtk", point_data={'y': y}) I'm not sure what's going on. |
Beta Was this translation helpful? Give feedback.
-
My bad I should have mentioned my system setup. Thanks for pointing me to the issue of serialization. I will take a look at it |
Beta Was this translation helpful? Give feedback.
-
Interesting. Thanks. For projecting a vector field should I do it component-by-component? Because once I have the solution vector say of shape mesh.save("filename.xdmf", point_data = {"u": du_proj[projected_basis.nodal_dofs].T}) run into a similar issue? I am thinking out loud here -- haven't tested this yet. Will do so later today. |
Beta Was this translation helpful? Give feedback.
-
I think these formats that support quadratic meshes expect one value per node. However, I think you can use projection with |
Beta Was this translation helpful? Give feedback.
-
The global The problem I was solving is that of an incompressible shell radially stretched outwards and simply plotting the projected displacement seemed to suggest that something isn't correct. See
|
Beta Was this translation helpful? Give feedback.
The global
L2
projection withElementVector
seems to be a bit off. I don't know why though. Will try component-wiseL2
projection.The problem I was solving is that of an incompressible shell radially stretched outwards and simply plotting the projected displacement seemed to suggest that something isn't correct. See
L2
projectionwhereas if I simply write the
nodal
andedge
dofs, it seems to work (at least I don't see anything noticeably off).Actual
nodal
andedge
dataThe dof-ordering is quite nice in the sense that I can simply stich together the nodal and edge dofs by