# How to apply component wise dirichlet boundary condition on a point

Hello Everyone, I am solving a plane strain periodic boundary condition problem for linear elastic analysis. In order to constraint the geometry, I am thinking of applying dirichlet on a point on the model.

Here, I want to constraint the point in the present example such it is constrained in x-direction, however, it is free to displace in y-direction. Such that during any applied load in y-direction, it is free to displace in y-direction, however, since x-direction is constrained which is at the x=0.05, my model is constraint enough for solution to run.

My question is how do I apply dirichlet boundary condition on a point in just x-direction.

Here, is my MWE for the problem.

Here is my mesh

``````cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {15};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
``````

Here is my code

``````import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import mesh, fem, io, common, default_scalar_type
import dolfinx.fem.petsc
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx.plot import vtk_mesh
from dolfinx.io import XDMFFile
from ufl import (FacetNormal, FiniteElement, as_matrix, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

# Mesh
mesh_path = "WPright-medium.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)

import ufl

gdim = 2

def strain(u, repr ="vectorial"):
if repr =="vectorial":
return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
elif repr =="tensorial":
return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
sigv = dot(C, strain(u))
if repr =="vectorial":
return sigv
elif repr =="tensorial":
return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))

#Cyclic Symmetry Condition
def periodic_relation(x):
out_x = np.zeros(x.shape)
out_x[0] = x[1]*(-0.05/0.1)
out_x[1] = x[1]
return out_x

#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

T = fem.Constant(domain, 1000000.0)

n = FacetNormal(domain)

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)

L_form = dot(T*n,u_) * ds(3)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()

ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)

# bc at point
def pointt(x):
return np.logical_and(np.isclose(x[0], 0.05), np.isclose(x[1], 0))
middof = fem.locate_dofs_geometrical(V, pointt)
u_center = fem.Function(V)
u_center.x.array[:] = 0

bc_center = fem.dirichletbc(u_center, middof)

bcs = [bc_center]

with common.Timer("~~Periodic: Compute mpc condition"):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V.sub(0), facet_markers, 2, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()

WPsolve = LinearProblem(a_form, L_form, mpc, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = WPsolve.solve()
uh.name = "u"

``````

I will illustrate this on a unit square (with a non-zero condition to show that it is applied).
In the following example I have constrained the x component at (0.5,0.5) to be 0.3:

``````from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (2,)))

def single_dof(x):
return np.isclose(x[0], 0.5) & np.isclose(x[1], 0.5)

constrained_vertex = dolfinx.mesh.locate_entities(mesh, 0, single_dof)
constrained_dof = dolfinx.fem.locate_dofs_topological(V.sub(0), 0, constrained_vertex)

u_bc = dolfinx.fem.Constant(mesh, 0.3)
bc_x = dolfinx.fem.dirichletbc(u_bc, constrained_dof, V.sub(0))

uh = dolfinx.fem.Function(V)
uh.x.array[:] = 0.0
dolfinx.fem.set_bc(uh.x.array, [bc_x])

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(uh)
``````
1 Like
``````dolfinx.fem.set_bc(uh.x.array, [bc_x])
``````

Thanks Jorgen for your response. Then, how do I call the bcs for my solving of the problem with mpc as shown in here

``````WPsolve = LinearProblem(a_form, L_form, mpc, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
``````

The whole point is that `bc_x` is a standard DirichletBC object, that you can send in to the linear problem as any other bc, simply define `bcs=[bc_x]` before creating `WPsolve`.

Thank you Jorgen for your prompt response. I was able to implement your suggestion, however, I am not able to constraint the model as I earlier thought for my analysis. As you know that I am solving a plane strain linear elastic problem with Periodic Boundary Condition and as per your suggestion, I could constraint the vertex in x-direction and in y-direction (together or separately). However, my results show large deflections when I constaint the vertex only in the x-direction.

I tried the same with COMSOL and it was successfull in constraining the model.
I show some images for problem definition and results from COMSOL and FENICSX.

I am attaching here some images to explain you my problem

By constraining only in â€śxâ€ť direction- your model can move arbitrarily in the y-direction, which I guess is what you are seeing (inspect the x and y components of `u` in Paraview rather than just looking at the magnitudes).

Hi Jorgen, I just plotted the displacement magnitudes for comparison. I have seen contour plots of u in x and y direction and for the 3rd case where only x is constrained for the vertex, the u displacements are 1e7 order of magnitude as if its free. However, if we see the results from COMSOL, the displacements are within of the O(1e-7). What other option I can try here? Any suggestionsâ€¦

Could you report how these figures above look like if you only consider the x-component and then only the y-component separately in the post-processing?
As I have already stated, if you only constrain in x-direction, you can translate in y direction without violating the PDE.

Dear Jorgen, please find attached results for x and y translation invidually for only x and y constrains respectively.

As you can clearly see by these pictures, it is not sufficient to only constrain one of the directions.

But then, how do you suggest I constrain the system without impacting the solution. As here, if I constraint the top and bottom boundary, it affects the solution. Also, as presented earlier, the results from COMSOL showed that even if with x-directional constraint, the system did not result in large displacements.

You can add nullspace handling to the solver, along the lines of:

Saying that COMSOL gave you a finite solution doesnâ€™t really help much, as COMSOL is a commercial package, which most people: