You can subscribe to this list here.
2003 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}
(2) 
_{Oct}
(2) 
_{Nov}
(27) 
_{Dec}
(31) 

2004 
_{Jan}
(6) 
_{Feb}
(15) 
_{Mar}
(33) 
_{Apr}
(10) 
_{May}
(46) 
_{Jun}
(11) 
_{Jul}
(21) 
_{Aug}
(15) 
_{Sep}
(13) 
_{Oct}
(23) 
_{Nov}
(1) 
_{Dec}
(8) 
2005 
_{Jan}
(27) 
_{Feb}
(57) 
_{Mar}
(86) 
_{Apr}
(23) 
_{May}
(37) 
_{Jun}
(34) 
_{Jul}
(24) 
_{Aug}
(17) 
_{Sep}
(50) 
_{Oct}
(24) 
_{Nov}
(10) 
_{Dec}
(60) 
2006 
_{Jan}
(47) 
_{Feb}
(46) 
_{Mar}
(127) 
_{Apr}
(19) 
_{May}
(26) 
_{Jun}
(62) 
_{Jul}
(47) 
_{Aug}
(51) 
_{Sep}
(61) 
_{Oct}
(42) 
_{Nov}
(50) 
_{Dec}
(33) 
2007 
_{Jan}
(60) 
_{Feb}
(55) 
_{Mar}
(77) 
_{Apr}
(102) 
_{May}
(82) 
_{Jun}
(102) 
_{Jul}
(169) 
_{Aug}
(117) 
_{Sep}
(80) 
_{Oct}
(37) 
_{Nov}
(51) 
_{Dec}
(43) 
2008 
_{Jan}
(71) 
_{Feb}
(94) 
_{Mar}
(98) 
_{Apr}
(125) 
_{May}
(54) 
_{Jun}
(119) 
_{Jul}
(60) 
_{Aug}
(111) 
_{Sep}
(118) 
_{Oct}
(125) 
_{Nov}
(119) 
_{Dec}
(94) 
2009 
_{Jan}
(109) 
_{Feb}
(38) 
_{Mar}
(93) 
_{Apr}
(88) 
_{May}
(29) 
_{Jun}
(57) 
_{Jul}
(53) 
_{Aug}
(48) 
_{Sep}
(68) 
_{Oct}
(151) 
_{Nov}
(23) 
_{Dec}
(35) 
2010 
_{Jan}
(84) 
_{Feb}
(60) 
_{Mar}
(184) 
_{Apr}
(112) 
_{May}
(60) 
_{Jun}
(90) 
_{Jul}
(23) 
_{Aug}
(70) 
_{Sep}
(119) 
_{Oct}
(27) 
_{Nov}
(47) 
_{Dec}
(54) 
2011 
_{Jan}
(22) 
_{Feb}
(19) 
_{Mar}
(92) 
_{Apr}
(93) 
_{May}
(35) 
_{Jun}
(91) 
_{Jul}
(32) 
_{Aug}
(61) 
_{Sep}
(7) 
_{Oct}
(69) 
_{Nov}
(81) 
_{Dec}
(23) 
2012 
_{Jan}
(64) 
_{Feb}
(95) 
_{Mar}
(35) 
_{Apr}
(36) 
_{May}
(63) 
_{Jun}
(98) 
_{Jul}
(70) 
_{Aug}
(171) 
_{Sep}
(149) 
_{Oct}
(64) 
_{Nov}
(67) 
_{Dec}
(126) 
2013 
_{Jan}
(108) 
_{Feb}
(104) 
_{Mar}
(171) 
_{Apr}
(133) 
_{May}
(108) 
_{Jun}
(100) 
_{Jul}
(93) 
_{Aug}
(126) 
_{Sep}
(74) 
_{Oct}
(59) 
_{Nov}
(145) 
_{Dec}
(93) 
2014 
_{Jan}
(38) 
_{Feb}
(45) 
_{Mar}
(26) 
_{Apr}
(41) 
_{May}
(125) 
_{Jun}
(70) 
_{Jul}
(61) 
_{Aug}
(66) 
_{Sep}
(60) 
_{Oct}
(110) 
_{Nov}
(27) 
_{Dec}
(30) 
2015 
_{Jan}
(43) 
_{Feb}
(67) 
_{Mar}
(71) 
_{Apr}
(92) 
_{May}
(39) 
_{Jun}
(15) 
_{Jul}
(46) 
_{Aug}
(63) 
_{Sep}
(84) 
_{Oct}
(82) 
_{Nov}
(69) 
_{Dec}
(45) 
2016 
_{Jan}
(92) 
_{Feb}
(91) 
_{Mar}
(148) 
_{Apr}
(43) 
_{May}
(58) 
_{Jun}
(117) 
_{Jul}
(92) 
_{Aug}
(140) 
_{Sep}
(49) 
_{Oct}
(33) 
_{Nov}
(85) 
_{Dec}
(40) 
2017 
_{Jan}
(41) 
_{Feb}
(36) 
_{Mar}
(49) 
_{Apr}
(41) 
_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 



