# AttributeError: 'numpy.float64' object has no attribute 'block\_variable'

**URL:** <https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500>\
**Category:** dolfin-adjoint\
**Created:** [October 19, 2023, 1:18am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500 "2023-10-19T01:18:58Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 19, 2023, 1:18am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/1 "2023-10-19T01:18:59Z")

</div>

Hi all,

I want my force to reach the values I want with the applied displacement. I determine these desired points in the polynomial. So I minimize my function to find the polynomial coefficients that will allow it to pass my target points in the force-displacement graph.  
However I am having an error like:  
**vector = vector.\_cpp\_object  
^^^^^^^^^^^^^^^^^^  
AttributeError: ‘float’ object has no attribute ‘\_cpp\_object’**

The error starts with minimization.  
I have divided the code into pieces and left my minimal example below. Even though I optimize only the first constant in the example, I will later ask it to optimize all constants and many polynomials will be included in my geometry.  
I can provide you more information if needed.  
Thank you!

**Geometry-Boundary-Function Definitions (No error here)**

```auto
from dolfin import *
from dolfin_adjoint import *
import gmsh
import numpy as np
import matplotlib.pyplot as plt
import meshio
import moola

a = 0.1
t = 1
h = 1
x_dir = 1.5
y_dir = h + a
scale = 1
N = 10
lc = 1/2
max_disp = -0.3
last_point_of_spline = 0.25
Traclist = []
Traclist2 = []
normal = Constant((0, 1, 0))
recu = []
J = 0

multiply = 10

x_valuesA = [0, 0.1, 0.225, last_point_of_spline]
y_valuesA = [0, multiply*0.32, multiply*0.1, multiply*0.15]

desired_response_x = [0, 0.1, 0.225, last_point_of_spline]
desired_response_y = [0, multiply*0.36, multiply*0.02, multiply*0.07]
def createGeometryAndMesh(a, t, lc):
    # Clear all models and create a new one
    gmsh.initialize()
    gmsh.clear()
    gmsh.model.add("t3")

    def geopoints(dx, dy, surf_point):
        P1 = gmsh.model.occ.addPoint(0 + dx, dy, 0, lc)
        P2 = gmsh.model.occ.addPoint(t + dx, dy, 0, lc)
        P3 = gmsh.model.occ.addPoint(dx, h + dy, 0, lc)
        P4 = gmsh.model.occ.addPoint(t + dx, h + dy, 0, lc)

        L1 = gmsh.model.occ.addLine(P1, P2)
        L2 = gmsh.model.occ.addLine(P2, P4)
        L3 = gmsh.model.occ.addLine(P4, P3)
        L4 = gmsh.model.occ.addLine(P3, P1)

        Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])
        Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 11 + surf_point)

        return Surface1

    # Entity creation operations
    SurfaceA = geopoints(0, 0, 2000)
    gmsh.model.occ.synchronize()
    # One group for each of the unit squares
    pg1 = gmsh.model.addPhysicalGroup(2, [SurfaceA])
    gmsh.model.setPhysicalName(2, pg1, "main_geometry")

    # Set mesh size at all points
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
    gmsh.model.occ.synchronize()

    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
    gmsh.model.mesh.generate()
    gmsh.write("./Result--" + str(commit) + "geoforstl" + ".geo_unrolled")
    gmsh.write("./Result--" + str(commit) + "geo.msh")
    # colors:
    gmsh.option.setNumber("Geometry.PointNumbers", 1)
    gmsh.option.setColor("Geometry.Points", 255, 165, 0)
    gmsh.option.setColor("General.Text", 255, 255, 255)
    gmsh.option.setColor("Mesh.Points", 255, 0, 0)
    rr, gg, bb, aaa = gmsh.option.getColor("Geometry.Points")
    gmsh.option.setColor("Geometry.Surfaces", rr, gg, bb, aaa)
    gmsh.onelab.set("""[
    {"type":"number",
    "name":"Parameters/Twisting angle",
    "values":[90],
    "min":0,
    "max":120,
    "step":1}]""")
    
    def checkForEvent():
        action = gmsh.onelab.getString("ONELAB/Action")
        if len(action) and action[0] == "check":
            gmsh.onelab.setString("ONELAB/Action", [""])
            createGeometryAndMesh()
            gmsh.graphics.draw()
        return True  

    # Refine the mesh
    gmsh.model.mesh.refine()
  
    gmsh.write("t3.opt");
    gmsh.finalize()

createGeometryAndMesh(a, t, lc)

# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("./Result--" + str(commit)+ "geo.msh")

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})

meshio.write("./Result--" + str(commit)+ "mesh.xdmf",
             triangle_mesh)
meshio.write("mesh.xdmf", triangle_mesh)

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True,
               "eliminate_zeros": True,
               "precompute_basis_const": True,
               "precompute_ip_const": True}

# Create mesh and define function space
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
print(mesh.topology().dim() - 1)

File("./Result--" + str(commit)+ "MSH.pvd") << mesh
File("./Result--"+ str(commit) + "/MSH.pvd") << mesh
File("./Result--"+ str(commit) + "MSH2.pvd") << mf
File("./Result--" + str(commit)+ "/MSH2" ".pvd") << mf
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)

b_c = Boundary()
b_c.mark(boundary_markers, 3)
File("./Result--" + str(commit)+ "MSH3.pvd") << boundary_markers

# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx_A = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
V = VectorFunctionSpace(mesh, "CG", 2) # Lagrange
W_tensor = TensorFunctionSpace(mesh, "CG", 2)
V_tensor = TensorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "R", 0)

def boundary_bot(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], 0, tol)

bc_bot = DirichletBC(V, [0, 0, 0], boundary_bot)

def boundary_top(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], (h), tol)

disp = Constant([0.0, max_disp, 0.0])
bc_top = DirichletBC(V, disp, boundary_top)
bcs = [bc_bot, bc_top]

# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
u0 = Function(V) # Initial displacement
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)
dss = ds(subdomain_data=boundary_subdomains)

```

**Calculation starts here**  
**Error is also here**

```auto
def kinematics(u):
    # Kinematics
    d = u.geometric_dimension() # Space dimension
    I = Identity(d) # Identity tensor
    F = variable(I + grad(u)) # Deformation gradient
    C = F.T * F # Right Cauchy-Green tensor 
    # Invariants of deformation tensors
    Ic = tr(C)
    J = det(F)
    E = (C-I)/2

    return Ic, J, F, C, I, E

def poly(x_values, y_values):
    # Fit a cubic polynomial (degree 3) to the given data points
    coeffs = np.polyfit(x_values, y_values, len(x_values) -1)
    print(coeffs)
    # print(coeffs.convert().coef)
    if len(coeffs) == 4:
        c1, c2, c3, c4= coeffs
    if len(coeffs) == 3:
        c1 = 0
        c2 = coeffs[0]
        c3 = coeffs[1]
        c4 = coeffs[2]
    if len(coeffs) == 2:
        c1 = 0
        c2 = 0
        c3 = coeffs[0]
        c4 = coeffs[1]
    
    return c1, c2, c3, c4

for i in np.linspace(0, max_disp, N):
    u.assign(u0)
    disp.assign(Constant([0.0, i, 0.0]))
    print("i", i)
    Ic, J, F, C, I, E = kinematics(u)
    # Define the custom strain energy density function
    def psi(x_values, y_values, num):
        aA, bA, cA, dA= poly(x_values, y_values) 
        aA = interpolate(Expression(("a_constant"), a_constant=aA, degree=1), W) # W = FunctionSpace(mesh, "R", 0)
        Ic, J, F, C, I, E = kinematics(u)
        psi = (aA((0,h,0))/ 4) * (E[1, 1]) **4 - (bA / 3) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) **2 - dA * E[1, 1] + 1E+4 * E[0,0]**2 + 1E+4 * E[0,1] **2
        stress = diff(psi, F)
        return psi, stress, aA, bA, cA, dA 
    
    psiA, stressA, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)   
    a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y) 
    target_y = (a_d) * (i) **3 + (b_d) * (i)** 2 + (c_d ) * (i) + d_d 
    psi = psiA 
    PiA = psiA * dx_A
    Pi = PiA
    dPiA = derivative(PiA, u, v)
    dPi = derivative(Pi, u, v)
    stress = stressA 
    solve(dPi == 0, u, bcs, 
          form_compiler_parameters=ffc_options)
    print("u", u((0, h, 0))[1])
    traction = dot(stress, normal)
    stored_traction = assemble((traction[1]) * dss(1)) # Traction in the y-direction
    Traclist.append(stored_traction)
    """
    ERROR STARTS HERE
    """
    J = (target_y - stored_traction) ** 2
    # aA, bA, cA, dA = Function(W), Function(W), Function(W),Function(W)
    # aA = interpolate(Expression(('a_constant'), a_constant=aA, degree=1), W)
    # bA = interpolate(Expression(('b_constant'), b_constant=bA, degree=1), W)
    # cA = interpolate(Expression(('c_constant'), c_constant=cA, degree=1), W)
    # dA = interpolate(Expression(('d_constant'), d_constant=dA, degree=1), W)
    control = Control(aA) # , Control(bA), Control(cA), Control(dA)

    rf = ReducedFunctional(J, control)
    problem = MoolaOptimizationProblem(rf)

    n_moola = moola.DolfinPrimalVector(aA)
    # d_moola = moola.DolfinPrimalVector(bA)
    # k_moola = moola.DolfinPrimalVector(cA)
    # l_moola = moola.DolfinPrimalVector(dA)
    m_moola = moola.DolfinPrimalVectorSet([n_moola]) # , d_moola, k_moola, l_moola
    solver = moola.BFGS(problem, m_moola, options={'jtol': 0,
                                                'gtol': 1e-9,
                                                'Hinit': "default",
                                                'maxiter': 10,
                                                'mem_lim': 10})    
    sol = solver.solve()
    m_opt = sol['control'].data
    print(m_opt)

```

---

<div class="post-metadata">

**Author:** ![spcvanschie](https://avatars.discourse-cdn.com/v4/letter/s/76d3ee/32.png) [@spcvanschie](https://fenicsproject.discourse.group/u/spcvanschie)\
**Post date:** [October 20, 2023, 4:52am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/2 "2023-10-20T04:52:31Z")

</div>

Hi Hailey,

I’m having trouble understanding what exactly is going wrong in your code. Could you directly copy the error trace that you encounter and specify in which exact line of your code it pops up?  
Right now the bar is pretty high for people to help you, since they’d need to run your code themselves. Supplying the exact error trace and line in which it happens lowers that bar.

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 20, 2023, 5:42pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/3 "2023-10-20T17:42:07Z")

</div>

Hello @spcvanschie , thank you a lot for your answer!  
I mentioned in the code where the error starts in

```auto
sol = solver.solve()

```

and also here is the full error:

```auto
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 7.678e-29 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton solver finished in 0 iterations and 0 linear solver iterations.
Traceback (most recent call last):
  File "/home/hailey/trial_opt.py", line 295, in <module>
    sol = solver.solve()
          ^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/moola/algorithms/bfgs.py", line 172, in solve
    J = objective(xk)
        ^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/optimization/moola_problem.py", line 50, in __call__
    self.latest_eval_eval.append(rf(x.data))
                                 ^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/reduced_functional.py", line 137, in __call__
    func_value = self.scale * self.functional.block_variable.checkpoint
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'numpy.float64' object has no attribute 'block_variable'

```

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 21, 2023, 4:42am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/4 "2023-10-21T04:42:39Z")

