Numerical values from ufl.SpatialCoordinate

Hi everyone!

I am implementing a Lubrication + Wear model in FeniCSx for my Masters Thesis.

While i was implementing further improvements i noticed, my previous approach will not be suitable to implement something i would need.

What I am trying to do is the following:

In my variational formulation i have a call to a height function, that returns the separation between the two geometries I am simulating.

Inside this function, i would want to obtain the value from a dictionary, which maps each node in the mesh to a certain value which changes over time.

However, i am unable to obtain the numerical coordinates (which I would need, as they are the key to the dictionary), as when solving, this height function is called with my ufl.SpatialCoordinate object.

There are probably many ways of doing such, do you have any suggestions? For context, the reason i am using a dictionary is because there is no analytic expression to calculate the value i map.

Many thanks!

Code:

from abc import ABC, abstractmethod ## Para clases abstractas
import ufl ## Parte de FeniCSx
from dolfinx import fem, io, mesh, plot, geometry ## Parte de FeniCSx
from mpi4py import MPI ## Parte de FeniCSx
from petsc4py.PETSc import ScalarType ## Parte de FeniCSx
import math, random ## Operaciones matem√°ticas b√°sicas
import numpy as np ## Numpy para poder operar con arrays
import matplotlib.pyplot as mpl ## Para pintamiento

class Problema: 

    ## Constructor
    def __init__(self, datos):
        self.datos = datos ## Diccionario que contiene los datos del problema
        self.deltat = 0.1 ## Placeholder 
        self.datos["speed"] = ((self.datos["Ut"] - self.datos["Ub"])**2 +  (self.datos["Vt"] - self.datos["Vb"])**2)**(1/2)
        self.BC = []  ## Array que contiene las condiciones de contorno
        ## Diccionario que contiene los valores de diversos parámetros para la iteración inicial. 
        self.timestepvals= {"Zh": 3, 
                                              "DeltaZh": 0,
                                              "DeltaZ0": 0, 
                                              "S1" : 0} 
    
    ## Setter para la variable h0 para el calculo de la altura. Necesario para poder variar la altura en calculos iterativos. 
    def set_h0(self,h0):
        self.datos["h0"] = h0
    
    ## Metodo abstracto para el calculo de la altura. Obligamos a las clases hijas a definirla de forma que represente el problema
    @abstractmethod
    def height(self):
        pass

    ## Clase para problemas tipo esfera - point contact
class ProblemaRodillo(Problema):

    ## Override del constructor para calcular el área del contacto de forma específica para cada caso. 
    def __init__(self, datos):
        Problema.__init__(self, datos) # Llamamos primero al constructor por defecto 
        self.datos["Area_nom"] = math.pi* self.datos["rRod"]**2 

    def height(self,x): 
        
        ## THIS VALUE SHOULD BE SOMETHING LIKE
        # worn = self.worn[(x[0], x[1], x[2])].
        worn = 0.001 

        return self.datos["h0"] + (x[0]**2)/(2*self.datos["rRod"]) # + worn
    
datosPruebaPsi={
    "rRod": 0.03, 
    "Load": 10e4,
    "h0min" : 5e-6,
    "h0max": 6e-6,
    "rho": 860,
    "mu": 0.1,
    "Ub":1,
    "Ut":0,
    "Vb":0,
    "Vt":0,
    
    "domain":{
        "xMin":-0.05,
        "xMax":0.05,
        "yMin":-0.5,
        "yMax":0.5,
        "numNodes":50
    }, 
}

Prueba = ProblemaRodillo(datosPruebaPsi)

xMin = Prueba.datos["domain"]["xMin"]
xMax = Prueba.datos["domain"]["xMax"]
yMin = Prueba.datos["domain"]["yMin"]
yMax = Prueba.datos["domain"]["yMax"]
numNodes = Prueba.datos["domain"]["numNodes"]
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                    points=((xMin, yMin), (xMax, yMax)), n=(numNodes, numNodes),
                    cell_type=mesh.CellType.triangle)

Prueba.FUNCSPACE = fem.FunctionSpace(msh, ("Lagrange", 2))
Prueba.PSPACE = ufl.TrialFunction(Prueba.FUNCSPACE)
Prueba.VSPACE = ufl.TestFunction(Prueba.FUNCSPACE)


# Definimos las fronteras para aplicar la condicion de contorno
facets = mesh.locate_entities_boundary(msh, dim=1, marker=lambda x: np.logical_or(np.logical_or(np.isclose(x[0],Prueba.datos["domain"]["xMin"]),np.isclose(x[0], Prueba.datos["domain"]["xMax"])),
                            np.logical_or(np.isclose(x[1], Prueba.datos["domain"]["yMin"]),np.isclose(x[1],Prueba.datos["domain"]["yMax"]))))
