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}
(27) 
_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 




1
(1) 
2
(4) 
3

4

5

6

7
(1) 
8

9

10

11

12

13
(2) 
14
(1) 
15
(6) 
16
(6) 
17

18

19
(2) 
20
(1) 
21
(2) 
22
(3) 
23
(8) 
24
(4) 
25

26
(2) 
27
(2) 
28
(1) 




From: Wout Ruijter <woutruijter@gm...>  20060228 16:22:38

ShengliXu // you forgot the all important line: dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices); before system.matrix>add_matrix (Ke, dof_indices); system.rhs>add_vector (Fe, dof_indices); So all the other explanations I came up with for the program not working are farfetched or nonsense. You can choose the solver on the command line, for example: time /usr/bin/mpirun np 1 $(target) ksp_type gmres pc_type jacobi >> log= .txt Now, one more time: If you apply a strain in xdirection the bounds are going to move apart, this CONFLICTS with the boundary conditions you just put up (as you're telling them they have the same displacement). As you will see, if you couple all dofs on all opposing sides as you propose you will only get some noise as an answer regardless of the loading. Read my earlier post on how to deal with that (by adding dofs to the system in a dummy element that represent a displacement of a whole side, keeping the sides in the same shape). Secondly, the code you put in to deal with corners does not do anything, as I said, the corners are in the sides and are dealt with automatically (because "couples" is a std::map and keys are unique in those). Cheers W 
From: Wout Ruijter <woutruijter@gm...>  20060227 11:13:41

ShengliXu First of all, you're trying to solve something that doesn't make sense (to me at least), you try to compress a square in xdirection and at the same time demand the left and right bound to have the same xdisplacement. You can't get a result that shows the effects of your periodic code because the penalty is 1e10 and therefore the DofConstraints (of order 1) are overshadowed. Always print out the convergence history of your solver, if you have trouble converging usually something like that is happening. If you want the bounds to have the same shape but allow them to move closer together as a whole, you have to add a "ghost" DOF to the system which represents this movement, and add it to all of the dofconstraintrows that you are adding for a certain dimension/dof combination. I needed this before and implemented it by having a loose element added to the system, which will give you a bunch of dofs to use for such purposes (multiple dimension periodicity, for example). Make sure you don't assemble that element though, just set the diagonal terms to nonzero for the dofs you don't use. Secondly, there must be some bug in your code as it only finds one set of opposite dofs, I wrote a periodic BC code in the file that I find works. It might be a matter of taste, but I think the preferred way is still to use element based looping and dof querying as opposed to node based (the smart people on the list might correct me on that ;). Cheers Wout 
From: Steffen Petersen <steffen.petersen@tu...>  20060227 08:33:35

I have just committed an example 17 including a generalized eigenvalue problem. Currently, the largest eigenvalues in the spectrum are computed. The code can be easily changed in order to compute the smallest values (or use the slepc comman line options). However, the Lanczos solver slepc offers with this option appears to be somewhat unstable when used to compute the smallest values (perhaps this is due to the zero eigenmode in this example). An option to work around this problem is to configure and install slepc with arpack and parpack (the arpack solver works fine for example 17). Steffen David Xu schrieb: > Hi Steffen and Ondrej, > > I read your previous message that you've extended the slepc interface > to support generalized eigen problems. I'd like to get the details > from you. I'm new to libmesh and I'm wondering if you can post some > new tutorial examples on how to solve eigen problems using libmesh > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > help. > > Ondrej, I read about your message about your implementation of libmesh > and pysparse's jacobidavidson solver. That sounds really good. Could > you please give me the details too? I really appreciate it. > > David 
From: Shengli Xu <shengli.xu.xu@gm...>  20060226 17:01:14