1
(5) 
2
(11) 
3
(3) 
4
(6) 
5
(1) 
6

7

8

9
(1) 
10
(4) 
11

12
(9) 
13
(2) 
14
(6) 
15
(3) 
16

17

18

19

20

21

22

23
(1) 
24
(2) 
25
(1) 
26

27

28

29
(8) 
30
(5) 



From: Roy Stogner <roystgnr@ic...>  20090902 23:10:03

If you're using a recent libMesh SVN head, now's the time for you to do an svn update. There was a really stupid bug I introduced a while back when adding ghosted vector support, which could lead to corrupted solutions for certain parallel partitioning / AMR combinations. Ironically, this isn't a ghosted vector bug; in fact without the new ghost vector assertions it would have been an order of magnitude harder to debug.  Roy 
From: Vijay S. Mahadevan <vijay.m@gm...>  20090902 21:37:55

John, good catch. I was wrong in my previous email and I apologize for confusing anybody (Ted) following the thread ! I made the mistake of subtituting u'(2) = 0.25 but then forgot about the x before this term which brings it back to 0.5. Its only the middle of the week and I'm already erratic. I blame the summer heat... On Wed, Sep 2, 2009 at 3:22 PM, John Peterson<jwpeterson@...> wrote: > On Wed, Sep 2, 2009 at 2:38 PM, Vijay S. Mahadevan<vijay.m@...> 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 > 
From: John Peterson <jwpeterson@gm...>  20090902 20:23:03

