Point data in high order meshes

[bhaveshshrimali]

How would one go about saving point_data for a 10-noded tetrahedral mesh for instance? I have right now, basis.nodal_dofs and basis.edge_dofs. I was using linear meshes up until now and hence basis.nodal_dofs sufficed even when having a high order approximation such as in #651

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 edge_dofs in saving point_data?

[gdmcbain]

Actually I don't get the ValueError here... I get an XDMF of which meshio-info reports

<meshio mesh object>
  Number of points: 29521
  Number of cells:
    tetra10: 19900
  Point data: u

and which ParaView is happy to open and render.

[gdmcbain]

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.

[gdmcbain]

General issue on serialization of meshes and solutions: #158

[kinnala]

I'd guess you need to provide one value per Mesh.doflocs. I'll try it out in a minute.

[kinnala]

If you use matching Element and Mesh, then you can simply pass the solution vector to point_data:

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.

[kinnala]

If your finite element space has higher order than P2, then L2-projection is well defined and you can use skfem.utils.projection to find the array which is passed to point_data.

[kinnala]

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 projection should be extended so that this lambda-hack is not required when interpolating?

[kinnala]

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.

[bhaveshshrimali]

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.

My bad I should have mentioned my system setup. meshio version is 4.4.5 and pygmsh 6.1.1 (I still feel that it is more intuitive and the newer versions can be avoided directly in favor of the gmsh python API.

Thanks for pointing me to the issue of serialization. I will take a look at it

[bhaveshshrimali]

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 projection should be extended so that this lambda-hack is not required when interpolating?

Interesting. Thanks. For projecting a vector field should I do it component-by-component? Because once I have the solution vector say of shape (3 * mesh.p.shape[1]) then wouldn't doing

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.

[kinnala]

I think these formats that support quadratic meshes expect one value per node. However, I think you can use projection with ElementVector. Just when passing data to point_data I think it must be component-by-component. Not sure though, I'm mostly using matplotlib for visualization.

[bhaveshshrimali]

The global L2 projection with ElementVector seems to be a bit off. I don't know why though. Will try component-wise L2 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 projection

whereas if I simply write the nodal and edge dofs, it seems to work (at least I don't see anything noticeably off).

Actual nodal and edge data

The dof-ordering is quite nice in the sense that I can simply stich together the nodal and edge dofs by

nodal_disps = du[basis["u"].nodal_dofs].T
edge_disps = du[basis["u"].edge_dofs].T

pointData = np.concatenate((
    nodal_disps, edge_disps
), axis=0)

I will close this since it seems I can get what I wanted to begin with. The complete code is basically example 36 but for different geometry and b.c's and now with ElementTetCCR

import numpy as np
from skfem.io import from_meshio
import os, meshio
from skfem.helpers import grad
from pygmsh import generate_mesh as gm
from pygmsh.opencascade import Geometry as geo
from numba import njit, prange
from scipy.sparse import bmat
from skfem.utils import projection
from skfem import * 

def generateMesh():
    geo_file = os.path.join(f"octant.geo")
    R_out = 1.0e-2
    fo = 0.05
    R_in = R_out * fo ** (1.0 / 3)
    geom = geo(R_out / 30.0, R_out / 10.)
    outer_ball = geom.add_ball([0.0, 0.0, 0.0], R_out)
    inner_ball = geom.add_ball([0.0, 0.0, 0.0], R_in)
    shell = geom.boolean_difference([outer_ball], [inner_ball])
    cube = geom.add_box([0.0, 0.0, 0.0], [R_out, R_out, R_out * 1.1])
    octant = geom.boolean_intersection([shell, cube])
    geom.add_raw_code("Mesh.ElementOrder=2;")
    geom.add_raw_code(f"Rout={R_out};")
    geom.add_raw_code(f"eps={R_out/1.e5};")
    geom.add_raw_code(
        "topline() = Curve In BoundingBox {-eps, - eps, -eps, eps, eps, Rout+eps};"
    )

    geom.add_raw_code('Physical Curve("topedge") = topline();')
    geom.add_raw_code('Physical Volume("entireDomain") = bo2[];')
    msh = gm(geom, geo_filename=geo_file)
    return msh

mesh = from_meshio(generateMesh())

@njit(nogil=True)
def Jstar(p):
    return (lmbda + p + np.sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda)

@njit(nogil=True)
def d2Jst_dp2(p):
    f = 2.*mu/(4.*lmbda*mu + (lmbda + p)**2.)**(3./2.)
    return f

@njit(nogil=True)
def dJst_dp(p):
    f = (lmbda + p + np.sqrt(4.*lmbda*mu + (lmbda + p)**2.))/(2.*lmbda*np.sqrt(4.*lmbda*mu + (lmbda + p)**2.))
    return f

@njit(nogil=True, parallel=True)
def vdet(F):
    J = np.zeros_like(F[0,0])
    for a in prange(J.shape[0]):
        for b in prange(J.shape[1]):
            J[a,b] += F[0, 0, a, b] * (F[1, 1, a, b] * F[2, 2, a, b] -
                                        F[1, 2, a, b] * F[2, 1, a, b]) -\
                    F[0, 1, a, b] * (F[1, 0, a, b] * F[2, 2, a, b] -
                                        F[1, 2, a, b] * F[2, 0, a, b]) +\
                    F[0, 2, a, b] * (F[1, 0,a ,b] * F[2, 1, a, b] -
                                        F[1, 1, a, b] * F[2, 0, a, b])
    return J


@njit(nogil=True, parallel=True)
def vinv(F):
    J = vdet(F)
    Finv = np.zeros_like(F)
    for a in prange(J.shape[0]):
        for b in prange(J.shape[1]):
            Finv[0, 0, a, b] += (-F[1, 2, a, b] * F[2, 1, a, b] +
                                F[1, 1, a, b] * F[2, 2, a, b]) / J[a, b]
            Finv[1, 0, a, b] += (F[1, 2, a, b] * F[2, 0, a, b] -
                                F[1, 0, a, b] * F[2, 2, a, b]) / J[a, b]
            Finv[2, 0, a, b] += (-F[1, 1, a, b] * F[2, 0, a, b] +
                                F[1, 0, a, b] * F[2, 1, a, b]) / J[a, b]
            Finv[0, 1, a, b] += (F[0, 2, a, b] * F[2, 1, a, b] -
                                F[0, 1, a, b] * F[2, 2, a, b]) / J[a, b]
            Finv[1, 1, a, b] += (-F[0, 2, a, b] * F[2, 0, a, b] +
                                F[0, 0, a, b] * F[2, 2, a, b]) / J[a, b]
            Finv[2, 1, a, b] += (F[0, 1, a, b] * F[2, 0, a, b] -
                                F[0, 0, a, b] * F[2, 1, a, b]) / J[a, b]
            Finv[0, 2, a, b] += (-F[0, 2, a, b] * F[1, 1, a, b] +
                                F[0, 1, a, b] * F[1, 2, a, b]) / J[a, b]
            Finv[1, 2, a, b] += (F[0, 2, a, b] * F[1, 0, a, b] -
                                F[0, 0, a, b] * F[1, 2, a, b]) / J[a, b]
            Finv[2, 2, a, b] += (-F[0, 1, a, b] * F[1, 0, a, b] +
                                F[0, 0, a, b] * F[1, 1, a, b]) / J[a, b]
    return Finv, J

@njit(nogil=True, parallel=True)
def F1(dv, du, p):
    
    out = np.zeros_like(du[0, 0])
    F = np.zeros_like(du)
    for a in prange(du.shape[2]):
        for b in prange(du.shape[3]):
            F[0, 0, a, b] += 1.
            F[1, 1, a, b] += 1.
            F[2, 2, a, b] += 1.
    
    F += du
    Finv, J = vinv(F)
    # p * J * transpose(Finv) + mu * F
    for i in range(du.shape[0]):
        for j in range(du.shape[1]):
            for a in prange(J.shape[0]):
                for b in prange(J.shape[1]):
                    out[a,b] += (mu * F[i, j, a, b] + p[a,b] * J[a,b] * Finv[j, i, a, b]) * dv[i, j, a, b]
    return out

@njit(nogil=True, parallel=True)
def F2(dv, du, p):
    # p = dp.value
    out = np.zeros_like(du[0,0])
    
    F = np.zeros_like(du)
    
    for a in prange(du.shape[2]):
        for b in prange(du.shape[3]):
            F[0, 0, a, b] += 1.
            F[1, 1, a, b] += 1.
            F[2, 2, a, b] += 1.
    
    F += du
    J = vdet(F)
   
    Js = Jstar(p)
    # J - (Js + (p + mu/Js - lmbda*(Js - 1))*dJst_dp(p))
    for a in prange(p.shape[0]):
        for b in prange(p.shape[1]):
            out[a, b] += (J[a,b] - (Js[a,b] + (p[a,b] + mu/Js[a,b] - lmbda*(Js[a,b] - 1.)) * dJst_dp(p)[a, b])) * dv[a,b]

    return out

@njit(nogil=True,parallel=True)
def A11(du, dv, dw, p):
    # p = dp.value
    out = np.zeros_like(p)
    F = np.zeros_like(dw)
    for a in prange(dw.shape[2]):
        for b in prange(dw.shape[3]):
            F[0, 0, a, b] += 1.
            F[1, 1, a, b] += 1.
            F[2, 2, a, b] += 1.
    
    F += dw
    Finv, J = vinv(F)
    kron = np.eye(dw.shape[0])
    
    for i in range(du.shape[0]):
        for j in range(du.shape[1]):
            for k in range(du.shape[0]):
                for l in range(du.shape[1]):
                    for a in range(J.shape[0]):
                        for b in range(J.shape[1]):
                            out[a, b] += (
                                p[a,b] * J[a,b] * Finv[l,k,a,b] * Finv[j, i, a, b] -\
                                    p[a,b] * J[a,b] * Finv[j,k,a,b] * Finv[l, i, a, b] +\
                                        mu * kron[i,k] * kron[j,l]
                            ) * du[i,j,a,b] * dv[k,l,a,b]

    return out

@njit(nogil=True, parallel=True)
def A12(du, dv, dw):
    out = np.zeros_like(dw[0,0])
    F = np.zeros_like(dw)
    for a in prange(dw.shape[2]):
        for b in prange(dw.shape[3]):
            F[0, 0, a, b] += 1.
            F[1, 1, a, b] += 1.
            F[2, 2, a, b] += 1.
    
    F += dw
    Finv, J = vinv(F)

    for i in range(dw.shape[0]):
        for j in range(dw.shape[1]):
            for a in prange(J.shape[0]):
                for b in prange(J.shape[1]):
                    out[a, b] += (
                        J[a,b] * Finv[j, i, a, b]
                    ) * dv[i, j, a, b] * du[a, b]
 
    return out

@njit(nogil=True, parallel=True)
def A22(du, dv, p):
    Js = Jstar(p)
    dJdp = dJst_dp(p)
    d2Jdp2 = d2Jst_dp2(p)
    
    out = np.zeros_like(p)
    for a in prange(p.shape[0]):
        for b in prange(p.shape[1]):
            out[a,b] += (-2. * dJdp[a, b] - p[a,b] * d2Jdp2[a,b] + mu/Js[a,b]**2 * dJdp[a,b]**2 - mu/Js[a,b] * d2Jdp2[a,b] + lmbda * (Js[a,b] - 1.) * d2Jdp2[a,b] + lmbda * dJdp[a,b]**2) * du[a,b] * dv[a,b]
    return out

@LinearForm(nthreads=8)
def a1(v, w):
    return F1(*(grad(v), grad(w["disp"]), w["press"].value))

@LinearForm(nthreads=8)
def a2(v, w):
    return F2(*(v.value, grad(w["disp"]), w["press"].value))

@BilinearForm(nthreads=8)
def b11(u, v, w):
    return A11(*(grad(u), grad(v), grad(w["disp"]), w["press"].value))

@BilinearForm(nthreads=8)
def b12(u, v, w):
    return A12(*(u.value, grad(v), grad(w["disp"])))

@BilinearForm(nthreads=8)
def b22(u, v, w):
    return A22(u.value, v.value, w["press"].value)

mu, lmbda = 0.01462, 14.62
stretch_ = 0.05

uelem = ElementVectorH1(ElementTetCCR())
pelem = ElementTetDG(ElementTetP1())

elems = {
    "u": uelem,
    "p": pelem
}
basis = {
    field: CellBasis(mesh, e, intorder=7)
    for field, e in elems.items()
}

du = basis["u"].zeros()
dp = basis["p"].zeros()

outer_nodes = lambda x: np.isclose(np.sqrt(x[0]**2 + x[1]**2 + x[2]**2)/r, 1., rtol=1.e-3)
dofs_outer_radius = basis["u"].find_dofs({
    "out": mesh.nodes_satisfying(outer_nodes)
})

u1Outer = lambda x,y,z: stretch_ * x
u2Outer = lambda x,y,z: stretch_ * y
u3Outer = lambda x,y,z: stretch_ * z


dofs_u1_outer = np.hstack((
    dofs_outer_radius["out"].nodal["u^1"],
    dofs_outer_radius["out"].edge["u^1"]
))
dofs_u2_outer = np.hstack((
    dofs_outer_radius["out"].nodal["u^2"],
    dofs_outer_radius["out"].edge["u^2"]
))
dofs_u3_outer = np.hstack((
    dofs_outer_radius["out"].nodal["u^3"],
    dofs_outer_radius["out"].edge["u^3"]
))

du[dofs_u1_outer] = u1Outer(*basis["u"].doflocs[:, dofs_u1_outer])
du[dofs_u2_outer] = u2Outer(*basis["u"].doflocs[:, dofs_u2_outer])
du[dofs_u3_outer] = u3Outer(*basis["u"].doflocs[:, dofs_u3_outer])

I = np.hstack((
    basis["u"].complement_dofs(dofs_outer_radius),
    basis["u"].N + np.arange(basis["p"].N)
))

for itr in range(50):
    uv = basis["u"].interpolate(du)
    pv = basis["p"].interpolate(dp) 

    K11 = asm(b11, basis["u"], basis["u"], disp=uv, press=pv)
    K12 = asm(b12, basis["p"], basis["u"], disp=uv, press=pv)
    K22 = asm(b22, basis["p"], basis["p"], disp=uv, press=pv)
    f = np.concatenate((
        asm(a1, basis["u"], disp=uv, press=pv),
        asm(a2, basis["p"], disp=uv, press=pv)
    ))
    K = bmat(
        [[K11, K12],
        [K12.T, K22]], "csr"
    )
    uvp = solve(*condense(K, -f, I=I), use_umfpack=True)
    delu, delp = np.split(uvp, [du.shape[0]])
    du += delu
    dp += delp
    normu = np.linalg.norm(delu)
    normp = np.linalg.norm(delp)
    print(f"{itr+1}, norm_du: {normu}, norm_dp: {normp}")
    if normu < 1.e-8 and normp < 1.e-8:
        break

nodal_disps = du[basis["u"].nodal_dofs].T
edge_disps = du[basis["u"].edge_dofs].T

pointData = np.concatenate((
    nodal_disps, edge_disps
), axis=0)
mesh.save("octant.xdmf", point_data={"u": pointData})