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}
(73) 
_{Jun}
(51) 
_{Jul}
(12) 
_{Aug}
(69) 
_{Sep}
(24) 
_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 







1

2

3

4

5
(10) 
6
(6) 
7
(6) 
8
(3) 
9

10

11
(19) 
12
(6) 
13
(21) 
14
(1) 
15

16

17
(6) 
18
(7) 
19
(5) 
20
(3) 
21
(2) 
22

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

27
(2) 
28
(2) 
29
(11) 
30
(1) 






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: 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...>  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: Maxime Debon <maxime.debon@et...>  20081129 21:10:46

Thank you for your answer. My explanation was incomplete. My problem comes from a second stage library, a library made on libmesh.so or libmesh.a [ ex: libfoo ] First I tried to make a library with example ex0 as in the tutorial : http://www.adpgmbh.ch/cpp/gcc/create_lib.html The classic compilation with the makefile works perfectly but when creating the library, either I obtain "ar: lz: No such file or directory" for a static library libex0.a or "segmentation fault (core dumped)" after compilation with shared libex0.so The segmentation fault seems to occur at instruction : libMesh::init(argc, argv); May be some libraries are not completely linked Maxime 
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: 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: 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: Benjamin Kirk <benjamin.kirk@na...>  20081129 17:57:38

> On the website, the section "Linking With Your Application" explains > how to build an executable with : > > c++ o foo foo.C `libmeshconfig cxxflags include ldflags` > > This is fine on my installation. > > My question, surely trivial, is : How to build a libmesh static > library to call foo functions in external programs ? I'm not sure why you need a static library to do this instead of a shared one, but at any rate $ ./configure disableshared ... will do what you want. Ben 
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: Maxime Debon <maxime.debon@et...>  20081129 12:52:24

Hi, On the website, the section "Linking With Your Application" explains how to build an executable with : c++ o foo foo.C `libmeshconfig cxxflags include ldflags` This is fine on my installation. My question, surely trivial, is : How to build a libmesh static library to call foo functions in external programs ? Thanks, Maxime 
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); } } 
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...>  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: 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: Vasilis Vavourakis <vasvav@gm...>  20081127 16:15:50

