From: jurak <jurak@ma...>  20081128 19:30:22

Dear all, I was testing ExactSolution class by a simple program that solves Poisson equation on a unite rectangle, using P1 elements on triangles, and calculates convergence rate of the method by using uniform refinement of the initial mesh. To my surprise I obtained the following errors and rates: E R R O R C O N V E R G E N C E R A T E H1 L2 Linf H1 L2 Linf ===== ===== ======= ===== ===== ====== 3.469577e01 1.363935e02 3.209919e02 1.744062e01 4.888679e03 1.110995e02 0.992308 1.480259 1.530684 8.753062e02 2.352227e03 5.084924e03 0.994592 1.055418 1.127554 4.406352e02 1.495102e03 3.113305e03 0.990203 0.653783 0.707779 The convergence in L2 norm is too slow (should be of second order). My guess is that I am missing here something fairly obvious, but I can't find what. The assembly function is taken from example 2, and the refining loop is given here: for(int refinement=0; refinement < max_refinement; ++ refinement) { system.solve(); ex_sol.compute_error("Poisson","p"); l2_error = ex_sol.l2_error("Poisson","p"); h1_error = ex_sol.h1_error("Poisson","p"); linf_error = ex_sol.l_inf_error("Poisson","p"); out_file << std::scientific << h1_error << " " << l2_error << " "<< linf_error<< " "; if(refinement >0) out_file << std::fixed << std::log2(h1_error_old/h1_error)<< " " << std::log2(l2_error_old/l2_error) << " " << std::log2(linf_error_old/linf_error); out_file <<std::endl; h1_error_old = h1_error; l2_error_old = l2_error; linf_error_old = linf_error; mesh_refinement.uniformly_refine(); es.reinit(); } I give the whole program in the attachment. Thanks for your help. Mladen Jurak <http://web.math.hr/%7Ejurak>; 
From: Benjamin Kirk <benjamin.kirk@na...>  20081128 21:02:09

I don't seem to see the attachment... But I have to ask: If the assembly function was taken from example 3 (2 does no assembly(?)), are you constraining the hanging degrees of freedom? That is, is there a line like // We have now built the element matrix and RHS vector in terms // of the element degrees of freedom. However, it is possible // that some of the element DOFs are constrained to enforce // solution continuity, i.e. they are not really "free". We need // to constrain those DOFs in terms of nonconstrained DOFs to // ensure a continuous solution. The // \p DofMap::constrain_element_matrix_and_vector() method does // just that. dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); (c.f. example 10) Ben On 11/28/08 1:32 PM, "jurak" <jurak@...> wrote: > Dear all, > > I was testing ExactSolution class by a simple program that solves > Poisson equation > on a unite rectangle, using P1 elements on triangles, and calculates > convergence > rate of the method by using uniform refinement of the initial mesh. > > To my surprise I obtained the following errors and rates: > > E R R O R C > O N V E R G E N C E R A T E > H1 L2 Linf > H1 L2 Linf > ===== ===== ======= ===== ===== ====== > 3.469577e01 1.363935e02 3.209919e02 > 1.744062e01 4.888679e03 1.110995e02 0.992308 1.480259 1.530684 > 8.753062e02 2.352227e03 5.084924e03 0.994592 1.055418 1.127554 > 4.406352e02 1.495102e03 3.113305e03 0.990203 0.653783 0.707779 > > The convergence in L2 norm is too slow (should be of second order). > > My guess is that I am missing here something fairly obvious, but I can't > find what. > The assembly function is taken from example 2, and the refining loop is > given here: > > for(int refinement=0; refinement < max_refinement; ++ refinement) > { > > system.solve(); > > ex_sol.compute_error("Poisson","p"); > > l2_error = ex_sol.l2_error("Poisson","p"); > h1_error = ex_sol.h1_error("Poisson","p"); > linf_error = ex_sol.l_inf_error("Poisson","p"); > > out_file << std::scientific > << h1_error << " " << l2_error << " "<< linf_error<< " "; > if(refinement >0) > out_file << std::fixed > << std::log2(h1_error_old/h1_error)<< " " > << std::log2(l2_error_old/l2_error) << " " > << std::log2(linf_error_old/linf_error); > out_file <<std::endl; > > h1_error_old = h1_error; > l2_error_old = l2_error; > linf_error_old = linf_error; > > mesh_refinement.uniformly_refine(); > es.reinit(); > } > > I give the whole program in the attachment. Thanks for your help. > > Mladen Jurak > > > > <http://web.math.hr/%7Ejurak>; 
From: Roy Stogner <roystgnr@ic...>  20081129 06:06:32

