Aligning dofs between subspace of mixed element and new function

General
1

Hi,
I’m working on some finite elasticity codes and I’ve ran into a small quirk and would love some and explanation.

In particular I’m curious about how Dofs of subspaces of mixed elements are ordered and how they are different than those created. Also if there is any reliability in the ordering of dofs. And lastly is it possible to align dofs between two functions in different function spaces?

To better explain what I’m talking about here is the following findings I saw

I was following the hyperplasticity demo where the mesh is warped by a vector but the coloring is done by a separate scalar.

This is my first attempt to plot the function:

pyvista.set_jupyter_backend('client')
plotter = pyvista.Plotter()

V , V_dofs = ME.sub(0).collapse() #Vector element subspace of ME
u_n = fem.Function(V)
u_ex = Expression(w.sub(0),V.element.interpolation_points())
u_n.interpolate(u_ex)

topology, cells, geometry = plot.create_vtk_mesh(V)

function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_n)] = u_n.x.array.reshape(geometry.shape[0], len(u_n))

function_grid["u"] = values
function_grid.set_active_vectors("u")

warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 67.5])

# Compute magnitude of displacement to visualize in GIF 
Vs = FunctionSpace(domain,("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(u_n[1], Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
#print(u_n.function_space.dofmap.list)
plotter.update_scalars(magnitude.x.array,render = False)
plotter.camera.position=[5,25,200]
plotter.camera.focal_point=[5,40,0]
plotter.render()
plotter.show()

Here is the result:

Screenshot 2023-09-10 at 8.17.52 PM
Screenshot 2023-09-10 at 8.17.52 PM1418×900 152 KB

It seems that the dofs between the Magnitude function and the u_n function are not in the same order.

I then created this example:

pyvista.set_jupyter_backend('client')
plotter = pyvista.Plotter()

V = fem.VectorFunctionSpace(domain,("Lagrange",2)) #HERE IS THE FIX JUST CREATE A NEW FUNCTION SPACE
u_n = fem.Function(V)
u_ex = Expression(w.sub(0),V.element.interpolation_points())
u_n.interpolate(u_ex)

topology, cells, geometry = plot.create_vtk_mesh(V)

function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_n)] = u_n.x.array.reshape(geometry.shape[0], len(u_n))

function_grid["u"] = values
function_grid.set_active_vectors("u")

warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 67.5])

# Compute magnitude of displacement to visualize in GIF 
Vs = FunctionSpace(domain,("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(u_n[1], Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
#print(u_n.function_space.dofmap.list)
plotter.update_scalars(magnitude.x.array,render = False)
plotter.camera.position=[5,25,200]
plotter.camera.focal_point=[5,40,0]
plotter.render()
plotter.show()

Screenshot 2023-09-10 at 8.18.53 PM
Screenshot 2023-09-10 at 8.18.53 PM1076×654 52.1 KB

This function seems to work correctly.

As seen above it seems that creating a new vector subspace solves the previous issue above issue and the DOFS between the two subspaces align. I’m curious if this is coincidence or if DOF ordering is assured.

Also is there anyway to get the first example to work? Maybe by aligning the dofs between the subspace and the new element?

Thanks for your help.

2

Please provide a minimal reproducible example, as it is unclear to me what function spaces you are using in your first example. (The code does not need to solve an actual elasticity problem to highlight any of these issues, just interpolate a function (say (x[0], x[1], x[2]) into the given function space.

3

Yup of course here are the quick two examples I’ve concocted

import numpy as np
import dolfinx

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot, log
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, Expression )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
import ufl
from ufl import (TestFunctions, TrialFunction, Identity, grad, det, div, dev, inv, tr, sqrt, conditional , gt, dx, inner, derivative, dot, ln, split)
from datetime import datetime
from dolfinx.plot import create_vtk_mesh

import pyvista
pyvista.set_jupyter_backend('client')

length = 10.0 # mm
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[length,length,length]],[2,2,2],mesh.CellType.hexahedron)

def Strech(x):
    return (-.5*x[0],x[1]*2,x[2])

U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2) # For displacement
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1) # For  pressure
TH = ufl.MixedElement([U2, P1])     # Taylor-Hood style mixed element
ME = FunctionSpace(domain, TH)    # Total space for all DOFs
w = Function(ME)