in case the attached file is not present, i copiedpasted the c++ code below: ************************************** SOURCE CODE STARTS HERE *********************************** //================================================================================================= #include <iostream> #include <algorithm> #include <cmath> #include "dense_submatrix.h" #include "dense_subvector.h" #include "gmsh_io.h" #include "vtk_io.h" #include "xdr_io.h" #include "plane.h" #include "boundary_info.h" #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "gmv_io.h" #include "gnuplot_io.h" #include "linear_implicit_system.h" #include "equation_systems.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 "perf_log.h" #include "elem.h" #include "string_to_enum.h" #include "getpot.h" //================================================================================================= // Global parameters: // Material parameters: const Real young =1.0e0, poisson =0.3; const Real cc0 =young/(1.0poisson*poisson), cc1 =young*poisson/(1.0poisson*poisson), cc2 =young/2.0/(1.0+poisson); // Normal pressure loading: const Real pressure =1.0; // Penalty factor for imposing essential BCs: const Real penalty =1.0e10; //================================================================================================= void assemble_navier_cauchy (EquationSystems& es, const std::string& system_name); //================================================================================================= int main (int argc, char* argv[]) { LibMeshInit init(argc, argv); // Construct a 2D square mesh: Mesh mesh(2); MeshTools::Generation::build_square(mesh, 16, 16, 0.0, 1.0, 0.0, 1.0, TRI3); mesh.print_info(); // PerfLog perf_log("square plate under uniform tension"); // Create a linear implicit system for the NavierCauchy PDE of 2D elastostatics: EquationSystems system(mesh); system.add_system<LinearImplicitSystem>("NavierCauchy"); // Add displacement components: system.get_system("NavierCauchy").add_variable("u0", FIRST, LAGRANGE); system.get_system("NavierCauchy").add_variable("u1", FIRST, LAGRANGE); // Attach assembly function: system.get_system("NavierCauchy").attach_assemble_function(assemble_navier_cauchy); system.init(); system.print_info(); // Solve linear system: system.get_system("NavierCauchy").solve(); // Print the results in Gmsh format: GmshIO(mesh).write_equation_systems("Kirsch2dElastostatics.msh", system); // Exit function normally: return 0; } //================================================================================================= void assemble_navier_cauchy (EquationSystems& es, const std::string& system_name) { libmesh_assert(system_name=="NavierCauchy"); PerfLog perf_log("matrix assembly"); const MeshBase& mesh =es.get_mesh(); const unsigned int dim = mesh.mesh_dimension(); LinearImplicitSystem& system =es.get_system<LinearImplicitSystem>("NavierCauchy"); QGauss quad_volu(dim, FIFTH); QGauss quad_surf(dim1, FIFTH); const unsigned int u0_var =system.variable_number("u0"); const unsigned int u1_var =system.variable_number("u1"); FEType fe_u_type =system.variable_type(u0_var); AutoPtr<FEBase> fe_u_volu(FEBase::build(dim, fe_u_type)); // FEM for volume integration. AutoPtr<FEBase> fe_u_surf(FEBase::build(dim, fe_u_type)); // FEM for surface integration. fe_u_volu>attach_quadrature_rule(&quad_volu); fe_u_surf>attach_quadrature_rule(&quad_surf); const std::vector<Real>& JxW_volu =fe_u_volu>get_JxW(); const std::vector<std::vector<Real> >& phi_u_volu =fe_u_volu>get_phi(); const std::vector<std::vector<RealGradient> >& dphi_u_volu =fe_u_volu>get_dphi(); const DofMap& dof_map =system.get_dof_map(); DenseMatrix<Number> Ke; DenseVector<Number> Fe; DenseSubMatrix<Number> K00(Ke), K01(Ke), K10(Ke), K11(Ke); DenseSubVector<Number> F0(Fe), F1(Fe); std::vector<unsigned int> dof_indices; std::vector<unsigned int> dof_indices_u0; std::vector<unsigned int> dof_indices_u1; for ( MeshBase::const_element_iterator citer_el=mesh.local_elements_begin(); citer_el!=mesh.local_elements_end(); citer_el++) { const Elem* elem =*citer_el; // Initialize corresponding arrays: perf_log.push("initialization"); dof_map.dof_indices(elem, dof_indices); dof_map.dof_indices(elem, dof_indices_u0, u0_var); dof_map.dof_indices(elem, dof_indices_u1, u1_var); const unsigned int n_dofs =dof_indices.size(); const unsigned int n_u0_dofs =dof_indices_u0.size(); const unsigned int n_u1_dofs =dof_indices_u1.size(); Ke.resize(n_dofs, n_dofs); Fe.resize(n_dofs); K00.reposition(u0_var*n_u0_dofs, u0_var*n_u0_dofs, n_u0_dofs, n_u0_dofs); K01.reposition(u0_var*n_u0_dofs, u1_var*n_u1_dofs, n_u0_dofs, n_u1_dofs); K10.reposition(u1_var*n_u1_dofs, u0_var*n_u0_dofs, n_u1_dofs, n_u0_dofs); K11.reposition(u1_var*n_u1_dofs, u1_var*n_u1_dofs, n_u1_dofs, n_u1_dofs); F0.reposition(u0_var*n_u0_dofs, n_u0_dofs); F1.reposition(u1_var*n_u1_dofs, n_u1_dofs); fe_u_volu>reinit(elem); perf_log.pop("initialization"); // Stiffness matrix: perf_log.push("stiffness matrix"); for (unsigned int qp=0; qp<quad_volu.n_points(); qp++) { for (unsigned int i=0; i<n_u0_dofs; i++) for (unsigned int j=0; j<n_u0_dofs; j++) K00(i,j) += JxW_volu[qp]*(dphi_u_volu[i][qp](0)*cc0*dphi_u_volu[j][qp](0)+ dphi_u_volu[i][qp](1)*cc2*dphi_u_volu[j][qp](1)); for (unsigned int i=0; i<n_u0_dofs; i++) for (unsigned int j=0; j<n_u1_dofs; j++) K01(i,j) += JxW_volu[qp]*(dphi_u_volu[i][qp](0)*cc1*dphi_u_volu[j][qp](1)+ dphi_u_volu[i][qp](1)*cc2*dphi_u_volu[j][qp](0)); for (unsigned int i=0; i<n_u1_dofs; i++) for (unsigned int j=0; j<n_u0_dofs; j++) K10(i,j) += JxW_volu[qp]*(dphi_u_volu[i][qp](1)*cc1*dphi_u_volu[j][qp](0)+ dphi_u_volu[i][qp](0)*cc2*dphi_u_volu[j][qp](1)); for (unsigned int i=0; i<n_u1_dofs; i++) for (unsigned int j=0; j<n_u1_dofs; j++) K11(i,j) += JxW_volu[qp]*(dphi_u_volu[i][qp](1)*cc0*dphi_u_volu[j][qp](1)+ dphi_u_volu[i][qp](0)*cc2*dphi_u_volu[j][qp](0)); } perf_log.pop("stiffness matrix"); // Bodyforce vector: perf_log.push("rhs vector"); // Nobody forces are present, that's is why the following loop is commented. /* for (unsigned int qp=0; qp<quad_volu.n_points(); qp++) { for (unsigned int i=0; i<n_u0_dofs; i++) F0(i) += JxW_volu[qp]*phi_u_volu[i][qp]*0.0; for (unsigned int i=0; i<n_u1_dofs; i++) F1(i) += JxW_volu[qp]*phi_u_volu[i][qp]*0.0; }*/ perf_log.pop("rhs vector"); // Boundary conditions: perf_log.push("insert BCs"); for (unsigned int s=0; s<elem>n_sides(); s++) { if (elem>neighbor(s)!=NULL) continue; const short int bottom =0, // LIBMESH square plate under uniform tension in 2D right =1, top =2, left =3; // Get the boundary ID for side 's'. These are set internally by build_square(): short int bc_id =mesh.boundary_info>boundary_id(elem,s); if (bc_id==BoundaryInfo::invalid_id) libmesh_error(); AutoPtr<Elem> side(elem>build_side(s)); const std::vector<Real>& JxW_surf =fe_u_surf>get_JxW(); const std::vector<std::vector<Real> >& phi_u_surf =fe_u_surf>get_phi(); const std::vector<Point>& eta =fe_u_surf>get_normals(); fe_u_surf>reinit(elem, s); if (bc_id==left) { // zero u0 Dirichlet BC: const Real u0_value =0.0; for (unsigned int ns=0; ns<side>n_nodes(); ns++) { for (unsigned int ne=0; ne<elem>n_nodes(); ne++) { if (elem>node(ne) == side>node(ns)) { K00(ne,ne) += penalty; F0(ne) += penalty * u0_value; } } } } else if (bc_id==bottom) { // zero u1 Dirichlet BC: const Real u1_value =0.0; for (unsigned int ns=0; ns<side>n_nodes(); ns++) { for (unsigned int ne=0; ne<elem>n_nodes(); ne++) { if (elem>node(ne) == side>node(ns)) { K11(ne,ne) += penalty; F1(ne) += penalty * u1_value; } } } } else if (bc_id==right) { // traction (normal pressure) load Neumann BC: for (unsigned int qp=0; qp<quad_surf.n_points(); qp++) { Real p_value =0.0; for (unsigned int j=0; j<phi_u_surf.size(); j++) { p_value += phi_u_surf[j][qp] * pressure; } for (unsigned int i=0; i<n_u0_dofs; i++) F0(i) += JxW_surf[qp] * phi_u_surf[i][qp] * eta[qp](0) * p_value; for (unsigned int i=0; i<n_u1_dofs; i++) F1(i) += JxW_surf[qp] * phi_u_surf[i][qp] * eta[qp](1) * p_value; } } } perf_log.pop("insert BCs"); // Insert array and vector in general system: perf_log.push("matrix insertion"); system.matrix>add_matrix(Ke, dof_indices); system.rhs>add_vector(Fe, dof_indices); perf_log.pop("matrix insertion"); } } //================================================================================================= ************************************** SOURCE CODE ENDS HERE *********************************** :) 
From: Vasilis Vavourakis <vasvav@gm...>  20081127 15:56:10

