Hi,
The general technique is to build a finite element on the "face" or
"edge" which is on the boundary and perform the integral (with a loop
over gauss points, i, and j on the face) which appears in your
variational statement. Here is a code fragment
// Handle boundary conditions
for (unsigned int s=0; s<elem>n_sides(); s++)
if (elem>neighbor(s) == NULL) // on the boundary
{
const unsigned int boundary_id =
mesh.boundary_info.boundary_id (elem, s);
AutoPtr<FEBase> fe_face (FEBase::build(dim,
flow_dof_map.variable_type (u_var)));
QGauss qface (dim1, Control::quadrature_order());
fe_face>attach_quadrature_rule (&qface);
fe_face>reinit (elem, s);
const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi();
const std::vector<std::vector<RealGradient> >& dphi_face = fe_face>get_dphi();
const std::vector<Real>& JxW_face = fe_face>get_JxW();
const std::vector<Point>& qface_point = fe_face>get_xyz();
const std::vector<Point>& normals = fe_face>get_normals();
// ...
// T Flux boundary values
else if (user>is_neumann_boundary(boundary_id,
((coupled) ? "rbm" : "thermal"), "T"))
{
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
// Extract the value which needs to be imposed for the flux
// at each Gauss point.
const Real val = user>neumann_boundary_value(qface_point[qp],
normals[qp],
time,
((coupled) ? "rbm" : "thermal"),
"T",
boundary_id);
// Compute the integration on the face
for (unsigned int i=0; i<phi_face.size(); i++)
Ft(i) += JxW_face[qp]*val*phi_face[i][qp];
}
}
}