</div>

@spcvanschie please let me know if I need to add more details

---

<div class="post-metadata">

**Author:** ![spcvanschie](https://avatars.discourse-cdn.com/v4/letter/s/76d3ee/32.png) [@spcvanschie](https://fenicsproject.discourse.group/u/spcvanschie)\
**Post date:** [October 21, 2023, 6:28am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/5 "2023-10-21T06:28:01Z")

</div>

So a disclaimer here is that I do not have much experience with Dolfin-adjoint. Nonetheless, I see an issue by comparing your code with [this](https://www.dolfin-adjoint.org/en/stable/documentation/poisson-mother/poisson-mother.html) example case. Your functional `J` seems to be incorrect. Right now you assemble a component of the `traction` form, then subtract `target_y` (which is a Numpy float) and then square that difference. I don’t know what exactly your goal functional is, but in any case it is currently defined as a Numpy float64, which is probably what’s causing your error to be thrown. You should formulate your functional `J` as a Dolfin object, so that Dolfin-adjoint can actually interpret it. You could try defining `J` as a Dolfin `Expression`.

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 21, 2023, 9:56pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/6 "2023-10-21T21:56:25Z")

</div>

Thank you so much! I tried that but it did not work. @spcvanschie @nate

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 22, 2023, 3:19pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/7 "2023-10-22T15:19:50Z")

