Interpolating numerical values (2D numpy array) over a mesh

General
1

Hi,

I am trying to interpolate a 2D-vector valued numpy array over a mesh. I have the numerical values of 2D-vectors at every node of the mesh as a numpy array. Now I would like to interpolate this numerical vector valued function over the mesh. Could you please let know how I can do that? Below I post my code.

### Linear elasticity example FEniCSx tutorial ###

######################################################################
###         importing packages and defining parameters starts      ### 
######################################################################
from __future__ import annotations

from pathlib import Path
from typing import Union

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx.fem as fem
import numpy as np
import scipy.sparse.linalg
from dolfinx import default_scalar_type
from dolfinx.common import Timer, TimingType, list_timings
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx,
                 exp, grad, inner, pi, sin, cos)

import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

import pandas as pd 
import pyvista
from dolfinx import mesh, fem, plot, io

from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
import ufl

### importing the subroutine to give analytical solution ###
import stokesTheoryForL2 
from stokesTheoryForL2 import anasolmain
from stokesTheoryForL2 import correctedsol

#L = 1
#W = 0.2
mu = 0.001 ### 0.001, 1
#rho = 1
#delta = W / L
#gamma = 0.4 * delta**2
#beta = 1.25
#lambda_ = beta
lambda_ = 1

print("Value of mu: {}, and lambda: {}".format(mu,lambda_))
#g = gamma

######################################################################
###      importing packages and defining parameters ends           ### 
######################################################################


######################################################################
###                          defining mesh starts                  ### 
######################################################################

# domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
#                          [20, 6, 6], cell_type=mesh.CellType.hexahedron) 

# nx, ny = 50, 50
# domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])], 
#                                [nx, ny], mesh.CellType.triangle) 


# Get PETSc int and scalar types
complex_mode = True if np.dtype(default_scalar_type).kind == 'c' else False

# Create mesh and finite element
NX = 50
NY = 50
mesh = create_unit_square(MPI.COMM_WORLD, NX, NY)
#V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))
V = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
tol = 250 * np.finfo(default_scalar_type).resolution

# Listing space points ### 
dof_coordinates = V.tabulate_dof_coordinates()
#print(dof_coordinates[:2,:2])
df = pd.DataFrame(dof_coordinates)
df.to_csv('outdata/planenodesPeriodic.csv', index = None, header=False) 

######################################################################
###                     defining mesh ends                         ### 
######################################################################


######################################################################
###                boundary conditions starts                      ###
######################################################################
bcs = []
#########################################################################

def periodic_boundary(x):
    #print('From periodic_boundary0 function, length of x: ', len(x))
    print('From periodic_boundary0 function, shape of x: ', x.shape)
    #print(np.isclose(x[0], 1, atol=tol))
    boollist = []
    for i in range(x.shape[1]):
        if np.isclose(x[0,i], 0, atol=tol) or np.isclose(x[1,i], 0, atol=tol):
            boollist = boollist + [1]
        else:
            boollist = boollist + [0]
    boollist = np.array(boollist, dtype='bool')
    return boollist

def mapping(x):
    left_constraint = x[1] > 1e-14
    #mapped_values = np.zeros((2, x.shape[1]))
    mapped_values = np.zeros_like(x)
    mapped_values[0] = left_constraint + x[0]
    mapped_values[1] = np.invert(left_constraint) + x[1]
    return mapped_values


with Timer("~PERIODIC: Initialize MPC"):
    mpc = MultiPointConstraint(V)
    mpc.create_periodic_constraint_geometrical(V, periodic_boundary, mapping, bcs)
    #mpc.create_periodic_constraint_geometrical(V, periodic_boundary0, periodic_relation0, bcs)
    #mpc.create_periodic_constraint_geometrical(V, periodic_boundary1, periodic_relation1, bcs)
    mpc.finalize()


### traction free BC starts on other faces ### 
T = fem.Constant(mesh, default_scalar_type((0, 0)))
### traction free BC ends on other faces ### 


### defining surface integral over the boundary starts ###
ds = ufl.Measure("ds", domain=mesh)
### defining surface integral over the boundary ends ###

######################################################################
###                  boundary conditions end                       ###
######################################################################

######################################################################
###              variational formulation starts                    ###
######################################################################