On Wed, Sep 2, 2009 at 2:38 PM, Vijay S. Mahadevan<vijay.m@...> 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 
From: Vijay S. Mahadevan <vijay.m@gm...>  20090902 19:39:11

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. Vijay On Wed, Sep 2, 2009 at 2:25 PM, Ted Kord<teddy.kord@...> wrote: > This is my assemble function. I can't see that I'm doing anything wrong but > then again, I'm new to the library. > > void assemble_poisson(EquationSystems& es, > const std::string& system_name) > { > libmesh_assert (system_name == "Poisson"); > > const MeshBase& mesh = es.get_mesh(); > > const unsigned int dim = mesh.mesh_dimension(); > > LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> > ("Poisson"); > > const DofMap& dof_map = system.get_dof_map(); > > FEType fe_type = dof_map.variable_type(0); > > AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); > > QGauss qrule (dim, FIFTH); > > fe>attach_quadrature_rule (&qrule); > > AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); > > QGauss qface(dim1, FIFTH); > > fe_face>attach_quadrature_rule (&qface); > > const std::vector<Real>& JxW = fe>get_JxW(); > > const std::vector<Point>& q_point = fe>get_xyz(); > > const std::vector<std::vector<Real> >& phi = fe>get_phi(); > > const std::vector<std::vector<RealGradient> >& dphi = fe>get_dphi(); > > DenseMatrix<Number> Ke; > DenseVector<Number> Fe; > > > std::vector<unsigned int> dof_indices; > > MeshBase::const_element_iterator el = > mesh.active_local_elements_begin(); > const MeshBase::const_element_iterator end_el = > mesh.active_local_elements_end(); > > for ( ; el != end_el ; ++el) > { > const Elem* elem = *el; > > dof_map.dof_indices (elem, dof_indices); > > fe>reinit (elem); > > Ke.resize (dof_indices.size(), > dof_indices.size()); > > Fe.resize (dof_indices.size()); > > for (unsigned int qp=0; qp<qrule.n_points(); qp++) > { > for (unsigned int i=0; i<phi.size(); i++) > for (unsigned int j=0; j<phi.size(); j++) > Ke(i,j) += > JxW[qp]*(q_point[qp](0)*dphi[i][qp]*dphi[j][qp]); > { > const Real x = q_point[qp](0); > //const Real y = q_point[qp](1); > //const Real eps = 1.e3; > > > const Real fxy = 2/(q_point[qp](0) * q_point[qp](0)); > for (unsigned int i=0; i<phi.size(); i++) > Fe(i) += JxW[qp]*fxy*phi[i][qp]; > } > } > > { > > for (unsigned int side=0; side<elem>n_sides(); side++){ > short int bc_id = mesh.boundary_info>boundary_id > (elem,side); > if (elem>neighbor(side) == NULL) > { > const std::vector<std::vector<Real> >& phi_face = > fe_face>get_phi(); > > const std::vector<Real>& JxW_face = fe_face>get_JxW(); > > const std::vector<Point >& qface_point = > fe_face>get_xyz(); > > fe_face>reinit(elem, side); > > for (unsigned int qp=0; qp<qface.n_points(); qp++) > { > const Real xf = qface_point[qp](0); > const Real yf = qface_point[qp](1); > > > const Real penalty = 1.e10; > > if (bc_id == 0) { //Dirichlet BC > // The boundary value. > const Real value = 2.0;//exact_solution(xf, yf); > > for (unsigned int i=0; i<phi_face.size(); i++) > for (unsigned int j=0; j<phi_face.size(); > j++) > Ke(i,j) += > JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp]; > > for (unsigned int i=0; i<phi_face.size(); i++) > Fe(i) += > JxW_face[qp]*penalty*value*phi_face[i][qp]; > } > else if (bc_id == 1) //Neumann BC > { > const Real value = 0.5;//neumann_value(xf, yf, > zf); > for (unsigned int i=0; i<phi.size(); i++) > Fe(i) += JxW[qp]*phi[i][qp]*value; > } > } > } > } > } > > dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); > > system.matrix>add_matrix (Ke, dof_indices); > system.rhs>add_vector (Fe, dof_indices); > } > } > > Regards > > Ted >  > Let Crystal Reports handle the reporting  Free Crystal Reports 2008 30Day > 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/bobjjuly > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 
From: Ted Kord <teddy.kord@go...>  20090902 19:25:27

This is my assemble function. I can't see that I'm doing anything wrong but then again, I'm new to the library. void assemble_poisson(EquationSystems& es, const std::string& system_name) { libmesh_assert (system_name == "Poisson"); const MeshBase& mesh = es.get_mesh(); const unsigned int dim = mesh.mesh_dimension(); LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson"); const DofMap& dof_map = system.get_dof_map(); FEType fe_type = dof_map.variable_type(0); AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); QGauss qrule (dim, FIFTH); fe>attach_quadrature_rule (&qrule); AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qface(dim1, FIFTH); fe_face>attach_quadrature_rule (&qface); const std::vector<Real>& JxW = fe>get_JxW(); const std::vector<Point>& q_point = fe>get_xyz(); const std::vector<std::vector<Real> >& phi = fe>get_phi(); const std::vector<std::vector<RealGradient> >& dphi = fe>get_dphi(); DenseMatrix<Number> Ke; DenseVector<Number> Fe; std::vector<unsigned int> dof_indices; MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); for ( ; el != end_el ; ++el) { const Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); fe>reinit (elem); Ke.resize (dof_indices.size(), dof_indices.size()); Fe.resize (dof_indices.size()); for (unsigned int qp=0; qp<qrule.n_points(); qp++) { for (unsigned int i=0; i<phi.size(); i++) for (unsigned int j=0; j<phi.size(); j++) Ke(i,j) += JxW[qp]*(q_point[qp](0)*dphi[i][qp]*dphi[j][qp]); { const Real x = q_point[qp](0); //const Real y = q_point[qp](1); //const Real eps = 1.e3; const Real fxy = 2/(q_point[qp](0) * q_point[qp](0)); for (unsigned int i=0; i<phi.size(); i++) Fe(i) += JxW[qp]*fxy*phi[i][qp]; } } { for (unsigned int side=0; side<elem>n_sides(); side++){ short int bc_id = mesh.boundary_info>boundary_id (elem,side); if (elem>neighbor(side) == NULL) { const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi(); const std::vector<Real>& JxW_face = fe_face>get_JxW(); const std::vector<Point >& qface_point = fe_face>get_xyz(); fe_face>reinit(elem, side); for (unsigned int qp=0; qp<qface.n_points(); qp++) { const Real xf = qface_point[qp](0); const Real yf = qface_point[qp](1); const Real penalty = 1.e10; if (bc_id == 0) { //Dirichlet BC // The boundary value. const Real value = 2.0;//exact_solution(xf, yf); for (unsigned int i=0; i<phi_face.size(); i++) for (unsigned int j=0; j<phi_face.size(); j++) Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp]; for (unsigned int i=0; i<phi_face.size(); i++) Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp]; } else if (bc_id == 1) //Neumann BC { const Real value = 0.5;//neumann_value(xf, yf, zf); for (unsigned int i=0; i<phi.size(); i++) Fe(i) += JxW[qp]*phi[i][qp]*value; } } } } } dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); system.matrix>add_matrix (Ke, dof_indices); system.rhs>add_vector (Fe, dof_indices); } } Regards Ted 
From: John Peterson <jwpeterson@gm...>  20090902 14:41:18

