2009/9/2 John Peterson <[email protected]> > On Wed, Sep 2, 2009 at 2:38 PM, Vijay S. Mahadevan<[email protected]> > wrote: > > On first look, I would say that the neumann_value is supposed to be > > -0.25 and not -0.5 at x=2. Try this change and see if it resolves the > > bug. Also like John suggested, see if your L2 and H1 errors give right > > convergence orders. > > I don't see anything obviously wrong... > > If the original problem is: d/dx (x * du/dx) = 2/x^2 > > The weak form is: > > -(x u',v') + x u'(2) v(2) - x u'(1) v(1) = (2/x^2,v) > > Bringing the right endpoint bc over to the rhs and assuming the > Dirichlet bc is handled, > > -(x u',v') = (2/x^2,v) - x u'(2) v(2) > > And replacing with the Neumann condition, (-x * du/dx) = 0.5 @ x = 2, we > get > > -(x u',v') = (2/x^2,v) + 0.5 v(2) > > Multiplying thru by neg. 1 as in the code, we get > > (x u',v') = -(2/x^2,v) - 0.5 v(2) > > If the theoretical convergence rates are off, I would probably start > by looking at the endpoint integral, though there may be something > else obvious I'm missing. > > -- > John >
The problem was that I was doing this: Fe(i) += JxW_face[qp] * phi_face[i][qp] * g; where g is the Neumann BC. Whereas, this IS the correct thing to do: Fe(i) += phi_face[i][qp] * g; The difference is there's no multiplication by JxW_face[qp] because applying the Neumann B.C does not require integration. Thanks all. I'm very excited and motivated. Ted ------------------------------------------------------------------------------ Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day trial. Simplify your report design, integration and deployment - and focus on what you do best, core application coding. Discover what's new with Crystal Reports now. http://p.sf.net/sfu/bobj-july _______________________________________________ Libmesh-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