def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = SpatialCoordinate(mesh)
#print('Printing x', x)
#f = as_vector((sin(2 * pi * x[0]), 0))
f = as_vector((1000 * cos(2 * pi * x[0] + 2 * pi * x[1]), 1000 * cos(2 * pi * x[0] + 2 * pi * x[1])))
#f = fem.Constant(mesh, default_scalar_type((0, 0)))
#f = as_vector((x[0] * sin(5.0 * pi * x[1]) + 1.0 * exp(-(dx_ * dx_ + dy_ * dy_) / 0.02), 0.3 * x[1]))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

######################################################################
###              variational formulation ends                      ###
######################################################################

######################################################################
###           solving variational problem starts                   ###
######################################################################
# Setup MPC system
with Timer("~PERIODIC: Initialize varitional problem"):
    problem = LinearProblem(a, L, mpc, bcs=bcs)

solver = problem.solver

# Give PETSc solver options a unique prefix
solver_prefix = "dolfinx_mpc_solve_{}".format(id(solver))
solver.setOptionsPrefix(solver_prefix)

petsc_options: dict[str, Union[str, int, float]]
if complex_mode or default_scalar_type == np.float32:
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
else:
    petsc_options = {"ksp_type": "cg", "ksp_rtol": 1e-6, "pc_type": "hypre", "pc_hypre_type": "boomeramg",
                     "pc_hypre_boomeramg_max_iter": 1, "pc_hypre_boomeramg_cycle_type": "v"  # ,
                     # "pc_hypre_boomeramg_print_statistics": 1
                     }

# Set PETSc options
opts = PETSc.Options()  # type: ignore
opts.prefixPush(solver_prefix)
if petsc_options is not None:
    for k, v in petsc_options.items():
        opts[k] = v
opts.prefixPop()
solver.setFromOptions()


with Timer("~PERIODIC: Assemble and solve MPC problem"):
    uh = problem.solve()
    # solver.view()
    it = solver.getIterationNumber()
    print("Constrained solver iterations {0:d}".format(it))

######################################################################
###             solving variational problem ends                   ###
######################################################################

######################################################################
###                  L^2 error computation starts                  ###
######################################################################

usol = uh.x.array
numpoints = int(len(usol)/2)
usol = np.reshape(usol, (numpoints,2))
df = pd.DataFrame(usol)
df.to_csv('outdata/planedisplacementsPeriodicNumerical.csv', index = None, header=False)


### interpolating the analytical solution over mesh using basis functions ###
V2 = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, ))) 
uex = fem.Function(V2) 
#uex.interpolate(lambda x: anasolmain(x[0],x[1]))
uex.interpolate(anasolmain)

### reshaping the 2D-planer analytical solution ###
#print(uex.x.array)
usolAna = uex.x.array
numpointsAna = int(len(usolAna)/2)
usolAna = np.reshape(usolAna, (numpointsAna,2))
#print(usol.tolist())

### saving the analytical solution into a file ###
df = pd.DataFrame(usolAna)
df.to_csv('outdata/planedisplacementsPeriodicAnaFenicsx.csv', index = None, header=False)

###=================================================================###
L2_error = fem.form(ufl.inner(uh - uex, uh - uex) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))
#print("L^2-norm of the error for {}x{} grid is: {}".format(NX, NY, error_L2))
###=================================================================###

###=================================================================###
L2_error = fem.form(uh[0] * ufl.dx)
error_local_uh0 = fem.assemble_scalar(L2_error)
#print("Integral of the first component of numerical solution: ", error_local_uh0)
###=================================================================###

###=================================================================###
L2_error = fem.form(uh[1] * ufl.dx)
error_local_uh1 = fem.assemble_scalar(L2_error)
#print("Integral of the second component of numerical solution: ", error_local_uh1)
###=================================================================###

###=================================================================###
L2_error = fem.form(uex[0] * ufl.dx)
error_local_uex0 = fem.assemble_scalar(L2_error)
#print("Integral of the first component of analytical solution: ", error_local_uex0)
###=================================================================###

###=================================================================###
L2_error = fem.form(uex[1] * ufl.dx)
error_local_uex1 = fem.assemble_scalar(L2_error)
#print("Integral of the second component of analytical solution: ", error_local_uex1)
###=================================================================###

