Skip to main content

Handling Dirichlet Bounding Conditions in easifem (Part 1)

· 6 min read
Vikas Sharma
Vikas Sharma
Assistant Professor, Kyoto University, Japan

Introduction​

To apply boundary condition in FEM computation, EASIFEM, provides a class called DirichletBC_.

info

DirichletBC_ is a subclass of AbstractBC.

To understand how DirichletBC works, lets consider an example of linear elasticity. Let's say we want to apply the following boundary condition.

u=U0, on Γ\mathbf{u} = \mathbf{U}_{0}, \text{ on } \Gamma

We may think that there is only one boundary condition. But in easifem this is not the case. Actually, u\mathbf{u}, has three components in 3D (and two components in 2D). Therefore, the above boundary condition is actually boundary condition for uxu_x, uyu_y, and uzu_z. So, we have three boundary condition on a given boundary Γ\Gamma.

The second point, which is quite obvious, is that every boundary condition has two things:

  • The boundary
  • The value (condition)

To define the boundary EASIFEM employs the MeshSelection class. The value can be specified in several ways as mentioned below in this section.

note

Several instances of DirichletBC can have same boundary but different condition.

Import from toml​

We can initiate an instance of DirichletBC_ by importing the information from a toml-file. To do so, we will use the method called ImportFromToml.

Let us consider the following toml file

[domain]
filename = "./mesh/square_3x3.h5"
group = "/"
totalMedium = 1

[u]
name = "u"
engine = "NATIVE_SERIAL"
fieldType = "Normal"
spaceCompo = 1
timeCompo = 1
fedofName = "space"
geofedofName = "geofedof"

[u.geofedof]
baseContinuity = "H1"
baseInterpolation = "Lagrange"
ipType = "Equidistance"
baseType = "Monomial"

[u.space]
baseContinuity = "H1"
baseInterpolation = "Hierarchical"
# ipType = "Equidistance"
# baseType = "Monomial"
order = 4
scaleForQuadOrder = 2
quadOptName = "quadOpt"

[u.space.quadOpt]
isHomogeneous = true
quadratureType = "GaussLegendre"
# order = 2

[dbc]
name = "left_bottom_fix"
idof = 1
nodalValueType = "Constant"
value = 10.0
[dbc.boundary]
isSelectionByMeshID = true
[dbc.boundary.meshID]
line = [1, 4]

We use the following code to import the boundary condition from the above toml file.

!> author: Vikas Sharma, Ph. D.
! date: 2025-10-23
! summary: Study of DirichletBC class

PROGRAM main
USE GlobalData, ONLY: DFP, I4B, LGT
USE ExceptionHandler_Class, ONLY: e, EXCEPTION_INFORMATION
USE DirichletBC_Class
USE FEDomain_Class
USE AbstractMesh_Class
USE FEDOF_Class
USE ScalarField_Class
USE Display_Method
USE ReallocateUtility

IMPLICIT NONE

!----------------------------------------------------------------------------
! Parameters
!----------------------------------------------------------------------------
CHARACTER(*), PARAMETER :: tomlFileName = "./test1.toml"

!----------------------------------------------------------------------------
! Types and variables
!----------------------------------------------------------------------------
TYPE(DirichletBC_) :: obj
TYPE(FEDomain_) :: dom
TYPE(ScalarField_) :: u
TYPE(FEDOF_) :: fedof, geofedof
INTEGER(I4B) :: tsize, iNodeOnNode, iNodeOnEdge, iNodeOnFace, nrow, ncol

!----------------------------------------------------------------------------
! Allocatables and pointers
!----------------------------------------------------------------------------
CLASS(AbstractMesh_), POINTER :: cellMesh, boundaryMesh
INTEGER(I4B), ALLOCATABLE :: nodeNum(:)
REAL(DFP), ALLOCATABLE :: nodalValue(:, :)

!----------------------------------------------------------------------------
! Setting the verbosity
! In Debug mode, there will many messages printed to the screen
! Following code suppress information-exception.
!----------------------------------------------------------------------------