</div>

> [@spcvanschie](#):
>
> Right now you assemble a component of the `traction` form, then subtract `target_y` (which is a Numpy float) and then square that difference. I don’t know what exactly your goal functional is, but in any case it is currently defined as a Numpy float64, which is probably what’s causing your error to be thrown.

This is indeed one issue. DOLFIN-adjoint relies on operator overloading, and has not overloaded all of numpy. Thus, you should wrap

> [@haileyfromde](#):
>
> ` J = (target_y - stored_traction) ** 2`

as  
`J = (stored_traction- AdjFloat(target_y)) ** 2` to ensure that the numpy float is overloaded

Secondly, your function `u` is not dependent of `aA`, as you are evaluating `aA` at a specific point, which only has very specific support in dolfin-adjoint:

```auto
        psi = (aA(np.array([0,h,0]))/ 4) * (E[1, 1]) **4 - (bA / 3.) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) **2 - dA * E[1, 1] + 1.E+4 * E[0,0]**2 + 1.E+4 * E[0,1] **2

```

i.e. wrapping the eval in a numpy array allows annotation.

As your example is far from minimal, I’m having issues with getting further in debugging your remaining issues. Please try to simplify your example.  
I’ve for instance hit a general issue with older versions of ufl ([Sanitise unicode input by dham · Pull Request #111 · FEniCS/ufl · GitHub](https://github.com/FEniCS/ufl/pull/111/files)) where non-linear solver output is not sanitized, so I cannot visualize the computational graph to see what is not overloaded correctly.

I believe that having a point evaluation of a function as the control variable doesn’t make much sense to me, as this would just be a scalar value (i.e a constant as control variable).

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 22, 2023, 6:44pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/8 "2023-10-22T18:44:30Z")

</div>

@dokken Thank you so much!  
The reason I use dolfin\_adjoint is that I will add hyperelastic materials later and my aim is to find the constants of the material that actually behaves like a polynomial. As you said, my Control value does not depend on u, but I do not know how I can achieve my goal by finding a value that depends on u.

I am still having errors and I minimized my code below:

```auto
def kinematics(u):
    # Kinematics
    d = u.geometric_dimension() # Space dimension
    I = Identity(d) # Identity tensor
    F = variable(I + grad(u)) # Deformation gradient
    C = F.T * F # Right Cauchy-Green tensor 
    # Invariants of deformation tensors
    Ic = tr(C)
    J = det(F)
    E = (C-I)/2

    return Ic, J, F, C, I, E

def poly(x_values, y_values):
    # Fit a cubic polynomial (degree 3) to the given data points
    coeffs = np.polyfit(x_values, y_values, len(x_values) -1)
    c1, c2, c3, c4= coeffs
    return c1, c2, c3, c4

disp.assign(Constant([0.0, max_disp, 0.0]))
Ic, J, F, C, I, E = kinematics(u)
# Define the custom strain energy density function
aA, bA, cA, dA= poly(x_valuesA, y_valuesA) 
aA = project(Expression(("a_constant"), a_constant=aA, degree=1), W) # W = FunctionSpace(mesh, "R", 0)
psi = (aA(np.array([0,h,0]))/ 4) * (E[1, 1]) **4 - (bA / 3.) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) **2 - dA * E[1, 1] + 1.E+4 * E[0,0]**2 + 1.E+4 * E[0,1] **2
stress = diff(psi, F)
traction = dot(stress, normal)
stored_traction = assemble((traction[1]) * dss(1)) 
Traclist.append(stored_traction)

"""
ERROR STARTS HERE
"""
target_y = 1e-2
J = (stored_traction- AdjFloat(target_y)) ** 2
control = Control(aA)

rf = ReducedFunctional(J, control)
problem = MoolaOptimizationProblem(rf)

n_moola = moola.DolfinPrimalVector(aA)
m_moola = moola.DolfinPrimalVectorSet([n_moola]) # , d_moola, k_moola, l_moola
solver = moola.BFGS(problem, m_moola, options={'jtol': 0,
                                            'gtol': 1e-9,
                                            'Hinit': "default",
                                            'maxiter': 10,
                                            'mem_lim': 10})    
sol = solver.solve()
m_opt = sol['control'].data
print(m_opt)
   

```

Please let me know if it is not clear enough.

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 23, 2023, 4:02pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/10 "2023-10-23T16:02:21Z")

</div>

> [@dokken](#):
>
> `aA(np.array([0,h,0]))`

Hello @dokken, I changed my code as below and now I do not have the error. However the result is not correct, it is giving me the initial guess back.

```auto
def kinematics(u):
    # Kinematics
    d = u.geometric_dimension() # Space dimension
    I = Identity(d) # Identity tensor
    F = variable(I + grad(u)) # Deformation gradient
    C = F.T * F # Right Cauchy-Green tensor 
    # Invariants of deformation tensors
    Ic = tr(C)
    J = det(F)
    E = (C-I)/2

    return Ic, J, F, C, I, E

def poly(x_values, y_values):
    # Fit a cubic polynomial (degree 3) to the given data points
    coeffs = np.polyfit(x_values, y_values, len(x_values) -1)
    if len(coeffs) == 4:
        c1, c2, c3, c4= coeffs
    return c1, c2, c3, c4

for i in np.linspace(0, -0.25, 6):
    disp.assign(Constant([0.0, i, 0.0]))
    Ic, J, F, C, I, E = kinematics(u)

    def psi(x_values, y_values, num):
        aA, bA, cA, dA= poly(x_values, y_values) 
        psi = (aA / 4) * (E[1, 1]) **4 - (bA / 3) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) ** 2 - dA * E[1, 1] 
        stress = diff(psi, F)
        return psi, stress, aA, bA, cA, dA 

    psiA, stress, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)   
    a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y) 
    target_y = [0,1,1,2,3,4]

    traction = dot(stress, normal)
    stored_traction = assemble((traction[1]) * dss(1))
    Traclist.append(stored_traction)

    
J = 0  
def objective(J, i):
    for i in range(len(Traclist)):
        J += (Traclist[i] - AdjFloat(target_y[i])) ** 2  
        print(J)
    return J

J = objective(J, i)

control = Control(Constant(aA))
control2 = Control(Constant(bA)) 
control3 = Control(Constant(cA)) 
control4 = Control(Constant(dA)) 
rf = ReducedFunctional(J, [control, control2, control3, control4])
# const = UFLInequalityConstraint(dJdx - 1e-5, control)

problem = MinimizationProblem(rf) # , constraints=const
parameters_ipopt = {"acceptable_tol": 1e-10, "maximum_iterations": 10}
solver = IPOPTSolver(problem, parameters=parameters_ipopt)
f_opt = solver.solve()
print("f_opt", f_opt[1]((0,h,0)))

```

---

<div class="post-metadata">

**Author:** ![spcvanschie](https://avatars.discourse-cdn.com/v4/letter/s/76d3ee/32.png) [@spcvanschie](https://fenicsproject.discourse.group/u/spcvanschie)\
**Post date:** [October 23, 2023, 5:18pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/11 "2023-10-23T17:18:13Z")

</div>

I can see two possible issues in this block of code:

```auto
for i in np.linspace(0, -0.25, 6):
    disp.assign(Constant([0.0, i, 0.0]))
    Ic, J, F, C, I, E = kinematics(u)

    def psi(x_values, y_values, num):
        aA, bA, cA, dA= poly(x_values, y_values) 
        psi = (aA / 4) * (E[1, 1]) **4 - (bA / 3) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) ** 2 - dA * E[1, 1] 
        stress = diff(psi, F)
        return psi, stress, aA, bA, cA, dA 

    psiA, stress, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)   
    a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y) 
    target_y = [0,1,1,2,3,4]

    traction = dot(stress, normal)
    stored_traction = assemble((traction[1]) * dss(1))
    Traclist.append(stored_traction)

    
J = 0  
def objective(J, i):
    for i in range(len(Traclist)):
        J += (Traclist[i] - AdjFloat(target_y[i])) ** 2  
        print(J)
    return J

J = objective(J, i)

```

Both are Python issues and have nothing to do with Dolfin:

- All lines after `J = 0` are not indented in the same way as the preceding lines, and so are not a part of your `for` loop.
- Within your function definition `objective()` you use `i` both as an input argument and as the argument of the `for` loop.

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 23, 2023, 5:19pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/12 "2023-10-23T17:19:37Z")

</div>

Thanks! I changed that and I am still having the same error. @spcvanschie

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 23, 2023, 6:29pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/13 "2023-10-23T18:29:17Z")

</div>

The issue is that you are defining the control as a constant at the end, you should define them as constants before sending then into the variational form, i.e

```python
    def psi(x_values, y_values, num):
        aA, bA, cA, dA= poly(x_values, y_values) 
        aA = Constant(aA)
        bA = Constant(bA)
        cA = Constant(cA)
        dA = Constant(dA)
        psi = (aA / 4) * (E[1, 1]) **4 - (bA / 3) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) ** 2 - dA * E[1, 1] 
        stress = diff(psi, F)
        return psi, stress, aA, bA, cA, dA 

    psiA, stress, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)   
    a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y) 
    target_y = [0,1,1,2,3,4]

    traction = dot(stress, normal)
    stored_traction = assemble((traction[1]) * dss(1))
    Traclist.append(stored_traction)

    
J = 0  
def objective(J, i):
    for i in range(len(Traclist)):
        J += (Traclist[i] - AdjFloat(target_y[i])) ** 2  
        print(J)
    return J

J = objective(J, i)

control = Control(aA)
control2 = Control(bA) 
control3 = Control(cA) 
control4 = Control(dA) 
rf = ReducedFunctional(J, [control, control2, control3, control4])

```

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 23, 2023, 10:26pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/14 "2023-10-23T22:26:20Z")

</div>

thank you so much! @dokken But then it gives me this error:

```auto
  File "/home/hailey/Desktop/hailey/research/datageneration/trial_opt.py", line 310, in <module>
    f_opt = solver.solve()
            ^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/optimization/ipopt_solver.py", line 199, in solve
    results = self.ipopt_problem.solve(guess)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "cyipopt/cython/ipopt_wrapper.pyx", line 642, in ipopt_wrapper.Problem.solve
  File "cyipopt/cython/ipopt_wrapper.pyx", line 676, in ipopt_wrapper.objective_cb
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/reduced_functional_numpy.py", line 36, in __call__
    return self.rf. __call__ (self.set_local(m_copies, m_array))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hailey/miniconda3/envs/fenicsproject/lib/python3.11/site-packages/pyadjoint/reduced_functional.py", line 137, in __call__
    func_value = self.scale * self.functional.block_variable.checkpoint
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'Sum' object has no attribute 'block_variable'

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 24, 2023, 4:41am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/15 "2023-10-24T04:41:05Z")

</div>

Please make a reproducible example, as the code I presented you ran on my system.

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 24, 2023, 2:54pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/16 "2023-10-24T14:54:25Z")

</div>

When I change the location of J = 0 (I showed in the code below), I am having the error.  
Also, I expect my optimized coefficients to change but they still the same after the optimization.  
Thanks!  
Here I shared it @dokken

```auto
def poly(x_values, y_values):
    # Fit a cubic polynomial (degree 3) to the given data points
    coeffs = np.polyfit(x_values, y_values, len(x_values) -1)
    c1, c2, c3, c4= coeffs

    return c1, c2, c3, c4

J = 0 # it does not work if I move it here
a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y)
for i in np.linspace(0, max_disp, N):
    disp.assign(Constant([0.0, i, 0.0]))
    Ic, J, F, C, I, E = kinematics(u)
    # Define the custom strain energy density function
    def psi(x_values, y_values, num):
            aA, bA, cA, dA= poly(x_values, y_values) 
            aA = Constant(aA)
            bA = Constant(bA)
            cA = Constant(cA)
            dA = Constant(dA)
            psi = (aA / 4) * (E[1, 1]) **4 - (bA / 3) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) ** 2 - dA * E[1, 1] 
            stress = diff(psi, F)
            return psi, stress, aA, bA, cA, dA 

    psiA, stress, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)   
    target_y = (a_d) * (-i) **3 + (b_d) * (-i)** 2 + (c_d ) * (-i) + d_d 
    traction = dot(stress, normal)
    stored_traction = assemble((traction[1]) * dss(1))
    Traclist.append(stored_traction)
    target_y_list.append(target_y)

