From: Ted K. <ted...@go...> - 2009-08-31 20:02:59
|
Hi How do you apply a non-homogeneous Neumann boundary condition? All the examples only deal with Dirichlet B.Cs. Thanks in advance. Ted Kord |
From: David K. <dknez@MIT.EDU> - 2009-08-31 20:09:10
|
Just add the appropriate boundary integral term in a similar way as you do a Dirichlet BC, e.g. something like Fe(i) += JxW_face[qp] * phi_face[i][qp] * g; where g is the Neumann BC. Boundary IDs are useful if you only want to impose the BC on a part of the boundary... Note that you don't include a penalty term in this case. - Dave Ted Kord wrote: > Hi > > How do you apply a non-homogeneous Neumann boundary condition? All the > examples only deal with Dirichlet B.Cs. > > Thanks in advance. > > Ted Kord > ------------------------------------------------------------------------------ > 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 > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: Ted K. <ted...@go...> - 2009-09-01 03:46:45
|
> > Well, to impose a Neumann BC on part of the boundary you should only add > Neumann BC terms to the entries of the rhs vector corresponding to dofs on > the relevant part of the boundary. So you'd need to detect which part of the > boundary you're on and add the Neumann BC terms accordingly. Have a look at > ex13 for an example usage of boundary IDs for this kind of purpose. For more > complicated domains, I typically generate a mesh using gmsh, which enables > you to specify boundary IDs (using a gmsh "physical"). > > In principle 1D is the same except that in that simpler case the Neumann > term is not an integral so you don't need a boundary quadrature rule etc... > Also there are only 2 boundary elements in a 1D mesh so you could easily > determine which boundary element you're on by testing the value of > elem->centroid() or something. > > If you can't get it to work, then just send your BC code fragment to the > list... > > - Dave > This is what I've tried. The Dirichlet BC works fine but the plot of the solution implies that something's terriby wrong with the Neumann B.C : for(unsigned int s=0; s<elem->n_sides(); s++) { { Ke(s,s) += penalty; short int bc_id = mesh.boundary_info->boundary_id (elem,s); if (bc_id == 0) //Left end : Dirichlet BC Fe(s) += 2*penalty; else if (bc_id == 1) //Right end : Neumann BC Fe(s) += - 0.5; } } What am I doing incorrectly? Thanks Ted |
From: Vijay S. M. <vi...@gm...> - 2009-09-01 05:32:39
|
Ted, I assume your dirichlet_value=2.0 and neumann_value=0.5 on your assembly routines for the boundary ? If yes, the boundary assembly is not correct. A generic code for bc specification would look like : for (unsigned int side=0; side<elem->n_sides(); side++) if (elem->neighbor(side) == NULL) { const Real penalty = 1.e10; 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 zf = qface_point[qp](2); if (bc_id == 0) { const Real value = exact_solution(xf, yf, zf); 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) { const Real value = neumann_value(xf, yf, zf); for (unsigned int i=0; i<phi_face.size(); i++) Fe(i) += JxW_face[qp]*value*phi_face[i][qp]; } } Most part of it is taken from ex4.C. Substitute in the neumann_value for the side and you should get the right solution. Just look at the variational form for your problem and this should be easy to see. Vijay On Mon, Aug 31, 2009 at 6:57 PM, Ted Kord<ted...@go...> wrote: >> >> Well, to impose a Neumann BC on part of the boundary you should only add >> Neumann BC terms to the entries of the rhs vector corresponding to dofs on >> the relevant part of the boundary. So you'd need to detect which part of the >> boundary you're on and add the Neumann BC terms accordingly. Have a look at >> ex13 for an example usage of boundary IDs for this kind of purpose. For more >> complicated domains, I typically generate a mesh using gmsh, which enables >> you to specify boundary IDs (using a gmsh "physical"). >> >> In principle 1D is the same except that in that simpler case the Neumann >> term is not an integral so you don't need a boundary quadrature rule etc... >> Also there are only 2 boundary elements in a 1D mesh so you could easily >> determine which boundary element you're on by testing the value of >> elem->centroid() or something. >> >> If you can't get it to work, then just send your BC code fragment to the >> list... >> >> - Dave >> > > > This is what I've tried. The Dirichlet BC works fine but the plot of the > solution implies that something's terriby wrong with the Neumann B.C : > > for(unsigned int s=0; s<elem->n_sides(); s++) > { > { > Ke(s,s) += penalty; > > short int bc_id = mesh.boundary_info->boundary_id (elem,s); > > if (bc_id == 0) //Left end : Dirichlet BC > Fe(s) += 2*penalty; > else if (bc_id == 1) //Right end : Neumann BC > Fe(s) += - 0.5; > > } > } > > What am I doing incorrectly? > > Thanks > > 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 > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: Ted K. <ted...@go...> - 2009-09-01 23:59:49
|
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 Could you possibly show me how to solve this in libmesh with exemplar code. This would really help me get through the boook I'm using to learn the FEM using libmesh. Thank you. Regards Ted |
From: John P. <pet...@cf...> - 2009-09-02 03:30:54
|
On Tue, Sep 1, 2009 at 6:59 PM, Ted Kord<ted...@go...> 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 zero-dimensional quadrature rule and using it to "integrate" that boundary term, as the folks above have suggested. -- John |
From: Ted K. <ted...@go...> - 2009-09-02 04:26:58
|
My apologies The suggestions you guys made actually worked. What I forgot to do the first time round was define the zero-dimensional quadrature rule. Thanks. Ted |
From: John P. <pet...@cf...> - 2009-09-02 04:29:23
|
On Tue, Sep 1, 2009 at 11:26 PM, Ted Kord<ted...@go...> wrote: > My apologies > > The suggestions you guys made actually worked. What I forgot to do the first > time round was define the zero-dimensional quadrature rule. Great, glad to hear it's working! -- John |
From: Ted K. <ted...@go...> - 2009-09-02 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: David K. <dknez@MIT.EDU> - 2009-09-02 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 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 > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: John P. <jwp...@gm...> - 2009-09-02 14:41:18
|
On Wed, Sep 2, 2009 at 9:04 AM, David Knezevic<dk...@mi...> 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: Ted K. <ted...@go...> - 2009-09-02 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(dim-1, 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.e-3; 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: Vijay S. M. <vi...@gm...> - 2009-09-02 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<ted...@go...> 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(dim-1, 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.e-3; > > > 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 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 > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: John P. <jwp...@gm...> - 2009-09-02 20:23:03
|
On Wed, Sep 2, 2009 at 2:38 PM, Vijay S. Mahadevan<vi...@gm...> 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. M. <vi...@gm...> - 2009-09-02 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<jwp...@gm...> wrote: > On Wed, Sep 2, 2009 at 2:38 PM, Vijay S. Mahadevan<vi...@gm...> 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: Ted K. <ted...@go...> - 2009-09-12 11:49:03
|
2009/9/2 John Peterson <jwp...@gm...> > On Wed, Sep 2, 2009 at 2:38 PM, Vijay S. Mahadevan<vi...@gm...> > 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 |
From: Roy S. <roy...@ic...> - 2009-09-12 14:12:25
|
On Sat, 12 Sep 2009, Ted Kord wrote: > Whereas, this IS the correct thing to do: > > Fe(i) += phi_face[i][qp] * g; No, it isn't. > The difference is there's no multiplication by JxW_face[qp] because applying > the Neumann B.C does not require integration. Yes, it does. You've covered up a bug, not found it, I'm afraid. This may work in 1-D (where JxW_face is supposed to simply be 1 on your single "integration point"), but it will give inaccurate results in 2D/3D. > Thanks all. I'm very excited and motivated. Sorry if I've squelched that any, but I didn't want to see your code break later when you move to higher dimensions. --- Roy |
From: Ted K. <ted...@go...> - 2009-09-12 20:56:05
|
2009/9/12 Roy Stogner <roy...@ic...> > > > On Sat, 12 Sep 2009, Ted Kord wrote: > > Whereas, this IS the correct thing to do: >> >> Fe(i) += phi_face[i][qp] * g; >> > > No, it isn't. > > The difference is there's no multiplication by JxW_face[qp] because >> applying >> the Neumann B.C does not require integration. >> > > Yes, it does. You've covered up a bug, not found it, I'm afraid. > This may work in 1-D (where JxW_face is supposed to simply be 1 on > your single "integration point"), but it will give inaccurate results > in 2D/3D. > > Thanks all. I'm very excited and motivated. >> > > Sorry if I've squelched that any, but I didn't want to see your code > break later when you move to higher dimensions. > --- > Roy > It's good you did that 'cause having gone over the code again, it turns out that even though Fe(i) += phi_face[i][qp] * g; works fine (because JxW_face is supposed to simply be 1 on the single "integration point" like you pointed out), the reason I was getting incorrect results before was 'cause I was using JxW[qp] instead of JxW_face[qp] like you mentioned. So, it all ends well, motivation's still there and I've learnt quite a bit. Thanks. Ted |