CALL e%setQuietMode(EXCEPTION_INFORMATION, .TRUE.)

!----------------------------------------------------------------------------
! Domain: Initiate FEDomain and print info
!----------------------------------------------------------------------------

CALL dom%ImportFromToml(tomlName="domain", filename=tomlFileName)
! CALL dom%DisplayDomainInfo(msg="DomainInfo: ")

!----------------------------------------------------------------------------
! ScalarField and FEDOF: Initiate scalarField and fedof
!----------------------------------------------------------------------------

CALL u%ImportFromToml(tomlName="u", fedof=fedof, geofedof=geofedof, dom=dom, &
                      filename=tomlFileName)
! CALL fedof%Display(msg="FEDOF info: ")

!----------------------------------------------------------------------------
! DirichletBC: Import from toml
!----------------------------------------------------------------------------

CALL obj%ImportFromToml(filename=tomlFileName, dom=dom, tomlName="dbc")
! CALL obj%Display("DirichletBC Info: ")

!----------------------------------------------------------------------------
! DirichletBC: Get the total node numbers
!----------------------------------------------------------------------------
tsize = obj%GetTotalNodeNum(fedof=fedof)
CALL Display(tsize, "Total Node Num: ")

!----------------------------------------------------------------------------
! DirichletBC: Get the node numbers
!----------------------------------------------------------------------------

CALL Reallocate(nodeNum, tsize, isExpand=.TRUE., expandFactor=2)
CALL obj%Get(fedof=fedof, nodeNum=nodeNum, tsize=tsize, &
             iNodeOnNode=iNodeOnNode, iNodeOnEdge=iNodeOnEdge, &
             iNodeOnFace=iNodeOnFace)
CALL Display(tsize, "tsize = ")
CALL Display(nodeNum(1:tsize), "nodeNum", full=.TRUE., orient="ROW")
CALL Display(iNodeOnNode, "iNodeOnNode = ")
CALL Display(iNodeOnFace, "iNodeOnFace = ")
CALL Display(iNodeOnEdge, "iNodeOnEdge = ")

!----------------------------------------------------------------------------
! DirichletBC: Get the nodal values
!----------------------------------------------------------------------------

CALL Reallocate(nodalValue, tsize, 2, isExpand=.TRUE., expandFactor=2)
CALL obj%Get(fedof=fedof, nodeNum=nodeNum, nodalValue=nodalValue, &
             nrow=nrow, ncol=ncol)
CALL Display(nodalValue(1:nrow, 1:ncol), "nodalValue: ", &
             full=.TRUE.)
CALL Display([nrow, ncol], "[nrow, ncol]: ")

!----------------------------------------------------------------------------
! Cleanup
!----------------------------------------------------------------------------
CALL dom%DEALLOCATE()
END PROGRAM main

Learn from example​

Let's consider the following example, in which we will specify the constant boundary condition.

Click here to see the example
PROGRAM main
USE easifemBase
USE easifemClasses
IMPLICIT NONE

TYPE(DirichletBC_) :: obj
TYPE(MeshSelection_) :: boundary
TYPE(ParameterList_) :: param
TYPE(Domain_) :: dom
TYPE(HDF5File_) :: domainfile
CHARACTER(*), PARAMETER :: domainfilename = "./mesh3D.h5"
INTEGER(I4B) :: bottom = 1, top = 2, left = 3, right = 4, &
                front = 5, behind = 6, nsd
INTEGER(I4B), ALLOCATABLE :: nodeNum(:)
REAL(DFP), ALLOCATABLE :: nodalValue(:, :)

CALL FPL_Init; CALL param%Initiate()
CALL domainfile%Initiate(filename=domainfilename, mode="READ")
CALL domainfile%OPEN()
CALL dom%Initiate(domainfile, group="")

nsd = dom%GetNSD()

