 |
v0.14.0 |
|
v0.14.0 |
The focus of our current tutorial is on solving Laplace-type partial differential equations (PDEs) Poisson Problem with non-homogenous Dirichlet boundary conditions. While the SCL-2: Poisson's equation (non-homogeneous BC) has designed for solving Dirichlet problems with homogeneous boundary conditions, this tutorial guides the methods to extend this to non-homogeneous cases and plot the solution.
We are also aiming to present MoFEM approaches for handling Dirichlet problems in addition to the least square approximation while establishing your necessary foundation in MoFEM as well. Prior to that, we would like to suggest understanding the operators(i.e. OpDomainLhsMatrixK and OpDomainRhsVectorF) in SCL-1: Poisson's equation (homogeneous BC) before proceeding to this tutorial.
We are demonstrating an approach to solve a simple Poisson's problem with nonhomogenous BC and help you grasp the approach better. For this problem, we denote the domain as Ω, and Γ represents the Dirichlet boundary, as illustrated in Figure 1.
Thus, we are solving a sample problem with nonhomogeneous BCs to find \(\mathbf{u}\). such that:
Domain Equation:
\[ \begin{align} -\nabla \cdot\nabla u = f(\mathbf{x}) \,\,\, in \,\,\, \Omega \end{align} \]
Boundary condition:
\[ \begin{align} u = g(\mathbf{x}) \,\,\, on \,\,\,\Gamma \end{align} \]
Where f is the source term of the body and g is denoted as boundary function. However, if g is a constant everywhere (i.e. 0), we speak of homogeneous boundary conditions but in this problem, we are considering a spatially varying function \(g(\mathbf{x})\).
In the context of solving this partial differential equation (PDE), Galerkin finite element methods with least square approximation have been executed to enforce nonhomogeneous boundary condition. We are also utilizing MoFEM libraries conforming finite element spaces with Least-Square approximation(LSA) for the boundary by minimizing the squared differences between the approximated values and the actual, similar to the following:
\[ \frac{1}{2}\int_{\Gamma} (u - g(\mathbf{x}))^2 \,d\Gamma \]
Remember, previously we derived a weak form for homogeneous BC in SCL-1: Poisson's equation (homogeneous BC). Following that approch, we can write similar weak form for our current problem for domain (Ω) and boundary(Γ) which also exists for \( \forall \,\delta u \) as follows:
\[ \begin{align} \int_{\Omega} \nabla \delta u \cdot \nabla u \,d\Omega = \int_{\Omega} \delta u \,f(\mathbf{x})\,d\Omega \quad\text{and} \end{align} \]
\[ \begin{align} \int_{\Gamma} \delta u \, \,d\Gamma = \int_{\Gamma} \delta u \,g(\mathbf{x})\,d\Gamma \end{align} \]
In the context of finite element method, the problem will be solved numerically by approximating the scalar space \( \mathbf{H^1} \) by the discrete counterparts associated with our FE mesh of the computational domain \(\Omega\). To do that, we are utilizing two key components to solve: (i) \(N_{\alpha} \) and \(N_{\beta} \) are the shape functions for the field variable and for test functions respectively, (ii) \(\bar{u}_{\alpha}\) and \(\delta\bar{u}_{\beta}\) represent for 'n' number of DoF within that element to obtain the solution 'u' as stated below:
\[ \begin{align} u &= \sum_{\alpha=0}^{n-1} \mathbf{N}_{\alpha} \bar{u}_\alpha \,\,\text{and,}\\ {\delta u} &= \sum_{\beta=0}^{n-1}\mathbf{{N}_{\beta}} \delta\bar{u}_{\beta} \end{align} \]
Further, to assemble the element stiffness matrix we are indicating \(\alpha\) for the collums and \(\beta\) for rows. A typical approximation can be understood well in tutorial FUN-2: Hierarchical approximation. So, to approximate the fields based on the weak forms in equation (5) and (6) can be written as follows:
\[ \begin{equation} \sum_{\alpha=0}^{n-1} \sum_{\beta=0}^{n-1} \delta\bar{u}_{\beta} \left[\int_{\Omega} \nabla \mathbf{N}_{\alpha} \cdot \nabla \mathbf{N}_{\beta} \,d\Omega\right] \bar{u}_{\alpha} \,\,= \,\, \sum_{\beta=0}^{n-1} \delta \bar{u}_{\beta} \left[\int_{\Omega} f(\mathbf{x})\,\, \mathbf{N}_{\beta} \,d\Omega\right] \quad ;\,\forall \,\delta \bar{u}_{\beta} \end{equation} \]
and
\[ \begin{equation} \sum_{\alpha=0}^{n-1} \sum_{\beta=0}^{n-1} \delta\bar{u}_{\beta} \left[\int_{\Gamma} \mathbf{N}_{\alpha} \,\,\mathbf{N}_{\beta} \,d\Gamma\right] \bar{u}_{\alpha} \,\,= \,\,\sum_{\beta=0}^{n-1} \delta \bar{u}_{\beta} \left[\int_{\Gamma} g(\mathbf{x})\,\, \mathbf{N}_{\beta} \,d\Gamma\right]\quad ;\forall \,\delta \bar{u}_{\beta} \end{equation} \]
For a better understanding the implementation of the aboves in MoFEM, we are explaining the forms in two parts: domain and boundary. Also, let's say, we are naming the individuals for respective UDOs as follows,
\[ \mathbf{A_k} = \int_{\Omega} \nabla \mathbf{N}_{\alpha} \cdot \nabla \mathbf{N}_{\beta} \,d\Omega, \quad \mathbf{B_f} = \int_{\Omega} f(\mathbf{x})\,\,\mathbf{N}_{\beta} \,d\Omega\,\,\ \text{ |\,For Domain}, \]
and
\[ \mathbf{C_k} = \int_{\Gamma} \mathbf{N}_{\alpha} \cdot \mathbf{N}_{\beta} \,d\Gamma, \quad \mathbf{D_f} = \int_{\Gamma} g(\mathbf{x})\,\, \mathbf{N}_{\beta} \,d\Gamma \,\,\text{ |\,For Boundary}. \]
Additionally, a similar discretization and implementation approach was referred in tutorial SCL-1: Poisson's equation (homogeneous BC).
Based on the above schemes, we are considering 4 operators in MoFEM. Let's say, we are naming those as follows:
| Respective Operators |
|---|
| \(\mathbf{A_k} \text{is implemented in } \text{OpDomainLhs}\) |
| \(\mathbf{B_f} \text{is implemented in } \text{OpDomainRhs}\) |
| \(\mathbf{C_k} \text{is implemented in } \text{OpBoundaryLhs}\) |
| \(\mathbf{D_f} \text{is implemented in } \text{OpBoundaryRhs}\) |
The numerical scheme has several similarities with SCL-1: Poisson's equation (homogeneous BC) in terms of coding in MOFEM. Therefore, we are mostly going to discuss the critical parts of it in the next section comparing with the SCL-1: Poisson's equation (homogeneous BC).
Before we explore the implementation details of the User-Defined Operators (UDOs) in "poisson_2d_nonhomogeneous.hpp", let's take a look at the initial lines of code in the file. It includes MoFEM libraries inclusions for finite element functionality and alias definitions for declaration and initialization as follows:
MoFEM has an in-built extensive library consisting of operators and to implement those operators users need to create some new UDOs according to the problems and PDEs. First, let's look at the structure of a UDO so that users can also learn the general structure of any UDOs implemented in MoFEM.
For example, we can notice that our OpDomainLhs is similar to the OpDomainLhsMatrixK UDO in SCL-1: Poisson's equation (homogeneous BC). So, let's identify the differences bellow and then the further will be discussed:
In the above, the class OpDomainLhs is inherited from OpFaceEle, which stands for the alias for the base class MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator. Within this class, there are two public functions: OpDomainLhs() and doWork(). Here, the doWork() does the same as iNtegrate as in SCL-1: Poisson's equation (homogeneous BC). However, the difference is that we are using a different Lambda function to calculate the body source term. Alongside, we have got two private members named locLhs and transLocLhs. This follows the concept of data encapsulation to hide values and objects inside the class as much as possible.
locLhs transLocLhs CHKERR MatSetValues The characterization of OpDomainRhs is straightforward and also aligns with OpDomainRhsVectorF in the context of SCL-1 for solving Poisson's equation with homogeneous BC. However, there is a slight difference in the body source term. The accompanying code snippet illustrates in the following:
In the above snippet, we can find a different body_source term which is calculated by calling a Lambda function sourceTermFunc() with coordinate(x, y, z) i.e. t_coords(0), t_coords(1), and t_coords(2) as arguments as follows:
OpBoundaryLhs is focusing on the boundary( \( \Gamma\)) instead of the domain. So we're using OpEdgeEle instead of OpFaceEle. You can realize the fact easily that for 2D problems any boundary element type should be edge(1D) and for 3D problems it will be face(2D) type element.This operator is quite similar to our first operator, OpDomainLhs. Therefore we are not describing it again. Further, the code snippet is here:
If we compare the OpBoundaryRhs operator, it also closely resembles the OpDomainRhsVectorF operator in its structure. Additionally, it is also similar to OpBoundaryRhs. However, we are now using boundary_term here instead of body_source as follows:
This main class Poisson2DNonhomogeneous contains functions each of which is responsible for a certain task as described in SCL-1. However, in this case, we need to declare two functions as a global variable of sourceTermFunction and boundaryFunction. These functions are applicable for the operators OpDomainRhs and OpBoundaryRhs respectively. In addition to that, we have declared a vector type shared pointer 'boundaryMarker'. We will use and discuss while applying the boundary condition. The code snippet is as follows:
Next, this code defines function in MoFEM to loads and reads your mesh file.
Since we are considering both domain and boundary, we need to add the domain and boundary field in the simpleInterface to setup the problem. Besides, we have selected the polynomial base function as AINSWORTH_BERNSTEIN_BEZIER_BASE and \( \mathbf{H^1}\) scalar space, where we are considering single(1) degree of freedom for each shape function has
Next, we are proceeding to discuss the detail for setting up boundary condition. This code involves marking degrees of freedom on boundary entities, merging markers, and obtaining entities for field splitting in the following:
The line bc_mng->pushMarkDOFsOnEntities(...) is used to collect the degrees of freedom (DOFs) of the entities only for marked block set ("BOUNDARY_CONDITION"). Then, the block sets are accumulated as in container and assembled in boundaryMarker. In addition, the `boundaryEntitiesForFieldsplit is to get the range/IDs for the blocks including vertex, edges or faces. It has been done to indicate the difference between the domain and boundary entities.
At this stage, we are all set to focus on assembling all we did before to the pipeline. However, we find some changes that happens here in codes compared with SCL-2: Poisson's equation (non-homogeneous BC). Before going to that, first we add operators to the OpDomainLhsPipeline using the AddHOOps function considering scalar space H1. Then we need to assemble to the pipeline for the domain using OpSetBc. Then push the main OpDomainLhs operator, and finally unset the assembly using OpUnSetBc. We will discuss it in detail after going through the code snippet in the following:
Here, the process involves using the OpSetBc, OpUnSetBc, and UDOs to push back to the RHS and LHS pipelines for assembling. Specifically, these operations are being carried out on the pipeline for the Domain and boundary in same way i.e. getOpDomainLhsPipeline(), OpBoundaryLhsPipeline. So we are only describing for only one set of assemble processes in the following:
The OpSetBc() object is responsible for assembling only for the boundary entities or everything except the boundary. In other words, when it is flagged 'true' it focuses exclusively on the domain entities and does not account for boundary entities tagged with the 'boundaryMarker.
For example, when we assemble for the boundary, first we use OpSetBc(), then add OpBoundaryLhs and OpUnSetBc to the getOpBoundaryLhsPipeline() and getOpBoundaryLhsPipeline. To be noted that it's important to set the flag as false here.
To recapitulate, when we are assembling for the domain,
OpSetBc()works for skipping(marked'true' to skip boundary) the boundary entities. On the other hand, it considers only the boundaries (marked 'false' to consider) when assembling only for boundary. And then push the UDO. Then OpUnSetBc() is used to set the index to the entities as it was initially.
Afterwards, the computation of LHS and RHS is defined in the setIntegrationRules. The instructions for the functions setIntegrationRules() is same as it was in SCL-1: Poisson's equation (homogeneous BC).
To solve the problem, we remember that we made our assembling systems ready using the set/unseting process in boundaryCondition() for domain and boundaries as shown in Figure: 2 where we can consider the subscripts 'D' is for domain and 'B' for boundary.
At this stage, we are not solving the system of equations directly using KSP solver as we did in SCL-1: Poisson's equation (homogeneous BC). Instead of that, we are proceeding to deploy a preconditioner to the solver named FIELDSPLIT. It splits the fields for domain and boundary as an option - yes/no if needed. The part of the code is as follows:
The algorithm for the above code is to check whether the indices in the field were marked for domain or boundary. It set up the field-split solver using PETSc libraries for boundary entities(boundaryEntitiesForFieldsplit) and the remaining will be for the domain. During the iteration process, IS considers solving the option for domain and boundary so that the index sets are then used to configure the PCFIELDSPLIT solver to solve the appropriate blockset. Finally, the code releases the memory associated with the boundary index set and goes for KSP to solve. And to accumulate all these we followed by VecGhostUpdateBegin(), VecGhostUpdateEnd() and DMoFEMMeshToLocalVector to make things ready before getting the results.
The rest of the code including outputResults() and main function does the same as SCL-2: Poisson's equation (non-homogeneous BC)
Then, run the analysis and get the newly created output file, namely out_result.h5m. Convert it to *.vtk format. Then open it in Paraview and use the filter WarpByScalar, you will be able to see the deformation as bellow:
Further requisite: https://en.wikiversity.org/wiki/Partial_differential_equations/Poisson_Equation.
The plain program for both the implementation of the UDOs (*.hpp) and the main program (*.cpp) are as follows
