### How to apply Dirichlet BC to some given degree of freedom

112
views
0
9 weeks ago by
Hi everyone,
I am solving a PDE problem on the domain Omega = (0,Lx) x (0, Ly).
-div (h^2 * grad(u)) = f in Omega.
u = 0 at x = 0 or x = 1.
suppose that h is given.
We consider Dirichlet boundary condition on the left, x=0, and the right, x = 1. We set the upper bound for this solution, say u <= ub = 0.2 where h = 1 and f = 2. After solving this PDE problem, there are some nodes of the mesh reaching the upper bound. Let us call the set UAS, which contains some degree of freedom of the problem.
Now I would like to solve the self-adjoint problem PDE problem i.e. consider the PDE with the same settings and  an additional Dirichlet boundary condition on the set UAS, u = 0 on UAS.
My question is that how can we apply the Dirichlet condition on some given  degree of freedom (elements of UAS)? Note that I use IndexSet to mark the UAS set.

Some lines of code that I am doing are attached. Thank you so much.

Best regards,
Nha
"""
FEniCS program:

-div(q(h,p)*grad(u)) = f   in (0, Lx) x (0, Ly)
u = uD on the boundary.

q(h,p) = h^p
i.e. p = 2 and mu =2
"""
# This programe is a modification from the code of Dr. Maurini to study a PDE problem.

# https://bitbucket.org/cmaurini/varfrac_for_cism
#==========================================================
#==========================================================
# Licensed under the GNU LGPL Version 3.
#
#==========================================================
#                      Import libraries
#==========================================================
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import numpy as np
from numpy import linalg as LA
from petsc4py import PETSc
try:
from petsc4py import PETSc
except:
print ("\n*** Warning: you need to have petsc4py installed for this module to run\n")
print ("\nExiting.\n")
exit()

if not has_petsc4py():
print ("\n*** Warning: Dolfin is not compiled with petsc4py support\n")
print ("\nExiting.\n")
exit()

#==========================================================
#==========================================================
#                   The main function
#==========================================================
#==========================================================
def main():
parameters["linear_algebra_backend"] = "PETSc"
#==========================================================
#                        Parameters
#==========================================================

# Some constants
mu = 2.; p = 2.; tol = 1e-14; zero = Constant(0.0);

ubVal = 0.2; lbVal = -100;

# Create mesh and define function space
Lx = 1.; Ly = 1.; Nx = 4; Ny = 4;
# The mesh
mesh = RectangleMesh(Point(0,0), Point(Lx, Ly), Nx, Ny, "right/left")

# Degree of the varialbes
u_degree = 2; grad_u_degree = 1; h_degree = 1;
# Define functional spaces
Vh = FunctionSpace(mesh, 'P', h_degree)
Vu = FunctionSpace(mesh, 'P', u_degree)

# Define Dirichlet boundary conditions Rec
def BoundaryX0(x, on_boundary):
return on_boundary and near(x[0], 0., tol)

def BoundaryX1(x, on_boundary):
return on_boundary and near(x[0], 1., tol)
# Apply B.C to Vu
bcX0 = DirichletBC(Vu, zero, BoundaryX0)
bcX1 = DirichletBC(Vu, zero, BoundaryX1)
bcs = [bcX0, bcX1]

# Define variational problem
u = TrialFunction(Vu)
usol = Function(Vu)
v = TestFunction(Vu)
h = Function(Vh)

# The RHS of the PDE
fVal = 2.
f = Expression('fVal', degree=1, fVal = fVal)

# Initial guess for h
h0Val = 1.
h0 = Expression('h0Val', degree=1, h0Val = h0Val)
h.assign(interpolate(h0, Vh))

L = f*v*dx

