diff --git a/demos/3D_soil_slope_limit_analysis/3D_soil_slope_limit_analysis.py b/demos/3D_soil_slope_limit_analysis/3D_soil_slope_limit_analysis.py new file mode 100644 index 0000000..0ff2e01 --- /dev/null +++ b/demos/3D_soil_slope_limit_analysis/3D_soil_slope_limit_analysis.py @@ -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)) +# %% diff --git a/dolfinx_optim/mosek_io.py b/dolfinx_optim/mosek_io.py index 7ad486c..80cc1e0 100755 --- a/dolfinx_optim/mosek_io.py +++ b/dolfinx_optim/mosek_io.py @@ -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 @@ -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 = {} @@ -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) @@ -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)) @@ -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) @@ -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) @@ -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) @@ -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() @@ -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