Micheal, Wout, Ben I am testing a code for the problem of 2D plane stress. I want to impose periodic boundary condition just like u1=3Du2. I mainly consider the post of Ben and Micheal. I don't konw where to call the subroutine for periodic bc. The result considering periodic bc is the same as not considering periodic bc. I don't konw how to do. Maybe The part of code considering periodic bc is wrong. I have mailed to Wout for this problem. Maybe I need to add some code to system.h and system.C. The source code is attached. The example is: 1x1 domain, 10x10 mesh, Quad4, fixed bc on x=3D0, p=3D1 on x=3D1, consider periodic bc of ux with node (0.5,0) and (0.8,0). // 2D elastic code /* The Next Great Finite Element Library. */ /* Copyright (C) 2003 Benjamin S. Kirk */ /* This library is free software; you can redistribute it and/or */ /* modify it under the terms of the GNU Lesser General Public */ /* License as published by the Free Software Foundation; either */ /* version 2.1 of the License, or (at your option) any later version. */ /* This library is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */ /* Lesser General Public License for more details. */ /* You should have received a copy of the GNU Lesser General Public */ /* License along with this library; if not, write to the Free Software */ /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 021111307 US= A */ #include <iostream> #include <fstream> #include <algorithm> #include <math.h> #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "gmv_io.h" #include "linear_implicit_system.h" #include "equation_systems.h" #include "fe.h" #include "quadrature_gauss.h" #include "sparse_matrix.h" #include "numeric_vector.h" #include "dense_matrix.h" #include "dense_vector.h" // add by shengli #include "elem.h" #include "dof_map.h" #include "perf_log.h" // for assembling the element matrix and vector on a componentbycomponent basis #include "dense_submatrix.h" //#include "dense_sub void assemble_elastic(EquationSystems& es, const std::string& system_name); // impose periodic boundary condition void apply_periodic_bc(EquationSystems &es, const std::string &system_name)= ; int main (int argc, char** argv) { libMesh::init (argc, argv); { std::cout << "Running " << argv[0]; for (int i=3D1; i<argc; i++) std::cout << " " << argv[i]; std::cout << std::endl << std::endl; PerfLog perf("Main Program"); perf.start_event("program init"); Mesh mesh (2); // Use the MeshTools::Generation mesh generator to create a uniform // grid on the square [0,1]^2. MeshTools::Generation::build_square(mesh,10,10,0.,1.,0.,1.,QUAD4); // MeshData mesh_data(mesh); //mesh_data.activate(); // mesh_data.enable_compatibility_mode(); // mesh.read("in.xda",&mesh_data); // mesh_data.read("data.xta"); mesh.print_info(); // EquationSystems equation_systems (mesh,&mesh_data); EquationSystems equation_systems (mesh); LinearImplicitSystem &system =3D equation_systems.add_system<LinearImplicitSystem>("elastic"); // equation_systems.add_system<LinearImplicitSystem> ("elastic"); //equation_systems("elastic").add_variable("u", SECOND); // equation_systems("elastic").add_variable("v", SECOND); // equation_systems("elastic").attach_assemble_function (assemble_elastic); system.add_variable("u",FIRST); system.add_variable("v",FIRST); system.attach_assemble_function(assemble_elastic); // define equation_systems variable equation_systems.parameters.set<Real>("mu")=3D0.3; equation_systems.parameters.set<Real>("E")=3D20000; equation_systems.init(); equation_systems.print_info(); perf.stop_event("program init"); perf.start_event("solve"); equation_systems("elastic").solve(); perf.stop_event("solve"); perf.start_event("saving"); GMVIO (mesh).write_equation_systems ("out.gmv", equation_systems); perf.stop_event("saving"); } return libMesh::close(); } void assemble_elastic(EquationSystems& es, const std::string& system_name) { assert (system_name =3D=3D "elastic"); PerfLog perf("Matrix Assembly"); const Mesh& mesh =3D es.get_mesh(); // const MeshData& mesh_data =3D es.get_mesh_data(); const unsigned int dim =3D mesh.mesh_dimension(); LinearImplicitSystem& system=3Des.get_system <LinearImplicitSystem>("elastic"); const Real mu =3D es.parameters.get<Real>("mu"); const Real E =3D es.parameters.get<Real>("E"); // Numeric ids corresponding to each variable in the system const unsigned int u_var=3Dsystem.variable_number("u"); const unsigned int v_var=3Dsystem.variable_number("v"); const DofMap& dof_map =3D system.get_dof_map(); FEType fe_type =3D dof_map.variable_type(0); AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); // QGauss qrule (dim, FIFTH); QGauss qrule(dim,fe_type.default_quadrature_order()); 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 =3D fe>get_JxW(); // const std::vector<Point>& q_point =3D fe>get_xyz(); // const std::vector<std::vector<Real> >& phi =3D fe>get_phi(); const std::vector<std::vector<RealGradient> >& dphi =3D fe>get_dphi(); DenseMatrix<Number> Ke; DenseVector<Number> Fe; DenseSubMatrix<Number> Kuu(Ke),Kuv(Ke), Kvu(Ke),Kvv(Ke); DenseSubVector<Number> Fu(Fe),Fv(Fe); std::vector<unsigned int> dof_indices; std::vector<unsigned int> dof_indices_u; std::vector<unsigned int> dof_indices_v; MeshBase::const_element_iterator el =3D mesh.elements_begin(); const MeshBase::const_element_iterator end_el =3D mesh.elements_end(); for ( ; el !=3D end_el ; ++el){ perf.start_event("elem init"); const Elem* elem =3D *el; dof_map.dof_indices (elem, dof_indices); dof_map.dof_indices(elem,dof_indices_u,u_var); dof_map.dof_indices(elem,dof_indices_v,v_var); // const unsigned int n_dofs=3Ddof_indices.size(); const unsigned int n_u_dofs=3Ddof_indices_u.size(); const unsigned int n_v_dofs=3Ddof_indices_v.size(); // Compute the elementspecific data for the current element. fe>reinit (elem); Ke.resize (n_dofs,n_dofs); Fe.resize (n_dofs); // DenseSubMatrix.repositon() unsigned int ioff,joff; ioff=3D0; joff=3D0; Kuu.reposition(ioff,joff,n_u_dofs,n_u_dofs); joff +=3D n_u_dofs; Kuv.reposition(ioff,joff,n_u_dofs,n_v_dofs); ioff +=3D n_u_dofs; joff =3D 0; Kvu.reposition(ioff,joff,n_v_dofs,n_u_dofs); joff +=3D n_u_dofs; Kvv.reposition(ioff,joff,n_v_dofs,n_v_dofs); ioff=3D0; Fu.reposition(ioff,n_u_dofs); ioff +=3D n_u_dofs; Fv.reposition(ioff,n_v_dofs); perf.stop_event("elem init"); perf.start_event("Ke"); for (unsigned int qp=3D0; qp<qrule.n_points(); qp++){ // Assemble Kuu for (unsigned int i=3D0; i<n_u_dofs;i++){ for (unsigned int j=3D0; j<n_u_dofs;j++){ Kuu(i,j) +=3D JxW[qp]*(dphi[i][qp](0)*dphi[j][qp](0)/(1mu*mu)+dphi[i][qp](1)*dphi[j][qp]= (1)/2./(1.+mu))*E; } } // Assemble Kuv for(unsigned int i=3D0;i<n_u_dofs;i++){ for(unsigned int j=3D0;j<n_v_dofs;j++){ Kuv(i,j) +=3D JxW[qp]*(dphi[i][qp](0)*dphi[j][qp](1)*mu/(1mu*mu)+dphi[i][qp](1)*dphi[j][= qp](0)/2./(1.+mu))*E; } } // Assemble Kvu for(unsigned int i=3D0;i<n_v_dofs;i++){ for(unsigned int j=3D0;j<n_u_dofs;j++){ Kvu(i,j) +=3D JxW[qp]*(dphi[i][qp](1)*dphi[j][qp](0)*mu/(1mu*mu)+dphi[i][qp](0)*dphi[j][= qp](1)/2./(1.+mu))*E; } } // Assemble Kvv for(unsigned int i=3D0;i<n_v_dofs;i++){ for(unsigned int j=3D0;j<n_v_dofs;j++){ Kvv(i,j) +=3D JxW[qp]*(dphi[i][qp](1)*dphi[j][qp](1)/(1mu*mu)+dphi[i][qp](0)*dphi[j][qp]= (0)/2./(1.+mu))*E; } } } perf.stop_event("Ke"); // output Ke std::cout<<"output stiffness matrix"<<std::endl; for(unsigned int i=3D0;i<n_dofs;i++){ for(unsigned int j=3D0;j<n_dofs;j++){ std::cout<<"("<<i<<","<<j<<")"<<Ke(i,j)<<", "; } std::cout<<std::endl; } perf.start_event("matrix insertion"); system.matrix>add_matrix (Ke, dof_indices); system.rhs>add_vector (Fe, dof_indices); perf.stop_event("matrix insertion"); } // assemble load and impose bc { assert(system_name =3D=3D "elastic"); LinearImplicitSystem &system =3D es.get_system <LinearImplicitSystem>("elastic"); // Get writable references to the overall matrix and vector SparseMatrix<Number> &matrix =3D *system.matrix; NumericVector<Number> &rhs =3D *system.rhs; // Get a const reference to the mesh object const Mesh &mesh =3D es.get_mesh(); // Apply load std::cout<<"Assemble Load"<<std::endl; unsigned int n_nodes =3D mesh.n_nodes(); // fy=3D100 at node (1,0) for(unsigned int n_cnt=3D0;n_cnt<n_nodes;n_cnt++){ const Node &curr_node=3Dmesh.node(n_cnt); const Real x=3Dcurr_node(0); const Real y=3Dcurr_node(1); /* if(fabs(y0.0)<TOLERANCE && fabs(x1.)<TOLERANCE){ unsigned int ndof =3D curr_node.dof_number(0,1,0); rhs.set(ndof,100.0); // fy=3D100 break; } */ //apply pressure=3D1 on x=3D1 if(fabs(y0.0)<TOLERANCE && fabs(x1.)<TOLERANCE){ unsigned int n =3D curr_node.dof_number(0,0,0); rhs.set(n,0.05); } if(fabs(y1.)<TOLERANCE && fabs(x1.)<TOLERANCE){ unsigned int n =3D curr_node.dof_number(0,0,0); rhs.set(n,0.05); } if(0.01<y && y< 0.99 && fabs(x1.)<TOLERANCE){ unsigned int n =3D curr_node.dof_number(0,0,0); rhs.set(n,0.1); } } // fixed bc at left side for(unsigned int n_cnt=3D0;n_cnt<n_nodes;n_cnt++){ const Node &curr_node=3Dmesh.node(n_cnt); const Real x=3Dcurr_node(0); const Real y=3Dcurr_node(1); const Real penalty =3D 1.e10; if(fabs(x0.0)<TOLERANCE){ unsigned int dofu =3D curr_node.dof_number(0,0,0); const Real u=3D0.0; rhs.set(dofu,u*penalty); std::cout<<"dofu=3D"<<dofu<<std::endl<<"matrix"<<matrix(dofu,dofu)<<std::en= dl; // matrix.set(dofu,dofu,penalty); matrix.add(dofu,dofu,penalty); /* unsigned int dofv =3D curr_node.dof_number(0,1,0); const Real v=3D0.0; rhs.set(dofv,v*penalty); matrix.add(dofv,dofv,penalty); */ if(fabs(y0.0)<TOLERANCE){ unsigned int n =3D curr_node.dof_number(0,1,0); const Real v=3D0.0; rhs.set(n,v*penalty); matrix.add(n,n,penalty); } } } } { apply_periodic_bc(es,system_name); } } void apply_periodic_bc(EquationSystems &es, const std::string &system_name)= { assert(system_name=3D=3D"elastic"); PerfLog perf_log("Periodic bc. Assembly",false); const Mesh &mesh=3Des.get_mesh(); LinearImplicitSystem &system=3Des.get_system <LinearImplicitSystem>("elastic"); DofMap &dof_map=3Dsystem.get_dof_map(); // impose ux(0.5,0.0)=3Dux(0.8,0.0) for(unsigned int n=3D0;n<mesh.n_nodes();n++){ const Node &node1=3Dmesh.node(n); const double x1=3Dnode1(0); if((fabs(node1(0)0.5)<TOLERANCE) && fabs(node1(1))<TOLERANCE){ for(unsigned int n1=3D0;n1<mesh.n_nodes();n1++){ const Node &node2=3Dmesh.node(n1); const double x2=3Dnode2(0); std::cout<<"x1=3D"<<x1<<",x2=3D"<<x2<<std::endl; if(fabs(node2(0)0.8)<TOLERANCE && fabs(node2(1))<TOLERANCE){ unsigned int d1=3Dnode1.dof_number(0,0,0); unsigned int d2=3Dnode2.dof_number(0,0,0); std::cout<<"n_left=3D"<<n<<",n_right=3D"<<n1<<std::endl; DofConstraintRow constraint; constraint[d1]=3D1; dof_map.add_constraint_row(d2,constraint); break; } } } } } 
From: Shengli Xu <shengli.xu.xu@gm...>  20060226 02:25:21