# Define the upper and lower bounds
ub = interpolate(Constant(ubVal), Vu)
lb = interpolate(Constant(lbVal), Vu)
#-----------------------------
# Solution with TAO (linear system formulation, working only for quadratic energies)
#-----------------------------
# Assemble the linear system
A, b = assemble_system(a, L, bcs)
# redefine the functional space for the displacement field, u
u = Function(Vu)

solver = TAOLinearBoundSolver("gpcg","gmres")
# Set some parameters
solver.parameters["monitor_convergence"]=True
solver.parameters["report"]=True

# Solve the problem
solver.solve(A, u.vector(), b , lb.vector(), ub.vector())
info(solver.parameters,True)

##==========================================================
##                  Find the Active set: Begin
##==========================================================
u_array = u.vector().get_local()
ub_array = ub.vector().get_local()

UAS_array = np.where(u_array[:] >(ub_array[:] - 1e-5))
UAS = PETSc.IS().createGeneral(UAS_array)
UAS.view()

#    ##==========================================================
#    ##                  Find the Active set: End
#    ##==========================================================

# The maximum of u
u_max = np.abs(u.vector().get_local()).max()
print('\n u_max = \n', u_max)

h_max = np.abs(h.vector().get_local()).max()
print('\n h_max = \n', h_max)

# plot the solution
plot(u, title= "Displacement")
plot(mesh)
vtkfile_u = File('Output/displacement.pvd')
vtkfile_u << u

if __name__ == '__main__':
main()
​
Community: FEniCS Project

3
9 weeks ago by
There may be a smarter way to do this through DOLFIN, but one solution is to just apply the BCs yourself with petsc4py.  Given the desired homogeneous Dirichlet BC DoFs as an index set, you would modify the LHS like
as_backend_type(A).mat().zeroRowsColumns(bc_dofs)
as_backend_type(A).mat().assemblyBegin()
as_backend_type(A).mat().assemblyEnd()​
where A is the LHS matrix and bc_dofs is an IS.  (Note that, despite the name of zeroRowsColumns(), its default behavior is to set the diagonal entry to 1, as desired.)  Then you can modify the RHS like
as_backend_type(b).vec().setValues(bc_dofs,zeros(bc_dofs.getLocalSize()))
as_backend_type(b).vec().assemblyBegin()
as_backend_type(b).vec().assemblyEnd()
where b is the RHS vector, and the zeros() function is imported from numpy.
Thank you so much for your help, Kamensky. I appreciate it.
I thought about that, however, I could not able to express my thoughts in terms of Fenics's language since Fenics is still new to me.
However, I did try another approach and it does work for now. What i did is to set upper bound and lower bound at these degrees of freedom to be zero. The results turn out to be good.
Thank you so much for your help.
Best wishes,
Nha
"""
FEniCS program: Optimization with PDE constraint

Min_h {J(h) = \int_\Omega f u dx + \mu * \int_\Omega h^p dx}
s.t
-div(q(h,p)*grad(u)) = f   in the unit square.
u = uD on the boundary.

q(h,p) = h^p
i.e. p = 2 and mu =2
"""
# This programe is a modification from the code of Dr. Maurini to study a PDE constrained optimization problem.

# https://bitbucket.org/cmaurini/varfrac_for_cism
#==========================================================
#==========================================================
# Licensed under the GNU LGPL Version 3.
#
#==========================================================
#                      Import libraries
#==========================================================
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
#import sympy
import numpy as np
from numpy import linalg as LA
from petsc4py import PETSc

try:
from petsc4py import PETSc
except:
print ("\n*** Warning: you need to have petsc4py installed for this module to run\n")
print ("\nExiting.\n")
exit()

if not has_petsc4py():
print ("\n*** Warning: Dolfin is not compiled with petsc4py support\n")
print ("\nExiting.\n")
exit()

#==========================================================
#==========================================================
#                   The main function
#==========================================================
#==========================================================
def main():
#    test_TAO_linear_bound_solver()
#    test_SNES_linear_bound_solver()
parameters["linear_algebra_backend"] = "PETSc"
#==========================================================
#                        Parameters
#==========================================================

