Essential BC for mixed problem

[markusrenoldner]

Dear community,

How do I implement essential (e.g. homogeneous Dirichlet) boundary conditions for a mixed finite element (involving multiple gridfunctions and FE-spaces) problem that involves a block matrix consisting of several submatrices?

Ideally, I would like to have a command, that sets the degrees of freedom of one FE-space (or rather gridfunction) on the boundary to a certain value, independent of how the resulting matrix will get added to the large system matrix.

I checked the example selection, the website and the documentation, but this was not very fruitful.
Can the resulting matrix from the procedure explained here: https://mfem.org/fem_bc/
just be added to a large blockmatrix?

Thank you in advance

[pazner]

Hi @markusrenoldner,

You can call Operator::FormLinearSystem on the BlockOperator, where the list of essential DOFs indexes into the monolithic block vector. For example, this modification of example 5 solves the Darcy system on a unit square mesh with essential boundary conditions on part of the boundary, and natural boundary conditions on the remainder.

You can run the example with ./ex5 -m ../data/inline-quad.mesh.

(Note you can also do the elimination in the 2x2 system "by hand", e.g. by using FormLinearSystem, FormSystemMatrix, and FormRectangularSystemMatrix on the relevant blocks that need to have DOFs eliminated. You will need to ensure that you modify the RHS appropriately).

diff --git a/examples/ex5.cpp b/examples/ex5.cpp
index a59dde087..955ca5562 100644
--- a/examples/ex5.cpp
+++ b/examples/ex5.cpp
@@ -128,6 +128,12 @@ int main(int argc, char *argv[])
    std::cout << "dim(R+W) = " << block_offsets.Last() << "\n";
    std::cout << "***********************************************************\n";
 
+   Array<int> natural_bdr({1, 1, 0, 0});
+   Array<int> ess_bdr({0, 0, 1, 1});
+
+   Array<int> ess_dofs;
+   R_space->GetEssentialTrueDofs(ess_bdr, ess_dofs);
+
    // 7. Define the coefficients, analytical solution, and rhs of the PDE.
    ConstantCoefficient k(1.0);
 
@@ -149,7 +155,8 @@ int main(int argc, char *argv[])
    LinearForm *fform(new LinearForm);
    fform->Update(R_space, rhs.GetBlock(0), 0);
    fform->AddDomainIntegrator(new VectorFEDomainLFIntegrator(fcoeff));
-   fform->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(fnatcoeff));
+   fform->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(fnatcoeff),
+                                natural_bdr);
    fform->Assemble();
    fform->SyncAliasMemory(rhs);
 
@@ -267,6 +274,12 @@ int main(int argc, char *argv[])
    darcyPrec.SetDiagonalBlock(0, invM);
    darcyPrec.SetDiagonalBlock(1, invS);
 
+   Operator *darcyOpConstrained;
+   Vector X, B;
+   x = 0.0;
+
+   darcyOp.FormLinearSystem(ess_dofs, x, rhs, darcyOpConstrained, X, B);
+
    // 11. Solve the linear system with MINRES.
    //     Check the norm of the unpreconditioned residual.
    int maxIter(1000);
@@ -279,14 +292,16 @@ int main(int argc, char *argv[])
    solver.SetAbsTol(atol);
    solver.SetRelTol(rtol);
    solver.SetMaxIter(maxIter);
-   solver.SetOperator(darcyOp);
+   solver.SetOperator(*darcyOpConstrained);
    solver.SetPreconditioner(darcyPrec);
    solver.SetPrintLevel(1);
-   x = 0.0;
-   solver.Mult(rhs, x);
+   X = 0.0;
+   solver.Mult(B, X);
    if (device.IsEnabled()) { x.HostRead(); }
    chrono.Stop();
 