Micheal, Wout, Ben I am testing a code for the problem of 2D plane stress. I want to impose periodic boundary condition just like u1=3Du2. I mainly consider the post of Ben and Micheal. I don't konw where to call the subroutine for periodic bc. The result considering periodic bc is the same as not considering periodic bc. I don't konw how to do. Maybe The part of code considering periodic bc is wrong. The source code is attached. The example is: 1x1 domain, 10x10 mesh, Quad4, fixed bc on x=3D0, p=3D1 on x=3D1, consider periodic bc of ux with node (0.5,0) and (0.8,0). 
From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20060224 17:46:09

I just installed petsc 2.3.1 and built libmesh with it. At least for the case of real numbers without SLEPc the current libMesh source required no modification. I am a little shocked... Anyway, I figured I'd document the Endeavour. See http://libmesh.sourceforge.net/wiki/index.php/Installation#libmesh_from_ cvs.2C_PETSc_2.3.1.2C_CentOS_4.2  Benjamin S. Kirk benjamin.kirk@...=20 NASA/JSC, Applied Aerosciences & CFD Branch 2814839491 
From: John Peterson <peterson@cf...>  20060224 16:53:16

jcch@... writes: > > I would like to understand the reason for including the neighbors > in the send_list. I'll give it a shot.... Sometimes (in continuous Galerkin) you need to compute interelement fluxjumps for the Kelly error indicator, which requires the neighbors. I'm not sure about DG. J 
From: <jcch@ma...>  20060224 16:22:46