u = w.sub(0)


u.interpolate(Strech)

pyvista.set_jupyter_backend('client')
plotter = pyvista.Plotter()

V , V_dofs = ME.sub(0).collapse() ## Difference 
u_n = fem.Function(V)
u_ex = Expression(w.sub(0),V.element.interpolation_points())
u_n.interpolate(u_ex)

topology, cells, geometry = plot.create_vtk_mesh(V)

function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_n)] = u_n.x.array.reshape(geometry.shape[0], len(u_n))

function_grid["u"] = values
function_grid.set_active_vectors("u")

warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])

# Compute magnitude of displacement to visualize in GIF 
Vs = FunctionSpace(domain,("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(u_n[1], Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
#print(u_n.function_space.dofmap.list)
plotter.update_scalars(magnitude.x.array,render = False)
plotter.camera.position=[5,25,200]
plotter.camera.focal_point=[5,40,0]
plotter.render()
plotter.show()

With the result

Screenshot 2023-09-11 at 10.53.09 AM
Screenshot 2023-09-11 at 10.53.09 AM932×988 127 KB

And

import numpy as np
import dolfinx

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot, log
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, Expression )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
import ufl
from ufl import (TestFunctions, TrialFunction, Identity, grad, det, div, dev, inv, tr, sqrt, conditional , gt, dx, inner, derivative, dot, ln, split)
from datetime import datetime
from dolfinx.plot import create_vtk_mesh

import pyvista
pyvista.set_jupyter_backend('client')

length = 10.0 # mm
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[length,length,length]],[2,2,2],mesh.CellType.hexahedron)

def Strech(x):
    return (-.5*x[0],x[1]*2,x[2])

U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2) # For displacement
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1) # For  pressure
TH = ufl.MixedElement([U2, P1])     # Taylor-Hood style mixed element
ME = FunctionSpace(domain, TH)    # Total space for all DOFs
w = Function(ME)



u = w.sub(0)


u.interpolate(Strech)

pyvista.set_jupyter_backend('client')
plotter = pyvista.Plotter()

V = fem.VectorFunctionSpace(domain,("Lagrange",2)) ## Difference 
u_n = fem.Function(V)
u_ex = Expression(w.sub(0),V.element.interpolation_points())
u_n.interpolate(u_ex)

topology, cells, geometry = plot.create_vtk_mesh(V)

function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_n)] = u_n.x.array.reshape(geometry.shape[0], len(u_n))

function_grid["u"] = values
function_grid.set_active_vectors("u")

warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 20])

# Compute magnitude of displacement to visualize in GIF 
Vs = FunctionSpace(domain,("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(u_n[1], Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
#print(u_n.function_space.dofmap.list)
plotter.update_scalars(magnitude.x.array,render = False)
plotter.camera.position=[5,25,200]
plotter.camera.focal_point=[5,40,0]
plotter.render()
plotter.show()

With the correct result:

Screenshot 2023-09-11 at 10.53.45 AM
Screenshot 2023-09-11 at 10.53.45 AM1054×1032 177 KB

The only difference is that in one we use
V , V_dofs = ME.sub(0).collapse() ## Difference
while in the other that produces correct results we use
V = fem.VectorFunctionSpace(domain,("Lagrange",2)) ## Difference

4

Use
Vs, _ = V.sub(0).collapse()
instead, as you then get the correct ordering for a component of the vector space.

5

Here is this example doing that. It does however still seem to have the incorrect ordering for the vector space. See the example below.

Could this just be bug in the current version of DOLFINx? I’m using DOLFINx 7.0 the main branch. Just to check I got the latest version and pulled a couple minutes ago.

import numpy as np
import dolfinx

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot, log
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, Expression )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
import ufl
from ufl import (TestFunctions, TrialFunction, Identity, grad, det, div, dev, inv, tr, sqrt, conditional , gt, dx, inner, derivative, dot, ln, split)
from datetime import datetime
from dolfinx.plot import vtk_mesh

import pyvista
pyvista.set_jupyter_backend('client')

length = 10.0 # mm
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[length,length,length]],[2,2,2],mesh.CellType.hexahedron)