Hello all :) I attached a source file that implements a displacementbased formulation for steadystate 2D elasticity (very simple example for someone to understand it) and this piece of code works FINE. For my practice i also implemented a mixed formulation (displacementsstresses) and i did manage to evaluate them correctly. If there is anyone that may need it i can post it as well!!! My problem in the displacement formulation (see source code file) is that after the numerical evaluation of displacements i couldn't find a way to calculate nodal stresses in the postprocess. Although the same question it was posted both from me and others in the past ( see http://sourceforge.net/mailarchive/forum.php?thread_name=.139.91.194.170.1226927623.squirrel%40www.mead.upatras.gr&forum_name=libmeshusers), i still cannot implement the ExplicitSystem class so as to do my postprocess! I red over & over again the documentation for this class but it was really hard for me to figure out a way to do it! :( Unfortunately i got very confused due to the fact that i'm bery new to LibMesh as well as its "vast" size of documentation for classes... I'd really appreciate if someone can help me out with this... Apart from useful advices, any handy source code is really really welcomed :) many thanks in advance ~:^) 
From: Yujie <recrusader@gm...>  20081125 00:02:39

Dear Roy: I am checking the codes about the output file. If using the method you provide, one can't simultaneously output the variable and vector in a file, right? Because the vectors are the copies of the solutions in my application, I am considering to refer to build_variable_names() and build_solution_vector() to write two functions (build_vector_names() and build_vector_vector()) to output the vector values. I feel it should be simple. could you give me some comments? thanks a lot. Regards, Yujie On Mon, Nov 24, 2008 at 3:51 PM, Roy Stogner <roystgnr@...>wrote: > On Mon, 24 Nov 2008, Yujie wrote: > > In addition, how to write the vectors added by add_vector() to the output >> files (such as gmv, tecplot formats)? >> > > Swap the solution and the added vector, write the solution, swap back. >  > Roy > 
From: Roy Stogner <roystgnr@ic...>  20081124 20:09:14