I am currently using libmesh for solving the 3D Euler equations using discontinuous finite elements. The dof_map members void DofMap::distribute_dofs_var_major (MeshBase& mesh) void DofMap::distribute_dofs_node_major (MeshBase& mesh) append to the send_list the DOFs from elements that live on neighboring processors (that are neighbors of the elements on the local processor) but the solution vector only allocates local space for the local DOFs. I would like to understand the reason for including the neighbors in the send_list. 
From: li pan <li76pan@ya...>  20060224 11:54:38

thanks for your explaination. I'm dealing with nonlinear problem with quadratic strain gradient in my M(i,j) . So I might use QGauss(SECOND). best regards pan Roy Stogner writes: > On Thu, 23 Feb 2006, li pan wrote: > > Ben, John, is our QGauss documentation wrong? It says, "Gauss > quadrature rules of order p have the property of integrating > polynomials of degree 2p1 exactly.", but that's not true. That > sounds like the theorem for a 1D Gauss rule with p points, but looking > in the code, we aren't using the Order argument to choose a number of > points, we're using it to choose what degree of polynomial to exactly > integrate. The documentation is wrong. According to the source code, it's something like this: order p = # of points 2p1 ===== =============== ==== CONSTANT,FIRST 1 1 SECOND 3 5 THIRD 4 7 I think it's based on computing entries in the "mass matrix" M_ij = \int (phi_i * phi_j) dx. So, I think a "SECOND" order rule will compute M_ij exactly if both phi_i and phi_j are quadratic polynomials. A "THIRD" order rule will compute M_ij exactly if both phi_i and phi_j are cubic polys. etc. __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com 
From: Roy Stogner <roystgnr@ic...>  20060223 21:39:09

