Skip to content

update to dolfinx 0.9.0 #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
cmaurini wants to merge 1 commit into
base: main
Choose a base branch
from
Open

update to dolfinx 0.9.0 #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 80 additions & 0 deletions demos/3D_soil_slope_limit_analysis/3D_soil_slope_limit_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
#%%
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import mesh, fem, io
from dolfinx_optim.mosek_io import MosekProblem
from dolfinx_optim.convex_function import ConvexTerm, QuadraticTerm
from dolfinx_optim.cones import SDP
from dolfinx_optim.utils import to_vect

# Define the box mesh
L, W, H = (1.2, 2.0, 1.0)
Nx, Ny, Nz = (10, 3, 20)
domain = mesh.create_box(MPI.COMM_WORLD, [(0, 0, 0), (L, W, H)], [Nx, Ny, Nz])

# Define the conic representation of the Mohr-Coulomb support function
c = fem.Constant(domain, 1.0)
phi = fem.Constant(domain, np.pi / 6.0)

class MohrCoulomb(ConvexTerm):
"""SDP implementation of Mohr-Coulomb criterion."""

def conic_repr(self, X):
Y1 = self.add_var((3,3), cone=SDP(3))
Y2 = self.add_var((3,3), cone=SDP(3))
a = (1 - ufl.sin(phi)) / (1 + ufl.sin(phi))
self.add_eq_constraint(X - to_vect(Y1) + to_vect(Y2))
self.add_eq_constraint(ufl.tr(Y2) - a * ufl.tr(Y1))
self.add_linear_term(2 * c * ufl.cos(phi) / (1 + ufl.sin(phi)) * ufl.tr(Y1))

# Set up the loading, function spaces and boundary conditions
gamma = 10.0
f = fem.Constant(domain, (0, 0, -gamma))

def border(x):
return np.isclose(x[0], L) | np.isclose(x[2], 0)

gdim = 3
deg_quad = 4


#%%

#V = fem.functionspace(domain, ("CR", 1, (gdim,)))
V = fem.functionspace(domain, ("CG", 2, (gdim,)))
bc_dofs = fem.locate_dofs_geometrical(V, border)
bcs = [fem.dirichletbc(np.zeros((gdim,)), bc_dofs, V)]

# Initiate the MosekProblem object and add the linear equality constraint
prob = MosekProblem(domain, "3D limit analysis")
u = prob.add_var(V, bc=bcs, name="Collapse mechanism")

# Add CR stabilization
#h = ufl.CellDiameter(domain)
#h_avg = (h("+") + h("-")) / 2.0
#n = ufl.FacetNormal(domain)
#jump_u = 1/h * ufl.jump(u)
#stabilization = QuadraticTerm(jump_u, 4, measure_type="ds")
#prob.add_convex_term(stabilization)

#%%
prob.add_eq_constraint(ufl.dot(f, u) * ufl.dx, b=1.0)

# Add the convex term corresponding to the support function
crit = MohrCoulomb(ufl.sym(ufl.grad(u)), 2)
prob.add_convex_term(crit)

# Solve the problem and export results to Paraview
pobj, dobj = prob.optimize()
# %%
V_plot = fem.functionspace(domain, ("CG", 2, (gdim,)))
u_plot = fem.Function(V_plot,name="u")
u_plot.interpolate(u)
with io.VTKFile(MPI.COMM_WORLD, "results.pvd", "w") as vtk:
vtk.write_function(u_plot)

# Check the solution compared with the exact solution
print("2D factor [Chen] (for phi=30°):", 6.69)
print("Computed factor:", pobj * float(gamma * H / c))
# %%
30 changes: 15 additions & 15 deletions dolfinx_optim/mosek_io.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ def _create_mass_matrix(V2, dx):
# mass matrix on V2
one = fem.Function(V2)
v = ufl.TestFunction(V2)
one.vector.set(1.0)
one.x.petsc_vec.set(1.0)
m_ufl = ufl.inner(one, v) * dx
return fem.assemble_vector(fem.form(m_ufl)).array

