Circular mean interpolation in Fenicsx?

General
1

Hello,
I was wondering if there is any way to perform a circular mean interpolation in Fenicsx, as the one shown in the screenshot?
My problem is expressed in terms of angular coordinates theta1 and theta2, which are defined over the intervals -pi<=theta2<=pi and 0<=theta1<=pi. Arithmetical interpolations lead to unphysical results if the nodal angles theta2 are close to the upper/lower limits of the definition range.
Can you please point me to an example where a user-defined interpolation is performed?
Many thanks!
Screenshot 2023-12-04 at 11.41.15

1 Like
2

Could you be more specific in your question and perhaps provide a minimal working example of what you’ve tried so far? The concept of a circular mean is likely unfamiliar to the FEniCS crowd; however, if the concepts are straightforward in the context of a FEniCS implementation we could provide help.

3

Hi Natan,
thank you for your kind reply.
The variational form of my problem is expressed in terms of the unknown vector theta=(theta1,theta2).
The domain of each of those functions is -pi<=theta2<=pi and -pi/2<=theta1<=pi/2
Referring to the equations in the screenshot above:

  • theta1^I and eta1^I respectively denote the values of the unknown and test function at the node I

  • theta1^h and eta1^h denote the corresponding interpolated values

  • N^I are the basis functions, for instance of Lagrange elements

In the following MWE, what I am doing under the hood is discretizing both theta1 and theta2 as:

theta1^h = sum_I (N_I theta1^I)
theta2^h = sum_I (N_I theta2^I)

(same for the corresponding test functions), while what I would like to do is expressing theta2 as a circular interpolation, as indicated in the screenshot above.

I am attaching a MWE, which I hope would allow me to clarify the problem.
Thank you for your help, I really appreciate all the work that you guys do.

import os
import sys
import numpy as np
import petsc4py
import ufl
from dolfinx import fem, mesh, la
from mpi4py import MPI
from petsc4py import PETSc
petsc4py.init(sys.argv)
os.system('clear')
# Create MPI communicator
comm = MPI.COMM_WORLD
# ----- Create the mesh
L = 1.0
N = 10
domain = mesh.create_box(comm=comm, points=[[0.0, 0.0, 0.0], [L, L, L]],\
                         n=[N, N, N], cell_type=mesh.CellType.tetrahedron)
metadata = {'quadrature_degree' : 4}
dv = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define the function space
V_theta = fem.VectorFunctionSpace(domain,("CG", 1), dim=2)
# Solution, Test, and Trial functions
theta = fem.Function(V_theta)
theta_test, theta_trial =  ufl.TestFunction(V_theta), ufl.TrialFunction(V_theta)
def initialize(x):
    d = x.shape[1]
    angles = [np.zeros(d), np.zeros(d)]
    # Create noise
    np.random.seed(123)
    noise = abs(np.random.normal(np.mean(np.ones(d)*0), 1/100, size = d))
    angles = [np.ones(d)*0.+noise, np.ones(d)*0.+noise]
    return angles
theta.interpolate(initialize)
class SNESProblem:
    def __init__(self, F, u, bcs, J=None):
        V = u.function_space
        du = ufl.TrialFunction(V)
        self.L = fem.form(F)
        if J is None:
            self.a = fem.form(ufl.derivative(F, u, du))
        else:
            self.a = fem.form(J)
        self.bcs = bcs
        self._F, self._J = None, None
        self.u = u
    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        with F.localForm() as f_local:
            f_local.set(0.0)
        fem.petsc.assemble_vector(F, self.L)
        fem.petsc.apply_lifting(F, [self.a], bcs=[self.bcs], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(F, self.bcs, x, -1.0)
    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        fem.petsc.assemble_matrix(J, self.a, bcs=self.bcs)
        J.assemble()
# --- Define lower bound and upper bound for theta
theta_min = fem.Function(V_theta)
theta_max = fem.Function(V_theta)
def eval_theta_lowerbound(x):
    d = x.shape[1]
    return [-np.pi*np.ones(d), -np.pi/2*np.ones(d)]
def eval_theta_upperbound(x):
    d = x.shape[1]
    return [np.pi*np.ones(d), np.pi/2*np.ones(d)]
theta_min.interpolate(eval_theta_lowerbound)
theta_max.interpolate(eval_theta_upperbound)


# convert m to Cartesian coordinates
def m(theta):
    return ufl.as_vector([ufl.cos(theta[0])*ufl.cos(theta[1]),
                          ufl.sin(theta[0])*ufl.cos(theta[1]),
                          ufl.sin(theta[1])])
# Exchange energy
W_ex = ufl.tr(ufl.grad(m(theta))*ufl.grad(m(theta)).T)
# Uniaxal anisotropy energy
anisotropy_axis = (1.0, 0.0, 0.0)
a_u = fem.Constant(domain,PETSc.ScalarType(anisotropy_axis))
W_ua = -ufl.dot(a_u, m(theta))**2
# Total energy
W = (W_ex + W_ua)*dv
W_residual_theta = ufl.derivative(W, theta, theta_test)
W_jacobian_theta = ufl.derivative(W_residual_theta, theta, theta_trial)
# --- Allocate the memory needed for TAO
dofs_domain, dofs_borders = V_theta.dofmap.index_map, V_theta.dofmap.index_map_bs
b = la.create_petsc_vector(dofs_domain, dofs_borders)
J = fem.petsc.create_matrix(fem.form(W_jacobian_theta))
# --- Create optimization problem
problem = SNESProblem(W_residual_theta, theta, [], J= W_jacobian_theta)
solver_snes = PETSc.SNES().create()
solver_snes.setType("vinewtonrsls")
solver_snes.setFunction(problem.F, b)
solver_snes.setJacobian(problem.J, J)
solver_snes.setTolerances(rtol=1.0e-9, max_it=50)
solver_snes.getKSP().setType("preonly")
solver_snes.getKSP().setTolerances(rtol=1.0e-9)
solver_snes.getKSP().getPC().setType("lu")
# We set the bound (Note: they are passed as reference and not as values)
solver_snes.setVariableBounds(theta_min.vector,theta_max.vector)
solver_snes.solve(None, theta.vector)
print("Done") (edited)