It's a bit late, but I have a very simple answer: there is nothing wrong with the code
I was simply missing the extra optimisation that Newman recommends in his paper. The results are perfectly in-line with what is reported in the paper without applying the extra optimisation step. I'll leave this up in case it helps someone implementing the algorithm.
You are approximating a definite integral with cumtrapz it won't give you the same result as the integrated equation unless you add a constant and plot with the given x coordinates:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as it
x = np.arange(-10,10, 0.01) # start,stop,step
f = x**2
As, $J2 = -21 meV$ is more dominant than $J1 = 2.3meV$, the system is antiferromagnetic in nature. The expected ground state energy per moelcule(NiO) is $-42meV$. When the state reaches equilibirum, two of the nearest neighbors are alike and two are opposite, cancelling the energy contirbutions of each other. The energy contribution is from the second ...
I think that the answer by @helloworld922 is misleading. The first image shown in the answer seems to be an effect of the Poisson effect, a contraction in one direction due to loads applied in the other direction.
If you want to obtain a state of constant stress in your simulation you need to change the boundary conditions that you are applying, namely:
Here is how you want to test this and you need only two elements in the mesh. You want to define your left BC so it will reproduce a constant stress state as follows: assuming $u$ is the displacement in the x-direction and $v$ the displacement in the y-direction, set $u=0$ at the two nodes on this edge and $v=0$ at the bottom node on this edge. The two nodes ...
The primary problem is that the CST approximation has a different displacement response depending on the orientation of mesh elements relative to the applied element loading (you're only allowed to applied forces on the nodes of triangles, so distributed loads must be approximated).
You can see the effect of this by looking at only a single triangle ...
I experienced the same issue with assumed-shape dummy arrays, and the approach by Jonatan Ostrom (here) worked for me also. Because f2py creates a fresh variable for dummy arguments with intent(out), an explicit shape information seems necessary for such arguments. This is in contrast to dummy arrays with intent(in) and intent(inout) or with no intent, for ...
CVXPY's norm atom won't accept a raw Python list as an argument; you need to pass it a CVXPY expression. Stack the list of scalars into a vector using the hstack atom, like so:
constraints = [cp.norm(
y_hat[col] - cp.trace(
np.transpose((B_hat_star[:,col][:,np.newaxis]*np.sqrt(L)*C_hat[col,:])) @ X)
for col in ...