On Thu, 23 Feb 2006, John Peterson wrote: > According to the source code, it's something like this: > > order p = # of points 2p1 > ===== =============== ==== > CONSTANT,FIRST 1 1 > SECOND 3 5 > THIRD 4 7 No, no, you're looking at triangle quadratures  the "2p1" rule only applies on line segments. Even in the code for triangles, it looks to me like the rules labeled "exact for quadratics", "exact for cubics", etc. should be exact when integrating quadratics, cubics, etc. not products of them. I'm committing a CVS change to the quadrature_gauss.h comments to say so.  Roy 
From: John Peterson <peterson@cf...>  20060223 21:24:51

Roy Stogner writes: > On Thu, 23 Feb 2006, li pan wrote: > > Ben, John, is our QGauss documentation wrong? It says, "Gauss > quadrature rules of order p have the property of integrating > polynomials of degree 2p1 exactly.", but that's not true. That > sounds like the theorem for a 1D Gauss rule with p points, but looking > in the code, we aren't using the Order argument to choose a number of > points, we're using it to choose what degree of polynomial to exactly > integrate. The documentation is wrong. According to the source code, it's something like this: order p = # of points 2p1 ===== =============== ==== CONSTANT,FIRST 1 1 SECOND 3 5 THIRD 4 7 I think it's based on computing entries in the "mass matrix" M_ij = \int (phi_i * phi_j) dx. So, I think a "SECOND" order rule will compute M_ij exactly if both phi_i and phi_j are quadratic polynomials. A "THIRD" order rule will compute M_ij exactly if both phi_i and phi_j are cubic polys. etc. J 
From: Roy Stogner <roystgnr@ic...>  20060223 20:22:04

