Constraint boundary nodes with assigning function values before variational solve

variational formulation
1

Hi FEniCS family,
I have few queries related to constraining boundary dof.
Can we give known values to w function (ex: for boundary nodes) as a constraint?

solve(a==L,w,[])

The problem is related to tapered body with non-periodic boundaries where PBC doesn’t work. For that, I wanted to give known/specific value to unknown function w at boundary, as a constraint before solving the variational formulation.

My thought was since solve command does A^{-1} b with command solve(a==L,w,[ ])
and assign all values to w, so how could it consider constraint ? Note, A is coefficient matrix from assembling form a and b is the vector from assembling form L.

Although, we can assign values to function w using w.vector[i] = constant, where i defines nodal dof on boundary.
Q2.Can you comment/ suggest for identifying boundary nodes and its corresponding dof from mesh?
Thank you in advance.

2

Q1. Can we give known values to w function (ex: for boundary nodes) as a constraint?

You are trying to assign a dirichlet boundary condition. See https://fenicsproject.org/olddocs/dolfin/latest/python/demos/poisson/demo_poisson.py.html for instance.

Q2.Can you comment/ suggest for identifying boundary nodes and its corresponding dof from mesh?

With the same notations as the demo linked above, with

dof_coordinates = V.tabulate_dof_coordinates()

you can get the coordinates of every dof. You also have a function boundary(x) that returns True or False if a node is on the boundary or not, so something like

boundary_dofs = dict()
for (i, x) in enumerate(dof_coordinates)
    if boundary(x):
        boundary_dofs[i] = x

should give you a good starting point.

1 Like
3

Thanks @francesco-ballarin for the previous response.

I am getting difficulty in applying the constraints using DirichletBC(). Please find the attached MWE where, I am not seeing any significant effect of applying it. Can you guide me the correct way of applying constraints?


from dolfin import *
import numpy as np
import meshio
mesh=UnitCubeMesh(5,5,5)

# Extracting boundary 2D mesh
dim=3
bdim = dim-1
bmesh = BoundaryMesh(mesh, "exterior")
mapping = bmesh.entity_map(bdim)
part_of_bot = MeshFunction("size_t", bmesh, bdim)
for cell in cells(bmesh):
    curr_facet_normal = Facet(mesh, mapping[cell.index()]).normal()
    if near(curr_facet_normal.x(), -1.0):  # On bot boundary-left
        part_of_bot[cell] = 1   
    elif near(curr_facet_normal.x(), 1.0): 
        part_of_bot[cell] = 2 
         
mesh_l,mesh_r = SubMesh(bmesh, part_of_bot, 1),SubMesh(bmesh, part_of_bot, 2) 

# Input functions
def eps(v):
    E1= as_vector([0,v[1].dx(1),v[1].dx(2),0,v[2].dx(2),0]) 
    return as_tensor([(E1[0],0.5*E1[5],0.5*E1[4]),(0.5*E1[5],E1[1],0.5*E1[3]),(0.5*E1[4],0.5*E1[3],E1[2])]),E1    
def sigma(v,Eps):
    E,nu=200e3,0.25
    G=E/(2*(1+nu))
    S=np.zeros((6,6))
    S[0,0], S[1,1], S[2,2]=1/E, 1/E, 1/E
    S[0,1], S[0,2]= -nu/E, -nu/E
    S[1,0], S[1,2]= -nu/E, -nu/E
    S[2,0], S[2,1]= -nu/E, -nu/E
    S[3,3], S[4,4], S[5,5]= 1/G, 1/G, 1/G    
    C=as_tensor(np.linalg.inv(S))
    s1= dot(C,eps(v)[1]+Eps)
    return as_tensor([(s1[0],s1[5],s1[4]),(s1[5],s1[1],s1[3]),(s1[4],s1[3],s1[2])]),C,s1
# FE 
Ve = VectorElement("CG", mesh.ufl_cell(), 1,dim=3) 
Re = VectorElement("R", mesh.ufl_cell(), 0,dim=3)
W = FunctionSpace(mesh, MixedElement([Ve, Re]))
V = FunctionSpace(mesh, Ve)

v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)

# Define boundary 
def boundary(x):
    return x[0] < -0.5+DOLFIN_EPS  or x[0] > 0.5-DOLFIN_EPS

# FE for boundary mesh

# Boundary mesh left
Ve_l = VectorElement("CG", mesh_l.ufl_cell(), 1,dim=3) 
Re_l = VectorElement("R", mesh_l.ufl_cell(), 0,dim=3)
W_l = FunctionSpace(mesh_l, MixedElement([Ve_l, Re_l]))
v_l,lamb_l = TestFunctions(W_l)
dvl, dlambl = TrialFunctions(W_l)
V_l = FunctionSpace(mesh_l, Ve_l)

# Boundary mesh right
Ve_r = VectorElement("CG", mesh_r.ufl_cell(), 1,dim=3) 
Re_r = VectorElement("R", mesh_r.ufl_cell(), 0,dim=3)
W_r = FunctionSpace(mesh_r, MixedElement([Ve_r, Re_r]))
V_r = FunctionSpace(mesh_r, Ve_r)
v_r,lamb_r = TestFunctions(W_r)
dvr, dlambr = TrialFunctions(W_r)
w_l,w_r=Function(W_l), Function(W_r) # Boundary function defined