# ######################################################################
# ###                  L^2 error computation ends                    ###
# ######################################################################

usolCorrected = np.zeros((numpoints,2))
### Subtracting the constant error from the numerical result
usolCorrected[:,0] = usol[:,0] - (error_local_uh0 - error_local_uex0)
usolCorrected[:,1] = usol[:,1] - (error_local_uh1 - error_local_uex1)

print("Printing corrected solution from outside any function len(usolCorrected[:,0])", len(usolCorrected[:,0]))
print("Printing corrected solution from outside any function", np.array([usolCorrected[:,0],usolCorrected[:,1]]))
df = pd.DataFrame(usolCorrected)
df.to_csv('outdata/planedisplacementsPeriodicCorrected.csv', index = None, header=False)

################################################################################################
###       Interpolating corrdected numerical solution over mesh using basis functions        ###
def usolcorrectSupply(x):
    x1 = x[0]
    x2 = x[1]

    v1anaCorrected = usolCorrected[:,0] #np.ones(len(x1)) #x1 #usolCorrected[:,0] #utest[:,0]
    v2anaCorrected = usolCorrected[:,1] #np.cos(2 * np.pi * x2) #x2 #usolCorrected[:,1] #utest[:,1]

    print("Length of v1anaCorrected (number of x[0] coordinates): ", len(x1))
    print("Printing corrected sol from inside of usolcorrectSupply function", np.array([v1anaCorrected, v2anaCorrected]))
    return np.array([v1anaCorrected, v2anaCorrected])

V2 = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
usolCorrectedMesh = fem.Function(V2)
usolCorrectedMesh.interpolate(usolcorrectSupply)
#usolCorrectedMesh.interpolate(anasolmain)
#usolCorrectedMesh.interpolate(correctedsol)


usolCorrectedPrint = usolCorrectedMesh.x.array
#print("Testing usolcurrectSupply function", usolCorrectedPrint)
numpoints = int(len(usolCorrectedPrint)/2)
usolCorrectedPrint = np.reshape(usolCorrectedPrint, (numpoints,2))
print("Printing interpolated corrected sol", np.array([usolCorrectedPrint[:,0],usolCorrectedPrint[:,1]]))

df = pd.DataFrame(usolCorrectedPrint)
#df = pd.DataFrame(usolCorrected)
df.to_csv('outdata/planedisplacementsPeriodicCorrectedMesh.csv', index = None, header=False)
#################################################################################################

###=================================================================###
L2_error = fem.form(usolCorrectedMesh[0] * ufl.dx)
error_local_ucorrect0 = fem.assemble_scalar(L2_error)
# error_L2 = np.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))
print("Integral of the first component of corrected solution: ", error_local_ucorrect0)
###=================================================================###

###=================================================================###
L2_error = fem.form(usolCorrectedMesh[1] * ufl.dx)
error_local_ucorrect1 = fem.assemble_scalar(L2_error)
# error_L2 = np.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))
print("Integral of the second component of corrected solution: ", error_local_ucorrect1)
###=================================================================###

print("Number of nodes {}".format(len(dof_coordinates)))
print("Number of sol points {}".format(len(usol)))

In the above code I am trying to interpolate the array “usolCorrected” (which is a numpy array of 2D-vectors evaluated at every node of the mesh) over the mesh.

Thanks.

1 Like
2

You typically want to search the forum before asking a question. Interpolating a function defined from an array is quite a common question, even though probably the case of a vector field is less common.

The very first search result I found is How to interpolate numpy array into VectorFunctionSpace? : it boils down to creating a LinearNDInterpolator from scipy. I won’t be searching further results for you.

Furthermore, taking the risk of sounding rude, I won’t be reading your code either. That is far from a minimal example: we don’t need to see your linear elasticity problem, which even uses dolfinx_mpc, for your to formulate the question “how can I interpolate a 2D field defined in numpy to a dolfinx Function”? Just take a built in mesh, prepare mock numpy data, and then ask how to interpolate those in a dolfinx Function.

3

This is far from a general example representative of your question:

Interpolating numerical values (2D numpy array) over a mesh

as the code involves dolfinx_mpc and solving a problem.
It is unclear why you need this for the problem at hand.

You should also consider the following:

and explain how your relation with the mesh nodes is assumed.

1 Like