def Strech(x):
    return (-.5*x[0],x[1]*2,x[2])

U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2) # For displacement
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1) # For  pressure
TH = ufl.MixedElement([U2, P1])     # Taylor-Hood style mixed element
ME = FunctionSpace(domain, TH)    # Total space for all DOFs
w = Function(ME)



u = w.sub(0)


u.interpolate(Strech)

pyvista.set_jupyter_backend('client')
plotter = pyvista.Plotter()

V , _ = ME.sub(0).collapse() ## Difference 
u_n = fem.Function(V)
u_ex = Expression(w.sub(0),V.element.interpolation_points())
u_n.interpolate(u_ex)

topology, cells, geometry = plot.vtk_mesh(V)

function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_n)] = u_n.x.array.reshape(geometry.shape[0], len(u_n))

function_grid["u"] = values
function_grid.set_active_vectors("u")

warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])

# Compute magnitude of displacement to visualize in GIF 
Vs = FunctionSpace(domain,("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(u_n[1], Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
#print(u_n.function_space.dofmap.list)
plotter.update_scalars(magnitude.x.array,render = False)
plotter.camera.position=[5,25,200]
plotter.camera.focal_point=[5,40,0]
plotter.render()
plotter.show()

Screenshot 2023-09-12 at 7.59.44 PM
Screenshot 2023-09-12 at 7.59.44 PM1128×1032 198 KB

6

This is not what I stated above.
I stated

7

Oh sorry I misunderstood.
As expected that works just fine.

I wanted to express my appreciation for your time not just helping me but really supporting a huge part of the community.

One last question for my curiosity, is why the dof ordering is between functions space and the collapsed subspaces of the mixed functions. Also is there anyway to reorder the dofs?

import numpy as np
import dolfinx

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot, log
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, Expression )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
import ufl
from ufl import (TestFunctions, TrialFunction, Identity, grad, det, div, dev, inv, tr, sqrt, conditional , gt, dx, inner, derivative, dot, ln, split)
from datetime import datetime
from dolfinx.plot import vtk_mesh

import pyvista
pyvista.set_jupyter_backend('client')

length = 10.0 # mm
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[length,length,length]],[2,2,2],mesh.CellType.hexahedron)

def Strech(x):
    return (-.5*x[0],x[1]*2,x[2])

U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2) # For displacement
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1) # For  pressure
TH = ufl.MixedElement([U2, P1])     # Taylor-Hood style mixed element
ME = FunctionSpace(domain, TH)    # Total space for all DOFs
w = Function(ME)



u = w.sub(0)


u.interpolate(Strech)

pyvista.set_jupyter_backend('client')
plotter = pyvista.Plotter()

V , _ = ME.sub(0).collapse() ## Difference 
u_n = fem.Function(V)
u_ex = Expression(w.sub(0),V.element.interpolation_points())
u_n.interpolate(u_ex)

topology, cells, geometry = plot.vtk_mesh(V)

function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_n)] = u_n.x.array.reshape(geometry.shape[0], len(u_n))

function_grid["u"] = values
function_grid.set_active_vectors("u")

warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])

# Compute magnitude of displacement to visualize in GIF 
Vs,_ = V.sub(0).collapse()
magnitude = fem.Function(Vs)
us = fem.Expression(u_n[1], Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
#print(u_n.function_space.dofmap.list)
plotter.update_scalars(magnitude.x.array,render = False)
plotter.camera.position=[5,25,200]
plotter.camera.focal_point=[5,40,0]
plotter.render()
plotter.show()

dolfinx.io.VTXWriter(MPI.COMM_WORLD,"test",magnitude)

Screenshot 2023-09-13 at 12.50.42 PM
Screenshot 2023-09-13 at 12.50.42 PM794×720 51.3 KB