# Some constants
mu = 2.; p = 2.; tol = 1e-14; zero = Constant(0.0);

ubVal = 0.2; lbVal = -100;

# Create mesh and define function space
Lx = 1.; Ly = 0.1; Nx = 100; Ny = 10;
# The mesh
#    mesh = UnitSquareMesh(Nx, Ny)

# Degree of the varialbes
u_degree = 2; grad_u_degree = 1; h_degree = 1;
# Define functional spaces
Vh = FunctionSpace(mesh, 'P', h_degree)
Vu = FunctionSpace(mesh, 'P', u_degree)

# Define Dirichlet boundary conditions Rec
def BoundaryX0(x, on_boundary):
return on_boundary and near(x[0], 0., tol)

def BoundaryX1(x, on_boundary):
return on_boundary and near(x[0], 1., tol)
# Apply B.C to Vu
bcX0 = DirichletBC(Vu, zero, BoundaryX0)
bcX1 = DirichletBC(Vu, zero, BoundaryX1)
bcs = [bcX0, bcX1]

#    # Define Dirichlet boundary conditions Unit Square
#    def boundary(x, on_boundary):
#        return on_boundary
#    bcs = DirichletBC(Vu, zero, boundary)

# Define variational problem
u = TrialFunction(Vu)
usol = Function(Vu)
v = TestFunction(Vu)
h = Function(Vh)

# The RHS of the PDE
fVal = 2.
f = Expression('fVal', degree=1, fVal = fVal)

# Initial guess for h
h0Val = 1.
h0 = Expression('h0Val', degree=1, h0Val = h0Val)
h.assign(interpolate(h0, Vh))

L = f*v*dx

# Define the upper and lower bounds
ub = interpolate(Constant(ubVal), Vu)
lb = interpolate(Constant(lbVal), Vu)
#-----------------------------
# Solution with TAO (linear system formulation, working only for quadratic energies)
#-----------------------------
# Assemble the linear system
A, b = assemble_system(a, L, bcs)
# redefine the functional space for the displacement field, u
u = Function(Vu)

solver = TAOLinearBoundSolver("gpcg","gmres")
# Set some parameters
solver.parameters["monitor_convergence"]=True
solver.parameters["report"]=True

# Solve the problem
solver.solve(A, u.vector(), b , lb.vector(), ub.vector())
#    info(solver.parameters,True)

# The maximum of u
u_max = np.abs(u.vector().get_local()).max()
print('\n u_max = \n', u_max)

h_max = np.abs(h.vector().get_local()).max()
print('\n h_max = \n', h_max)

# plot the solution
plot(u, title= "Displacement")
plot(mesh)
vtkfile_u = File('Output/displacement.pvd')
vtkfile_u << u

#==========================================================
#    Find the Active set and solve the adjoint equation
#==========================================================
u_array = u.vector().get_local()
ub_array = ub.vector().get_local()
# Find the active set associated to the upper bound
ActiveSet = list(u_array[:] >(ub_array[:] - 1e-5))

pinfty = np.inf; ninfty = -np.inf;

# Define the upper and lower bounds for adjoint

# Get upper and lower bound arrays

#Set upper and lower bound to be zero on the ActiveSet
for atDof, isActivePoint in enumerate(ActiveSet):
if isActivePoint == 1 :
#            print(atDof, isActivePoint) #List the active points if needed

#Assign values
# Print the upper and lower arrays to see

# Assemble the linear system

# Define the functional space for the displacement field of the adjoint solution

solver = TAOLinearBoundSolver("gpcg","gmres")
# Set some parameters
solver.parameters["monitor_convergence"]=True
solver.parameters["report"]=True

# Solve the problem

h_max = np.abs(h.vector().get_local()).max()
print('\n h_max = \n', h_max)

#Print the solution of the adjoint equation
# plot the solution
​