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<[email protected]> 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 > [email protected] > https://lists.sourceforge.net/lists/listinfo/libmesh-users > ------------------------------------------------------------------------------ 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