# J = 0 # it works when J=0 is here
def objective(J):
    for k in range(len(Traclist)):
        J += (Traclist[k] - AdjFloat(target_y_list[k])) ** 2  
        print(J)
    return J

J = objective(J)
control = Control(aA)
control2 = Control(bA) 
control3 = Control(cA) 
control4 = Control(dA) 
rf = ReducedFunctional(J, [control, control2, control3, control4])
problem = MinimizationProblem(rf)  
parameters_ipopt = {"acceptable_tol": 1e-10, "maximum_iterations": 10}
solver = IPOPTSolver(problem, parameters=parameters_ipopt)
f_opt = solver.solve()
print("f_opt a", f_opt[0]((0,h,0)))
print("f_opt b", f_opt[1]((0,h,0)))
print("f_opt c", f_opt[2]((0,h,0)))
print("f_opt d", f_opt[3]((0,h,0)))

```

---

<div class="post-metadata">

**Author:** ![spcvanschie](https://avatars.discourse-cdn.com/v4/letter/s/76d3ee/32.png) [@spcvanschie](https://fenicsproject.discourse.group/u/spcvanschie)\
**Post date:** [October 24, 2023, 5:52pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/17 "2023-10-24T17:52:34Z")

</div>

> When I change the location of J = 0 (I showed in the code below), I am having the error.

Well yes. That makes sense, since you define another variable with the name `J` four lines further down:

```auto
J = 0 # it does not work if I move it here
a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y)
for i in np.linspace(0, max_disp, N):
    disp.assign(Constant([0.0, i, 0.0]))
    Ic, J, F, C, I, E = kinematics(u)