Expand Down Expand Up @@ -130,7 +130,7 @@ def __init__(self, domain: dolfinx.mesh.Mesh, name: str):
self.c = []
self.parameters = self._default_parameters()
self.M = mf.Model(name)
# self.vectors_dict = {}
# self.x.petsc_vecs_dict = {}
self.vectors = []
self.objectives = []
self.constraints = {}
Expand All @@ -151,10 +151,10 @@ def _default_parameters(self):
def _create_variable_vector(self, var, name, ux, lx, cone):
if cone is not None:
n = cone.dim
m = len(var.vector.array) // n
m = len(var.x.petsc_vec.array) // n
domain = mosek_cone_domain(cone)
if cone.type == "sdp":
m = len(var.vector.array) // n**2
m = len(var.x.petsc_vec.array) // n**2
domain = mf.Domain.inPSDCone(n, m)
if name is not None:
vec = self.M.variable(name, domain)
Expand Down Expand Up @@ -193,11 +193,11 @@ def _create_variable_vector(self, var, name, ux, lx, cone):
def _add_boundary_conditions(self, variable, vector, bcs):
V = variable.function_space
u_bc = fem.Function(V)
u_bc.vector.set(np.inf)
u_bc.x.petsc_vec.set(np.inf)

fem.set_bc(u_bc.vector, to_list(bcs))
dof_indices = np.where(u_bc.vector.array < np.inf)[0].astype(np.int32)
bc_values = u_bc.vector.array[dof_indices]
fem.set_bc(u_bc.x.petsc_vec, to_list(bcs))
dof_indices = np.where(u_bc.x.petsc_vec.array < np.inf)[0].astype(np.int32)
bc_values = u_bc.x.petsc_vec.array[dof_indices]

self.M.constraint(vector.pick(dof_indices), mf.Domain.equalsTo(bc_values))

Expand Down Expand Up @@ -436,7 +436,7 @@ def add_convex_term(self, conv_fun: ConvexTerm):
cone = conv_fun.cones[i]

vector = self._create_variable_vector(var, name, ux, lx, cone)
assert len(var.vector.array) == vector.getSize()
assert len(var.x.petsc_vec.array) == vector.getSize()
self.variables.append(var)
self.vectors.append(vector)

Expand Down Expand Up @@ -516,11 +516,11 @@ def _apply_linear_constraints(self, conv_fun):
if cons["bu"] is not None:
if isinstance(cons["bu"], float):
bu = fem.Function(cons["V"])
xbu = bu.vector.array
xbu = bu.x.petsc_vec.array
xbu[:] = cons["bu"]
elif cons["bu"] == 0:
bu = fem.Function(cons["V"])
xbu = bu.vector.array
xbu = bu.x.petsc_vec.array
else:
bu = fem.assemble_vector(
fem.form(ufl.inner(lamb_, cons["bu"]) * conv_fun.dx)
Expand All @@ -531,11 +531,11 @@ def _apply_linear_constraints(self, conv_fun):
if cons["bl"] is not None:
if isinstance(cons["bl"], float):
bl = fem.Function(cons["V"])
xbl = bl.vector.array
xbl = bl.x.petsc_vec.array
xbl[:] = cons["bl"]
elif cons["bl"] == 0:
bl = fem.Function(cons["V"])
xbl = bl.vector.array
xbl = bl.x.petsc_vec.array
else:
bl = fem.assemble_vector(
fem.form(ufl.inner(lamb_, cons["bl"]) * conv_fun.dx)
Expand Down Expand Up @@ -715,7 +715,7 @@ def optimize(self, sense="min", dump=False):
# Retrieve solution and save to file
for var, vec in zip(self.variables, self.vectors):
if isinstance(var, fem.Function):
var.vector.array[:] = vec.level()
var.x.petsc_vec.array[:] = vec.level()
else:
var.value = vec.level()

Expand All @@ -735,5 +735,5 @@ def get_lagrange_multiplier(self, name):
"""Retrieves Lagrange multiplier function associated with constraint `name`."""
constraint, V_cons = self.constraints[name]
lag = fem.Function(V_cons, name=name)
lag.vector.array[:] = constraint.dual()
lag.x.petsc_vec.array[:] = constraint.dual()
return lag