On Mon, 24 Nov 2008, Yujie wrote: > I am sorry to further bother you. Howeve, I have read the patch you provided > to me. I think I didn't demonstrate my problem detailedly. You demonstrated your problem, correctly; you just didn't understand the "diff" format of the patch I provided. The lines with "" in front of them are lines from the existing code which the patch *removes*. See the "patch" command for details. If you can't get it to work, I'll send you a complete .C file when I'm done traveling.  Roy 
From: Yujie <recrusader@gm...>  20081124 17:42:44

Dear Roy: I am sorry to further bother you. Howeve, I have read the patch you provided to me. I think I didn't demonstrate my problem detailedly. I need to dynamically add vectors after System::init(). However, regarding the patch, I think it only is for before System::init(), right? Because if System::init() is implemented, the program will generate errors in the following codes to my knowledge "  // We can only add new vectors before initializing...  if (!_can_add_vectors)  {  std::cerr << "ERROR: Too late. Cannot add vectors to the system after initialization"  << std::endl  << " any more. You should have done this earlier."  << std::endl;  libmesh_error();  } " thanks a lot. Regards, Yujie On Fri, Nov 21, 2008 at 1:05 PM, Roy Stogner <roystgnr@...>wrote: > > On Thu, 20 Nov 2008, Yujie wrote: > > I notice that if i need to add vectors in LinearImplicitSystem using >> System::add_vector() it must be done before System::init(). However, I >> can't >> confirm the vector number at the beginning. How to dynamically add >> vectors? >> > > Try the following patch, and let me know if it works; if so I'll commit it > to SVN. > > > Index: src/solvers/system.C > =================================================================== >  src/solvers/system.C (revision 3160) > +++ src/solvers/system.C (working copy) > @@ 493,21 +493,17 @@ > return *(_vectors[vec_name]); > } > >  // We can only add new vectors before initializing... >  if (!_can_add_vectors) >  { >  std::cerr << "ERROR: Too late. Cannot add vectors to the system > after initialization" >  << std::endl >  << " any more. You should have done this earlier." >  << std::endl; >  libmesh_error(); >  } >  >  // Otherwise build the vector and return it. > + // Otherwise build the vector > NumericVector<Number>* buf = NumericVector<Number>::build().release(); > _vectors.insert (std::make_pair (vec_name, buf)); > _vector_projections.insert (std::make_pair (vec_name, projections)); > > + // Initialize it if necessary > + if (!_can_add_vectors) > + { > + buf>init (this>n_dofs(), this>n_local_dofs()); > + } > + > return *buf; > } > > > In addition, In the same LinearImplicitSystem, the vectors added by >> add_vector() have similar characteristics with the variable added by >> add_variable(), such as parallel distribution, same Dofs? thanks a lot. >> > > Every vector added to a system has dofs for every variable in that system. > All > vectors have the same indexing and the same parallel distribution. (In > particular they're all parallelized like solution, and if you want ghost > dof > coefficients you'll need to create a local vector to hold them) >  > Roy > 
From: John Peterson <jwpeterson@gm...>  20081124 15:20:01