```

> Also, I expect my optimized coefficients to change but they still the same after the optimization.

What initial values are you using? Did you already try using a different set of initial values to see whether those are changed in any way by the optimization process?

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 25, 2023, 12:12am UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/18 "2023-10-25T00:12:42Z")

</div>

Thanks! @spcvanschie Yes you are right, I totally missed this!  
I tried with different initial values, with the code below, I am getting the results as:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/f05fdc7c1fbcbd3b5a6a512576cbab087274e108.png)

I think it is only optimizing the last control value.

```auto
x_spline = np.linspace(0, 0.25, N, endpoint=True)  
y_spline = f_opt[0]((0,h,0)) * x_spline **3 + f_opt[1]((0,h,0)) * x_spline** 2 + f_opt[2]((0,h,0)) * x_spline + f_opt[3]((0,h,0))
y_spline_target = a_d * x_spline **3 + b_d * x_spline** 2 + c_d * x_spline + d_d
# Plot the data points and the spline
plt.plot(x_spline, y_spline, 'o-', label='Optimized Curve')
plt.plot(x_spline, y_spline_target, label='Target Curve')
aA, bA, cA, dA= poly(x_valuesA, y_valuesA) 
y_spline_initial_guess = aA * x_spline **3 + bA * x_spline** 2 + cA * x_spline + dA
plt.plot(x_spline, y_spline_initial_guess, label='Initial Guess Curve')
plt.legend()
plt.grid(True)
plt.xlabel("Displacement (mm)")
plt.ylabel("Force (N)")
plt.title("Force-Displacement Curve")