! We call Set SetAbstractBCParam to set the parameter for boundary condition
CALL SetAbstractBCParam(param=param, prefix=obj%GetPrefix(),  &
  & name="ZeroBC", idof=1, nodalValueType=Constant)

! We call SetMeshSelectionParam to set the parameter for boundary condition
CALL SetMeshSelectionParam(param=param, prefix=boundary%GetPrefix(),  &
  & isSelectionByMeshID=.TRUE.)

CALL boundary%Initiate(param)

CALL boundary%Add(dom=dom, dim=nsd - 1, meshID=[top])
CALL boundary%Set()

CALL obj%Initiate(param=param, boundary=boundary, dom=dom)

CALL obj%Set(constantNodalValue=0.0_DFP)

CALL obj%Get(nodeNum=nodeNum, nodalValue=nodalValue)

CALL Display(nodeNum, "nodeNum", advance="NO")
CALL Display(nodalValue, "nodalValue", advance="YES")

CALL domainfile%DEALLOCATE()
CALL dom%DEALLOCATE()
CALL param%DEALLOCATE(); CALL FPL_Finalize
END PROGRAM main

In the above code, to define the boundary condition, we follow the steps given below.

Step 1: Set the properties of the DirichletBC​

We set the properties of DirichletBC_ by using the method called SetAbstractBCParam.

CALL SetAbstractBCParam(param=param, prefix=obj%GetPrefix(),  &
  & name="ZeroBC", idof=1, nodalValueType=Constant)
note

Because we are setting constant boundary condition, we used nodalValueType=Constant.

You can learn more about this method here.

Step 2: Define a boundary​

To define a boundary we will use the MeshSelection. In the above code, we select the boundary by specifing the meshID.

CALL SetMeshSelectionParam(param=param, prefix=boundary%GetPrefix(),  &
  & isSelectionByMeshID=.TRUE.)

After setting the boundary parameter we call Initiate method.

CALL boundary%Initiate(param)

Subsequently, we call Add method to add the information of meshID.

CALL boundary%Add(dom=dom, dim=nsd - 1, meshID=[top])
CALL boundary%Set()
info

After adding the information of meshID we should call Set method, which means that we are done adding information to the boundary.

You can learn more about SetMeshSelectionParam here

Step 3: Initiate instance of DirichletBC​

After initiating the boundary, call Initiate. To initiate an instance of DirichletBC_ we need to pass the boundary, paramters, and domain.

CALL obj%Initiate(param=param, boundary=boundary, dom=dom)

Step 4: Set the boundary condition​

After initiating an instance of DirichletBC_, next step is to set the boundary condition. To do so, we will use the method Set.

While setting the value we should respect the configuration used while calling SetAbstractBCParam. For example, in the above example we configure boundary condition for nodalValueType=Constant. Therefore, we should set the constantNodalValue while calling the set method.

CALL obj%Set(constantNodalValue=0.0_DFP)

Step 5: Get the value of boundary condition​

To get the boundary condition we will use the method Get. The Get function can take two arguments nodeNum(:) and nodalValue(:,:). The nodeNum(:) is compulsory, whereas nodalValue can be optional.

CALL obj%Get(nodeNum=nodeNum, nodalValue=nodalValue)
note

On return, the size of nodeNum and SIZE(nodalValue, 1) is same.

The columns in nodalValue denotes the boundary condition at different times. You can read more about this subroutine here.

Further reading​

There is more to DirichletBC_, and you can learn about them from following pages. (Here DBC stands for DirichletBC_)

Constant user function
This example shows how to initiate DBC by using a constant UserFunction $u=\alpha$
Spatial user function
This example shows how to initiate DBC by using a space dependent UserFunction, $u=g(x, y, z)$
Toml file constant value
This example shows how to initiate DBC by reading a toml-file. Constant boundary condition, $u=\alpha$
Toml file constant function
This example shows how to initiate DBC by reading a toml-file. Constant boundary condition, $u=\alpha$
Toml file space function
This example shows how to initiate DBC by reading a toml-file. Space boundary condition, $u=g(x,y,z)$