Equivalent of set_local() in Dolfinx

dolfinx
1

Hi All,

I am trying to port this fenics notebook to Dolfinx.
I somehow navigated through the syntax but I am unable to “translate” a legacy some_function.vector().set_local()) command. Searching I found the setValueLocal command (e.g. here, but it does require an index.
Looking at the original fenics code, the line theta.vector().set_local(np.maximum(np.minimum(1, thetav), thetamin)) populates a Function with certain values present in another Function, dof by dof (both are defined on the same VectorSpace.
How could I achieve the same on Dolfinx?
I attach the code I put together so far, it does run indefinitely as the update_theta function is not doing what it should, the way I am trying to “translate” the set_local command might well be total nonsense.

########## DOLFINX ##############
%matplotlib notebook
from dolfinx import *
import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
from dolfinx.fem import FunctionSpace, Function
from dolfinx.fem import VectorFunctionSpace
from dolfinx import fem, io
from dolfinx import mesh
from dolfinx.mesh import CellType, create_rectangle

import ufl
# Algorithmic parameters
niternp = 20 # number of non-penalized iterations
niter = 80 # total number of iterations
pmax = 4        # maximum SIMP exponent
exponent_update_frequency = 4 # minimum number of steps between exponent update
tol_mass = 1e-4 # tolerance on mass when finding Lagrange multiplier
thetamin = 0.001 # minimum density modeling void

L = 5

domain = create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[3*L, L]]), [30,30], cell_type=CellType.triangle)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)


# Problem parameters
thetamoy = 0.4 # target average material density
E = fem.Constant(domain, ScalarType(1))
nu = fem.Constant(domain, ScalarType(0.3))
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
f = fem.Constant(domain, ScalarType([0,-1])) # vertical downwards force
body = fem.Constant(domain, ScalarType((0, 0.1)))

# Function space for density field
V0 = fem.FunctionSpace(domain, ("DG", 0))
# Function space for displacement
V2 = fem.VectorFunctionSpace(domain, ("CG", 2))
def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0,0], dtype=ScalarType)

bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V2, fdim, boundary_facets), V2)



def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], 3*L)

fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])




metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
p = fem.Constant(domain, ScalarType(1)) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = fem.Constant(domain, ScalarType(1))# Lagrange multiplier for volume constraint

thetaold = Function(V0, name="Density")
thetaold.x.set(thetamoy)
coeff = thetaold**p
theta = Function(V0)

volume  = fem.assemble_scalar(fem.form(thetaold*ufl.dx))
avg_density_0 = fem.assemble_scalar(fem.form(thetaold*ufl.dx))/volume
avg_density = 0.
def eps(v):
    return ufl.sym(ufl.grad(v))
def sigma(v):
    return coeff*(lamda*ufl.nabla_div(v)*ufl.Identity(2)+2*mu*eps(v))
#   return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)
def energy_density(u, v):
    return ufl.inner(sigma(u), eps(v))

# Inhomogeneous elastic variational problem
u_ = ufl.TestFunction(V2)
du = ufl.TrialFunction(V2)
a = ufl.inner(sigma(u_), eps(du)) * ufl.dx
L = ufl.dot(f, u_)*ds(2)
def local_project(u,V, mesh):
#    ut = ufl.TrialFunction(V)
#    v = ufl.TestFunction(V)
    uh = Function(V)
#    solve(inner(ut,v)*dx==inner(u,v)*dx,uh)
    dx = ufl.Measure("dx", domain=mesh)
    #problem_int = fem.petsc.LinearProblem(ufl.inner(ut,v)*dx,ufl.inner(u,v)*dx,
    #                                  petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    u_0 = u
    uh = problem_int.solve()
    #uh = fem.solve(ufl.inner(ut,v)*dx, ufl.inner(u,v)*dx,uh)
    return uh


def update_theta():
    theta_new = Function(V0)
    theta.x.array[:] = local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0, domain).x.array
    thetav = theta.vector.getArray()
    ### The next line is the issue  ####
    theta.x.array[:] = (np.maximum(np.minimum(1, thetav), thetamin))

    avg_density = fem.assemble_scalar(fem.form(theta*ufl.dx))/volume
    return avg_density