On Thu, 23 Feb 2006, li pan wrote: > thanks for your suggestion. I think I can try to do > the integration like with hex mesh. Another simple > question. If I'm using Tet4, and considering u v w > components (variables), and all the components have > order one, can I use QGauss(SECOND)? If you're integrating products of two linear functions, then any rule that exactly integrates quadratics will work for you  QGauss(SECOND) is the cheapest such rule. Of course, I don't know exactly what you're integrating  if you're just solving a Laplacian problem and only need to integrate products of gradients, then all your integrals except the boundary integrals are elementwise constant, and QGauss(CONSTANT) will give you the same results with one quadrature point instead of four. If you're solving a nonlinear problem or a problem with variable coefficients, then getting exact integrals from a quadrature rule may be impossible, and even if QGauss(SECOND) doesn't ruin your O(h^2) accuracy you may want to experiment with higher rules to get a better constant in front of that h^2. Ben, John, is our QGauss documentation wrong? It says, "Gauss quadrature rules of order p have the property of integrating polynomials of degree 2p1 exactly.", but that's not true. That sounds like the theorem for a 1D Gauss rule with p points, but looking in the code, we aren't using the Order argument to choose a number of points, we're using it to choose what degree of polynomial to exactly integrate.  Roy 
From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20060223 15:23:36

I'd be happy to help here too, but there is another option you could consider. PETSc says they can interface directly with UMFPACK. In that case, you simply need to build PETSc such that it uses UMFPACK internally and then you should have access to (at leas a subset of) the UMFPACK sparse directly. The standard PETSc allows for "pc_type direct" which does a dense direct LU factorization and then uses it as the "preconditioner" in an iterative method (which will then converge in 1 iteration). This only works on serial problems, mind you, but it looks like UMFPACK is serial as well. If you install PETSc with UMFPACK support I would guess you'd get something like "pc_type umfpack" which would then build the "preconditioner" using UMFPACK's Unsymmetric sparse multifrontal method. The only reason I mention it is because it may be the path of least resistance, especially if UMFPACK is relegated to serial problems only. Adding a new solver interface is fairly straightforward but time consuming, and it is a lot of work to go through for a solver which will only work on serial problems. Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of John Peterson Sent: Thursday, February 23, 2006 8:47 AM To: Blome Mark Cc: Libmeshusers@... Subject: [Libmeshusers] direct solver Blome Mark writes: > > I would like to use a direct sparse matrix solver package within > libmesh (for example umfpack, pardiso or thelike). > Am I right that I have to implement an interface for the solver > package by deriving from the class LinearSolver (linear_solver.h)? > Are there any pitfalls I should be aware of ? How do I include the > solver package into the libmesh Makefile so it will get linked > correctly? Did anybody alreay implement a direct solver interface ? You might get some hints from the Petsc and Laspack implementations of the LinearSolver interface. Some things you will probably need are: 1.) UMF_Matrix derived from SparseMatrix 2.) UMF_Vector derived from NumericVector 3.) UMF_Interface as you mentioned This sounds like an interesting project. If UMF is open source, and small enough, we can distribute it directly in libmesh's contrib directory. Otherwise, you will need 4.) configure tests, etc. I think we can help with 1.)  4.) if you do end up trying to do this. J  This SF.Net email is sponsored by xPML, a groundbreaking scripting language that extends applications into web and mobile media. Attend the live webcast and join the prime developer group breaking into this new coding territory! http://sel.asus.falkag.net/sel?cmd=3Dlnk&kid=3D110944&bid=3D241720&dat=3D= 121642 _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: John Peterson <peterson@cf...>  20060223 14:47:22

