Newton and script dependant Neumann Boundary condition

dolfinx
,
1

Hi, i would like to implement a script dependant Neumann BC on Fenicsx 0.7.2.
For information the non linear system to solve is the same as :

Let’s say i have a python script that take values of x and values of u(x) as input from the Neumann BC, and that give me the corresponding value of g : script_python(x_Neum_BC,u_Neum_BC) = g

how can i implement this g, while preserving the ufl formalism and the newton solver that go with it.
Exemple of how i think it should/could look like :

# exemple algorithm
    g = script_python ( retrieve_value_of_x_at_Neum_BC , retrieve_value_of_u_at_Neum_BC)   

#g is now a list or numpy array

    F+= inner(    rebuild ufl form of g with values calculated from script     ,  test_fct ) * ds   

    self.problem = NonlinearProblem(F,self.uh, bcs=self.bcs)

i thought of doing smth like this, but i don’t know how to get x values and u values at the boundary in a numpy array, when the solution is not yet calculated.
Also i don’t know if the Newton solver will call the script for each iterations to have a new “g” each time, as i want it to do?

is there any way to have a Scripted Neumann BC with a Newton Solver in Fenicsx?
Thanks fo reading

2

If your input u is defined everywhere, then just interpolate it in the FE space and use the interpolated field in ufl.

If your input u is only defined at the boundary, then you may want to restrict interpolation to the cells which are attached to the boundary. See for instance Interpolation and IO — DOLFINx 0.7.2 documentation for an example on how to restrict interpolation to a specific set of cells.

3

Thanks for answering Francesco,
i will implement the interpolation you spoke of this week, and see if it does work as i want.
i’ll come back here then

4

Hello,
So I tried my best to implement what I needed, and after some thought I don’t think I need interpolation.
• I first need to take numerical values of uh at the Neumann boundary
• send them with corresponding numerical x values to my script
• then change numerical values in uh at the boundary with the new calculated ones
• do the Neumann BC fct

For testing purpose the script using numerical value of uh is actually doing the same calculation as the “g” I use before : “g = - constant * u² *sgn(u)”, it’s just that now it’s a function using numerical values of uh.
I manage to get, with this Neumann(script(uh)) BC, the same residual as the 1st prgrm I made for the first Newton iteration. But alas it’s always this same residual for 50 iterations. I don’t know, it’s like uh is not updated in Newton.

Part of my program :
Add neumann to F :

g2d =  self.fct_for_uh()
F += self.set_neumann_bc(self.markers_1d["top_mat2d"],g2d,self.v)

Self.set_neumann_bc() is the same as my previous program, you can see it in my previous post by following the link.

Function for recalculating uh at BC and recreating uh ufl

def fct_for_uh(self):
    tdim = self.mesh_geom.topology.dim-2# got the Dofs
    cells1 = locate_entities_boundary(self.mesh_geom, tdim, lambda m: m[1] >= 59.9)  ## is exactly the indexes where x[1]>=59.9
    coords = self.V.tabulate_dof_coordinates()
    xx = coords[cells1]# coordinate of these dofs
    uu = self.uh.eval(xx,cells1)# value of u at these dofs

    uu = fct_dens_chrg(xx[0],uu) ## the new script, using uh values at BC to calculate g values
    self.uh.vector[cells1]=uu

    return self.uh

Function to calculate g, using numerical values of uh at BC :

import numpy as np
def fct_dens_chrg(x,u):
    u=-13.295*pow(u,2)*np.sign(u)   # x will be used later
    return u

self.problem = NonlinearProblem(F,self.uh, bcs=self.bcs) 
or
self.problem = NonlinearProblem(F, self.fct_for_uh(), bcs=self.bcs)

above, either one gives the same results : 1st iterations seems ok, but residual never change

so i think i’m close to the solution, i just need to make the newton solver actually use my script for uh, and not just for the 1st iterations. if anyone has an idea??
Thanks for reading

5

We can’t help anymore unless we have a complete code, i.e. a code we can copy and paste and run. For the sake of increasing the chances of people volunteering their time to read it, please make it as minimal (but still complete) as possible.

6

ok this is the code for the mesh, just launch it to create the file.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 22 15:28:10 2024

