From: John Peterson <peterson@cf...>  20040504 01:24:30

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]; } } } 