+   x.Vector::operator=(X);
+
    if (solver.GetConverged())
    {
       std::cout << "MINRES converged in " << solver.GetNumIterations()

[markusrenoldner]

@pazner

Dear Will,
Thank you for your answer.

My problem is to find u $\in$ X and w $\in$ Y st.

$$\begin{bmatrix}A & B \\ B^T & 0\end{bmatrix} \cdot \begin{bmatrix} u^{k+1} \\ w^{k+1} \end{bmatrix} = \begin{bmatrix} A u^{k} + B w^{k} \\ 0 \end{bmatrix}$$

which is computed iteratively (k denote the time steps). I am adapting a working code for the periodic-cube.mesh to a homogeneous Dirichlet scenario on the ref-cube.mesh.

Hence, trying to mimic your 2nd suggestion ("you can also do the elimination in the 2x2 system by hand [...]") - i did the following:

mfem::Array<int> X_essdof, Y_essdof;
mfem::Array<int> ess_bdr(mesh.bdr_attributes.Max()); 
ess_bdr = 1;
X.GetEssentialTrueDofs(ess_bdr, x_essdof); 
Y.GetEssentialTrueDofs(ess_bdr, Y_essdof);

this should set the boundary values to zero:

mfem::GridFunction u(&X);
mfem::GridFunction w(&Y);
u = 0.0;
w = 0.0;

submatrix assembly:

mfem::BilinearForm blf_A(&X);
mfem::SparseMatrix mat_A;
blf_A.AddDomainIntegrator(new mfem::VectorFEMassIntegrator());
blf_A.Assemble();
blf_A.FormSystemMatrix(X_essdof,mat_A);
mat_A.Finalize();

(and similar for mat_B and its transposed mat_B_T)
assembly of big matrix:

mfem::Array<int> offsets (3);
offsets[0] = 0;
offsets[1] = u.Size();
offsets[2] = w.Size();
offsets.PartialSum();
mfem::BlockMatrix M(offsets);
M.SetBlock(0,0, &mat_A);
M.SetBlock(0,1, &mat_B);
M.SetBlock(1,0, &mat_B_T);

Right hand side assembly:

mfem::Vector RHS (u.Size() + w.Size());
mfem::Vector RHS_subvec (u.Size());
mat_A.Mult(u, RHS_subvec);
mat_B.AddMult(w, RHS_subvec);
RHS = 0.0;
RHS.AddSubVector(RHS_subvec, 0);

and then i call MINRES

mfem::MINRES(M, RHS, x, 0, iter, tol, tol);

To summarize: I am using the ess_tdof arrays when calling FormSystemMatrix() on the submatrices. Nothing else was changed in comparison to the implementation on the periodic mesh.

Generally speaking, does this procedure suffice? Does MFEM process everything correctly, or do I also need to eliminate DOFs in the RHS, for example?

[pazner]

FormSystemMatrix performs the elimination on the matrix/operator, whereas FormLinearSystem also performs the elimination of the RHS.

If you use only FormSystemMatrix, then it is your responsibility to eliminate the BCs from the RHS.

Since it looks like your essential boundary conditions are homogeneous, you don't need to move any contribution to the RHS (it would be zero). In the case that you are using inhomogeneous conditions, then you need to move the contribution from the essential DOFs to the RHS.

You also need to modify the values of the RHS vector so that you obtain the correct values at the essential BCs. This will depends on the DiagonalPolicy. In your case, you can just zero out the essential DOFs in the right-hand side.

You can look at the code in ConstrainedOperator::EliminateRHS to see how this elimination is done in the case of DiagonalPolicy::DIAG_ONE (i.e. where the eliminated portion of the matrix is replaced with identity):

mfem/linalg/operator.cpp

Lines 520 to 545 in 2f6eb88

	void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
	{
	w = 0.0;
	const int csz = constraint_list.Size();
	auto idx = constraint_list.Read();
	auto d_x = x.Read();
	// Use read+write access - we are modifying sub-vector of w
	auto d_w = w.ReadWrite();
	MFEM_FORALL(i, csz,
	{
	const int id = idx[i];
	d_w[id] = d_x[id];
	});
	// A.AddMult(w, b, -1.0); // if available to all Operators
	A->Mult(w, z);
	b -= z;
	// Use read+write access - we are modifying sub-vector of b
	auto d_b = b.ReadWrite();
	MFEM_FORALL(i, csz,
	{
	const int id = idx[i];
	d_b[id] = d_x[id];
	});
	}

The best way to be sure that you have achieved the result you're looking for is by testing with a problem with a known solution (e.g. manufactured solution).

[markusrenoldner]

Many thanks for the qick answer @pazner!

One final question. Lets consider the 1D poisson problem with Dirichlet BC (say with bound. values $a$ and $b$).
Either you put the BC to the RHS (and get a smaller system, but you mess with the unknown vector and the matrix)

$$\begin{pmatrix} -2 & 1 & & & \\\ 1 & -2 &1 & & \\\ & &... & & \\\ & & 1& -2&1 \\\ & & & 1&-2 \end{pmatrix} \cdot \begin{pmatrix} u_1 \\\ u_2 \\\ ...\\\ ...\\\ u_n \end{pmatrix} = \begin{pmatrix} f_1 -a\\\ f_2 \\\ ...\\\ ...\\\ f_n-b \end{pmatrix}$$

Or you can change the matrix and keep all DOFs in the unknown vector (larger system, but possibly easier to implement):

$$\begin{pmatrix} 1 & & & & & \\\ 1 & -2 &1 & & &\\\ & 1 &-2 &1 & &\\\ & & &... && \\\ & & & &...& \\\ & & & 1&-2&1 \\ & & & &&1 \end{pmatrix} \cdot \begin{pmatrix} u_0 \\\ u_1 \\\ u_2 \\\ ...\\\ ...\\\ u_n \\\ u_{n+1} \end{pmatrix} = \begin{pmatrix} a\\\ f_1\\\ f_2 \\\ ...\\\ ...\\\ f_n\\\ b \end{pmatrix}$$

Which option does MFEM chose? Is there a difference between FormLinearSystem and FormSystemMatrix concerning the BC implementation?

[pazner]

MFEM implements the larger system, except the off-diagonal entries in the columns associated with the essential BCs are also set to zero (so that the matrix is symmetric), requiring that the right-hand side be modified as in your first example.

[markusrenoldner]

Thank you!