Blome Mark writes: > > I would like to use a direct sparse matrix solver package within > libmesh (for example umfpack, pardiso or thelike). > Am I right that I have to implement an interface for the solver > package by deriving from the class LinearSolver (linear_solver.h)? > Are there any pitfalls I should be aware of ? How do I include the > solver package into the libmesh Makefile so it will get linked > correctly? Did anybody alreay implement a direct solver interface ? You might get some hints from the Petsc and Laspack implementations of the LinearSolver interface. Some things you will probably need are: 1.) UMF_Matrix derived from SparseMatrix 2.) UMF_Vector derived from NumericVector 3.) UMF_Interface as you mentioned This sounds like an interesting project. If UMF is open source, and small enough, we can distribute it directly in libmesh's contrib directory. Otherwise, you will need 4.) configure tests, etc. I think we can help with 1.)  4.) if you do end up trying to do this. J 
From: li pan <li76pan@ya...>  20060223 13:27:55

Hi Roy and John, thanks for your suggestion. I think I can try to do the integration like with hex mesh. Another simple question. If I'm using Tet4, and considering u v w components (variables), and all the components have order one, can I use QGauss(SECOND)? best regards pan On Wed, 22 Feb 2006, li pan wrote: > I realized all the example which include an assemble() > function are written for square or hex. Has anyone got > a example for tetrahedron? It's not clear for me how > to do volume integration and surface integration in > the case of tetrahedron. It's exactly the same as area integration and edge integration for the 2D cases. I've got an altered 2D+3D example 15 that I could send you if you want (or I could finally get around to cleaning it up and committing it to CVS...), but it's not going to tell you much  the assemble function is exactly the same no matter whether the mesh is CloughTocher triangles, Hermite rectangles, or Hermite cubes.  Roy __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com 
From: Blome Mark <blome@au...>  20060223 12:49:30

Hi everybody, I would like to use a direct sparse matrix solver package within libmesh = (for example umfpack, pardiso or thelike).=20 Am I right that I have to implement an interface for the solver package = by deriving from the class LinearSolver (linear_solver.h)? Are there any = pitfalls I should be aware of ? How do I include the solver package into = the libmesh Makefile so it will get linked correctly? Did anybody alreay = implement a direct solver interface ?=20 Thanks for any help on this, Mark=20 
From: Shengli Xu <shengli.xu.xu@gm...>  20060223 10:09:43

Dear all, I need to consider the pressure load. But I don't know how to form the inpu= t file such as .unv or .xda considering the pressure. The example laplace.tar.bz2 supported by Ondrej Certik on the libmesh wiki. But there i= s an error when run ./mesh. the error is :  ... Info : Writing mesh file 't1.msh' Info : Mesh 2D complete (0.028995 s) Info : Writing mesh file 't1.msh' Info : 411 nodes Info : 1537 elements Info : Wrote mesh file 't1.msh' reading... Traceback ( most recent call last): File "./msh2libmesh", line 28, in ? m.readmsh("t1.msh",b, False,False) TypeError: readmsh() takes at most 4 arguments (5 given) I doubt that "m.readmsh("t1.msh",b,False,False)" just has 4 argumets. why i= t display "5 given"? best ShengliXu 
From: John Peterson <peterson@cf...>  20060222 16:49:30