"""

import sys

try:
    import gmsh
except ImportError:
    print("This demo requires gmsh to be installed")
    sys.exit(0)

from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI

# WARNING: unit length is 1.0 nm

c=50## coté carré
esp=20## espacement carré avec 2dmat
ld=600## long domaine
hd=300## haut domaine

mesh_res_r_min = 1
mesh_res_r_max = 20

gmsh.initialize()
gmsh.model.add("syst")

# Geometry

#### le carré à 20nm du 2dmat
p1 = gmsh.model.occ.addPoint(-c/2, esp, 0.0, mesh_res_r_min)
p2 = gmsh.model.occ.addPoint(c/2, esp, 0.0, mesh_res_r_min)
p3 = gmsh.model.occ.addPoint(c/2, esp+c, 0.0, mesh_res_r_min)
p4 = gmsh.model.occ.addPoint(-c/2, esp+c, 0.0, mesh_res_r_min)

#### le domaine entier 600 x 300 nm
p5 = gmsh.model.occ.addPoint(-ld/2, 0.0, 0.0, mesh_res_r_min)
p6 = gmsh.model.occ.addPoint(+ld/2, 0.0, 0.0, mesh_res_r_min)
p7 = gmsh.model.occ.addPoint(+ld/2, hd, 0.0, mesh_res_r_max)
p8 = gmsh.model.occ.addPoint(-ld/2, hd, 0.0, mesh_res_r_max)

### carré
line1 = gmsh.model.occ.addLine(p1, p2)
line2 = gmsh.model.occ.addLine(p2, p3)
line3 = gmsh.model.occ.addLine(p3, p4)
line4 = gmsh.model.occ.addLine(p4, p1)
### domaine
line5 = gmsh.model.occ.addLine(p5, p6)
line6 = gmsh.model.occ.addLine(p6, p7)
line7 = gmsh.model.occ.addLine(p7, p8)
line8 = gmsh.model.occ.addLine(p8, p5)

carre = gmsh.model.occ.addCurveLoop([line1, line2, line3, line4])
surf_carre        = gmsh.model.occ.addPlaneSurface([carre])

domaine = gmsh.model.occ.addCurveLoop([line5, line6, line7, line8])
surf_dom          = gmsh.model.occ.addPlaneSurface([domaine])

surf_poisson = gmsh.model.occ.fragment([(2, surf_dom)], [(2, surf_carre)], removeTool=False)
# #####
gmsh.model.occ.synchronize()

bottom_line = [gmsh.model.getBoundary(surf_poisson[1][0])[0]] ## après multiple modif

gmsh.model.mesh.setSize(bottom_line, mesh_res_r_min)   ###  bonne syntax mais ne change pas la resol du mesh
####  /!\   /!\   ####   ne marchait pas car il fallait ajouter que les points au coin du bas avait eux aussi cette resolution
##############################################################################
boundary_dirichlet   = gmsh.model.addPhysicalGroup(1, [line1, line2,line3,line4])
boundary_neumann_air_haut = gmsh.model.addPhysicalGroup(1, [line7])
boundary_neumann_air_gauche = gmsh.model.addPhysicalGroup(1, [line8])
boundary_neumann_air_droite = gmsh.model.addPhysicalGroup(1, [line6])
boundary_neumann_2dmat = gmsh.model.addPhysicalGroup(1, [line5])

gmsh.model.setPhysicalName(1, boundary_neumann_2dmat, "boundary_2dmat")
##############################################################################

###  ajouter Permittivité relative au matériaux   
dielectric_poisson = gmsh.model.addPhysicalGroup(2,surf_poisson[0][0]) #####  le [0][0] correspond à la surface 2-1 dom-carré,  [1][0] c'est les lignes je crois
gmsh.model.setPhysicalName(2, dielectric_poisson, "dielectric_poisson")
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate()

gmsh.write("mesh2d_liu_test.msh")

and this is the program, i can’t really make it smaller, most important part is from line 200 to 249 : " Solver Poisson".
line 243 is what i want to make work. g = script (value of u at BC)
in this mesh the Neumann BC is at y=0

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu May 23 15:18:00 2024

"""
####     environnement fenew
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import dolfinx.fem.petsc
import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import meshio
from petsc4py import PETSc
###### ajouté depuis https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import vtk_mesh
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
                 div, dot, dx, grad, inner, lhs, rhs)
import pyvista
from dolfinx import mesh, fem, io, nls, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import CellType, create_rectangle, locate_entities, locate_entities_boundary

##############################################################################

def fct_dens_chrg(x,u):
    u=-13.295*pow(u,2)*np.sign(u)
    return u

##############################################################################
class PoissonSolverBase:
  """
  Class that defines basic quantities useful for all
  Poisson solvers
  """
  def __init__(self):
    self.proc_id = MPI.COMM_WORLD.rank
    self.markers_1d = {}
    self.markers_2d = {}
    self.markers_3d = {}
    self.bcs = []

  def __create_meshio(self, mesh, cell_type, prune_z=False):
    """
    Select parts of an existing meshio mesh and create an existing meshio mesh
    """
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

  def reads_gmsh_mesh(self, name_of_file):
    """
    Reads gmsh mesh
    """
    if self.proc_id == 0:
      # Read in mesh
      self.mesh_meshio = meshio.read(name_of_file)

  def convert_to_xdmf(self, name_out_file, element_type, prune_z):
    """
    Convert gmsh mesh to xdmf
    """
    if self.proc_id == 0:
      # Create and save one file for the mesh, and one file for the facets
      mesh_tmp =self. __create_meshio(self.mesh_meshio, element_type, prune_z=prune_z)
      meshio.write(name_out_file, mesh_tmp)

  def create_mesh(self, name_of_file):
    """
    Read the H5-XDMF files and generate the mesh
    """
    with XDMFFile(MPI.COMM_WORLD, name_of_file, "r") as xdmf:
      self.mesh_geom       = xdmf.read_mesh(name="Grid")
      self.mesh_geom_tags  = xdmf.read_meshtags(self.mesh_geom, name="Grid")
    self.mesh_geom.topology.create_connectivity(self.mesh_geom.topology.dim, self.mesh_geom.topology.dim-1)

  def create_mesh_boundaries(self, name_of_file):
    """
    Read the meshtags for boundaries
    """
    with XDMFFile(MPI.COMM_WORLD, name_of_file, "r") as xdmf:
      self.mesh_boundaries_tags = xdmf.read_meshtags(self.mesh_geom, name="Grid")

  def set_dielectric_constant(self, eps_val_markers):
    """
    Set the dielectric constant function
    """
    self.dielectric_cons = fem.Function(self.Q)
    for values in eps_val_markers:
      cells_tmp = self.mesh_geom_tags.find(values[0])
      self.dielectric_cons.x.array[cells_tmp] = np.full_like(cells_tmp, values[1], dtype=ScalarType)

  def set_dirichlet_bc(self, marker, value):
    """
    Set Dirichlet boundary conditions from marker
    """
    facets = self.mesh_boundaries_tags.find(marker)
    dofs   = fem.locate_dofs_topological(self.V, self.mesh_geom.topology.dim-1, facets)
    self.bcs.append(fem.dirichletbc(ScalarType(value), dofs, self.V))  ##  je crois c'est ça à changer pour neumann value ce sra la value de g, et on veut g=0

