how to realize hybrid boundary conditions ( Dirichlet and Neumann B.C.)？

[jsfeng-fudan]

Ｈｉ，　Ｉ　ａｍ　ａ　ｎｅｗｅｒ　ｔｏ　ｍｆｅｍ　ｌｉｂ．　Ｉ　ｗａｎｔ　ｔｏ　ｕｓｅ　ｔｗｏ　Ｂ．Ｃ．ｓ　ｉｎ　１Ｄ　ｄｉｆｆｅｒｅｎｔｉａｌ　ｐｒｏｂｌｅｍ．　Ｈｏｗｅｖｅｒ，　Ｉ　ｃａｎ　ｎｏｔ　ｇｅｔ　ｔｈｅ　ｒｉｇｈｔ　ｓｏｌｕｔｉｏｎ．

Ｉ　ｐｏｓｔ　ｍｙ　ｐｒｏｂｌｅｍ　ａｎｄ　ｃｏｄｅ．　Ｃａｎ　ａｎｙｏｎｅ　ｈｅｌｐ　ｍｅ　ｔｏ　ｒｅｓｏｌｖｅ　ｔｈｅ　ｐｒｏｂｌｅｍ？
// Problem: (1) -y''[x] = 1, y[0]=1, y'[1]= -1; two boundary ( Drichlet and Neumann B.C. )

#include "mfem.hpp"
#include
#include

using namespace std;
using namespace mfem;

int main(int argc, char *argv[])
{

int order = 1;

double sx = 1.0;
Mesh mesh = Mesh::MakeCartesian1D(10, sx);

H1_FECollection fec(order, mesh.Dimension());
cout << fec.GetBasisType() << endl;
cout <<"order : " << fec.GetOrder() << endl;
FiniteElementSpace fespace(&mesh, &fec);

cout << "Number of finite element unknowns: " << fespace.GetTrueVSize() << endl;

Array ess_tdof_list;

if( mesh.bdr_attributes.Size())
{
Array ess_bdr(mesh.bdr_attributes.Max());
ess_bdr = 1;
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}

LinearForm b(&fespace);
LinearForm b1(&fespace);

ConstantCoefficient one(1.0);

b.AddDomainIntegrator(new DomainLFIntegrator(one));
b.Assemble();
b.Print();

cout << "after adding boundary condition .... " << endl;
int dim = 1;
VectorArrayCoefficient f(dim);
for (int i = 0; i < dim-1; i++)
{
f.Set(i, new ConstantCoefficient(0.0));
}
{
Vector pull_force(mesh.bdr_attributes.Max());
cout << "bdr_attributes.Max " <<mesh.bdr_attributes.Max() << endl;
pull_force = 0.0;
pull_force(1) = -1.0;
f.Set(dim-1, new PWConstCoefficient(pull_force));
}
b1.AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));

b1.Assemble();
b1.Print();
b.Add(1,b1).Print();
b=b.Add(1,b1);

GridFunction x(&fespace); // need to satisfy BC
x = 0.0;
x[ess_tdof_list[0]] = 1.0; // left B.C. x[0] = 1.0;

BilinearForm a(&fespace);
a.AddDomainIntegrator( new DiffusionIntegrator(one));
//a.Print();
cout << a.Size() << endl;
//cout << "After assemble (mass matrix) : " << endl;
a.Assemble();

OperatorPtr A;
Vector B, X;
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);

cout << "Size of linear system : " << A->Height() << endl;

GSSmoother M( (SparseMatrix&) (*A));

PCG(*A,M,B,X,1, 200, 1e-12,0.0);
cout << "recover the solution on nodes ...." << endl;
a.RecoverFEMSolution(X,b,x);
x.Save("sol.gf");
mesh.Save("mesh.mesh");

return 0;
}

[psocratis]

Hi @jsfeng-fudan, in order to specify different kinds of BC on different parts of the boundary you need to use bdr_attributes.
In the 1D case, you have the left vertex (which has to be marked as essential) and the right one for which you have an additional Bdr integrator.

Please see below for the solution to the minimal example you posted above

#include "mfem.hpp"
using namespace std;
using namespace mfem;

double y_exact(const Vector & x)
{
   return -x(0)*x(0)/2.0+1;
}