#Solve w_r, for right Boundary
x,dx = SpatialCoordinate(mesh_r), Measure('dx')(domain=mesh_r)
Eps= as_vector((0,0,0,0,x[1],-x[2]))
F2 = sum([inner(sigma(dvr, Eps)[0], eps(v_r)[0])*dx])
a2r=lhs(F2)+dot(lamb_r,dvr)*dx + dot(dlambr,v_r)*dx
L2r=rhs(F2)
solve(a2r == L2r, w_r,[])

#Solve w_l, for LEFT Boundary
x,dx = SpatialCoordinate(mesh_l), Measure('dx')(domain=mesh_l)
Eps= as_vector((0,0,0,0,x[1],-x[2]))
F2 = sum([inner(sigma(dvl, Eps)[0], eps(v_l)[0])*dx])
a2l=lhs(F2)+dot(lamb_l,dvl)*dx + dot(dlambl,v_l)*dx 
L2l=rhs(F2)
solve(a2l == L2l, w_l,[])

# Constraint apply to 3D mesh
w2a=Function(W)
xx=V.tabulate_dof_coordinates()
xx2,yy2=V_l.tabulate_dof_coordinates(), V_r.tabulate_dof_coordinates()

# equating solution at boundary dof
for i in range(int(len(xx)/3)):
    X,Y=round(xx[3*i][1],6),round(xx[3*i][2],6) # XY cross sec. coordinates (3D)
    if near(xx[3*i][0],-0.5): #left boundary
        for j in range(int(len(xx2)/3)): 
            xl,yl=round(xx2[3*j][1],6),round(xx2[3*j][2],6)
            if near(xl,X) and near(yl,Y): 
                for k in range(3):  # 3 dof at node
                    w2a.vector()[3*i+k] = w_l.vector()[3*j+k]

    elif near(xx[3*i][0],0.5): #right boundary
        for j in range(int(len(yy2)/3)): 
            xr,yr=round(yy2[3*j][1],6),round(yy2[3*j][2],6)
            if near(xr,X) and near(yr,Y):
                for k in range(3):
                    w2a.vector()[3*i+k] = w_r.vector()[3*j+k]
#   Solve
x, dx=SpatialCoordinate(mesh), Measure('dx')(domain=mesh)
Eps= as_vector((0,0,0,0,x[1],-x[2]))
F2 = sum([inner(sigma(dv,Eps)[0], eps(v_)[0])*dx])
a2=lhs(F2)+dot(lamb_,dv)*dx + dot(dlamb,v_)*dx 
w2 = Function(W) 
L2=rhs(F2)
#
  1. My query is the method of using w2a (which is obtained from taking boundary solution), is correctly implemented or not. I want all dof at boundary of 3D mesh is constrained using corresponding solution obtained through w_l and w_r. What is the correct way of applying this constraint?
  2. The dof related to w2 function are not reduced when we are constraining the boundary, which I don’t know should decrease (like in perodic boundary) or remain as it is.
bc= DirichletBC(W, w2a,boundary)   [Sphereincube_1mat.msh](https://drive.google.com/file/d/1vZLgJq8u8M_BitNVPFyxMizQw_ig-8U2/view?usp=sharing)
solve(a2 == L2, w2,bc)       
Dhe= assemble(L2).get_local() # Dhe form
(w2.vector().get_local()) # V0 form  # Issue: bc not implemented
  1. Does bc= DirichletBC(W, w2a,boundary) put assigned w2a boundary for left and right boundary as constraint. The output of w2.vector() and Dhe does not align with expected result after constraint. Pleae find the attached .msh file.

I am getting difficulty for mapping of given constraints in original mesh. Any help is appreciated.

4

There’s too much code to sort through, please simplify it as much as you can, and use a builtin mesh.

You seem to use mixed dimensional capabilities in legacy DOLFIN. I don’t have much experience with them, so you may have to wait for someone else with more experience than me to comment. However, please help them by providing a more compact example.

5

@francesco-ballarin, Thanks for your response. I have added the user defined mesh. I am sorry for the long code.

In summary, my problem is constraining the domain alternative to periodic boundary condition, where
PeriodicBoundary(Subdomain) can be used to map the x coordinates on opposite faces. However, for this case, I am constraining the domain but, instead of PBC where mapping occurs to opposite x coordinates on boundary, I am giving function value to dofs of boundary nodes.
@dokken, Can you please give your comments in FEniCS legacy problem.
Asking for general guidance, can we use something like

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and x[0] < -0.5+DOLFIN_EPS ) 
    def(map):
          # Lost

inside
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=Left(w2a))
where, w2a has corresponding output of boundary solved function as:

Q2. What can we use as constrained domain with Left(Subdomain) function, so that boundary solved function, w2a can be given as constraint in solving 3D mesh?

Just adding, the Dhe (below quote) from assembling rhs vector of weak form of 3D mesh, should have got affected like in PBC constraints. Therefore, the DirichiletBC is not working for this case. Something like FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=Left(w2a)) could apply constraint.

Thank you for your help.

6

There is too much going on in this code for me to have the bandwidth to parse it.
I would start by simply applying boundary conditions to a function in the function space, i.e.

# Create mesh here

W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=Left(w2a))

bc = DirichletBc(W.sub(0), .....)

Please note that there are no Dirichlet BCs in your code atm, so it is really hard to understand what you are trying to do.

1 Like
7

Thanks @dokken for the response. The W.sub(0) was helpful and new thing for me to learn.

Can you please share some existing resource to constrain function space W with DirchiletBC? I could only find Periodic Boundary (SubDomain), but not Dirichilet(Subdomain). My aim is to constrain the dofs of function space.
I apologize for asking silly questions.