def update_lagrange_multiplier(avg_density):
    print("made it to Lagrange update proc")
    avg_density1 = avg_density
    # Initial bracketing of Lagrange multiplier
    if (avg_density1 < avg_density_0):
        lagmin = float(lagrange)
        while (avg_density < avg_density_0):
        #    lagrange.assign(fem.Constant(domain, ScalarType(lagrange/2)))
            #lagrange.value(fem.Constant(domain, ScalarType(lagrange/2)))
            print(f"First while loop, lagrange is {lagrange.value}")
            lagrange.value= (lagrange.value/2)
            avg_density = update_theta()
        lagmax = float(lagrange)
    elif (avg_density1 > avg_density_0):
        lagmax = float(lagrange)
        while (avg_density > avg_density_0):
        #    lagrange.assign(fem.Constant(domain, ScalarType(lagrange*2)))
            #lagrange.value(fem.Constant(domain, ScalarType(lagrange*2)))
            print(f"Second while loop, lagrange is {lagrange.value}")
            lagrange.value = lagrange.value*2
            avg_density = update_theta()
        lagmin = float(lagrange)
    else:
        lagmin = float(lagrange)
        lagmax = float(lagrange)

    # Dichotomy on Lagrange multiplier
    inddico=0
    print("made it to Lagrange dichotomy")
    while ((abs(1.-avg_density/avg_density_0)) > tol_mass):
    #    lagrange.assign(Constant((lagmax+lagmin)/2))
        #lagrange.value(fem.Constant((lagmax+lagmin)/2))
        lagrange.value= (lagmax+lagmin)/2
        avg_density = update_theta()
        inddico += 1;
        if (avg_density < avg_density_0):
            lagmin = float(lagrange)
        else:
            lagmax = float(lagrange)
    print("   Dichotomy iterations:", inddico)
def update_exponent(exponent_counter):
    exponent_counter += 1
    if (i < niternp):
        p = fem.Constant(domain, ScalarType(1))
    elif (i >= niternp):
        if i == niternp:
            print("\n Starting penalized iterations\n")
        if ((abs(compliance-old_compliance) < 0.01*compliance_history[0]) and
            (exponent_counter > exponent_update_frequency) ):
            # average gray level
            gray_level = fem.assemble_scalar((theta-thetamin)*(1.-theta)*dx)*4/volume
 #           p.assign(Constant(min(float(p)*(1+0.3**(1.+gray_level/2)), pmax)))
            p.value = (min(float(p.value)*(1+0.3**(1.+gray_level/2)), pmax))
            p.assign(fem.Constant(domain, ScalarType(1)))
           
            exponent_counter = 0
            print("   Updated SIMP exponent p = ", float(p))
    return exponent_counter
u = Function(V2, name="Displacement")
ut = ufl.TrialFunction(V0)
u_0 = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem_int = fem.petsc.LinearProblem(ufl.inner(ut,v)*dx,ufl.inner(u_0,v)*dx,
                                      petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

old_compliance = 1e30
compliance_history = []
for i in range(niter):
#    problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()

    print("Made it past uh")

    compliance = fem.assemble_scalar(fem.form(ufl.action(L, uh)) )

    compliance_history.append(compliance)
    print("Iteration {}: compliance =".format(i), compliance)

    avg_density = update_theta()

    print("Made it past update theta")

    update_lagrange_multiplier(avg_density)

    print("Made it past Lagrange update")

    exponent_counter = update_exponent(exponent_counter)

    # Update theta field and compliance
    thetaold.x.array[:]= theta.x.array
    old_compliance = compliance
2

Could you please elaborate on what the issue is with this line?

What I can see as an issue is that you mix using dolfinx.fem.Function.vector and dolfinx.fem.Function.x.

Using the vector command gives you a petsc4py vector, which you are only accessing the local entries (those owned by the process). Using .x you Get the dolfinx.la.vector object, that has both owned entities and ghosted entities.

Please choose to use either of them.

3

Thank you this was useful. I do struggle in getting the differences between objects like dolfinx.fem.Function.vector and dolfinx.fem.Function.x.

Could you please elaborate on what the issue is with this line?

Well the code is not doing what it should do just around that line so i got dubious, and in any case I was wondering if there were better methods to assign a field to a function, other than modifying its .x.array attrubute. Thanks a lot. I found what I think is the real source of error, I will post it separately

4

The fastest way to modify data of the same length is by accessing the array. However, as pointed out above, the vector.array and x.array is of different length