# Localizamos los grados de libertad en las fronteras                                                                  
dofs = fem.locate_dofs_topological(V=Prueba.FUNCSPACE, entity_dim=1, entities=facets)
# Aplicamos condiciones de contorno tipo dirichlet
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=Prueba.FUNCSPACE)
## Introducimos en el array que contiene las condiciones de contorno
Prueba.set_h0(Prueba.datos["h0max"])
rho = Prueba.datos["rho"]
mu = Prueba.datos["mu"]
Um = 0.5*(Prueba.datos["Ub"]+Prueba.datos["Ut"])
Vm = 0.5*(Prueba.datos["Vb"]+Prueba.datos["Vt"])

x = ufl.SpatialCoordinate(msh)
p = Prueba.PSPACE
v = Prueba.VSPACE

## Formulación variacional
a = (rho*(Prueba.height(x)**3)/(12*mu))*ufl.inner(ufl.grad(p),ufl.grad(v))*ufl.dx
L = rho*Prueba.height(x)*ufl.inner(ufl.tensors.as_vector([Um,Vm]),ufl.grad(v))*ufl.dx 

## Asignamos el problema a nuestra clase
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "ibcgs", "pc_type":"lu"})
Prueba.solution = problem.solve()

points = np.zeros((3, len(msh.geometry.x[:,0])))
points[0] = msh.geometry.x[:,0]

bb_tree = geometry.BoundingBoxTree(msh, msh.topology.dim)
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
colliding_cells = geometry.compute_colliding_cells(msh, cell_candidates, points.T)
cells = [colliding_cells.links(i)[0] for i, point in enumerate(points.T)]

p_values = Prueba.solution.eval(points.T, cells)

mpl.scatter(points[0,:], p_values, label="Numerical solution")

ufl.SpatialCoordinate does not imply the ‚Äúmesh-nodes‚ÄĚ (meaning the vertices for linear meshes). The x=ufl.SpatialCoordinate(mesh) represents any quadrature point used during integration in its physical space.

What I would suggest to do, is to interpolate your function (non-analytical) into a suitable function space. Say, if you want to actually use the mesh nodes (vertices), you can create a first order Lagrange function and interpolate into that, i.e.

V  = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(V)

def your_expression(x):
     # map your values here
     return value
u.interpolate(your_expression)

If you instead you want to use the quadrature points (as described by x=ufl.SpatialCoordinate(mesh), I suggest something like the following:

from IPython import embed
import dolfinx
from mpi4py import MPI
import ufl
import basix
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)


def height(x):
    # Add whatever function you want to call here, for now choose identity
    return x[0]


quadrature_degree = 4
dx = ufl.Measure("dx", domain=mesh, metadata={
                 "quadrature_degree": quadrature_degree})
quadrature_points, wts = basix.make_quadrature(
    basix.cell.string_to_type(mesh.ufl_cell().cellname()), quadrature_degree)
Q_element = ufl.FiniteElement(
    "Quadrature", mesh.ufl_cell(), quadrature_degree, quad_scheme="default")
Q = dolfinx.fem.FunctionSpace(mesh, Q_element)
height_function = dolfinx.fem.Function(Q)

x = ufl.SpatialCoordinate(mesh)
x_expr = dolfinx.fem.Expression(x, quadrature_points)
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
integration_points = x_expr.eval(np.arange(num_cells, dtype=np.int32))

num_points_per_cell = quadrature_points.shape[0]
gdim = mesh.geometry.dim
assert(integration_points.shape[1] == num_points_per_cell*gdim)
assert(integration_points.shape[0] == num_cells)
for i in range(num_cells):
    for j in range(num_points_per_cell):
        coord = integration_points[i, gdim*j:gdim*(j+1)]
        value = height(coord)
        height_function.x.array[i*num_points_per_cell +
                                j:i*num_points_per_cell + (j+1)] = value

for i in range(mesh.geometry.dim):
    # verify that height_function maps to x[0]
    print(dolfinx.fem.assemble_scalar(
        dolfinx.fem.form((height_function-x[0])**2*dx)))


# Use heigh function in general()
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = dolfinx.fem.form(height_function*ufl.inner(ufl.grad(u), ufl.grad(v))*dx)
A = dolfinx.fem.petsc.assemble_matrix(a)
1 Like

Good evening!

Many thanks for the response, @dokken! With your hints I have been able to achieve what I wanted.

I tried the interpolation route, but i was getting nonsensical values ( e-311, i suppose because of the mapping i was doing)

While testing this, i tried to manually assign the fem.Function.vector.array, and realized that that is exactly the structure that i need to save the data. I need a way to assign a value to each DOF, thus the vector of a Function fits perfectly

It is probably not the intended usage, but after some quick tests i am confident I will be able to work like that.

Again, thanks a lot!!