how to install all dependencies at a time because when one thing working on the same time the second module not found came.
I am working on “scordelis_lo.py”
This post does not contain enough information for us to give any guidance.
For instance
this is not a file from FEniCS as far as I am aware.
You also haven’t listed what dependencies you are talking about.
file is not uploaded here that’s why i am unable to send here.
One cannot upload files to the forum. You can add code, by encapsulating it with 3x`,
i.e.
```python
#add code here
```
tell me the solution of that
Scordelis-Lo Roof Shell Benchmark
import os
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import Index
You need to be more clear what your problem actually is.
The first set of your packages are standard python packages. The second set of your packages come from legacy FEniCS.
I have already stated how you should format your code in:
# Partly Clamped Hyperbolic Paraboloid
# ====================================
#
import os
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import Index
parameters["form_compiler"]["quadrature_degree"] = 4
output_dir = "output/"
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# We consider an hyperbolic paraboloid shell made of a linear elastic isotropic
# homogeneous material with Young modulus ``Y`` and Poisson ratio ``nu``; ``mu``
# and ``lb`` denote the shear modulus :math:`\mu = Y/(2 (1 + \nu))` and the Lamé
# constant :math:`\lambda = 2 \mu \nu/(1 - 2 \nu)`. The (uniform) shell
# thickness is denoted by ``t``. :
Y = 1.
t_ = 0.1
gravity = 1.
factor = 1.
t = Constant(t_)
nu = 0.3
mu = Y/(2.0*(1.0 + nu))
lb = 2.0*mu*nu/(1.0 - 2.0*nu)
export = True
# The midplane of the initial (stress-free) configuration :math:`\mathcal{S}_I`
# of the shell is given in the form of an analytical expression
#
# .. math:: y_I:x\in\mathcal{M}\subset R^2\to y_I(x)\in\mathcal{S}_I\subset \mathcal R^3 \qquad \mathcal{S}_I = \{x_1, x_2, x_1^2-x_2^2\}
#
# in terms of the curvilinear coordinates :math:`x`. Hence, we mesh the
# two-dimensional domain :math:`\mathcal{M}\equiv [-L/2,L/2]\times [-L/2,L/2]`.
L = 1.
P1, P2 = Point(-L/2, -L/2), Point(L/2, L/2)
ndiv = 80
mesh = RectangleMesh(P1,P2, ndiv, ndiv)
# We provide the analytical expression of the initial shape as an
# ``Expression`` that we represent on a suitable ``FunctionSpace`` (here
# :math:`P_2`, but other are choices are possible)::
initial_shape = Expression(('x[0]','x[1]','0*(x[0]*x[0] - x[1]*x[1])'), degree = 4)
V_y = FunctionSpace(mesh, VectorElement('P', triangle, degree=2, dim=3))
yI = project(initial_shape, V_y)
# We compute the covariant and contravariant component of the
# metric tensor, ``aI``, ``aI_contra``; the covariant base vectors ``g0``, ``g1``;
# the contravariant base vectors ``g0_c``, ``g1_c``. ::
aI = grad(yI).T*grad(yI)
aI_contra, jI = inv(aI), det(aI)
g0, g1 = yI.dx(0), yI.dx(1)
g0_c, g1_c = aI_contra[0,0]*g0 + aI_contra[0,1]*g1, aI_contra[1,0]*g0 + aI_contra[1,1]*g1
# Given the midplane, we define the corresponding unit normal as below and
# project on a suitable function space (here :math:`P_1` but other choices
# are possible)::
def normal(y):
n = cross(y.dx(0), y.dx(1))
return n/sqrt(inner(n,n))/factor
V_normal = FunctionSpace(mesh, VectorElement("P", triangle, degree = 1, dim = 3))
nI = project(normal(yI), V_normal)
# We can visualize the shell shape and its normal with this
# utility function ::
def plot_shell(y,n=None):
y_0, y_1, y_2 = y.split(deepcopy=True)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(y_0.compute_vertex_values(),
y_1.compute_vertex_values(),
y_2.compute_vertex_values(),
triangles=y.function_space().mesh().cells(),
linewidth=1, antialiased=True, shade = False)
if n:
n_0, n_1, n_2 = n.split(deepcopy=True)
ax.quiver(y_0.compute_vertex_values(),
y_1.compute_vertex_values(),
y_2.compute_vertex_values(),
n_0.compute_vertex_values(),
n_1.compute_vertex_values(),
n_2.compute_vertex_values(),
length = .2, color = "r")
ax.view_init(elev=50, azim=30)
ax.set_xlim(-0.5, 0.5)
ax.set_ylim(-0.5, 0.5)
ax.set_zlim(-0.5, 0.5)
ax.set_xlabel(r"$x_0$")
ax.set_ylabel(r"$x_1$")
ax.set_zlabel(r"$x_2$")
return ax
plot_shell(yI)
plt.savefig('output/initial_configuration.png')
# In our 5-parameter Naghdi shell model the configuration of the shell is
# assigned by
#
# - the 3-component vector field ``u_`` representing the (small) displacement
# with respect to the initial configuration ``yI``
#
# - the 2-component vector field ``theta_`` representing the (small) rotation
# of fibers orthogonal to the middle surface.
#
# Following [2, 3], we use a ``P2+bubble`` element for ``y_`` and a ``P2``
# element for ``theta_``, and collect them in the state vector
# ``z_=[u_, theta_]``. We further define ``Function``, ``TestFunction``, and
# ``TrialFucntion`` and their different split views, which are useful for
# expressing the variational formulation. ::
P2 = FiniteElement("P", triangle, degree = 2)
bubble = FiniteElement("B", triangle, degree = 3)
Z = FunctionSpace(mesh, MixedElement(3*[P2 + bubble] + 2*[P2]))
z_ = Function(Z)
z, zt = TrialFunction(Z), TestFunction(Z)
u0_, u1_, u2_, th0_, th1_ = split(z_)
u0t, u1t, u2t, th0t, th1t = split(zt)
u0, u1, u2, th0, th1 = split(z)
# We define the displacement vector and the rotation vector, with this latter
# tangent to the middle surface, :math:`\theta = \theta_\sigma \, g^\sigma`,
# :math:`\sigma = 0, 1` ::
u_, u, ut = as_vector([u0_, u1_, u2_]), as_vector([u0, u1, u2]), as_vector([u0t, u1t, u2t])
theta_, theta, thetat = th0_*g0_c + th1_*g1_c, th0*g0_c + th1*g1_c, th0t*g0_c + th1t*g1_c
# The extensional, ``e_naghdi``, bending, ``k_naghdi``, and shearing, ``g_naghdi``, strains
# in the linear Naghdi model are defined by ::
e_naghdi = lambda v: 0.5*(grad(yI).T*grad(v) + grad(v).T*grad(yI))
k_naghdi = lambda v, t: -0.5*(grad(yI).T*grad(t) + grad(t).T*grad(yI)) - 0.5*(grad(nI).T*grad(v) + grad(v).T*grad(nI))
g_naghdi = lambda v, t: grad(yI).T*t + grad(v).T*nI
# Using curvilinear coordinates, the constitutive equations are written in terms
# of the matrix ``A_hooke`` below, representing the contravariant components of
# the constitutive tensor for isotropic elasticity in plane stress, see *e.g.*
# [4]. We use the index notation offered by ``UFL`` to express operations
# between tensors ::
i, j, k, l = Index(), Index(), Index(), Index()
A_hooke = as_tensor((((2.0*lb*mu)/(lb + 2.0*mu))*aI_contra[i,j]*aI_contra[k,l] + 1.0*mu*(aI_contra[i,k]*aI_contra[j,l] + aI_contra[i,l]*aI_contra[j,k])),[i,j,k,l])
# The membrane stress and bending moment tensors, ``N`` and ``M``, and shear
# stress vector, ``T``, are ::
N = as_tensor((A_hooke[i,j,k,l]*e_naghdi(u_)[k,l]),[i, j])
M = as_tensor(((1./12.0)*A_hooke[i,j,k,l]*k_naghdi(u_, theta_)[k,l]),[i, j])
T = as_tensor(((5./6.)*mu*aI_contra[i,j]*g_naghdi(u_, theta_)[j]), [i])
# Hence, the contributions to the elastic energy densities due to membrane,
# ``psi_m``, bending, ``psi_b``, and shear, ``psi_s``, are ::
psi_m = .5*inner(N, e_naghdi(u_))
psi_b = .5*inner(M, k_naghdi(u_, theta_))
psi_s = .5*inner(T, g_naghdi(u_, theta_))
memb_strain = inner(e_naghdi(u_),e_naghdi(u_))
bend_strain = inner(k_naghdi(u_, theta_),k_naghdi(u_, theta_))
shear_strain = inner(g_naghdi(u_, theta_),g_naghdi(u_, theta_))
# Shear and membrane locking are treated using the PSRI proposed by Arnold and
# Brezzi, see [2, 3]. In this approach, shear and membrane energies are splitted
# as a sum of two weighted contributions, one of which is computed with a reduced
# integration. Thus, shear and membrane energies have the form
#
# .. math:: (i=m,s) \qquad \alpha \int_\mathcal{M} \psi_i \, \sqrt{j_I} \mathrm{d}x + (\kappa - \alpha) \int_\mathcal{M} \psi_i \, \sqrt{j_I} \mathrm{d}x_h, \qquad \text{where} \, \kappa \propto t^{-2}
#
# While [2, 3] suggest a 1-point reduced integration, we observed that this leads
# to spurious modes in the present case.
# We use then :math:`2\times 2`-points Gauss integration for the
# portion :math:`\kappa - \alpha` of the energy, whilst the rest is integrated with a :math:`4\times 4` scheme.
# As suggested in [3], we adopt an optimized weighting factor :math:`\alpha = 1` ::
dx_h = dx(metadata={'quadrature_degree': 2})
alpha = 1.0
kappa = 1.0/t**2
shear_energy = alpha*psi_s*sqrt(jI)*dx + (kappa - alpha)*psi_s*sqrt(jI)*dx_h
membrane_energy = alpha*psi_m*sqrt(jI)*dx + (kappa - alpha)*psi_m*sqrt(jI)*dx_h
bending_energy = psi_b*sqrt(jI)*dx
# Then, the elastic energy is ::
elastic_energy = (t**3)*(bending_energy + membrane_energy + shear_energy)
# The shell is subjected to a constant vertical load. Thus, the external work is ::
f = Constant(gravity)
external_work = f*t*u2_*sqrt(jI)*dx
# We now compute the total potential energy with its first and second derivatives ::
Pi_total = elastic_energy - external_work
residual = derivative(Pi_total, z_, zt)
hessian = derivative(residual, z_, z)
# The boundary conditions prescribe a full clamping on the :math:`x_0 = -L/2` boundary,
# while the other sides are left free ::
left_boundary = lambda x, on_boundary: abs(x[0] + L/2) <= DOLFIN_EPS and on_boundary
right_boundary = lambda x, on_boundary: abs(x[0] - L/2) <= DOLFIN_EPS and on_boundary
axial_boundary = lambda x, on_boundary: abs(x[0] + L/2) <= DOLFIN_EPS or abs(x[0] - L/2) <= DOLFIN_EPS and on_boundary
full_boundary = lambda x, on_boundary: abs(x[0] + L/2) <= DOLFIN_EPS or abs(x[1] + L/2) <= DOLFIN_EPS or abs(x[0] - L/2) <= DOLFIN_EPS or abs(x[1] - L/2) <= DOLFIN_EPS and on_boundary
clamp = DirichletBC(Z, project(Expression(("0.0", "0.0", "0.0", "0.0", "0.0"), degree = 1), Z), left_boundary)
bcs = [clamp]
# We now solve the linear system of equations ::
A, b = assemble_system(hessian, residual, bcs=bcs)
solve(A, z_.vector(), b, 'lu')
u0_h, u1_h, u2_h, th0_h, th1_h = z_.split(deepcopy=True)
# Finally, we can plot the final configuration of the shell ::
scale_factor = 1.
plot_shell(project(scale_factor*u_ + yI, V_y))
plt.savefig("output/linear_naghdi/finalconfiguration.png")
XDMFFile("output/linear_naghdi/z_.xdmf").write(u1_h)
file_phi = File(output_dir + "linear_naghdi/configuration.pvd")
phi = project(scale_factor*u_+ yI, V_y)
phi.rename("phi", "phi")
file_phi << (phi)
file_energy = File(output_dir + "linear_naghdi/energy.pvd")
en_function = project((t**3)*((alpha*psi_m*sqrt(jI) + (kappa - alpha)*psi_m*sqrt(jI)) + psi_b*sqrt(jI) + (alpha*psi_s*sqrt(jI) + (kappa - alpha)*psi_s*sqrt(jI))), FunctionSpace(mesh, 'Lagrange', 1))
work_function = project(f*t**3*u2_*sqrt(jI), FunctionSpace(mesh, 'Lagrange', 1))
ben_share_function = project(psi_b*sqrt(jI)/(psi_s*sqrt(jI) + psi_m*sqrt(jI) + psi_b*sqrt(jI)), FunctionSpace(mesh, 'Lagrange', 1))
memb_en_function = project((t**3)*(psi_m*sqrt(jI)), FunctionSpace(mesh, 'Lagrange', 1))
bend_en_function = project((t**3)*(psi_b*sqrt(jI)), FunctionSpace(mesh, 'Lagrange', 1))
shear_en_function = project((t**3)*(psi_s*sqrt(jI)), FunctionSpace(mesh, 'Lagrange', 1))
en_function.rename("Elastic Energy", "Elastic Energy")
file_energy << (en_function)
memb_strain_function = project(memb_strain, FunctionSpace(mesh, 'Lagrange', 1))
bend_strain_function = project(bend_strain, FunctionSpace(mesh, 'Lagrange', 1))
shear_strain_function = project(shear_strain, FunctionSpace(mesh, 'Lagrange', 1))
center = (0.0, 0.0)
left = (0.5, -0.5)
right = (0.5, 0.5)
test = (0.25,0.25)
# output
if export:
# for L2-error
grid_samples = 100
x = (np.linspace(0, 1, grid_samples)-0.5)*L
y = (np.linspace(0, 1, grid_samples)-0.5)*L
xv, yv = np.meshgrid(x, y)
path = 'sobol_8192.csv'
coords = np.genfromtxt(path, delimiter=',')-0.5
u = np.zeros(len(coords))
u1 = np.zeros(len(coords))
u2 = np.zeros(len(coords))
th1 = np.zeros(len(coords))
th2 = np.zeros(len(coords))
for i in np.arange(0, len(coords)):
u[i] = u0_h(coords[i])
u1[i] = u1_h(coords[i])
u2[i] = u2_h(coords[i])
th1[i] = th0_h(coords[i])
th2[i] = th1_h(coords[i])
pred = np.stack([coords[:,0],coords[:,1],u,u1,u2,th1,th2],axis=1)
np.savetxt("fenics_pred.csv", pred, delimiter=",", header="xi_1,xi_2,u_x,u_y,u_z,th_1,th_2", comments="")
print("Displacement at (0.5,0) in z: ", u2_h((0.5, 0.0)))
print("E: {:.2e}, g: {:.2e}, t:{:.2e}".format(Y,gravity,t_))
print("Membrane/Bending/Shear energy contribution: {:.2f}/{:.2f}/{:.2f}".format(assemble(t**3*membrane_energy)/assemble(elastic_energy), assemble(t**3*bending_energy)/assemble(elastic_energy), assemble(t**3*shear_energy)/assemble(elastic_energy)))
print("Membrane energy: {:.2e}".format(assemble(t**3*membrane_energy)))
print("Bending energy: {:.2e}".format(assemble(t**3*bending_energy)))
print("Shear energy: {:.2e}".format(assemble(t**3*shear_energy)))
print("Inner energy: {:.2e}".format(assemble(elastic_energy)))
print("External work: {:.2e}".format(-assemble(external_work)))
Now you have just added a code, with no context. What modules do you not have installed?
What error messages are you getting?
What system are you using?
I am install fenics by following command: at site https://fenicsproject.org/download/
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt update
sudo apt install fenicsx
after installation see the following commands:
vm236@vm236:~ cat /etc/debian_version
bookworm/sid
vm236@vm236:~ uname -a
Linux vm236 5.15.0-122-generic #132-Ubuntu SMP Thu Aug 29 13:45:52 UTC 2024 x86_64 x86_64 x86_64 GNU/Linux
vm236@vm236:~$ python3 -c “import dolfinx; print(dolfinx.version); print(dolfinx.path)”
Authorization required, but no authorization protocol specified
0.8.0
[‘/usr/lib/petsc/lib/python3/dist-packages/dolfinx’]
You haven’t let us know the secret of what your error is. That makes it a little difficult to help you.
when I run the above code:
i get the following Errors:
ModuleNotFoundError Traceback (most recent call last)
/tmp/ipykernel_31910/3164710742.py in
12 import matplotlib.pyplot as plt
13
—> 14 from dolfin import *
15 from mshr import *
16 from ufl import Index
ModuleNotFoundError: No module named ‘dolfin’
That must be the problem. fenicsx installs modern FEniCS (FEniCS-X), dolfinx not dolfin.
dolfin is legacy FEniCS
sudo apt install fenics
1 Like
after that following error came
ModuleNotFoundError Traceback (most recent call last)
/tmp/ipykernel_31910/2007352019.py in
13
14 from dolfinx import *
—> 15 from mshr import *
16 from ufl import Index
17
ModuleNotFoundError: No module named ‘mshr’
Are you sure you installed fenics?
dpkg -l fenics
Note that dolfinx does not use mshr. You need to make up your mind if you’re using dolfin or dolfinx.
Hello,
Installed all dependencies of dolfinx can be critical, we spend a lot of time about this. Now we use a guix cook doing it and ensuring repeatability of software configuration. This cook is available here : https://olivier.bonnefon.pages.mia.inra.fr/mse/index.html .
regards,