On Fri, 28 Nov 2008, Benjamin Kirk wrote: > I don't seem to see the attachment... But I have to ask: > > If the assembly function was taken from example 3 (2 does no assembly(?)), > are you constraining the hanging degrees of freedom? That shouldn't be the problem, if he's just doing uniform refinement.  Roy 
From: Benjamin Kirk <benjamin.kirk@na...>  20081129 17:13:07

>> I don't seem to see the attachment... But I have to ask: >> >> If the assembly function was taken from example 3 (2 does no assembly(?)), >> are you constraining the hanging degrees of freedom? > > That shouldn't be the problem, if he's just doing uniform refinement. Roy's right, that's not the problem, but I think this is: MeshBase::const_element_iterator el = mesh.local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.local_elements_end(); for ( ; el != end_el; ++el) { ... What happens if you change the loop to MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); Note that in the code you sent you are assembling the system for *all* the elements, including the coarser mesh levels. What does that do? Ben 
From: Roy Stogner <roystgnr@ic...>  20081129 19:42:31

On Sat, 29 Nov 2008, Benjamin Kirk wrote: > Roy's right, that's not the problem, but I think this is: > > MeshBase::const_element_iterator el = > mesh.local_elements_begin(); > What happens if you change the loop to > > MeshBase::const_element_iterator el = > mesh.active_local_elements_begin(); That'd do it. Didn't David or someone do this once too? It's an easy mistake to make. We ought to add a warning message somehow... maybe when an FE object is reinit'ed using a grandparent or higher element?  Roy 
From: jurak <jurak@ma...>  20081129 19:48:24