########################    Neumann     BC       ###########################
   
  def set_neumann_bc(self, markers, value, test_fct):
       """
       Set Neumann boundary conditions from marker
       """
       if type(markers)==int :
           markers=[markers]
   
       facet_indices, facet_markers = [], []
       for marker in markers:
         facets = self.mesh_boundaries_tags.find(marker)
         facet_indices.append(facets)
         facet_markers.append(np.full_like(facets, marker))

       facet_indices = np.hstack(facet_indices).astype(np.int32)
       facet_markers = np.hstack(facet_markers).astype(np.int32)
       sorted_facets = np.argsort(facet_indices)
       facet_tag = mesh.meshtags(self.mesh_geom,  self.mesh_geom.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
       dss = ufl.Measure("ds", domain=self.mesh_geom, subdomain_id = facet_markers[0], subdomain_data = facet_tag)  # avec subdomain_data = facet_tag on a assertion error
       Formul = -1*inner(value, test_fct) * dss
       return Formul

##############################################################################

  def define_measure_on_boundary(self, markers):

    facet_indices, facet_markers = [], []
    for marker in markers:
      facets = self.mesh_boundaries_tags.find(marker)
      facet_indices.append(facets)
      facet_markers.append(np.full_like(facets, marker))
    facet_indices = np.hstack(facet_indices).astype(np.int32)
    facet_markers = np.hstack(facet_markers).astype(np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = mesh.meshtags(self.mesh_geom,  self.mesh_geom.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
    return ufl.Measure("dS", domain=self.mesh_geom, subdomain_data=facet_tag)

  def solve(self):
    """
    Solve the linear problem
    """
    self.uh = self.problem.solve()

  def save_solution_xdmf(self, name_of_file):
    """
    Save the solution to xdfm format
    """
    with io.XDMFFile(self.mesh_geom.comm, name_of_file, "w") as file:
      file.write_mesh(self.mesh_geom)
      file.write_function(self.uh)

  def visualize_solution(self):
    """
    Visualize the solution with pyvista
    """
    try:
      import pyvista
      cells, types, x = plot.vtk_mesh(self.V)   ### plot.create_vtk obsolete^
      grid = pyvista.UnstructuredGrid(cells, types, x)
      grid.point_data["u"] = self.uh.x.array.real
      grid.set_active_scalars("u")
      plotter = pyvista.Plotter()
      plotter.add_mesh(grid, show_edges=True)
      warped = grid.warp_by_scalar()
      plotter.add_mesh(warped)
      plotter.show()
    except ModuleNotFoundError:
      print("'pyvista' is required to visualise the solution")
      print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
###############################################################################
###############################################################################

  def fct_for_uh(self):

       tdim = self.mesh_geom.topology.dim-2# got the Dofs
       cells1 = locate_entities_boundary(self.mesh_geom, tdim, lambda m: m[1] == 0)  ## is exactly the indexes where x[1]>=59.9
       coords= self.V.tabulate_dof_coordinates()
       xx=coords[cells1]# coordinate of these dofs
       uu=self.uh.eval(xx,cells1)# value of u at these dofs
       uu=fct_dens_chrg(xx[0],uu)## pour pybinding il faudra r soit x[0]
       self.uh.vector[cells1]=uu
       
       return self.uh
   
###############################################################################
###################             Solver  Poisson              #################

class system_to_solve(PoissonSolverBase):
  """
  Class that solves the Poisson equation for Liu non linear system
  """
  def __init__(self, name_of_gmsh_file, eps, bias):
    super().__init__()

    self.reads_gmsh_mesh(name_of_gmsh_file)
    self.convert_to_xdmf("mesh_system.xdmf", "triangle", True)
    self.convert_to_xdmf("mesh_system_boundaries.xdmf", "line", True)
    self.create_mesh("mesh_system.xdmf")
    self.create_mesh_boundaries("mesh_system_boundaries.xdmf")

    # Define markers  
    self.markers_1d["square"]    = 1 # carre
    self.markers_1d["left"]  = 4 # bord gauche
    self.markers_1d["right"]  = 2 # bord droite
    self.markers_1d["top"]    = 3 # bord haut
    self.markers_1d["mat2d"]  = 5 ## bord bas ici
    self.markers_2d["air"] = 6 # le milieu air

    # Create function space
    self.Q = fem.functionspace(self.mesh_geom, ("DG", 0))  ## c'est pour .set_dielectric_constant
    self.V = fem.functionspace(self.mesh_geom, ("CG", 1))  ## 
   
    ####################      u    v  #####   pour NEWTON nonlinear :  
    self.uh = fem.Function(self.V)
    self.v = ufl.TestFunction(self.V)
    # Create a function dielectric_cons for the dielectric constant with different value in every region (one region here)
    eps_markers = [(self.markers_2d["air"], eps[0])]
    self.set_dielectric_constant(eps_markers)

##########################         F         ##############################

    F= self.dielectric_cons * inner(grad(self.uh), grad(self.v)) * dx  

    ###Apply Dirichlet boundary conditions where needed using markers
    self.set_dirichlet_bc(self.markers_1d["square"], bias)

    ####### Neumann border ext + 2dmat
    
    # g2d = -13.295*pow(self.uh,2)*ufl.sign(self.uh)##  explicit neumann BC ( WORKING)
    g2d =  self.fct_for_uh()   #  that's my objective  script dependant Neumann BC not working (not far?)
    F+=self.set_neumann_bc(self.markers_1d["mat2d"],g2d,self.v)
   
    self.problem = NonlinearProblem(F,self.uh, bcs=self.bcs)

##########################################################################

bias = 1.5  # potential at the square
eps = [1.0]  #  dielec constant, air

to_solve = system_to_solve("mesh2d_liu_test.msh", eps, bias)

######   non linear Newton
solver = NewtonSolver(MPI.COMM_WORLD, to_solve.problem)
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()

opts[f"{option_prefix}setInitialGuessNonzero"] = "True"

ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(to_solve.uh)
assert (converged)
print(f"Number of iterations: {n:d}")





###########    SAVE AND FIGURE
to_solve.save_solution_xdmf("sol_poisson_liu_system.xdmf")
to_solve.visualize_solution()

cells, types, x = plot.vtk_mesh(to_solve.V) 
UH = to_solve.uh.x.array.real
y=x[:,1]

## les x en y = 0
ind=np.where(y<1E-10)[0] # 601 values, i took 1nm step for 2dmat for x = -300 : +300nm, so looks good
## UH for y=0
UHX = UH[ind]
### we take x at y=0
X = x[ind,0]

sort=np.argsort(X) # index in increasing order for X
X=X[sort]### applied to X
UHX=UHX[sort] ### applied to UHX

plt.plot(X,UHX,'-')
plt.title('Potentiel(x) y=0')

hope it helps

7

Sorry, that’s not true. There are a lot of things you could do to make it smaller: for instance, (1) we don’t need to see your actual mesh, just use a builtin mesh, (2) we don’t need to have the part related to the visualization.

This is a forum of volunteers belonging to academia and/or companies, who reply to messages on top of everything else they have going on in their professional life.
I don’t mean to be rude, but I won’t be reading a 300 lines code. The earlier you understand what is the best way to interact with this community, the highest the probability someone will reply to your messages.

8

And, just to avoid any misunderstanding, I am not making these comments up out of my own mind: see Read before posting: How do I get my question answered?

1 Like
9

thanks for replying i understand totally, i will create a new one.

10

As I should have done sooner, i created a minimal working exemple :
it’s a square domain, with +1V on top border, and Neumann BC on bottom that is function of u.

from dolfinx import default_scalar_type, log
from dolfinx.fem import (Constant, Function, functionspace,
                          assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem,NonlinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags, locate_entities
from dolfinx.plot import vtk_mesh
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner, Measure
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import numpy as np
import ufl

def fct_dens_chrg(x,u):   #  THE g using numerical Values of uh at Neum bottom boundary
    g=-13.295*pow(u,2)*np.sign(u)     # same as before but with uh values at Neumann bc
    return g
def fct_for_uh():  # this function use fct_dens_chrg(x,u) to change values in ufl uh at Neum BC
     tdim = mesh.topology.dim-2# got the Dofs
     cells1 = locate_entities_boundary(mesh, tdim, lambda m: m[1] == 0) # bottom boundary for Neumann
     coords= V.tabulate_dof_coordinates()
     xx=coords[cells1]# coordinate of these dofs
     uu=uh.eval(xx,cells1)# values of uh at these dofs
     uu=fct_dens_chrg(xx[0],uu) 
     uh.vector[cells1]=uu
     return uh
#####################   Actual  Code    #####################################
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(mesh, ("Lagrange", 1))
uh = Function(V)
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
######  Neumann BC , g explicit uh :  Works if decommented    #########
# g = -13.295 *pow(uh,2) *ufl.sign(uh) 
####  Neumann g BC ### using values of uh at bottom border :
g =  fct_for_uh()     # stagnant residual if decommented (same residual as 1st iteration of above g)

#######   ds fabrication     #####
boundaries = [(1, lambda x: x[1]==0)]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

L = inner(g, v) * ds
###  Dirichlet BC on top +1V  #####
dofs_up = locate_dofs_geometrical(V, lambda x: np.isclose(x[1], 1))
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]

F=a-L
problem = NonlinearProblem(F ,uh, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.rtol = 1e-6
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}setInitialGuessNonzero"] = "True"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(uh)
assert (converged)
print(f"Number of iterations: {n:d}")

with g = -13.295 *pow(uh,2) *ufl.sign(uh), it works. with g=fct_for_uh() it doesn’t converge.
i think i didn’t implement a way to tell Newton solver to update uh at Neumann BC at each iterations with my function fct_for_uh(), however i don’t really see how to do that. Tried with problem = NonlinearProblem(F ,g, bcs=bcs) but same result. Thanks for your insight.

11

This function is unsafe (as you have an explicit dependency on V and uh, and you are trying to relate the dofs of your mesh to the cell indices.
What I mean by this is that:
uh.vector[cells1]=uu does not make sense, as uh is in a first order Lagrange space, therefore being located at vertices. It seems like you are mixing these.

It is unclear to me what you are trying to do in this code.
Are you trying to assign -13.292 * pow(uh, 2)*np.sign(u) to your g?
If so, have you considered using `dolfinx.fem.expression + interpolate for this, i.e. see: Projection and interpolation — FEniCS Tutorial @ Sorbonne

12

thanks for the reply,
About what i am trying to do with this code :
I have made a program with g explicitly using uh, “g= mathematical function of uh”, and it did converge using Newton solver.

What I want is to, instead of using an explicit mathematical function, use another program to compute the values of g(uh) and apply them at the boundary for each Newton iteration.

This mentionned program works with numpy array input of coordinate x on x axis and numpy array input of u. so in final u(x).
So I need to get the value of uh and x at the Neumann boundary,
hence :
tdim = mesh.topology.dim - 2 # I want the points of facets that are in the boundary, i.e all point of the mesh where x[1]=0
cells1 = locate_entities_boundary(mesh, tdim, lambda m: m[1] == 0) # bottom boundary for Neumann , i call it cells1 but it’s not a good name i imagine since tdim=0
coords= V.tabulate_dof_coordinates() ## coordinates (x,y) from V functionspace(“lagrange”,1)
xx=coords[cells1]# coordinate of these dofs, so points in this V space
uu=uh.eval(xx,cells1) # values of uh at these points , uh is also using V functionspace so there shouldn’t be problem?
the mentionned program that i want to use will then give me the numerical value of g, g(uh(x)). for now i use the same function as the mathematical g i used before, but with numpy array. Later i will have a much bigger program :
uu=fct_dens_chrg(xx[0],uu)

Finally i will use these new g values :
uh.vector[cells1]=uu
return uh
in “F+=inner(g,v)*ds” that goes inside the Newton solver.
So that at each iteration of the Newton solver the current uh_iteration should give a new g(uh_iteration) to use in “inner(g,v)*ds”.
i hope my objectiv is clearer, it’s really about using an other program to calculate values of g and using them in an iterativ way with the Newton solver.

“uh.vector[cells1]=uu does not make sense” idk for me it was, i must be very confused about what vertices is ? a facet is a triangle and a dof or vertices should be one of the 3 points of a facet? i thought that since uh=Function(V) and cells1 are vertices on the boundary following V coordinates they can be linked?, i have no error messages of at least shape or size, even on more complex mesh.

I tried to use interpolate to get values of uh at Neumann BC for calculating g(uh) but since i’ve managed to do it with eval i stop trying, i’ll look at it maybe it will help.
Thanks

13

Ok, so lets split this into multiple steps.
What you want to do is to find the value of some dolfinx.fem.Function named uh at the boundary. When you say at the boundary, I assume you mean the degrees of freedom associated with the facet.

This would give you edges, not facets.

Here you are mixing mesh entities (edges/facets) and degrees of freedom.
Consider the following code to get the degrees of freedom at the boundary

from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))

# Create a dummy functon with some data n
def f(x):
    return np.sin(x[0])*np.sin(0.1+x[1])

uh = dolfinx.fem.Function(V)
uh.interpolate(f)


# Next we define what boundary we want to locate dofs on. I will consider the same function as you used, i.e. y=0
fdim = mesh.topology.dim - 1
facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda m: m[1] == 0)

# With this, we can find the degrees of freedom associated (i.e. being defined on the facet, edges or vertices) with these facets
dofs = dolfinx.fem.locate_dofs_topological(V, fdim, facets)

# We now verify our extraction
np.set_printoptions( formatter={'all':lambda x: f"{x:.3f}"})
x_coord = V.tabulate_dof_coordinates()
for dof in dofs:
    print(f"Coordinates of dof {dof}: {x_coord[dof]}, uh: {uh.x.array[dof]:.5e} Expected{f(x_coord[dof]):.5e}")

yielding:

Coordinates of dof 0: [0.900 -0.000 0.000], uh: 7.82022e-02 Expected7.82022e-02
Coordinates of dof 1: [1.000 -0.000 0.000], uh: 8.40069e-02 Expected8.40069e-02
Coordinates of dof 5: [0.950 0.000 0.000], uh: 8.12060e-02 Expected8.12060e-02
Coordinates of dof 9: [0.800 -0.000 0.000], uh: 7.16161e-02 Expected7.16161e-02
Coordinates of dof 11: [0.850 0.000 0.000], uh: 7.50029e-02 Expected7.50029e-02
Coordinates of dof 21: [0.700 -0.000 0.000], uh: 6.43145e-02 Expected6.43145e-02
Coordinates of dof 23: [0.750 0.000 0.000], uh: 6.80503e-02 Expected6.80503e-02
Coordinates of dof 37: [0.600 -0.000 0.000], uh: 5.63702e-02 Expected5.63702e-02
Coordinates of dof 39: [0.650 0.000 0.000], uh: 6.04178e-02 Expected6.04178e-02
Coordinates of dof 57: [0.500 -0.000 0.000], uh: 4.78627e-02 Expected4.78627e-02
Coordinates of dof 59: [0.550 0.000 0.000], uh: 5.21817e-02 Expected5.21817e-02
Coordinates of dof 81: [0.400 -0.000 0.000], uh: 3.88770e-02 Expected3.88770e-02
Coordinates of dof 83: [0.450 0.000 0.000], uh: 4.34241e-02 Expected4.34241e-02
Coordinates of dof 109: [0.300 -0.000 0.000], uh: 2.95028e-02 Expected2.95028e-02
Coordinates of dof 111: [0.350 0.000 0.000], uh: 3.42327e-02 Expected3.42327e-02
Coordinates of dof 141: [0.200 -0.000 0.000], uh: 1.98338e-02 Expected1.98338e-02
Coordinates of dof 143: [0.250 0.000 0.000], uh: 2.46992e-02 Expected2.46992e-02
Coordinates of dof 177: [0.100 -0.000 0.000], uh: 9.96671e-03 Expected9.96671e-03
Coordinates of dof 179: [0.150 0.000 0.000], uh: 1.49189e-02 Expected1.49189e-02
Coordinates of dof 217: [0.000 0.000 0.000], uh: 2.15064e-19 Expected2.15064e-19
Coordinates of dof 219: [0.050 0.000 0.000], uh: 4.98959e-03 Expected4.98959e-03
1 Like
14

Thanks for the reply,
i’ve changed the part that takes the dofs at the bottom boundary.
for reminder, my objective is to have g calculated from a function that takes for input uh_neumann (values of uh in numpy array) and the corresponding x(numpy array) at Neumann boundary, and also this g must be updated for each newton iterations.

I still have the same problem of non convergence, relative residual is constant for 50 iterations, i tried lots of stuff these last days to try to make it work, without success :
First i have put g = Function(V) to initialize g instead of g=uh which was causing bugs. now relative residual goes to 8e-16 after one iterations.
After some testing, I think that my g(uh values at boundary) is not well updated with the uh found after each Newton iterations.
so I did this :
i put solver.rtol at 1e10, that way Newton always converge and do one step of iteration,
i put the calculus of g and the call of this newton solver inside a for loop with 6 iterations (the iteration number that was needed with g(uh) explicit to reach residual <= 1e-8)
the aim was to get the newly calculated uh from Newton and use it in g(uh) before re-entering Newton.
With this change, using g(uh) explicit still give me the right converged answer.
And using g(uh values at boundary) make the residual changing across iterations, but it’s diverging.
Looking at the values of uh at the boundary after one newton iteration, i can see that they are expressed by ~ -13.295 * uh(prev iterations) *sign(uh prev iterations).
so indeed now g is updated by uh values at boundary :

uh at dof[0] at y=0, for g(uh) Working EXPLICITE formula :
iteration = 0
uh bndry (0) = 0.0 ( Before entering Newton iteration, at dofs[0] of Neumann boundary)
uh after Newton bndry (0) = 0.9999999999999906 ( After Newton iteration)

iteration = 1
uh bndry (0) = 0.9999999999999906
uh after Newton bndry (0) = 0.5026490711024902

iteration = 2
uh bndry (0) = 0.5026490711024902
uh after Newton bndry (0) = 0.2612376401526086

iteration = 3
uh bndry (0) = 0.2612376401526086
uh after Newton bndry (0) = 0.15760262089435506
so this is from the working g(uh) for comparison with the non worjing one below :

uh at dof[0] at y=0, for g(uh) Not working numpy formula iteration = 0
uh bndry (0) = 0.0
uh after Newton bndry (0) = 0.9999999999999906 same as above after iteration 0

iteration = 1
uh bndry (0) = 0.9999999999999906
uh after Newton bndry (0) = -12.803907652525295 ~ -13.295 * (0.999)² *sign(0.999) diverged compared to explicit g(uh)

iteration = 2
uh bndry (0) = -12.803907652525295
uh after Newton bndry (0) = 2252.3460503325764 ~ -13.295 * (-12.8)² *sign(-12.8)

iteration = 3
uh bndry (0) = 2252.3460503325764
uh after Newton bndry (0) = -69405356.56087281 ~ -13.295 * (2252)² *sign(2252)

I think uh is internally updated inside the newton solver, and i need to get this uh updated to apply it to my function to get a right g(uh)
since i can have indices of uh corresponding to uh in my Neumann border (dofs), i thought that maybe inside the code of Newton solver
i can add after uh is updated : g.vector[dofs] = -13.295*pow(uh.x.array[dofs],2)*np.sign(uh.x.array[dofs]), that way g get the latest uh value.

i don’t know where and exactly how the newton solver is implemented, but i don’t think i can have a g, dependant of a function of values of uh, defined outside the newton solver and expect it to work fine with it.
if you have insight or idea or question, i’ll be gratefull, Thanks for reading
the updated MWE :

from dolfinx import default_scalar_type, log
from dolfinx.fem import (Constant, Function, functionspace,
                          assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem,NonlinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags, locate_entities
from dolfinx.plot import vtk_mesh
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from ufl import * #SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner, Measure
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import numpy as np
import pyvista
import ufl

def fct_dens_chrg(x,potential_vec):   #  THE g using numerical Values of u at Neum bottom boundary
    ggg=-13.295*pow(potential_vec,2)*np.sign(potential_vec)
    return ggg
def fct_for_uh(gg,potential,xx):  # use fct_dens_chrg(x,u) to change values in ufl uh
      fdim = mesh.topology.dim - 1
      facets = locate_entities_boundary(mesh, fdim, lambda m: m[1] <= 0.01)
      # With this, we can find the degrees of freedom associated (i.e. being defined on the facet, edges or vertices) with these facets
      dofs = locate_dofs_topological(V, fdim, facets)
      x_coord = V.tabulate_dof_coordinates()
      xx=x_coord[dofs]
      gg.vector[dofs] = fct_dens_chrg(xx[0],potential.x.array[dofs])
      return gg
#####################   Actual  Code    #####################################
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(mesh, ("Lagrange", 2))
uh = Function(V)
g = Function(V) 
# g = uh  # does not diverge as g = Function(V) does, instead uh_n+1 - uh_n = 0 at the boundary across all Newtons iterations 
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
#######   ds fabrication     #####
boundaries = [(1, lambda x: x[1]<=0.001)]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
###  Dirichlet BC on top +1V  #####
dofs_up = locate_dofs_geometrical(V, lambda x: x[1]>=0.99)
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]
# DOFs on bottom boundary, Dokken# With this, we can find the degrees of freedom associated (i.e. being defined on the facet, edges or vertices) with these facets
fdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, fdim, lambda m: m[1] <=0.001)
dofs = locate_dofs_topological(V, fdim, facets)
x_coord = V.tabulate_dof_coordinates() 

###########  Newton iterations (one at a time)
for iter in range(6): # let's do 6 Newton iterations
    bndry = uh.x.array[dofs]
    print(' ',' ','iteration = ',iter),print('uh before Newt at bndry = ',uh.x.array[dofs[0]]),print(' ')
    
    ####  Working Neumann g BC, explicit g(uh) :       ################################
    # g = -13.295 *pow(uh,2) *ufl.sign(uh)   # Works if decommented, g is explicit function of uh
    ####  Neumann g BC ### using values of uh at bottom border :
    # g =  fct_for_uh(g,uh,x)     
    ##### Exact same as above, but  without functions called ###########
    g.vector[dofs] = -13.295*pow(uh.x.array[dofs],2)*np.sign(uh.x.array[dofs]) 
    # print('g at boundary = ',g.x.array[dofs]),print(' ')
    ###########################################################################
    
    L = inner(g,v) * ds
    F=a-L
    problem = NonlinearProblem(F ,uh, bcs=bcs) 
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.rtol = 1e10      ### allowing one step at a time with rtol=1e10
    n, converged = solver.solve(uh)
    solver.report = True
    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}setInitialGuessNonzero"] = "True"
    ksp.setFromOptions()
    log.set_log_level(log.LogLevel.INFO)
    n, converged = solver.solve(uh)
    assert (converged)
    # print(f"Number of iterations: {n:d}"),print(' ')
    bndry_old = bndry 
    bndry = uh.x.array[dofs]
    print('difference uh new - old at boundary:', bndry - bndry_old),print(' ')
    print('uh after Newt at bndry (0) = ',bndry),print(' ')
15

Then you should either:

  1. Use a manual implementation of the Newton iteration where you perform this update step, see: Custom Newton solvers — FEniCSx tutorial
  2. Use @Andrash’s external operator code: GitHub - a-latyshev/dolfinx-external-operator: Extension of DOLFINx implementing the concept of external operator which is well-documented at: External operators within FEniCSx/DOLFINx — External operators within FEniCSx/DOLFINx
1 Like
16

Thank you very much for the reply,
I think number 2 proposition is what i was looking for. And ultimately i can create a custom Newton solver. Have a good day

17

Hello,
i advanced towards my goal, but i’m facing a difference in the solution found betwenn g(UFL) and g(Numpy).

I changed my FEniCSx version, I’m using FEniCSx 0.8.0,

My system is a square in 2d planar, where top border has uh = 1 (Dirichlet BC) and bottom border is Neumann BC with g = -uh². Using g in ufl format works fine, but i need to use a g from a numpy array.

I used the links provided by Dokken, thanks again they were usefull, to have a custom Newton solver and to use the external operator package.
With this 2 combined i was finally able to have convergence when i use g from numpy!
Unfortunately i have a non negligeable difference in the results after convergence, also the convergence rate is not the same see fig below :
diff_convergence
i think there is something that i didn’t implement, but i really can’t see what! i would like advice on that, g (UFL) is working fine but how is it different from what i did with g (numpy)?
the MWE with g numpy and commmented g ufl:

Thanks for reading

"""
Created on Tue Jun 18 13:02:24 2024
In 2d planar square system, solve laplacian(uh) = 0 in Omega, uh = 1 in d_omega_D (top border), g = -uh² in d_omega_N (Neumann bottom border)
g use external operator package
"""
import dolfinx
from dolfinx import default_scalar_type, log, fem
from dolfinx.fem import (Constant, Function, functionspace,petsc,
                          assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological)
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags, locate_entities
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
import ufl
from ufl import * #SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner, Measure
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import numpy as np
import pyvista
import matplotlib.pyplot as plt
from dolfinx_external_operator import (
    FEMExternalOperator,
    evaluate_external_operators,
    evaluate_operands,
    replace_external_operators)
import basix

def g_impl(uuh):
    output = -1*pow(uuh,2)     # NUMPY
    return output.reshape(-1)# The output must be returned flattened to one dimension  

def dgdu_impl(uuh):
    aa = -2*uuh             # NUMPY
    return aa.reshape(-1) # The output must be returned flattened to one dimension

def g_external(derivatives):
    if derivatives == (0,):  # no derivation, the function itself
        return g_impl
    elif derivatives == (1,):  # the derivative with respect to the operand `uh`
        return dgdu_impl
    else:
        return NotImplementedError
#####################   Code    #######################################
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(mesh, ("CG", 1))
uh = Function(V)
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
#####################  ds fabrication for Neumann BC, bottom of square ####
boundaries = [(1, lambda x: x[1]<=0.001)]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
#####################  Dirichlet BC on top +1V  #####
dofs_up = locate_dofs_geometrical(V, lambda x: x[1]>=0.999)
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+ 1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]
#####   dofs for NEUMANN BC
fdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, fdim, lambda m: m[1] <=0.001)
dofs = locate_dofs_topological(V, fdim, facets)

#######      g(uh ufl)  :
# g = -1*pow(uh,2)  ###   Full UFL formulation for comparaison purpose
##########################################################
#######      g with external operator package :
quadrature_degree = 2
Qe = basix.ufl.quadrature_element(mesh.topology.cell_name(), degree=quadrature_degree, value_shape=())
Q = fem.functionspace(mesh, Qe)
dx = Measure("dx", metadata={"quadrature_scheme": "default", "quadrature_degree": quadrature_degree})
g = FEMExternalOperator(uh, function_space=Q)  # g is now Symbolic not numpy involved
g.external_function = g_external
##########################################################
du = dolfinx.fem.Function(V)
maxiter = 10
i=0
correct_list=[]
L1 = inner(g,v) * ds  # g symbolic used,   this should be Neumann BC
F=a-L1
while i < maxiter:
    #########  Below work with g external operator
    ######### Assemble Jacobian and residual
    F_replaced, F_external_operators = replace_external_operators(F)
    evaluated_operands = evaluate_operands(F_external_operators)
    _ = evaluate_external_operators(F_external_operators, evaluated_operands)
    F_compiled = fem.form(F_replaced)  # Residual
    ####################    Jacobian     ###############################
    J = derivative(F, uh)# define the derivative
    J_expanded = ufl.algorithms.expand_derivatives(J)# actual value calculations of the derivative
    J_replaced, J_external_operators = replace_external_operators(J_expanded)
    _ = evaluate_external_operators(J_external_operators, evaluated_operands)
    J_compiled = fem.form(J_replaced)  # Jacobian var of the prgrm that works, jacobian = dolfinx.fem.form(J)
    
    b_vector = fem.assemble_vector(F_compiled)
    A_matrix = fem.assemble_matrix(J_compiled)
    A = dolfinx.fem.petsc.create_matrix(J_compiled)
    L = dolfinx.fem.petsc.create_vector(F_compiled)
    solver = PETSc.KSP().create(mesh.comm)
    solver.setOperators(A)

    # # # dolfinx.fem.petsc.assemble_matrix(A, jacobian)  ## replaced by :
    dolfinx.fem.petsc.assemble_matrix(A, J_compiled, bcs=bcs)
    A.assemble()  ### important
    # # # dolfinx.fem.petsc.assemble_vector(L, residual)   ## replaced by :
    dolfinx.fem.petsc.assemble_vector(L, F_compiled)
    L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    L.scale(-1)
    
    # Compute b - J(u_D-u_(i-1))
    dolfinx.fem.petsc.apply_lifting(L, [J_compiled], [bcs], x0=[uh.vector], scale=1)
    # Set du|_bc = u_{i-1}-u_D
    dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
    L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
    # Solve linear problem
    solver.solve(L, du.vector)  ##  error if A = create matrix but not if A = assemble
    du.x.scatter_forward()
    # Update u_{i+1} = u_i + delta u_i
    uh.x.array[:] += du.x.array
    # Compute norm of update
    correction_norm = du.vector.norm(0)
    correct_list.append(correction_norm)
    
    ########   Below work with g UFL , to use comment above until while loop, uncomment g UFL and fully comment g external operator
    # L1 = inner(g,v) * ds
    # F=a-L1
    # residual = dolfinx.fem.form(F)
    # J = ufl.derivative(F, uh)
    # jacobian = dolfinx.fem.form(J)
    # A = dolfinx.fem.petsc.create_matrix(jacobian)
    # L = dolfinx.fem.petsc.create_vector(residual)
    # solver = PETSc.KSP().create(mesh.comm)
    # solver.setOperators(A)
    # dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bcs)
    # A.assemble()
    # dolfinx.fem.petsc.assemble_vector(L, residual)
    # L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    # L.scale(-1)
    # # Compute b - J(u_D-u_(i-1))
    # dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[uh.vector], scale=1)
    # # Set du|_bc = u_{i-1}-u_D
    # dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
    # L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
    # # Solve linear problem
    # solver.solve(L, du.vector)
    # du.x.scatter_forward()
    # # Update u_{i+1} = u_i + delta u_i
    # uh_old=uh.copy()
    # uh.x.array[:] += du.x.array
    # # Compute norm of update
    # correction_norm = du.vector.norm(0)
    # correct_list.append(correction_norm)
    ##########################################################################

    print(f"Iteration {i}: Correction norm {correction_norm}")
    i += 1
    if correction_norm < 1e-8:
        break
    
plt.figure(1)
plt.semilogy(range(i),correct_list,'-+')
plt.title('correction norm (iteration)')
18

I cannot get your code to run at all, as g lives on a boundary, while the quadrature space is defined on the cells.
I’ve made a PR that adds support for submeshes that would capture the boundary external operator g appropriately:

1 Like
19

Thanks Dokken for the reply,

i’m really sorry about my last post, as i forgot to precise that i’m coding in fenicsx 0.8 and not anymore in 0.7.2. I changed my python environment since i couldn’t install the external operator package with fenicsx 0.7.2. i will update the post. i also retested the code of my previous post and it did work with fenicsx 0.8.

thanks for pulling a request on the github of external operator package, i will continue for now my work with g ufl, and wait for it to be updated.
Have a nice day!

20

Dear @katkit,

I am pleased to inform you that the @dokken 's PR with external operators adapted for boundaries has been merged into the main branch.

In particular, you must be interested in the test dolfinx-external-operator/tests/test_codim_external_operator.py at main · a-latyshev/dolfinx-external-operator · GitHub where the external operator g exists in the functional space based on the external boundary of the domain.

I draw your attention that the new feature was tasted on the recent nightly built of dolfinx and it does not work for v0.8.0, so please run a container with dolfinx nightly built first and then install dolfinx-external-operator from its main branch.

Do not hesitate to share any struggles with the application of the package in your work as well as a successful implementation.

2 Likes