On Wed, Sep 2, 2009 at 9:04 AM, David Knezevic<dknez@...> wrote: > The Neumann condition is only imposed weakly, so you might have to > refine the mesh further to get a more accurate solution... if that > doesn't work then there's probably still a bug in your code. Yes, please also check that the theoretical convergence rates of the L2 and H1 norms under uniform refinement are correct. (Actually compute the integral norms, not the discrete l2/h1 norms, of the error.)  John 
From: David Knezevic <dknez@MIT.EDU>  20090902 14:04:58

The Neumann condition is only imposed weakly, so you might have to refine the mesh further to get a more accurate solution... if that doesn't work then there's probably still a bug in your code.  Dave Ted Kord wrote: > Having said that it works, the error compared to the exact solution > (0.5ln(x) + 2/x) is still significant though particularly towards the end > where the Neumann BC was imposed. > > Ted >  > Let Crystal Reports handle the reporting  Free Crystal Reports 2008 30Day > 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/bobjjuly > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 
From: Ted Kord <teddy.kord@go...>  20090902 05:04:05

Having said that it works, the error compared to the exact solution (0.5ln(x) + 2/x) is still significant though particularly towards the end where the Neumann BC was imposed. Ted 
From: John Peterson <peterson@cf...>  20090902 04:29:23

On Tue, Sep 1, 2009 at 11:26 PM, Ted Kord<teddy.kord@...> wrote: > My apologies > > The suggestions you guys made actually worked. What I forgot to do the first > time round was define the zerodimensional quadrature rule. Great, glad to hear it's working!  John 
From: Ted Kord <teddy.kord@go...>  20090902 04:26:58

My apologies The suggestions you guys made actually worked. What I forgot to do the first time round was define the zerodimensional quadrature rule. Thanks. Ted 
From: John Peterson <peterson@cf...>  20090902 03:30:54

On Tue, Sep 1, 2009 at 6:59 PM, Ted Kord<teddy.kord@...> wrote: > Thanks guys but the alterations in code that you gave are not working > either. > > This is the actual equation I'm trying to solve: > > d/dx (x * du/dx) = 2/x^2 on the domain 1 < x < 2 > > subject to these boundary conditions: > > u(1) = 2 > (x * du/dx) = 0.5 @ x = 2 > The boundary terms are x u'(2) v(2)  x u'(1) v(1) where v is a test function. At the x=1 end we enforce the Dirichlet condition via the penalty method. At the x=2 end, We replace x u'(2) by 0.5. There is no integral in the usual sense since this is 1D, so all you really need to do is subtract 0.5 from the equation corresponding to the rightmost endpoint to implement that BC. You can also do the same thing by creating a zerodimensional quadrature rule and using it to "integrate" that boundary term, as the folks above have suggested.  John 