Benjamin Kirk wrote: > Roy's right, that's not the problem, but I think this is: > > MeshBase::const_element_iterator el = > mesh.local_elements_begin(); > const MeshBase::const_element_iterator end_el = > mesh.local_elements_end(); > > for ( ; el != end_el; ++el) > { > ... > > What happens if you change the loop to > > MeshBase::const_element_iterator el = > mesh.active_local_elements_begin(); > const MeshBase::const_element_iterator end_el = > mesh.active_local_elements_end(); > > > Note that in the code you sent you are assembling the system for *all* the > elements, including the coarser mesh levels. > > What does that do? > > Ben > That's it. With /active_local iterators /I get second order convergence in L2 norm. Thank's for that. But, it seems that Example 3 has the same problem since (current version) uses *MeshBase*::const_element_iterator el = mesh.elements_begin(); *const* MeshBase::const_element_iterator end_el = mesh.elements_end() in order to iterate through the elements. Mladen 
From: Benjamin Kirk <benjamin.kirk@na...>  20081129 19:52:33

>> What happens if you change the loop to >> >> MeshBase::const_element_iterator el = >> mesh.active_local_elements_begin(); >> const MeshBase::const_element_iterator end_el = >> mesh.active_local_elements_end(); >> >> >> Note that in the code you sent you are assembling the system for *all* the >> elements, including the coarser mesh levels. >> >> What does that do? >> > That's it. With /active_local iterators /I get second order convergence > in L2 norm. > Thank's for that. > > But, it seems that Example 3 has the same problem since (current > version) uses ... Strictly speaking, ex3 is immune to the problem because we do not introduce adaptivity until ex9 (uniform refinement) and ex10 (adaptive refinement). In retrospect, perhaps we should introduce the active element iterators immediately in ex3  it is not like they are a complicated notion or anything. Ben 
From: John Peterson <jwpeterson@gm...>  20081129 21:49:33

On Sat, Nov 29, 2008 at 1:42 PM, Roy Stogner <roystgnr@...> wrote: > > On Sat, 29 Nov 2008, Benjamin Kirk wrote: > >> Roy's right, that's not the problem, but I think this is: >> >> MeshBase::const_element_iterator el = >> mesh.local_elements_begin(); > >> What happens if you change the loop to >> >> MeshBase::const_element_iterator el = >> mesh.active_local_elements_begin(); > > That'd do it. Didn't David or someone do this once too? It's an easy > mistake to make. We ought to add a warning message somehow... maybe > when an FE object is reinit'ed using a grandparent or higher element? I've wondered about what we could do about this problem as well. At the time of instantiation of a regular element iterator, in debug mode, we might be able to check the mesh for the presence of ancestor elements, and warn that you may be iterating over an improper subset of elements.  John 
From: Roy Stogner <roystgnr@ic...>  20081129 22:16:45

On Sat, 29 Nov 2008, John Peterson wrote: > On Sat, Nov 29, 2008 at 1:42 PM, Roy Stogner <roystgnr@...> wrote: > >> That'd do it. Didn't David or someone do this once too? It's an easy >> mistake to make. We ought to add a warning message somehow... maybe >> when an FE object is reinit'ed using a grandparent or higher element? > > I've wondered about what we could do about this problem as well. At > the time of instantiation of a regular element iterator, in debug > mode, we might be able to check the mesh for the presence of ancestor > elements, and warn that you may be iterating over an improper subset > of elements. That would probably go too far  I think there's a lot of legit library level code that needs to iterate over all elements. My first thought was to warn on reinit() or dof_indices() of nonactive elements, but that might get falsely triggered by library projection code. I think grandfather elements would be safe, but that wouldn't issue a warning until after a second level of refinement was done.  Roy 
From: John Peterson <jwpeterson@gm...>  20081130 02:56:38

On Sat, Nov 29, 2008 at 4:16 PM, Roy Stogner <roystgnr@...> wrote: > > > On Sat, 29 Nov 2008, John Peterson wrote: > >> On Sat, Nov 29, 2008 at 1:42 PM, Roy Stogner <roystgnr@...> >> wrote: >> >>> That'd do it. Didn't David or someone do this once too? It's an easy >>> mistake to make. We ought to add a warning message somehow... maybe >>> when an FE object is reinit'ed using a grandparent or higher element? >> >> I've wondered about what we could do about this problem as well. At >> the time of instantiation of a regular element iterator, in debug >> mode, we might be able to check the mesh for the presence of ancestor >> elements, and warn that you may be iterating over an improper subset >> of elements. > > That would probably go too far  I think there's a lot of legit > library level code that needs to iterate over all elements. > > My first thought was to warn on reinit() or dof_indices() of > nonactive elements, but that might get falsely triggered by > library projection code. How about a flag in FEBase we can toggle to turn on/off warnings about calling reinit on nonactive elements? By default it would always warn, but we could flip the flag during library routines where we really need to do reinitialization of nonactive elements.  John 
From: jurak <jurak@ma...>  20081129 09:07:38

Benjamin Kirk wrote: > I don't seem to see the attachment... But I have to ask: > > If the assembly function was taken from example 3 (2 does no assembly(?)), > are you constraining the hanging degrees of freedom? That is, is there a > line like > > // We have now built the element matrix and RHS vector in terms > // of the element degrees of freedom. However, it is possible > // that some of the element DOFs are constrained to enforce > // solution continuity, i.e. they are not really "free". We need > // to constrain those DOFs in terms of nonconstrained DOFs to > // ensure a continuous solution. The > // \p DofMap::constrain_element_matrix_and_vector() method does > // just that. > dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); > > (c.f. example 10) > > Ben > It is example 3; since the attachment didn't pass through i give the whole program below. I tried also constraining element matrices but it doesn't make any difference. Here is the program: Mladen #include <iostream> #include <algorithm> #include <cmath> #include <cassert> #include <cstdlib> #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "mesh_refinement.h" #include "linear_implicit_system.h" #include "equation_systems.h" #include "vector_value.h" #include "exact_solution.h" #include "fe.h" #include "quadrature_gauss.h" #include "dof_map.h" #include "sparse_matrix.h" #include "numeric_vector.h" #include "dense_matrix.h" #include "dense_vector.h" #include "elem.h" #include "string_to_enum.h" void assemble_poisson(EquationSystems& es, const std::string& system_name); Number es_ptr(const Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name) { const Real pi = libMesh::pi; double x=p(0), y=p(1); return std::sin( pi * x) * std::sin( pi * y); } Gradient ges_ptr(const Point &p, const Parameters ¶meters, const std::string &sys_name, const std::string &unknown_name) { const Real pi = libMesh::pi; double x=p(0), y=p(1); Gradient g; g(0) = pi * std::cos( pi * x) * std::sin( pi * y); g(1) = pi * std::sin( pi * x) * std::cos( pi * y); g(2) = 0.0; return g; } //============================================================== Real exact(double x, double y, double z) { const Real pi = libMesh::pi; return std::sin( pi * x) * std::sin( pi * y); } double exact(Point p) { const Real pi = libMesh::pi; double x=p(0), y=p(1); return std::sin( pi * x) * std::sin( pi * y); } //============================================================== int main (int argc, char** argv) { using std::cout; using std::endl; LibMeshInit init (argc, argv); Mesh mesh (2); ElemType el_type = Utility::string_to_enum<ElemType>("TRI3"); unsigned int n_refinement = 4; MeshTools::Generation::build_square (mesh, 10, 10, 0.0, 1.0, 0.0, 1.0, el_type); MeshRefinement mesh_refinement(mesh); EquationSystems es(mesh); LinearImplicitSystem & system = es.add_system<LinearImplicitSystem>("Poisson"); system.add_variable ("p", Utility::string_to_enum<Order> ("FIRST")); system.attach_assemble_function (assemble_poisson); es.parameters.set<Real>("linear solver tolerance") =1.0E15; es.init(); ExactSolution ex_sol(es); ex_sol.attach_exact_value(es_ptr); ex_sol.attach_exact_deriv(ges_ptr); ex_sol.extra_quadrature_order(5); Real l2_error=0.0, l2_error_old=0.0; Real linf_error=0.0, linf_error_old=0.0; Real h1_error=0.0, h1_error_old=0.0; std::ofstream out_file("results.txt"); // output file // out_file.open("results.txt"); out_file << " E R R O R C O N V E R G E N C E R A T E \n"; out_file << " H1 L2 Linf H1 L2 Linf \n"; out_file << "========== ========= ========== ========= ========= ==========\n"; int max_refinement=n_refinement; for(int refinement=0; refinement < max_refinement; ++ refinement) { system.solve(); ex_sol.compute_error("Poisson","p"); l2_error = ex_sol.l2_error("Poisson","p"); h1_error = ex_sol.h1_error("Poisson","p"); linf_error = ex_sol.l_inf_error("Poisson","p"); out_file << std::scientific << h1_error << " " << l2_error << " "<< linf_error<< " "; if(refinement >0) out_file << std::fixed << std::log2(h1_error_old/h1_error)<< " " << std::log2(l2_error_old/l2_error) << " " << std::log2(linf_error_old/linf_error); out_file <<std::endl; h1_error_old = h1_error; l2_error_old = l2_error; linf_error_old = linf_error; mesh_refinement.uniformly_refine(); es.reinit(); } out_file.close(); return 0; } void assemble_poisson(EquationSystems& es, const std::string& system_name) { 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); // 0 = var_number 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.local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.local_elements_end(); for ( ; el != end_el; ++el) { const Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); unsigned int n_dofs = dof_indices.size(); fe>reinit (elem); Ke.resize (n_dofs, n_dofs); Fe.resize (n_dofs); 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]*(dphi[i][qp]*dphi[j][qp]); for (unsigned int qp=0; qp<qrule.n_points(); qp++) { const Real x = q_point[qp](0); const Real y = q_point[qp](1); const Real z = q_point[qp](2); const Real eps = 1.e5; const Real uxx = (exact(xeps,y,z) + exact(x+eps,y,z) 2.*exact(x,y,z))/eps/eps; const Real uyy = (exact(x,yeps,z) + exact(x,y+eps,z) 2.*exact(x,y,z))/eps/eps; Real fxy; fxy =  (uxx + uyy); 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++) if (elem>neighbor(side) == NULL) { const Real penalty = 1.e20; 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); const Real value = exact(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]; } } } system.matrix>add_matrix (Ke, dof_indices); system.rhs>add_vector (Fe, dof_indices); } } 