On Mon, Nov 24, 2008 at 7:37 AM, Tim Kroeger <tim.kroeger@...> wrote: > Dear all, > > On Mon, 24 Nov 2008, Tim Kroeger wrote: > >> src/geom/elem_refinement.C: In member function 'unsigned int >> Elem::_cast_node_address_to_unsigned_int(unsigned int)': >> src/geom/elem_refinement.C:160: error: cast from 'Node*' to 'unsigned int' >> loses precision > > I now manually modified the following three things, which seems to solve my > problem: > > > 1. I added > > libmesh_CXXFLAGS += m32 > libmesh_CFLAGS += m32 > > at the end of Make.common. > > > 2. In the Tecplot section of Make.common, I adjusted the path to tecio.a > from "x86_64unknownlinuxgnu" to "i686pclinuxgnu". > > > 3. In contrib/tetgen/Makefile, I again added > > libmesh_CXXFLAGS += m32 > > after line 10. > > > It seems to work now, but anyway, I would appreciate if these things > happened automatically. I think adding m32, or whatever this option is called on other compilers, to libmesh_CXXFLAGS should be the way to go. Then tetgen et al. should use the libmesh_CXXFLAGS as well. The latest Intel compilers also appear to support m32. Under PGI it seems to get a little more complicated, with "tp barcelona tp barcelona32" required on the Barcelona chips. It's probably worthwhile writing our own autoconf test which sets the correct flag based on the compiler. Unless autoconf already has this...I'll take a look.  John 
From: Tim Kroeger <tim.kroeger@ce...>  20081124 13:37:54

Dear all, On Mon, 24 Nov 2008, Tim Kroeger wrote: > src/geom/elem_refinement.C: In member function ‘unsigned int > Elem::_cast_node_address_to_unsigned_int(unsigned int)’: > src/geom/elem_refinement.C:160: error: cast from ‘Node*’ to ‘unsigned int’ > loses precision I now manually modified the following three things, which seems to solve my problem: 1. I added libmesh_CXXFLAGS += m32 libmesh_CFLAGS += m32 at the end of Make.common. 2. In the Tecplot section of Make.common, I adjusted the path to tecio.a from "x86_64unknownlinuxgnu" to "i686pclinuxgnu". 3. In contrib/tetgen/Makefile, I again added libmesh_CXXFLAGS += m32 after line 10. It seems to work now, but anyway, I would appreciate if these things happened automatically. Best Regards, Tim  Dr. Tim Kroeger Phone +494212187710 tim.kroeger@..., tim.kroeger@... Fax +494212184236 MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany Amtsgericht Bremen HRB 16222 Geschaeftsfuehrer: Prof. Dr. H.O. Peitgen 
From: Tim Kroeger <tim.kroeger@ce...>  20081124 11:41:43