li pan writes: > dear all, > I realized all the example which include an assemble() > function are written for square or hex. Has anyone got > a example for tetrahedron? It's not clear for me how > to do volume integration and surface integration in > the case of tetrahedron. It doesn't make any difference. The code is the same for hexes or tets. See example 4 if you want to see a 3D example. J 
From: Roy Stogner <roystgnr@ic...>  20060222 16:46:47

On Wed, 22 Feb 2006, li pan wrote: > I realized all the example which include an assemble() > function are written for square or hex. Has anyone got > a example for tetrahedron? It's not clear for me how > to do volume integration and surface integration in > the case of tetrahedron. It's exactly the same as area integration and edge integration for the 2D cases. I've got an altered 2D+3D example 15 that I could send you if you want (or I could finally get around to cleaning it up and committing it to CVS...), but it's not going to tell you much  the assemble function is exactly the same no matter whether the mesh is CloughTocher triangles, Hermite rectangles, or Hermite cubes.  Roy 
From: li pan <li76pan@ya...>  20060222 16:23:05

dear all, I realized all the example which include an assemble() function are written for square or hex. Has anyone got a example for tetrahedron? It's not clear for me how to do volume integration and surface integration in the case of tetrahedron. best pan __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com 
From: Steffen Petersen <steffen.petersen@tu...>  20060221 18:28:34

Hello David, I will commit an adjusted example within the next days. Note that some rather useful functionality is only available with the latest slepc version (e.g. computation of the smallest eigenvalues within the spectrum), so I recommend to use the 2.3.0 versions of petsc and slepc. Steffen > Hi Steffen and Ondrej, > > I read your previous message that you've extended the slepc interface > to support generalized eigen problems. I'd like to get the details > from you. I'm new to libmesh and I'm wondering if you can post some > new tutorial examples on how to solve eigen problems using libmesh > (eg. solve 1D or 2D simple harmonic oscillator)? That will a great > help. > > Ondrej, I read about your message about your implementation of libmesh > and pysparse's jacobidavidson solver. That sounds really good. Could > you please give me the details too? I really appreciate it. > > David > 
From: Shengli Xu <shengli.xu.xu@gm...>  20060221 07:28:31

Hello, everyone It is something about how to impose boundary condition. The question is: Firstly, dof A is fixed by the penalty method. dof B is also needed to fixed. But dof B is imposed the periodic boundary condition to dof A. I just want to know is it right in libMesh. Another question: Are there any other method in libMesh to impose fixed boundary condition beside the penalty meshod? Best Regards! ShengliXu 
From: David Xu <dxu@my...>  20060220 22:48:43

Hi Steffen and Ondrej, I read your previous message that you've extended the slepc interface to support generalized eigen problems. I'd like to get the details from you. I'm new to libmesh and I'm wondering if you can post some new tutorial examples on how to solve eigen problems using libmesh (eg. solve 1D or 2D simple harmonic oscillator)? That will a great help. Ondrej, I read about your message about your implementation of libmesh and pysparse's jacobidavidson solver. That sounds really good. Could you please give me the details too? I really appreciate it. David 
From: Roy Stogner <roystgnr@ic...>  20060219 16:55:32

On Sun, 19 Feb 2006, Shengli Xu wrote: > Quad9 is used in my program. Why the element type in the gmv result file is > Quad4? Is it correct? It's not correct, but it may be the best we can do. GMV doesn't support Quad9 elements, so for GMV output we refine each of them into four Quad4 elements. This usually looks okay (just a little less smooth than it ought to) for uniform meshes, but it can make solutions on adaptive meshes look discontinuous when they really aren't.  Roy 
From: Shengli Xu <shengli.xu.xu@gm...>  20060219 05:46:27

Hello, Quad9 is used in my program. Why the element type in the gmv result file is Quad4? Is it correct? Best Regards! ShengliXu 