```

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 25, 2023, 1:46pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/19 "2023-10-25T13:46:00Z")

</div>

it is only changing the 1 control value. @dokken @spcvanschie

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 25, 2023, 2:03pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/20 "2023-10-25T14:03:50Z")

</div>

As stated above in: [AttributeError: 'numpy.float64' object has no attribute 'block\_variable' - #17 by spcvanschie](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/17)  
you are overwriting `J` inside your loop

---

<div class="post-metadata">

**Author:** ![haileyfromde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/haileyfromde/32/3643_2.png) [@haileyfromde](https://fenicsproject.discourse.group/u/haileyfromde)\
**Post date:** [October 25, 2023, 2:10pm UTC](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500/21 "2023-10-25T14:10:26Z")

</div>

I fixed it @dokken, now it is working but the result is just a shifted version of the initial guess. So it is only optimizing the Control(dA). Thanks!

```auto
objective_J = 0
a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y)
for i in np.linspace(0, max_disp, N):
    u.assign(u0)
    disp.assign(Constant([0.0, i, 0.0]))
    Ic, J, F, C, I, E = kinematics(u)
    def psi(x_values, y_values, num):
        aAA, bAA, cAA, dAA = poly(x_values, y_values) 
        aA = Constant(aAA) 
        bA = Constant(bAA)
        cA = Constant(cAA)
        dA = Constant(dAA)
        psi = (aA/ 4) * (E[1, 1]) **4 - (bA/ 3) * (E[1, 1])** 3 + (cA / 2) * (E[1, 1]) **2 - dA * E[1, 1]** 1
        stress = diff(psi, F)
        return psi, stress, aA, bA, cA, dA 
    
    psiA, stressA, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)    
    target_y = (a_d) * (-i) **3 + (b_d) * (-i)** 2 + (c_d ) * (-i) + d_d 
    stress = stressA 
    solve(dPi == 0, u, bcs, 
          form_compiler_parameters=ffc_options)
    print("u", u((0, h, 0))[1])
    traction = dot(stress, normal)
    stored_traction = assemble((traction[1]) * dss(1)) # Traction in the y-direction
    Traclist.append(stored_traction)
    target_y_list.append(target_y)
    
    objective_J += (stored_traction - AdjFloat(target_y)) ** 2 
    
control = [Control(aA), Control(bA), Control(cA), Control(dA)]

rf = ReducedFunctional(objective_J, control)
problem = MinimizationProblem(rf) # , constraints=const
parameters_ipopt = {"acceptable_tol": 1e-10, "maximum_iterations": 10}
solver = IPOPTSolver(problem, parameters=parameters_ipopt)
f_opt = solver.solve()
print("f_opt a", f_opt[0]((0,h,0)))
print("f_opt b", f_opt[1]((0,h,0)))
print("f_opt c", f_opt[2]((0,h,0)))
print("f_opt d", f_opt[3]((0,h,0)))

```

[Next page](https://fenicsproject.discourse.group/t/attributeerror-numpy-float64-object-has-no-attribute-block-variable/12500.md?page=2)