Dear libMesh developers, For certain reasons, I need to compile a 32bit libMesh on a 64bit machine. What I did was: 1. I compiled mpich2 in 32 bit mode, i.e. by saying export CFLAGS=m32 export LDFLAGS=m32 export CPPFLAGS=m32 export FC=gfortran export F90=gfortran export FFLAGS=m32 export F90FLAGS=m32 prior to configuring mpich2 (without problems). 2. I compiled PETSc2.3.3p11 against that mpich2 (without problems). 3. I set the above variables again, before calling libMesh's "./configure". 4. I said "make" in the libMesh directory. What I get is: src/geom/elem_refinement.C: In member function ‘unsigned int Elem::_cast_node_address_to_unsigned_int(unsigned int)’: src/geom/elem_refinement.C:160: error: cast from ‘Node*’ to ‘unsigned int’ loses precision It seems to me that libMesh configured for 32bit mode (in particular determining LIBMESH_SIZEOF_INT and LIBMESH_SIZEOF_VOID_P to be equal), but then compiled without the "m32" option. This is confirmed by saying "make n" and observing that the command to be executed does not contain "m32". What did I do wrong? Best Regards, Tim  Dr. Tim Kroeger Phone +494212187710 tim.kroeger@..., tim.kroeger@... Fax +494212184236 MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany Amtsgericht Bremen HRB 16222 Geschaeftsfuehrer: Prof. Dr. H.O. Peitgen 
From: Roy Stogner <roystgnr@ic...>  20081123 19:06:19

On Sat, 22 Nov 2008, Yujie wrote: > I don't sorry that I don't understand the second problem. To my knowledge, > if I use add_variable() to add the variables, these variables will add dofs, > right? Right. > If I use add_vector() to add the vectors, I don't specify the > purpose of the vectors and they don't add dofs, right? Right; each new vector has a coefficient corresponding to each existing dof. > Now, I want to copy the solution (variables) to the vectors(added by > add_vector()), why do I need to create the local vector? You don't need to create anything to copy the solution to the new vectors. But if you want to copy the current_local_solution to the new vectors, it won't work, because the new vectors only have storage for local dof coefficients, not for ghost dof coefficients.  Roy 
From: Yujie <recrusader@gm...>  20081123 06:44:24

Dear Roy: I don't sorry that I don't understand the second problem. To my knowledge, if I use add_variable() to add the variables, these variables will add dofs, right? If I use add_vector() to add the vectors, I don't specify the purpose of the vectors and they don't add dofs, right? Now, I want to copy the solution (variables) to the vectors(added by add_vector()), why do I need to create the local vector? why not just use the vectors added by add_vector()？ it doesn't work? thanks a lot. Regards, Yujie > > > > In addition, In the same LinearImplicitSystem, the vectors added by >> add_vector() have similar characteristics with the variable added by >> add_variable(), such as parallel distribution, same Dofs? thanks a lot. >> > > Every vector added to a system has dofs for every variable in that system. > All > vectors have the same indexing and the same parallel distribution. (In > particular they're all parallelized like solution, and if you want ghost > dof > coefficients you'll need to create a local vector to hold them) >  > Roy > 
From: Roy Stogner <roystgnr@ic...>  20081121 21:06:00

On Thu, 20 Nov 2008, Yujie wrote: > I notice that if i need to add vectors in LinearImplicitSystem using > System::add_vector() it must be done before System::init(). However, I can't > confirm the vector number at the beginning. How to dynamically add vectors? Try the following patch, and let me know if it works; if so I'll commit it to SVN. Index: src/solvers/system.C ===================================================================  src/solvers/system.C (revision 3160) +++ src/solvers/system.C (working copy) @@ 493,21 +493,17 @@ return *(_vectors[vec_name]); }  // We can only add new vectors before initializing...  if (!_can_add_vectors)  {  std::cerr << "ERROR: Too late. Cannot add vectors to the system after initialization"  << std::endl  << " any more. You should have done this earlier."  << std::endl;  libmesh_error();  }   // Otherwise build the vector and return it. + // Otherwise build the vector NumericVector<Number>* buf = NumericVector<Number>::build().release(); _vectors.insert (std::make_pair (vec_name, buf)); _vector_projections.insert (std::make_pair (vec_name, projections)); + // Initialize it if necessary + if (!_can_add_vectors) + { + buf>init (this>n_dofs(), this>n_local_dofs()); + } + return *buf; } > In addition, In the same LinearImplicitSystem, the vectors added by > add_vector() have similar characteristics with the variable added by > add_variable(), such as parallel distribution, same Dofs? thanks a lot. Every vector added to a system has dofs for every variable in that system. All vectors have the same indexing and the same parallel distribution. (In particular they're all parallelized like solution, and if you want ghost dof coefficients you'll need to create a local vector to hold them)  Roy 