int main(int argc, char *argv[])
{

   int order = 1;
   double sx = 1.0;
   Mesh mesh = Mesh::MakeCartesian1D(10, sx);

   H1_FECollection fec(order, mesh.Dimension());
   FiniteElementSpace fespace(&mesh, &fec);

   cout << "Number of finite element unknowns: " << fespace.GetTrueVSize() << endl;

   Array<int> ess_tdof_list;
   Array<int> ess_bdr;
   Array<int> newmann_bdr;
   if( mesh.bdr_attributes.Size())
   {
      ess_bdr.SetSize(mesh.bdr_attributes.Max());
      newmann_bdr.SetSize(mesh.bdr_attributes.Max());
      ess_bdr[0] = 1;
      ess_bdr[1] = 0;
      newmann_bdr[0] = 0;
      newmann_bdr[1] = 1;
      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
   }

   LinearForm b(&fespace);
   ConstantCoefficient one(1.0);
   ConstantCoefficient negone(-1.0);

   b.AddDomainIntegrator(new DomainLFIntegrator(one));
   b.AddBoundaryIntegrator(new BoundaryLFIntegrator(negone),newmann_bdr);
   b.Assemble();

   GridFunction x(&fespace); 
   x = 0.0;
   x.ProjectBdrCoefficient(one,ess_bdr);

   BilinearForm a(&fespace);
   a.AddDomainIntegrator(new DiffusionIntegrator(one));
   a.Assemble();

   OperatorPtr A;
   Vector B, X;
   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);

   cout << "Size of linear system : " << A->Height() << endl;

   GSSmoother M( (SparseMatrix&) (*A));

   PCG(*A,M,B,X,1, 200, 1e-12,0.0);
   a.RecoverFEMSolution(X,b,x);

   char vishost[] = "localhost";
   int  visport   = 19916;
   socketstream sol_sock(vishost, visport);
   sol_sock.precision(8);
   sol_sock << "solution\n" << mesh << x << flush;


   // exact solution
   GridFunction y_ex(&fespace);
   FunctionCoefficient Ex(y_exact);
   y_ex.ProjectCoefficient(Ex);

   socketstream exact_sock(vishost, visport);
   exact_sock.precision(8);
   exact_sock << "solution\n" << mesh << y_ex << flush;

   return 0;
}

I hope this helps.

Best,

Socratis

[jsfeng-fudan]

Thank you every much @psocratis . I have another problem, which is similar to this problem, but now there is no Dirichlet B.C.
The problem is :
-y' ' [x] = 3 in domain [0,1], with the Robin BCs : y'[0]-y[0] = -1 and y'[1]+y[1] = 1;
my solution is

correct solution is :

I attach my key code for enforcing the Robin B.C.

rbc_bdr = 1;
fespace.GetEssentialTrueDofs(rbc_bdr, ess_tdof_list);

double rbc_a_val1 = -1.0; // for left endpoint
double rbc_b_val1 = -1.0; 
double rbc_a_val2 = 1.0;  // for right endpoint
double rbc_b_val2 = 1.0; 
ConstantCoefficient rbc_a_coeff1(rbc_a_val1);// entering a matrix
ConstantCoefficient rbc_a_coeff2(rbc_a_val2);// entering a matrix
ConstantCoefficient rbc_b_coeff1(rbc_b_val1); // entering b linear form
ConstantCoefficient rbc_b_coeff2(rbc_b_val2); // entering b linear form           
b.AddDomainIntegrator(new DomainLFIntegrator(f_src));
// key steps: add B.C. into linear system
// Add Robin BCs n.(matCoef Grad u) + rbcACoef u = rbcBCoef on the boundary
// // marked in the rbc_marker array.
// left
rbc_bdr = 0;
rbc_bdr[0] =1 ;
b.AddBoundaryIntegrator(new BoundaryLFIntegrator(rbc_b_coeff1), rbc_bdr);
rbc_bdr = 0;
rbc_bdr[1] =1 ;
b.AddBoundaryIntegrator(new BoundaryLFIntegrator(rbc_b_coeff2), rbc_bdr);                                   
b.Assemble();
//getchar();
a.AddDomainIntegrator( new DiffusionIntegrator(one));
// // Add Robin BCs n.(matCoef grad u) + rbcACoef u = rbcBCoef on the boundary
// // marked in the rbc_marker array.
rbc_bdr = 0;
rbc_bdr[0] = 1;
a.AddBoundaryIntegrator(new MassIntegrator(rbc_a_coeff1), rbc_bdr);
rbc_bdr = 0;
rbc_bdr[1] = 1;
a.AddBoundaryIntegrator(new MassIntegrator(rbc_a_coeff2), rbc_bdr);       
cout << a.Size() << endl;
a.Assemble();

can you help to solve this problem ?

[psocratis]

please see #3121